IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
May 21, 2021, 10:03:06 AM (5 years ago)
Author:
eugene
Message:

major revision of the photometric flat-field correction / grid analysis; strip out code which uses astrometric information to define the grid

File:
1 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/relphot.20210521/src/MosaicOps.c

    r41559 r41603  
    5050# undef COMPARE
    5151
    52 }
    53 
    54 /* find mosaic frames (unique time periods & photcode name matches mosaic) */
    55 void initMosaics_old (Image *image, off_t Nimage) {
    56 
    57   off_t i, j, status, found, NMOSAIC, *MosaicN_IMAGE;
    58   unsigned int start, stop;
    59   char *pname;
    60 
    61   if (!MOSAIC_ZEROPT) return;
    62 
    63   Nmosaic = 0;
    64   NMOSAIC = 10;
    65   ALLOCATE (mosaic, Mosaic, NMOSAIC);
    66 
    67   ALLOCATE (MosaicToImage, off_t *, NMOSAIC);
    68   ALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
    69   ALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
    70 
    71   ALLOCATE (ImageToMosaic, off_t, Nimage); // mosaic to which image belongs
    72 
    73   /* a 'mosaic' in relphot is (unlike relastro) a virtual concept: there is no
    74    * entry in the image table that represents this mosaic.  Instead, it is an
    75    * internal construct that defines a group of related images
    76    */
    77 
    78   /* generate list of unique mosaics */
    79   for (i = 0; i < Nimage; i++) {
    80     ImageToMosaic[i] = -1;
    81 
    82     /* select valid mosaic images by photcode */
    83     pname = GetPhotcodeNamebyCode (image[i].photcode);
    84     status = strncmp (pname, MOSAICNAME, strlen (MOSAICNAME));
    85     if (status) continue;
    86 
    87     /* set image time range */
    88     start = image[i].tzero - MAX(0.01*image[i].trate*image[i].NY, 1);
    89     stop  = image[i].tzero + MAX(1.01*image[i].trate*image[i].NY, 1);
    90 
    91     /* find existing mosaic with this time range */
    92     found = FALSE;
    93     for (j = 0; !found && (j < Nmosaic); j++) {
    94       if (stop  < mosaic[j].start) continue;
    95       if (start > mosaic[j].stop)  continue;
    96       found = TRUE;
    97 
    98       // add reference from image to mosaic
    99       ImageToMosaic[i] = j;
    100 
    101       /* add image to mosaic image list */
    102       MosaicToImage[j][MosaicN_Image[j]] = i;
    103       MosaicN_Image[j] ++;
    104       if (MosaicN_Image[j] == MosaicN_IMAGE[j]) {
    105         MosaicN_IMAGE[j] += 10;
    106         REALLOCATE (MosaicToImage[j], off_t, MosaicN_IMAGE[j]);
    107       }
    108 
    109     }
    110     if (found) continue;
    111    
    112     /* a new mosaic, define ranges -- preserve the original values incase this image is not used */
    113     mosaic[Nmosaic].start     = start;
    114     mosaic[Nmosaic].stop      = stop;
    115     mosaic[Nmosaic].McalPSF   = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
    116     mosaic[Nmosaic].McalAPER  = 0.0; // note : mosaic stores only offsets relative to the original image values
    117     mosaic[Nmosaic].dMcal     = 0.0; // note : at the end, mosaic.Mcal is added back to the input images
    118     mosaic[Nmosaic].dMsys     = 0.0;
    119     mosaic[Nmosaic].McalChiSq = 0.0;// NAN or 0.0?
    120     mosaic[Nmosaic].flags     = image[i].flags;
    121     mosaic[Nmosaic].secz      = image[i].secz;
    122     mosaic[Nmosaic].photcode  = GetPhotcodeEquivCodebyCode (image[i].photcode);
    123 
    124     // XXX do we need to do something about flag consistency across a mosaic?
    125 
    126     /* add image to mosaic image list */
    127     MosaicN_IMAGE[Nmosaic] = 10;
    128     MosaicN_Image[Nmosaic] = 1;
    129     ALLOCATE (MosaicToImage[Nmosaic], off_t, MosaicN_IMAGE[Nmosaic]);
    130     MosaicToImage[Nmosaic][0] = i;
    131 
    132     // add reference from image to mosaic
    133     ImageToMosaic[i] = Nmosaic;
    134    
    135     Nmosaic ++;
    136     if (Nmosaic == NMOSAIC) {
    137       NMOSAIC += 10;
    138       REALLOCATE (mosaic,  Mosaic,  NMOSAIC);
    139       REALLOCATE (MosaicToImage,  off_t *, NMOSAIC);
    140       REALLOCATE (MosaicN_Image, off_t,   NMOSAIC);
    141       REALLOCATE (MosaicN_IMAGE, off_t,   NMOSAIC);
    142     }
    143   }
    144 
    145   // free this or not?
    146   free (MosaicN_IMAGE);
    147 
    148   initMosaicGrid (image, Nimage);
    149 
    150   fprintf (stderr, "matched %d images to %d mosaics\n", (int) Nimage, (int) Nmosaic);
    151   return;
    15252}
    15353
    … …  
    318218  free (startTimesMosaic);
    319219
    320   initMosaicGrid (subset, Nsubset);
     220  initMosaicMcal (subset, Nsubset);
    321221
    322222  fprintf (stderr, "matched %d images to %d mosaics, %d simple chips not matched to mosaics\n", (int) (Nsubset - Nsimple), (int) Nmosaic, (int) Nsimple);
    … …  
    460360
    461361  if (mergeMcal) {
    462     initMosaicGrid (image, Nimage);
     362    initMosaicMcal (image, Nimage);
    463363  }
    464364
    … …  
    511411}
    512412
     413// this function sets mosaic->coords to the median of the individual chips.  This
     414// coordinate frame is used by the parallel region analysis to assign exposures (mosaics)
     415// to specific machines by a single center (rather than individual chips)
    513416void setMosaicCenters (Image *image, off_t Nimage) {
    514417  OHANA_UNUSED_PARAM(Nimage);
    … …  
    599502    mosaic[i].coords.cdelt1 = 1.0 / 3600.0;
    600503    mosaic[i].coords.cdelt2 = 1.0 / 3600.0;
    601 
    602     mosaic[i].McalPSF   = 0.0;
    603     mosaic[i].McalAPER  = 0.0;
    604     mosaic[i].dMcal     = 0.0;
    605     mosaic[i].dMsys     = 0.0;
    606     mosaic[i].McalChiSq = 0.0;
    607504  }
    608505  return;
    609506}
    610507
    611 void initMosaicGrid (Image *image, off_t Nimage) {
     508void initMosaicMcal (Image *image, off_t Nimage) {
    612509  OHANA_UNUSED_PARAM(Nimage);
    613510
    614   /* find max dR, dD range for all mosaics */
    615   /* define mosaic.coords to cover dR, dD */
    616   /* send results to initGridBins */
    617 
    618   off_t i, j, m, NX, NY;
    619   int dXmax, dYmax;
    620   double dS, dX, dY;
    621   double R, D, Rmin, Rmax, Dmin, Dmax;
    622511  double McalPSF, McalAPER, dMcal, McalChiSq;
    623512
    624513  fprintf (stderr, "*** moving Mcal from image.Mcal to mosaic.Mcal ***\n");
    625514
    626   dXmax = dYmax = 0.0;
    627   for (i = 0; i < Nmosaic; i++) {
    628     Dmin = Rmin =  360.0;
    629     Dmax = Rmax = -360.0;
    630     dS = 0.0;
     515  for (off_t i = 0; i < Nmosaic; i++) {
    631516    McalPSF = McalAPER = dMcal = McalChiSq = 0;
    632     for (j = 0; j < MosaicN_Image[i]; j++) {
    633       m = MosaicToImage[i][j];
    634       NX = image[m].NX;
    635       NY = image[m].NY;
    636       dS += hypot (image[m].coords.cdelt1*image[m].coords.pc1_1, image[m].coords.cdelt1*image[m].coords.pc2_1);
    637 
    638       OhanaProjection proj = GetProjection (image[m].coords.ctype);
    639       if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
    640         XY_to_LM (&R, &D, 0.0, 0.0, &image[m].coords);
    641       } else {
    642         XY_to_RD (&R, &D, 0.0, 0.0, &image[m].coords);
    643       }
    644       Rmin = MIN (Rmin, R);
    645       Rmax = MAX (Rmax, R);
    646       Dmin = MIN (Dmin, D);
    647       Dmax = MAX (Dmax, D);
    648 
    649       if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
    650         XY_to_LM (&R, &D, (double) NX, 0.0, &image[m].coords);
    651       } else {
    652         XY_to_RD (&R, &D, (double) NX, 0.0, &image[m].coords);
    653       }
    654       Rmin = MIN (Rmin, R);
    655       Rmax = MAX (Rmax, R);
    656       Dmin = MIN (Dmin, D);
    657       Dmax = MAX (Dmax, D);
    658 
    659       if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
    660         XY_to_LM (&R, &D, (double) NX, (double) NY, &image[m].coords);
    661       } else {
    662         XY_to_RD (&R, &D, (double) NX, (double) NY, &image[m].coords);
    663       }
    664       Rmin = MIN (Rmin, R);
    665       Rmax = MAX (Rmax, R);
    666       Dmin = MIN (Dmin, D);
    667       Dmax = MAX (Dmax, D);
    668 
    669       if ((proj == PROJ_WRP) && !image[m].coords.mosaic) {
    670         XY_to_LM (&R, &D, 0.0, (double) NY, &image[m].coords);
    671       } else {
    672         XY_to_RD (&R, &D, 0.0, (double) NY, &image[m].coords);
    673       }
    674       Rmin = MIN (Rmin, R);
    675       Rmax = MAX (Rmax, R);
    676       Dmin = MIN (Dmin, D);
    677       Dmax = MAX (Dmax, D);
    678 
    679       // XXX : this probably does not handle mosaics at RA = 0,360 well
    680 
    681       /* we are using mosaic.Mcal, not image.Mcal. reset image.Mcal */
    682 
    683       // XXX: how does this work with UBERCAL?  We want to keep the Mcal values supplied by ubercal, but
    684       // solve for a single offset for each exposure (Mosaic.Mcal). 
    685       // we also want to keep the flat-field terms for each exposure (regardless of ubercal or not)
    686       // if it helps, note that ubercal uses a single zp per exposure, so the mean of those values is the same as the value
     517    for (off_t j = 0; j < MosaicN_Image[i]; j++) {
     518      off_t m = MosaicToImage[i][j];
    687519
    688520      if (!isfinite(image[m].McalPSF)) {
    689         image[m].McalPSF = 0.0;
    690         image[m].McalAPER = 0.0;
    691         image[m].dMcal = 0.0;
     521        image[m].McalPSF   = 0.0;
     522        image[m].McalAPER  = 0.0;
     523        image[m].dMcal     = 0.0;
    692524        image[m].McalChiSq = 0.0;
    693525        fprintf (stderr, "warning: resetting NAN value for Mcal %s\n", image[m].name);
    … …  
    703535      image[m].dMcal     = NAN;
    704536      image[m].McalChiSq = NAN;
    705 
    706     }
    707     dS /= MosaicN_Image[i];
    708 
    709     InitCoords (&mosaic[i].coords, "DEC--TAN");
    710     mosaic[i].coords.crval1 = Rmin;
    711     mosaic[i].coords.crval2 = Dmin;
    712     mosaic[i].coords.cdelt1 = dS;
    713     mosaic[i].coords.cdelt2 = dS;
    714 
    715     RD_to_XY (&dX, &dY, Rmax, Dmax, &mosaic[i].coords);
     537    }
    716538
    717539    mosaic[i].McalPSF   = McalPSF / MosaicN_Image[i];
    … …  
    720542    mosaic[i].McalChiSq = McalChiSq / MosaicN_Image[i];
    721543  }
    722   if (!USE_GRID) return;
    723 
    724   dXmax = MAX (dXmax, dX);
    725   dYmax = MAX (dYmax, dY);
    726   initGrid (dXmax, dYmax);
    727544  return;
    728545}
    … …  
    12311048      continue;
    12321049    }
    1233     float Mgrid = getMgrid (m, c);
     1050    float Mgrid = getMgridTiny (&catalog[c].measureT[m]);
    12341051    if (isnan(Mgrid)) {
    12351052      info->Ngrid ++;
    … …  
    13821199    float Mcal     = getMcal  (m, c, MAG_CLASS_PSF);
    13831200    // float Mgrp     = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
    1384     float Mgrid    = getMgrid (m, c);
     1201    float Mgrid    = getMgridTiny (&catalog[c].measureT[m]);
    13851202    float MrelPSF  = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
    13861203    float Mflat    = getMflat (m, c, flatcorr, catalog);
    … …  
    14811298           summary.Nskip);
    14821299  free (threadinfo);
     1300  // XXX SetMmosInfoFree (&summary);
    14831301
    14841302  npass_output ++;
    … …  
    18111629      Mgrp = getMgrp  (m, c, catalog[c].measureT[m].airmass, NULL);
    18121630      if (isnan(Mgrp)) continue;
    1813       Mgrid = getMgrid (m, c);
     1631      Mgrid = getMgridTiny (&catalog[c].measureT[m]);
    18141632      if (isnan(Mgrid)) continue;
    18151633      Mrel = getMrel  (catalog, m, c, MAG_CLASS_PSF, MAG_SRC_CHP);
Note: See TracChangeset for help on using the changeset viewer.