IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Nov 16, 2021, 8:28:36 AM (5 years ago)
Author:
eugene
Message:

clean up imfit, fix imfit-trail, add imfit tests

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/ipp-20211108/Ohana/src/opihi/cmd.astro/imfit.c

    r41666 r41916  
    33int imfit (int argc, char **argv) {
    44
    5   int i, j, N, Npts, Save, VERBOSE;
    6   int sx, sy, nx, ny, Nx, Ny;
    7   float chisq, ochisq, dchisq, Gain, RDnoise, SatThreshold;
    8   opihi_flt *x, *y, *z, *dz;
    9   float *V;
     5  int N;
    106  Buffer *buf;
    117
    12   Save = FALSE;
     8  char *Save = NULL;
    139  if ((N = get_argument (argc, argv, "-save"))) {
    1410    remove_argument (N, &argc, argv);
    15     Save = TRUE;
    16   }
    17 
    18   // int ShapeVariation = FALSE;
    19   // if ((N = get_argument (argc, argv, "-shapes"))) {
    20   //   remove_argument (N, &argc, argv);
    21   //   ShapeVariation = TRUE;
    22   // }
    23 
    24   SatThreshold = 0xffff;
     11    Save = strcreate (argv[N]);
     12    remove_argument (N, &argc, argv);
     13  }
     14
     15  int Insert = FALSE;
     16  if ((N = get_argument (argc, argv, "-insert"))) {
     17    remove_argument (N, &argc, argv);
     18    Insert = TRUE;
     19    if (Save) { gprint (GP_ERR, "-save and -insert are mutually exclusive\n"); free (Save); return (FALSE); }
     20  }
     21
     22  int SatThreshold = 0xffff;
    2523  if ((N = get_argument (argc, argv, "-sat"))) {
    2624    remove_argument (N, &argc, argv);
    … …  
    3028
    3129  /* Gain in e/DN */
    32   Gain = 1.0;
     30  float Gain = 1.0;
    3331  if ((N = get_argument (argc, argv, "-gain"))) {
    3432    remove_argument (N, &argc, argv);
    … …  
    3836
    3937  /* RD noise in DN */
    40   RDnoise = 0.0;
     38  float RDnoise = 0.0;
    4139  if ((N = get_argument (argc, argv, "-rdnoise"))) {
    4240    remove_argument (N, &argc, argv);
    … …  
    4543  }
    4644
    47   VERBOSE = FALSE;
     45  int VERBOSE = FALSE;
    4846  if ((N = get_argument (argc, argv, "-v"))) {
    4947    remove_argument (N, &argc, argv);
    … …  
    5149  }
    5250
    53   /* set fitting function */
     51  /* set fitting function : defines par, Npar, fitfunc, etc globals (imfit.h) */
    5452  fgauss_setup ("fgauss");
    5553  if ((N = get_argument (argc, argv, "-func"))) {
    5654    fitfunc = NULL;
    5755    remove_argument (N, &argc, argv);
    58     fgauss_setup (argv[N]);
    59     pgauss_setup (argv[N]);
    60     pgauss_psf_setup (argv[N]);
     56    fgauss_setup (argv[N]); // OK
     57    pgauss_setup (argv[N]); // OK
    6158    sgauss_setup (argv[N]);
    62     sgauss_psf_setup (argv[N]);
    63     qgauss_setup (argv[N]);
    64     qgauss_psf_setup (argv[N]);
    65     qfgauss_setup (argv[N]);
     59    qgauss_setup (argv[N]); // OK
     60    qfgauss_setup (argv[N]); // OK
    6661    qrgauss_setup (argv[N]);
    6762    trail_setup (argv[N]);
     63    pgauss_psf_setup (argv[N]);
     64    sgauss_psf_setup (argv[N]);
     65    qgauss_psf_setup (argv[N]);
    6866    if (fitfunc == NULL) {
    6967      gprint (GP_ERR, "unknown function %s\n", argv[N]);
     68      FREE (Save);
    7069      return (FALSE);
    7170    }
    … …  
    7473
    7574  if (argc != 6) {
    76     gprint (GP_ERR, "USAGE: imfit <buffer> sx sy nx ny\n");
     75    gprint (GP_ERR, "USAGE: imfit <buffer> Xo Yo dX dY\n");
     76    gprint (GP_ERR, "options: [-save buffer] [-insert] [-sat value] [-gain value] [-rdnoise value] [-v] [-func option]\n");
     77    gprint (GP_ERR, "   (Xo,Yo) : center\n");
     78    gprint (GP_ERR, "   (dX,dY) : window size\n");
     79    FREE (Save);
    7780    return (FALSE);
    7881  }
    7982
    8083  /* non-optional arguments */
    81   if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) return (FALSE);
    82   sx = atof (argv[2]);
    83   sy = atof (argv[3]);
    84   nx = atof (argv[4]);
    85   ny = atof (argv[5]);
    86   Nx = buf[0].matrix.Naxis[0];
    87   Ny = buf[0].matrix.Naxis[1];
    88 
    89   /* check if region is valid */
    90   if (sx + 0.5*nx < 0) goto range;
    91   if (sy + 0.5*ny < 0) goto range;
    92   if (sx + 0.5*nx >= Nx) goto range;
    93   if (sy + 0.5*ny >= Ny) goto range;
     84  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) { FREE (Save); return (FALSE); }
     85  int Xo = atof (argv[2]);
     86  int Yo = atof (argv[3]);
     87  int dX = atof (argv[4]);
     88  int dY = atof (argv[5]);
     89  int Nx = buf[0].matrix.Naxis[0];
     90  int Ny = buf[0].matrix.Naxis[1];
     91
     92  int sx = Xo - dX/2;
     93  int sy = Yo - dY/2;
     94
     95  /* check if region is valid (center must be in range of image pixels) */
     96  if (Xo < 0) goto range;
     97  if (Yo < 0) goto range;
     98  if (Xo >= Nx) goto range;
     99  if (Yo >= Ny) goto range;
     100
     101  // image value in DN
     102  // rdnoise in DN
     103  // sigma_DN^2 = sigma_e^2 / gain^2
     104  // sigma_e^2  = Ne
     105  // sigma_e^2  = DN * gain
     106  // sigma_DN^2 = DN * gain / gain^2 = DN / gain
     107
     108  if (Insert) {
     109    float *Vi = (float *)buf[0].matrix.buffer;
     110    for (int j = 0; j < dY; j++) {
     111      for (int i = 0; i < dX; i++) {
     112        float vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL);
     113        Vi[(i+sx)+(j+sy)*Nx] += vf;
     114      }
     115    }
     116    return TRUE;
     117  }
    94118
    95119  /* convert array z[x,y] to x[i], y[i], z[i] */
    96120  N = 0;
    97   Npts = nx*ny;
    98   ALLOCATE (x,  opihi_flt, 2*Npts);
    99   ALLOCATE (y,  opihi_flt, 2*Npts);
    100   ALLOCATE (z,  opihi_flt, 2*Npts);
    101   ALLOCATE (dz, opihi_flt, 2*Npts);
    102   for (j = 0; j < ny; j++) {
     121  int Npts = dX*dY;
     122  ALLOCATE_PTR (x,  opihi_flt, 2*Npts);
     123  ALLOCATE_PTR (y,  opihi_flt, 2*Npts);
     124  ALLOCATE_PTR (z,  opihi_flt, 2*Npts);
     125  ALLOCATE_PTR (dz, opihi_flt, 2*Npts);
     126  for (int j = 0; j < dY; j++) {
    103127    if (j + sy < 0) continue;
    104128    if (j + sy >= Ny) continue;
    105     V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
    106     for (i = 0; i < nx; i++) {
     129    float *V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
     130    for (int i = 0; i < dX; i++) {
    107131      if (i + sx < 0) continue;
    108132      if (i + sx >= Nx) continue;
    109133      if (*V > SatThreshold) goto next;
    110       dz[N] = (SQ(RDnoise) + *V/Gain);
     134      dz[N] = (SQ(RDnoise) + MAX(0.0, *V/Gain)); // treat negative pixels as pure read noise
    111135      if (dz[N] <= 0) goto next;
    112136      dz[N] = 1.0 / dz[N];
    … …  
    122146
    123147  /* run fit routine */
    124   ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
    125   dchisq = ochisq;
    126   chisq  = ochisq;
    127 
    128 //for (i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
    129 
    130   for (i = 0; (i < 25); i++) {
     148  float ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
     149  float dchisq = ochisq;
     150  float chisq  = ochisq;
     151
     152  int Niter = 0;
     153  // for (int i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
     154
     155  // keep iterating if chisq is increasing or
     156  for (Niter = 0; (Niter < 25) && ((dchisq <= 0.0) || (dchisq > 0.1*(Npts - Npar))); Niter++) {
    131157    chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
    132158    dchisq = ochisq - chisq;
     159    // fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq);
    133160    ochisq = chisq;
    134     fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq);
    135161  } 
    136   set_int_variable ("Niter",  i);
     162  set_int_variable ("Niter",  Niter);
    137163
    138164  /** create output image (keep in sky) **/
    … …  
    141167    float *Vi, *Vo, vr, vf;
    142168
    143     if ((out = SelectBuffer ("out",   ANYBUFFER, TRUE)) == NULL) return (FALSE);
     169    if ((out = SelectBuffer (Save, ANYBUFFER, TRUE)) == NULL) { free (Save); return (FALSE); }
    144170    free (out[0].header.buffer);
    145171    free (out[0].matrix.buffer);
    146172
    147173    strcpy (out[0].file, "(empty)");
    148     if (!CreateBuffer (out, 2*nx, 2*ny, -32, 0.0, 1.0)) return FALSE;
    149 
    150     /* four panels: 1) raw image. 2) fit  3) raw - fit   4) ?? */
     174    if (!CreateBuffer (out, 2*dX, 2*dY, -32, 0.0, 1.0)) { free (Save); return FALSE; }
     175
     176    /* four panels: 1) raw image. 2) fit  3) raw - fit  4) absolute deviation */
    151177    Vi = (float *)buf[0].matrix.buffer;
    152178    Vo = (float *)out[0].matrix.buffer;
    153     for (j = 0; j < ny; j++) {
    154       for (i = 0; i < nx; i++) {
     179    for (int j = 0; j < dY; j++) {
     180      for (int i = 0; i < dX; i++) {
    155181        vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL);
    156182        vr = Vi[(i+sx)+(j+sy)*Nx];
    157         Vo[(i   )+(j   )*2*nx] = vr;
    158         Vo[(i+nx)+(j   )*2*nx] = vf;
    159         Vo[(i   )+(j+ny)*2*nx] = vr - vf + *sky;
    160         Vo[(i+nx)+(j+ny)*2*nx] = fabs(vr-vf) + *sky;
     183        Vo[(i   )+(j   )*2*dX] = vr;
     184        Vo[(i+dX)+(j   )*2*dX] = vf;
     185        Vo[(i   )+(j+dY)*2*dX] = vr - vf + *sky;
     186        Vo[(i+dX)+(j+dY)*2*dX] = fabs(vr-vf) + *sky;
    161187      }
    162188    }
     189    free (Save);
    163190  }
    164191
    … …  
    169196
    170197  if (VERBOSE) {
    171     for (i = 0; i < Npar; i++) {
     198    for (int i = 0; i < Npar; i++) {
    172199      gprint (GP_ERR, "%g ", par[i]);
    173200    }
Note: See TracChangeset for help on using the changeset viewer.