Index: /trunk/ppBackground/src/Makefile.am =================================================================== --- /trunk/ppBackground/src/Makefile.am (revision 36579) +++ /trunk/ppBackground/src/Makefile.am (revision 36580) @@ -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: /trunk/ppBackground/src/ppBackgroundStack.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStack.c (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStack.c (revision 36580) @@ -18,18 +18,18 @@ ppBackgroundErrorRegister(); - ppBackgroundData *data = ppBackgroundDataInit(&argc, argv); + ppBackgroundStackData *data = ppBackgroundStackDataInit(&argc, argv); if (!data) { goto DIE; } - if (!ppBackgroundArguments(data, argc, argv)) { + if (!ppBackgroundStackArguments(data, argc, argv)) { goto DIE; } - if (!ppBackgroundCamera(data)) { + if (!ppBackgroundStackCamera(data)) { goto DIE; } - if (!ppBackgroundLoop(data)) { + if (!ppBackgroundStackLoop(data)) { goto DIE; } @@ -39,18 +39,18 @@ 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 && 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) { Index: /trunk/ppBackground/src/ppBackgroundStack.h =================================================================== --- /trunk/ppBackground/src/ppBackgroundStack.h (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStack.h (revision 36580) @@ -7,5 +7,5 @@ #include "ppBackgroundErrorCodes.h" -#define PPBACKGROUND_RECIPE "PPBACKGROUND" // Recipe name +#define PPBACKGROUND_RECIPE "PPBACKGROUND_STACK" // Recipe name // Data for processing @@ -30,5 +30,5 @@ /// Initialise data for processing -ppBackgroundData *ppBackgroundDataInit(int *argc, char *argv[] // Command-line arguments +ppBackgroundStackData *ppBackgroundStackDataInit(int *argc, char *argv[] // Command-line arguments ); @@ -39,29 +39,16 @@ /// Parse camera configurations -bool ppBackgroundStackCamera(ppBackgroundData *data // Data for processing +bool ppBackgroundStackCamera(ppBackgroundStackData *data // Data for processing ); /// Loop over input data, processing -bool ppBackgroundStackLoop(ppBackgroundData *data // Data for processing +bool ppBackgroundStackLoop(ppBackgroundStackData *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 - ); +bool ppBackgroundStackModelFitOTASolution(ppBackgroundStackData *data); +bool ppBackgroundStackDataModelFit(ppBackgroundStackData *data); +bool ppBackgroundStackCalibApply(ppBackgroundStackData *data); +bool ppBackgroundStackModelFit(ppBackgroundStackData *data); - -/// 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 Index: /trunk/ppBackground/src/ppBackgroundStackArguments.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackArguments.c (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStackArguments.c (revision 36580) @@ -48,6 +48,6 @@ 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, ""); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-fitOTAs", 0, "", false); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-OTApath", 0, "", NULL); if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) { @@ -56,6 +56,6 @@ unsigned int numBad = 0; // Number of bad lines - data->contents = psMetadataConfigRead(NULL, &numBad,psMetadataLookupStr(NULL, arguments, "-input")); - if (!inputs || numBad > 0) { + 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); Index: /trunk/ppBackground/src/ppBackgroundStackCamera.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackCamera.c (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStackCamera.c (revision 36580) @@ -31,11 +31,23 @@ { bool status; // Status of file definition - + pmConfig *config = data->config; // Because I'm reusing code. + int u,v; + // FIX Figure out what stacks we have to deal with. + psString stackName = data->stacks->data[0]; + fileArguments("IMAGE", stackName, "Input image", data->config); + pmFPAfile *stack = pmFPAfileDefineFromArgs(&status, data->config, "PPBACKGROUND.STACK", + "IMAGE"); + if (!status || !stack) { + psError(psErrorCodeLast(), false, "Failed to build file from PPBACKGROUND.STACK"); + return false; + } + + // Read over the input background models. psMetadataIterator *iter = psMetadataIteratorAlloc(data->contents, PS_LIST_HEAD, NULL); // Iterator psMetadataItem *item; // Item from iteration - int i = 0; + // int i = 0; while ((item = psMetadataGetAndIncrement(iter))) { if (item->type != PS_DATA_METADATA) { - psError(PPBACKGROUND_ERR_ARGUMENTS, true + psError(PPBACKGROUND_ERR_ARGUMENTS, true, "Component %s of the input metadata is not of type METADATA", item->name); psFree(iter); @@ -47,13 +59,9 @@ 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, + // Read the smf file from this item 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)) { @@ -95,5 +103,5 @@ } // read WCS data from the corresponding header - pmHDU *hdu = pmFPAviewThisHDU (view, astromFile->fpa); + pmHDU *hdu = pmFPAviewThisHDU (view, smfFile->fpa); if (bilevelAstrometry) { if (!pmAstromReadBilevelChip (chip, hdu->header)) { @@ -108,5 +116,5 @@ // 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)) { + if (!pmAstromReadWCS(smfFile->fpa, chip, hdu->header, 1.0)) { psError(psErrorCodeLast(), false, "Unable to read WCS astrometry for input FPA."); psFree(view); @@ -116,13 +124,17 @@ // FIX Load model data for this chip - psString modelFileName = psMetadataLookupStr(NULL, input->models, chipName); + psString modelFileName = psMetadataLookupStr(NULL, input, chipName); psFits *modelFits = psFitsOpen(modelFileName,"r"); - psImage *image = psFitsReadImage(modelFits,NULL,0); + int nC, nR; + psElemType *type; + psRegion region; + psFitsImageSize(&nC,&nR,type,modelFits,region); + psImage *image = psFitsReadImage(modelFits,region,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); + 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 @@ -142,26 +154,23 @@ psSphere *sky = psSphereAlloc(); // Sky coordinates - for (v = 0; v < image->nrows; v++) { + for (v = 0; v < image->numRows; v++) { pix->y = (v - yoffset) * ybin; - for (u = 0; u < image->ncols; u++) { + 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, 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; + 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; // 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; } + 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; } } } @@ -187,5 +196,5 @@ // 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); + 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); } @@ -194,5 +203,6 @@ psStringAppend(&solutionFileName, ".%s.fits",chipName); psFits *solutionFits = psFitsOpen(solutionFileName,"r"); - psImage *solution = psFitsReadImage(solutionFits,NULL,0); + psImage *solution = psFitsReadImage(solutionFits,region,0); + psMetadataAddPtr(data->OTA_solutions,PS_LIST_TAIL, chipName, PS_DATA_UNKNOWN | PS_META_REPLACE, "OTA solution element", solution); } } @@ -215,4 +225,15 @@ } + // 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->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; + + data->modelMap = psImageMapNoImageAlloc( binning,stats); + return true; } Index: /trunk/ppBackground/src/ppBackgroundStackData.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackData.c (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStackData.c (revision 36580) @@ -36,6 +36,6 @@ data->smfs = psArrayAlloc(0); - data->OTA_solutions = psMetadataAlloc(0); - data->fitOTAs = false; + data->OTA_solutions = psMetadataAlloc(); + data->fit_OTAS = false; data->OTApath = NULL; @@ -49,6 +49,4 @@ data->outRoot = NULL; data->config = NULL; - - data-> return data; Index: /trunk/ppBackground/src/ppBackgroundStackLoop.c =================================================================== --- /trunk/ppBackground/src/ppBackgroundStackLoop.c (revision 36579) +++ /trunk/ppBackground/src/ppBackgroundStackLoop.c (revision 36580) @@ -13,11 +13,11 @@ ) { - pmConfig *config = data->config; // Configuration 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->fitOTAs) { + if (data->fit_OTAS) { if (!ppBackgroundStackModelFitOTASolution(data)) { // Currently can't fail. @@ -67,5 +67,5 @@ // 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 36580) @@ -27,6 +27,6 @@ bool ppBackgroundStackModelFitOTASolution(ppBackgroundStackData *data) { long i; - - pmFPAview *view = pmFPAviewAllow(0); + int u,v; + pmFPAview *view = pmFPAviewAlloc(0); psMetadataIterator *iter = psMetadataIteratorAlloc(data->OTA_solutions, PS_LIST_HEAD, NULL); // Iterate over all chips. psMetadataItem *item; @@ -38,14 +38,14 @@ continue; } - psString *workingChip = item->name; + psString workingChip = item->name; psImage *solution = item->data.V; - for (v = 0; v < solution->nrows; v++) { - for (u = 0; u < solution->ncols; u++) { + 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]; + pmFPAfile *smfFile = data->smfs->data[i]; pmChip *chip; // Chip from FPA @@ -84,5 +84,5 @@ for (i = 0; i < data->smfs->n; i++) { pmFPAviewReset(view); - pmFPAfile *smfFile = data->smfs.data[i]; + pmFPAfile *smfFile = data->smfs->data[i]; pmChip *chip; // Chip from FPA @@ -95,20 +95,21 @@ 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]); + psVector *obs = psVectorAlloc(image->numRows * image->numCols, PS_TYPE_F32); + psVector *model= psVectorAlloc(image->numRows * image->numCols, PS_TYPE_F32); + + for (v = 0; v < image->numRows; v++) { + for (u = 0; u < image->numCols; u++) { + obs->data.F32[v * image->numCols + u] = image->data.F32[v][u]; + model->data.F32[v * image->numCols + 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); + 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); } // End OTA loop } // End smf/exp loop @@ -127,8 +128,8 @@ pmFPAview *view = pmFPAviewAlloc(0); - stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + // psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); for (i = 0; i < data->smfs->n; i++) { pmFPAviewReset(view); - pmFPAfile *smfFile = data->smfs.data[i]; + pmFPAfile *smfFile = data->smfs->data[i]; pmChip *chip; // Chip from FPA @@ -146,6 +147,6 @@ psF32 scale = psMetadataLookupF32(NULL,chip->concepts,"bkg scale"); - for (v = 0; v < image->nrows; v++) { - for (u = 0; u < image->ncols; u++) { + for (v = 0; v < image->numRows; v++) { + for (u = 0; u < image->numCols; u++) { model->data.F32[v][u] = scale * image->data.F32[v][u] - offset - camera->data.F32[v][u]; } @@ -175,5 +176,5 @@ for (i = 0; i < data->smfs->n; i++) { pmFPAviewReset(view); - pmFPAfile *smfFile = data->smfs.data[i]; + pmFPAfile *smfFile = data->smfs->data[i]; pmChip *chip; // Chip from FPA @@ -182,10 +183,10 @@ continue; } - const char *chipName = psMetadataLookupStr(NULL, chip->concepts, "CHIP.NAME"); // Name of chip + // 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++) { + for (v = 0; v < calib->numRows; v++) { + for (u = 0; u < calib->numCols; u++) { X->data.F32[j] = ra->data.F32[v][u]; Y->data.F32[j] = dec->data.F32[v][u]; @@ -197,5 +198,7 @@ } // End smfs - status = psImageMapClipFit(&fitStatus,data->modelMap,stats, NULL, 0, X, Y, Z, E); + bool fitStatus; + bool status = psImageMapClipFit(&fitStatus,data->modelMap,stats, NULL, 0, X, Y, Z, E); + psFree(X); psFree(Y); @@ -204,3 +207,4 @@ psFree(stats); psFree(view); -} + return(status); +}