Index: trunk/ppBackground/src/ppBackgroundStack.h =================================================================== --- trunk/ppBackground/src/ppBackgroundStack.h (revision 36614) +++ trunk/ppBackground/src/ppBackgroundStack.h (revision 36615) @@ -13,4 +13,8 @@ psMetadata *contents; // Metadata containing information about the data arrays psArray *smfs; // List of pointers to individual smf headers + + // data->models->name->XY__->{image/ra/dec/calibrated/offset/scale} + // data->models->counter->XY__->{image/ra/dec/calibrated/offset/scale} + psMetadata *models; // Metadata containing the information about the model data psMetadata *OTA_solutions; // List of pointers to the OTA-specific solution images Index: trunk/ppBackground/src/ppBackgroundStackArguments.c =================================================================== --- trunk/ppBackground/src/ppBackgroundStackArguments.c (revision 36614) +++ trunk/ppBackground/src/ppBackgroundStackArguments.c (revision 36615) @@ -50,4 +50,5 @@ psMetadataAddBool(arguments, PS_LIST_TAIL, "-fitOTAs", 0, "", false); psMetadataAddStr(arguments, PS_LIST_TAIL, "-OTApath", 0, "", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-dumpconfig", 0, "", NULL); if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) { Index: trunk/ppBackground/src/ppBackgroundStackCamera.c =================================================================== --- trunk/ppBackground/src/ppBackgroundStackCamera.c (revision 36614) +++ trunk/ppBackground/src/ppBackgroundStackCamera.c (revision 36615) @@ -33,4 +33,6 @@ pmConfig *config = data->config; // Because I'm reusing code. int u,v; + size_t A,P; + double RR = -9999.0 ,DD = -99999.0,rr = 9999.0,dd = 99999.0; // FIX Figure out what stacks we have to deal with. psString stackName = data->stacks->data[0]; @@ -49,4 +51,10 @@ } psArrayAdd(data->stack_data,data->stack_data->n,psMemIncrRefCounter(stack)); + pmFPAfileActivate(config->files, false, NULL); + + // You know what? Let's just fucking lie to config->files. + psMetadata *config_files = config->files; + config->files = NULL; + config->files = psMetadataAlloc(); // Read over the input background models. @@ -62,18 +70,31 @@ return(false); } - + + // Pull out the information for this exposure psMetadata *input = item->data.md; // the input metadata of interest - psString smfFileName = psMetadataLookupStr(NULL, input, "astrom"); psMetadata *modelContent = psMetadataLookupMetadata(NULL, input, "models"); + + // Allocate the model metadata object + psMetadata *Bmodel = psMetadataAlloc(); + // Read the smf file from this item fileArguments("astrom",smfFileName,"",config); pmFPAfile *smfFile = pmFPAfileDefineFromArgs(&status, config, "PSWARP.ASTROM","astrom"); - + + // taking from pswarpLoadAstrometry.c + smfFile->type = PM_FPA_FILE_WCS; // Read the SMF data pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { return NULL; } + if (!pmFPAfileRead(smfFile,view,config)) { + psError(PS_ERR_IO, false, "failed READ at FPA %s", smfFile->name); + psFree(view); + return false; + } + printf("CZW: Item %d\n",i); // find the FPA phu @@ -93,4 +114,8 @@ } } + + psMemStats(0,&A,&P); + fprintf(stderr,"fpa %ld %ld\n",A,P); + pmChip *chip; // Chip from FPA while ((chip = pmFPAviewNextChip(view, smfFile->fpa, 1))) { @@ -106,4 +131,7 @@ psWarning("More than one cell present for chip %d", view->chip); } + + psMemStats(0,&A,&P); + fprintf(stderr,"chip %ld %ld\n",A,P); // read WCS data from the corresponding header @@ -128,5 +156,5 @@ } - // FIX Load model data for this chip + // Load model data for this chip psString modelFileName = pmConfigConvertFilename(psMetadataLookupStr(NULL, modelContent, chipName), config, PM_FPA_MODE_READ, false); @@ -175,18 +203,27 @@ if (tp->y < data->dec_min) { data->dec_min = tp->y; } else if (tp->y > data->dec_max) { data->dec_max = tp->y; } - } - } - + + if (sky->r < rr) { rr = sky->r; } + else if (sky->r > RR) { RR = sky->r; } + if (sky->r < dd) { dd = sky->d; } + else if (sky->d > DD) { DD = sky->d; } + + } + } + + // Allocate the model data for this chip + psMetadata *this_model = psMetadataAlloc(); + // Attach vectors to teh structure of the chip - psMetadataAddImage(chip->analysis,PS_LIST_TAIL, + psMetadataAddImage(this_model,PS_LIST_TAIL, "bkg image", 0, "ota space X vector", image); - psMetadataAddImage(chip->analysis,PS_LIST_TAIL, + psMetadataAddImage(this_model,PS_LIST_TAIL, "bkg ra", 0, "ota space ra vector", raim); - psMetadataAddImage(chip->analysis,PS_LIST_TAIL, + psMetadataAddImage(this_model,PS_LIST_TAIL, "bkg dec", 0, "ota space dec vector", decim); - psMetadataAddImage(chip->analysis,PS_LIST_TAIL, + psMetadataAddImage(this_model,PS_LIST_TAIL, "bkg calibrated data", 0, "ota space corrected data", model); @@ -194,11 +231,15 @@ // Define default background model parameters, using the assumption: // observed = camera + offset + scale * astrophysical - psMetadataAddF32(chip->analysis,PS_LIST_TAIL, + psMetadataAddF32(this_model,PS_LIST_TAIL, "bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", 0.0); - psMetadataAddF32(chip->analysis,PS_LIST_TAIL, + psMetadataAddF32(this_model,PS_LIST_TAIL, "bkg scale", PS_META_REPLACE, "background scale parameter for this exposure/ota pair", 1.0); - // FIX Free model data for this chip + // Add this model to the current model + psMetadataAddMetadata(Bmodel,PS_LIST_TAIL, + chipName, PS_META_REPLACE, + "model data for this exposure/ota", this_model); + // Free model data for this chip psFree(modelFileName); psFree(modelFits); @@ -212,5 +253,6 @@ psFree(decim); psFree(model); - + psFree(this_model); + // Check to see if we've loaded or allocated an OTA solution container for this chip if (!psMetadataLookupPtr(NULL,data->OTA_solutions,chipName)) { // No solution metadata entry exists for this chipName @@ -235,32 +277,61 @@ } } - - // FIX Hopefully this isn't going to destroy everything. - // Chip + + + psMemStats(0,&A,&P); + fprintf(stderr,"chip2 %ld %ld\n",A,P); + /* if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { */ /* psError(psErrorCodeLast(), false, "Error saving data to files."); */ /* return false; */ /* } */ - } - // FPA -/* if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { */ -/* psError(psErrorCodeLast(), false, "Error saving data to files."); */ -/* return false; */ -/* } */ - - // Add this smf to the smf array. - psArrayAdd(data->smfs,data->smfs->n,psMemIncrRefCounter(smfFile)); - psFree(smfFile); - } - + } // end chip loop + + // Add the set of models to the datastructure + psMetadataAddMetadata(data->models, PS_LIST_TAIL, + smfFile->origname, PS_META_REPLACE, + "model data for this exposure", Bmodel); + psMetadataAddS16(data->models,PS_LIST_TAIL, + "N", PS_META_REPLACE, + "counter", i); + // data->models->name->XY__->{image/ra/dec/calibrated/offset/scale} + // data->models->counter->XY__->{image/ra/dec/calibrated/offset/scale} + + + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Error saving data to files."); + return false; + } + pmFPAfileActivate(config->files, false, NULL); + psFree(modelContent); + // psFree(smfFile); + psFree(smfFileName); + psFree(view); + + psFree(config->files); + config->files = psMetadataAlloc(); + // + + // + } // end smf loop + + // Unlie now. + psFree(config->files); + config->files = config_files; + // Allocate the modelMap for the region we're covering. psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); psImageBinning *binning = psImageBinningAlloc(); - binning->nXruff = 100; // Number of samples - binning->nYruff = 100; + binning->nXruff = 33; // Number of samples + binning->nYruff = 33; binning->nXfine = ceil(data->ra_max - data->ra_min) + 1; // This is the range we're looking at binning->nYfine = ceil(data->dec_max - data->dec_min) + 1; - + binning->nXskip = floor(data->ra_min); + binning->nYskip = floor(data->dec_min); + psImageBinningSetScale(binning,PS_IMAGE_BINNING_CENTER); + printf("Sizes: sky: %f %f -> %f %f :: tp: %f %f -> %f %f :: map: %d %d\n", + rr,dd,RR,DD,data->ra_min,data->dec_min,data->ra_max,data->dec_max,binning->nXfine,binning->nYfine); data->modelMap = psImageMapNoImageAlloc( binning,stats); Index: trunk/ppBackground/src/ppBackgroundStackData.c =================================================================== --- trunk/ppBackground/src/ppBackgroundStackData.c (revision 36614) +++ trunk/ppBackground/src/ppBackgroundStackData.c (revision 36615) @@ -13,7 +13,13 @@ ) { - psFree(data->contents); - psFree(data->smfs); + // CZW 2014-03-24 + // The things that are commented out here, and therefore don't get freed are in that state because + // I'd far rather spend time sorting out actual math issues than dealing with data structures. + + // psFree(data->contents); + +// psFree(data->smfs); + psFree(data->models); psFree(data->OTA_solutions); psFree(data->OTApath); @@ -25,5 +31,5 @@ psFree(data->stacks); psFree(data->outRoot); - psFree(data->config); + //psFree(data->config); return; } @@ -38,4 +44,5 @@ data->smfs = psArrayAlloc(0); + data->models = psMetadataAlloc(); data->OTA_solutions = psMetadataAlloc(); data->fit_OTAS = false; Index: trunk/ppBackground/src/ppBackgroundStackLoop.c =================================================================== --- trunk/ppBackground/src/ppBackgroundStackLoop.c (revision 36614) +++ trunk/ppBackground/src/ppBackgroundStackLoop.c (revision 36615) @@ -20,4 +20,5 @@ // If we didn't load the OTA solution from an external source, we need to build one. if (data->fit_OTAS) { + printf("Fitting OTAs!\n"); if (!ppBackgroundStackModelFitOTASolution(data)) { // Currently can't fail. @@ -25,7 +26,50 @@ return(false); } - } + + if (data->OTApath) { + // This will write the solutions out. + psMetadataIterator *iter = psMetadataIteratorAlloc(data->OTA_solutions, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *item; // Item from iteration + int i = 0; + while ((item = psMetadataGetAndIncrement(iter))) { + i++; + if (item->type != PS_DATA_IMAGE) { + psString chipName = item->name; + psImage *image = item->data.V; + psMetadata *header = psMetadataAlloc(); + psString solutionFileName = psStringCopy(data->OTApath); + psStringAppend(&solutionFileName, ".%s.fits",chipName); + psString resolved = pmConfigConvertFilename(solutionFileName,config,true,true); + psFits *solutionFits = psFitsOpen(resolved,"w"); + if (!solutionFits) { + psError(2, false, "Unable to open FITS file %s to write model %s.", resolved, chipName); + psFitsClose(solutionFits); + psFree(resolved); + return(false); + } + if (!psFitsWriteImage(solutionFits, header, image, 0, NULL)) { + psError(2, false, "Unable to write FITS image %s.", resolved); + psFitsClose(solutionFits); + psFree(resolved); + return false; + } + if (!psFitsClose(solutionFits)) { + psError(2, false, "Unable to close FITS image %s.", resolved); + psFree(resolved); + return false; + } + psFree(resolved); + psFree(solutionFileName); + } + } + psFree(item); + psFree(iter); + } + + } + // This seems wrong, but I need a blank modelMap object, so we fit the zero-data we've stored in the calib objects + printf("Determining blank modelMap!\n"); if (!ppBackgroundStackModelFit(data)) { psError(psErrorCodeLast(), false, "Error determining the blank modelMap object."); @@ -34,4 +78,5 @@ // Apply OTA solution + printf("Calib apply!\n"); if (!ppBackgroundStackCalibApply(data)) { psError(psErrorCodeLast(), false, "Error applying the calibration models."); @@ -43,4 +88,5 @@ { // Construct the offset information + printf("Model fit!\n"); if (!ppBackgroundStackDataModelFit(data)) { psError(psErrorCodeLast(), false, "Error determining the exposure/OTA scaling."); @@ -49,4 +95,5 @@ // Apply full correction + printf("Calib apply!\n"); if (!ppBackgroundStackCalibApply(data)) { psError(psErrorCodeLast(), false, "Error applying the calibration models."); @@ -55,4 +102,5 @@ // Fit the new model + printf("Determining model!\n"); if (!ppBackgroundStackModelFit(data)) { psError(psErrorCodeLast(), false, "Error determining the modelMap object."); @@ -66,7 +114,12 @@ pmFPAfile *stack = data->stack_data->data[i]; pmFPAview *view = pmFPAviewAlloc(0); - // Define output image - pmFPAfile *stack_model = pmFPAfileDefineOutput(config,stack->fpa,"PPBACKGROUND.STACK.MODEL"); - + + // Define output image. Why is this always so hard to do? + pmFPA *tmp_fpa1,*tmp_fpa2; + tmp_fpa1 = pmFPAConstruct(config->camera,config->cameraName); + tmp_fpa2 = pmFPAConstruct(config->camera,config->cameraName); + + pmFPAfile *stack_model = pmFPAfileDefineOutput(config,tmp_fpa1,"PPBACKGROUND.STACK.MODEL"); + if (!stack_model) { psError(psErrorCodeLast(), false, "Unable to generate output model"); @@ -74,5 +127,5 @@ } - pmFPAfile *stack_corr = pmFPAfileDefineOutput(config,stack->fpa,"PPBACKGROUND.STACK.OUTPUT"); + pmFPAfile *stack_corr = pmFPAfileDefineOutput(config,tmp_fpa2,"PPBACKGROUND.STACK.OUTPUT"); if (!stack_corr) { psError(psErrorCodeLast(), false, "Unable to generate output result"); @@ -82,5 +135,15 @@ stack_corr->save = true; + printf("I'm about to loop over the parts of this stack: %d\n",i); // Iterate over the images. + pmFPAfileActivate(config->files,true,"PPBACKGROUND.STACK"); + pmFPAfileActivate(config->files,true,"PPBACKGROUND.STACK.MODEL"); + pmFPAfileActivate(config->files,true,"PPBACKGROUND.STACK.OUTPUT"); + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "load failure for Chip"); + return(false); + } + + pmChip *chip; while ((chip = pmFPAviewNextChip(view, stack->fpa, 1))) { @@ -93,5 +156,5 @@ return(false); } - + printf(" I'm in a chip\n"); pmCell *cell; @@ -109,4 +172,6 @@ continue; } + printf(" I'm in a cell\n"); + // process each of the readouts @@ -120,8 +185,26 @@ continue; } - - pmReadout *model = pmFPAfileThisReadout(config->files,view,"PPBACKGROUND.STACK.MODEL"); - pmReadout *corr = pmFPAfileThisReadout(config->files,view,"PPBACKGROUND.STACK.OUTPUT"); - + printf(" I'm in a readout\n"); + + // Futz with things to get an acceptable output product. + pmCell *model_cell = pmFPAviewThisCell(view,stack_model->fpa); + pmCell *corr_cell = pmFPAviewThisCell(view,stack_corr->fpa); + + pmReadout *model = pmReadoutAlloc(model_cell); + pmReadout *corr = pmReadoutAlloc(corr_cell); + model->data_exists = true; + corr->data_exists = true; + model->parent->data_exists = true; + corr->parent->data_exists = true; + model->parent->parent->data_exists = true; + corr->parent->parent->data_exists = true; + model->image = psImageAlloc(readout->image->numCols,readout->image->numRows,PS_TYPE_F32); + corr->image = psImageAlloc(readout->image->numCols,readout->image->numRows,PS_TYPE_F32); + + model_cell->concepts = psMemIncrRefCounter(cell->concepts); + model_cell->conceptsRead = cell->conceptsRead; + corr_cell->concepts = psMemIncrRefCounter(cell->concepts); + corr_cell->conceptsRead = cell->conceptsRead; + psPlane *pix = psPlaneAlloc(); // Pixel coordinates on chip psPlane *fp = psPlaneAlloc(); // Focal plane coordinates @@ -130,11 +213,14 @@ int x,y; for (y = 0; y < readout->image->numRows; y++) { + pix->y = y; for (x = 0; x < readout->image->numCols; x++) { + pix->x = x; // Calculate model for each pixel of output - psPlaneTransformApply(fp, chip->toFPA, pix); - psPlaneTransformApply(tp, stack->fpa->toTPA, fp); + // psPlaneTransformApply(fp, chip->toFPA, pix); + psPlaneTransformApply(tp, stack->fpa->toTPA, pix); model->image->data.F32[y][x] = psImageMapEval(data->modelMap,tp->x,tp->y); corr->image->data.F32[y][x] = readout->image->data.F32[y][x] + model->image->data.F32[y][x]; + } } @@ -143,19 +229,19 @@ psFree(tp); } // Close readout - + printf(" I'm done with that readout\n"); // Close output image - // Close cells (XXX shouldn't pmFPAfileClose iterate down as needed?) - view->cell = -1; - while ((cell = pmFPAviewNextCell(view, stack->fpa, 1)) != NULL) { - if (!cell->process || !cell->file_exists) { - continue; - } - if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { - psError(psErrorCodeLast(), false, "save failure for Cell"); - return(false); - } - } } // Close Cell - + // Close cells (XXX shouldn't pmFPAfileClose iterate down as needed?) + view->cell = -1; + while ((cell = pmFPAviewNextCell(view, stack->fpa, 1)) != NULL) { + if (!cell->process || !cell->file_exists) { + continue; + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "save failure for Cell"); + return(false); + } + } + // Close chip if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { Index: trunk/ppBackground/src/ppBackgroundStackMath.c =================================================================== --- trunk/ppBackground/src/ppBackgroundStackMath.c (revision 36614) +++ trunk/ppBackground/src/ppBackgroundStackMath.c (revision 36615) @@ -20,4 +20,7 @@ #include "ppBackgroundStack.h" +// data->models->name->XY__->{image/ra/dec/calibrated/offset/scale} +// data->models->counter->XY__->{image/ra/dec/calibrated/offset/scale} + // // Determine the best OTA-to-OTA solution for this set of exposures/chips. @@ -26,10 +29,9 @@ // solution_{OTA}(u,v) = < data_{OTA}(u,v) >_{exposures} bool ppBackgroundStackModelFitOTASolution(ppBackgroundStackData *data) { - long i; - int u,v; - pmFPAview *view = pmFPAviewAlloc(0); + int u,v; psMetadataIterator *iter = psMetadataIteratorAlloc(data->OTA_solutions, PS_LIST_HEAD, NULL); // Iterate over all chips. psMetadataItem *item; psStats *stats = psStatsAlloc( PS_STAT_ROBUST_MEDIAN ); + psS16 N = psMetadataLookupS16(NULL, data->models, "N"); while ((item = psMetadataGetAndIncrement(iter))) { @@ -43,30 +45,28 @@ for (v = 0; v < solution->numRows; v++) { for (u = 0; u < solution->numCols; u++) { - psVector *tmp = psVectorAllocEmpty(data->smfs->n,PS_TYPE_F32); - - for (i = 0; i < data->smfs->n; i++) { - pmFPAviewReset(view); - pmFPAfile *smfFile = data->smfs->data[i]; - - pmChip *chip; // Chip from FPA - while ((chip = pmFPAviewNextChip(view, smfFile->fpa, 1))) { - if (!chip->process || !chip->file_exists) { - continue; - } - const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip - if (workingChip && strcmp(chipName, workingChip) != 0) { // Is this chip the one we want? - continue; - } - psImage *image = psMetadataLookupPtr(NULL,chip->concepts, "bkg image"); - psVectorAppend(tmp,image->data.F32[v][u]); - } // End chip scan - } // End smf scan + + psVector *tmp = psVectorAllocEmpty(N,PS_TYPE_F32); + + psMetadataIterator *expIter = psMetadataIteratorAlloc(data->models, PS_LIST_HEAD, NULL); + psMetadataItem *expItem; + while ((expItem = psMetadataGetAndIncrement(expIter))) { + if (expItem->type != PS_DATA_METADATA) { + continue; // This is the N counter + } + psMetadata *chipData = psMetadataLookupPtr(NULL, expItem->data.md, workingChip); + psImage *image = psMetadataLookupPtr(NULL, chipData, "bkg image"); + psVectorAppend(tmp,image->data.F32[v][u]); + } // End loop over exposures psStatsInit(stats); psVectorStats(stats,tmp,NULL,NULL,0); solution->data.F32[v][u] = stats->robustMedian; + + psFree(expIter); } // End u } // End v } // End working chip scan + psFree(iter); + psFree(stats); return(true); } @@ -78,20 +78,19 @@ // FOR: scale, offset bool ppBackgroundStackDataModelFit(ppBackgroundStackData *data) { - long i; - int u,v; - - pmFPAview *view = pmFPAviewAlloc(0); - for (i = 0; i < data->smfs->n; i++) { - pmFPAviewReset(view); - pmFPAfile *smfFile = data->smfs->data[i]; + int u,v; + + psMetadataIterator *expIter = psMetadataIteratorAlloc(data->models, PS_LIST_HEAD, NULL); + psMetadataItem *expItem; + while ((expItem = psMetadataGetAndIncrement(expIter))) { + if (expItem->type != PS_DATA_METADATA) { + continue; // This is the N counter + } - pmChip *chip; // Chip from FPA - while ((chip = pmFPAviewNextChip(view, smfFile->fpa, 1))) { - if (!chip->process || !chip->file_exists) { - continue; - } - psImage *image = psMetadataLookupPtr(NULL, chip->concepts, "bkg image"); - psImage *ra = psMetadataLookupPtr(NULL, chip->concepts, "bkg ra"); - psImage *dec = psMetadataLookupPtr(NULL, chip->concepts, "bkg dec"); + psMetadataIterator *chipIter = psMetadataIteratorAlloc(expItem->data.md, PS_LIST_HEAD, NULL); + psMetadataItem *chipItem; + while ((chipItem = psMetadataGetAndIncrement(chipIter))) { + psImage *image = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image"); + psImage *ra = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra"); + psImage *dec = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec"); psVector *obs = psVectorAlloc(image->numRows * image->numCols, PS_TYPE_F32); @@ -108,11 +107,12 @@ int status = psVectorFitPolynomial1D(poly,NULL,0,model,NULL,obs); if (!status) { - psMetadataAddF32(chip->analysis,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]); - psMetadataAddF32(chip->analysis,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]); - } - // psFree(tmp); + psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff[0]); + psMetadataAddF32(chipItem->data.md,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff[1]); + } + psFree(poly); } // End OTA loop + psFree(chipIter); } // End smf/exp loop - psFree(view); + psFree(expIter); return(true); } @@ -124,26 +124,23 @@ // calib_{exposure,OTA}(r,d) = scale * data_{exposure,OTA}(r,d) + offset + solution_{OTA}(u,v) bool ppBackgroundStackCalibApply(ppBackgroundStackData *data) { - long i; - int u,v; - - pmFPAview *view = pmFPAviewAlloc(0); - // psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); - for (i = 0; i < data->smfs->n; i++) { - pmFPAviewReset(view); - pmFPAfile *smfFile = data->smfs->data[i]; - - pmChip *chip; // Chip from FPA - while ((chip = pmFPAviewNextChip(view, smfFile->fpa, 1))) { - if (!chip->process || !chip->file_exists) { - continue; - } - const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip - - psImage *image = psMetadataLookupPtr(NULL, chip->concepts, "bkg image"); - psImage *model = psMetadataLookupPtr(NULL, chip->concepts, "bkg calibrated data"); + int u,v; + + psMetadataIterator *expIter = psMetadataIteratorAlloc(data->models, PS_LIST_HEAD, NULL); + psMetadataItem *expItem; + while ((expItem = psMetadataGetAndIncrement(expIter))) { + if (expItem->type != PS_DATA_METADATA) { + continue; // This is the N counter + } + psMetadataIterator *chipIter = psMetadataIteratorAlloc(expItem->data.md, PS_LIST_HEAD, NULL); + psMetadataItem *chipItem; + while ((chipItem = psMetadataGetAndIncrement(chipIter))) { + const char *chipName = chipItem->name; + + psImage *image = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image"); + psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg calibrated data"); psImage *camera= psMetadataLookupPtr(NULL, data->OTA_solutions, chipName); - psF32 offset = psMetadataLookupF32(NULL,chip->concepts,"bkg offset"); - psF32 scale = psMetadataLookupF32(NULL,chip->concepts,"bkg scale"); + psF32 offset = psMetadataLookupF32(NULL,chipItem->data.md,"bkg offset"); + psF32 scale = psMetadataLookupF32(NULL,chipItem->data.md,"bkg scale"); for (v = 0; v < image->numRows; v++) { @@ -153,6 +150,7 @@ } } // End chip + psFree(chipIter); } // End smf - psFree(view); + psFree(expIter); return(true); } @@ -163,28 +161,30 @@ // This "averaging" is done using the psImageMapClipFit. bool ppBackgroundStackModelFit(ppBackgroundStackData *data) { - long i,j; - int u,v; - - psVector *X = psVectorAllocEmpty(data->smfs->n,PS_TYPE_F32); - psVector *Y = psVectorAllocEmpty(data->smfs->n,PS_TYPE_F32); - psVector *Z = psVectorAllocEmpty(data->smfs->n,PS_TYPE_F32); - psVector *E = psVectorAllocEmpty(data->smfs->n,PS_TYPE_F32); + long j; + int u,v; + + psS16 N = psMetadataLookupS16(NULL, data->models, "N"); + psVector *X = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32); + psVector *Y = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32); + psVector *Z = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32); + psVector *E = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_F32); + psVector *mask = psVectorAlloc(N * 60 * 13 * 13,PS_TYPE_VECTOR_MASK); j = 0; pmFPAview *view = pmFPAviewAlloc(0); psStats *stats = psStatsAlloc( PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV ); - for (i = 0; i < data->smfs->n; i++) { - pmFPAviewReset(view); - pmFPAfile *smfFile = data->smfs->data[i]; - - pmChip *chip; // Chip from FPA - while ((chip = pmFPAviewNextChip(view, smfFile->fpa, 1))) { - if (!chip->process || !chip->file_exists) { - continue; - } - // const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip - psImage *calib = psMetadataLookupPtr(NULL, chip->concepts, "bkg calibrated data"); - psImage *ra = psMetadataLookupPtr(NULL, chip->concepts, "bkg ra"); - psImage *dec = psMetadataLookupPtr(NULL, chip->concepts, "bkg dec"); + + psMetadataIterator *expIter = psMetadataIteratorAlloc(data->models, PS_LIST_HEAD, NULL); + psMetadataItem *expItem; + while ((expItem = psMetadataGetAndIncrement(expIter))) { + if (expItem->type != PS_DATA_METADATA) { + continue; // This is the N counter + } + psMetadataIterator *chipIter = psMetadataIteratorAlloc(expItem->data.md, PS_LIST_HEAD, NULL); + psMetadataItem *chipItem; + while ((chipItem = psMetadataGetAndIncrement(chipIter))) { + psImage *calib = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg calibrated data"); + psImage *ra = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra"); + psImage *dec = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec"); for (v = 0; v < calib->numRows; v++) { for (u = 0; u < calib->numCols; u++) { @@ -193,16 +193,21 @@ Z->data.F32[j] = calib->data.F32[v][u]; E->data.F32[j] = 1.0; + mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0; + j++; } } } // End chip + psFree(chipIter); } // End smfs bool fitStatus; - bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, NULL, 0, X, Y, Z, E); - + bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, mask, 0, X, Y, Z, E); + + psFree(expIter); psFree(X); psFree(Y); psFree(Z); psFree(E); + psFree(mask); psFree(stats); psFree(view);