IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Feb 28, 2022, 12:10:28 PM (4 years ago)
Author:
eugene
Message:

add Molweide projection; improvements to source fitting functions (imfit) for testing

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/opihi/cmd.astro/imfit.c

    r41666 r42078  
    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 minIter = 5;
     23  if ((N = get_argument (argc, argv, "-min-iter"))) {
     24    remove_argument (N, &argc, argv);
     25    minIter = atoi (argv[N]);
     26    remove_argument (N, &argc, argv);
     27  }
     28
     29  int SatThreshold = 0xffff;
    2530  if ((N = get_argument (argc, argv, "-sat"))) {
    2631    remove_argument (N, &argc, argv);
    … …  
    3035
    3136  /* Gain in e/DN */
    32   Gain = 1.0;
     37  float Gain = 1.0;
    3338  if ((N = get_argument (argc, argv, "-gain"))) {
    3439    remove_argument (N, &argc, argv);
    … …  
    3843
    3944  /* RD noise in DN */
    40   RDnoise = 0.0;
     45  float RDnoise = 0.0;
    4146  if ((N = get_argument (argc, argv, "-rdnoise"))) {
    4247    remove_argument (N, &argc, argv);
    … …  
    4550  }
    4651
    47   VERBOSE = FALSE;
     52  int VERBOSE = FALSE;
    4853  if ((N = get_argument (argc, argv, "-v"))) {
    4954    remove_argument (N, &argc, argv);
    … …  
    5156  }
    5257
    53   /* set fitting function */
     58  /* set fitting function : defines par, Npar, fitfunc, etc globals (imfit.h) */
    5459  fgauss_setup ("fgauss");
    5560  if ((N = get_argument (argc, argv, "-func"))) {
    5661    fitfunc = NULL;
    5762    remove_argument (N, &argc, argv);
    58     fgauss_setup (argv[N]);
    59     pgauss_setup (argv[N]);
    60     pgauss_psf_setup (argv[N]);
     63    fgauss_setup (argv[N]); // OK
     64    pgauss_setup (argv[N]); // OK
    6165    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]);
     66    qgauss_setup (argv[N]); // OK
     67    rgauss_setup (argv[N]);
     68    qfgauss_setup (argv[N]); // OK
    6669    qrgauss_setup (argv[N]);
    6770    trail_setup (argv[N]);
     71    pgauss_psf_setup (argv[N]);
     72    sgauss_psf_setup (argv[N]);
     73    qgauss_psf_setup (argv[N]);
     74
     75    fgauss_pol_setup (argv[N]);
     76    rgauss_pol_setup (argv[N]);
     77
    6878    if (fitfunc == NULL) {
    6979      gprint (GP_ERR, "unknown function %s\n", argv[N]);
     80      FREE (Save);
    7081      return (FALSE);
    7182    }
    … …  
    7485
    7586  if (argc != 6) {
    76     gprint (GP_ERR, "USAGE: imfit <buffer> sx sy nx ny\n");
     87    gprint (GP_ERR, "USAGE: imfit <buffer> Xo Yo dX dY\n");
     88    gprint (GP_ERR, "options: [-save buffer] [-insert] [-sat value] [-gain value] [-rdnoise value] [-v] [-func option]\n");
     89    gprint (GP_ERR, "   (Xo,Yo) : center\n");
     90    gprint (GP_ERR, "   (dX,dY) : window size\n");
     91    FREE (Save);
    7792    return (FALSE);
    7893  }
    7994
    8095  /* 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;
     96  if ((buf = SelectBuffer (argv[1], OLDBUFFER, TRUE)) == NULL) { FREE (Save); return (FALSE); }
     97  int Xo = atof (argv[2]);
     98  int Yo = atof (argv[3]);
     99  int dX = atof (argv[4]);
     100  int dY = atof (argv[5]);
     101  int Nx = buf[0].matrix.Naxis[0];
     102  int Ny = buf[0].matrix.Naxis[1];
     103
     104  int sx = Xo - dX/2;
     105  int sy = Yo - dY/2;
     106
     107  /* check if region is valid (center must be in range of image pixels) */
     108  if (Xo < 0) goto range;
     109  if (Yo < 0) goto range;
     110  if (Xo >= Nx) goto range;
     111  if (Yo >= Ny) goto range;
     112
     113  // image value in DN
     114  // rdnoise in DN
     115  // sigma_DN^2 = sigma_e^2 / gain^2
     116  // sigma_e^2  = Ne
     117  // sigma_e^2  = DN * gain
     118  // sigma_DN^2 = DN * gain / gain^2 = DN / gain
     119
     120  if (Insert) {
     121    float *Vi = (float *)buf[0].matrix.buffer;
     122    for (int j = 0; j < dY; j++) {
     123      for (int i = 0; i < dX; i++) {
     124        float vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL);
     125        Vi[(i+sx)+(j+sy)*Nx] += vf;
     126      }
     127    }
     128    return TRUE;
     129  }
    94130
    95131  /* convert array z[x,y] to x[i], y[i], z[i] */
    96132  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++) {
     133  int Npts = dX*dY;
     134  ALLOCATE_PTR (x,  opihi_flt, 2*Npts);
     135  ALLOCATE_PTR (y,  opihi_flt, 2*Npts);
     136  ALLOCATE_PTR (z,  opihi_flt, 2*Npts);
     137  ALLOCATE_PTR (dz, opihi_flt, 2*Npts);
     138  for (int j = 0; j < dY; j++) {
    103139    if (j + sy < 0) continue;
    104140    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++) {
     141    float *V = (float *)(buf[0].matrix.buffer) + (j+sy)*buf[0].matrix.Naxis[0] + sx;
     142    for (int i = 0; i < dX; i++) {
    107143      if (i + sx < 0) continue;
    108144      if (i + sx >= Nx) continue;
    109       if (*V > SatThreshold) goto next;
    110       dz[N] = (SQ(RDnoise) + *V/Gain);
     145      if (*V > SatThreshold) goto next; // skip pixels above threshold
     146      if (!isfinite(*V)) goto next; // skip nan pixels
     147      dz[N] = (SQ(RDnoise) + MAX(0.0, *V/Gain)); // treat negative pixels as pure read noise
    111148      if (dz[N] <= 0) goto next;
    112149      dz[N] = 1.0 / dz[N];
    … …  
    122159
    123160  /* 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++) {
     161  float ochisq = mrq2dinit (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
     162  float dchisq = ochisq;
     163  float chisq  = ochisq;
     164
     165  int Niter = 0;
     166  // for (int i = 0; (i < 25) && ((dchisq <= 0.0) || (dchisq > 0.01*(Npts - Npar))); i++) {
     167
     168  // keep iterating if chisq is increasing or
     169  for (Niter = 0; (Niter < 25) && ((Niter < minIter) || ((dchisq <= 0.0) || (dchisq > 0.1*(Npts - Npar)))); Niter++) {
    131170    chisq = mrq2dmin (x, y, z, dz, Npts, par, Npar, fitfunc, VERBOSE);
    132171    dchisq = ochisq - chisq;
     172    // fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq);
    133173    ochisq = chisq;
    134     fprintf (stderr, "%f -> %f : %f\n", ochisq, chisq, dchisq);
    135174  } 
    136   set_int_variable ("Niter",  i);
     175  set_int_variable ("Niter",  Niter);
    137176
    138177  /** create output image (keep in sky) **/
    … …  
    141180    float *Vi, *Vo, vr, vf;
    142181
    143     if ((out = SelectBuffer ("out",   ANYBUFFER, TRUE)) == NULL) return (FALSE);
     182    if ((out = SelectBuffer (Save, ANYBUFFER, TRUE)) == NULL) { free (Save); return (FALSE); }
    144183    free (out[0].header.buffer);
    145184    free (out[0].matrix.buffer);
    146185
    147186    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) ?? */
     187    if (!CreateBuffer (out, 2*dX, 2*dY, -32, 0.0, 1.0)) { free (Save); return FALSE; }
     188
     189    /* four panels: 1) raw image. 2) fit  3) raw - fit  4) absolute deviation */
    151190    Vi = (float *)buf[0].matrix.buffer;
    152191    Vo = (float *)out[0].matrix.buffer;
    153     for (j = 0; j < ny; j++) {
    154       for (i = 0; i < nx; i++) {
     192    for (int j = 0; j < dY; j++) {
     193      for (int i = 0; i < dX; i++) {
    155194        vf = fitfunc ((float)(i+sx), (float)(j+sy), par, Npar, NULL);
    156195        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;
     196        Vo[(i   )+(j   )*2*dX] = vr;
     197        Vo[(i+dX)+(j   )*2*dX] = vf;
     198        Vo[(i   )+(j+dY)*2*dX] = vr - vf + *sky;
     199        Vo[(i+dX)+(j+dY)*2*dX] = fabs(vr-vf) + *sky;
    161200      }
    162201    }
     202    free (Save);
    163203  }
    164204
    … …  
    169209
    170210  if (VERBOSE) {
    171     for (i = 0; i < Npar; i++) {
     211    for (int i = 0; i < Npar; i++) {
    172212      gprint (GP_ERR, "%g ", par[i]);
    173213    }
Note: See TracChangeset for help on using the changeset viewer.