IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 41615


Ignore:
Timestamp:
May 25, 2021, 5:00:52 PM (5 years ago)
Author:
eugene
Message:

save grid correction into well-defined location; pass to relphot_client if needed

Location:
branches/eam_branches/relphot.20210521
Files:
7 edited

Legend:

Unmodified
Added
Removed

  • branches/eam_branches/relphot.20210521/include/relphot.h

    r41611 r41615  
    6666  PS1_w = 6,
    6767} PS1filters;
     68
     69typedef enum {
     70  GRID_MEAN = 0,
     71  GRID_STDEV = 1,
     72  GRID_NPTS = 2,
     73} GridOutputMode;
    6874
    6975typedef struct {
    … …  
    374380int    SKY_DEPTH;  /** XXX EAM : depth of catalog tables, fix usage */
    375381char  *SYNTH_ZERO_POINTS;
     382char  *GRID_MEANFILE;
    376383
    377384// globals for parallel region operations
    … …  
    574581void          setMflatFromGrid    PROTO((Catalog *catalog));
    575582void          initGridBins        PROTO(());
    576 GridCorrectionType *getGridCorr   PROTO((int *Nlast));
    577 int           GridCorrectionSave  PROTO((char *filename));
     583GridCorrectionType *getGridCorrByCode PROTO((int code));
     584GridCorrectionType *getGridCorrNext   PROTO((int *Nlast));
     585int           GridCorrectionSave  PROTO((void));
     586int           GridCorrectionSaveFile  PROTO((char *filename, GridOutputMode mode));
     587void          GridCorrectionLoad  PROTO((char *filename));
    578588void          initImageBins       PROTO((Catalog *catalog, int Ncatalog, int doImageList));
    579589void          initImagesSubset    PROTO((ImageSubset *input, off_t *line_number, off_t N));

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

    r41610 r41615  
    11# include "relphot.h"
    22
    3 int GridCorrectionSave (char *filename) {
     3int GridCorrectionSave () {
     4
     5  if (!GRID_ZEROPT) return TRUE;
     6
     7  char filename[1024], uniquer[12];
     8  int TIME = time(NULL);
     9  int PID = getpid();
     10
     11  snprintf_nowarn (uniquer, 12, "%05d.%05d", PID, TIME % 100000);
     12
     13  snprintf (filename, 1024, "%s/gridcorr.mean.%s.fits", CATDIR, uniquer);
     14  if (!GridCorrectionSaveFile (filename, GRID_MEAN)) return FALSE;
     15  GRID_MEANFILE = strcreate (filename);  // save in the global so we can pass to e.g., reload_catalogs
     16
     17  snprintf (filename, 1024, "%s/gridcorr.stdev.%s.fits", CATDIR, uniquer);
     18  if (!GridCorrectionSaveFile (filename, GRID_STDEV)) return FALSE;
     19
     20  snprintf (filename, 1024, "%s/gridcorr.npts.%s.fits", CATDIR, uniquer);
     21  if (!GridCorrectionSaveFile (filename, GRID_NPTS)) return FALSE;
     22
     23  return TRUE;
     24}
     25
     26int GridCorrectionSaveFile (char *filename, GridOutputMode mode) {
    427
    528  /* open file for input */
    … …  
    2750
    2851  while (TRUE) {
    29     // getGridCorr() takes a photcode and returns the *next* valid grid correction pointer (start at -1)
    30     GridCorrectionType *GridCorr = getGridCorr (&Nlast);
     52    // getGridCorreNext() takes a photcode and returns the *next* valid grid correction pointer (start at -1)
     53    GridCorrectionType *GridCorr = getGridCorrNext (&Nlast);
    3154    if (!GridCorr) break;
    3255
    3356    int filter = GetPhotcodeEquivCodebyCode (GridCorr->photcode);
     57    char *filtname = GetPhotcodeNamebyCode (filter); // reference, do not free
     58    char *photname = GetPhotcodeNamebyCode (GridCorr->photcode); // reference, do not free
    3459
    3560    Header header;
    … …  
    3762
    3863    gfits_init_header (&header);
    39     header.bitpix   = -32;
    40     header.Naxes    =   2;
     64    header.bitpix   = (mode == GRID_NPTS) ? 32 : -32;
     65    header.Naxes    = 2;
    4166    header.Naxis[0] = GridCorr->Nx;
    4267    header.Naxis[1] = GridCorr->Ny;
    … …  
    4873    // create header & matrix
    4974    if (!gfits_modify (&header, "PHOTCODE", "%d", 1, GridCorr->photcode)) myAbort ("failed to set PHOTCODE");
    50     if (!gfits_modify (&header, "FILTER",   "%d", 1, filter))             myAbort ("failed to set FILTER");
     75    if (!gfits_modify (&header, "FILTER",   "%s", 1, filtname))           myAbort ("failed to set FILTER");
    5176    if (!gfits_modify (&header, "SEASON",   "%d", 1, 0))                  myAbort ("failed to set SEASON"); // XXX hard-wired for now
    5277    if (!gfits_modify (&header, "X_CHIP",   "%d", 1, -1))                 myAbort ("failed to set X_CHIP");
    … …  
    5883
    5984    int Nx = GridCorr->Nx;
    60     float *value = (float *) matrix.buffer;
    6185
     86    float *fvalue = (float *) matrix.buffer;
     87    int   *ivalue = (int   *) matrix.buffer;
     88   
    6289    for (int iy = 0; iy < GridCorr->Ny; iy++) {
    6390      for (int ix = 0; ix < GridCorr->Nx; ix++) {
    64         value[ix + iy*Nx] = GridCorr->Mgrid[ix][iy];
     91        switch (mode) {
     92          case GRID_MEAN:  fvalue[ix + iy*Nx] = GridCorr-> Mgrid[ix][iy]; break;
     93          case GRID_STDEV: fvalue[ix + iy*Nx] = GridCorr->dMgrid[ix][iy]; break;
     94          case GRID_NPTS:  ivalue[ix + iy*Nx] = GridCorr->nMgrid[ix][iy]; break;
     95        }
    6596      }
    6697    }
    6798   
     99    char extname[85];
     100    switch (mode) {
     101      // XXX If we add seasons, add here:
     102      case GRID_MEAN:  snprintf (extname, 85, "%s.%s", photname, "MEAN");  break;
     103      case GRID_STDEV: snprintf (extname, 85, "%s.%s", photname, "STDEV"); break;
     104      case GRID_NPTS:  snprintf (extname, 85, "%s.%s", photname, "NPTS");  break;
     105    }
     106    gfits_modify (&header, "EXTNAME", "%s", 1, extname);
     107
    68108    gfits_primary_to_extended (&header, "IMAGE", "Image Extension");
    69109    gfits_fwrite_header  (f, &header);
    … …  
    78118}
    79119
    80 # if (0)
    81 CamPhotomCorrection *CamPhotomCorrectionLoad (char *filename) {
     120void GridCorrectionLoad (char *filename) {
    82121
    83   int i, ix, iy, season, filter;
     122  int photcode;
    84123
    85   Header *phu = gfits_alloc_header();
     124  Header phu;
     125  Header header;
     126  Matrix matrix;
     127
     128  if (!GRID_ZEROPT) return;
     129
     130  // allocate the basic containing structures:
     131  initGridBins ();
    86132
    87133  /* open file for input */
    … …  
    89135  if (f == (FILE *) NULL) {
    90136    gprint (GP_ERR, "can't open file for read : %s\n", filename);
    91     return FALSE;
     137    return;
    92138  }
    93139
    94140  // read the PHU header
    95   if (!gfits_load_header (f, phu)) return FALSE;
     141  if (!gfits_fread_header (f, &phu)) return;
    96142 
    97   int Nbytes = gfits_data_size (phu);
     143  // move to the first extension
     144  int Nbytes = gfits_data_size (&phu);
    98145  fseeko (f, Nbytes, SEEK_CUR);
    99146
    100   int Nseason, Nfilter, Nx, Ny, dX, dY, NxCCD, NyCCD;
    101   if (!gfits_scan (phu, "NSEASON", "%d", 1, &Nseason)) return FALSE;
    102   if (!gfits_scan (phu, "NFILTER", "%d", 1, &Nfilter)) return FALSE;
    103   if (!gfits_scan (phu, "NX",      "%d", 1, &Nx))      return FALSE;
    104   if (!gfits_scan (phu, "NY",      "%d", 1, &Ny))      return FALSE;
    105   if (!gfits_scan (phu, "DX",      "%d", 1, &dX))      return FALSE;
    106   if (!gfits_scan (phu, "DY",      "%d", 1, &dY))      return FALSE;
    107   if (!gfits_scan (phu, "NX_CHIP", "%d", 1, &NxCCD))   return FALSE;
    108   if (!gfits_scan (phu, "NY_CHIP", "%d", 1, &NyCCD))   return FALSE;
     147  gfits_free_header (&phu);
    109148
    110   CamPhotomCorrection *cam = CamPhotomCorrectionAlloc (Nx, Ny, Nfilter, Nseason);
    111   cam->phu = phu;
     149  while (gfits_fread_header (f, &header)) {
     150    if (!gfits_fread_matrix (f, &matrix, &header))                         myAbort ("failed to read matrix");
    112151
    113   cam->dX      = dX;       // superpixel sampling, x-dir
    114   cam->dY      = dY;       // superpixel sampling, y-dir
     152    if (!gfits_scan (&header, "PHOTCODE", "%d", 1, &photcode))    myAbort ("failed to get PHOTCODE");
     153    GridCorrectionType *GridCorr = getGridCorrByCode (photcode);
     154    GridCorr->photcode = photcode;
    115155
    116   cam->NxCCD   = NxCCD;    // pixels per chip
    117   cam->NyCCD   = NyCCD;    // pixels per chip
     156    if (!gfits_scan (&header, "DX_CHIP",  "%f", 1, &GridCorr->dX)) myAbort ("failed to get DX_CHIP");
     157    if (!gfits_scan (&header, "DY_CHIP",  "%f", 1, &GridCorr->dY)) myAbort ("failed to get DY_CHIP");
     158    if (!gfits_scan (&header, "NX_CHIP",  "%d", 1, &GridCorr->Nx)) myAbort ("failed to get NX_CHIP");
     159    if (!gfits_scan (&header, "NY_CHIP",  "%d", 1, &GridCorr->Ny)) myAbort ("failed to get NY_CHIP");
     160   
     161    GridCorr->dMgrid = NULL;
     162    GridCorr->nMgrid = NULL;
     163    ALLOCATE (GridCorr->Mgrid, float *, GridCorr->Nx);
     164    for (int ix = 0; ix < GridCorr->Nx; ix++) {
     165      ALLOCATE (GridCorr->Mgrid[ix], float, GridCorr->Ny);
     166    }     
    118167
    119   for (i = 0; i < Nseason; i++) {
    120     char field[256];
    121     double value;
     168    int Nx = GridCorr->Nx;
     169
     170    float *fvalue = (float *) matrix.buffer;
    122171   
    123     snprintf (field, 256, "TS_%03d", i);
    124     if (!gfits_scan (phu, field, "%lf", 1, &value)) return FALSE;
    125     cam->tstart[i] = ohana_mjd_to_sec(value);
    126    
    127     snprintf (field, 256, "TE_%03d", i);
    128     if (!gfits_scan (phu, field, "%lf", 1, &value)) return FALSE;
    129     cam->tstop[i] = ohana_mjd_to_sec(value);
     172    for (int iy = 0; iy < GridCorr->Ny; iy++) {
     173      for (int ix = 0; ix < GridCorr->Nx; ix++) {
     174        GridCorr-> Mgrid[ix][iy] = fvalue[ix + iy*Nx];
     175      }
     176    }
     177    gfits_free_header (&header);
     178    gfits_free_matrix (&matrix);
    130179  }
    131 
    132   while (TRUE) {
    133     // load data for this header : if not found, assume we hit the end of the file
    134     Header *header = gfits_alloc_header();
    135     if (!gfits_load_header (f, header)) {
    136       free (header);
    137       break;
    138     }
    139    
    140     if (!gfits_scan (header, "FILTER", "%d", 1, &filter)) myAbort ("failed to find FILTER");
    141     if (!gfits_scan (header, "SEASON", "%d", 1, &season)) myAbort ("failed to find SEASON");
    142     if (!gfits_scan (header, "X_CHIP", "%d", 1, &ix))     myAbort ("failed to find X_CHIP");
    143     if (!gfits_scan (header, "Y_CHIP", "%d", 1, &iy))     myAbort ("failed to find Y_CHIP");
    144     // XXX NOTE: astroflat.20150209.fits had ix and iy backwards in header
    145     // double-check that the new flats are OK
    146 
    147     Matrix *matrix = gfits_alloc_matrix();
    148     if (!gfits_load_matrix (f, matrix, header)) myAbort ("failed to read matrix");
    149 
    150     int index = ix + iy*cam->Nx + filter*cam->Nchips + season*cam->Nflats;
    151     myAssert (index >= 0, "index too small");
    152     myAssert (index < cam->Nvalues, "index too big");
    153 
    154     myAssert (!cam->matrix[index], "entry already assigned?");
    155 
    156     cam->matrix[index] = matrix;
    157     cam->header[index] = header;
    158   }
    159   return cam;
     180  gfits_free_header (&header);
    160181}
    161 # endif
    162 

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

    r41610 r41615  
    6868    for (int ix = 0; ix < GridCorr[code]->Nx; ix++) {
    6969      FREE (GridCorr[code]-> Mgrid[ix]);
    70       FREE (GridCorr[code]->dMgrid[ix]);
    71       FREE (GridCorr[code]->nMgrid[ix]);
     70      if (GridCorr[code]->dMgrid) FREE (GridCorr[code]->dMgrid[ix]);
     71      if (GridCorr[code]->nMgrid) FREE (GridCorr[code]->nMgrid[ix]);
    7272    }     
    7373
    … …  
    8585}
    8686
    87 GridCorrectionType *getGridCorr (int *Nlast) {
     87GridCorrectionType *getGridCorrNext (int *Nlast) {
    8888
    8989  if (GridCorr == NULL) return NULL;
    … …  
    9999    *Nlast = i;
    100100  }
     101  return result;
     102}
     103
     104GridCorrectionType *getGridCorrByCode (int code) {
     105
     106  if (GridCorr == NULL) return NULL;
     107
     108  if (code >= NGridCorr) return NULL;
     109  if (code <          0) return NULL;
     110
     111  ALLOCATE(GridCorr[code], GridCorrectionType, 1)
     112  GridCorrectionType *result = GridCorr[code];
    101113  return result;
    102114}
    … …  
    135147    ALLOCATE (GridCorr[code]->nMgrid,   int *, NX_BIN);
    136148    for (int ix = 0; ix < NX_BIN; ix++) {
    137       ALLOCATE (GridCorr[code]-> Mgrid[ix], float, NX_BIN);
    138       ALLOCATE (GridCorr[code]->dMgrid[ix], float, NX_BIN);
    139       ALLOCATE (GridCorr[code]->nMgrid[ix],   int, NX_BIN);
     149      ALLOCATE (GridCorr[code]-> Mgrid[ix], float, NY_BIN);
     150      ALLOCATE (GridCorr[code]->dMgrid[ix], float, NY_BIN);
     151      ALLOCATE (GridCorr[code]->nMgrid[ix],   int, NY_BIN);
    140152    }     
    141153  }

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

    r41603 r41615  
    382382  // GRID_ZEROPT is not valid for all cases, probably should be its own mode...
    383383  GRID_ZEROPT = FALSE;
     384  GRID_MEANFILE = NULL; // this is set by GridCorrectionSave()
    384385  if ((N = get_argument (argc, argv, "-grid"))) {
    385386    remove_argument (N, &argc, argv);
    … …  
    602603  FREE (REGION_FILE);
    603604  FREE (SYNTH_ZERO_POINTS);
     605  FREE (GRID_MEANFILE);
    604606  FREE (IMAGE_TABLE);
    605607  FREE (CATDIR);
    … …  
    858860  }
    859861
     862  // load the GridCorrection file
     863  GRID_ZEROPT = FALSE;
     864  GRID_MEANFILE = NULL;
     865  if ((N = get_argument (argc, argv, "-grid"))) {
     866    GRID_ZEROPT = TRUE;
     867    remove_argument (N, &argc, argv);
     868    GRID_MEANFILE = strcreate (argv[N]);
     869    remove_argument (N, &argc, argv);
     870  }
     871
    860872  ImagSelect = FALSE;
    861873  if ((N = get_argument (argc, argv, "-instmag"))) {
    … …  
    916928
    917929  FREE (SYNTH_ZERO_POINTS);
     930  FREE (GRID_MEANFILE);
    918931  FREE (BOUNDARY_TREE);
    919932  FREE (UPDATE_CATFORMAT);

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

    r41607 r41615  
    264264    if (BOUNDARY_TREE)     { strextend (&command, "-boundary-tree %s", BOUNDARY_TREE); }
    265265    if (SYNTH_ZERO_POINTS) { strextend (&command, "-synthphot-zpts %s", SYNTH_ZERO_POINTS); }
     266    if (GRID_ZEROPT)       { strextend (&command, "-grid %s", GRID_MEANFILE); }
    266267    if (USE_BASIC_CHECK)   {  strextend (&command, "-basic-image-search"); }
    267268    if (USE_MCAL_PSF_FOR_STACK_APER) { strextend (&command, "-use-mcal-psf-for-stack-aper"); }

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

    r41606 r41615  
    8080      initImagesSubset (image, NULL, Nimage);
    8181
     82      // load grid corrections here
     83      GridCorrectionLoad (GRID_MEANFILE);
     84
    8285      reload_catalogs (skylist, HOST_ID, HOSTDIR);
    8386      freeImages ((char *)image);
    8487      free (image);
     88      freeGridBins ();
    8589      client_logger_message ("updated catalogs\n");
    8690      relphot_client_free (sky, skylist);

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

    r41610 r41615  
    144144    freeTGroupBins (Ncatalog);
    145145
    146     GridCorrectionSave ("test.grid.fits");
     146    GridCorrectionSave ();
     147
    147148    // end of if (NLOOP > 0) block : this loop determines the offsets per chip
    148149  } else {
Note: See TracChangeset for help on using the changeset viewer.