IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36635


Ignore:
Timestamp:
Apr 2, 2014, 1:13:06 PM (12 years ago)
Author:
watersc1
Message:

This version appears to complete and generate the expected output products. Now to test and tune the parameters.

Location:
trunk/ppBackground/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppBackground/src/ppBackgroundStack.h

    r36615 r36635  
    2323
    2424  psImageMap *modelMap;
     25  // These are the full extent of the input data
    2526  psF32 ra_min;
    2627  psF32 ra_max;
    2728  psF32 dec_min;
    2829  psF32 dec_max;
     30
     31  // Because the binning code can't handle subsections sanely, these hold the working subsection.
     32  psF32 x_min;
     33  psF32 x_max;
     34  psF32 y_min;
     35  psF32 y_max;
    2936
    3037  psArray *stack_data;

  • trunk/ppBackground/src/ppBackgroundStackCamera.c

    r36615 r36635  
    3333    pmConfig *config = data->config;    // Because I'm reusing code.
    3434    int u,v;
    35     size_t A,P;
     35    //    size_t A,P;
    3636    double RR = -9999.0 ,DD = -99999.0,rr = 9999.0,dd = 99999.0;
    3737    // FIX Figure out what stacks we have to deal with.
    … …  
    115115      }
    116116
    117       psMemStats(0,&A,&P);
    118       fprintf(stderr,"fpa %ld %ld\n",A,P);
     117/*       psMemStats(0,&A,&P); */
     118/*       fprintf(stderr,"fpa %ld %ld\n",A,P); */
    119119
    120120      pmChip *chip;                       // Chip from FPA
    … …  
    132132        }
    133133
    134         psMemStats(0,&A,&P);
    135         fprintf(stderr,"chip %ld %ld\n",A,P);
     134/*      psMemStats(0,&A,&P); */
     135/*      fprintf(stderr,"chip %ld %ld\n",A,P); */
    136136
    137137        // read WCS data from the corresponding header
    … …  
    279279
    280280       
    281         psMemStats(0,&A,&P);
    282         fprintf(stderr,"chip2 %ld %ld\n",A,P);
    283        
    284 /*      if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { */
    285 /*        psError(psErrorCodeLast(), false, "Error saving data to files."); */
    286 /*        return false; */
    287 /*      } */
     281/*      psMemStats(0,&A,&P); */
     282/*      fprintf(stderr,"chip2 %ld %ld\n",A,P); */
    288283      } // end chip loop
    289284
    … …  
    306301      pmFPAfileActivate(config->files, false, NULL);
    307302      psFree(modelContent);
    308       //      psFree(smfFile);
     303      //       psFree(smfFile);       // This gets freed when we trash the config->files structure.  I think. 
    309304      psFree(smfFileName);
    310305      psFree(view);
    … …  
    313308      config->files = psMetadataAlloc();
    314309      //
    315      
     310
    316311      //
    317312    } // end smf loop
    … …  
    321316    config->files = config_files;
    322317   
    323     // Allocate the modelMap for the region we're covering.
    324 
    325     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    326     psImageBinning *binning = psImageBinningAlloc();
    327     binning->nXruff = 33; // Number of samples
    328     binning->nYruff = 33;
    329     binning->nXfine = ceil(data->ra_max - data->ra_min) + 1; // This is the range we're looking at
    330     binning->nYfine = ceil(data->dec_max - data->dec_min) + 1;
    331     binning->nXskip = floor(data->ra_min);
    332     binning->nYskip = floor(data->dec_min);
    333     psImageBinningSetScale(binning,PS_IMAGE_BINNING_CENTER);
    334     printf("Sizes: sky: %f %f -> %f %f :: tp: %f %f -> %f %f :: map: %d %d\n",
    335            rr,dd,RR,DD,data->ra_min,data->dec_min,data->ra_max,data->dec_max,binning->nXfine,binning->nYfine);
    336     data->modelMap = psImageMapNoImageAlloc( binning,stats);
     318    // remove the {ra|dec}_min values from everything so we don't confuse the binning code.
     319/*     psMetadataIterator *expIter = psMetadataIteratorAlloc(data->models, PS_LIST_HEAD, NULL); */
     320/*     psMetadataItem *expItem; */
     321/*     while ((expItem = psMetadataGetAndIncrement(expIter))) { */
     322/*       if (expItem->type != PS_DATA_METADATA) { */
     323/*      continue; // This is the N counter */
     324/*       } */
     325/*       psMetadataIterator *chipIter = psMetadataIteratorAlloc(expItem->data.md, PS_LIST_HEAD, NULL); */
     326/*       psMetadataItem *chipItem; */
     327/*       while ((chipItem = psMetadataGetAndIncrement(chipIter))) { */
     328/*      psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra"); */
     329/*      psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec"); */
     330/*      for (v = 0; v < ra->numRows; v++) { */
     331/*        for (u = 0; u < ra->numCols; u++) { */
     332/*          ra->data.F32[v][u] -= data->ra_min; */
     333/*          dec->data.F32[v][u] -= data->dec_min; */
     334/*        } */
     335/*      } */
     336/*       } // End chip */
     337/*       psFree(chipIter); */
     338/*     } // End smfs */
     339/*     psFree(expIter); */
     340/*     // And from the structure objects. */
     341/*     data->x_min -= data->ra_min; */
     342/*     data->y_min -= data->dec_min; */
     343/*     data->x_max -= data->ra_min; */
     344/*     data->y_max -= data->dec_min; */
     345/* /\*     data->ra_max -= data->ra_min; *\/ */
     346/* /\*     data->dec_max -= data->dec_min; *\/ */
     347/* /\*     data->ra_min = 0.0; *\/ */
     348/* /\*     data->dec_min = 0.0; *\/ */
     349
    337350   
    338351    return true;

  • trunk/ppBackground/src/ppBackgroundStackLoop.c

    r36615 r36635  
    6868       
    6969      }
    70      
    71      
     70    } // End initialization block.     
     71
     72
     73    // Loop over the input images, and apply the models to construct the restored versions.
     74    for (i = 0; i < data->stack_data->n; i++) {
     75      pmFPAfile *stack = data->stack_data->data[i];
     76      pmFPAview *view = pmFPAviewAlloc(0);
     77
     78      // PART 1:
     79      // Determine the extent of the model map for this stack
     80      data->x_min = 99e99; data->x_max = -99e99;
     81      data->y_min = 99e99; data->y_max = -99e99;
     82      // Allocate the modelMap for the region we're covering.
     83      psPlane *pix = psPlaneAlloc();   // Pixel coordinates on chip
     84      psPlane *tp = psPlaneAlloc();    // Focal plane coordinates
     85     
     86      pix->x = 0; pix->y = 0;
     87      psPlaneTransformApply(tp, stack->fpa->toTPA, pix);
     88      printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y);
     89      if (tp->x < data->x_min) { data->x_min = tp->x; }
     90      if (tp->x > data->x_max) { data->x_max = tp->x; }
     91      if (tp->y < data->y_min) { data->y_min = tp->y; }
     92      if (tp->y > data->y_max) { data->y_max = tp->y; }
     93
     94      pix->x = 6240; pix->y = 0;
     95      psPlaneTransformApply(tp, stack->fpa->toTPA, pix);
     96      printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y);
     97      if (tp->x < data->x_min) { data->x_min = tp->x; }
     98      if (tp->x > data->x_max) { data->x_max = tp->x; }
     99      if (tp->y < data->y_min) { data->y_min = tp->y; }
     100      if (tp->y > data->y_max) { data->y_max = tp->y; }
     101
     102      pix->x = 6240; pix->y = 6243;
     103      psPlaneTransformApply(tp, stack->fpa->toTPA, pix);
     104      printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y);
     105      if (tp->x < data->x_min) { data->x_min = tp->x; }
     106      if (tp->x > data->x_max) { data->x_max = tp->x; }
     107      if (tp->y < data->y_min) { data->y_min = tp->y; }
     108      if (tp->y > data->y_max) { data->y_max = tp->y; }
     109
     110      pix->x = 0; pix->y = 6243;
     111      psPlaneTransformApply(tp, stack->fpa->toTPA, pix);
     112      printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y);
     113      if (tp->x < data->x_min) { data->x_min = tp->x; }
     114      if (tp->x > data->x_max) { data->x_max = tp->x; }
     115      if (tp->y < data->y_min) { data->y_min = tp->y; }
     116      if (tp->y > data->y_max) { data->y_max = tp->y; }
     117
     118/*       data->x_min -= data->ra_min; */
     119/*       data->x_max -= data->ra_min; */
     120/*       data->y_min -= data->dec_min; */
     121/*       data->y_max -= data->dec_min; */
     122     
     123     
     124      psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
     125      psImageBinning *binning = psImageBinningAlloc();
     126      binning->nXruff = 15; // Number of samples
     127      binning->nYruff = 15;
     128      //    binning->nXfine = ceil(data->ra_max - data->ra_min) + 1; // This is the range we're looking at
     129      //    binning->nYfine = ceil(data->dec_max - data->dec_min) + 1;
     130      binning->nXfine = ceil(data->x_max - data->x_min) + 1;
     131      binning->nYfine = ceil(data->y_max - data->y_min) + 1;
     132      binning->nXskip = floor(data->x_min - data->ra_min) + 1;
     133      binning->nYskip = floor(data->y_min - data->dec_min) + 1;
     134      psImageBinningSetScale(binning,PS_IMAGE_BINNING_CENTER);
     135      printf("Sizes: sky: %f %f -> %f %f :: tp: %f %f -> %f %f :: map: %d %d\n",
     136             0.0,0.0,0.0,0.0,data->ra_min,data->dec_min,data->ra_max,data->dec_max,binning->nXfine,binning->nYfine);
     137      printf("Sizes: corners: %f %f -> %f %f map: %d %d\n",
     138             data->x_min,data->y_min,data->x_max,data->y_max,
     139             binning->nXfine,binning->nYfine);
     140      psImage *sizeImage = psImageAlloc(binning->nXfine,binning->nYfine,PS_TYPE_F32);
     141      P_PSIMAGE_SET_COL0(sizeImage, data->x_min);
     142      P_PSIMAGE_SET_ROW0(sizeImage, data->y_min);
     143      data->modelMap = psImageMapAlloc(sizeImage,binning,stats);
     144/*       psFree(sizeImage); */
     145/*       data->modelMap = psImageMapNoImageAlloc( binning,stats); */
     146/*       P_PSIMAGE_SET_COL0(data->modelMap->map, (data->x_min - binning->nXskip) / binning->nXbin); */
     147/*       P_PSIMAGE_SET_ROW0(data->modelMap->map, (data->y_min - binning->nYskip) / binning->nYbin); */
     148
     149      // force col0/row0
     150/*       data->modelMap->map->col0 = data->modelMap->binning->nXskip; */
     151/*       data->modelMap->map->row0 = data->modelMap->binning->nYskip; */
     152     
     153      // PART 2:
     154      // Solve the data into a model for this region of the sky.
     155
    72156      // This seems wrong, but I need a blank modelMap object, so we fit the zero-data we've stored in the calib objects
    73157      printf("Determining blank modelMap!\n");
    … …  
    83167        return(false);
    84168      }
    85     } // End initialization block.
    86 
    87     // This is where a loop would likely start.
    88     {
    89       // Construct the offset information
    90       printf("Model fit!\n");
    91       if (!ppBackgroundStackDataModelFit(data)) {
    92         psError(psErrorCodeLast(), false, "Error determining the exposure/OTA scaling.");
    93         return(false);
    94       }
    95      
    96       // Apply full correction
    97       printf("Calib apply!\n");
    98       if (!ppBackgroundStackCalibApply(data)) {
    99         psError(psErrorCodeLast(), false, "Error applying the calibration models.");
    100         return(false);
    101       }     
    102      
    103       // Fit the new model
    104       printf("Determining model!\n");
    105       if (!ppBackgroundStackModelFit(data)) {
    106         psError(psErrorCodeLast(), false, "Error determining the modelMap object.");
    107         return(false);
    108       }
    109     } // End loop
    110 
    111     // Loop over the input images, and apply the models to construct the restored versions.
    112 
    113     for (i = 0; i < data->stack_data->n; i++) {
    114       pmFPAfile *stack = data->stack_data->data[i];
    115       pmFPAview *view = pmFPAviewAlloc(0);
    116 
     169
     170      // This is where an iterative solution loop would likely start.
     171      {
     172        // Construct the offset information
     173        printf("Model fit!\n");
     174        if (!ppBackgroundStackDataModelFit(data)) {
     175          psError(psErrorCodeLast(), false, "Error determining the exposure/OTA scaling.");
     176          return(false);
     177        }
     178       
     179        // Apply full correction
     180        printf("Calib apply!\n");
     181        if (!ppBackgroundStackCalibApply(data)) {
     182          psError(psErrorCodeLast(), false, "Error applying the calibration models.");
     183          return(false);
     184        }     
     185       
     186        // Fit the new model
     187        printf("Determining model!\n");
     188        if (!ppBackgroundStackModelFit(data)) {
     189          psError(psErrorCodeLast(), false, "Error determining the modelMap object.");
     190          return(false);
     191        }
     192      } // End loop
     193
     194      // PART 3:
     195      // Define output products
    117196      // Define output image.  Why is this always so hard to do?
    118197      pmFPA *tmp_fpa1,*tmp_fpa2;
    119198      tmp_fpa1 = pmFPAConstruct(config->camera,config->cameraName);
    120199      tmp_fpa2 = pmFPAConstruct(config->camera,config->cameraName);
    121 
     200     
    122201      pmFPAfile *stack_model = pmFPAfileDefineOutput(config,tmp_fpa1,"PPBACKGROUND.STACK.MODEL");
    123202     
    … …  
    134213      stack_model->save = true;
    135214      stack_corr->save = true;
    136 
     215     
    137216      printf("I'm about to loop over the parts of this stack: %d\n",i);
    138217      // Iterate over the images.
    … …  
    144223        return(false);
    145224      }
    146 
    147225     
    148226      pmChip *chip;
    … …  
    151229          continue;
    152230        }
    153 
     231       
    154232        if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    155233          psError(psErrorCodeLast(), false, "load failure for Chip");
    … …  
    256334      }
    257335      psFree(view);
    258     }     
     336      psFree(data->modelMap);
     337    }
    259338               
    260339   

  • trunk/ppBackground/src/ppBackgroundStackMath.c

    r36615 r36635  
    9494      psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
    9595     
    96       psVector *obs  = psVectorAlloc(image->numRows * image->numCols, PS_TYPE_F32);
    97       psVector *model= psVectorAlloc(image->numRows * image->numCols, PS_TYPE_F32);
     96      psVector *obs  = psVectorAllocEmpty(image->numRows, PS_TYPE_F32);
     97      psVector *model= psVectorAllocEmpty(image->numRows, PS_TYPE_F32);
    9898     
    9999      for (v = 0; v < image->numRows; v++) {
    100100        for (u = 0; u < image->numCols; u++) {
    101           obs->data.F32[v * image->numCols + u] = image->data.F32[v][u];
    102           model->data.F32[v * image->numCols + u] = psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u]);
     101          if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)||
     102              (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; }
     103          psVectorAppend(obs,image->data.F32[v][u]);
     104          psVectorAppend(model, psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u]));
     105          //      obs->data.F32[v * image->numCols + u] = image->data.F32[v][u];
     106          //      model->data.F32[v * image->numCols + u] = psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u]);
    103107        }
    104108      }
    … …  
    140144      psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg calibrated data");
    141145      psImage *camera= psMetadataLookupPtr(NULL, data->OTA_solutions, chipName);
     146      psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
     147      psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
    142148     
    143149      psF32 offset = psMetadataLookupF32(NULL,chipItem->data.md,"bkg offset");
    … …  
    146152      for (v = 0; v < image->numRows; v++) {
    147153        for (u = 0; u < image->numCols; u++) {
     154          if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)||
     155              (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; }
     156
    148157          model->data.F32[v][u] = scale * image->data.F32[v][u] - offset - camera->data.F32[v][u];
    149158        }
    … …  
    164173  int u,v;
    165174
     175  long used;
    166176  psS16 N = psMetadataLookupS16(NULL, data->models, "N");
    167   psVector *X = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32);
    168   psVector *Y = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32);
    169   psVector *Z = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32);
    170   psVector *E = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32);
    171   psVector *mask = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_VECTOR_MASK);
     177  psVector *X = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32);
     178  psVector *Y = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32);
     179  psVector *Z = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32);
     180  psVector *E = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32);
     181  psVector *mask = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_VECTOR_MASK);
    172182  j = 0;
    173183 
    … …  
    187197      psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
    188198      psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
     199      //      psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
    189200      for (v = 0; v < calib->numRows; v++) {
    190201        for (u = 0; u < calib->numCols; u++) {
    191           X->data.F32[j] = ra->data.F32[v][u];
    192           Y->data.F32[j] = dec->data.F32[v][u];
    193           Z->data.F32[j] = calib->data.F32[v][u];
    194           E->data.F32[j] = 1.0;
    195           mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0;
     202          if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)||
     203              (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) {
     204            j++;           
     205            continue; }
     206          psVectorAppend(X,ra->data.F32[v][u]);
     207          psVectorAppend(Y,dec->data.F32[v][u]);
     208          psVectorAppend(Z,calib->data.F32[v][u]);
     209          psVectorAppend(E,1.0);
     210          psVectorAppend(mask,0);
     211/*        X->data.F32[j] = ra->data.F32[v][u]; */
     212/*        Y->data.F32[j] = dec->data.F32[v][u]; */
     213/*        Z->data.F32[j] = calib->data.F32[v][u]; */
     214/*        E->data.F32[j] = 1.0; */
     215/*        mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0; */
     216          used++;
    196217          j++;
     218         
     219/*        printf("DATA %ld %ld %f %f %f %f\n", */
     220/*               j,used, */
     221/*               ra->data.F32[v][u], */
     222/*               dec->data.F32[v][u], */
     223/*               calib->data.F32[v][u], */
     224/*               model->data.F32[v][u]); */
     225                 
    197226        }
    198227      }
    … …  
    200229    psFree(chipIter);
    201230  } // End smfs
    202 
     231  printf("%ld %ld\n", j,used);
    203232  bool fitStatus;
    204   bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, mask, 0, X, Y, Z, E);
     233  bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, mask, 1, X, Y, Z, E);
    205234
    206235  psFree(expIter);
Note: See TracChangeset for help on using the changeset viewer.