IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36649 for trunk


Ignore:
Timestamp:
Apr 3, 2014, 6:11:18 PM (12 years ago)
Author:
watersc1
Message:

I'm reasonably happy with the results it's making now.

Location:
trunk/ppBackground/src
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppBackground/src/ppBackgroundStack.h

    r36635 r36649  
    2323
    2424  psImageMap *modelMap;
     25  psS32 model_iteration;
    2526  // These are the full extent of the input data
    2627  psF32 ra_min;

  • trunk/ppBackground/src/ppBackgroundStackData.c

    r36615 r36649  
    5050
    5151    data->modelMap = NULL;
     52    data->model_iteration = 0;
    5253    data->ra_min   = 1e9;
    5354    data->ra_max   = -1e9;

  • trunk/ppBackground/src/ppBackgroundStackLoop.c

    r36642 r36649  
    99
    1010#include "ppBackgroundStack.h"
     11
     12#define WCS_TOLERANCE 0.001             // Tolerance for WCS
     13
    1114
    1215bool ppBackgroundStackLoop(ppBackgroundStackData *data // Run-time data
    … …  
    120123      if (tp->y > data->y_max) { data->y_max = tp->y; }
    121124
    122 /*       data->x_min -= data->ra_min; */
    123 /*       data->x_max -= data->ra_min; */
    124 /*       data->y_min -= data->dec_min; */
    125 /*       data->y_max -= data->dec_min; */
    126      
    127      
    128125      psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    129126      psImageBinning *binning = psImageBinningAlloc();
    130127      binning->nXruff = 13; // Number of samples
    131128      binning->nYruff = 13;
    132       //    binning->nXfine = ceil(data->ra_max - data->ra_min) + 1; // This is the range we're looking at
    133       //    binning->nYfine = ceil(data->dec_max - data->dec_min) + 1;
    134129      binning->nXfine = ceil(data->x_max - data->x_min) + 1;
    135130      binning->nYfine = ceil(data->y_max - data->y_min) + 1;
    … …  
    146141      P_PSIMAGE_SET_ROW0(sizeImage, data->y_min);
    147142      data->modelMap = psImageMapAlloc(sizeImage,binning,stats);
    148 /*       psFree(sizeImage); */
    149 /*       data->modelMap = psImageMapNoImageAlloc( binning,stats); */
    150 /*       P_PSIMAGE_SET_COL0(data->modelMap->map, (data->x_min - binning->nXskip) / binning->nXbin); */
    151 /*       P_PSIMAGE_SET_ROW0(data->modelMap->map, (data->y_min - binning->nYskip) / binning->nYbin); */
    152 
    153       // force col0/row0
    154 /*       data->modelMap->map->col0 = data->modelMap->binning->nXskip; */
    155 /*       data->modelMap->map->row0 = data->modelMap->binning->nYskip; */
    156143     
    157144      // PART 2:
    … …  
    173160
    174161      // This is where an iterative solution loop would likely start.
    175       for (int iterator = 0; iterator < 2; iterator++) {
     162      for (int iterator = 0; iterator < 4; iterator++) {
    176163        // Construct the offset information
    177164        printf("Model fit!\n");
    … …  
    310297            psFree(fp);
    311298            psFree(tp);
     299
     300            // Copy WCS (from ppStackUpdateHeader)
     301            pmHDU *inHDU = pmHDUFromCell(readout->parent);
     302            model->parent->hdu = pmHDUAlloc(NULL);
     303            corr->parent->hdu = pmHDUAlloc(NULL);
     304            pmHDU *modHDU= pmHDUFromCell(model->parent);
     305            pmHDU *corHDU= pmHDUFromCell(corr->parent);
     306
     307            if (!modHDU || !inHDU) {
     308              psWarning("Unable to find HDU at FPA level to copy wcs!");
     309            }
     310            else {
     311              if (!pmAstromReadWCS(stack_model->fpa,model_cell->parent,inHDU->header,1.0)) {
     312                psErrorClear();
     313                psWarning("Unable to read WCS astrometry from input FPA!");
     314              }
     315              else {
     316                if (!modHDU->header) {
     317                  modHDU->header = psMetadataAlloc();
     318                }
     319                if (!pmAstromWriteWCS(modHDU->header, stack_model->fpa,model_cell->parent, WCS_TOLERANCE)) {
     320                  psErrorClear();
     321                  psWarning("Unable to read WCS astrometry from input FPA!");
     322                }
     323              }
     324              if (!pmAstromReadWCS(stack_corr->fpa,corr_cell->parent,inHDU->header,1.0)) {
     325                psErrorClear();
     326                psWarning("Unable to read WCS astrometry from input FPA!");
     327              }
     328              else {
     329                if (!corHDU->header) {
     330                  corHDU->header = psMetadataAlloc();
     331                }
     332                if (!pmAstromWriteWCS(corHDU->header, stack_corr->fpa,corr_cell->parent, WCS_TOLERANCE)) {
     333                  psErrorClear();
     334                  psWarning("Unable to read WCS astrometry from input FPA!");
     335                }
     336              }
     337            } // End WCS saving.
     338
     339           
    312340          } // Close readout
    313341          printf("    I'm done with that readout\n");
    … …  
    339367      psFree(view);
    340368      psFree(data->modelMap);
     369      psFree(sizeImage);
    341370    }
    342371               

  • trunk/ppBackground/src/ppBackgroundStackMath.c

    r36642 r36649  
    5555          }
    5656          psMetadata *chipData = psMetadataLookupPtr(NULL, expItem->data.md, workingChip);
     57          if (!chipData) { continue; }
    5758          psImage *image      = psMetadataLookupPtr(NULL, chipData, "bkg image");
     59          if (!image) { continue; }
    5860          psVectorAppend(tmp,image->data.F32[v][u]);
    5961        } // End loop over exposures
    … …  
    6365
    6466        psFree(expIter);
     67        psFree(tmp);
    6568      } // End u
    6669    } // End v
     70    // Remove the median value from this data.  We just want the tilts, not the offsets
     71    psStatsInit(stats);
     72    psImageStats(stats,solution,NULL,0);
     73    for (v = 0; v < solution->numRows; v++) {
     74      for (u = 0; u < solution->numCols; u++) {
     75        solution->data.F32[v][u] -= stats->robustMedian;
     76      }
     77    }
     78   
     79   
    6780  } // End working chip scan
    6881
    … …  
    90103    psMetadataItem *chipItem;
    91104    while ((chipItem = psMetadataGetAndIncrement(chipIter))) {
     105      //      const char *chipName = chipItem->name;
     106     
    92107      psImage *image = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
    93108      psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
    94109      psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
    95      
     110      //      psImage *camera= psMetadataLookupPtr(NULL, data->OTA_solutions, chipName);     
    96111      psVector *obs  = psVectorAllocEmpty(image->numRows, PS_TYPE_F32);
    97112      psVector *model= psVectorAllocEmpty(image->numRows, PS_TYPE_F32);
    98      
     113
     114      int j = 0;
     115      int used = 0;
    99116      for (v = 0; v < image->numRows; v++) {
    100117        for (u = 0; u < image->numCols; u++) {
    101118          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]);
     119              (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { j++; continue; }
     120          psVectorAppend(obs,image->data.F32[v][u]);// - camera->data.F32[v][u]);
    104121          psVectorAppend(model, psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u]));
    105         }
    106       }
    107 
    108       psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,1);
    109       int status = psVectorFitPolynomial1D(poly,NULL,0,model,NULL,obs);
    110       if (!status) {
    111         psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]);
    112         psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]);
    113       }
    114       psFree(poly);
     122          j++;
     123          used++;
     124        }
     125      }
     126
     127      if (used > 0) {
     128
     129        psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,1);
     130        int status = psVectorFitPolynomial1D(poly,NULL,0,model,NULL,obs);
     131        printf("in model fit loop: %d %d %d %f %f\n",status,j,used,poly->coeff[0],poly->coeff[1]);
     132        if (status && (poly->coeff[1] != 0.0)) {
     133          psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]);
     134          psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]);
     135        }
     136        psFree(poly);
     137      }
    115138    } // End OTA loop
    116139    psFree(chipIter);
    … …  
    153176              (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; }
    154177
    155           model->data.F32[v][u] = scale * image->data.F32[v][u] - offset - camera->data.F32[v][u];
     178          model->data.F32[v][u] = scale * image->data.F32[v][u] + offset - camera->data.F32[v][u];
    156179        }
    157180      }
    … …  
    168191// This "averaging" is done using the psImageMapClipFit.
    169192bool ppBackgroundStackModelFit(ppBackgroundStackData *data) {
    170   long j;
    171   int u,v;
    172 
    173   long used;
     193  long j = 0;
     194  int u,v;
     195
     196  long used = 0;
    174197  psS16 N = psMetadataLookupS16(NULL, data->models, "N");
    175198  psVector *X = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32);
    … …  
    195218      psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
    196219      psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
    197       //      psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
     220#define DUMP_DATA 0
     221#if DUMP_DATA
     222      psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
     223
     224      psF32 offset = psMetadataLookupF32(NULL,chipItem->data.md,"bkg offset");
     225      psF32 scale  = psMetadataLookupF32(NULL,chipItem->data.md,"bkg scale");
     226
     227#endif
    198228      for (v = 0; v < calib->numRows; v++) {
    199229        for (u = 0; u < calib->numCols; u++) {
    … …  
    214244          used++;
    215245          j++;
    216          
    217 /*        printf("DATA %ld %ld %f %f %f %f\n", */
    218 /*               j,used, */
    219 /*               ra->data.F32[v][u], */
    220 /*               dec->data.F32[v][u], */
    221 /*               calib->data.F32[v][u], */
    222 /*               model->data.F32[v][u]); */
    223                  
     246#if DUMP_DATA
     247          printf("DATA %d %ld %ld %f %f %f %f %f %f\n",
     248                 data->model_iteration,j,used,
     249                 offset,scale,
     250                 ra->data.F32[v][u],
     251                 dec->data.F32[v][u],
     252                 calib->data.F32[v][u],
     253                 model->data.F32[v][u]);
     254#endif
    224255        }
    225256      }
    … …  
    230261  bool fitStatus;
    231262  bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, mask, 1, X, Y, Z, E);
    232 
     263  data->model_iteration++;
    233264  psFree(expIter);
    234265  psFree(X);
Note: See TracChangeset for help on using the changeset viewer.