IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 20323


Ignore:
Timestamp:
Oct 21, 2008, 6:11:09 PM (18 years ago)
Author:
eugene
Message:

remove old GRID_V1 code; allow grid analysis to slowly anneal

Location:
trunk/Ohana/src/relphot
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/Ohana/src/relphot/include/relphot.h

    r20185 r20323  
    6363
    6464int    NLOOP;
     65int    NGRID;
    6566int    RESET;
    6667int    UPDATE;
    … …  
    7374int    MEAS_BAD;
    7475int    STAR_TOOFEW;
     76int    GRID_TOOFEW;
    7577int    IMAGE_TOOFEW;
    7678double IMAGE_GOOD_FRACTION;
    7779int    IMAGE_BAD;
    7880int    FREEZE_IMAGES;
     81int    FREEZE_MOSAICS;
    7982int    USE_GRID;
    8083char  *OUTROOT;

  • trunk/Ohana/src/relphot/src/ConfigInit.c

    r19897 r20323  
    2525
    2626  GetConfig (config, "STAR_CHISQ",             "%lf", 0, &STAR_CHISQ);
     27  GetConfig (config, "GRID_TOOFEW",            "%d",  0, &GRID_TOOFEW);
    2728  GetConfig (config, "STAR_TOOFEW",            "%d",  0, &STAR_TOOFEW);
    2829  GetConfig (config, "IMAGE_TOOFEW",           "%d",  0, &IMAGE_TOOFEW);

  • trunk/Ohana/src/relphot/src/GridOps.c

    r20186 r20323  
    11# include "relphot.h"
    22
     3typedef enum {
     4  GRID_FITTED,
     5  GRID_FROZEN,
     6  GRID_REFERENCE,
     7};
     8
    39static int     Ngrid;
    4 static float   *gridM;
    5 static float   *gridS;
    6 static int     *gridN;
     10static float   *gridM; // magnitude offset for this grid cell
     11static float   *gridS; // stdev of the magnitude offset for this grid cell
     12static int     *gridN; // number of stars used to measure the magnitude offset for this grid cell
     13static int     *gridV; // data mode for this cell: fitted, frozen, reference
    714static int      gridX;
    815static int      gridY;
    … …  
    3037static char *config;
    3138
    32 # ifdef GRID_V1
    3339void initGrid (int dX, int dY) {
    3440
    35   int i;
    36 
    37   /* define mosaic 2d correction grid */
    38   gridX = dX / RELPHOT_GRID_BINNING + 1;
    39   gridY = dY / RELPHOT_GRID_BINNING + 1;
    40   Ngrid = gridX * gridY;
    41 
    42   ALLOCATE (gridM, float, Ngrid);
    43   ALLOCATE (gridS, float, Ngrid);
    44   ALLOCATE (gridN, int,   Ngrid);
    45   bzero (gridM, Ngrid*sizeof(float));
    46   bzero (gridS, Ngrid*sizeof(float));
    47   bzero (gridN, Ngrid*sizeof(int));
    48 
    49 }
    50 # endif
    51 
    52 # ifdef GRID_V2
    53 void initGrid (int dX, int dY) {
    54 
    55   int i, N, ccdnum_max;
     41  int i, N, ccdnum_max, refX, refY, refBin;
    5642  char *p, field[64], line[256];
    5743  int *Fx, *Fy;  /* chip flip */
    … …  
    7258  ScanConfig (config, "NAXIS1", "%d", 1, &camera.Nx);
    7359  ScanConfig (config, "NAXIS2", "%d", 1, &camera.Ny);
     60
     61  ScanConfig (config, "REFCELL.X", "%d", 1, &refX);
     62  ScanConfig (config, "REFCELL.Y", "%d", 1, &refY);
    7463
    7564  // temporary storage
    … …  
    143132  ALLOCATE (gridS, float, Ngrid);
    144133  ALLOCATE (gridN, int,   Ngrid);
    145   bzero (gridM, Ngrid*sizeof(float));
    146   bzero (gridS, Ngrid*sizeof(float));
    147   bzero (gridN, Ngrid*sizeof(int));
    148 
    149 }
    150 # endif
     134  ALLOCATE (gridV, int,   Ngrid);
     135
     136  // the grid bins may have one of three possible states: fitted, frozen, reference
     137  // set the initial values to indicate that the bins are frozen
     138  for (i = 0; i < Ngrid; i++) {
     139    gridM[i] = 0.0;
     140    gridS[i] = 0.0;
     141    gridN[i] = 0;
     142    gridV[i] = GRID_FROZEN;
     143  }
     144
     145  // refBin is the index of the grid cell which is kept at 0.0 (all others are relative to this)
     146  refBin = refX + refY*gridX;
     147  gridV[refBin] = GRID_REFERENCE;
     148}
    151149
    152150void initGridBins (Catalog *catalog, int Ncatalog) {
    … …  
    208206}
    209207
    210 # ifdef GRID_V1
    211 int setGridMeasure (int meas, int cat, double X, double Y) {
    212 
    213   int ix, iy, i;
    214 
    215   ix = X / RELPHOT_GRID_BINNING;
    216   iy = Y / RELPHOT_GRID_BINNING;
    217   if (ix < 0) goto escape;
    218   if (iy < 0) goto escape;
    219   if (ix >= gridX) goto escape;
    220   if (iy >= gridY) goto escape;
     208int setGridMeasure (int meas, int cat, double X, double Y, int ccdnum) {
     209
     210  int ix, iy, Cx, Cy, i;
     211  double x, y;
     212
     213  /* X, Y are chip coords on chip ccdnum */
     214
     215  /* normalize X & Y */
     216  x = X;
     217  if (camera.Fx[ccdnum]) x = camera.Nx - X;
     218  y = Y;
     219  if (camera.Fy[ccdnum]) y = camera.Ny - Y;
     220
     221  /* grid coords on the chip */
     222  Cx = MIN (MAX ((x / camera.Nx) * RELPHOT_GRID_X, 0), RELPHOT_GRID_X - 1);
     223  Cy = MIN (MAX ((y / camera.Ny) * RELPHOT_GRID_Y, 0), RELPHOT_GRID_Y - 1);
     224 
     225  /* coordinates in the grid */
     226  ix = Cx + camera.Ox[ccdnum]*RELPHOT_GRID_X;
     227  iy = Cy + camera.Oy[ccdnum]*RELPHOT_GRID_Y;
    221228
    222229  i = ix + iy*gridX;
    223230
    224231  bin[cat][meas] = i;
    225   Xmeas[cat][meas] = X;
    226   Ymeas[cat][meas] = Y;
     232  Xmeas[cat][meas] = x + camera.Ox[ccdnum]*camera.Nx;
     233  Ymeas[cat][meas] = y + camera.Oy[ccdnum]*camera.Ny;
    227234  clist[i][Nlist[i]] = cat;
    228235  mlist[i][Nlist[i]] = meas;
    … …  
    236243  return (TRUE);
    237244
    238 escape:
    239245  fprintf (stderr, "error: star out of grid\n");
    240246  exit (1);
    241247}
    242 # endif
    243 
    244 # ifdef GRID_V2
    245 int setGridMeasure (int meas, int cat, double X, double Y, int ccdnum) {
    246 
    247   int ix, iy, Cx, Cy, i;
    248   double x, y;
    249 
    250   /* X, Y are chip coords on chip ccdnum */
    251 
    252   /* normalize X & Y */
    253   x = X;
    254   if (camera.Fx[ccdnum]) x = camera.Nx - X;
    255   y = Y;
    256   if (camera.Fy[ccdnum]) y = camera.Ny - Y;
    257 
    258   /* grid coords on the chip */
    259   Cx = MIN (MAX ((x / camera.Nx) * RELPHOT_GRID_X, 0), RELPHOT_GRID_X - 1);
    260   Cy = MIN (MAX ((y / camera.Ny) * RELPHOT_GRID_Y, 0), RELPHOT_GRID_Y - 1);
    261  
    262   /* coordinates in the grid */
    263   ix = Cx + camera.Ox[ccdnum]*RELPHOT_GRID_X;
    264   iy = Cy + camera.Oy[ccdnum]*RELPHOT_GRID_Y;
    265 
    266   i = ix + iy*gridX;
    267 
    268   bin[cat][meas] = i;
    269   Xmeas[cat][meas] = x + camera.Ox[ccdnum]*camera.Nx;
    270   Ymeas[cat][meas] = y + camera.Oy[ccdnum]*camera.Ny;
    271   clist[i][Nlist[i]] = cat;
    272   mlist[i][Nlist[i]] = meas;
    273 
    274   Nlist[i] ++;
    275   if (Nlist[i] == NLIST[i]) {
    276     NLIST[i] += 100;
    277     REALLOCATE (clist[i], int, NLIST[i]);
    278     REALLOCATE (mlist[i], int, NLIST[i]);
    279   }     
    280   return (TRUE);
    281 
    282   fprintf (stderr, "error: star out of grid\n");
    283   exit (1);
    284 }
    285 # endif
    286248
    287249float getMgrid (int meas, int cat) {
    … …  
    294256  if (i == -1) return (NAN);
    295257
     258  // during the grid annealing process, we skip over grid cells until they have enough
     259  // valid stars to be fitted
     260  if (gridV[i] == GRID_FROZEN) return (NAN);
     261
    296262  value = gridM[i];
    297263  return (value);
    … …  
    301267void setMgrid (Catalog *catalog) {
    302268
    303   int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ncal, Nrel, Nsys;
     269  int i, j, m, c, n, N, Nmax, Nbad, Nmos, Ncal, Nrel, Nsys, Nfit;
    304270  double *list, *dlist;
    305271  float Msys, Mrel, Mcal, Mmos;
    … …  
    315281  ALLOCATE (dlist, double, Nmax);
    316282
    317   Nbad = Ncal = Nmos = Nrel = Nsys = 0;
     283  Nbad = Ncal = Nmos = Nrel = Nsys = Nfit = 0;
    318284
    319285  for (i = 0; i < Ngrid; i++) {
    … …  
    356322    }
    357323
     324    // the reference Cell is forced to have a value of 0.0, and is never changed to GRID_FITTED
     325    if (gridV[i] == GRID_REFERENCE) {
     326      gridM[i] = 0.0;
     327      gridS[i] = 0.0;
     328      gridN[i] = N;
     329      continue;
     330    }
     331
     332    // until we have enough valid measurements on this grid cell, skip it
     333    if (N < GRID_TOOFEW) {
     334      gridV[i] = GRID_FROZEN;
     335      continue;
     336    }
     337
    358338    liststats (list, dlist, N, &stats);
    359339    gridM[i] = stats.mean;
    360340    gridS[i] = stats.sigma;
    361341    gridN[i] = N;
    362   }
    363 
    364   fprintf (stderr, "grid cells having too few measurements (Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Nbad, Nmos, Ncal, Nrel, Nsys);
     342    gridV[i] = GRID_FITTED;
     343    Nfit++;
     344  }
     345
     346  fprintf (stderr, "%d of %d grid cells fitted (+ reference cell) (Nbad: %d, Nmos: %d, Ncal: %d, Nrel: %d, Nsys: %d)\n", Nfit, Ngrid, Nbad, Nmos, Ncal, Nrel, Nsys);
    365347
    366348  free (list);
    … …  
    531513  gfits_free_matrix (&matrix);
    532514
    533 # ifdef GRID_V1
    534   /* calculate pixel values for each CCD pixel, write out CCD images */
    535   /* grid pixels are in RA,DEC coords, transform to image and interpolate */
    536   for (i = 0; i < Nimage; i++) {
    537     image = getimage (imlist[i]);
    538     pname = GetPhotcodeNamebyCode (image[0].photcode);
    539 
    540     /* this is kind of bogus... */
    541     /* pname is CAMERA.FILTER.CCD, grab the CCD */
    542     p = strrchr (pname, '.');
    543     if (p == (char *) NULL) {
    544       fprintf (stderr, "error parsing photcode %s\n", pname);
    545       exit (2);
    546     }
    547     p++;
    548     sprintf (ccdname, "ccd%s", p);
    549 
    550     gfits_modify (&theader, "EXTNAME", "%s", 1, ccdname);
    551     gfits_modify (&theader, "FILTER", "%s", 1, photcode[0].name);
    552     gfits_modify (&theader, "PHOTCODE", "%s", 1, pname);
    553     gfits_modify (&theader, "NX", "%d", 1, image[i].NX);
    554     gfits_modify (&theader, "NY", "%d", 1, image[i].NY);
    555     write_coords (&theader, &image[0].coords);
    556 
    557     Nx = 2 * image[i].NX / RELPHOT_GRID_BINNING;
    558     Ny = 2 * image[i].NY / RELPHOT_GRID_BINNING;
    559     theader.Naxis[0] = Nx;
    560     theader.Naxis[1] = Ny;
    561     gfits_modify (&theader, "NAXIS1", "%d", 1, Nx);
    562     gfits_modify (&theader, "NAXIS2", "%d", 1, Ny);
    563     gfits_create_matrix  (&theader, &matrix);
    564 
    565     InterpolateGrid ((float *)matrix.buffer, Nx, Ny, &image[0].coords, &refmosaic[0].coords);
    566     gfits_fwrite_header (f, &theader);
    567     gfits_fwrite_matrix (f, &matrix);
    568     gfits_free_matrix (&matrix);
    569   }
    570 # endif
    571 
    572 # ifdef GRID_V2
    573515  /* calculate value for each CCD pixel, write out CCD images */
    574516  /* grid pixels are tied to detector pixels, but are flipped to match focal plane */
    … …  
    610552    gfits_free_matrix (&matrix);
    611553  }
    612 # endif
    613554
    614555  free (filename);

  • trunk/Ohana/src/relphot/src/ImageOps.c

    r20187 r20323  
    178178    if (measure[0].t > stop[i]) continue;
    179179   
    180 # ifdef GRID_V2 /* this section is added to support GridOps.v2.c */
    181180    if (USE_GRID) {
    182181
    … …  
    199198         constructed, there will be null values for undefined ccdnums */
    200199
    201 # if (0)
    202       /* add this measurement to the grid cell for this chip */
    203       ave = measure[0].averef;
    204       ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
    205       dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
    206        
    207       /* X,Y always positive-definite in range 0,0 - dX, dY */
    208       RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
    209 # endif
     200      // old code to add this measurement to the grid cell for this chip
     201      // ave = measure[0].averef;
     202      // ra  = catalog[cat].average[ave].R - measure[0].dR / 3600.0;
     203      // dec = catalog[cat].average[ave].D - measure[0].dD / 3600.0;
     204      // RD_to_XY (&X, &Y, ra, dec, &image[i].coords);
    210205
    211206      // XXX we can now use these values (but need to be careful about old formats)
    … …  
    214209      setGridMeasure (meas, cat, X, Y, ccdnum);
    215210    }
    216 # endif
    217211
    218212    bin[cat][meas] = i;

  • trunk/Ohana/src/relphot/src/MosaicOps.c

    r20192 r20323  
    320320
    321321  if (!MOSAIC_ZEROPT) return (FALSE);
     322  if (FREEZE_MOSAICS) return (FALSE);
    322323
    323324  image = getimages (&N);
    … …  
    453454  bzero (&stats, sizeof (StatType));
    454455  if (!MOSAIC_ZEROPT) return (stats);
     456  if (FREEZE_MOSAICS) return (stats);
    455457
    456458  ALLOCATE (list, double, Nmosaic);
    … …  
    481483  bzero (&stats, sizeof (StatType));
    482484  if (!MOSAIC_ZEROPT) return (stats);
     485  if (FREEZE_MOSAICS) return (stats);
    483486
    484487  ALLOCATE (list, double, Nmosaic);
    … …  
    523526  bzero (&stats, sizeof (StatType));
    524527  if (!MOSAIC_ZEROPT) return (stats);
     528  if (FREEZE_MOSAICS) return (stats);
    525529
    526530  ALLOCATE (list, double, Nmosaic);
    … …  
    551555
    552556  if (!MOSAIC_ZEROPT) return;
     557  if (FREEZE_MOSAICS) return;
    553558
    554559  if (VERBOSE) fprintf (stderr, "marking poor mosaics\n");
    … …  
    645650
    646651  if (!MOSAIC_ZEROPT) return;
     652  if (FREEZE_MOSAICS) return;
    647653
    648654  ALLOCATE (xlist, double, Nmosaic);

  • trunk/Ohana/src/relphot/src/StarOps.c

    r20192 r20323  
    8585        N++;
    8686      }
    87       if (N < STAR_TOOFEW) { /* too few measurements */
     87
     88      // when performing the grid analysis, STAR_TOOFEW will be set to 1;
     89
     90      if (N <= STAR_TOOFEW) { /* too few measurements */
    8891        catalog[i].average[j].code |= ID_STAR_FEW;
    8992        Nfew ++;
    … …  
    367370        N++;
    368371      }
    369       if (N < TOOFEW) continue;
     372      if (N <= TOOFEW) continue;
    370373
    371374      /* 3-sigma clip based on stats of inner 50% */

  • trunk/Ohana/src/relphot/src/args.c

    r20188 r20323  
    9393    remove_argument (N, &argc, argv);
    9494  }
     95  if ((N = get_argument (argc, argv, "-nloop"))) {
     96    remove_argument (N, &argc, argv);
     97    NLOOP = atof (argv[N]);
     98    remove_argument (N, &argc, argv);
     99  }
     100
     101  NGRID = 8;
     102  if ((N = get_argument (argc, argv, "-ngrid"))) {
     103    remove_argument (N, &argc, argv);
     104    NGRID = atof (argv[N]);
     105    remove_argument (N, &argc, argv);
     106  }
    95107
    96108  RESET = FALSE;
    … …  
    139151    remove_argument (N, &argc, argv);
    140152    FREEZE_IMAGES = TRUE;
     153  }
     154
     155  FREEZE_MOSAICS = FALSE;
     156  if ((N = get_argument (argc, argv, "-mosfreeze"))) {
     157    remove_argument (N, &argc, argv);
     158    FREEZE_MOSAICS = TRUE;
    141159  }
    142160

  • trunk/Ohana/src/relphot/src/bcatalog.c

    r20189 r20323  
    9090
    9191    // XXXX test : what checks do I need to make elsewhere to avoid problems here?
    92     if (Nm < STAR_TOOFEW) { /* enough measurements in band? */
     92    if (Nm <= STAR_TOOFEW) { /* enough measurements in band? */
    9393      Nmeasure -= Nm;
    9494      continue;

  • trunk/Ohana/src/relphot/src/relphot.c

    r17242 r20323  
    7070  }
    7171
     72  // if we are measuring the flat-field correction grid, we need to perform a number of iterations first:
     73  if (USE_GRID) {
     74      int star_toofew;
     75
     76      // until we finish the grid analysis, do not reject stars out-of-hand based on ID_STAR_FEW
     77      // XXX this is kind of poor: need to have a better distinctions about STAR_BAD in setMrel vs getMrel
     78      star_toofew = STAR_TOOFEW;
     79      for (i = 0; i < NGRID; i++) {
     80          STAR_BAD = ID_STAR_POOR;
     81          setMrel  (catalog, Ncatalog);
     82          STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
     83          setMgrid (catalog);
     84      }
     85      STAR_BAD  = ID_STAR_POOR | ID_STAR_FEW;
     86      STAR_TOOFEW = star_toofew;
     87  }
     88
    7289  /* determine fit values */
    7390  for (i = 0; i < NLOOP; i++) {
Note: See TracChangeset for help on using the changeset viewer.