Index: /trunk/ppBackground/src/ppBackgroundStack.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStack.c (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStack.c (revision 36579) @@ -0,0 +1,72 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppBackgroundStack.h" + +int main(int argc, char *argv[]) +{ + ppBackgroundVersionPrint(); + + pmErrorRegister(); + psphotErrorRegister(); + ppBackgroundErrorRegister(); + + ppBackgroundData *data = ppBackgroundDataInit(&argc, argv); + if (!data) { + goto DIE; + } + + if (!ppBackgroundArguments(data, argc, argv)) { + goto DIE; + } + + if (!ppBackgroundCamera(data)) { + goto DIE; + } + + if (!ppBackgroundLoop(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: /trunk/ppBackground/src/ppBackgroundStack.h =================================================================== --- /trunk/ppBackground/src/ppBackgroundStack.h (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStack.h (revision 36579) @@ -0,0 +1,80 @@ +#ifndef PP_BACKGROUNDSTACK_H +#define PP_BACKGROUNDSTACK_H + +#include +#include + +#include "ppBackgroundErrorCodes.h" + +#define PPBACKGROUND_RECIPE "PPBACKGROUND" // 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 + + 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; + psF32 ra_min; + psF32 ra_max; + psF32 dec_min; + psF32 dec_max; + + psArray *stacks; // List of stacks to be corrected. + psString outRoot; // Output root name + pmConfig *config; // Configuration +} ppBackgroundStackData; + +/// Initialise data for processing +ppBackgroundData *ppBackgroundDataInit(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(ppBackgroundData *data // Data for processing + ); + +/// Loop over input data, processing +bool ppBackgroundStackLoop(ppBackgroundData *data // Data for processing + ); + +/// Determine the binning from the recipe if available. +psImageBinning *ppBackgroundStackBinningByRecipe(const psImage *image, // Image for which to generate a bg model + const pmConfig *config, // Configuration + psString recipe_name, + psString Xbin_name, + psString Ybin_name + ); + + +/// Restore the background to an image +bool ppBackgroundStackRestore( + pmChip *chip, // Chip to correct + const pmChip *background, // Chip with background model + const pmChip *pattern, // Chip with pattern + const pmFPAview *view, // View to data + pmConfig *config, // Configuration + psImageMaskType maskBad // value to use for bad pixels + ); + +/// 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: /trunk/ppBackground/src/ppBackgroundStackArguments.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackArguments.c (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStackArguments.c (revision 36579) @@ -0,0 +1,90 @@ +/** @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, ""); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-OTApath", 0, ""); + + 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")); + if (!inputs || 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: /trunk/ppBackground/src/ppBackgroundStackCamera.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackCamera.c (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStackCamera.c (revision 36579) @@ -0,0 +1,218 @@ +#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 + + psMetadataIterator *iter = psMetadataIteratorAlloc(data->contents, PS_LIST_HEAD, NULL); // Iterator + psMetadataItem *item; // Item from iteration + int i = 0; + while ((item = psMetadataGetAndIncrement(iter))) { + 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); + } + + psMetadata *input = item->data.md; // the input metadata of interest + + psString smfFileName = psMetadataLookupStr(NULL, input, "astrom"); + + // FIX Figure out how to calculate the smf name + // psString cam_file = pmFPAfileNameFromRule("PSASTRO.OUTPUT"); + // pmFPAfile *smfFile = pmFPAfileBindFromArgs(&status, NULL, config); + // pmFPAfile *smfFile = pmFPAfileDefineInput(config, + fileArguments("astrom",smfFileName,"",config); + pmFPAfile *smfFile = pmFPAfileDefineFromArgs(&status, config, "PSWARP.ASTROM","astrom"); + + // Read the SMF data + + pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + return NULL; + } + + // 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; + } + } + 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 (data->chipName && strcmp(chipName, data->chipName) != 0) { */ +/* continue; */ +/* } */ + 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); + } + // read WCS data from the corresponding header + pmHDU *hdu = pmFPAviewThisHDU (view, astromFile->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(astromFile->fpa, chip, hdu->header, 1.0)) { + psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA."); + psFree(view); + return false; + } + } + + // FIX Load model data for this chip + psString modelFileName = psMetadataLookupStr(NULL, input->models, chipName); + psFits *modelFits = psFitsOpen(modelFileName,"r"); + psImage *image = psFitsReadImage(modelFits,NULL,0); + psMetadata *header = psFitsReadHeader(NULL,modelFits); + + // Allocate the data structures for this chip + psImage *raim = psImageAlloc(image->cols,image->rows,PS_TYPE_F32); + psImage *decim = psImageAlloc(image->cols,image->rows,PS_TYPE_F32); + psImage *model = psImageAlloc(image->cols,image->rows,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->nrows; v++) { + pix->y = (v - yoffset) * ybin; + for (u = 0; u < image->ncols; u++) { + pix->x = (u - xoffset) * xbin; + psPlaneTransformApply(fp, chip->toFPA, pix); + psPlaneTransformApply(tp, smfFile->fpa->toTPA, fp); + psDeproject(sky, tp, astromFile->fpa->toSky); + + // FIX project sky to stack projection cell grid. + if (!data->radians) { + sky->r *= 180.0 / M_PI; + sky->d *= 180.0 / M_PI; + } + raim->data.F32[v][u] = sky->r; + decim->data.F32[v][u] = sky->d; + model->data.F32[v][u] = 0.0; + + // Check the bounds so we'll know how large of an area to model in the map + if (sky->r < data->ra_min) { data->ra_min = sky->r; } + else if (sky->r > data->ra_max) { data->ra_max = sky->r; } + if (sky->d < data->dec_min) { data->dec_min = sky->d; } + else if (sky->d > data->dec_max) { data->dec_max = sky->d; } + } + } + + // Attach vectors to teh structure of the chip + psMetadataAddPtr(chip->analysis,PS_LIST_TAIL,"bkg image", PS_DATA_UNKNOWN | PS_META_REPLACE, "ota space X vector", image); + psMetadataAddPtr(chip->analysis,PS_LIST_TAIL,"bkg ra", PS_DATA_UNKNOWN | PS_META_REPLACE, "ota space X vector", raim); + psMetadataAddPtr(chip->analysis,PS_LIST_TAIL,"bkg dec", PS_DATA_UNKNOWN | PS_META_REPLACE, "ota space X vector", decim); + psMetadataAddPtr(chip->analysis,PS_LIST_TAIL,"bkg calibrated data", PS_DATA_UNKNOWN | PS_META_REPLACE, "ota space corrected data", model); + + // Define default background model parameters, using the assumption: + // observed = camera + offset + scale * astrophysical + psMetadataAddF32(chip->analysis,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", 0.0); + psMetadataAddF32(chip->analysis,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 + psFree(header); + psFree(modelFits); + psFree(modelFileName); + + + // 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->cols,image->rows,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); + psFits *solutionFits = psFitsOpen(solutionFileName,"r"); + psImage *solution = psFitsReadImage(solutionFits,NULL,0); + } + } + + // FIX Hopefully this isn't going to destroy everything. + // Chip + 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)); + } + + return true; +} Index: /trunk/ppBackground/src/ppBackgroundStackData.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackData.c (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStackData.c (revision 36579) @@ -0,0 +1,70 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "ppBackgroundStack.h" + +// Destructor +static void backgroundDataFree(ppBackgroundStackData *data // Data to free + ) +{ + psFree(data->contents); + psFree(data->smfs); + + psFree(data->OTA_solutions); + psFree(data->OTApath); + + psFree(data->modelMap); + + 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->OTA_solutions = psMetadataAlloc(0); + data->fitOTAs = false; + data->OTApath = NULL; + + data->modelMap = NULL; + data->ra_min = 1e9; + data->ra_max = -1e9; + data->dec_min = 1e9; + data->dec_max = -1e9; + + data->stacks = psArrayAlloc(0); + data->outRoot = NULL; + data->config = NULL; + + data-> + + 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: /trunk/ppBackground/src/ppBackgroundStackLoop.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackLoop.c (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStackLoop.c (revision 36579) @@ -0,0 +1,75 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "ppBackgroundStack.h" + +bool ppBackgroundStackLoop(ppBackgroundStackData *data // Run-time data + ) +{ + pmConfig *config = data->config; // Configuration data + + // + // 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->fitOTAs) { + if (!ppBackgroundStackModelFitOTASolution(data)) { + // Currently can't fail. + psError(psErrorCodeLast(), false, "Error calculating the OTA solution"); + return(false); + } + } + + // This seems wrong, but I need a blank modelMap object, so we fit the zero-data we've stored in the calib objects + if (!ppBackgroundStackModelFit(data)) { + psError(psErrorCodeLast(), false, "Error determining the blank modelMap object."); + return(false); + } + + // Apply OTA solution + if (!ppBackgroundStackCalibApply(data)) { + psError(psErrorCodeLast(), false, "Error applying the calibration models."); + return(false); + } + } // End initialization block. + + // This is where a loop would likely start. + { + // Construct the offset information + if (!ppBackgroundStackDataModelFit(data)) { + psError(psErrorCodeLast(), false, "Error determining the exposure/OTA scaling."); + return(false); + } + + // Apply full correction + if (!ppBackgroundStackCalibApply(data)) { + psError(psErrorCodeLast(), false, "Error applying the calibration models."); + return(false); + } + + // Fit the new model + if (!ppBackgroundStackModelFit(data)) { + psError(psErrorCodeLast(), false, "Error determining the modelMap object."); + return(false); + } + } // End loop + + // Loop over the input images, and apply the models to construct the restored versions. + + for (i = 0; i < data->stacks->n; i++) { + // Define output image + // Calculate model for each pixel of output + // Close output image + + + return(true); +} + + + Index: /trunk/ppBackground/src/ppBackgroundStackMath.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackMath.c (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStackMath.c (revision 36579) @@ -0,0 +1,206 @@ +/** @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" + +// +// 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) { + long i; + + pmFPAview *view = pmFPAviewAllow(0); + psMetadataIterator *iter = psMetadataIteratorAlloc(data->OTA_solutions, PS_LIST_HEAD, NULL); // Iterate over all chips. + psMetadataItem *item; + psStats *stats = psStatsAlloc( PS_STAT_ROBUST_MEDIAN ); + + 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->nrows; v++) { + for (u = 0; u < solution->ncols; 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 + psStatsInit(stats); + psVectorStats(stats,tmp,NULL,NULL,0); + solution->data.F32[v][u] = stats->robustMedian; + } // End u + } // End v + } // End working chip scan + + 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) { + 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]; + + 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"); + + psVector *obs = psVectorAlloc(image->nrows * image->ncols, PS_TYPE_F32); + psVector *model= psVectorAlloc(image->nrows * image->ncols, PS_TYPE_F32); + + for (v = 0; v < image->nrows; v++) { + for (u = 0; u < image->ncols; u++) { + obs->data.F32[v * image->ncols + u] = image->data.F32[v][u]; + model->data.F32[v * image->ncols + u] = psImageMapEval(data->modelMap,ra->data.F32[v][u],dec->data.F32[v][u]); + } + } + + psPolynomial1D *poly = psPolynomial1DAlloc(PS_POLYNOMIAL_ORD,1); + status = psVectorFitPolynomial1d(poly,NULL,0,model,NULL,obs); + psMetadataAddF32(chip->analysis,PS_LIST_TAIL,"bkg offset", PS_META_REPLACE, "background offset for this exposure/ota pair", poly->coeff.F64[0]); + psMetadataAddF32(chip->analysis,PS_LIST_TAIL,"bkg scale", PS_META_REPLACE, "background scale for this exposure/ota pair", poly->coeff.F64[1]); + + psFree(tmp); + } // End OTA loop + } // End smf/exp loop + psFree(view); + 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) { + long i; + int u,v; + + pmFPAview *view = pmFPAviewAlloc(0); + 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"); + psImage *camera= psMetadataLookupPtr(NULL, data->OTA_solutions, chipName); + + psF32 offset = psMetadataLookupF32(NULL,chip->concepts,"bkg offset"); + psF32 scale = psMetadataLookupF32(NULL,chip->concepts,"bkg scale"); + + for (v = 0; v < image->nrows; v++) { + for (u = 0; u < image->ncols; u++) { + model->data.F32[v][u] = scale * image->data.F32[v][u] - offset - camera->data.F32[v][u]; + } + } + } // End chip + } // End smf + psFree(view); + 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 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); + 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"); + for (v = 0; v < image->nrows; v++) { + for (u = 0; u < image->ncols; u++) { + 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; + } + } + } // End chip + } // End smfs + + status = psImageMapClipFit(&fitStatus,data->modelMap,stats, NULL, 0, X, Y, Z, E); + psFree(X); + psFree(Y); + psFree(Z); + psFree(E); + psFree(stats); + psFree(view); +}