IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Ignore:
Timestamp:
Mar 27, 2014, 4:21:00 PM (12 years ago)
Author:
watersc1
Message:

Working with metadata structures is far easier than pmFPAfiles. First cleanly running-to-completion version. Now to sort out the resolution issues.

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppBackground/src/ppBackgroundStackMath.c

    r36580 r36615  
    2020#include "ppBackgroundStack.h"
    2121
     22// data->models->name->XY__->{image/ra/dec/calibrated/offset/scale}
     23// data->models->counter->XY__->{image/ra/dec/calibrated/offset/scale}
     24
    2225//
    2326// Determine the best OTA-to-OTA solution for this set of exposures/chips.
    … …  
    2629// solution_{OTA}(u,v) = < data_{OTA}(u,v) >_{exposures}
    2730bool ppBackgroundStackModelFitOTASolution(ppBackgroundStackData *data) {
    28   long i;
    29   int u,v;
    30   pmFPAview *view    = pmFPAviewAlloc(0);
     31  int u,v;
    3132  psMetadataIterator *iter = psMetadataIteratorAlloc(data->OTA_solutions, PS_LIST_HEAD, NULL); // Iterate over all chips.
    3233  psMetadataItem *item;
    3334  psStats *stats = psStatsAlloc( PS_STAT_ROBUST_MEDIAN );
     35  psS16 N = psMetadataLookupS16(NULL, data->models, "N");
    3436 
    3537  while ((item = psMetadataGetAndIncrement(iter))) {
    … …  
    4345    for (v = 0; v < solution->numRows; v++) {
    4446      for (u = 0; u < solution->numCols; u++) {
    45         psVector *tmp = psVectorAllocEmpty(data->smfs->n,PS_TYPE_F32); 
    46        
    47         for (i = 0; i < data->smfs->n; i++) {
    48           pmFPAviewReset(view);
    49           pmFPAfile *smfFile = data->smfs->data[i];
    50          
    51           pmChip *chip;                       // Chip from FPA
    52           while ((chip = pmFPAviewNextChip(view, smfFile->fpa, 1))) {
    53             if (!chip->process || !chip->file_exists) {
    54               continue;
    55             }
    56             const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip
    57             if (workingChip && strcmp(chipName, workingChip) != 0) { // Is this chip the one we want?
    58               continue;
    59             }
    60             psImage *image = psMetadataLookupPtr(NULL,chip->concepts, "bkg image");
    61             psVectorAppend(tmp,image->data.F32[v][u]);
    62           } // End chip scan
    63         } // End smf scan
     47
     48        psVector *tmp = psVectorAllocEmpty(N,PS_TYPE_F32);     
     49
     50        psMetadataIterator *expIter = psMetadataIteratorAlloc(data->models, PS_LIST_HEAD, NULL);
     51        psMetadataItem *expItem;
     52        while ((expItem = psMetadataGetAndIncrement(expIter))) {
     53          if (expItem->type != PS_DATA_METADATA) {
     54            continue; // This is the N counter
     55          }
     56          psMetadata *chipData = psMetadataLookupPtr(NULL, expItem->data.md, workingChip);
     57          psImage *image      = psMetadataLookupPtr(NULL, chipData, "bkg image");
     58          psVectorAppend(tmp,image->data.F32[v][u]);
     59        } // End loop over exposures
    6460        psStatsInit(stats);
    6561        psVectorStats(stats,tmp,NULL,NULL,0);
    6662        solution->data.F32[v][u] = stats->robustMedian;
     63
     64        psFree(expIter);
    6765      } // End u
    6866    } // End v
    6967  } // End working chip scan
    7068
     69  psFree(iter);
     70  psFree(stats);
    7171  return(true);
    7272}
    … …  
    7878// FOR: scale, offset
    7979bool ppBackgroundStackDataModelFit(ppBackgroundStackData *data) {
    80   long i;
    81   int u,v;
    82  
    83   pmFPAview *view    = pmFPAviewAlloc(0);
    84   for (i = 0; i < data->smfs->n; i++) {
    85     pmFPAviewReset(view);
    86     pmFPAfile *smfFile = data->smfs->data[i];
     80  int u,v;
     81
     82  psMetadataIterator *expIter = psMetadataIteratorAlloc(data->models, PS_LIST_HEAD, NULL);
     83  psMetadataItem *expItem;
     84  while ((expItem = psMetadataGetAndIncrement(expIter))) {
     85    if (expItem->type != PS_DATA_METADATA) {
     86      continue; // This is the N counter
     87    }
    8788   
    88     pmChip *chip;                       // Chip from FPA
    89     while ((chip = pmFPAviewNextChip(view, smfFile->fpa, 1))) {
    90       if (!chip->process || !chip->file_exists) {
    91         continue;
    92       }
    93       psImage *image = psMetadataLookupPtr(NULL, chip->concepts, "bkg image");
    94       psImage *ra    = psMetadataLookupPtr(NULL, chip->concepts, "bkg ra");
    95       psImage *dec   = psMetadataLookupPtr(NULL, chip->concepts, "bkg dec");
     89    psMetadataIterator *chipIter = psMetadataIteratorAlloc(expItem->data.md, PS_LIST_HEAD, NULL);
     90    psMetadataItem *chipItem;
     91    while ((chipItem = psMetadataGetAndIncrement(chipIter))) {
     92      psImage *image = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
     93      psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
     94      psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
    9695     
    9796      psVector *obs  = psVectorAlloc(image->numRows * image->numCols, PS_TYPE_F32);
    … …  
    108107      int status = psVectorFitPolynomial1D(poly,NULL,0,model,NULL,obs);
    109108      if (!status) {
    110         psMetadataAddF32(chip->analysis,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]);
    111         psMetadataAddF32(chip->analysis,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]);
    112       }
    113       //      psFree(tmp);
     109        psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]);
     110        psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]);
     111      }
     112      psFree(poly);
    114113    } // End OTA loop
     114    psFree(chipIter);
    115115  } // End smf/exp loop
    116   psFree(view);
     116  psFree(expIter);
    117117  return(true);
    118118}
    … …  
    124124// calib_{exposure,OTA}(r,d) = scale * data_{exposure,OTA}(r,d) + offset + solution_{OTA}(u,v)
    125125bool ppBackgroundStackCalibApply(ppBackgroundStackData *data) {
    126   long i;
    127   int u,v;
    128  
    129   pmFPAview *view    = pmFPAviewAlloc(0);
    130   //  psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN);
    131   for (i = 0; i < data->smfs->n; i++) {
    132     pmFPAviewReset(view);
    133     pmFPAfile *smfFile = data->smfs->data[i];
    134    
    135     pmChip *chip;                       // Chip from FPA
    136     while ((chip = pmFPAviewNextChip(view, smfFile->fpa, 1))) {
    137       if (!chip->process || !chip->file_exists) {
    138         continue;
    139       }
    140       const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip
    141      
    142       psImage *image = psMetadataLookupPtr(NULL, chip->concepts, "bkg image");
    143       psImage *model = psMetadataLookupPtr(NULL, chip->concepts, "bkg calibrated data");
     126  int u,v;
     127
     128  psMetadataIterator *expIter = psMetadataIteratorAlloc(data->models, PS_LIST_HEAD, NULL);
     129  psMetadataItem *expItem;
     130  while ((expItem = psMetadataGetAndIncrement(expIter))) {
     131    if (expItem->type != PS_DATA_METADATA) {
     132      continue; // This is the N counter
     133    }
     134    psMetadataIterator *chipIter = psMetadataIteratorAlloc(expItem->data.md, PS_LIST_HEAD, NULL);
     135    psMetadataItem *chipItem;
     136    while ((chipItem = psMetadataGetAndIncrement(chipIter))) {
     137      const char *chipName = chipItem->name;
     138     
     139      psImage *image = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image");
     140      psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg calibrated data");
    144141      psImage *camera= psMetadataLookupPtr(NULL, data->OTA_solutions, chipName);
    145142     
    146       psF32 offset = psMetadataLookupF32(NULL,chip->concepts,"bkg offset");
    147       psF32 scale  = psMetadataLookupF32(NULL,chip->concepts,"bkg scale");
     143      psF32 offset = psMetadataLookupF32(NULL,chipItem->data.md,"bkg offset");
     144      psF32 scale  = psMetadataLookupF32(NULL,chipItem->data.md,"bkg scale");
    148145
    149146      for (v = 0; v < image->numRows; v++) {
    … …  
    153150      }
    154151    } // End chip
     152    psFree(chipIter);
    155153  } // End smf
    156   psFree(view);
     154  psFree(expIter);
    157155  return(true);
    158156}
    … …  
    163161// This "averaging" is done using the psImageMapClipFit.
    164162bool ppBackgroundStackModelFit(ppBackgroundStackData *data) {
    165   long i,j;
    166   int u,v;
    167 
    168   psVector *X = psVectorAllocEmpty(data->smfs->n,PS_TYPE_F32);
    169   psVector *Y = psVectorAllocEmpty(data->smfs->n,PS_TYPE_F32);
    170   psVector *Z = psVectorAllocEmpty(data->smfs->n,PS_TYPE_F32);
    171   psVector *E = psVectorAllocEmpty(data->smfs->n,PS_TYPE_F32);
     163  long j;
     164  int u,v;
     165
     166  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);
    172172  j = 0;
    173173 
    174174  pmFPAview *view    = pmFPAviewAlloc(0);
    175175  psStats *stats     = psStatsAlloc( PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV );
    176   for (i = 0; i < data->smfs->n; i++) {
    177     pmFPAviewReset(view);
    178     pmFPAfile *smfFile = data->smfs->data[i];
    179    
    180     pmChip *chip;                       // Chip from FPA
    181     while ((chip = pmFPAviewNextChip(view, smfFile->fpa, 1))) {
    182       if (!chip->process || !chip->file_exists) {
    183         continue;
    184       }
    185       //      const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip
    186       psImage *calib = psMetadataLookupPtr(NULL, chip->concepts, "bkg calibrated data");
    187       psImage *ra    = psMetadataLookupPtr(NULL, chip->concepts, "bkg ra");
    188       psImage *dec    = psMetadataLookupPtr(NULL, chip->concepts, "bkg dec");
     176
     177  psMetadataIterator *expIter = psMetadataIteratorAlloc(data->models, PS_LIST_HEAD, NULL);
     178  psMetadataItem *expItem;
     179  while ((expItem = psMetadataGetAndIncrement(expIter))) {
     180    if (expItem->type != PS_DATA_METADATA) {
     181      continue; // This is the N counter
     182    }
     183    psMetadataIterator *chipIter = psMetadataIteratorAlloc(expItem->data.md, PS_LIST_HEAD, NULL);
     184    psMetadataItem *chipItem;
     185    while ((chipItem = psMetadataGetAndIncrement(chipIter))) {
     186      psImage *calib = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg calibrated data");
     187      psImage *ra    = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra");
     188      psImage *dec   = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec");
    189189      for (v = 0; v < calib->numRows; v++) {
    190190        for (u = 0; u < calib->numCols; u++) {
    … …  
    193193          Z->data.F32[j] = calib->data.F32[v][u];
    194194          E->data.F32[j] = 1.0;
     195          mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0;
     196          j++;
    195197        }
    196198      }
    197199    } // End chip
     200    psFree(chipIter);
    198201  } // End smfs
    199202
    200203  bool fitStatus;
    201   bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, NULL, 0, X, Y, Z, E);
    202  
     204  bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, mask, 0, X, Y, Z, E);
     205
     206  psFree(expIter);
    203207  psFree(X);
    204208  psFree(Y);
    205209  psFree(Z);
    206210  psFree(E);
     211  psFree(mask);
    207212  psFree(stats);
    208213  psFree(view);
Note: See TracChangeset for help on using the changeset viewer.