Index: branches/eam_branches/ipp-ops-20130712/ppBackground/src/Makefile.am =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppBackground/src/Makefile.am (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppBackground/src/Makefile.am (revision 37068) @@ -1,3 +1,3 @@ -bin_PROGRAMS = ppBackground +bin_PROGRAMS = ppBackground ppBackgroundStack if HAVE_SVNVERSION @@ -28,4 +28,7 @@ ppBackground_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PPSTATS_LIBS) $(PSPHOT_LIBS) $(PPBACKGROUND_LIBS) +ppBackgroundStack_CPPFLAGS = $(ppBackground_CPPFLAGS) +ppBackgroundStack_LDFLAGS = $(ppBackground_LDFLAGS) + ppBackground_SOURCES = \ ppBackground.c \ @@ -38,4 +41,16 @@ ppBackgroundErrorCodes.c \ ppBackgroundExit.c + +ppBackgroundStack_SOURCES = \ + ppBackgroundStack.c \ + ppBackgroundStackArguments.c \ + ppBackgroundStackCamera.c \ + ppBackgroundStackData.c \ + ppBackgroundStackLoop.c \ + ppBackgroundStackMath.c \ + ppBackgroundVersion.c \ + ppBackgroundErrorCodes.c \ + ppBackgroundExit.c + noinst_HEADERS = \ Index: branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackground.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackground.c (revision 37064) +++ branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackground.c (revision 37068) @@ -18,17 +18,19 @@ ppBackgroundErrorRegister(); + printf("Initializing data\n"); ppBackgroundData *data = ppBackgroundDataInit(&argc, argv); if (!data) { goto DIE; } - + printf("Reading Arguments\n"); if (!ppBackgroundArguments(data, argc, argv)) { goto DIE; } - + printf("Setting up Camera\n"); if (!ppBackgroundCamera(data)) { goto DIE; } + printf("Looping over data!\n"); if (!ppBackgroundLoop(data)) { goto DIE; Index: branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStack.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStack.c (revision 37068) +++ branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStack.c (revision 37068) @@ -0,0 +1,75 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppBackgroundStack.h" + +int main(int argc, char *argv[]) +{ + ppBackgroundVersionPrint(); + + pmErrorRegister(); + psphotErrorRegister(); + ppBackgroundErrorRegister(); + + printf("Data init\n"); + ppBackgroundStackData *data = ppBackgroundStackDataInit(&argc, argv); + if (!data) { + goto DIE; + } + printf("Read arguments\n"); + if (!ppBackgroundStackArguments(data, argc, argv)) { + goto DIE; + } + + printf("Setup camera\n"); + if (!ppBackgroundStackCamera(data)) { + goto DIE; + } + + printf("Do Loop\n"); + if (!ppBackgroundStackLoop(data)) { + goto DIE; + } + + DIE: + ; // Empty statement to satisy compiler + psExit exitValue = ppBackgroundExitCode(PS_EXIT_SUCCESS); // Exit code + +/* if (data && data->stats && data->statsFile) { */ +/* psString stats = psMetadataConfigFormat(data->stats); // Statistics to output */ +/* if (!stats || strlen(stats) == 0) { */ +/* psError(PPBACKGROUND_ERR_IO, false, "Unable to format statistics file"); */ +/* } else if (fprintf(data->statsFile, "%s", stats) != strlen(stats)) { */ +/* psError(PPBACKGROUND_ERR_IO, true, "Unable to write statistics file"); */ +/* } */ +/* psFree(stats); */ +/* if (fclose(data->statsFile) == EOF) { */ +/* psError(PPBACKGROUND_ERR_IO, true, "Unable to close statistics file"); */ +/* } */ +/* data->statsFile = NULL; */ +/* exitValue = ppBackgroundExitCode(exitValue); */ +/* } */ + + if (data) { + psString dump_file = psMetadataLookupStr(NULL, data->config->arguments, "-dumpconfig"); + if (dump_file) { + if (!pmConfigDump(data->config, dump_file)) { + psError(psErrorCodeLast(), false, "Unable to dump configuration."); + exitValue = ppBackgroundExitCode(exitValue); + } + } + psFree(data); + } + + pmConfigDone(); + psLibFinalize(); + + return ppBackgroundExitCode(exitValue); +} + Index: branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStack.h =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStack.h (revision 37068) +++ branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStack.h (revision 37068) @@ -0,0 +1,80 @@ +#ifndef PP_BACKGROUNDSTACK_H +#define PP_BACKGROUNDSTACK_H + +#include +#include + +#include "ppBackgroundErrorCodes.h" + +#define PPBACKGROUND_RECIPE "PPBACKGROUND_STACK" // Recipe name + +// Data for processing +typedef struct { + 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 + bool fit_OTAS; // Calculate OTA solutions based on these inputs. + psString OTApath; // Location for pre-solved OTA solutions. + + psImageMap *modelMap; + psS32 model_iteration; + // These are the full extent of the input data + psF32 ra_min; + psF32 ra_max; + psF32 dec_min; + psF32 dec_max; + + // Because the binning code can't handle subsections sanely, these hold the working subsection. + psF32 x_min; + psF32 x_max; + psF32 y_min; + psF32 y_max; + + psArray *stack_data; + psArray *stacks; // List of stacks to be corrected. + psString outRoot; // Output root name + pmConfig *config; // Configuration +} ppBackgroundStackData; + +/// Initialise data for processing +ppBackgroundStackData *ppBackgroundStackDataInit(int *argc, char *argv[] // Command-line arguments + ); + +/// Parse command-line arguments +bool ppBackgroundStackArguments(ppBackgroundStackData *data, // Data for processing + int argc, char *argv[] // Command-line arguments + ); + +/// Parse camera configurations +bool ppBackgroundStackCamera(ppBackgroundStackData *data // Data for processing + ); + +/// Loop over input data, processing +bool ppBackgroundStackLoop(ppBackgroundStackData *data // Data for processing + ); + +bool ppBackgroundStackModelFitOTASolution(ppBackgroundStackData *data); +bool ppBackgroundStackDataModelFit(ppBackgroundStackData *data); +bool ppBackgroundStackCalibApply(ppBackgroundStackData *data); +bool ppBackgroundStackModelFit(ppBackgroundStackData *data); + + +/// Determine exit code +psExit ppBackgroundExitCode( + psExit exitValue // Current exit code + ); + +/// Add version information to header +bool ppBackgroundVersionHeader( + psMetadata *header // Header to supplement + ); + +/// Print version information to stdout +void ppBackgroundVersionPrint(void); + +#endif Index: branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackArguments.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackArguments.c (revision 37068) +++ branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackArguments.c (revision 37068) @@ -0,0 +1,91 @@ +/** @file ppBackgroundStackArguments.c + * + * @brief + * + * @ingroup ppBackgroundStack + * + * @author Paul Price + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppBackgroundStack.h" + +/// Print usage information and die +static void usage(const char *program, // Name of the program + psMetadata *arguments, // Command-line arguments + ppBackgroundStackData *data // Run-time data + ) +{ + fprintf(stderr, "\nPan-STARRS background replacement\n\n"); + fprintf(stderr, "Usage: %s OUTPUT_ROOT\n\n", program); + fprintf(stderr, "\n"); + psArgumentHelp(arguments); + psFree(data); + + pmConfigDone(); + psLibFinalize(); + + exit(PS_EXIT_CONFIG_ERROR); +} + + +bool ppBackgroundStackArguments(ppBackgroundStackData *data, int argc, char *argv[]) +{ + assert(data); + assert(data->config); + + psMetadata *arguments = data->config->arguments; + psMetadataAddStr(arguments, PS_LIST_TAIL, "-input", 0, "Filename of input metadata", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-image", 0, "Filename of image (required)", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask", 0, "Filename of mask", NULL); + 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) { + usage(argv[0], arguments, data); + } + + unsigned int numBad = 0; // Number of bad lines + data->contents = psMetadataConfigRead(NULL, &numBad,psMetadataLookupStr(NULL, arguments, "-input"),false); + if (!(data->contents) || numBad > 0) { + psError(PPBACKGROUND_ERR_CONFIG, false, "Unable to cleanly read MDC file with inputs."); + return(false); + } + + // Add any images to generate for to the list + // If we're going to do this for a full projection cell, this probably needs to be a metadata object. + psArrayAdd(data->stacks, data->stacks->n, psMemIncrRefCounter(psMetadataLookupStr(NULL, arguments, "-image"))); + + // Determine what we're doing with the OTA solution + if ((!psMetadataLookupBool(NULL,arguments,"-fitOTAs"))&& + (!psMetadataLookupStr(NULL,arguments,"-OTApath"))) { + psError(PPBACKGROUND_ERR_CONFIG, false, "No OTA solution path (-OTApath) provided, and no request to model this (-fitOTAs)."); + return(false); + } + data->fit_OTAS = psMetadataLookupBool(NULL, arguments, "-fitOTAs"); + data->OTApath = psMetadataLookupStr(NULL, arguments, "-OTApath"); + + + // This is the output base + data->outRoot = psStringCopy(argv[1]); + psMetadataAddStr(arguments, PS_LIST_TAIL, "OUTPUT", 0, "Output root name", data->outRoot); + + psTrace("ppBackgroundStack", 1, "Done reading command-line arguments\n"); + + + PS_ASSERT_STRING_NON_EMPTY(data->outRoot, false); + + return true; +} + + Index: branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackCamera.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackCamera.c (revision 37068) +++ branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackCamera.c (revision 37068) @@ -0,0 +1,345 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppBackgroundStack.h" + +/// Add a single filename to the arguments as an array, so that it can be used with pmFPAfileBindFromArgs, etc +static void fileArguments(const char *file, // The symbolic name for the file + const char *name, // The name of the file + const char *comment, // Description of the file + pmConfig *config // Configuration + ) +{ + psArray *files = psArrayAlloc(1); // Array with file names + files->data[0] = psStringCopy(name); + if (psMetadataLookup(config->arguments, file)) { + psMetadataRemoveKey(config->arguments, file); + } + psMetadataAddArray(config->arguments, PS_LIST_TAIL, file, 0, comment, files); + psFree(files); + return; +} + + +bool ppBackgroundStackCamera(ppBackgroundStackData *data // Run-time data + ) +{ + bool status; // Status of file definition + 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]; + fileArguments("IMAGE", stackName, "Input image", data->config); + psFree(stackName); + pmFPAfile *stack = pmFPAfileDefineFromArgs(&status, data->config, + "PPBACKGROUND.STACK", "IMAGE"); + if (!status || !stack) { + psError(psErrorCodeLast(), false, "Failed to build file from PPBACKGROUND.STACK"); + return false; + } + + if (!pmAstromReadBilevelMosaic(stack->fpa,stack->fpa->hdu->header)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input stack."); + return false; + } + 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. + psMetadataIterator *iter = psMetadataIteratorAlloc(data->contents, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *item; // Item from iteration + int i = 0; + while ((item = psMetadataGetAndIncrement(iter))) { + i++; + if (item->type != PS_DATA_METADATA) { + psError(PPBACKGROUND_ERR_ARGUMENTS, true, + "Component %s of the input metadata is not of type METADATA", item->name); + psFree(iter); + 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 + bool bilevelAstrometry = false; + pmHDU *phu = pmFPAviewThisPHU(view, smfFile->fpa); + if (phu) { + char *ctype = psMetadataLookupStr(NULL, phu->header, "CTYPE1"); + if (ctype) { + bilevelAstrometry = !strcmp (&ctype[4], "-DIS"); + } + } + if (bilevelAstrometry) { + if (!pmAstromReadBilevelMosaic(smfFile->fpa, phu->header)) { + psError(psErrorCodeLast(), false, "Unable to read bilevel mosaic astrometry for input FPA."); + psFree(view); + return false; + } + } + + psF32 exptime = 1.0; + exptime = psMetadataLookupF32(NULL, phu->header, "EXPTIME"); + + 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 (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "Error loading data from files."); + return false; + } + if (chip->cells->n != 1) { + 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 + pmHDU *hdu = pmFPAviewThisHDU (view, smfFile->fpa); + if (bilevelAstrometry) { + if (!pmAstromReadBilevelChip (chip, hdu->header)) { + psWarning("Unable to read bilevel chip astrometry for chip %s.", chipName); + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "Error saving data to files."); + return false; + } + continue; + } + } else { + // we use a default FPA pixel scale of 1.0 + psWarning("Reading WCS astrometry for chip %s.", chipName); + if (!pmAstromReadWCS(smfFile->fpa, chip, hdu->header, 1.0)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA."); + psFree(view); + return false; + } + } + + // Load model data for this chip + psString modelFileName = pmConfigConvertFilename(psMetadataLookupStr(NULL, modelContent, chipName), + config, PM_FPA_MODE_READ, false); + psFits *modelFits = psFitsOpen(modelFileName,"r"); + psImage *image = psFitsReadImage(modelFits,psRegionSet(0,0,0,0),0); + psMetadata *header = psFitsReadHeader(NULL,modelFits); + + // Allocate the data structures for this chip + psImage *raim = psImageAlloc(image->numCols,image->numRows,PS_TYPE_F32); + psImage *decim = psImageAlloc(image->numCols,image->numRows,PS_TYPE_F32); + psImage *model = psImageAlloc(image->numCols,image->numRows,PS_TYPE_F32); + + // Read header and construct original positions + psS32 naxis1 = psMetadataLookupS32(NULL, header, "NAXIS1"); + psS32 naxis2 = psMetadataLookupS32(NULL, header, "NAXIS2"); + psS32 imnaxis1 = psMetadataLookupS32(NULL, header, "IMNAXIS1"); + psS32 imnaxis2 = psMetadataLookupS32(NULL, header, "IMNAXIS2"); + psString ccdsum = psMetadataLookupStr(NULL, header, "CCDSUM"); + psS32 xbin = atoi(strtok(ccdsum," ")); + psS32 ybin = atoi(strtok(NULL, " ")); + psS32 xoffset = (naxis1 * xbin - imnaxis1) / (2 * xbin); + psS32 yoffset = (naxis2 * ybin - imnaxis2) / (2 * ybin); + + psPlane *pix = psPlaneAlloc(); // Pixel coordinates on chip + psPlane *fp = psPlaneAlloc(); // Focal plane coordinates + psPlane *tp = psPlaneAlloc(); // Tangent plane coordinates + psSphere *sky = psSphereAlloc(); // Sky coordinates + + for (v = 0; v < image->numRows; v++) { + pix->y = (v - yoffset) * ybin; + for (u = 0; u < image->numCols; u++) { + pix->x = (u - xoffset) * xbin; + psPlaneTransformApply(fp, chip->toFPA, pix); + psPlaneTransformApply(tp, smfFile->fpa->toTPA, fp); + psDeproject(sky, tp, smfFile->fpa->toSky); + + psProject(tp,sky,stack->fpa->toSky); + + raim->data.F32[v][u] = tp->x; + decim->data.F32[v][u] = tp->y; + model->data.F32[v][u] = 0.0; + image->data.F32[v][u] /= exptime; + // Check the bounds so we'll know how large of an area to model in the map + if (tp->x < data->ra_min) { data->ra_min = tp->x; } + else if (tp->x > data->ra_max) { data->ra_max = tp->x; } + 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(this_model,PS_LIST_TAIL, + "bkg image", 0, + "ota space X vector", image); + psMetadataAddImage(this_model,PS_LIST_TAIL, + "bkg ra", 0, + "ota space ra vector", raim); + psMetadataAddImage(this_model,PS_LIST_TAIL, + "bkg dec", 0, + "ota space dec vector", decim); + psMetadataAddImage(this_model,PS_LIST_TAIL, + "bkg calibrated data", 0, + "ota space corrected data", model); + + // Define default background model parameters, using the assumption: + // observed = camera + offset + scale * astrophysical + psMetadataAddF32(this_model,PS_LIST_TAIL, + "bkg offset", PS_META_REPLACE, + "background offset for this exposure/ota pair", 0.0); + psMetadataAddF32(this_model,PS_LIST_TAIL, + "bkg scale", PS_META_REPLACE, + "background scale parameter for this exposure/ota pair", 1.0); + // 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); + psFree(header); + psFree(pix); + psFree(fp); + psFree(tp); + psFree(sky); + psFree(image); + psFree(raim); + 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 + // FIX this should try to find the imagefile. + if (data->fit_OTAS) { // We are fitting OTAs, so allocate a new image + psImage *solution = psImageAlloc(image->numCols,image->numRows,PS_TYPE_F32); + psMetadataAddPtr(data->OTA_solutions,PS_LIST_TAIL, + chipName, PS_DATA_UNKNOWN | PS_META_REPLACE, + "OTA solution element", solution); + } + else { // We are not fitting OTAs, so read the one that should be saved on OTApath. + psString solutionFileName = psStringCopy(data->OTApath); + psStringAppend(&solutionFileName, ".%s.fits",chipName); + psString resolvedFileName = pmConfigConvertFilename(solutionFileName,config,false,false); + psFits *solutionFits = psFitsOpen(resolvedFileName,"r"); + psImage *solution = psFitsReadImage(solutionFits,psRegionSet(0,0,0,0),0); + psMetadataAddImage(data->OTA_solutions,PS_LIST_TAIL, + chipName, 0, + "OTA solution element", solution); + psFree(solutionFits); + psFree(solutionFileName); + psFree(solution); + } + } + } // 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); + + 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); // This gets freed when we trash the config->files structure. I think. + psFree(smfFileName); + psFree(view); + + psFree(config->files); + config->files = psMetadataAlloc(); + // + + // + } // end smf loop + + // Unlie now. + psFree(config->files); + config->files = config_files; + + // remove the {ra|dec}_min values from everything so we don't confuse the binning code. +/* 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 *ra = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra"); */ +/* psImage *dec = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec"); */ +/* for (v = 0; v < ra->numRows; v++) { */ +/* for (u = 0; u < ra->numCols; u++) { */ +/* ra->data.F32[v][u] -= data->ra_min; */ +/* dec->data.F32[v][u] -= data->dec_min; */ +/* } */ +/* } */ +/* } // End chip */ +/* psFree(chipIter); */ +/* } // End smfs */ +/* psFree(expIter); */ +/* // And from the structure objects. */ +/* data->x_min -= data->ra_min; */ +/* data->y_min -= data->dec_min; */ +/* data->x_max -= data->ra_min; */ +/* data->y_max -= data->dec_min; */ +/* /\* data->ra_max -= data->ra_min; *\/ */ +/* /\* data->dec_max -= data->dec_min; *\/ */ +/* /\* data->ra_min = 0.0; *\/ */ +/* /\* data->dec_min = 0.0; *\/ */ + + + return true; +} Index: branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackData.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackData.c (revision 37068) +++ branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackData.c (revision 37068) @@ -0,0 +1,79 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppBackgroundStack.h" + +// Destructor +static void backgroundDataFree(ppBackgroundStackData *data // Data to free + ) +{ + // 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); + + psFree(data->modelMap); + + + psFree(data->stack_data); + psFree(data->stacks); + psFree(data->outRoot); + //psFree(data->config); + return; +} + + +ppBackgroundStackData *ppBackgroundStackDataAlloc(void) +{ + ppBackgroundStackData *data = psAlloc(sizeof(ppBackgroundStackData)); // Processing data, to return + psMemSetDeallocator(data, (psFreeFunc)backgroundDataFree); + + data->contents = NULL; + data->smfs = psArrayAlloc(0); + + data->models = psMetadataAlloc(); + data->OTA_solutions = psMetadataAlloc(); + data->fit_OTAS = false; + data->OTApath = NULL; + + data->modelMap = NULL; + data->model_iteration = 0; + data->ra_min = 1e9; + data->ra_max = -1e9; + data->dec_min = 1e9; + data->dec_max = -1e9; + + data->stack_data = psArrayAlloc(0); + data->stacks = psArrayAlloc(0); + data->outRoot = NULL; + data->config = NULL; + + return data; +} + + +ppBackgroundStackData *ppBackgroundStackDataInit(int *argc, char **argv) +{ + PS_ASSERT_PTR_NON_NULL(argc, NULL); + PS_ASSERT_PTR_NON_NULL(argv, NULL); + + ppBackgroundStackData *data = ppBackgroundStackDataAlloc(); // Processing data, to return + data->config = pmConfigRead(argc, argv, PPBACKGROUND_RECIPE); + if (!data->config) { + psError(psErrorCodeLast(), false, "Unable to read configuration."); + return NULL; + } + return data; +} Index: branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackLoop.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackLoop.c (revision 37068) +++ branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackLoop.c (revision 37068) @@ -0,0 +1,377 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppBackgroundStack.h" + +#define WCS_TOLERANCE 0.001 // Tolerance for WCS + + +bool ppBackgroundStackLoop(ppBackgroundStackData *data // Run-time data + ) +{ + pmConfig *config = data->config; // Configuration data + int i; + // + // Solve the data into a consistent model. + { // This block does the initialization + // 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. + psError(psErrorCodeLast(), false, "Error calculating the OTA solution"); + 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); + } + + } + } // End initialization block. + + + // Loop over the input images, and apply the models to construct the restored versions. + for (i = 0; i < data->stack_data->n; i++) { + pmFPAfile *stack = data->stack_data->data[i]; + pmFPAview *view = pmFPAviewAlloc(0); + + // pmHDU *phu = pmFPAviewThisPHU(view, stack->fpa); + psF32 exptime = 1.0; + exptime = psMetadataLookupF32(NULL, stack->fpa->hdu->header, "EXPTIME"); + + // PART 1: + // Determine the extent of the model map for this stack + data->x_min = 99e99; data->x_max = -99e99; + data->y_min = 99e99; data->y_max = -99e99; + // Allocate the modelMap for the region we're covering. + psPlane *pix = psPlaneAlloc(); // Pixel coordinates on chip + psPlane *tp = psPlaneAlloc(); // Focal plane coordinates + + pix->x = 0; pix->y = 0; + psPlaneTransformApply(tp, stack->fpa->toTPA, pix); + printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y); + if (tp->x < data->x_min) { data->x_min = tp->x; } + if (tp->x > data->x_max) { data->x_max = tp->x; } + if (tp->y < data->y_min) { data->y_min = tp->y; } + if (tp->y > data->y_max) { data->y_max = tp->y; } + + pix->x = 6240; pix->y = 0; + psPlaneTransformApply(tp, stack->fpa->toTPA, pix); + printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y); + if (tp->x < data->x_min) { data->x_min = tp->x; } + if (tp->x > data->x_max) { data->x_max = tp->x; } + if (tp->y < data->y_min) { data->y_min = tp->y; } + if (tp->y > data->y_max) { data->y_max = tp->y; } + + pix->x = 6240; pix->y = 6243; + psPlaneTransformApply(tp, stack->fpa->toTPA, pix); + printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y); + if (tp->x < data->x_min) { data->x_min = tp->x; } + if (tp->x > data->x_max) { data->x_max = tp->x; } + if (tp->y < data->y_min) { data->y_min = tp->y; } + if (tp->y > data->y_max) { data->y_max = tp->y; } + + pix->x = 0; pix->y = 6243; + psPlaneTransformApply(tp, stack->fpa->toTPA, pix); + printf("%f %f -> %f %f\n",pix->x,pix->y,tp->x,tp->y); + if (tp->x < data->x_min) { data->x_min = tp->x; } + if (tp->x > data->x_max) { data->x_max = tp->x; } + if (tp->y < data->y_min) { data->y_min = tp->y; } + if (tp->y > data->y_max) { data->y_max = tp->y; } + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); + psImageBinning *binning = psImageBinningAlloc(); + binning->nXruff = 13; // Number of samples + binning->nYruff = 13; + binning->nXfine = ceil(data->x_max - data->x_min) + 1; + binning->nYfine = ceil(data->y_max - data->y_min) + 1; + binning->nXskip = floor(data->x_min - data->ra_min) + 1; + binning->nYskip = floor(data->y_min - data->dec_min) + 1; + psImageBinningSetScale(binning,PS_IMAGE_BINNING_CENTER); + printf("Sizes: sky: %f %f -> %f %f :: tp: %f %f -> %f %f :: map: %d %d\n", + 0.0,0.0,0.0,0.0,data->ra_min,data->dec_min,data->ra_max,data->dec_max,binning->nXfine,binning->nYfine); + printf("Sizes: corners: %f %f -> %f %f map: %d %d\n", + data->x_min,data->y_min,data->x_max,data->y_max, + binning->nXfine,binning->nYfine); + psImage *sizeImage = psImageAlloc(binning->nXfine,binning->nYfine,PS_TYPE_F32); + P_PSIMAGE_SET_COL0(sizeImage, data->x_min); + P_PSIMAGE_SET_ROW0(sizeImage, data->y_min); + data->modelMap = psImageMapAlloc(sizeImage,binning,stats); + + // PART 2: + // Solve the data into a model for this region of the sky. + + // 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."); + return(false); + } + + // Apply OTA solution + printf("Calib apply!\n"); + if (!ppBackgroundStackCalibApply(data)) { + psError(psErrorCodeLast(), false, "Error applying the calibration models."); + return(false); + } + + // This is where an iterative solution loop would likely start. + for (int iterator = 0; iterator < 4; iterator++) { + // Construct the offset information + printf("Model fit!\n"); + if (!ppBackgroundStackDataModelFit(data)) { + psError(psErrorCodeLast(), false, "Error determining the exposure/OTA scaling."); + return(false); + } + + // Apply full correction + printf("Calib apply!\n"); + if (!ppBackgroundStackCalibApply(data)) { + psError(psErrorCodeLast(), false, "Error applying the calibration models."); + return(false); + } + + // Fit the new model + printf("Determining model!\n"); + if (!ppBackgroundStackModelFit(data)) { + psError(psErrorCodeLast(), false, "Error determining the modelMap object."); + return(false); + } + } // End loop + + // PART 3: + // Define output products + // 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"); + return (false); + } + + pmFPAfile *stack_corr = pmFPAfileDefineOutput(config,tmp_fpa2,"PPBACKGROUND.STACK.OUTPUT"); + if (!stack_corr) { + psError(psErrorCodeLast(), false, "Unable to generate output result"); + return (false); + } + stack_model->save = true; + 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))) { + if (!chip->process || !chip->file_exists) { + continue; + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "load failure for Chip"); + return(false); + } + printf(" I'm in a chip\n"); + pmCell *cell; + + while ((cell = pmFPAviewNextCell(view, stack->fpa, 1)) != NULL) { + psLogMsg ("ppImageLoop", 5, "Cell %d: %x %x\n", view->cell, cell->file_exists, cell->process); + if (!cell->process || !cell->file_exists) { + continue; + } + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "load failure for Cell"); + return(false); + } + if (cell->readouts->n > 1) { + psWarning ("Skipping Video Cell for ppImageDetrendReadout"); + continue; + } + printf(" I'm in a cell\n"); + + + // process each of the readouts + pmReadout *readout; // Readout from cell + while ((readout = pmFPAviewNextReadout (view, stack->fpa, 1)) != NULL) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "load failure for Readout"); + return(false); + } + if (!readout->data_exists) { + continue; + } + 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 + psPlane *tp = psPlaneAlloc(); // Tangent plane coordinates + + 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, pix); + + model->image->data.F32[y][x] = exptime * 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]; + + } + } + psFree(pix); + psFree(fp); + psFree(tp); + + // Copy WCS (from ppStackUpdateHeader) + pmHDU *inHDU = pmHDUFromCell(readout->parent); + model->parent->hdu = pmHDUAlloc(NULL); + corr->parent->hdu = pmHDUAlloc(NULL); + pmHDU *modHDU= pmHDUFromCell(model->parent); + pmHDU *corHDU= pmHDUFromCell(corr->parent); + + if (!modHDU || !inHDU) { + psWarning("Unable to find HDU at FPA level to copy wcs!"); + } + else { + if (!pmAstromReadWCS(stack_model->fpa,model_cell->parent,inHDU->header,1.0)) { + psErrorClear(); + psWarning("Unable to read WCS astrometry from input FPA!"); + } + else { + if (!modHDU->header) { + modHDU->header = psMetadataAlloc(); + } + if (!pmAstromWriteWCS(modHDU->header, stack_model->fpa,model_cell->parent, WCS_TOLERANCE)) { + psErrorClear(); + psWarning("Unable to read WCS astrometry from input FPA!"); + } + } + if (!pmAstromReadWCS(stack_corr->fpa,corr_cell->parent,inHDU->header,1.0)) { + psErrorClear(); + psWarning("Unable to read WCS astrometry from input FPA!"); + } + else { + if (!corHDU->header) { + corHDU->header = psMetadataAlloc(); + } + if (!pmAstromWriteWCS(corHDU->header, stack_corr->fpa,corr_cell->parent, WCS_TOLERANCE)) { + psErrorClear(); + psWarning("Unable to read WCS astrometry from input FPA!"); + } + } + } // End WCS saving. + + + } // Close readout + printf(" I'm done with that readout\n"); + // Close output image + } // 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)) { + psError(psErrorCodeLast(), false, "save failure for Chip"); + return(false); + } + } // Close chip. + // Output and Close FPA + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "save failure for FPA"); + return(false); + } + psFree(view); + psFree(data->modelMap); + psFree(sizeImage); + } + + + return(true); +} + + + Index: branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackMath.c =================================================================== --- branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackMath.c (revision 37068) +++ branches/eam_branches/ipp-ops-20130712/ppBackground/src/ppBackgroundStackMath.c (revision 37068) @@ -0,0 +1,273 @@ +/** @file ppBackgroundStackArguments.c + * + * @brief + * + * @ingroup ppBackgroundStack + * + * @author Paul Price + * Copyright 2009 Institute for Astronomy, University of Hawaii + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#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. +// This is done by finding the median value of the set of (ota,u,v) values for +// all exposures. +// solution_{OTA}(u,v) = < data_{OTA}(u,v) >_{exposures} +bool ppBackgroundStackModelFitOTASolution(ppBackgroundStackData *data) { + 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))) { + if (item->type != PS_DATA_UNKNOWN) { + psWarning("This seems like I have a not-data, when I expect a data."); + continue; + } + psString workingChip = item->name; + psImage *solution = item->data.V; + + for (v = 0; v < solution->numRows; v++) { + for (u = 0; u < solution->numCols; u++) { + + 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); + if (!chipData) { continue; } + psImage *image = psMetadataLookupPtr(NULL, chipData, "bkg image"); + if (!image) { continue; } + 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); + psFree(tmp); + } // End u + } // End v + // Remove the median value from this data. We just want the tilts, not the offsets + psStatsInit(stats); + psImageStats(stats,solution,NULL,0); + for (v = 0; v < solution->numRows; v++) { + for (u = 0; u < solution->numCols; u++) { + solution->data.F32[v][u] -= stats->robustMedian; + } + } + + + } // End working chip scan + + psFree(iter); + psFree(stats); + return(true); +} + +// +// Determine the relative offset of this exposure/ota relative to the current estimate of the sky at this point. +// Note: the initial assumption of the model is a flat grid of zero, "the sky contains nothing real." +// SOLVE: model(r,d) = scale * data(r,d) + offset + solution_{OTA}(u,v) +// FOR: scale, offset +bool ppBackgroundStackDataModelFit(ppBackgroundStackData *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 *ra = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra"); + psImage *dec = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec"); + // psImage *camera= psMetadataLookupPtr(NULL, data->OTA_solutions, chipName); + psVector *obs = psVectorAllocEmpty(image->numRows, PS_TYPE_F32); + psVector *model= psVectorAllocEmpty(image->numRows, PS_TYPE_F32); + + int j = 0; + int used = 0; + for (v = 0; v < image->numRows; v++) { + for (u = 0; u < image->numCols; u++) { + if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)|| + (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { j++; continue; } + psVectorAppend(obs,image->data.F32[v][u]);// - camera->data.F32[v][u]); + psVectorAppend(model, psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u])); + j++; + used++; + } + } + + if (used > 0) { + + psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,1); + int status = psVectorFitPolynomial1D(poly,NULL,0,model,NULL,obs); + printf("in model fit loop: %d %d %d %f %f\n",status,j,used,poly->coeff[0],poly->coeff[1]); + if (status && (poly->coeff[1] != 0.0)) { + 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(expIter); + return(true); +} + + +// +// Apply the corrections to attempt to put all exposures, and all OTAs onto a common level, +// such that they should all match the true sky +// calib_{exposure,OTA}(r,d) = scale * data_{exposure,OTA}(r,d) + offset + solution_{OTA}(u,v) +bool ppBackgroundStackCalibApply(ppBackgroundStackData *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); + psImage *ra = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg ra"); + psImage *dec = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg dec"); + + 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++) { + for (u = 0; u < image->numCols; u++) { + if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)|| + (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { continue; } + + model->data.F32[v][u] = scale * image->data.F32[v][u] + offset - camera->data.F32[v][u]; + } + } + } // End chip + psFree(chipIter); + } // End smf + psFree(expIter); + return(true); +} + +// +// Determine the best estimate of the true sky from the ensemble of calibrated samples: +// model(r,d) = < calib(r,d) >_{exposures,OTAs} +// This "averaging" is done using the psImageMapClipFit. +bool ppBackgroundStackModelFit(ppBackgroundStackData *data) { + long j = 0; + int u,v; + + long used = 0; + psS16 N = psMetadataLookupS16(NULL, data->models, "N"); + psVector *X = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32); + psVector *Y = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32); + psVector *Z = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32); + psVector *E = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_F32); + psVector *mask = psVectorAllocEmpty(N * 13 * 13,PS_TYPE_VECTOR_MASK); + j = 0; + + pmFPAview *view = pmFPAviewAlloc(0); + psStats *stats = psStatsAlloc( PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV ); + + 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"); +#define DUMP_DATA 0 +#if DUMP_DATA + psImage *model = psMetadataLookupPtr(NULL, chipItem->data.md, "bkg image"); + + psF32 offset = psMetadataLookupF32(NULL,chipItem->data.md,"bkg offset"); + psF32 scale = psMetadataLookupF32(NULL,chipItem->data.md,"bkg scale"); + +#endif + for (v = 0; v < calib->numRows; v++) { + for (u = 0; u < calib->numCols; u++) { + if ((ra->data.F32[v][u] < data->x_min)||(ra->data.F32[v][u] > data->x_max)|| + (dec->data.F32[v][u] < data->y_min)||(dec->data.F32[v][u] > data->y_max)) { + j++; + continue; } + psVectorAppend(X,ra->data.F32[v][u]); + psVectorAppend(Y,dec->data.F32[v][u]); + psVectorAppend(Z,calib->data.F32[v][u]); + psVectorAppend(E,1.0); + psVectorAppend(mask,0); +/* X->data.F32[j] = ra->data.F32[v][u]; */ +/* Y->data.F32[j] = dec->data.F32[v][u]; */ +/* Z->data.F32[j] = calib->data.F32[v][u]; */ +/* E->data.F32[j] = 1.0; */ +/* mask->data.PS_TYPE_VECTOR_MASK_DATA[j] = 0; */ + used++; + j++; +#if DUMP_DATA + printf("DATA %d %ld %ld %f %f %f %f %f %f\n", + data->model_iteration,j,used, + offset,scale, + ra->data.F32[v][u], + dec->data.F32[v][u], + calib->data.F32[v][u], + model->data.F32[v][u]); +#endif + } + } + } // End chip + psFree(chipIter); + } // End smfs + printf("%ld %ld\n", j,used); + bool fitStatus; + bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, mask, 1, X, Y, Z, E); + data->model_iteration++; + psFree(expIter); + psFree(X); + psFree(Y); + psFree(Z); + psFree(E); + psFree(mask); + psFree(stats); + psFree(view); + return(status); +}