IPP Software Navigation Tools IPP Links Communication Pan-STARRS Links

Changeset 36615 for trunk


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.

Location:
trunk/ppBackground/src
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/ppBackground/src/ppBackgroundStack.h

    r36582 r36615  
    1313  psMetadata *contents;   // Metadata containing information about the data arrays
    1414  psArray *smfs;          // List of pointers to individual smf headers
     15 
     16  // data->models->name->XY__->{image/ra/dec/calibrated/offset/scale}
     17  // data->models->counter->XY__->{image/ra/dec/calibrated/offset/scale}
     18  psMetadata *models;     // Metadata containing the information about the model data
    1519 
    1620  psMetadata *OTA_solutions; // List of pointers to the OTA-specific solution images

  • trunk/ppBackground/src/ppBackgroundStackArguments.c

    r36580 r36615  
    5050    psMetadataAddBool(arguments, PS_LIST_TAIL, "-fitOTAs", 0, "", false);
    5151    psMetadataAddStr(arguments, PS_LIST_TAIL, "-OTApath", 0, "", NULL);
     52    psMetadataAddStr(arguments, PS_LIST_TAIL, "-dumpconfig", 0, "", NULL);
    5253
    5354    if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) {

  • trunk/ppBackground/src/ppBackgroundStackCamera.c

    r36582 r36615  
    3333    pmConfig *config = data->config;    // Because I'm reusing code.
    3434    int u,v;
     35    size_t A,P;
     36    double RR = -9999.0 ,DD = -99999.0,rr = 9999.0,dd = 99999.0;
    3537    // FIX Figure out what stacks we have to deal with.
    3638    psString stackName = data->stacks->data[0];
    … …  
    4951    }
    5052    psArrayAdd(data->stack_data,data->stack_data->n,psMemIncrRefCounter(stack));
     53    pmFPAfileActivate(config->files, false, NULL);
     54
     55    // You know what? Let's just fucking lie to config->files.
     56    psMetadata *config_files = config->files;
     57    config->files = NULL;
     58    config->files = psMetadataAlloc();
    5159   
    5260    // Read over the input background models.   
    … …  
    6270        return(false);
    6371      }
    64 
     72     
     73      // Pull out the information for this exposure
    6574      psMetadata *input = item->data.md; // the input metadata of interest
    66 
    6775      psString smfFileName = psMetadataLookupStr(NULL, input, "astrom");
    6876      psMetadata *modelContent = psMetadataLookupMetadata(NULL, input, "models");
     77
     78      // Allocate the model metadata object
     79      psMetadata *Bmodel = psMetadataAlloc();
     80     
    6981      // Read the smf file from this item
    7082      fileArguments("astrom",smfFileName,"",config);
    7183      pmFPAfile *smfFile = pmFPAfileDefineFromArgs(&status, config, "PSWARP.ASTROM","astrom");
    72 
     84     
     85      // taking from pswarpLoadAstrometry.c
     86      smfFile->type = PM_FPA_FILE_WCS;
    7387      // Read the SMF data
    7488      pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy
     89
    7590      if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
    7691        return NULL;
    7792      }
     93      if (!pmFPAfileRead(smfFile,view,config)) {
     94        psError(PS_ERR_IO, false, "failed READ at FPA %s", smfFile->name);
     95        psFree(view);
     96        return false;
     97      }
     98     
    7899      printf("CZW: Item %d\n",i);
    79100      // find the FPA phu
    … …  
    93114        }
    94115      }
     116
     117      psMemStats(0,&A,&P);
     118      fprintf(stderr,"fpa %ld %ld\n",A,P);
     119
    95120      pmChip *chip;                       // Chip from FPA
    96121      while ((chip = pmFPAviewNextChip(view, smfFile->fpa, 1))) {
    … …  
    106131          psWarning("More than one cell present for chip %d", view->chip);
    107132        }
     133
     134        psMemStats(0,&A,&P);
     135        fprintf(stderr,"chip %ld %ld\n",A,P);
    108136
    109137        // read WCS data from the corresponding header
    … …  
    128156        }
    129157
    130         // FIX Load model data for this chip
     158        // Load model data for this chip
    131159        psString modelFileName = pmConfigConvertFilename(psMetadataLookupStr(NULL, modelContent, chipName),
    132160                                                         config, PM_FPA_MODE_READ, false);
    … …  
    175203            if (tp->y < data->dec_min) { data->dec_min = tp->y; }
    176204            else if (tp->y > data->dec_max) { data->dec_max = tp->y; }
    177           }
    178         }
    179 
     205
     206            if (sky->r < rr) { rr = sky->r; }
     207            else if (sky->r > RR) { RR = sky->r; }
     208            if (sky->r < dd) { dd = sky->d; }
     209            else if (sky->d > DD) { DD = sky->d; }
     210           
     211          }
     212        }
     213
     214        // Allocate the model data for this chip
     215        psMetadata *this_model = psMetadataAlloc();
     216               
    180217        // Attach vectors to teh structure of the chip
    181         psMetadataAddImage(chip->analysis,PS_LIST_TAIL,
     218        psMetadataAddImage(this_model,PS_LIST_TAIL,
    182219                           "bkg image", 0,
    183220                           "ota space X vector", image);
    184         psMetadataAddImage(chip->analysis,PS_LIST_TAIL,
     221        psMetadataAddImage(this_model,PS_LIST_TAIL,
    185222                           "bkg ra", 0,
    186223                           "ota space ra vector", raim);
    187         psMetadataAddImage(chip->analysis,PS_LIST_TAIL,
     224        psMetadataAddImage(this_model,PS_LIST_TAIL,
    188225                           "bkg dec", 0,
    189226                           "ota space dec vector", decim);
    190         psMetadataAddImage(chip->analysis,PS_LIST_TAIL,
     227        psMetadataAddImage(this_model,PS_LIST_TAIL,
    191228                           "bkg calibrated data", 0,
    192229                           "ota space corrected data", model);
    … …  
    194231        // Define default background model parameters, using the assumption:
    195232        // observed = camera + offset + scale * astrophysical
    196         psMetadataAddF32(chip->analysis,PS_LIST_TAIL,
     233        psMetadataAddF32(this_model,PS_LIST_TAIL,
    197234                         "bkg offset", PS_META_REPLACE,
    198235                         "background offset for this exposure/ota pair", 0.0);
    199         psMetadataAddF32(chip->analysis,PS_LIST_TAIL,
     236        psMetadataAddF32(this_model,PS_LIST_TAIL,
    200237                         "bkg scale", PS_META_REPLACE,
    201238                         "background scale parameter for this exposure/ota pair", 1.0);
    202         // FIX Free model data for this chip
     239        // Add this model to the current model
     240        psMetadataAddMetadata(Bmodel,PS_LIST_TAIL,
     241                              chipName, PS_META_REPLACE,
     242                              "model data for this exposure/ota", this_model);
     243        // Free model data for this chip
    203244        psFree(modelFileName);
    204245        psFree(modelFits);
    … …  
    212253        psFree(decim);
    213254        psFree(model);
    214 
     255        psFree(this_model);
     256       
    215257        // Check to see if we've loaded or allocated an OTA solution container for this chip
    216258        if (!psMetadataLookupPtr(NULL,data->OTA_solutions,chipName)) { // No solution metadata entry exists for this chipName
    … …  
    235277          }
    236278        }
    237        
    238         // FIX Hopefully this isn't going to destroy everything.
    239         // Chip
     279
     280       
     281        psMemStats(0,&A,&P);
     282        fprintf(stderr,"chip2 %ld %ld\n",A,P);
     283       
    240284/*      if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { */
    241285/*        psError(psErrorCodeLast(), false, "Error saving data to files."); */
    242286/*        return false; */
    243287/*      } */
    244       }
    245       // FPA
    246 /*       if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { */
    247 /*      psError(psErrorCodeLast(), false, "Error saving data to files."); */
    248 /*      return false; */
    249 /*       } */
    250 
    251       // Add this smf to the smf array.
    252       psArrayAdd(data->smfs,data->smfs->n,psMemIncrRefCounter(smfFile));
    253       psFree(smfFile);
    254     }
    255 
     288      } // end chip loop
     289
     290      // Add the set of models to the datastructure
     291      psMetadataAddMetadata(data->models, PS_LIST_TAIL,
     292                            smfFile->origname, PS_META_REPLACE,
     293                            "model data for this exposure", Bmodel);
     294      psMetadataAddS16(data->models,PS_LIST_TAIL,
     295                       "N", PS_META_REPLACE,
     296                       "counter", i);
     297      // data->models->name->XY__->{image/ra/dec/calibrated/offset/scale}
     298      // data->models->counter->XY__->{image/ra/dec/calibrated/offset/scale}
     299
     300     
     301     
     302      if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     303        psError(psErrorCodeLast(), false, "Error saving data to files.");
     304        return false;
     305      }
     306      pmFPAfileActivate(config->files, false, NULL);
     307      psFree(modelContent);
     308      //      psFree(smfFile);
     309      psFree(smfFileName);
     310      psFree(view);
     311
     312      psFree(config->files);
     313      config->files = psMetadataAlloc();
     314      //
     315     
     316      //
     317    } // end smf loop
     318
     319    // Unlie now.
     320    psFree(config->files);
     321    config->files = config_files;
     322   
    256323    // Allocate the modelMap for the region we're covering.
    257324
    258325    psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
    259326    psImageBinning *binning = psImageBinningAlloc();
    260     binning->nXruff = 100; // Number of samples
    261     binning->nYruff = 100;
     327    binning->nXruff = 33; // Number of samples
     328    binning->nYruff = 33;
    262329    binning->nXfine = ceil(data->ra_max - data->ra_min) + 1; // This is the range we're looking at
    263330    binning->nYfine = ceil(data->dec_max - data->dec_min) + 1;
    264 
     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);
    265336    data->modelMap = psImageMapNoImageAlloc( binning,stats);
    266337   

  • trunk/ppBackground/src/ppBackgroundStackData.c

    r36582 r36615  
    1313    )
    1414{
    15   psFree(data->contents);
    16   psFree(data->smfs);
     15  // CZW 2014-03-24
     16  // The things that are commented out here, and therefore don't get freed are in that state because
     17  // I'd far rather spend time sorting out actual math issues than dealing with data structures.
     18 
     19  //  psFree(data->contents);
     20 
     21//  psFree(data->smfs);
    1722
     23  psFree(data->models);
    1824  psFree(data->OTA_solutions);
    1925  psFree(data->OTApath);
    … …  
    2531  psFree(data->stacks);
    2632  psFree(data->outRoot);
    27   psFree(data->config);
     33  //psFree(data->config);
    2834  return;
    2935}
    … …  
    3844    data->smfs = psArrayAlloc(0);
    3945
     46    data->models = psMetadataAlloc();
    4047    data->OTA_solutions = psMetadataAlloc();
    4148    data->fit_OTAS = false;

  • trunk/ppBackground/src/ppBackgroundStackLoop.c

    r36582 r36615  
    2020      // If we didn't load the OTA solution from an external source, we need to build one.
    2121      if (data->fit_OTAS) {
     22        printf("Fitting OTAs!\n");
    2223        if (!ppBackgroundStackModelFitOTASolution(data)) {
    2324          // Currently can't fail.
    … …  
    2526          return(false);
    2627        }
    27       }
     28
     29        if (data->OTApath) {
     30          // This will write the solutions out.
     31          psMetadataIterator *iter = psMetadataIteratorAlloc(data->OTA_solutions, PS_LIST_HEAD, NULL); // Iterator
     32          psMetadataItem *item; // Item from iteration
     33          int i = 0;
     34          while ((item = psMetadataGetAndIncrement(iter))) {
     35            i++;
     36            if (item->type != PS_DATA_IMAGE) {
     37              psString chipName = item->name;
     38              psImage  *image    = item->data.V;
     39              psMetadata *header = psMetadataAlloc();
     40              psString solutionFileName = psStringCopy(data->OTApath);
     41              psStringAppend(&solutionFileName, ".%s.fits",chipName);
     42              psString resolved = pmConfigConvertFilename(solutionFileName,config,true,true);
     43              psFits *solutionFits = psFitsOpen(resolved,"w");
     44              if (!solutionFits) {
     45                psError(2, false, "Unable to open FITS file %s to write model %s.", resolved, chipName);
     46                psFitsClose(solutionFits);
     47                psFree(resolved);
     48                return(false);
     49              }
     50              if (!psFitsWriteImage(solutionFits, header, image, 0, NULL)) {
     51                psError(2, false, "Unable to write FITS image %s.", resolved);
     52                psFitsClose(solutionFits);
     53                psFree(resolved);
     54                return false;
     55              }
     56              if (!psFitsClose(solutionFits)) {
     57                psError(2, false, "Unable to close FITS image %s.", resolved);
     58                psFree(resolved);
     59                return false;
     60              }
     61              psFree(resolved);
     62              psFree(solutionFileName);
     63            }
     64          }
     65          psFree(item);
     66          psFree(iter);
     67        }
     68       
     69      }
     70     
    2871     
    2972      // This seems wrong, but I need a blank modelMap object, so we fit the zero-data we've stored in the calib objects
     73      printf("Determining blank modelMap!\n");
    3074      if (!ppBackgroundStackModelFit(data)) {
    3175        psError(psErrorCodeLast(), false, "Error determining the blank modelMap object.");
    … …  
    3478     
    3579      // Apply OTA solution
     80      printf("Calib apply!\n");
    3681      if (!ppBackgroundStackCalibApply(data)) {
    3782        psError(psErrorCodeLast(), false, "Error applying the calibration models.");
    … …  
    4388    {
    4489      // Construct the offset information
     90      printf("Model fit!\n");
    4591      if (!ppBackgroundStackDataModelFit(data)) {
    4692        psError(psErrorCodeLast(), false, "Error determining the exposure/OTA scaling.");
    … …  
    4995     
    5096      // Apply full correction
     97      printf("Calib apply!\n");
    5198      if (!ppBackgroundStackCalibApply(data)) {
    5299        psError(psErrorCodeLast(), false, "Error applying the calibration models.");
    … …  
    55102     
    56103      // Fit the new model
     104      printf("Determining model!\n");
    57105      if (!ppBackgroundStackModelFit(data)) {
    58106        psError(psErrorCodeLast(), false, "Error determining the modelMap object.");
    … …  
    66114      pmFPAfile *stack = data->stack_data->data[i];
    67115      pmFPAview *view = pmFPAviewAlloc(0);
    68       // Define output image
    69       pmFPAfile *stack_model = pmFPAfileDefineOutput(config,stack->fpa,"PPBACKGROUND.STACK.MODEL");
    70 
     116
     117      // Define output image.  Why is this always so hard to do?
     118      pmFPA *tmp_fpa1,*tmp_fpa2;
     119      tmp_fpa1 = pmFPAConstruct(config->camera,config->cameraName);
     120      tmp_fpa2 = pmFPAConstruct(config->camera,config->cameraName);
     121
     122      pmFPAfile *stack_model = pmFPAfileDefineOutput(config,tmp_fpa1,"PPBACKGROUND.STACK.MODEL");
     123     
    71124      if (!stack_model) {
    72125        psError(psErrorCodeLast(), false, "Unable to generate output model");
    … …  
    74127      }
    75128     
    76       pmFPAfile *stack_corr  = pmFPAfileDefineOutput(config,stack->fpa,"PPBACKGROUND.STACK.OUTPUT");
     129      pmFPAfile *stack_corr  = pmFPAfileDefineOutput(config,tmp_fpa2,"PPBACKGROUND.STACK.OUTPUT");
    77130      if (!stack_corr) {
    78131        psError(psErrorCodeLast(), false, "Unable to generate output result");
    … …  
    82135      stack_corr->save = true;
    83136
     137      printf("I'm about to loop over the parts of this stack: %d\n",i);
    84138      // Iterate over the images.
     139      pmFPAfileActivate(config->files,true,"PPBACKGROUND.STACK");
     140      pmFPAfileActivate(config->files,true,"PPBACKGROUND.STACK.MODEL");
     141      pmFPAfileActivate(config->files,true,"PPBACKGROUND.STACK.OUTPUT");
     142      if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) {
     143        psError(psErrorCodeLast(), false, "load failure for Chip");
     144        return(false);
     145      }
     146
     147     
    85148      pmChip *chip;
    86149      while ((chip = pmFPAviewNextChip(view, stack->fpa, 1))) {
    … …  
    93156          return(false);
    94157        }
    95 
     158        printf("  I'm in a chip\n");
    96159        pmCell *cell;
    97160       
    … …  
    109172            continue;
    110173          }
     174          printf("    I'm in a cell\n");
     175
    111176
    112177          // process each of the readouts
    … …  
    120185              continue;
    121186            }
    122 
    123             pmReadout *model = pmFPAfileThisReadout(config->files,view,"PPBACKGROUND.STACK.MODEL");
    124             pmReadout *corr  = pmFPAfileThisReadout(config->files,view,"PPBACKGROUND.STACK.OUTPUT");
    125 
     187            printf("      I'm in a readout\n");
     188
     189            // Futz with things to get an acceptable output product.
     190            pmCell *model_cell = pmFPAviewThisCell(view,stack_model->fpa);
     191            pmCell *corr_cell  = pmFPAviewThisCell(view,stack_corr->fpa);
     192           
     193            pmReadout *model = pmReadoutAlloc(model_cell);
     194            pmReadout *corr  = pmReadoutAlloc(corr_cell);
     195            model->data_exists = true;
     196            corr->data_exists = true;
     197            model->parent->data_exists = true;
     198            corr->parent->data_exists = true;
     199            model->parent->parent->data_exists = true;
     200            corr->parent->parent->data_exists = true;
     201            model->image = psImageAlloc(readout->image->numCols,readout->image->numRows,PS_TYPE_F32);
     202            corr->image = psImageAlloc(readout->image->numCols,readout->image->numRows,PS_TYPE_F32);
     203
     204            model_cell->concepts = psMemIncrRefCounter(cell->concepts);
     205            model_cell->conceptsRead = cell->conceptsRead;
     206            corr_cell->concepts = psMemIncrRefCounter(cell->concepts);
     207            corr_cell->conceptsRead = cell->conceptsRead;
     208           
    126209            psPlane *pix = psPlaneAlloc();   // Pixel coordinates on chip
    127210            psPlane *fp = psPlaneAlloc();    // Focal plane coordinates
    … …  
    130213            int x,y;
    131214            for (y = 0; y < readout->image->numRows; y++) {
     215              pix->y = y;
    132216              for (x = 0; x < readout->image->numCols; x++) {
     217                pix->x = x;
    133218                // Calculate model for each pixel of output
    134                 psPlaneTransformApply(fp, chip->toFPA, pix);
    135                 psPlaneTransformApply(tp, stack->fpa->toTPA, fp);
     219                //              psPlaneTransformApply(fp, chip->toFPA, pix);
     220                psPlaneTransformApply(tp, stack->fpa->toTPA, pix);
    136221               
    137222                model->image->data.F32[y][x] = psImageMapEval(data->modelMap,tp->x,tp->y);
    138223                corr->image->data.F32[y][x] = readout->image->data.F32[y][x] + model->image->data.F32[y][x];
     224               
    139225              }
    140226            }
    … …  
    143229            psFree(tp);
    144230          } // Close readout
    145 
     231          printf("    I'm done with that readout\n");
    146232          // Close output image
    147           // Close cells (XXX shouldn't pmFPAfileClose iterate down as needed?)
    148           view->cell = -1;
    149           while ((cell = pmFPAviewNextCell(view, stack->fpa, 1)) != NULL) {
    150             if (!cell->process || !cell->file_exists) {
    151               continue;
    152             }
    153             if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
    154                 psError(psErrorCodeLast(), false, "save failure for Cell");
    155                 return(false);
    156             }
    157           }
    158233        } // Close Cell
    159        
     234        // Close cells (XXX shouldn't pmFPAfileClose iterate down as needed?)
     235        view->cell = -1;
     236        while ((cell = pmFPAviewNextCell(view, stack->fpa, 1)) != NULL) {
     237          if (!cell->process || !cell->file_exists) {
     238            continue;
     239          }
     240          if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {
     241            psError(psErrorCodeLast(), false, "save failure for Cell");
     242            return(false);
     243          }
     244        }
     245
    160246        // Close chip
    161247        if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) {

  • 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.