Index: /tags/pap_root_080214/psModules/src/imcombine/.cvsignore =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/.cvsignore (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/.cvsignore (revision 16477) @@ -0,0 +1,6 @@ +.deps +.libs +Makefile +Makefile.in +*.la +*.lo Index: /tags/pap_root_080214/psModules/src/imcombine/Makefile.am =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/Makefile.am (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/Makefile.am (revision 16477) @@ -0,0 +1,33 @@ +noinst_LTLIBRARIES = libpsmodulesimcombine.la + +libpsmodulesimcombine_la_CPPFLAGS = $(SRCINC) $(PSMODULES_CFLAGS) +libpsmodulesimcombine_la_LDFLAGS = -release $(PACKAGE_VERSION) + +libpsmodulesimcombine_la_SOURCES = \ + pmReadoutCombine.c \ + pmStack.c \ + pmStackReject.c \ + pmSubtraction.c \ + pmSubtractionEquation.c \ + pmSubtractionKernels.c \ + pmSubtractionMask.c \ + pmSubtractionMatch.c \ + pmSubtractionParams.c \ + pmSubtractionStamps.c \ + pmPSFEnvelope.c + +pkginclude_HEADERS = \ + pmImageCombine.h \ + pmReadoutCombine.h \ + pmStack.h \ + pmStackReject.h \ + pmSubtraction.h \ + pmSubtractionEquation.h \ + pmSubtractionKernels.h \ + pmSubtractionMask.h \ + pmSubtractionMatch.h \ + pmSubtractionParams.h \ + pmSubtractionStamps.h \ + pmPSFEnvelope.h + +CLEANFILES = *~ Index: /tags/pap_root_080214/psModules/src/imcombine/demo_psftool.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/demo_psftool.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/demo_psftool.c (revision 16477) @@ -0,0 +1,194 @@ +// This is a very quick and dirty test program for pmPSFEnvelope +// +// gcc -g demo_psftool.c -o demo_psftool `psmodules-config --cflags --libs` --std=gnu99 -Wall +// +// ./demo_psftool test MOPS.skycell.0198767.wrp*.psf +// + +#include +#include +#include + +// Add a single filename to the arguments as an array, so that it can be used with pmFPAfileBindFromArgs, etc +static void fileList(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); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, file, 0, comment, files); + psFree(files); + return; +} + + + +int main(int argc, char *argv[]) +{ + psLibInit(NULL); + pmConfig *config = pmConfigRead(&argc, argv, "PSF"); + if (!config) { + psErrorStackPrint(stderr, "Error reading configuration."); + exit(PS_EXIT_CONFIG_ERROR); + } + + psTraceSetLevel("psModules.imcombine", 5); + psTraceSetLevel("psModules.objects", 0); + psTraceSetLevel("psLib.math", 0); + + pmModelClassInit(); + + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output", argv[1]); + + psArray *files = psArrayAlloc(argc - 2); + for (int i = 2; i < argc; i++) { + psString name = NULL; // Name of file list + psStringAppend(&name, "INPUT_%d", i); + + fileList(name, argv[i], "Input PSF", config); + + pmFPAfile *file = pmFPAfileDefineFromArgs(NULL, config, "PSPHOT.PSF.LOAD", name); + psFree(name); + if (!file) { + psErrorStackPrint(stderr, "Can't define PSF file from %s --> %s", name, argv[i]); + psFree(files); + psFree(config); + exit(PS_EXIT_SYS_ERROR); + } + + files->data[i - 2] = psMemIncrRefCounter(file); + } + + pmFPAfile *outFile = pmFPAfileDefineOutput(config, NULL, "PSPHOT.PSF.SAVE"); + if (!outFile) { + psErrorStackPrint(stderr, "Can't define output PSF file"); + psFree(files); + psFree(config); + exit(PS_EXIT_SYS_ERROR); + } + outFile->save = true; + + // XXX This is a bit dodgy; should be more rigorous for a real system + { + pmFPAview *phuView = pmFPAviewAlloc(0); + phuView->chip = 0; + if (!pmFPAAddSourceFromView(outFile->fpa, "Envelope PSF", phuView, outFile->format)) { + psError(PS_ERR_UNKNOWN, false, "Unable to add PHU to output."); + psFree(phuView); + return false; + } + psFree(phuView); + } + + pmFPAview *view = pmFPAviewAlloc(0); + + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psErrorStackPrint(stderr, "Problem in I/O"); + psFree(view); + psFree(files); + psFree(config); + exit(PS_EXIT_SYS_ERROR); + } + + pmChip *chip; + while ((chip = pmFPAviewNextChip(view, outFile->fpa, 1)) != NULL) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psErrorStackPrint(stderr, "Problem in I/O"); + psFree(view); + psFree(files); + psFree(config); + exit(PS_EXIT_SYS_ERROR); + } + +#if 0 + pmFPAfileActivate(config->files, "PSPHOT.PSF.SAVE", false); + pmCell *cell; + while ((cell = pmFPAviewNextCell(view, outFile->fpa, 1)) != NULL) { + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psErrorStackPrint(stderr, "Problem in I/O"); + psFree(view); + psFree(files); + psFree(config); + exit(PS_EXIT_SYS_ERROR); + } + } +#endif + + psArray *inputs = psArrayAlloc(files->n); + int numCols = 4501, numRows = 4751; + for (int i = 0; i < files->n; i++) { + pmFPAfile *file = files->data[i]; + pmChip *chip = pmFPAviewThisChip(view, file->fpa); + + pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); + if (!psf) { + psErrorStackPrint(stderr, "Can't find PSF in file %d", i); + psFree(inputs); + psFree(files); + psFree(config); + exit(PS_EXIT_PROG_ERROR); + } + +#if 0 + pmHDU *hdu = pmHDUGetLowest(file->fpa, chip, NULL); + int imaxis1 = psMetadataLookupS32(NULL, hdu->header, "IMAXIS1"); + int imaxis2 = psMetadataLookupS32(NULL, hdu->header, "IMAXIS2"); + if (imaxis1 == 0 || imaxis2 == 0) { + psErrorStackPrint(stderr, "Size of image %d can't be determined.", i); + psFree(inputs); + psFree(files); + psFree(config); + exit(PS_EXIT_SYS_ERROR); + } + if (numCols == 0 && numRows == 0) { + numCols = imaxis1; + numRows = imaxis2; + } else if (imaxis1 != numCols || imaxis2 != numRows) { + psErrorStackPrint(stderr, "Image %d differs in size: %dx%d vs %dx%d", + i, imaxis1, imaxis2, numCols, numRows); + psFree(inputs); + psFree(files); + psFree(config); + exit(PS_EXIT_SYS_ERROR); + } +#endif + + inputs->data[i] = psMemIncrRefCounter(psf); + } + + pmPSF *psf = pmPSFEnvelope(numCols, numRows, inputs, 5, 20, "PS_MODEL_RGAUSS", 3, 3); + psFree(inputs); + if (!psf) { + psErrorStackPrint(stderr, "Can't generate envelope PSF."); + psFree(config); + exit(PS_EXIT_SYS_ERROR); + } + + psMetadataAddPtr(chip->analysis, PS_LIST_TAIL, "PSPHOT.PSF", PS_DATA_UNKNOWN, "Envelope PSF", psf); + psFree(psf); + chip->data_exists = true; + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psErrorStackPrint(stderr, "Problem in I/O"); + psFree(view); + psFree(files); + psFree(config); + exit(PS_EXIT_SYS_ERROR); + return false; + } + } + + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psErrorStackPrint(stderr, "Problem in I/O"); + psFree(view); + psFree(files); + psFree(config); + exit(PS_EXIT_SYS_ERROR); + } + + psFree(config); + + exit(PS_EXIT_SUCCESS); +} Index: /tags/pap_root_080214/psModules/src/imcombine/pmImageCombine.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmImageCombine.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmImageCombine.c (revision 16477) @@ -0,0 +1,683 @@ +/** @file pmImageCombine.c + * + * This file will perform image combination of several images of the + * same field, produce a list of questionable pixels, then tag some + * of those pixels as cosmic rays. + * + * @author Paul Price, IfA (original prototype) + * @author GLG, MHPCC + * + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-04-04 22:42:48 $ + * + * XXX: pmRejectPixels() has a known bug with the pmImageTransform() call. + * + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include "pslib.h" + +#define PIXEL_LIST_BUFFER 100 // Size of the pixel list buffer + +// Data structure for use as a buffer in combining pixels +typedef struct +{ + psVector *pixels; // Pixel values + psVector *masks; // Pixel masks + psVector *errors; // Pixel errors + psStats *stats; // Statistics to use with combination +} +combineBuffer; + +void combineBufferFree(combineBuffer *buffer) +{ + psFree(buffer->pixels); + psFree(buffer->masks); + psFree(buffer->errors); + psFree(buffer->stats); +} + +combineBuffer *combineBufferAlloc(long numImages // Number of images that will be combined + ) +{ + combineBuffer *buffer = psAlloc(sizeof(combineBuffer)); + psMemSetDeallocator(buffer, (psFreeFunc)combineBufferFree); + + buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32); + buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK); + buffer->errors = psVectorAlloc(numImages, PS_TYPE_F32); + buffer->stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + + return buffer; +} + + +static bool combinePixels(psImage *combine, // Combined image, for output + psArray *questionablePixels, // Array of rejection masks + int x, int y, // Position in the images + const psArray *images, // Array of input images + const psArray *errors, // Array of input error images + const psArray *masks, // Array of input masks + psU32 maskVal, // Mask value + psS32 numIter, // Number of rejection iterations + psF32 sigmaClip, // Number of standard deviations at which to reject + combineBuffer *buffer // Buffer for combination; to avoid multiple allocations + ) +{ + assert(combine); + assert(x >= 0 && x < combine->numCols); + assert(y >= 0 && y < combine->numRows); + assert(images); + int numImages = images->n; // Number of images to combine + if (masks) { + assert(masks->n == numImages); + } + if (errors) { + assert(errors->n == numImages); + } + assert(numIter >= 0); + assert(sigmaClip > 0); + assert(!questionablePixels || (questionablePixels && questionablePixels->n == numImages)); + + if (buffer) { + psMemIncrRefCounter(buffer); + } else { + buffer = combineBufferAlloc(numImages); + } + + psVector *pixelData = buffer->pixels; // Values for the pixel of interest + psVector *pixelErrors = buffer->errors; // Errors for the pixel of interest + psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest + psStats *stats = buffer->stats; // Statistics for combination + + // + // Loop through each image, extract the pixel/mask/error data into psVectors. + // + if (!masks) { + psVectorInit(pixelMasks, 0); + } + if (!errors) { + pixelErrors = NULL; + } + for (int i = 0; i < numImages; i++) { + // Set the pixel data + psImage *image = images->data[i]; // Image of interest + pixelData->data.F32[i] = image->data.F32[y][x]; + // Set the pixel mask data, if necessary + if (masks) { + psImage *mask = masks->data[i]; // Mask of interest + pixelMasks->data.U8[i] = mask->data.U8[y][x]; + } // Set the pixel error data, if necessary + if (errors) { + psImage *error = errors->data[i]; // Error image of interest + pixelErrors->data.F32[i] = error->data.F32[y][x]; + } + } + + // + // Iterate on the pixels, rejecting outliers + // + for (int iter = 0; iter < numIter; iter++) { + // Combine all the pixels, using the specified stat. + if (!psVectorStats(stats, pixelData, pixelErrors, pixelMasks, maskVal)) { + combine->data.F32[y][x] = NAN; + psFree(buffer); + return false; + } + float combinedPixel = stats->sampleMean; // Value of the combination + + if (iter == 0) { + // Save the value produced with no rejection, since it may be useful later + // (if the rejection turns out to be unnecessary) + combine->data.F32[y][x] = combinedPixel; + } + + // + // Reject all pixels that lie more that sigmaClip standard deviations from + // the combined pixel value. + // + int numRejects = 0; // Number of rejections + float stdev = stats->sampleStdev; + for (int i = 0; i < numImages; i++) { + if (!(pixelMasks->data.U8[i] & maskVal) && + fabs(pixelData->data.F32[i] - combinedPixel) > sigmaClip * stdev) { + // Reject pixel as questionable + numRejects++; + pixelMasks->data.U8[i] = maskVal; + if (questionablePixels) { + // Mark the pixel as questionable + psPixels *qp = questionablePixels->data[i]; // Questionable pixels for this image + int qpNum = qp->n; // Number of QPs in the image of interest + if (qpNum >= qp->nalloc) { + // Grow dynamically, if required + qp = psPixelsRealloc(qp, qp->nalloc + PIXEL_LIST_BUFFER); + questionablePixels->data[i] = qp; + } + qp->data[qpNum].x = x; + qp->data[qpNum].y = y; + qp->n++; + } + } + } + + // + // If the number of rejected pixels is zero, then there's no point continuing the loop. + // + if (numRejects == 0) { + break; + } + + // + // XXX: Is it possible to have all pixels rejected? If so, we should exit the loop. + // + } + + psFree(buffer); + return true; +} + +psImage *pmCombineImages( + psImage *combine, ///< Combined image (output) + psArray **questionablePixels, ///< Array of rejection masks + const psArray *images, ///< Array of input images + const psArray *errors, ///< Array of input error images + const psArray *masks, ///< Array of input masks + psU32 maskVal, ///< Mask value + const psPixels *pixels, ///< Pixels to combine + psS32 numIter, ///< Number of rejection iterations + psF32 sigmaClip ///< Number of standard deviations at which to reject +) +{ + psTrace("psModules.imcombine", 3, "Calling pmCombineImages(%ld)\n", images->n); + + PS_ASSERT_ARRAY_NON_NULL(images, NULL); + PS_ASSERT_INT_POSITIVE(images->n, NULL); + long numImages = images->n; // Number of images + int numCols = ((psImage*)images->data[0])->numCols; // Size in x + int numRows = ((psImage*)images->data[0])->numRows; // Size in y + + if (combine) { + PS_ASSERT_IMAGE_NON_NULL(combine, NULL); + PS_ASSERT_IMAGE_SIZE(combine, numCols, numRows, NULL); + } + if (questionablePixels && !*questionablePixels) { + PS_ASSERT_ARRAY_NON_NULL(*questionablePixels, NULL); + PS_ASSERT_ARRAY_SIZE(*questionablePixels, numImages, 0); + } + for (int i = 1; i < images->n; i++) { + psImage *image = images->data[i]; // Image of interest + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); + PS_ASSERT_IMAGE_SIZE(image, numCols, numRows, NULL); + } + if (errors) { + PS_ASSERT_ARRAYS_SIZE_EQUAL(images, errors, NULL); + for (int i = 0; i < images->n; i++) { + psImage *error = errors->data[i]; + PS_ASSERT_IMAGE_SIZE(error, numCols, numRows, NULL); + PS_ASSERT_IMAGE_TYPE(error, PS_TYPE_F32, NULL); + } + } + if (masks) { + PS_ASSERT_ARRAYS_SIZE_EQUAL(images, masks, NULL); + for (int i = 0; i < images->n; i++) { + psImage *mask = masks->data[i]; + PS_ASSERT_IMAGE_SIZE(mask, numCols, numRows, NULL); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, NULL); + } + } + PS_ASSERT_INT_POSITIVE(numIter, NULL); + PS_ASSERT_FLOAT_LARGER_THAN(sigmaClip, 0.0, NULL); + + // Allocate and initialize the combined image, if necessary. + if (!combine) { + combine = psImageAlloc(numCols, numRows, PS_TYPE_F32); + if (pixels) { + // Set everything we're not working on to NAN + psImageInit(combine, NAN); + } + } + + // + // Allocate the questionablePixels psArray, if necesssary, then create a psPixels + // struct for each image. + // + if (questionablePixels) { + if (*questionablePixels == NULL) { + *questionablePixels = psArrayAlloc(numImages); + } else if ((*questionablePixels)->n != numImages) { + *questionablePixels = psArrayRealloc(*questionablePixels, numImages); + } + for (int i = 0; i < numImages; i++) { + psFree((*questionablePixels)->data[i]); + (*questionablePixels)->data[i] = psPixelsAlloc(PIXEL_LIST_BUFFER); + } + } + + combineBuffer *buffer = combineBufferAlloc(numImages); // Buffer for combination + + if (pixels) { + // Only those specified pixels should be combined. + + for (int p = 0; p < pixels->n; p++) { + int x = pixels->data[p].x; // Column of interest + int y = pixels->data[p].y; // Row of interest + + if (!combinePixels(combine, questionablePixels ? *questionablePixels : NULL, x, y, + images, errors, masks, maskVal, numIter, sigmaClip, buffer)) { + // Bad pixel --- no big deal + psErrorClear(); + } + } + } else { + // + // We get here if there is a NULL list of pixels to combine. + // Therefore, we combine all pixels in all images. + // + + // + // Loop over all pixels in all images, set the appropriate data, mask, + // error vectors, call psVectorStats(), and set the result in the + // combine image. + // + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (!combinePixels(combine, questionablePixels ? *questionablePixels : NULL, x, y, + images, errors, masks, maskVal, numIter, sigmaClip, buffer)) { + // Bad pixel --- no big deal + psErrorClear(); + } + } + } + } + + psFree(buffer); + + psTrace("psModules.imcombine", 3, "Exiting pmCombineImages(%ld)\n", images->n); + return combine; +} + + +/****************************************************************************** +XXX: Directly from Paul Price + *****************************************************************************/ +static psF32 CalcGradient( + psImage *image, + psImage *imageMask, + psS32 x, + psS32 y +) +{ + psTrace("psModules.imcombine", 4, "Calling CalcGradient(%d, %d)\n", x, y); + int num = 0; + psVector *pixels = psVectorAlloc(8, PS_TYPE_F32); // Array of pixels + psVector *mask = psVectorAlloc(8, PS_TYPE_U8); // Corresponding mask + + // Get limits + int xMin = PS_MAX(x - 1, 0); + int xMax = PS_MIN(x + 1, image->numCols - 1); + int yMin = PS_MAX(y - 1, 0); + int yMax = PS_MIN(y + 1, image->numRows - 1); + if (imageMask != NULL) { + for (int j = yMin; j <= yMax; j++) { + for (int i = xMin; i <= xMax; i++) { + if ((i != x) && (j != y) && (0 == imageMask->data.U8[j][i])) { + pixels->data.F32[num] = image->data.F32[j][i]; + mask->data.U8[num] = 0; + num++; + } else { + mask->data.U8[num] = 1; + } + } + } + } else { + // + // This code is simply the previous loop without the imageMask. + // XXX: Consider restructuring this. + // + for (int j = yMin; j <= yMax; j++) { + for (int i = xMin; i <= xMax; i++) { + if ((i != x) && (j != y)) { + pixels->data.F32[num] = image->data.F32[j][i]; + mask->data.U8[num] = 0; + num++; + } else { + mask->data.U8[num] = 1; + } + } + } + } + + pixels->n = num; + mask->n = num; + + // Get the median + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); + psVectorStats(stats, pixels, NULL, mask, 1); + float median = stats->sampleMedian; + psFree(stats); + psFree(pixels); + psFree(mask); + + psTrace("psModules.imcombine", 4, "Exiting CalcGradient(%d, %d)\n", x, y); + return(median / image->data.F32[y][x]); +} + +/****************************************************************************** +DetermineRegion(image, myOutToIn): for a psImage and a psPlaneTransform to that +image, this routine determines the size of the input image which maps to that +image, and returns the result in a psRegion struct. + +XXX: Basically, this routine is only guaranteed to work if the transform is +linear. + +XXX: Shouldn't this functionality be part of psImageTransform()? + *****************************************************************************/ +static psRegion DetermineRegion(psImage *image, + psPlaneTransform *myOutToIn) +{ + psTrace("psModules.imcombine", 4, "Calling DetermineRegion()\n"); + psRegion myRegion; + myRegion.x0 = PS_MAX_F32; + myRegion.x1 = PS_MIN_F32; + myRegion.y0 = PS_MAX_F32; + myRegion.y1 = PS_MIN_F32; + psPlane in; + psPlane out; + + in.x = 0.0; + in.y = 0.0; + + psPlaneTransformApply(&out, myOutToIn, &in); + if (out.x < myRegion.x0) { + myRegion.x0 = out.x; + } + if (out.x > myRegion.x1) { + myRegion.x1 = out.x; + } + if (out.y < myRegion.y0) { + myRegion.y0 = out.y; + } + if (out.y > myRegion.y1) { + myRegion.y1 = out.y; + } + + in.x = (psF32) (image->numCols); + in.y = 0.0; + psPlaneTransformApply(&out, myOutToIn, &in); + if (out.x < myRegion.x0) { + myRegion.x0 = out.x; + } + if (out.x > myRegion.x1) { + myRegion.x1 = out.x; + } + if (out.y < myRegion.y0) { + myRegion.y0 = out.y; + } + if (out.y > myRegion.y1) { + myRegion.y1 = out.y; + } + + in.x = (psF32) (image->numCols); + ; + in.y = 0.0; + psPlaneTransformApply(&out, myOutToIn, &in); + if (out.x < myRegion.x0) { + myRegion.x0 = out.x; + } + if (out.x > myRegion.x1) { + myRegion.x1 = out.x; + } + if (out.y < myRegion.y0) { + myRegion.y0 = out.y; + } + if (out.y > myRegion.y1) { + myRegion.y1 = out.y; + } + + in.x = (psF32) (image->numCols); + in.y = (psF32) (image->numRows); + psPlaneTransformApply(&out, myOutToIn, &in); + if (out.x < myRegion.x0) { + myRegion.x0 = out.x; + } + if (out.x > myRegion.x1) { + myRegion.x1 = out.x; + } + if (out.y < myRegion.y0) { + myRegion.y0 = out.y; + } + if (out.y > myRegion.y1) { + myRegion.y1 = out.y; + } + + psTrace("psModules.imcombine", 4, "Exiting DetermineRegion()\n"); + return(myRegion); +} + +/****************************************************************************** +XXX: Don't we have a psLib function for this? + *****************************************************************************/ +static psImage *ImageConvertF32(psImage *image) +{ + psTrace("psModules.imcombine", 4, "Calling ImageConvertF32()\n"); + psImage *imgF32 = psImageAlloc(image->numCols, image->numRows, PS_TYPE_F32); + + for (psS32 i = 0 ; i < image->numRows ; i++) { + for (psS32 j = 0 ; j < image->numCols ; j++) { + imgF32->data.F32[i][j] = (psF32) image->data.U8[i][j]; + } + } + + psTrace("psModules.imcombine", 4, "Exiting ImageConvertF32()\n"); + return(imgF32); +} + + +// +// The following macros define how big the initial pixel list will be, and +// how much it should be incremented when realloc'ed. +// +#define PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH 100 +#define PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC 100 +/****************************************************************************** +pmRejectPixels(images, errors, inToOut, outToIn, rejThreshold, +gradLimit) + +XXX: Optimization: we don't need to transform the entire mask image. +XXX: The inToOut and outToIn transforms are confusing. Verify that what + I think they mean syncs with PWP. + *****************************************************************************/ +psArray *pmRejectPixels( + const psArray *images, ///< Array of input images + const psArray *masks, ///< Array of input image masks + const psArray *errors, ///< The pixels which were rejected in the combination + const psArray *inToOut, ///< Transformation from input to output system + const psArray *outToIn, ///< Transformation from output to input system + psF32 rejThreshold, ///< Rejection threshold + psF32 gradLimit ///< Gradient limit +) +{ + psTrace("psModules.imcombine", 3, "Calling pmRejectPixels()\n"); + PS_ASSERT_PTR_NON_NULL(images, NULL); + for (psS32 im = 0 ; im < images->n ; im++) { + psImage *tmpImage = (psImage *) images->data[im]; + PS_ASSERT_IMAGE_NON_NULL(tmpImage, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(tmpImage, NULL); + PS_ASSERT_IMAGE_TYPE(tmpImage, PS_TYPE_F32, NULL); + if (masks != NULL) { + PS_ASSERT_INT_EQUAL(images->n, masks->n, NULL); + psImage *tmpMask = (psImage *) masks->data[im]; + PS_ASSERT_IMAGE_NON_NULL(tmpMask, NULL); + PS_ASSERT_IMAGE_NON_EMPTY(tmpMask, NULL); + PS_ASSERT_IMAGE_TYPE(tmpMask, PS_TYPE_F32, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(tmpImage, tmpMask, NULL); + } + PS_ASSERT_IMAGES_SIZE_EQUAL(((psImage *) images->data[0]), tmpImage, NULL); + } + PS_ASSERT_PTR_NON_NULL(errors, NULL); + PS_ASSERT_PTR_NON_NULL(inToOut, NULL); + PS_ASSERT_PTR_NON_NULL(outToIn, NULL); + // Ensure that the psArray parameters have an element for each image. + psS32 numImages = images->n; + PS_ASSERT_INT_EQUAL(numImages, errors->n, NULL); + PS_ASSERT_INT_EQUAL(numImages, inToOut->n, NULL); + PS_ASSERT_INT_EQUAL(numImages, outToIn->n, NULL); + + // + // Create the psArray of psPixelLists, one for each image, for rejected pixels. + // + psArray *rejects = psArrayAlloc(numImages); + for (psS32 im = 0 ; im < numImages ; im++) { + rejects->data[im] = (psPtr *) psPixelsAlloc(PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH); + ((psPixels *)(rejects->data[im]))->n = ((psPixels *)(rejects->data[im]))->nalloc; + psPixels *pixels = (psPixels *) rejects->data[im]; + pixels->n = 0; + } + // + // rPtr is used to maintain a count of the questionable pixels for each image. + // + psVector *rPtr = psVectorAlloc(numImages, PS_TYPE_S32); + psVectorInit(rPtr, 0); + + psS32 numCols = ((psImage *) images->data[0])->numCols; + psS32 numRows = ((psImage *) images->data[0])->numRows; + psRegion myRegion = psRegionSet(0, numCols-1, 0, numRows-1); + psU32 maskVal = 1; // XXX: Is this appropriate? + + psPlane *inCoords = psAlloc(sizeof(psPlane)); + psPlane *outCoords = psAlloc(sizeof(psPlane)); + + for (psS32 im = 0 ; im < numImages ; im++) { + // + // Extract data from psArrays. + // + psPixels *pixelList = (psPixels *) errors->data[im]; + + psImage *currImage = (psImage *) images->data[im]; + myRegion.x0 = 0; + myRegion.x1 = currImage->numCols; + myRegion.y0 = 0; + myRegion.y1 = currImage->numRows; + psPlaneTransform *myInToOut = (psPlaneTransform *) inToOut->data[im]; + psPlaneTransform *myOutToIn = (psPlaneTransform *) outToIn->data[im]; + + // + // Create a psU8 mask image from the list of cosmic pixels. + // + psImage *maskImage = NULL; + maskImage = psPixelsToMask(maskImage, pixelList, myRegion, maskVal); + psImage *maskImageF32 = ImageConvertF32(maskImage); + + // + // Transform that mask image into detector coordinate space + // + psRegion myRegionXForm = DetermineRegion(maskImageF32, myOutToIn); + psImage *transformedImage = psImageTransform(NULL, NULL, maskImageF32, NULL, + 0, myOutToIn, myRegionXForm, NULL, + PS_INTERPOLATE_BILINEAR, 0); + + // + // Loop over all cosmic pixels. Transform their coords to detector space. + // If the value of the transformed mask is larger than rejThreshold, then + // this might be a cosmic ray pixel. We then calculate the mean gradient + // in other images. + // + + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, + transformedImage, NULL, NULL, + 0, 0.0, 0.0, 0, 0, 0.0); + + for (psS32 p = 0 ; p < pixelList->n ; p++) { + inCoords->x = 0.5 + (psF32) (pixelList->data[p]).x; + inCoords->y = 0.5 + (psF32) (pixelList->data[p]).y; + psPlaneTransformApply(outCoords, myInToOut, inCoords); + double maskVal; + if (!psImageInterpolate(&maskVal, NULL, NULL, outCoords->x, outCoords->y, interp)) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate image."); + psFree(interp); + psFree(maskImage); + psFree(maskImageF32); + psFree(transformedImage); + psFree(inCoords); + psFree(outCoords); + psFree(rejects); + return NULL; + } + if (maskVal > rejThreshold) { + + // This is a possible cosmic array pixel. We must calculate the gradient + // at this location in all input images. + psF32 meanGrads = 0.0; + psS32 numGrads = 0; + // + // Loop through all other images, calculate their mean gradient. + // + for (psS32 otherImg = 0 ; otherImg < numImages ; otherImg++) { + if (im != otherImg) { + // Map the outCoords to inCoords that for otherImg space. + psImage *tmpMask = NULL; + if (masks != NULL) { + tmpMask = masks->data[otherImg]; + } + psPlaneTransformApply(inCoords, + (psPlaneTransform * )outToIn->data[otherImg], + outCoords); + psS32 xPix = (int)(inCoords->x + 0.5); + psS32 yPix = (int)(inCoords->y + 0.5); + if ((xPix >= 0) && (xPix <= ((psImage*)(images->data[otherImg]))->numCols - 1) && + (yPix >= 0) && (yPix <= ((psImage*)(images->data[otherImg]))->numRows - 1)) { + meanGrads += CalcGradient(images->data[otherImg], tmpMask, xPix, yPix); + numGrads++; + } + } + } + if (numGrads > 0) { + meanGrads /= (psF32) numGrads; + } else { + // XXX: my idea. Verify with PWP: + meanGrads = 1.0 + gradLimit; + } + + // XXX: The SDRS and the prototype code differ significantly here: + // if (CalcGradient(inputs->data[i], pixelList->data.x, pixelList->data.y) < (gradLimit * meanGrads)) { + if (meanGrads < gradLimit) { + // + // Add this to the list of questionable pixels. We must ensure that the + // pixelList is large enough; if not, we realloc() + // + psS32 ptr = rPtr->data.S32[im]; + psPixels *pixelListPtr = (psPixels *) rejects->data[im]; + if (ptr >= pixelListPtr->nalloc) { + rejects->data[im] = (psPtr *) psPixelsRealloc(((psPixels *) rejects->data[im]), + ((((psPixels *) rejects->data[im])->nalloc) + PS_REJECT_PIXEL_INITIAL_PIXEL_LIST_LENGTH_INC)); + // XXX: Can the realloc() fail? Must we check for NULL? + } + + ((psPixels *) rejects->data[im])->data[ptr].x = (pixelList->data[p]).x; + ((psPixels *) rejects->data[im])->data[ptr].y = (pixelList->data[p]).y; + (rPtr->data.S32[im])++; + // XXX: this pixel ->n increment is wierd + (((psPixels *) rejects->data[im])->n)++; + } + } + } + + psFree(interp); + psFree(maskImage); + psFree(maskImageF32); + psFree(transformedImage); + } + + psFree(inCoords); + psFree(outCoords); + psTrace("psModules.imcombine", 3, "Exiting pmRejectPixels()\n"); + return(rejects); +} Index: /tags/pap_root_080214/psModules/src/imcombine/pmImageCombine.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmImageCombine.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmImageCombine.h (revision 16477) @@ -0,0 +1,44 @@ +/* @file pmImageCombine.h + * + * This file will perform image combination of several images of the + * same field, produce a list of questionable pixels, then tag some + * of those pixels as cosmic rays. + * + * @author Paul Price, IfA (original prototype) + * @author GLG, MHPCC + * + * @version $Revision: 1.6 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-03-30 21:12:56 $ + * Copyright 2004 Maui High Performance Computing Center, University of Hawaii + */ + +#ifndef PM_IMAGE_COMBINE_H +#define PM_IMAGE_COMBINE_H + +/// @addtogroup imcombine Image Combinations +/// @{ + +psImage *pmCombineImages( + psImage *combine, ///< Combined image (output) + psArray **questionablePixels, ///< Array of rejection masks + const psArray *images, ///< Array of input images + const psArray *errors, ///< Array of input error images + const psArray *masks, ///< Array of input masks + psU32 maskVal, ///< Mask value + const psPixels *pixels, ///< Pixels to combine + psS32 numIter, ///< Number of rejection iterations + psF32 sigmaClip ///< Number of standard deviations at which to reject +); + +psArray *pmRejectPixels( + const psArray *images, ///< Array of input images + const psArray *masks, ///< Array of input image masks + const psArray *errors, ///< The pixels which were rejected in the combination + const psArray *inToOut, ///< Transformation from input to output system + const psArray *outToIn, ///< Transformation from output to input system + psF32 rejThreshold, ///< Rejection threshold + psF32 gradLimit ///< Gradient limit +); + +/// @} +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmPSFEnvelope.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmPSFEnvelope.c (revision 16477) @@ -0,0 +1,302 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmReadoutFake.h" + +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmModelClass.h" +#include "pmModelUtils.h" +#include "pmPeaks.h" +#include "pmSource.h" +#include "pmSourceUtils.h" +#include "pmSourceFitModel.h" +#include "pmPSFtry.h" + + +#include "pmPSFEnvelope.h" + + + +//#define TESTING // Enable test output +#define PEAK_FLUX 1.0e3 // Peak flux for each source +#define SKY_VALUE 0.0e0 // Sky value for fake image +#define WEIGHT_VAL 1.0e0 // Weighting for image +#define PSF_STATS PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV // Statistics options for measuring PSF + + + +// XXX To do: +// +// * PSF variation when only a portion of the image is present (e.g., the edge of an FPA overlapping a +// skycell) may mean a disastrously weird PSF in the missing regions. To counter this, get a region of +// validity for each PSF (perhaps from an associated mask, or have the user work it out), and taper the PSF +// when outside this region (perhaps multiply the peak flux by a Gaussian whose arguments are the distance +// from the valid region, and a width that the user supplies). + + +// We deliberately do not include the calculation of and storing of residuals (data - model) for the PSF +// model, because (1) there is no code in psModules to do this, and we're not going to implement it here; and +// (2) this is intended to generate "nice" or "ideal" PSFs to feed into pmSubtraction (PSF matching code), so +// any residuals will hopefully be dealt with by that. + + +pmPSF *pmPSFEnvelope(int numCols, int numRows, // Size of original image + const psArray *inputs, // Input PSF models + int instances, // Number of instances per dimension + int radius, // Radius of each PSF + const char *modelName,// Name of PSF model to use + int xOrder, int yOrder // Order for PSF variation fit + ) +{ + PS_ASSERT_INT_POSITIVE(numCols, NULL); + PS_ASSERT_INT_POSITIVE(numRows, NULL); + PS_ASSERT_ARRAY_NON_NULL(inputs, NULL); + PS_ASSERT_INT_POSITIVE(instances, NULL); + PS_ASSERT_INT_POSITIVE(radius, NULL); + PS_ASSERT_STRING_NON_EMPTY(modelName, NULL); + PS_ASSERT_INT_NONNEGATIVE(xOrder, NULL); + PS_ASSERT_INT_NONNEGATIVE(yOrder, NULL); + + float xOrigSpacing = (float)(numCols - 2 * radius) / (float)(instances - 1); // Spacing between instances + float yOrigSpacing = (float)(numRows - 2 * radius) / (float)(instances - 1); // Spacing between instances + int fakeSpacing = 2 * radius + 1; // Spacing between instances (x and y) in the fake image + int fakeSize = instances * fakeSpacing; // Size of fake image + + // Generate list of fake sources (instances of the PSF) + int numFakes = PS_SQR(instances); // Number of fake sources + psArray *fakes = psArrayAlloc(numFakes); // Fake sources + psVector *xOffset = psVectorAlloc(numFakes, PS_TYPE_S32); // X offset from fake position to image + psVector *yOffset = psVectorAlloc(numFakes, PS_TYPE_S32); // Y offset from fake position to image + for (int j = 0, index = 0; j < instances; j++) { + float yOrig = j * yOrigSpacing + radius; // Source position in original image + float yFake = j * fakeSpacing + radius; // Position in fake image + int dy = yFake - yOrig; // Difference between fake and original position + + for (int i = 0; i < instances; i++, index++) { + float xOrig = i * xOrigSpacing + radius; // Source position in original image + float xFake = i * fakeSpacing + radius; // Position in fake image + int dx = xFake - xOrig; // Difference between fake and original position + + pmSource *fake = pmSourceAlloc(); // Fake source + fake->peak = pmPeakAlloc(xFake - dx, yFake - dy, PEAK_FLUX, PM_PEAK_LONE); + fake->type = PM_SOURCE_TYPE_STAR; + fake->psfMag = -2.5 * log10(PEAK_FLUX); + + psTrace("psModules.imcombine", 5, "Source %d: %.2f,%.2f\n", + index, xOrig, yOrig); + + fakes->data[index] = fake; + xOffset->data.S32[index] = dx; + yOffset->data.S32[index] = dy; + } + } + + // Generate fake images with each PSF, and take the envelope + psImage *envelope = psImageAlloc(fakeSize, fakeSize, PS_TYPE_F32); // Image with envelope of PSFs + psImageInit(envelope, SKY_VALUE); + pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout + for (int i = 0; i < inputs->n; i++) { + pmPSF *psf = inputs->data[i]; // PSF of interest + pmResiduals *resid = psf->residuals;// PSF residuals + psf->residuals = NULL; + if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, xOffset, yOffset, psf, + NAN, radius, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout."); + psFree(envelope); + psFree(yOffset); + psFree(xOffset); + psFree(fakes); + psf->residuals = resid; + return NULL; + } + psf->residuals = resid; + + // Need to renormalise sources so they all have the same peak. You would think they do have the same + // peak already, but it seems that the residual map messes things up by adding extra flux + for (int j = 0; j < numFakes; j++) { + pmSource *source = fakes->data[j]; // Fake source + float x = source->peak->xf + xOffset->data.S32[j]; // x coordinate of source + float y = source->peak->yf + yOffset->data.S32[j]; // y coordinate of source + + double flux = fakeRO->image->data.F32[(int)y][(int)x]; + float norm = PEAK_FLUX / flux; // Normalisation for source + psRegion region = psRegionSet(x - radius, x + radius, y - radius, y + radius); // PSF region + psImage *subImage = psImageSubset(fakeRO->image, region); // Subimage of fake PSF + psImage *subEnv = psImageSubset(envelope, region); // Subimage of envelope + psBinaryOp(subImage, subImage, "*", psScalarAlloc(norm, PS_TYPE_F32)); + psBinaryOp(subEnv, subEnv, "MAX", subImage); + psFree(subImage); + psFree(subEnv); + } + +#ifdef TESTING + { + // Write out the PSF field + psString name = NULL; + psStringAppend(&name, "psf_field_%03d.fits", i); + psFits *fits = psFitsOpen(name, "w"); + psFitsWriteImage(fits, NULL, fakeRO->image, 0, NULL); + psFitsClose(fits); + psFree(name); + } +#endif + + } + psFree(fakeRO); + +#ifdef TESTING + { + // Write out the envelope + psFits *fits = psFitsOpen("psf_field_envelope.fits", "w"); + psFitsWriteImage(fits, NULL, envelope, 0, NULL); + psFitsClose(fits); + } +#endif + + // Put the fake sources onto a full-size image + pmReadout *readout = pmReadoutAlloc(NULL); // Readout to contain envelope pixels + readout->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(readout->image, 0.0); + for (int i = 0; i < numFakes; i++) { + pmSource *source = fakes->data[i]; // Fake source + // Position of source on fake image + int xFake = source->peak->x + xOffset->data.S32[i]; + int yFake = source->peak->y + yOffset->data.S32[i]; + psRegion region = psRegionSet(xFake - radius, xFake + radius, + yFake - radius, yFake + radius); // PSF region + psImage *subImage = psImageSubset(envelope, region); // Subimage of fake PSF + + // Position of source on "real" image + int x0 = source->peak->x - radius; + int y0 = source->peak->y - radius; + + if (!psImageOverlaySection(readout->image, subImage, x0, y0, "=")) { + psError(PS_ERR_UNKNOWN, false, "Unable to overlay PSF"); + psFree(readout); + psFree(xOffset); + psFree(yOffset); + psFree(fakes); + return NULL; + } + } + psFree(xOffset); + psFree(yOffset); + + readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "+", psScalarAlloc(1.0, PS_TYPE_F32)); + readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); + psImageInit(readout->mask, 0); + +#ifdef TESTING + { + // Write out the envelope + psFits *fits = psFitsOpen("psf_field_full.fits", "w"); + psFitsWriteImage(fits, NULL, readout->image, 0, NULL); + psFitsClose(fits); + } +#endif + + // Reset the sources to point to the new pixels, and measure the moments in preparation for PSF fitting + for (int i = 0; i < numFakes; i++) { + pmSource *source = fakes->data[i]; // Fake source + float x = source->peak->xf; // x coordinates of source + float y = source->peak->yf; // y coordinates of source + + psFree(source->pixels); + psFree(source->weight); + psFree(source->maskView); + psFree(source->maskObj); + source->pixels = NULL; + source->weight = NULL; + source->maskView = NULL; + source->maskObj = NULL; + + if (!pmSourceDefinePixels(source, readout, x, y, (float)radius)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define pixels for source."); + psFree(readout); + psFree(fakes); + return NULL; + } + + if (!pmSourceMoments(source, radius)) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure moments for source."); + psFree(readout); + psFree(fakes); + return NULL; + } + } + + // Don't assume Poisson errors + pmPSFOptions *options = pmPSFOptionsAlloc(); // Options for fitting a PSF + options->poissonErrorsPhotLMM = true; + options->poissonErrorsPhotLin = false; + options->poissonErrorsParams = false; + options->stats = psStatsAlloc(PSF_STATS); + options->radius = radius; + options->psfTrendMode = PM_TREND_MAP; + options->psfTrendNx = xOrder; + options->psfTrendNy = yOrder; + options->psfFieldNx = numCols; + options->psfFieldNy = numRows; + options->psfFieldXo = 0; + options->psfFieldYo = 0; + + pmSourceFitModelInit (15, 0.01, WEIGHT_VAL, options->poissonErrorsPhotLMM); + + pmPSFtry *try = pmPSFtryModel(fakes, modelName, options, 0, 0xff); + psFree(options); + if (!try) { + psError(PS_ERR_UNKNOWN, false, "Unable to fit PSF model to PSF envelope."); + psFree(readout); + psFree(fakes); + return NULL; + } + + pmPSF *psf = psMemIncrRefCounter(try->psf); // Output PSF + psFree(try); + +#ifdef TESTING + { + // Need to translate peak flux --> integrated flux + pmModel *fakeModel = pmModelFromPSFforXY(psf, (float)numCols / 2.0, (float)numRows / 2.0, + 1.0); // Fake model, with central intensity of 1.0 + float flux0 = fakeModel->modelFlux(fakeModel->params); // Flux for central intensity of 1.0 + for (int i = 0; i < numFakes; i++) { + pmSource *source = fakes->data[i]; // Fake source + source->psfMag -= 2.5 * log10(flux0); + } + + pmReadout *generated = pmReadoutAlloc(NULL); // Generated image + pmReadoutFakeFromSources(generated, numCols, numRows, fakes, NULL, NULL, psf, NAN, radius, + false); + { + psFits *fits = psFitsOpen("psf_field_model.fits", "w"); + psFitsWriteImage(fits, NULL, generated->image, 0, NULL); + psFitsClose(fits); + } + psBinaryOp(generated->image, generated->image, "-", readout->image); + { + psFits *fits = psFitsOpen("psf_field_resid.fits", "w"); + psFitsWriteImage(fits, NULL, generated->image, 0, NULL); + psFitsClose(fits); + } + psFree(generated); + } +#endif + + psFree(readout); + + return psf; +} Index: /tags/pap_root_080214/psModules/src/imcombine/pmPSFEnvelope.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmPSFEnvelope.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmPSFEnvelope.h (revision 16477) @@ -0,0 +1,22 @@ +#ifndef PM_PSF_ENVELOPE_H +#define PM_PSF_ENVELOPE_H + +#include +#include +#include +#include +#include +#include + +/// Generate a PSF which is the envelope of an array of PSFs +/// +/// Generates multiple instances of the PSFs (distributed over an image) +pmPSF *pmPSFEnvelope(int numCols, int numRows, // Size of original image + const psArray *inputs, // Input PSF models + int instances, // Number of instances per dimension + int radius, // Radius of each PSF + const char *modelName, // Name of PSF model to use + int xOrder, int yOrder // Order for PSF variation + ); + +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmReadoutCombine.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmReadoutCombine.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmReadoutCombine.c (revision 16477) @@ -0,0 +1,366 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmHDUUtils.h" +#include "pmFPAMaskWeight.h" +#include "pmConceptsAverage.h" +#include "pmReadoutStack.h" + +#include "pmReadoutCombine.h" + +//#define SHOW_BUSY 1 // Show that the function is busy + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Allocator for pmCombineParams +pmCombineParams *pmCombineParamsAlloc(psStatsOptions combine) +{ + pmCombineParams *params = psAlloc(sizeof(pmCombineParams)); + + params->combine = combine; + params->maskVal = 0; + params->blank = 0; + params->nKeep = 0; + params->fracHigh = 0.0; + params->fracHigh = 0.0; + params->iter = 1; + params->rej = INFINITY; + params->weights = false; + + return params; +} + + +// XXX: Maybe add support for S16 and S32 types. Currently, only F32 supported. +bool pmReadoutCombine(pmReadout *output, const psArray *inputs, const psVector *zero, const psVector *scale, + const pmCombineParams *params) +{ + // Check inputs + PS_ASSERT_PTR_NON_NULL(output, false); + PS_ASSERT_ARRAY_NON_NULL(inputs, false); + PS_ASSERT_PTR_NON_NULL(params, false); + if (zero) { + PS_ASSERT_VECTOR_TYPE(zero, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_SIZE(zero, inputs->n, false); + } + if (scale) { + PS_ASSERT_VECTOR_TYPE(scale, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_SIZE(scale, inputs->n, false); + } + PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracLow, 0.0, 1.0, false); + PS_ASSERT_FLOAT_WITHIN_RANGE(params->fracHigh, 0.0, 1.0, false); + if (params->combine != PS_STAT_SAMPLE_MEAN && params->combine != PS_STAT_SAMPLE_MEDIAN && + params->combine != PS_STAT_ROBUST_MEDIAN && params->combine != PS_STAT_FITTED_MEAN && + params->combine != PS_STAT_CLIPPED_MEAN) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Combination method is not SAMPLE_MEAN, SAMPLE_MEDIAN, " + "ROBUST_MEDIAN, FITTED_MEAN or CLIPPED_MEAN.\n"); + return false; + } + for (int i = 0; i < inputs->n; i++) { + pmReadout *readout = inputs->data[i]; // Readout of interest + if (params->weights && !readout->weight) { + psError(PS_ERR_UNEXPECTED_NULL, true, + "Rejection based on weights requested, but no weights supplied for image %d.\n", i); + return false; + } + } + + bool first = !output->image; // First pass through? + + pmHDU *hdu = pmHDUFromReadout(output); // Output HDU + if (!hdu) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find HDU for readout.\n"); + return false; + } + + if (first) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Combining using statistic: %x", params->combine); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + + psStats *stats = psStatsAlloc(params->combine); // The statistics to use in the combination + if (params->combine == PS_STAT_CLIPPED_MEAN) { + stats->clipSigma = params->rej; + stats->clipIter = params->iter; + + if (first) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Combination clipping: %d iterations, rejection at %f sigma", + params->iter, params->rej); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + } + + int minInputCols, maxInputCols, minInputRows, maxInputRows; // Smallest and largest values to combine + int xSize, ySize; // Size of the output image + if (!pmReadoutStackValidate(&minInputCols, &maxInputCols, &minInputRows, &maxInputRows, &xSize, &ySize, + inputs)) { + psError(PS_ERR_UNKNOWN, false, "No valid input readouts."); + return false; + } + + pmReadoutUpdateSize(output, minInputCols, minInputRows, xSize, ySize, true); + psTrace("psModules.imcombine", 7, "Output minimum: %d,%d\n", output->col0, output->row0); + + psStatsOptions combineStdev = 0; // Statistics option for weights + if (params->weights) { + + if (!output->weight) { + output->weight = psImageAlloc(xSize, ySize, PS_TYPE_F32); + } + if (output->weight->numCols < xSize || output->weight->numRows < ySize) { + psImage *newWeight = psImageAlloc(xSize, ySize, PS_TYPE_F32); + psImageInit(newWeight, 0.0); + psImageOverlaySection(newWeight, output->weight, output->col0, output->row0, "="); + psFree(output->weight); + output->weight = newWeight; + } + + if (first) { + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, + "Using input weights to combine images", ""); + } + + // Get the correct statistics option for weights + switch (params->combine) { + case PS_STAT_SAMPLE_MEAN: + case PS_STAT_SAMPLE_MEDIAN: + combineStdev = PS_STAT_SAMPLE_STDEV; + break; + case PS_STAT_ROBUST_MEDIAN: + combineStdev = PS_STAT_ROBUST_STDEV; + break; + case PS_STAT_FITTED_MEAN: + combineStdev = PS_STAT_FITTED_STDEV; + break; + case PS_STAT_CLIPPED_MEAN: + combineStdev = PS_STAT_CLIPPED_STDEV; + break; + default: + psAbort("Should never get here --- checked params->combine before.\n"); + } + stats->options |= combineStdev; + } + + // We loop through each pixel in the output image. We loop through each input readout. We determine if + // that output pixel is contained in the image from that readout. If so, we save it in psVector pixels. + // If not, we set a mask for that element in pixels. Then, we mask off pixels not between fracLow and + // fracHigh. Then we call the vector stats routine on those pixels/mask. Then we set the output pixel + // value to the result of the stats call. + + psVector *pixels = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of pixels + psF32 *pixelsData = pixels->data.F32; // Dereference pixels + + psVector *mask = psVectorAlloc(inputs->n, PS_TYPE_U8); // Mask for stack + psU8 *maskData = mask->data.U8; // Dereference mask + + psVector *weights = NULL; // Stack of weights + psVector *errors = NULL; // Stack of errors (sqrt of variance/weights), for psVectorStats + psF32 *weightsData = NULL; // Dereference weights + if (params->weights) { + weights = psVectorAlloc(inputs->n, PS_TYPE_F32); // Stack of weights + weightsData = weights->data.F32; + } + psVector *index = NULL; // The indices to sort the pixels + + float keepFrac = 1.0 - params->fracLow - params->fracHigh; // Fraction of pixels to keep + if (keepFrac != 1.0 && first) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Min/max rejection: %f high, %f low, keep %d", + params->fracHigh, params->fracLow, params->nKeep); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + + psMaskType maskVal = params->maskVal; // The mask value + if (maskVal && first) { + psString comment = NULL; // Comment to add to header + psStringAppend(&comment, "Mask for combination: %x", maskVal); + psMetadataAddStr(hdu->header, PS_LIST_TAIL, "HISTORY", PS_META_DUPLICATE_OK, comment, ""); + psFree(comment); + } + + #ifndef PS_NO_TRACE + psTrace("psModules.imcombine", 3, "Iterating output: %d --> %d, %d --> %d\n", + minInputCols - output->col0, maxInputCols - output->col0, + minInputRows - output->row0, maxInputRows - output->row0); + if (psTraceGetLevel("psModules.imcombine") >= 3) { + for (int r = 0; r < inputs->n; r++) { + pmReadout *readout = inputs->data[r]; // Input readout + psTrace("psModules.imcombine", 3, "Iterating input %d: %d --> %d, %d --> %d\n", r, + minInputCols - readout->col0, maxInputCols - readout->col0, + minInputRows - readout->row0, maxInputRows - readout->row0); + } + } + #endif + + // Dereference output products + psF32 **outputImage = output->image->data.F32; // Output image + psU8 **outputMask = output->mask->data.U8; // Output mask + psF32 **outputWeight = NULL; // Output weight map + if (output->weight) { + outputWeight = output->weight->data.F32; + } + + psVector *invScale = NULL; // Inverse scale; pre-calculated for efficiency + if (scale) { + invScale = (psVector*)psBinaryOp(NULL, psScalarAlloc(1.0, PS_TYPE_F32), "/", (const psPtr)scale); + } + + for (int i = minInputRows; i < maxInputRows; i++) { + int yOut = i - output->row0; // y position on output readout + #ifdef SHOW_BUSY + + if (psTraceGetLevel("psModules.imcombine") > 9) { + printf("Processing row %d\r", i); + fflush(stdout); + } + #endif + for (int j = minInputCols; j < maxInputCols; j++) { + int xOut = j - output->col0; // x position on output readout + + int numValid = 0; // Number of valid pixels in the stack + memset(maskData, 0, mask->n * sizeof(psU8)); // Reset the mask + for (int r = 0; r < inputs->n; r++) { + pmReadout *readout = inputs->data[r]; // Input readout + int yIn = i - readout->row0; // y position on input readout + int xIn = j - readout->col0; // x position on input readout + psImage *image = readout->image; // The readout image + + #if 0 // This should have been taken care of already: + // Check bounds + if (xIn < 0 || xIn >= image->numCols || yIn < 0 || yIn >= image->numRows) { + continue; + } + #endif + + pixelsData[r] = image->data.F32[yIn][xIn]; + if (!isfinite(pixelsData[r])) { + maskData[r] = 1; + continue; + } + + // Check mask + psImage *roMask = readout->mask; // The mask image + if (roMask && roMask->data.U8[yIn][xIn] & maskVal) { + maskData[r] = 1; + continue; + } + + if (params->weights) { + weightsData[r] = readout->weight->data.F32[yIn][xIn]; + } + + if (zero) { + pixelsData[r] -= zero->data.F32[r]; + } + if (scale) { + pixelsData[r] *= invScale->data.F32[r]; + if (params->weights) { + weightsData[r] *= invScale->data.F32[r] * invScale->data.F32[r]; + } + } + + numValid++; + } + + if (numValid == 0) { + outputMask[yOut][xOut] = params->blank; + outputImage[yOut][xOut] = NAN; + continue; + } + + // Apply fracLow,fracHigh if there are enough pixels + if (numValid * keepFrac >= params->nKeep && keepFrac != 1.0) { + index = psVectorSortIndex(index, pixels); + int numLow = numValid * params->fracLow; // Number of low pixels to clip + int numHigh = numValid * params->fracHigh; // Number of high pixels to clip + // Low pixels + psS32 *indexData = index->data.S32; // Dereference index + for (int k = 0, numMasked = 0; numMasked < numLow && k < index->n; k++) { + // Don't count the ones that are already masked + if (!maskData[indexData[k]]) { + maskData[indexData[k]] = 1; + numMasked++; + } + } + // High pixels + for (int k = pixels->n - 1, numMasked = 0; numMasked < numHigh && k >= 0; k--) { + // Don't count the ones that are already masked + if (! maskData[indexData[k]]) { + maskData[indexData[k]] = 1; + numMasked++; + } + } + } + + // XXXXX this step probably is very expensive : convert errors to variance everywhere? + if (params->weights) { + errors = (psVector*)psUnaryOp(errors, weights, "sqrt"); + } + + // Combination + if (!psVectorStats(stats, pixels, errors, mask, 1)) { + // Can't do much about it, but it's not worth worrying about + psErrorClear(); + outputImage[yOut][xOut] = NAN; + outputMask[yOut][xOut] = params->blank; + if (params->weights) { + outputWeight[yOut][xOut] = NAN; + } + } else { + outputImage[yOut][xOut] = psStatsGetValue(stats, params->combine); + if (!isfinite(outputImage[yOut][xOut])) { + outputMask[yOut][xOut] = params->blank; + } + if (params->weights) { + float stdev = psStatsGetValue(stats, combineStdev); + outputWeight[yOut][xOut] = PS_SQR(stdev); // Variance + // XXXX this is not the correct formal error. + // also, the weighted mean is not obviously the correct thing here + } + } + } + } + #ifdef SHOW_BUSY + if (psTraceGetLevel("psModules.imcombine") > 9) { + printf("\n"); + } + #endif + psFree(index); + psFree(pixels); + psFree(mask); + psFree(weights); + psFree(errors); + psFree(stats); + psFree(invScale); + + // Update the "concepts" + psList *inputCells = psListAlloc(NULL); // List of cells + for (long i = 0; i < inputs->n; i++) { + pmReadout *readout = inputs->data[i]; // Readout of interest + psListAdd(inputCells, PS_LIST_TAIL, readout->parent); + } + bool success = pmConceptsAverageCells(output->parent, inputCells, NULL, NULL, true); + psFree(inputCells); + + output->data_exists = true; + output->parent->data_exists = true; + output->parent->parent->data_exists = true; + + return success; +} + Index: /tags/pap_root_080214/psModules/src/imcombine/pmReadoutCombine.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmReadoutCombine.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmReadoutCombine.h (revision 16477) @@ -0,0 +1,49 @@ +/* @file pmReadoutCombine.h + * @brief Combine multiple readouts + * + * @author George Gusciora, MHPCC + * @author Paul Price, IfA + * + * @version $Revision: 1.12 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-06-02 03:51:03 $ + * Copyright 2004-2006 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_READOUT_COMBINE_H +#define PM_READOUT_COMBINE_H + +/// @addtogroup imcombine Image Combinations +/// @{ + +/// Combination parameters for pmReadoutCombine. +/// +/// These values define how the combination is performed, and should not vary by detector, so that it can be +/// re-used for multiple combinations. +typedef struct +{ + psStatsOptions combine; ///< Statistic to use when performing the combination + psMaskType maskVal; ///< Mask value + psMaskType blank; ///< Mask value to give blank (i.e., no data) pixels + int nKeep; ///< Mimimum number of pixels to keep + float fracHigh; ///< Fraction of high pixels to immediately throw + float fracLow; ///< Fraction of low pixels to immediately throw + int iter; ///< Number of iterations for clipping (for CLIPPED_MEAN only) + float rej; ///< Rejection threshould for clipping (for CLIPPED_MEAN only) + bool weights; ///< Use the supplied weights (instead of calculated stdev)? +} +pmCombineParams; + +// Allocator for pmCombineParams +pmCombineParams *pmCombineParamsAlloc(psStatsOptions statsOptions ///< Statistic to use for combination + ); + +/// Combine multiple readouts, applying zero and scale, with optional minmax clipping +bool pmReadoutCombine(pmReadout *output,///< Output readout; altered and returned + const psArray *inputs, ///< Array of input readouts (F32 image and weight, U8 mask) + const psVector *zero, ///< Zero corrections to subtract from input, or NULL + const psVector *scale, ///< Scale corrections to divide into input, or NULL + const pmCombineParams *params ///< Combination parameters + ); + +/// @} +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmStack.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmStack.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmStack.c (revision 16477) @@ -0,0 +1,737 @@ +/** @file pmStack.c + * + * This file will perform image combination of several images of the + * same field, produce a list of questionable pixels, then tag some + * of those pixels as cosmic rays. + * + * @author Paul Price, IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.16 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-13 02:55:33 $ + * Copyright 2004-2007 Institute for Astronomy, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmConceptsAverage.h" + +#include "pmStack.h" + +#define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time +#define PIXEL_MAP_BUFFER 2 // Number of entries to add to pixel map at a time + +// Internal values for masks +#define MASK_BAD 0x01 // Mask value for pixels that have formerly been masked as bad +#define MASK_SUSPECT 0x02 // Mask value for pixels that are suspected bad (by pmStackCombine) + + +// Data structure for use as a buffer when combining pixels +// Use of this structure means we don't have to do an allocation in the combination function for each pixel +typedef struct { + psVector *pixels; // Pixel values + psVector *masks; // Pixel masks + psVector *weights; // Pixel weights + psVector *sort; // Buffer for sorting (to get a robust estimator of the standard dev) +} combineBuffer; + +static void combineBufferFree(combineBuffer *buffer) +{ + psFree(buffer->pixels); + psFree(buffer->masks); + psFree(buffer->weights); + psFree(buffer->sort); + return; +} + +static combineBuffer *combineBufferAlloc(long numImages, // Number of images that will be combined + psStatsOptions stat // Statistic to use + ) +{ + combineBuffer *buffer = psAlloc(sizeof(combineBuffer)); + psMemSetDeallocator(buffer, (psFreeFunc)combineBufferFree); + + buffer->pixels = psVectorAlloc(numImages, PS_TYPE_F32); + buffer->masks = psVectorAlloc(numImages, PS_TYPE_MASK); + buffer->weights = psVectorAlloc(numImages, PS_TYPE_F32); + buffer->sort = psVectorAlloc(numImages, PS_TYPE_F32); + + return buffer; +} + + +// Deallocator for the stack data +static void stackDataFree(pmStackData *data) +{ + psFree(data->readout); + psFree(data->pixels); + return; +} + + +// Generate a mean value for the combination +// Not using psVectorStats because it assumes that the weights are errors, and weights by 1/error^2 +static bool combinationMean(float *mean, // Mean value, to return + const psVector *values, // Values to combine + const psVector *weights, // Weights to apply + const psVector *masks, // Mask to apply + psMaskType maskVal // Value to mask + ) +{ + assert(values && weights && masks); + assert(values->n == weights->n); + assert(values->n == masks->n); + assert(values->type.type == PS_TYPE_F32); + assert(weights->type.type == PS_TYPE_F32); + assert(masks->type.type == PS_TYPE_MASK); + + float sumValueWeight = 0.0; // Sum of the value multiplied by the weight + float sumWeight = 0.0; // Sum of the weights + for (int i = 0; i < values->n; i++) { + if (!(masks->data.PS_TYPE_MASK_DATA[i] & maskVal)) { + sumValueWeight += values->data.F32[i] * weights->data.F32[i]; + sumWeight += weights->data.F32[i]; + } + } + if (sumWeight <= 0) { + return false; + } + + *mean = sumValueWeight / sumWeight; + return true; +} + +// Generate a standard deviation for the combination +// Not using psVectorStats because it has additional allocations which slow things down +static bool combinationStdev(float *median, // Median value, to return + float *stdev, // Standard deviation value, to return + const psVector *values, // Values to combine + const psVector *masks, // Mask to apply + psMaskType maskVal, // Value to mask + psVector *sortBuffer // Buffer for sorting + ) +{ + assert(values && masks); + assert(values->n == masks->n); + assert(values->type.type == PS_TYPE_F32); + assert(masks->type.type == PS_TYPE_MASK); + assert(sortBuffer && sortBuffer->nalloc >= values->n && sortBuffer->type.type == PS_TYPE_F32); + + int num = 0; // Number of valid values + for (int i = 0; i < values->n; i++) { + if (!(masks->data.PS_TYPE_MASK_DATA[i] & maskVal)) { + sortBuffer->data.F32[num++] = values->data.F32[i]; + } + } + sortBuffer->n = num; + if (!psVectorSortInPlace(sortBuffer)) { + return false; + } + *median = num % 2 ? (sortBuffer->data.F32[num / 2] + sortBuffer->data.F32[num / 2 + 1]) / 2.0 : + sortBuffer->data.F32[num / 2]; + *stdev = 0.74 * (sortBuffer->data.F32[(int)(0.75 * num)] - sortBuffer->data.F32[(int)(0.25 * num)]); + + return true; +} + +// Generate a variance value for the combination +static bool combinationVariance(float *variance, // Variance value, to return + const psVector *variances, // Pixel variances to combine + const psVector *weights, // Image weights to apply + const psVector *masks, // Mask to apply + psMaskType maskVal // Value to mask + ) +{ + assert(variances && weights && masks); + assert(variances->n == weights->n); + assert(variances->n == masks->n); + assert(variances->type.type == PS_TYPE_F32); + assert(weights->type.type == PS_TYPE_F32); + assert(masks->type.type == PS_TYPE_MASK); + + // Get the variance in the combination. We're not using the input pixel variances to generate a + // weighted average for the pixel flux (because that introduces systematic biases), so the variance + // of the output pixel value should simply be: + // simga^2 = sum(weight_i^2 * sigma_i^2) / (sum(weight_i))^2 + // This reduces, when the weights are all identically unity, to: + // variance_combination = sum(variance_i) / N^2 + // and if the variances are all equal: + // variance_combination = variance_individual / N + // which makes sense --- the standard deviation of the combination is reduced by a factor of sqrt(N). + + float sumWeights = 0.0; // Sum of the global weights + float sumVarianceWeights = 0.0; // Sum of the pixel variances multiplied by the global weights + for (int i = 0; i < variances->n; i++) { + if (!(masks->data.PS_TYPE_MASK_DATA[i] & maskVal)) { + sumVarianceWeights += variances->data.F32[i] * PS_SQR(weights->data.F32[i]); + sumWeights += weights->data.F32[i]; + } + } + + if (sumWeights <= 0) { + return false; + } + + *variance = sumVarianceWeights / PS_SQR(sumWeights); + return true; +} + +// Given a stack of images, combine with optional rejection. +// Pixels in the stack that are rejected are marked for subsequent +static bool combinePixels(psImage *image, // Combined image, for output + psImage *mask, // Combined mask, for output + psImage *weight, // Combined variance map, for output + const psArray *inputs, // Stack data + const psVector *weights, // Global (single value) weights for data, or NULL + const psVector *reject, // Indices of pixels to reject, or NULL + int x, int y, // Coordinates of interest + psMaskType maskVal, // Value to mask + psMaskType bad, // Value to give bad pixels + int numIter, // Number of rejection iterations + float rej, // Number of standard deviations at which to reject + combineBuffer *buffer // Buffer for combination; to avoid multiple allocations + ) +{ + // Rudimentary error checking + assert(image); + assert(mask); + assert(inputs); + assert(numIter >= 0); + assert(rej > 0); + assert(buffer); + + int num = inputs->n; // Number of images to combine + + psVector *pixelData = buffer->pixels; // Values for the pixel of interest + psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest + psVector *pixelWeights = buffer->weights; // Weights for the pixel of interest + psVector *sort = buffer->sort; // Sort buffer + + // Extract the pixel and mask data + int numBad = 0; // Number of bad pixels + for (int i = 0; i < num; i++) { + pmStackData *data = inputs->data[i]; // Stack data of interest + if (!data) { + continue; + } + psImage *image = data->readout->image; // Image of interest + psImage *weight = data->readout->weight; // Weight map of interest + psImage *mask = data->readout->mask; // Mask of interest + pixelData->data.F32[i] = image->data.F32[y][x]; + if (weight) { + pixelWeights->data.F32[i] = weight->data.F32[y][x]; + } + pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[y][x]; + if (pixelMasks->data.PS_TYPE_MASK_DATA[i] & maskVal) { + numBad++; + } + } + if (numBad == num) { + // Catch the case where everything is masked + image->data.F32[y][x] = NAN; + if (weight) { + weight->data.F32[y][x] = NAN; + } + mask->data.PS_TYPE_MASK_DATA[y][x] = bad; + return false; + } + + if (reject) { + for (int i = 0; i < reject->n; i++) { + pixelMasks->data.PS_TYPE_MASK_DATA[reject->data.U16[i]] |= maskVal; + } + } + + // Record the value derived with no clipping, because pixels rejected using the harsh clipping applied in + // the first pass might later be accepted. + float mean, variance; // Mean and variance of the combination + if (!combinationMean(&mean, pixelData, weights, pixelMasks, maskVal) || + (weight && !combinationVariance(&variance, pixelWeights, weights, pixelMasks, maskVal))) { + image->data.F32[y][x] = NAN; + mask->data.PS_TYPE_MASK_DATA[y][x] = bad; + if (weight) { + weight->data.F32[y][x] = NAN; + } + return false; + } + + image->data.F32[y][x] = mean; + if (weight) { + weight->data.F32[y][x] = variance; + } + mask->data.PS_TYPE_MASK_DATA[y][x] = 0; + + // The clipping that follows is solely to identify suspect pixels. + // These suspect pixels will be inspected in more detail by other functions. + long numClipped = LONG_MAX; // Number of pixels clipped + psMaskType ignore = maskVal | bad; // Ignore values with this mask value + for (int i = 0; i < numIter && numClipped > 0; i++) { + float median, stdev; // Median and stdev of the combination, for rejection + // Only get the mean if it's not the first iteration (because we got the mean earlier) + if (!combinationStdev(&median, &stdev, pixelData, pixelMasks, maskVal, sort)) { + image->data.F32[y][x] = NAN; + mask->data.PS_TYPE_MASK_DATA[y][x] = bad; + weight->data.F32[y][x] = NAN; + return false; + } + + float limit = rej * stdev; // Rejection limit + numClipped = 0; + for (int j = 0; j < num; j++) { + if (pixelMasks->data.PS_TYPE_MASK_DATA[j] & ignore) { + continue; + } + float diff = pixelData->data.F32[j] - median; + if (fabsf(diff) > limit) { + // Add the pixel as one to inspect + pmStackData *data = inputs->data[j]; // Stack data of interest + if (!data) { + continue; + } + data->pixels = psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y); + // Mask it so it's not considered in other iterations within this function + pixelMasks->data.PS_TYPE_MASK_DATA[j] |= bad; + numClipped++; + } + } + } + + return true; +} + + +// Ensure the input array of pmStackData is valid, and get some details out of it +static bool validateInputData(bool *haveWeights, // Do we have weights in the sky images? + bool *havePixels, // Do we have lists of pixels? + int *num, // Number of inputs + int *numCols, int *numRows, // Size of (sky) images + psArray *input // Input array of pmStackData to validate + ) +{ + PS_ASSERT_ARRAY_NON_NULL(input, false); + *num = input->n; + + pmStackData *data = NULL; // First image off the rank, used as a template + for (int i = 0; !data && i < input->n; i++) { + data = input->data[i]; + } + PS_ASSERT_PTR_NON_NULL(data, false); + assert(psMemGetDeallocator(data) == (psFreeFunc)stackDataFree); // Ensure it's the right type + *haveWeights = false; + PS_ASSERT_IMAGE_NON_NULL(data->readout->image, false); + PS_ASSERT_IMAGE_TYPE(data->readout->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(data->readout->mask, false); + PS_ASSERT_IMAGE_TYPE(data->readout->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->mask, false); + *numCols = data->readout->image->numCols; + *numRows = data->readout->image->numRows; + if (data->readout->weight) { + *haveWeights = true; + PS_ASSERT_IMAGE_NON_NULL(data->readout->weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->weight, false); + PS_ASSERT_IMAGE_TYPE(data->readout->weight, PS_TYPE_F32, false); + } + *havePixels = (data->pixels != NULL); + + // Make sure the rest correspond with the first + for (int i = 1; i < *num; i++) { + pmStackData *data = input->data[i]; // Stack data for this input + if (!data) { + continue; + } + assert(psMemGetDeallocator(data) == (psFreeFunc)stackDataFree); // Ensure it's the right type + if (!data->readout) { + psError(PS_ERR_UNEXPECTED_NULL, true, "The readout is specified in some but not all inputs."); + return false; + } + if ((*havePixels && !data->pixels) || (data->pixels && !*havePixels)) { + psError(PS_ERR_UNEXPECTED_NULL, true, "The pixels are specified in some but not all inputs."); + return false; + } + PS_ASSERT_IMAGE_NON_NULL(data->readout->image, false); + PS_ASSERT_IMAGE_NON_NULL(data->readout->mask, false); + PS_ASSERT_IMAGE_TYPE(data->readout->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_TYPE(data->readout->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGE_SIZE(data->readout->image, *numCols, *numRows, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->mask, false); + if (*haveWeights) { + PS_ASSERT_IMAGE_NON_NULL(data->readout->weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->readout->image, data->readout->weight, false); + PS_ASSERT_IMAGE_TYPE(data->readout->weight, PS_TYPE_F32, false); + } + } + + return true; +} + + +#if 0 +// Examine a pixel carefully to see if it should be rejected in the stack by convolving all the input images +// to the same seeing, and clipping there. +static bool convolveTest(psVector *reject, // Image indices to reject for this pixel (output) + psVector *inspect, // Images indices to inspect for this pixel (input) + psArray *input, // The input pmStackData + int xPix, int yPix, // Pixel coordinates of interest + psVector *seeing, // Seeing convolution for each image + psMaskType maskVal, // Value to mask + float extent, // Extent of convolution (Gaussian sigmas) + float threshold// Rejection threshold (standard deviations) + ) +{ + assert(reject && reject->type.type == PS_TYPE_U16); + assert(input); + assert(seeing); + assert(seeing->n == input->n); + + int num = input->n; // Number of input images + psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); // Statistics + psVector *values = psVectorAlloc(num, PS_TYPE_F32); // Values from convolving the image + psVector *valuesMask = psVectorAlloc(num, PS_TYPE_MASK); // Mask for convolution values + + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stacking data + int radius = extent * seeing->data.F32[i]; // How much to convolve + psImage *image = data->readout->image; // Image to convolve + psImage *mask = data->readout->mask; // Image mask + + int xMin = PS_MAX(xPix - radius, 0); + int xMax = PS_MIN(xPix + radius, image->numCols - 1); + int yMin = PS_MAX(yPix - radius, 0); + int yMax = PS_MIN(yPix + radius, image->numRows - 1); + + float sum = 0.0; // Sum of pixels in convolution + float sumKernel = 0.0; // Sum of kernel + for (int y = yMin; y <= yMax; y++) { + int v2 = PS_SQR(y - yPix); + for (int x = xMin; x <= xMax; x++) { + int u2 = PS_SQR(x - xPix); + float kernel = expf(-(u2 + v2) / seeing->data.F32[i]); + if (!(mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal)) { + sum += image->data.F32[y][x] * kernel; + sumKernel += kernel; + } + } + } + + if (sumKernel != 0.0) { + values->data.F32[i] = sum / sumKernel; + valuesMask->data.PS_TYPE_MASK_DATA[i] = 0; + } else { + // Everything's masked out; don't include this pixel in the analysis + values->data.F32[i] = NAN; + valuesMask->data.PS_TYPE_MASK_DATA[i] = MASK_BAD; + } + } + + // Mask the ones we're inspecting + for (int i = 0; i < inspect->n; i++) { + valuesMask->data.PS_TYPE_MASK_DATA[inspect->data.U16[i]] |= MASK_SUSPECT; + } + + if (!psVectorStats(stats, values, NULL, valuesMask, MASK_BAD | MASK_SUSPECT)) { + psFree(stats); + psFree(values); + psFree(valuesMask); + return false; + } + + // Inspect the pixels of interest + float mean = stats->sampleMean; // Mean value + float limit = threshold * stats->sampleStdev; // Rejection threshold + for (int j = 0; j < inspect->n; j++) { + int index = inspect->data.U16[j]; // Index of interest + if (fabsf(values->data.F32[index] - mean) > limit) { + reject->data.U16[reject->n++] = index; + } + } + + psFree(stats); + psFree(values); + psFree(valuesMask); + + return true; +} +#endif + + +// Generate a "pixel map". +// +// A "pixel map" is an image-like structure containing a vector that contains the indices of images. The idea +// is to provide a reverse lookup for an array of pixel lists, so that the image for which a pixel is flagged +// can be identified easily. +static psArray *pixelMapGenerate(const psArray *input, // Data to stack + int numCols, int numRows // Size of (sky) images + ) +{ + psArray *map = psArrayAlloc(numRows); // The pixel map + for (int y = 0; y < numRows; y++) { + map->data[y] = psArrayAlloc(numCols); + } + + for (int i = 0; i < input->n; i++) { + pmStackData *data = input->data[i]; + if (!data) { + continue; + } + assert(data->pixels); + psPixels *pixels = data->pixels;// The pixels of interest + for (int j = 0; j < pixels->n; j++) { + int x = pixels->data[j].x, y = pixels->data[j].y; // Coordinates of interest + if (x < 0 || x >= numCols || y < 0 || y >= numRows) { + psWarning("Bad pixel coordinate: %d,%d --- ignored.", x, y); + continue; + } + psArray *columns = map->data[y]; // The columns for that row + psVector *images = columns->data[x]; // The images for that column + if (!images) { + images = columns->data[x] = psVectorAllocEmpty(PIXEL_MAP_BUFFER, PS_TYPE_U16); + } + int size = images->n; // Element number at which to add + columns->data[x] = psVectorExtend(images, PIXEL_MAP_BUFFER, 1); + images->data.U16[size] = i; + } + } + + return map; +} + +// Query a "pixel map", by returning the list of image indices for a particular pixel. +static psVector *pixelMapQuery(const psArray *map, // Pixel map + int x, int y // Coordinates of interest + ) +{ + assert(y >= 0 && y < map->n); + psArray *colMap = map->data[y]; // Columns for that row + assert(x >= 0 && x < colMap->n); + return colMap->data[x]; +} + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +/// Constructor +pmStackData *pmStackDataAlloc(pmReadout *readout, float weight) +{ + pmStackData *data = psAlloc(sizeof(pmStackData)); // Stack data, to return + psMemSetDeallocator(data, (psFreeFunc)stackDataFree); + + data->readout = psMemIncrRefCounter(readout); + data->pixels = NULL; + data->weight = weight; + + return data; +} + +/// Stack input images +bool pmStackCombine(pmReadout *combined, psArray *input, psMaskType maskVal, psMaskType bad, + int numIter, float rej) +{ + PS_ASSERT_PTR_NON_NULL(combined, false); + bool haveWeights; // Do we have the weight maps? + bool havePixels; // Do we have lists of pixels? + int num; // Number of inputs + int numCols, numRows; // Size of (sky) images + if (!validateInputData(&haveWeights, &havePixels, &num, &numCols, &numRows, input)) { + return false; + } + PS_ASSERT_INT_POSITIVE(bad, false); + PS_ASSERT_INT_NONNEGATIVE(numIter, false); + PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); + if (havePixels) { + // This is a subsequent combination, so expect that the image and mask already exist + PS_ASSERT_IMAGE_NON_NULL(combined->image, false); + PS_ASSERT_IMAGE_TYPE(combined->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_NON_NULL(combined->mask, false); + PS_ASSERT_IMAGE_TYPE(combined->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(combined->image, combined->mask, false); + } + + // Pull out the weightings + psVector *weights = psVectorAlloc(num, PS_TYPE_F32); + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stack data for this input + if (!data) { + weights->data.F32[i] = 0.0; + } + weights->data.F32[i] = data->weight; + } + + // Buffer for combination + combineBuffer *buffer = combineBufferAlloc(num, numIter == 0 ? PS_STAT_SAMPLE_MEAN : + PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN); + + if (havePixels) { + // Only combine select pixels + psImage *combinedImage = combined->image; // Combined image + psImage *combinedMask = combined->mask; // Combined mask + psImage *combinedWeight = combined->weight; // Combined mask + + psArray *pixelMap = pixelMapGenerate(input, numCols, numRows); // Map of pixels to source + psPixels *pixels = NULL; // Total list of pixels, with no duplicates + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels + if (!data) { + continue; + } + pixels = psPixelsConcatenate(pixels, data->pixels); + data->pixels = psPixelsRealloc(data->pixels, PIXEL_LIST_BUFFER); // Just in case more rejection + } + for (int i = 0; i < pixels->n; i++) { + int x = pixels->data[i].x, y = pixels->data[i].y; // Coordinates of interest + psVector *reject = pixelMapQuery(pixelMap, x, y); // Inspect these images closely + combinePixels(combinedImage, combinedMask, combinedWeight, input, weights, reject, x, y, + maskVal, bad, numIter, rej, buffer); + } + psTrace("psModules.imcombine", 5, "Additional %ld pixels fixed.\n", pixels->n); + psFree(pixels); + psFree(pixelMap); + } else { + // Pull the products out, allocate if necessary + psImage *combinedImage = combined->image; // Combined image + if (!combinedImage) { + combined->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + combinedImage = combined->image; + } + psImage *combinedMask = combined->mask; // Combined mask + if (!combinedMask) { + combined->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); + combinedMask = combined->mask; + } + + psImage *combinedWeights = combined->weight; // Combined weight map + if (haveWeights && !combinedWeights) { + combined->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + combinedWeights = combined->weight; + } + + // Generate the pixel lists in which to place the rejected pixels + if (numIter != 0) { + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stack data for this input + if (!data) { + continue; + } + data->pixels = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); + } + } + + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + combinePixels(combinedImage, combinedMask, combinedWeights, input, weights, NULL, x, y, + maskVal, bad, numIter, rej, buffer); + } + } + +#ifndef PS_NO_TRACE + if (psTraceGetLevel("psModules.imcombine") >= 5) { + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stack data for this input + if (!data) { + continue; + } + psTrace("psModules.imcombine", 5, "Image %d: %ld pixels to inspect.\n", i, data->pixels->n); + } + } +#endif + } + + psFree(weights); + psFree(buffer); + + return true; +} + +#if 0 +bool pmStackReject(psArray *input, psMaskType maskVal, float extent, float threshold) +{ + bool haveWeights; // Do we have the sky weight maps? + bool havePixels; // Do we have lists of pixels? + int num; // Number of inputs + int numCols, numRows; // Size of (sky) images + if (!validateInputData(&haveWeights, &havePixels, &num, &numCols, &numRows, input)) { + return false; + } + if (!havePixels) { + psError(PS_ERR_UNEXPECTED_NULL, true, "No pixels have been flagged for possible rejection."); + return false; + } + PS_ASSERT_FLOAT_LARGER_THAN(extent, 0.0, false); + + // Want to concatenate the list of pixels, but also want to know where each came from + psArray *pixelMap = pixelMapGenerate(input, numCols, numRows); // Map of pixels to source + psPixels *pixels = NULL; // Total list of pixels, with no duplicates + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stacking data; contains the list of pixels + pixels = psPixelsConcatenate(pixels, data->pixels); + // Throw out the old in preparation for the new + data->pixels = psPixelsRealloc(data->pixels, 0); + } + + // Get the convolution widths + psVector *seeing = psVectorAlloc(num, PS_TYPE_F32); // Seeing for each image + const float fwhm2sigma = 1.0 / (2.0*sqrt((double)2.0*(log((double)2.0)))); // Conversion FWHM --> sigma + float seeingMax = -INFINITY; // Maximum FWHM + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stacking data + seeing->data.F32[i] = data->seeing * fwhm2sigma; + if (data->seeing > seeingMax) { + seeingMax = data->seeing; + } + } + (void)psBinaryOp(seeing, seeing, "*", seeing); + (void)psBinaryOp(seeing, psScalarAlloc(PS_SQR(seeingMax), PS_TYPE_F32), "-", seeing); + (void)psUnaryOp(seeing, seeing, "sqrt"); + + + psVector *reject = psVectorAllocEmpty(num, PS_TYPE_U16); // Image indices to reject for a given pixel + + for (int i = 0; i < pixels->n; i++) { + int x = pixels->data[i].x; // x coordinate of interest + int y = pixels->data[i].y; // y coordinate of interest + psVector *inspect = pixelMapQuery(pixelMap, x, y); // Inspect these images closely + + // Can daisy-chain multiple tests here + reject->n = 0; + + if (!convolveTest(reject, inspect, input, x, y, seeing, maskVal, extent, threshold)) { + // Some sort of error message here; NAN and mask all the pixels + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stacking data + psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y); + } + continue; + } + + // Add the rejected pixels to the official rejection list + for (int i = 0; i < reject->n; i++) { + pmStackData *data = input->data[reject->data.U16[i]]; // Stacking data + psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y); + } + } + psFree(reject); + psFree(seeing); + psFree(pixels); + psFree(pixelMap); + +#ifndef PS_NO_TRACE + if (psTraceGetLevel("psModules.imcombine") >= 5) { + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stack data for this input + psTrace("psModules.imcombine", 5, "Image %d: %ld pixels rejected.\n", i, data->pixels->n); + } + } +#endif + + return true; +} +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmStack.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmStack.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmStack.h (revision 16477) @@ -0,0 +1,83 @@ +/* @file pmStack.h + * + * This file will perform image combination of several images of the + * same field, produce a list of questionable pixels, then tag some + * of those pixels as defects. + * + * @author Paul Price, IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.4 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-08-23 23:43:12 $ + * + * Copyright 2004-2007 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_STACK_H +#define PM_STACK_H + +#include +#include +#include + +/// @addtogroup imcombine Image Combinations +/// @{ + + +/// Container for input image +typedef struct { + pmReadout *readout; ///< Warped readout (sky cell) + psPixels *pixels; ///< Pixels to inspect or reject + float weight; ///< Weight to apply +} pmStackData; + +/// Constructor +pmStackData *pmStackDataAlloc(pmReadout *readout, ///< Warped readout (sky cell) + float weight ///< Weight to apply + ); + +/// Stack input images +bool pmStackCombine(pmReadout *combined,///< Combined readout (output) + psArray *input, ///< Input array of pmStackData + psMaskType maskVal, ///< Mask value of bad pixels + psMaskType bad, ///< Mask value to give rejected pixels + int numIter, ///< Number of iterations + float rej ///< Rejection limit (standard deviations) + ); + +#if 0 +psImage *pmStackCombine(psImage *combine, ///< Combined image (output) + psArray **questionablePixels, ///< Array of rejection masks + const psArray *images, ///< Array of input images + const psArray *errors, ///< Array of input error images + const psArray *masks, ///< Array of input masks + psU32 maskVal, ///< Mask value + const psPixels *pixels, ///< Pixels to combine + psS32 numIter, ///< Number of rejection iterations + psF32 sigmaClip ///< Number of standard deviations at which to reject +); +#endif + +#if 0 +/// Reject pixles in input images +bool pmStackReject(psArray *input, ///< Input array of pmStackData + psMaskType maskVal, ///< Value to mask + float extent, ///< Gaussian convolution extent, in Gaussian sigmas + float threshold ///< Rejection threshold, in standard deviations + ); +#endif + +#if 0 +psArray *pmStackReject( + const psArray *images, ///< Array of input images + const psArray *masks, ///< Array of input image masks + const psArray *errors, ///< The pixels which were rejected in the combination + const psArray *inToOut, ///< Transformation from input to output system + const psArray *outToIn, ///< Transformation from output to input system + psF32 rejThreshold, ///< Rejection threshold + psF32 gradLimit ///< Gradient limit +); +#endif + +/// @} +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmStackReject.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmStackReject.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmStackReject.c (revision 16477) @@ -0,0 +1,160 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmSubtraction.h" +#include "pmSubtractionKernels.h" + +#define PIXEL_LIST_BUFFER 100 // Number of pixels to add to list at a time + +psPixels *pmStackReject(const psPixels *in, const psRegion *valid, float threshold, + const psArray *subRegions, const psArray *kernels) +{ + PS_ASSERT_PIXELS_NON_NULL(in, NULL); + PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(threshold, 0.0, NULL); + PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(threshold, 1.0, NULL); + PS_ASSERT_ARRAY_NON_NULL(subRegions, NULL); + PS_ASSERT_ARRAY_NON_NULL(kernels, NULL); + PS_ASSERT_ARRAYS_SIZE_EQUAL(subRegions, kernels, NULL); + + // Get the original image size + int numRegions = subRegions->n; // Number of regions + int numCols = 0, numRows = 0; // Size of original image + int minCols = INT_MAX, minRows = INT_MAX; // Minimum coordinate for image --- should be 0,0 + int size = 0; // Size of kernel + for (int i = 0; i < numRegions; i++) { + psRegion *subRegion = subRegions->data[i]; // Region of interest + if (subRegion->x0 < minCols) { + minCols = subRegion->x0; + } + if (subRegion->y0 < minRows) { + minRows = subRegion->y0; + } + if (subRegion->x1 > numCols) { + numCols = subRegion->x1; + } + if (subRegion->y1 > numRows) { + numRows = subRegion->y1; + } + + pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest + if (size == 0) { + size = kernel->size; + } else if (kernel->size != size) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Kernel sizes are not identical: %d vs %d", + size, kernel->size); + return NULL; + } + } + + // Adjust the size for the size of the subimage + if (valid) { + minCols = PS_MAX(valid->x0, minCols); + minRows = PS_MAX(valid->y0, minRows); + numCols = PS_MIN(valid->x1, numCols); + numRows = PS_MIN(valid->y1, numRows); + } + + psImage *mask = psPixelsToMask(NULL, in, psRegionSet(minCols, numCols - 1, minRows, numRows - 1), + 0x01); // Mask + psImage *image = psImageCopy(NULL, mask, PS_TYPE_F32); // Floating-point version, so we can convolve + psFree(mask); + + // Convolve the image with the kernel --- we're basically applying a matched filter and then thresholding + pmReadout *convRO = pmReadoutAlloc(NULL); // Readout with convolved image + pmReadout *inRO = pmReadoutAlloc(NULL); // Readout with input image + inRO->image = image; + inRO->col0 = minCols; + inRO->row0 = minRows; + for (int i = 0; i < numRegions; i++) { + psRegion *region = subRegions->data[i]; // Region of interest + if (valid && (region->x0 > valid->x1 || region->x1 < valid->x0 || + region->y0 > valid->y1 || region->y1 < valid->y0)) { + // Region is outside of our sub-image + continue; + } + pmSubtractionKernels *kernel = kernels->data[i]; // Kernel of interest + if (!pmSubtractionConvolve(convRO, NULL, inRO, NULL, NULL, 0, region, kernel, true)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve mask image in region %d.", i); + psFree(convRO); + psFree(inRO); + return NULL; + } + + // Need to adjust the thresholding level for the normalisation of the kernel --- the application of + // the kernel may scale the unit level that we've inserted. + + // Image of the kernel at the centre of the region + float xNorm = (region->x0 + 0.5 * (region->x1 - region->x0) - kernel->numCols/2.0) / + (float)kernel->numCols; + float yNorm = (region->y0 + 0.5 * (region->y1 - region->y0) - kernel->numRows/2.0) / + (float)kernel->numRows; + psImage *image = pmSubtractionKernelImage(kernel, xNorm, yNorm, false); + if (!image) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); + psFree(convRO); + psFree(inRO); + return NULL; + } + float sum = 0.0; + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + sum += image->data.F32[y][x]; + } + } + psFree(image); + + // Range for normalisation + int yMin = PS_MAX(minRows, region->y0) - inRO->row0; + int yMax = PS_MIN(numRows - 1, region->y1) - inRO->row0; + int xMin = PS_MAX(minCols, region->x0) - inRO->col0; + int xMax = PS_MIN(numCols - 1, region->x1) - inRO->col0; + psTrace("psModules.imcombine", 2, "Normalising convolved mask image by %f over %d:%d,%d:%d\n", + sum, xMin, xMax, yMin, yMax); + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + convRO->image->data.F32[y][x] /= sum; + } + } + } + psFree(inRO); + psImage *convolved = psMemIncrRefCounter(convRO->image); + psFree(convRO); + + // Threshold the convolved image + psPixels *bad = psPixelsAllocEmpty(PIXEL_LIST_BUFFER); // List of pixels that should be masked + for (int y = 0; y < convolved->numRows; y++) { + for (int x = 0; x < convolved->numCols; x++) { + if (convolved->data.F32[y][x] > threshold) { + bad = psPixelsAdd(bad, PIXEL_LIST_BUFFER, x, y); + } + } + } + psFree(convolved); + + // Now, grow the mask to include everything that touches a bad pixel in the convolution + mask = psPixelsToMask(NULL, bad, psRegionSet(0, numCols - 1, 0, numRows - 1), 0xff); + assert(mask->numCols == numCols && mask->numRows == numRows); + for (int i = 0; i < bad->n; i++) { + int xPix = bad->data[i].x, yPix = bad->data[i].y; // Coordinates of interest + // Convolution limits + int xMin = PS_MAX(xPix - size, 0); + int xMax = PS_MIN(xPix + size, numCols - 1); + int yMin = PS_MAX(yPix - size, 0); + int yMax = PS_MIN(yPix + size, numRows - 1); + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + assert(x < mask->numCols && y < mask->numRows); + mask->data.PS_TYPE_MASK_DATA[y][x] = 0xff; + } + } + } + + bad = psPixelsFromMask(bad, mask, 0xff); + psFree(mask); + + return bad; +} Index: /tags/pap_root_080214/psModules/src/imcombine/pmStackReject.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmStackReject.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmStackReject.h (revision 16477) @@ -0,0 +1,19 @@ +#ifndef PM_STACK_REJECT_H +#define PM_STACK_REJECT_H + +#include +#include + +/// Given a list of pixels from the convolved image, find the corresponding (smaller subset of) pixels in the +/// original image, and then convolve those to get the list of all pixels which should be rejected +/// +/// We apply a matched filter to the corresponding mask image, and threshold to find the original pixels +psPixels *pmStackReject(const psPixels *in, ///< List of pixels in the convolved image + const psRegion *valid, ///< Valid region to consider + float threshold, ///< Threshold on convolved image, 0..1 + const psArray *regions, ///< Array of image regions for image + const psArray *kernels ///< Array of kernel parameters for each region + ); + + +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtraction.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtraction.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtraction.c (revision 16477) @@ -0,0 +1,840 @@ +/** @file pmSubtraction.c + * + * @author Paul Price, IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.79 $ $Name: not supported by cvs2svn $ + * @date $Date: 2008-02-14 23:18:34 $ + * + * Copyright 2004-2007 Institute for Astronomy, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmHDU.h" // Required for pmFPA.h +#include "pmFPA.h" +#include "pmSubtractionStamps.h" +#include "pmSubtractionEquation.h" + +#include "pmSubtraction.h" + +//#define TESTING + +#define PIXEL_LIST_BUFFER 100 // Number of entries to add to pixel list at a time + + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private (file-static) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Generate the kernel to apply to the variance from the normal kernel +static psKernel *varianceKernel(psKernel *out, // Output kernel + psKernel *normalKernel // Normal kernel + ) +{ + // Kernel range + int xMin = normalKernel->xMin, xMax = normalKernel->xMax; + int yMin = normalKernel->yMin, yMax = normalKernel->yMax; + + if (!out) { + out = psKernelAlloc(xMin, xMax, yMin, yMax); + } + + // Take the square of the normal kernel + for (int v = yMin; v <= yMax; v++) { + for (int u = xMin; u <= xMax; u++) { + out->kernel[v][u] = PS_SQR(normalKernel->kernel[v][u]); + } + } + + return out; +} + +// Generate an image of the solved kernel +static psKernel *solvedKernel(psKernel *kernel, // Kernel, to return + const pmSubtractionKernels *kernels, // Kernel basis functions + const psImage *polyValues, // Spatial polynomial values + bool wantDual // Want the dual (second) kernel? + ) +{ + assert(kernels); + assert(polyValues); + + int numKernels = kernels->num; // Number of kernel basis functions + int size = kernels->size; // Kernel half-size + if (!kernel) { + kernel = psKernelAlloc(-size, size, -size, size); + } + psImageInit(kernel->image, 0.0); + + for (int i = 0; i < numKernels; i++) { + double value = p_pmSubtractionSolutionCoeff(kernels, polyValues, i, wantDual); // Polynomial value + + switch (kernels->type) { + case PM_SUBTRACTION_KERNEL_POIS: { + int u = kernels->u->data.S32[i]; // Offset in x + int v = kernels->v->data.S32[i]; // Offset in y + kernel->kernel[v][u] += value; + kernel->kernel[0][0] -= value; + break; + } + /* SPAM and FRIES use the same method */ + case PM_SUBTRACTION_KERNEL_SPAM: + case PM_SUBTRACTION_KERNEL_FRIES: { + int uStart = kernels->u->data.S32[i]; + int uStop = kernels->uStop->data.S32[i]; + int vStart = kernels->v->data.S32[i]; + int vStop = kernels->vStop->data.S32[i]; + + // Normalising sum of kernel component to unity + float norm = 1.0 / (float)((uStop - uStart + 1) * (vStop - vStart + 1)); + + for (int v = vStart; v <= vStop; v++) { + for (int u = uStart; u <= uStop; u++) { + kernel->kernel[v][u] += norm * value; + kernel->kernel[0][0] -= value; + } + } + break; + } + case PM_SUBTRACTION_KERNEL_GUNK: { + if (i < kernels->inner) { + // Using pre-calculated function + psKernel *preCalc = kernels->preCalc->data[i]; // Precalculated values + // Iterating over the kernel + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + kernel->kernel[v][u] += preCalc->kernel[v][u] * value; + // Photometric scaling is built into the preCalc kernel --- no subtraction! + } + } + } else { + // Using delta function + int u = kernels->u->data.S32[i]; // Offset in x + int v = kernels->v->data.S32[i]; // Offset in y + kernel->kernel[v][u] += value; + kernel->kernel[0][0] -= value; + } + break; + } + case PM_SUBTRACTION_KERNEL_ISIS: { + psKernel *preCalc = kernels->preCalc->data[i]; // Precalculated values + // Iterating over the kernel + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + kernel->kernel[v][u] += preCalc->kernel[v][u] * value; + // Photometric scaling is built into the preCalc kernel --- no subtraction! + } + } + break; + } + case PM_SUBTRACTION_KERNEL_RINGS: { + psArray *preCalc = kernels->preCalc->data[i]; // Precalculated data + psVector *uCoords = preCalc->data[0]; // u coordinates + psVector *vCoords = preCalc->data[1]; // v coordinates + psVector *poly = preCalc->data[2]; // Polynomial values + int num = uCoords->n; // Number of pixels + + for (int j = 0; j < num; j++) { + int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates + kernel->kernel[v][u] += poly->data.F32[j] * value; + } + // Photometric scaling is built into the kernel --- no subtraction! + break; + } + default: + psAbort("Should never get here."); + } + } + + // Put in the normalisation component + kernel->kernel[0][0] += (wantDual ? 1.0 : p_pmSubtractionSolutionNorm(kernels)); + + return kernel; +} + +// Subtract the (0,0) element to preserve photometric scaling +static void convolveSub(psKernel *convolved, // Convolved image + const psKernel *image, // Image being convolved + int footprint // Size of region of interest + ) +{ + // Can't use psBinaryOp because the images are of different size + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + convolved->kernel[y][x] -= image->kernel[y][x]; + } + } + return; +} + +// Generate the convolution given some offset +static psKernel *convolveOffset(const psKernel *image, // Image to convolve (a kernel for convenience) + int u, int v, // Offset to apply + int footprint // Size of region of interest + ) +{ + psKernel *convolved = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Convolved image + int numBytes = (2 * footprint + 1) * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy + for (int y = -footprint; y <= footprint; y++) { + // Convolution with a delta function is just the value specified by the offset + memcpy(&convolved->kernel[y][-footprint], &image->kernel[y - v][-footprint - u], numBytes); + } + return convolved; +} + + +// Convolve an image using FFT +static void convolveFFT(psImage *target,// Place the result in here + const psImage *image, // Image to convolve + const psImage *mask, // Mask image + psMaskType maskVal, // Value to mask + const psKernel *kernel, // Kernel by which to convolve + psRegion region,// Region of interest + float background, // Background to add + int size // Size of (square) kernel + ) +{ + psRegion border = psRegionSet(region.x0 - size, region.x1 + size, + region.y0 - size, region.y1 + size); // Add a border + + // Casting away const so psImageSubset can add the child + psImage *subImage = psImageSubset((psImage*)image, border); // Subimage to convolve + psImage *subMask = mask ? psImageSubset((psImage*)mask, border) : NULL; // Subimage mask + psImage *convolved = psImageConvolveFFT(subImage, subMask, maskVal, kernel, 0.0); // Convolution + psFree(subImage); + psFree(subMask); + + // Now, we have to chop off the borders, and stick it in where it belongs + psImage *subConv = psImageSubset(convolved, psRegionSet(size, -size, size, -size)); // Cut off the edges + psImage *subTarget = psImageSubset(target, region); // Target for this subregion + if (background != 0.0) { + psBinaryOp(subTarget, subConv, "+", psScalarAlloc(background, PS_TYPE_F32)); + } else { + int numBytes = subTarget->numCols * PSELEMTYPE_SIZEOF(PS_TYPE_F32); // Number of bytes to copy + for (int y = 0; y < subTarget->numRows; y++) { + memcpy(subTarget->data.F32[y], subConv->data.F32[y], numBytes); + } + } + psFree(subConv); + psFree(convolved); + psFree(subTarget); + return; +} + +// Convolve an image directly +static void convolveDirect(psImage *target, // Put the result here + const psImage *image, // Image to convolve + const psKernel *kernel, // Kernel by which to convolve + psRegion region,// Region of interest + float background, // Background to add + int size // Size of (square) kernel + ) +{ + for (int y = region.y0; y < region.y1; y++) { + for (int x = region.x0; x < region.x1; x++) { + target->data.F32[y][x] = background; + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + target->data.F32[y][x] += kernel->kernel[v][u] * + image->data.F32[y - v][x - u]; + } + } + } + } + return; +} + +// Convolve a region of an image +static inline void convolveRegion(psImage *convImage, // Convolved image (output) + psImage *convWeight, // Convolved weight map (output), or NULL + psKernel **kernelImage, // Convolution kernel for the image + psKernel **kernelWeight, // Convolution kernel for the weight map, or NULL + const psImage *image, // Image to convolve + const psImage *weight, // Weight map to convolve, or NULL + const psImage *subMask, // Subtraction mask + psMaskType maskVal, // Mask value to apply in convolution + const pmSubtractionKernels *kernels, // Kernels + psImage *polyValues, // Polynomial values + float background, // Background value to apply + psRegion region, // Region to convolve + bool useFFT, // Use FFT to convolve? + bool wantDual // Want the dual convolution? + ) +{ + *kernelImage = solvedKernel(*kernelImage, kernels, polyValues, wantDual); + if (weight) { + *kernelWeight = varianceKernel(*kernelWeight, *kernelImage); + } + + if (useFFT) { + // Use Fast Fourier Transform to do the convolution + // This provides a big speed-up for large kernels + convolveFFT(convImage, image, subMask, maskVal, *kernelImage, region, background, kernels->size); + if (weight) { + convolveFFT(convWeight, weight, subMask, maskVal, *kernelWeight, region, 0.0, kernels->size); + } + } else { + convolveDirect(convImage, image, *kernelImage, region, background, kernels->size); + if (weight) { + convolveDirect(convWeight, weight, *kernelWeight, region, 0.0, kernels->size); + } + } + + return; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Semi-public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, const psKernel *kernel) +{ + PS_ASSERT_KERNEL_NON_NULL(image, NULL); + PS_ASSERT_KERNEL_NON_NULL(kernel, NULL); + + psImage *conv = psImageConvolveFFT(image->image, NULL, 0, kernel, 0.0); // Convolved image + int x0 = - image->xMin, y0 = - image->yMin; // Position of centre of convolved image + psKernel *convolved = psKernelAllocFromImage(conv, x0, y0); // Kernel version + psFree(conv); + return convolved; +} + +psImage *p_pmSubtractionPolynomial(psImage *output, int spatialOrder, float x, float y) +{ + assert(spatialOrder >= 0); + assert(x >= -1 && x <= 1); + assert(y >= -1 && y <= 1); + + output = psImageRecycle(output, spatialOrder + 1, spatialOrder + 1, PS_TYPE_F64); + output->data.F64[0][0] = 1.0; + + double value = 1.0; + for (int i = 1; i <= spatialOrder; i++) { + value *= x; + output->data.F64[0][i] = value; + } + + value = 1.0; + for (int j = 1; j <= spatialOrder; j++) { + value *= y; + output->data.F64[j][0] = value; + } + + for (int j = 1; j <= spatialOrder; j++) { + for (int i = 1; i <= spatialOrder - j; i++) { + output->data.F64[j][i] = output->data.F64[j][0] * output->data.F64[0][i]; + } + } + + return output; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + +// Convolve a stamp by a single kernel basis function +static psKernel *convolveStampSingle(const pmSubtractionKernels *kernels, // Kernel basis functions + int index, // Kernel basis function index + const psKernel *image, // Image to convolve (a kernel for convenience) + int footprint // Size of region of interest + ) +{ +#if 0 + if (index == kernels->num) { + // Normalisation component + return convolveOffset(image, 0, 0, footprint); + } +#endif + + switch (kernels->type) { + case PM_SUBTRACTION_KERNEL_POIS: { + int u = kernels->u->data.S32[index]; // Offset in x + int v = kernels->v->data.S32[index]; // Offset in y + psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image + convolveSub(convolved, image, footprint); + return convolved; + } + // Method for SPAM and FRIES is the same + case PM_SUBTRACTION_KERNEL_SPAM: + case PM_SUBTRACTION_KERNEL_FRIES: { + psKernel *convolved = psKernelAlloc(-footprint, footprint, + -footprint, footprint); // Convolved image + int uStart = kernels->u->data.S32[index]; + int uStop = kernels->uStop->data.S32[index]; + int vStart = kernels->v->data.S32[index]; + int vStop = kernels->vStop->data.S32[index]; + float norm = 1.0 / (uStop - uStart + 1) * (vStop - vStart + 1); // Normalisation + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + double sum = 0.0; + for (int v = vStart; v <= vStop; v++) { + for (int u = uStart; u <= uStop; u++) { + sum += image->kernel[y - v][x - u]; + } + } + convolved->kernel[y][x] = norm * sum; + } + } + convolveSub(convolved, image, footprint); + return convolved; + } + case PM_SUBTRACTION_KERNEL_GUNK: { + if (index < kernels->inner) { + // Photometric scaling is already built in to the precalculated kernel + return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]); + } + // Using delta function + int u = kernels->u->data.S32[index]; // Offset in x + int v = kernels->v->data.S32[index]; // Offset in y + psKernel *convolved = convolveOffset(image, u, v, footprint); // Convolved image + convolveSub(convolved, image, footprint); + return convolved; + } + case PM_SUBTRACTION_KERNEL_ISIS: { + // Photometric scaling is already built in to the precalculated kernel + return p_pmSubtractionConvolveStampPrecalc(image, kernels->preCalc->data[index]); + } + case PM_SUBTRACTION_KERNEL_RINGS: { + psKernel *convolved = psKernelAlloc(-footprint, footprint, + -footprint, footprint); // Convolved image + psArray *preCalc = kernels->preCalc->data[index]; // Precalculated data + psVector *uCoords = preCalc->data[0]; // u coordinates + psVector *vCoords = preCalc->data[1]; // v coordinates + psVector *poly = preCalc->data[2]; // Polynomial values + int num = uCoords->n; // Number of pixels + for (int y = -footprint; y <= footprint; y++) { + for (int x = -footprint; x <= footprint; x++) { + double sum = 0.0; // Accumulated sum from convolution + for (int j = 0; j < num; j++) { + int u = uCoords->data.S32[j], v = vCoords->data.S32[j]; // Kernel coordinates + sum += image->kernel[y - v][x - u] * poly->data.F32[j]; + } + convolved->kernel[y][x] = sum; + // Photometric scaling is built into the kernel --- no subtraction! + } + } + return convolved; + } + default: + psAbort("Should never get here."); + } + return NULL; +} + +// Convolve the stamp by each of the kernel basis functions +static psArray *convolveStamp(psArray *convolutions, // The convolutions + const psKernel *image, // Image to convolve + const pmSubtractionKernels *kernels, // Kernel basis functions + int footprint // Stamp half-size + ) +{ + assert(image); + assert(kernels); + assert(footprint >= 0); + + if (convolutions) { + // Already done + return convolutions; + } + + int numKernels = kernels->num; // Number of kernels + convolutions = psArrayAlloc(numKernels); + + for (int i = 0; i < numKernels; i++) { + convolutions->data[i] = convolveStampSingle(kernels, i, image, footprint); + } + + return convolutions; +} + + +bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, const pmSubtractionKernels *kernels, int footprint) +{ + PS_ASSERT_PTR_NON_NULL(stamp, false); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + PS_ASSERT_INT_NONNEGATIVE(footprint, false); + + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Stamp not marked for calculation."); + return false; + } + + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_TARGET: + case PM_SUBTRACTION_MODE_1: + stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint); + break; + case PM_SUBTRACTION_MODE_2: + stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint); + break; + case PM_SUBTRACTION_MODE_UNSURE: + case PM_SUBTRACTION_MODE_DUAL: + stamp->convolutions1 = convolveStamp(stamp->convolutions1, stamp->image1, kernels, footprint); + stamp->convolutions2 = convolveStamp(stamp->convolutions2, stamp->image2, kernels, footprint); + break; + default: + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + return true; +} + + + + +int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, const psVector *deviations, + psImage *subMask, float sigmaRej, int footprint) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, -1); + PS_ASSERT_VECTOR_NON_NULL(deviations, -1); + PS_ASSERT_VECTOR_TYPE(deviations, PS_TYPE_F32, -1); + PS_ASSERT_IMAGE_NON_EMPTY(subMask, -1); + PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, -1); + + double totalSquareDev = 0.0; // Total square deviation from zero + int numStamps = 0; // Number of used stamps + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + continue; + } + totalSquareDev += PS_SQR(deviations->data.F32[i]); + numStamps++; + } + + if (numStamps == 0) { + psError(PS_ERR_UNKNOWN, true, "No good stamps found."); + return -1; + } + + if (!isfinite(sigmaRej) || sigmaRej <= 0.0) { + // User just wanted to calculate and record the RMS for posterity + psLogMsg("psModules.imcombine", PS_LOG_INFO, "RMS deviation: %f", sqrt(totalSquareDev / numStamps)); + return 0; + } + + float limit = sigmaRej * sqrt(totalSquareDev / (double)numStamps); // Limit on maximum deviation + psTrace("psModules.imcombine", 1, "Deviation limit: %f\n", limit); + + int numRejected = 0; // Number of stamps rejected + int numGood = 0; // Number of good stamps + double newSquareDev = 0.0; // New square deviation + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + if (deviations->data.F32[i] > limit) { + // Mask out the stamp in the image so you it's not found again + psTrace("psModules.imcombine", 3, "Rejecting stamp %d (%d,%d)\n", i, + (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); + numRejected++; + for (int y = stamp->y - footprint; y <= stamp->y + footprint; y++) { + for (int x = stamp->x - footprint; x <= stamp->x + footprint; x++) { + subMask->data.PS_TYPE_MASK_DATA[y][x] |= PM_SUBTRACTION_MASK_REJ; + } + } + + // Set stamp for replacement + stamp->x = 0; + stamp->y = 0; + stamp->xNorm = NAN; + stamp->yNorm = NAN; + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + // Recalculate convolutions + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + stamp->convolutions1 = stamp->convolutions2 = NULL; + psFree(stamp->image1); + psFree(stamp->image2); + psFree(stamp->weight); + stamp->image1 = stamp->image2 = stamp->weight = NULL; + psFree(stamp->matrix1); + psFree(stamp->matrix2); + psFree(stamp->matrixX); + stamp->matrix1 = stamp->matrix2 = stamp->matrixX = NULL; + psFree(stamp->vector1); + psFree(stamp->vector2); + stamp->vector1 = stamp->vector2 = NULL; + } else { + numGood++; + newSquareDev += PS_SQR(deviations->data.F32[i]); + } + } + } + + if (numRejected > 0) { + psLogMsg("psModules.imcombine", PS_LOG_INFO, + "%d good stamps; %d rejected.\nRMS deviation: %f --> %f\n", + numGood, numRejected, sqrt(totalSquareDev / numStamps), + sqrt(newSquareDev / (double)numGood)); + } else { + psLogMsg("psModules.imcombine", PS_LOG_INFO, + "%d good stamps; 0 rejected.\nRMS deviation: %f\n", + numGood, sqrt(totalSquareDev / numStamps)); + } + + return numRejected; +} + +psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); + + // Precalulate polynomial values + psImage *polyValues = p_pmSubtractionPolynomial(NULL, kernels->spatialOrder, x, y); + + // The appropriate kernel + psKernel *kernel = solvedKernel(NULL, kernels, polyValues, wantDual); + + psFree(polyValues); + + psImage *image = psMemIncrRefCounter(kernel->image); // Image of the kernel + psFree(kernel); + + return image; +} + +#if 0 +psArray *pmSubtractionKernelSolutions(const pmSubtractionKernels *kernels, float x, float y, bool wantDual) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(x, -1.0, 1.0, NULL); + PS_ASSERT_FLOAT_WITHIN_RANGE(y, -1.0, 1.0, NULL); + + psArray *images = psArrayAlloc(solution->n - 1); // Images of each kernel to return + psVector *fakeSolution = psVectorAlloc(solution->n, PS_TYPE_F64); // Fake solution vector + psVectorInit(fakeSolution, 0.0); + + for (int i = 0; i < solution->n - 1; i++) { + fakeSolution->data.F64[i] = solution->data.F64[i]; + images->data[i] = pmSubtractionKernelImage(kernels, x, y, wantDual); + fakeSolution->data.F64[i] = 0.0; + } + + psFree(fakeSolution); + + return images; +} +#endif + + + +bool pmSubtractionConvolve(pmReadout *out1, pmReadout *out2, const pmReadout *ro1, const pmReadout *ro2, + const psImage *subMask, psMaskType blank, const psRegion *region, + const pmSubtractionKernels *kernels, bool useFFT) +{ + PS_ASSERT_PTR_NON_NULL(out1, false); + PS_ASSERT_PTR_NON_NULL(ro1, false); + PS_ASSERT_IMAGE_NON_NULL(ro1->image, false); + PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false); + if (ro1->weight) { + PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false); + PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false); + } + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, false); + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + PS_ASSERT_PTR_NON_NULL(out2, false); + PS_ASSERT_PTR_NON_NULL(ro2, false); + PS_ASSERT_IMAGE_NON_NULL(ro2->image, false); + PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false); + if (ro2->weight) { + PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false); + PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false); + } + } + if (subMask) { + PS_ASSERT_IMAGE_NON_NULL(subMask, false); + PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(subMask, ro1->image, false); + } + if (out1->image) { + PS_ASSERT_IMAGE_NON_NULL(out1->image, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(out1->image, ro1->image, false); + PS_ASSERT_IMAGE_TYPE(out1->image, PS_TYPE_F32, false); + } + if (out1->mask) { + PS_ASSERT_IMAGE_NON_NULL(out1->mask, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(out1->mask, ro1->image, false); + PS_ASSERT_IMAGE_TYPE(out1->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGE_NON_NULL(subMask, false); + } + if (out1->weight) { + PS_ASSERT_IMAGE_NON_NULL(out1->weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(out1->weight, ro1->image, false); + PS_ASSERT_IMAGE_TYPE(out1->weight, PS_TYPE_F32, false); + } + if (region && psRegionIsNaN(*region)) { + psString string = psRegionToString(*region); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string); + psFree(string); + return false; + } + + psImage *image, *weight; // Image and weight map to convolve + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_TARGET: + case PM_SUBTRACTION_MODE_1: + case PM_SUBTRACTION_MODE_DUAL: + image = ro1->image; + weight = ro1->weight; + break; + case PM_SUBTRACTION_MODE_2: + image = ro2->image; + weight = ro2->weight; + break; + default: + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions + int x0 = ro1->col0, y0 = ro1->row0; // Image offset + + psImage *convImage1 = out1->image; // Convolved image + if (!convImage1) { + convImage1 = out1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImage *convMask = NULL; // Convolved mask image + if (subMask) { + if (!out1->mask) { + out1->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); + } + convMask = out1->mask; + psImageInit(convMask, 0); + } + psImage *convWeight1 = NULL; // Convolved weight (variance) image + if (weight) { + if (!out1->weight) { + out1->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + convWeight1 = out1->weight; + psImageInit(convWeight1, 0.0); + } + + psImage *convImage2 = NULL; + psImage *convWeight2 = NULL; // Convolved products for dual mode + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + convImage2 = out2->image; + if (!convImage2) { + convImage2 = out2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + convWeight2 = NULL; // Convolved weight (variance) image + if (ro2->weight) { + if (!out2->weight) { + out2->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + convWeight2 = out2->weight; + psImageInit(convWeight2, 0.0); + } + if (subMask) { + // Copying mask --- they're the same! + psFree(out2->mask); + out2->mask = psMemIncrRefCounter(convMask); + } + } + + int size = kernels->size; // Half-size of kernel + int fullSize = 2 * size + 1; // Full size of kernel + + // Get region for convolution: [xMin:xMax,yMin:yMax] + int xMin = size, xMax = numCols - size; + int yMin = size, yMax = numRows - size; + if (region) { + xMin = PS_MAX(region->x0, xMin); + xMax = PS_MIN(region->x1, xMax); + yMin = PS_MAX(region->y0, yMin); + yMax = PS_MIN(region->y1, yMax); + } + + // Size to use when calculating normalised coordinates (different from actual size when convolving + // subimage) + int xNormSize = (kernels->numCols > 0 ? kernels->numCols : numCols); + int yNormSize = (kernels->numRows > 0 ? kernels->numRows : numRows); + + psMaskType maskSource; // Mask these pixels when convolving + psMaskType maskTarget; // Mark these pixels as bad when propagating the subtractionMask + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_TARGET: + case PM_SUBTRACTION_MODE_1: + maskSource = PM_SUBTRACTION_MASK_BAD_1; + maskTarget = PM_SUBTRACTION_MASK_BAD_2 | PM_SUBTRACTION_MASK_CONVOLVE_1; + break; + case PM_SUBTRACTION_MODE_2: + maskSource = PM_SUBTRACTION_MASK_BAD_2; + maskTarget = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_CONVOLVE_2; + break; + case PM_SUBTRACTION_MODE_DUAL: + maskSource = PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2; + maskTarget = PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_CONVOLVE_2; + break; + default: + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + psImage *polyValues = NULL; // Pre-calculated polynomial values + psKernel *kernelImage = NULL; // Kernel for the images + psKernel *kernelWeight = NULL; // Kernel for the weight maps + + for (int j = yMin; j < yMax; j += fullSize) { + int ySubMax = PS_MIN(j + fullSize, yMax); // Range for subregion of interest + float yNorm = 2.0 * (float)(j + y0 + size + 1 - yNormSize/2.0) / + (float)yNormSize; // Normalised coordinate + for (int i = xMin; i < xMax; i += fullSize) { + int xSubMax = PS_MIN(i + fullSize, xMax); // Range for subregion of interest + float xNorm = 2.0 * (float)(i + x0 + size + 1 - xNormSize/2.0) / + (float)xNormSize; // Normalised coordinate + + // Only generate polynomial values every kernel footprint, since we have already assumed + // (with the stamps) that it does not vary rapidly on this scale. + polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, xNorm, yNorm); + float background = p_pmSubtractionSolutionBackground(kernels, polyValues); // Background term + + psRegion subRegion = psRegionSet(i, xSubMax, j, ySubMax); // Sub-region to convolve + + convolveRegion(convImage1, convWeight1, &kernelImage, &kernelWeight, image, weight, + subMask, maskSource, kernels, polyValues, background, subRegion, useFFT, false); + + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + convolveRegion(convImage2, convWeight2, &kernelImage, &kernelWeight, ro2->image, ro2->weight, + subMask, maskSource, kernels, polyValues, background, subRegion, useFFT, true); + } + + // Propagate the mask + if (subMask) { + for (int y = j; y < ySubMax; y++) { + for (int x = i; x < xSubMax; x++) { + if (subMask->data.PS_TYPE_MASK_DATA[y][x] & maskTarget) { + convMask->data.PS_TYPE_MASK_DATA[y][x] |= blank; + convImage1->data.F32[y][x] = NAN; + if (weight) { + convWeight1->data.F32[y][x] = NAN; + } + } + } + } + } + + } + } + psFree(kernelImage); + psFree(kernelWeight); + psFree(polyValues); + + return true; +} Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtraction.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtraction.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtraction.h (revision 16477) @@ -0,0 +1,111 @@ +/* @file pmSubtraction.h + * + * PSF-matched image subtraction, based on the Alard & Lupton (1998) and Alard (2000) methods. + * + * @author Paul Price, IfA + * @author GLG, MHPCC + * + * @version $Revision: 1.22 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-12-07 01:57:15 $ + * Copyright 2004-207 Institute for Astronomy, University of Hawaii + */ + +#ifndef PM_SUBTRACTION_H +#define PM_SUBTRACTION_H + +#include + +#include +#include +#include +#include + +/// @addtogroup imcombine Image Combinations +/// @{ + +/// Mask values for the subtraction mask +typedef enum { + PM_SUBTRACTION_MASK_CLEAR = 0x00, // No masking + PM_SUBTRACTION_MASK_BAD_1 = 0x01, // Image 1 is bad + PM_SUBTRACTION_MASK_BAD_2 = 0x02, // Image 2 is bad + PM_SUBTRACTION_MASK_CONVOLVE_1 = 0x04, // If image 1 is convolved, would be bad + PM_SUBTRACTION_MASK_CONVOLVE_2 = 0x08, // If image 2 is convolved, would be bad + PM_SUBTRACTION_MASK_FOOTPRINT_1 = 0x10, // Bad pixel within the stamp footprint of image 1 + PM_SUBTRACTION_MASK_FOOTPRINT_2 = 0x20, // Bad pixel within the stamp footprint of image 2 + PM_SUBTRACTION_MASK_BORDER = 0x40, // Image border + PM_SUBTRACTION_MASK_REJ = 0x80, // Previously tried as a stamp, and rejected +} pmSubtractionMasks; + + +/// Number of terms in a polynomial +#define PM_SUBTRACTION_POLYTERMS(ORDER) (((ORDER) + 1) * ((ORDER) + 2) / 2) + +/// Set the indices for the normalisation and background terms +#define PM_SUBTRACTION_INDICES(NORM,BG,KERNELS) { \ + int numSpatial = PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder); /* Number of spatial terms */ \ + NORM = (KERNELS)->num * numSpatial; \ + BG = NORM + 1; \ +} + +/// Return the index for the start of the normalisation terms +#define PM_SUBTRACTION_INDEX_NORM(KERNELS) \ + ((KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder)) + +/// Return the index for the start of the background terms +#define PM_SUBTRACTION_INDEX_BG(KERNELS) \ + (((KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder)) + 1) + + +/// Convolve the reference stamp with the kernel components +bool pmSubtractionConvolveStamp(pmSubtractionStamp *stamp, ///< Stamp to convolve + const pmSubtractionKernels *kernels, ///< Kernel parameters + int footprint ///< Half-size of region over which to calculate equation + ); + +/// Reject stamps +int pmSubtractionRejectStamps(pmSubtractionStampList *stamps, ///< Stamps + const psVector *deviations, ///< Deviations for each stamp + psImage *subMask, ///< Subtraction mask + float sigmaRej, ///< Number of RMS deviations above zero at which to reject + int footprint ///< Half-size of stamp + ); + +/// Generate an image of the convolution kernel +psImage *pmSubtractionKernelImage(const pmSubtractionKernels *kernels, ///< Kernel parameters + float x, float y,///< Normalised position [-1,1] for which to generate image + bool wantDual ///< Calculate for the dual kernel? + ); + +/// Generate images of the convolution kernel elements +psArray *pmSubtractionKernelSolutions(const psVector *solution, ///< Solution vector + const pmSubtractionKernels *kernels, ///< Kernel parameters + float x, float y ///< Normalised position [-1,1] for images + ); + +/// Convolve image in preparation for subtraction +bool pmSubtractionConvolve(pmReadout *out1, ///< Output image 1 + pmReadout *out2, ///< Output image 2 (DUAL mode only) + const pmReadout *ro1, // Input image 1 + const pmReadout *ro2, // Input image 2 + const psImage *subMask, ///< Subtraction mask (or NULL) + psMaskType blank, ///< Mask value for blank regions + const psRegion *region, ///< Region to convolve (or NULL) + const pmSubtractionKernels *kernels, ///< Kernel parameters + bool useFFT ///< Use Fast Fourier Transform for the convolution? + ); + +/// Generate the convolution of an image, given a precalculated kernel +/// +/// The 'image' is a kernel for convenience --- intended to be a stamp +psKernel *p_pmSubtractionConvolveStampPrecalc(const psKernel *image, ///< Image to convolve + const psKernel *kernel ///< Kernel by which to convolve + ); + +/// Given (normalised) coordinates (x,y), generate a matrix where the elements (i,j) are x^i * y^j +psImage *p_pmSubtractionPolynomial(psImage *output, ///< Output matrix, or NULL + int spatialOrder, ///< Maximum spatial polynomial order + float x, float y ///< Normalised position of interest, [-1,1] + ); + +/// @} +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionEquation.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionEquation.c (revision 16477) @@ -0,0 +1,1052 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#include "pmSubtraction.h" +#include "pmSubtractionKernels.h" +#include "pmSubtractionStamps.h" + +#include "pmSubtractionEquation.h" + + +//#define TESTING + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private (file-static) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Calculate the sum over a stamp product +static inline double calculateSumProduct(const psKernel *image1, // First image in multiplication + const psKernel *image2, // Second image in multiplication + const psKernel *weight, // Weight image + int footprint // (Half-)Size of stamp + ) +{ + double sum = 0.0; // Sum of the image products + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + sum += image1->kernel[y][x] * image2->kernel[y][x] / weight->kernel[y][x]; + } + } + return sum; +} + +// Calculate a single element of the least-squares matrix, with the polynomial expansions in one direction +static inline bool calculateMatrixElement1(psImage *matrix, // Matrix to calculate + int i, int j, // Coordinates of element + const psKernel *image1, // First image in multiplication + const psKernel *image2, // Second image in multiplication + const psKernel *weight, // Weight image + const psImage *polyValues, // Spatial polynomial values + int numKernels, // Number of kernel basis functions + int footprint, // (Half-)Size of stamp + int spatialOrder, // Maximum order of spatial variation + bool symmetric // Is the matrix symmetric? + ) +{ + double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products + if (!isfinite(sum)) { + return false; + } + + // Generate the pseudo-convolutions from the spatial polynomial terms + for (int iyOrder = 0, iIndex = i; iyOrder <= spatialOrder; iyOrder++) { + for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex += numKernels) { + double convPoly = sum * polyValues->data.F64[iyOrder][ixOrder]; + + assert(iIndex < matrix->numRows && j < matrix->numCols); + + matrix->data.F64[iIndex][j] = convPoly; + if (symmetric) { + + assert(iIndex < matrix->numCols && j < matrix->numRows); + + matrix->data.F64[j][iIndex] = convPoly; + } + } + } + return true; +} + +// Calculate a single element of the least-squares matrix, with the polynomial expansions in both directions +static inline bool calculateMatrixElement2(psImage *matrix, // Matrix to calculate + int i, int j, // Coordinates of element + const psKernel *image1, // First image in multiplication + const psKernel *image2, // Second image in multiplication + const psKernel *weight, // Weight image + const psImage *polyValues, // Spatial polynomial values + int numKernels, // Number of kernel basis functions + int footprint, // (Half-)Size of stamp + int spatialOrder, // Maximum order of spatial variation + bool symmetric // Is the matrix symmetric? + ) +{ + double sum = calculateSumProduct(image1, image2, weight, footprint); // Sum of the image products + if (!isfinite(sum)) { + return false; + } + + // Generate the pseudo-convolutions from the spatial polynomial terms + for (int iyOrder = 0, iIndex = i; iyOrder <= spatialOrder; iyOrder++) { + for (int ixOrder = 0; ixOrder <= spatialOrder - iyOrder; ixOrder++, iIndex += numKernels) { + double iPoly = polyValues->data.F64[iyOrder][ixOrder]; // Value of polynomial + for (int jyOrder = 0, jIndex = j; jyOrder <= spatialOrder; jyOrder++) { + for (int jxOrder = 0; jxOrder <= spatialOrder - jyOrder; jxOrder++, jIndex += numKernels) { + double convPoly = sum * iPoly * polyValues->data.F64[jyOrder][jxOrder]; + + assert(iIndex < matrix->numRows && jIndex < matrix->numCols); + + matrix->data.F64[iIndex][jIndex] = convPoly; + if (symmetric) { + + assert(iIndex < matrix->numCols && jIndex < matrix->numRows); + + matrix->data.F64[jIndex][iIndex] = convPoly; + } + } + } + } + } + return true; +} + +// Calculate the square part of the matrix derived from multiplying convolutions +static bool calculateMatrixSquare(psImage *matrix, // Matrix to calculate + const psArray *convolutions1, // Convolutions for element 1 + const psArray *convolutions2, // Convolutions for element 2 + const psKernel *weight, // Weight image + const psImage *polyValues, // Polynomial values + int numKernels, // Number of kernel basis functions + int spatialOrder, // Order of spatial variation + int footprint // Half-size of stamp + ) +{ + bool symmetric = (convolutions1 == convolutions2 ? true : false); // Is matrix symmetric? + + for (int i = 0; i < numKernels; i++) { + psKernel *iConv = convolutions1->data[i]; // Convolution for i-th element + + for (int j = (symmetric ? i : 0); j < numKernels; j++) { + psKernel *jConv = convolutions2->data[j]; // Convolution for j-th element + + if (!calculateMatrixElement2(matrix, i, j, iConv, jConv, weight, polyValues, numKernels, + footprint, spatialOrder, symmetric)) { + psTrace("psModules.imcombine", 2, "Bad sumCC at %d, %d", i, j); + return false; + } + } + } + + return true; +} + +// Calculate least-squares matrix and vector +static bool calculateMatrix(psImage *matrix, // Matrix to calculate + const pmSubtractionKernels *kernels, // Kernel components + const psArray *convolutions, // Convolutions of source with kernels + const psKernel *input, // Input stamp, or NULL + const psKernel *weight, // Weight stamp + const psImage *polyValues, // Spatial polynomial values + int footprint, // (Half-)Size of stamp + bool normAndBG // Calculate normalisation and background terms? + ) +{ + int numKernels = kernels->num; // Number of kernel components + int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variation terms + int bgOrder = kernels->bgOrder; // Maximum order of background fit + int numBackground = normAndBG ? PM_SUBTRACTION_POLYTERMS(bgOrder) : 0; // Number of background terms + int numTerms = numKernels * numSpatial + (normAndBG ? 1 + numBackground : 0); // Total number of terms + assert(matrix); + assert(matrix->numCols == matrix->numRows); + assert(matrix->numCols == numTerms); + assert(convolutions && convolutions->n == numKernels); + assert(polyValues); + assert(!normAndBG || input); // If we want the normalisation and BG, then we need the input image + + // Square part of the matrix (convolution-convolution products) + if (!calculateMatrixSquare(matrix, convolutions, convolutions, weight, polyValues, numKernels, + spatialOrder, footprint)) { + return false; + } + + // XXX To support higher-order background model than simply constant, the below code needs to be updated. + if (normAndBG) { + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index in matrix for background + + for (int i = 0; i < numKernels; i++) { + psKernel *conv = convolutions->data[i]; // Convolution for i-th element + + // Normalisation-convolution terms + if (!calculateMatrixElement1(matrix, i, normIndex, conv, input, weight, polyValues, numKernels, + footprint, spatialOrder, true)) { + psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); + return false; + } + + // Background-convolution terms + double sumC = 0.0; // Sum of the convolution + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + sumC += conv->kernel[y][x] / weight->kernel[y][x]; + } + } + if (!isfinite(sumC)) { + psTrace("psModules.imcombine", 2, "Bad sumC at %d", i); + return false; + } + + for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { + for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { + double value = sumC * polyValues->data.F64[yOrder][xOrder]; + matrix->data.F64[index][bgIndex] = value; + matrix->data.F64[bgIndex][index] = value; + } + } + } + + // Background only, normalisation only, and background-normalisation terms + double sum1 = 0.0; // Sum of the weighting + double sumI = 0.0; // Sum of the input + double sumII = 0.0; // Sum of the input squared + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double invNoise2 = 1.0 / weight->kernel[y][x]; + double value = input->kernel[y][x] * invNoise2; + sumI += value; + sumII += value * input->kernel[y][x]; + sum1 += invNoise2; + } + } + if (!isfinite(sumI)) { + psTrace("psModules.imcombine", 2, "Bad sumI detected"); + return false; + } + if (!isfinite(sumII)) { + psTrace("psModules.imcombine", 2, "Bad sumII detected"); + return false; + } + if (!isfinite(sum1)) { + psTrace("psModules.imcombine", 2, "Bad sum1 detected"); + return false; + } + matrix->data.F64[normIndex][normIndex] = sumII; + matrix->data.F64[bgIndex][bgIndex] = sum1; + matrix->data.F64[normIndex][bgIndex] = sumI; + matrix->data.F64[bgIndex][normIndex] = sumI; + } + + return true; +} + + +// Calculate least-squares matrix and vector +static bool calculateVector(psVector *vector, // Vector to calculate, or NULL + const pmSubtractionKernels *kernels, // Kernel components + const psArray *convolutions, // Convolutions of source with kernels + const psKernel *input, // Input stamp, or NULL if !normAndBG + const psKernel *target, // Target stamp + const psKernel *weight, // Weight stamp + const psImage *polyValues, // Spatial polynomial values + int footprint, // (Half-)Size of stamp + bool normAndBG // Calculate normalisation and background terms? + ) +{ + int numKernels = kernels->num; // Number of kernel components + int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variation terms + int bgOrder = kernels->bgOrder; // Maximum order of background fit + int numBackground = normAndBG ? PM_SUBTRACTION_POLYTERMS(bgOrder) : 0; // Number of background terms + int numTerms = numKernels * numSpatial + (normAndBG ? 1 + numBackground : 0); // Total number of terms + assert(vector && vector->n == numTerms); + assert(convolutions && convolutions->n == numKernels); + assert(target); + assert(polyValues); + assert(!normAndBG || input); // If we want the normalisation and BG, then we need the input image + + // Convolution terms + for (int i = 0; i < numKernels; i++) { + psKernel *conv = convolutions->data[i]; // Convolution for i-th element + double sumTC = 0.0; // Sum of the target and convolution + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + sumTC += target->kernel[y][x] * conv->kernel[y][x] / weight->kernel[y][x]; + } + } + if (!isfinite(sumTC)) { + psTrace("psModules.imcombine", 2, "Bad sumTC at %d", i); + return false; + } + for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { + for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { + vector->data.F64[index] = sumTC * polyValues->data.F64[yOrder][xOrder]; + } + } + } + + if (normAndBG) { + // Background terms + double sumT = 0.0; // Sum of the target + double sumIT = 0.0; // Sum of the input-target product + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + float value = target->kernel[y][x] / weight->kernel[y][x]; + sumIT += value * input->kernel[y][x]; + sumT += value; + } + } + if (!isfinite(sumT)) { + psTrace("psModules.imcombine", 2, "Bad sumI detected"); + return false; + } + if (!isfinite(sumIT)) { + psTrace("psModules.imcombine", 2, "Bad sumIT detected"); + return false; + } + + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation term + vector->data.F64[normIndex] = sumIT; + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background term + vector->data.F64[bgIndex] = sumT; + } + + return true; +} + + + +// Calculate the cross-matrix, composed of convolutions of each image +// Note that the cross-matrix is NOT square +static bool calculateMatrixCross(psImage *matrix, // Matrix to calculate + const pmSubtractionKernels *kernels, // Kernel components + const psArray *convolutions1, // Convolutions of image 1 + const psArray *convolutions2, // Convolutions of image 2 + const psKernel *image1, // Image 1 stamp + const psKernel *weight, // Weight stamp + const psImage *polyValues, // Spatial polynomial values + int footprint // (Half-)Size of stamp + ) +{ + assert(matrix); + int numKernels = kernels->num; // Number of kernel components + int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial polynomial terms + int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + int numCols = numKernels * numSpatial + 1 + numBackground; // Number of columns + int numRows = numKernels * numSpatial; // Number of rows + assert(matrix->numCols == numCols && matrix->numRows == numRows); + assert(convolutions1 && convolutions1->n == numKernels); + assert(convolutions2 && convolutions2->n == numKernels); + + int normIndex, bgIndex; // Indices in matrix for normalisation and background terms + PM_SUBTRACTION_INDICES(normIndex, bgIndex, kernels); + + if (!calculateMatrixSquare(matrix, convolutions1, convolutions2, weight, polyValues, numKernels, + spatialOrder, footprint)) { + return false; + } + + for (int i = 0; i < numKernels; i++) { + // Normalisation + psKernel *conv = convolutions2->data[i]; // Convolution + if (!calculateMatrixElement1(matrix, i, normIndex, conv, image1, weight, polyValues, numKernels, + footprint, spatialOrder, false)) { + psTrace("psModules.imcombine", 2, "Bad sumIC at %d", i); + return false; + } + + // Background + double sumC = 0.0; // Sum of the weighting + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + sumC += conv->kernel[y][x] / weight->kernel[y][x]; + } + } + if (!isfinite(sumC)) { + psTrace("psModules.imcombine", 2, "Bad sumC detected at %d", i); + return false; + } + for (int yOrder = 0, index = i; yOrder <= spatialOrder; yOrder++) { + for (int xOrder = 0; xOrder <= spatialOrder - yOrder; xOrder++, index += numKernels) { + matrix->data.F64[index][bgIndex] = sumC * polyValues->data.F64[yOrder][xOrder]; + } + } + } + + return true; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Semi-public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Calculate the value of a polynomial, specified by coefficients and polynomial values +inline double p_pmSubtractionCalculatePolynomial(const psVector *coeff, // Coefficients + const psImage *polyValues, // Polynomial values + int order, // Order of polynomials + int index, // Index at which to begin + int step // Step between subsequent indices + ) +{ + double sum = 0.0; // Value of the polynomial sum + for (int yOrder = 0; yOrder <= order; yOrder++) { + for (int xOrder = 0; xOrder <= order - yOrder; xOrder++, index += step) { + + assert(index < coeff->n); + + sum += coeff->data.F64[index] * polyValues->data.F64[yOrder][xOrder]; + } + } + return sum; +} + +inline double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, const psImage *polyValues, + int index, bool wantDual) +{ +#if 0 + // This is probably in a tight loop, so don't check inputs + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); + PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); + PS_ASSERT_INT_POSITIVE(index, NAN); +#endif + + psVector *solution = wantDual ? kernels->solution2 : kernels->solution1; // Solution vector + return p_pmSubtractionCalculatePolynomial(solution, polyValues, kernels->spatialOrder, index, + kernels->num); +} + +inline double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels) +{ +#if 0 + // This is probably in a tight loop, so don't check inputs + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); + PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); +#endif + + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation + return kernels->solution1->data.F64[normIndex]; +} + +inline double p_pmSubtractionSolutionBackground(const pmSubtractionKernels *kernels, + const psImage *polyValues) +{ +#if 0 + // This is probably in a tight loop, so don't check inputs + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NAN); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NAN); + PS_ASSERT_IMAGE_NON_NULL(polyValues, NAN); +#endif + + int bgIndex = PM_SUBTRACTION_INDEX_BG(kernels); // Index for background + return p_pmSubtractionCalculatePolynomial(kernels->solution1, polyValues, kernels->bgOrder, bgIndex, 1); +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, const pmSubtractionKernels *kernels) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + + int footprint = stamps->footprint; // Half-size of stamps + int spatialOrder = kernels->spatialOrder; // Maximum order of spatial variation + int numKernels = kernels->num; // Number of kernel basis functions + int numSpatial = PM_SUBTRACTION_POLYTERMS(spatialOrder); // Number of spatial variations + int numBackground = PM_SUBTRACTION_POLYTERMS(kernels->bgOrder); // Number of background terms + + // Total number of parameters to solve for: coefficient of each kernel basis function, multipled by the + // number of coefficients for the spatial polynomial, normalisation and a constant background offset. + int numParams = numKernels * numSpatial + 1 + numBackground; + + psImage *polyValues = NULL; // Polynomial terms + + // We iterate over each stamp, allocate the matrix and vectors if + // necessary, and then calculate those matrix/vectors. + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE) { + continue; + } + + // Generate convolutions + if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i); + psFree(polyValues); + return NULL; + } + +#ifdef TESTING + for (int j = 0; j < numKernels; j++) { + if (stamp->convolutions1) { + psString convName = NULL; + psStringAppend(&convName, "conv1_%03d_%03d.fits", i, j); + psFits *fits = psFitsOpen(convName, "w"); + psFree(convName); + psKernel *conv = stamp->convolutions1->data[j]; + psFitsWriteImage(fits, NULL, conv->image, 0, NULL); + psFitsClose(fits); + } + + if (stamp->convolutions2) { + psString convName = NULL; + psStringAppend(&convName, "conv2_%03d_%03d.fits", i, j); + psFits *fits = psFitsOpen(convName, "w"); + psFree(convName); + psKernel *conv = stamp->convolutions2->data[j]; + psFitsWriteImage(fits, NULL, conv->image, 0, NULL); + psFitsClose(fits); + } + } +#endif + + polyValues = p_pmSubtractionPolynomial(polyValues, spatialOrder, stamp->xNorm, stamp->yNorm); + + + stamp->matrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); + stamp->vector1 = psVectorAlloc(numParams, PS_TYPE_F64); +#ifdef TESTING + psImageInit(stamp->matrix1, NAN); + psVectorInit(stamp->vector1, NAN); +#endif + + bool status; // Status of least-squares matrix/vector calculation + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_TARGET: + case PM_SUBTRACTION_MODE_1: + status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, + stamp->weight, polyValues, footprint, true); + status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, + stamp->image2, stamp->weight, polyValues, footprint, true); + break; + case PM_SUBTRACTION_MODE_2: + status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions2, stamp->image2, + stamp->weight, polyValues, footprint, true); + status &= calculateVector(stamp->vector1, kernels, stamp->convolutions2, stamp->image2, + stamp->image1, stamp->weight, polyValues, footprint, true); + break; + case PM_SUBTRACTION_MODE_DUAL: + stamp->matrix2 = psImageAlloc(numKernels * numSpatial, numKernels * numSpatial, PS_TYPE_F64); + stamp->matrixX = psImageAlloc(numParams, numKernels * numSpatial, PS_TYPE_F64); + stamp->vector2 = psVectorAlloc(numKernels * numSpatial, PS_TYPE_F64); +#ifdef TESTING + psImageInit(stamp->matrix2, NAN); + psImageInit(stamp->matrixX, NAN); + psVectorInit(stamp->vector2, NAN); +#endif + status = calculateMatrix(stamp->matrix1, kernels, stamp->convolutions1, stamp->image1, + stamp->weight, polyValues, footprint, true); + status &= calculateMatrix(stamp->matrix2, kernels, stamp->convolutions2, NULL, + stamp->weight, polyValues, footprint, false); + status &= calculateMatrixCross(stamp->matrixX, kernels, stamp->convolutions1, + stamp->convolutions2, stamp->image1, stamp->weight, polyValues, + footprint); + status &= calculateVector(stamp->vector1, kernels, stamp->convolutions1, stamp->image1, + stamp->image2, stamp->weight, polyValues, footprint, true); + status &= calculateVector(stamp->vector2, kernels, stamp->convolutions2, NULL, + stamp->image2, stamp->weight, polyValues, footprint, false); + break; + default: + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + if (!status) { + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + psWarning("Rejecting stamp %d (%d,%d) because of bad equation", + i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); + } else { + stamp->status = PM_SUBTRACTION_STAMP_USED; + } + +#ifdef TESTING + if (psTraceGetLevel("psModules.imcombine.equation") >= 10) { + psString matrixName = NULL; + psStringAppend(&matrixName, "matrix1_%d.fits", i); + psFits *matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, stamp->matrix1, 0, NULL); + psFitsClose(matrixFile); + + matrixName = NULL; + psStringAppend(&matrixName, "vector1_%d.fits", i); + psImage *dummy = psImageAlloc(stamp->vector1->n, 1, PS_TYPE_F64); + memcpy(dummy->data.F64[0], stamp->vector1->data.F64, + PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector1->n); + matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL); + psFree(dummy); + psFitsClose(matrixFile); + + if (stamp->vector2) { + matrixName = NULL; + psStringAppend(&matrixName, "vector2_%d.fits", i); + dummy = psImageAlloc(stamp->vector2->n, 1, PS_TYPE_F64); + memcpy(dummy->data.F64[0], stamp->vector2->data.F64, + PSELEMTYPE_SIZEOF(PS_TYPE_F64) * stamp->vector2->n); + matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, dummy, 0, NULL); + psFree(dummy); + psFitsClose(matrixFile); + } + + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + matrixName = NULL; + psStringAppend(&matrixName, "matrix2_%d.fits", i); + matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, stamp->matrix2, 0, NULL); + psFitsClose(matrixFile); + + matrixName = NULL; + psStringAppend(&matrixName, "matrixX_%d.fits", i); + matrixFile = psFitsOpen(matrixName, "w"); + psFree(matrixName); + psFitsWriteImage(matrixFile, NULL, stamp->matrixX, 0, NULL); + psFitsClose(matrixFile); + } + } +#endif + + } + psFree(polyValues); + + return true; +} + +bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, const pmSubtractionStampList *stamps) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + + // Check inputs + int numParams = -1; // Number of parameters + int numParams2 = 0; // Number of parameters for part solution (DUAL mode) + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + PS_ASSERT_PTR_NON_NULL(stamp, false); + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + continue; + } + + PS_ASSERT_VECTOR_NON_NULL(stamp->vector1, false); + if (numParams == -1) { + numParams = stamp->vector1->n; + } + PS_ASSERT_VECTOR_SIZE(stamp->vector1, (long)numParams, false); + PS_ASSERT_VECTOR_TYPE(stamp->vector1, PS_TYPE_F64, false); + PS_ASSERT_IMAGE_NON_NULL(stamp->matrix1, false); + PS_ASSERT_IMAGE_SIZE(stamp->matrix1, numParams, numParams, false); + PS_ASSERT_IMAGE_TYPE(stamp->matrix1, PS_TYPE_F64, false); + + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + PS_ASSERT_IMAGE_NON_NULL(stamp->matrix2, false); + PS_ASSERT_IMAGE_NON_NULL(stamp->matrixX, false); + if (numParams2 == 0) { + numParams2 = stamp->matrix2->numCols; + } + PS_ASSERT_IMAGE_SIZE(stamp->matrix2, numParams2, numParams2, false); + PS_ASSERT_IMAGE_SIZE(stamp->matrixX, numParams, numParams2, false); + PS_ASSERT_IMAGE_TYPE(stamp->matrix2, PS_TYPE_F64, false); + PS_ASSERT_IMAGE_TYPE(stamp->matrixX, PS_TYPE_F64, false); + PS_ASSERT_VECTOR_NON_NULL(stamp->vector2, false); + PS_ASSERT_VECTOR_SIZE(stamp->vector2, (long)numParams2, false); + PS_ASSERT_VECTOR_TYPE(stamp->vector2, PS_TYPE_F64, false); + } + } + if (numParams == -1) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "No suitable stamps found."); + return NULL; + } + + if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { + // Accumulate the least-squares matricies and vectors + psImage *sumMatrix = psImageAlloc(numParams, numParams, PS_TYPE_F64); // Combined matrix + psVector *sumVector = psVectorAlloc(numParams, PS_TYPE_F64); // Combined vector + psVectorInit(sumVector, 0.0); + psImageInit(sumMatrix, 0.0); + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + (void)psBinaryOp(sumMatrix, sumMatrix, "+", stamp->matrix1); + (void)psBinaryOp(sumVector, sumVector, "+", stamp->vector1); + } + } + + psVector *permutation = NULL; // Permutation vector, required for LU decomposition + psImage *luMatrix = psMatrixLUD(NULL, &permutation, sumMatrix); + psFree(sumMatrix); + if (!luMatrix) { + psError(PS_ERR_UNKNOWN, true, "LU Decomposition of least-squares matrix failed.\n"); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + return NULL; + } + kernels->solution1 = psMatrixLUSolve(kernels->solution1, luMatrix, sumVector, permutation); + psFree(sumVector); + psFree(luMatrix); + psFree(permutation); + if (!kernels->solution1) { + psError(PS_ERR_UNKNOWN, true, "Failed to solve the least-squares system.\n"); + return NULL; + } + } else { + // Dual convolution solution + + // Accumulation of stamp matrices/vectors + psImage *sumMatrix1 = psImageAlloc(numParams, numParams, PS_TYPE_F64); + psImage *sumMatrix2 = psImageAlloc(numParams2, numParams2, PS_TYPE_F64); + psImage *sumMatrixX = psImageAlloc(numParams, numParams2, PS_TYPE_F64); + psVector *sumVector1 = psVectorAlloc(numParams, PS_TYPE_F64); + psVector *sumVector2 = psVectorAlloc(numParams, PS_TYPE_F64); + psImageInit(sumMatrix1, 0.0); + psImageInit(sumMatrix2, 0.0); + psImageInit(sumMatrixX, 0.0); + psVectorInit(sumVector1, 0.0); + psVectorInit(sumVector2, 0.0); + + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status == PM_SUBTRACTION_STAMP_USED) { + (void)psBinaryOp(sumMatrix1, sumMatrix1, "+", stamp->matrix1); + (void)psBinaryOp(sumMatrix2, sumMatrix2, "+", stamp->matrix2); + (void)psBinaryOp(sumMatrixX, sumMatrixX, "+", stamp->matrixX); + (void)psBinaryOp(sumVector1, sumVector1, "+", stamp->vector1); + (void)psBinaryOp(sumVector2, sumVector2, "+", stamp->vector2); + } + } + +#if 0 + // Apply weighting to maximise the difference between solution coefficients for the same functions + double fudge = PS_SQR(2 * stamps->footprint + 1); // Fudge factor + for (int i = 0; i < kernels->num; i++) { + sumMatrix1->data.F64[i][i] -= fudge; + sumMatrix2->data.F64[i][i] += fudge; + } +#endif + + // Pure matrix operations + + // A * a = Ct * b + d + // C * a = B * b + e + // + // a = (Ct * Bi * C - A)i (Ct * Bi * e - d) + // b = Bi * (C * a - e) + psVector *a = psVectorRecycle(kernels->solution1, numParams, PS_TYPE_F64); + psVector *b = psVectorRecycle(kernels->solution2, numParams2, PS_TYPE_F64); +#ifdef TESTING + psVectorInit(a, NAN); + psVectorInit(b, NAN); +#endif + psImage *A = sumMatrix1; + psImage *B = sumMatrix2; + psImage *C = sumMatrixX; + psVector *d = sumVector1; + psVector *e = sumVector2; + + assert(a->n == numParams); + assert(b->n == numParams2); + assert(A->numRows == numParams && A->numCols == numParams); + assert(B->numRows == numParams2 && B->numCols == numParams2); + assert(C->numRows == numParams2 && C->numCols == numParams); + assert(d->n == numParams); + assert(e->n == numParams2); + + psImage *Bi = psMatrixInvert(NULL, B, NULL); + assert(Bi->numRows == numParams2 && Bi->numCols == numParams2); + psImage *Ct = psMatrixTranspose(NULL, C); + assert(Ct->numRows == numParams && Ct->numCols == numParams2); + + psImage *BiC = psMatrixMultiply(NULL, Bi, C); + assert(BiC->numRows == numParams2 && BiC->numCols == numParams); + psImage *CtBi = psMatrixMultiply(NULL, Ct, Bi); + assert(CtBi->numRows == numParams && CtBi->numCols == numParams2); + + psImage *CtBiC = psMatrixMultiply(NULL, Ct, BiC); + assert(CtBiC->numRows == numParams && CtBiC->numCols == numParams); + + psImage *F = (psImage*)psBinaryOp(NULL, CtBiC, "-", A); + assert(F->numRows == numParams && F->numCols == numParams); + float det = NAN; + psImage *Fi = psMatrixInvert(NULL, F, &det); + assert(Fi->numRows == numParams && Fi->numCols == numParams); + psTrace("psModules.imcombine", 4, "Determinant of F: %f\n", det); + + psVector *g = psVectorAlloc(numParams, PS_TYPE_F64); +#ifdef TESTING + psVectorInit(g, NAN); +#endif + assert(CtBi->numRows == numParams && CtBi->numCols == numParams2); + assert(e->n == numParams2); + assert(d->n == numParams); + for (int i = 0; i < numParams; i++) { + double value = 0.0; + for (int j = 0; j < numParams2; j++) { + value += CtBi->data.F64[i][j] * e->data.F64[j]; + } + g->data.F64[i] = value - d->data.F64[i]; + } + + assert(Fi->numRows == numParams && Fi->numCols == numParams); + assert(g->n == numParams); + for (int i = 0; i < numParams; i++) { + double value = 0.0; + for (int j = 0; j < numParams; j++) { + value += Fi->data.F64[i][j] * g->data.F64[j]; + } + a->data.F64[i] = value; + } + + psVector *h = psVectorAlloc(numParams2, PS_TYPE_F64); +#ifdef TESTING + psVectorInit(h, NAN); +#endif + assert(C->numRows == numParams2 && C->numCols == numParams); + assert(a->n == numParams); + assert(e->n == numParams2); + for (int i = 0; i < numParams2; i++) { + double value = 0.0; + for (int j = 0; j < numParams; j++) { + value += C->data.F64[i][j] * a->data.F64[j]; + } + h->data.F64[i] = value - e->data.F64[i]; + } + + assert(Bi->numRows == numParams2 && Bi->numCols == numParams2); + assert(h->n == numParams2); + for (int i = 0; i < numParams2; i++) { + double value = 0.0; + for (int j = 0; j < numParams2; j++) { + value += Bi->data.F64[i][j] * h->data.F64[j]; + } + b->data.F64[i] = value; + } + + +#if 0 + for (int i = 0; i < numParams; i++) { + double aVal1 = 0.0, bVal1 = 0.0; + for (int j = 0; j < numParams2; j++) { + aVal1 += A->data.F64[i][j] * a->data.F64[j]; + bVal1 += Ct->data.F64[i][j] * b->data.F64[j]; + } + bVal1 += d->data.F64[i]; + for (int j = numParams2; j < numParams; j++) { + aVal1 += A->data.F64[i][j] * a->data.F64[j]; + } + printf("%d: %lf\n", i, aVal1 - bVal1); + } + + for (int i = 0; i < numParams2; i++) { + double aVal2 = 0.0, bVal2 = 0.0; + for (int j = 0; j < numParams2; j++) { + aVal2 += C->data.F64[i][j] * a->data.F64[j]; + bVal2 += B->data.F64[i][j] * b->data.F64[j]; + } + bVal2 += e->data.F64[i]; + for (int j = numParams2; j < numParams; j++) { + aVal2 += C->data.F64[i][j] * a->data.F64[j]; + } + printf("%d: %lf\n", i, aVal2 - bVal2); + } +#endif + + { + psFits *fits = psFitsOpen("sumMatrix1.fits", "w"); + psFitsWriteImage(fits, NULL, sumMatrix1, 0, NULL); + psFitsClose(fits); + } + { + psFits *fits = psFitsOpen("sumMatrix2.fits", "w"); + psFitsWriteImage(fits, NULL, sumMatrix2, 0, NULL); + psFitsClose(fits); + } + { + psFits *fits = psFitsOpen("sumMatrixX.fits", "w"); + psFitsWriteImage(fits, NULL, sumMatrixX, 0, NULL); + psFitsClose(fits); + } + { + psFits *fits = psFitsOpen("sumFinverse.fits", "w"); + psFitsWriteImage(fits, NULL, Fi, 0, NULL); + psFitsClose(fits); + } + + + kernels->solution1 = a; + kernels->solution2 = b; + + // XXXXX Free temporary matrices and vectors + + } + + if (psTraceGetLevel("psModules.imcombine") >= 7) { + for (int i = 0; i < kernels->solution1->n; i++) { + psTrace("psModules.imcombine", 7, "Solution 1 %d: %f\n", i, kernels->solution1->data.F64[i]); + } + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { + for (int i = 0; i < kernels->solution2->n; i++) { + psTrace("psModules.imcombine", 7, "Solution 2 %d: %f\n", i, kernels->solution2->data.F64[i]); + } + } + } + + return true; +} + +psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, + const pmSubtractionKernels *kernels) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, NULL); + PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(kernels, NULL); + + psVector *deviations = psVectorAlloc(stamps->num, PS_TYPE_F32); // Mean deviation for stamps + int footprint = stamps->footprint; // Half-size of stamps + long numPixels = PS_SQR(2 * footprint + 1); // Number of pixels in footprint + double devNorm = 1.0 / (double)numPixels; // Normalisation for deviations + int numKernels = kernels->num; // Number of kernels + + psImage *polyValues = NULL; // Polynomial values + psKernel *residual = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Residual image + + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // The stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_USED) { + deviations->data.F32[i] = NAN; + continue; + } + + // Calculate coefficients of the kernel basis functions + polyValues = p_pmSubtractionPolynomial(polyValues, kernels->spatialOrder, stamp->xNorm, stamp->yNorm); + double norm = p_pmSubtractionSolutionNorm(kernels); // Normalisation + double background = p_pmSubtractionSolutionBackground(kernels, polyValues);// Difference in background + + // Calculate residuals + psKernel *weight = stamp->weight; // Weight postage stamp + psImageInit(residual->image, 0.0); + if (kernels->mode != PM_SUBTRACTION_MODE_DUAL) { + psKernel *target; // Target postage stamp + psKernel *source; // Source postage stamp + psArray *convolutions; // Convolution postage stamps for each kernel basis function + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_TARGET: + case PM_SUBTRACTION_MODE_1: + target = stamp->image2; + source = stamp->image1; + convolutions = stamp->convolutions1; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + source = stamp->image2; + convolutions = stamp->convolutions2; + break; + default: + psAbort("Unsupported subtraction mode: %x", kernels->mode); + } + + for (int j = 0; j < numKernels; j++) { + psKernel *convolution = convolutions->data[j]; // Convolution + double coefficient = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, + false); // Coefficient + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] -= convolution->kernel[y][x] * coefficient; + } + } + } + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += target->kernel[y][x] - background - source->kernel[y][x] * norm; + } + } + } else { + // Dual convolution + psArray *convolutions1 = stamp->convolutions1; // Convolutions of the first image + psArray *convolutions2 = stamp->convolutions2; // Convolutions of the second image + psKernel *image1 = stamp->image1; // The first image + psKernel *image2 = stamp->image2; // The second image + + for (int j = 0; j < numKernels; j++) { + psKernel *conv1 = convolutions1->data[j]; // Convolution of first image + psKernel *conv2 = convolutions2->data[j]; // Convolution of second image + double coeff1 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, false); // Coefficient 1 + double coeff2 = p_pmSubtractionSolutionCoeff(kernels, polyValues, j, true); // Coefficient 2 + + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += conv2->kernel[y][x] * coeff2 - conv1->kernel[y][x] * coeff1; + } + } + } + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + residual->kernel[y][x] += image2->kernel[y][x] - background - image1->kernel[y][x] * norm; + } + } + } + + double deviation = 0.0; // Sum of differences + for (int y = - footprint; y <= footprint; y++) { + for (int x = - footprint; x <= footprint; x++) { + double dev = PS_SQR(residual->kernel[y][x]) / weight->kernel[y][x]; + deviation += dev; +#ifdef TESTING + residual->kernel[y][x] = dev; +#endif + } + } + deviations->data.F32[i] = sqrtf(devNorm * deviation); + psTrace("psModules.imcombine", 5, "Deviation for stamp %d (%d,%d): %f\n", + i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5), deviations->data.F32[i]); + if (!isfinite(deviations->data.F32[i])) { + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + psTrace("psModules.imcombine", 5, + "Rejecting stamp %d (%d,%d) because of non-finite deviation\n", + i, (int)(stamp->x + 0.5), (int)(stamp->y + 0.5)); + continue; + } + +#ifdef TESTING + { + psString filename = NULL; + psStringAppend(&filename, "resid_%03d.fits", i); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, residual->image, 0, NULL); + psFitsClose(fits); + } + { + psString filename = NULL; + psStringAppend(&filename, "stamp_image_%03d.fits", i); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, stamp->image1->image, 0, NULL); + psFitsClose(fits); + } + { + psString filename = NULL; + psStringAppend(&filename, "stamp_weight_%03d.fits", i); + psFits *fits = psFitsOpen(filename, "w"); + psFree(filename); + psFitsWriteImage(fits, NULL, stamp->weight->image, 0, NULL); + psFitsClose(fits); + } +#endif + + } + psFree(residual); + psFree(polyValues); + + return deviations; +} Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionEquation.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionEquation.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionEquation.h (revision 16477) @@ -0,0 +1,46 @@ +#ifndef PM_SUBTRACTION_EQUATION_H +#define PM_SUBTRACTION_EQUATION_H + +#include "pmSubtractionStamps.h" +#include "pmSubtractionKernels.h" + +/// Calculate the least-squares equation to match the image quality +bool pmSubtractionCalculateEquation(pmSubtractionStampList *stamps, ///< Stamps + const pmSubtractionKernels *kernels ///< Kernel parameters + ); + +/// Solve the least-squares equation to match the image quality +bool pmSubtractionSolveEquation(pmSubtractionKernels *kernels, ///< Kernel parameters + const pmSubtractionStampList *stamps ///< Stamps + ); + +/// Calculate deviations +psVector *pmSubtractionCalculateDeviations(pmSubtractionStampList *stamps, ///< Stamps + const pmSubtractionKernels *kernels ///< Kernel parameters + ); + +/// Calculate the value of a polynomial, specified by coefficients and polynomial values +inline double p_pmSubtractionCalculatePolynomial(const psVector *coeff, ///< Coefficients + const psImage *polyValues, ///< Polynomial values + int order, ///< Order of polynomials + int index, ///< Index at which to begin + int step ///< Step between subsequent indices + ); + +/// Return the specified coefficient in the solution +inline double p_pmSubtractionSolutionCoeff(const pmSubtractionKernels *kernels, ///< Kernel parameters + const psImage *polyValues, ///< Polynomial values + int index, ///< Coefficient index to calculate + bool wantDual ///< Calculate the coefficient for the dual solution? + ); + +/// Return the normalisation in the solution +inline double p_pmSubtractionSolutionNorm(const pmSubtractionKernels *kernels ///< Kernel parameters + ); + +/// Return the background (difference) in the solution +inline double p_pmSubtractionSolutionBackground(const pmSubtractionKernels *kernels, ///< Kernel parameters + const psImage *polyValues ///< Polynomial values + ); + +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionKernels.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionKernels.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionKernels.c (revision 16477) @@ -0,0 +1,637 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmSubtraction.h" +#include "pmSubtractionKernels.h" + +#define RINGS_BUFFER 10 // Buffer size for RINGS data + + +// Free function for pmSubtractionKernels +static void subtractionKernelsFree(pmSubtractionKernels *kernels) +{ + psFree(kernels->description); + psFree(kernels->u); + psFree(kernels->v); + psFree(kernels->widths); + psFree(kernels->uStop); + psFree(kernels->vStop); + psFree(kernels->preCalc); + psFree(kernels->solution1); + psFree(kernels->solution2); +} + +// Raise an integer to an integer power +static inline long power(int value, // Value + int exp // Exponent + ) +{ + if (exp == 0) { + return 1.0; + } + long result = value; // Result to return + for (int i = 2; i <= exp; i++) { + result *= value; + } + return result; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Semi-public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, int start, int size) +{ + PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(kernels, false); + PS_ASSERT_INT_NONNEGATIVE(start, false); + PS_ASSERT_INT_NONNEGATIVE(size, false); + + int numNew = PS_SQR(2 * size + 1) - 1; // Number of new kernel parameters to add + + // Ensure the sizes match + kernels->widths = psVectorRealloc(kernels->widths, start + numNew); + kernels->u = psVectorRealloc(kernels->u, start + numNew); + kernels->v = psVectorRealloc(kernels->v, start + numNew); + kernels->preCalc = psArrayRealloc(kernels->preCalc, start + numNew); + kernels->inner = start; + + // Generate a set of kernels for each (u,v) + for (int v = - size, index = start; v <= size; v++) { + for (int u = - size; u <= size; u++, index++) { + if (v == 0 && u == 0) { + // Skip normalisation component: added explicitly + index--; + continue; + } + kernels->widths->data.F32[index] = NAN; + kernels->u->data.S32[index] = u; + kernels->v->data.S32[index] = v; + kernels->preCalc->data[index] = NULL; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d\n", index, u, v); + } + } + + return true; +} + +pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, int spatialOrder, + const psVector *fwhms, const psVector *orders, + pmSubtractionMode mode) +{ + PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); + PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(orders, NULL); + PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(fwhms, orders, NULL); + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + + int numGaussians = fwhms->n; // Number of Gaussians + + int num = 0; // Number of basis functions + psString params = NULL; // List of parameters + for (int i = 0; i < numGaussians; i++) { + int gaussOrder = orders->data.S32[i]; // Polynomial order to apply to Gaussian + psStringAppend(¶ms, "(%.1f,%d)", fwhms->data.F32[i], orders->data.S32[i]); + num += (gaussOrder + 1) * (gaussOrder + 2) / 2; + } + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_ISIS, + size, spatialOrder, mode); // The kernels + psStringAppend(&kernels->description, "ISIS(%d,%s,%d)", size, params, spatialOrder); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "ISIS kernel: %s,%d --> %d elements", + params, spatialOrder, num); + psFree(params); + + // Set the kernel parameters + for (int i = 0, index = 0; i < numGaussians; i++) { + float sigma = fwhms->data.F32[i] / (2.0 * sqrtf(2.0 * logf(2.0))); // Gaussian sigma + float norm = 1.0 / (M_2_PI * sqrtf(sigma)); // Normalisation for Gaussian + // Iterate over (u,v) order + for (int uOrder = 0; uOrder <= orders->data.S32[i]; uOrder++) { + for (int vOrder = 0; vOrder <= orders->data.S32[i] - uOrder; vOrder++, index++) { + // Set the pre-calculated kernel + psKernel *preCalc = psKernelAlloc(-size, size, -size, size); + double sum = 0.0; // Normalisation + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + sum += preCalc->kernel[v][u] = norm * power(u, uOrder) * power(v, vOrder) * + expf(-0.5 * (PS_SQR(u) + PS_SQR(v)) / PS_SQR(sigma)); + } + } + + // Normalise sum of kernel component to unity for even functions + if (uOrder % 2 == 0 && vOrder % 2 == 0) { + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + preCalc->kernel[v][u] = preCalc->kernel[v][u] / sum; + } + } + preCalc->kernel[0][0] -= 1.0; + } + + kernels->widths->data.F32[index] = fwhms->data.F32[i]; + kernels->u->data.S32[index] = uOrder; + kernels->v->data.S32[index] = vOrder; + if (kernels->preCalc->data[index]) { + psFree(kernels->preCalc->data[index]); + } + kernels->preCalc->data[index] = preCalc; + + psTrace("psModules.imcombine", 7, "Kernel %d: %f %d %d\n", index, + fwhms->data.F32[i], uOrder, vOrder); + } + } + } + + return kernels; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, pmSubtractionKernelsType type, + int size, int spatialOrder, pmSubtractionMode mode) +{ + pmSubtractionKernels *kernels = psAlloc(sizeof(pmSubtractionKernels)); // Kernels, to return + psMemSetDeallocator(kernels, (psFreeFunc)subtractionKernelsFree); + + kernels->type = type; + kernels->description = NULL; + kernels->num = numBasisFunctions; + kernels->u = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); + kernels->v = psVectorAlloc(numBasisFunctions, PS_TYPE_S32); + kernels->widths = psVectorAlloc(numBasisFunctions, PS_TYPE_F32); + kernels->preCalc = psArrayAlloc(numBasisFunctions); + kernels->uStop = NULL; + kernels->vStop = NULL; + kernels->size = size; + kernels->inner = 0; + kernels->spatialOrder = spatialOrder; + kernels->bgOrder = 0; + kernels->mode = mode; + kernels->numCols = 0; + kernels->numRows = 0; + kernels->solution1 = NULL; + kernels->solution2 = NULL; + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, int spatialOrder, pmSubtractionMode mode) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + + int num = PS_SQR(2 * size + 1) - 1; // Number of basis functions + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_POIS, + size, spatialOrder, mode); // The kernels + psStringAppend(&kernels->description, "POIS(%d,%d)", size, spatialOrder); + psLogMsg("psModules.imcombine", PS_LOG_INFO, "POIS kernel: %d,%d --> %d elements", + size, spatialOrder, num); + + if (!p_pmSubtractionKernelsAddGrid(kernels, 0, size)) { + psAbort("Should never get here."); + } + + return kernels; +} + + +pmSubtractionKernels *pmSubtractionKernelsISIS(int size, int spatialOrder, + const psVector *fwhms, const psVector *orders, + pmSubtractionMode mode) +{ + pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, + fwhms, orders, mode); // Kernels + if (!kernels) { + return NULL; + } + + if (psTraceGetLevel("psModules.imcombine.kernel") >= 10) { + for (int i = 0; i < kernels->num; i++) { + psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest + psString kernelName = NULL; + psStringAppend(&kernelName, "kernel%03d.fits", i); + psFits *kernelFile = psFitsOpen(kernelName, "w"); + psFree(kernelName); + psFitsWriteImage(kernelFile, NULL, kernel->image, 0, NULL); + psFitsClose(kernelFile); + double sum = 0.0; + for (int y = 0; y < kernel->image->numRows; y++) { + for (int x = 0; x < kernel->image->numCols; x++) { + sum += kernel->image->data.F32[y][x]; + } + } + psTrace("psModules.imcombine.kernel", 10, "Kernel %d sum: %le\n", i, sum); + } + } + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, int spatialOrder, int inner, int binning, + pmSubtractionMode mode) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); + PS_ASSERT_INT_POSITIVE(binning, NULL); + + // The outer region should be divisible by the "binning"; otherwise allocate remainder to the inner region + int numOuter = (size - inner) / binning; // Number of summed pixels in the outer region + int numInner = inner + (size - inner) % binning; // Number of pixels in the inner region + assert(numOuter * binning + numInner == size); + int numTotal = numOuter + numInner; // Total number of summed pixels + + psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); + + int num = PS_SQR(2 * numTotal + 1) - 1; // Number of basis functions + + psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_SPAM, + size, spatialOrder, mode); // The kernels + psStringAppend(&kernels->description, "SPAM(%d,%d,%d,%d)", size, inner, binning, spatialOrder); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "SPAM kernel: %d,%d,%d,%d --> %d elements", + size, inner, binning, spatialOrder, num); + + kernels->uStop = psVectorAlloc(num, PS_TYPE_S32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_S32); + + psVector *locations = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Locations for each kernel element + psVector *widths = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); // Widths for each kernel element + locations->data.S32[numTotal] = 0; + widths->data.S32[numTotal] = 0; + for (int i = 1; i <= numInner; i++) { + locations->data.S32[numTotal + i] = i; + widths->data.S32[numTotal + i] = 0; + locations->data.S32[numTotal - i] = - i; + widths->data.S32[numTotal - i] = 0; + } + for (int i = numInner + 1; i <= numTotal; i++) { + locations->data.S32[numTotal + i] = locations->data.S32[numTotal + i - 1] + + widths->data.S32[numTotal + i - 1] + 1; + widths->data.S32[numTotal + i] = binning - 1; + locations->data.S32[numTotal - i] = locations->data.S32[numTotal - i + 1] - binning; + widths->data.S32[numTotal - i] = binning - 1; + } + + if (psTraceGetLevel("psModules.imcombine") >= 10) { + for (int i = 0; i < 2 * numTotal + 1; i++) { + psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, locations->data.S32[i], + locations->data.S32[i] + widths->data.S32[i]); + } + } + + // Set the kernel parameters + for (int i = - numTotal, index = 0; i <= numTotal; i++) { + int u = locations->data.S32[numTotal + i]; // Location of pixel + int uStop = u + widths->data.S32[numTotal + i]; // Width of pixel + + for (int j = - numTotal; j <= numTotal; j++, index++) { + if (i == 0 && j == 0) { + // Skip normalisation component: added explicitly + index--; + continue; + } + int v = locations->data.S32[numTotal + j]; // Location of pixel + int vStop = v + widths->data.S32[numTotal + j]; // Width of pixel + + kernels->u->data.S32[index] = u; + kernels->v->data.S32[index] = v; + kernels->uStop->data.S32[index] = uStop; + kernels->vStop->data.S32[index] = vStop; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, + u, uStop, v, vStop); + } + } + + psFree(locations); + psFree(widths); + + return kernels; +} + + +pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, int spatialOrder, int inner, pmSubtractionMode mode) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LARGER_THAN(size, inner, NULL); + + int fibNum = 0; // Number of Fibonacci values + int fibLast = 1, fibTotal = 2; // Fibonacci sequence + while (fibTotal < size - inner) { + int temp = fibTotal; + fibTotal += fibLast; + fibLast = temp; + fibNum++; + } + + int numInner = inner; // Number of pixels in the inner region + int numOuter = fibNum; // Number of summed pixels in the outer region + int numTotal = numOuter + numInner; // Total number of summed pixels + + psTrace("psModules.imcombine", 3, "Inner: %d Outer: %d\n", numInner, numOuter); + + int num = PS_SQR(2 * numTotal + 1) - 1; // Number of basis functions + + psTrace("psModules.imcombine", 3, "Number of basis functions: %d\n", num); + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_FRIES, + size, spatialOrder, mode); // The kernels + psStringAppend(&kernels->description, "FRIES(%d,%d,%d)", size, inner, spatialOrder); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "FRIES kernel: %d,%d,%d --> %d elements", + size, inner, spatialOrder, num); + + kernels->uStop = psVectorAlloc(num, PS_TYPE_S32); + kernels->vStop = psVectorAlloc(num, PS_TYPE_S32); + + psVector *start = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); + psVector *stop = psVectorAlloc(2 * numTotal + 1, PS_TYPE_S32); + start->data.S32[numTotal] = 0; + stop->data.S32[numTotal] = 0; + for (int i = 1; i <= numInner; i++) { + start->data.S32[numTotal + i] = i; + stop->data.S32[numTotal + i] = i; + start->data.S32[numTotal - i] = -i; + stop->data.S32[numTotal - i] = -i; + } + for (int i = numInner + 1, fibLast = 1, fib = 2, temp; i <= numTotal; + i++, fib = (temp = fib) + fibLast, fibLast = temp) { + start->data.S32[numTotal + i] = stop->data.S32[numTotal + i - 1] + 1; + stop->data.S32[numTotal + i] = PS_MIN(start->data.S32[numTotal + i] + fib - 1, size); + start->data.S32[numTotal - i] = - stop->data.S32[numTotal + i]; + stop->data.S32[numTotal - i] = - start->data.S32[numTotal + i]; + } + + if (psTraceGetLevel("psModules.imcombine") >= 10) { + for (int i = 0; i < 2 * numTotal + 1; i++) { + psTrace("psModules.imcombine", 10, "%d: %d -> %d\n", i, start->data.S32[i], stop->data.S32[i]); + } + } + + // Set the kernel parameters + for (int i = - numTotal, index = 0; i <= numTotal; i++) { + int u = start->data.S32[numTotal + i]; // Location of pixel + int uStop = stop->data.S32[numTotal + i]; // Width of pixel + for (int j = - numTotal; j <= numTotal; j++, index++) { + if (i == 0 && j == 0) { + // Skip normalisation component: added explicitly + index--; + continue; + } + int v = start->data.S32[numTotal + j]; // Location of pixel + int vStop = stop->data.S32[numTotal + j]; // Width of pixel + + kernels->u->data.S32[index] = u; + kernels->v->data.S32[index] = v; + kernels->uStop->data.S32[index] = uStop; + kernels->vStop->data.S32[index] = vStop; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d %d\n", index, + u, uStop, v, vStop); + } + } + + psFree(start); + psFree(stop); + + return kernels; +} + +// Grid United with Normal Kernel +pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, int spatialOrder, const psVector *fwhms, + const psVector *orders, int inner, pmSubtractionMode mode) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); + PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(orders, NULL); + PS_ASSERT_VECTOR_TYPE(orders, PS_TYPE_S32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(fwhms, orders, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LESS_THAN(inner, size, NULL); + + pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, + fwhms, orders, mode); // Kernels + psStringPrepend(&kernels->description, "GUNK="); + psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); + + int numISIS = kernels->num; // Number of ISIS kernels + + if (!p_pmSubtractionKernelsAddGrid(kernels, numISIS, inner)) { + psAbort("Should never get here."); + } + + return kernels; +} + +// RINGS --- just what it says +pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, int spatialOrder, int inner, int ringsOrder, + pmSubtractionMode mode) +{ + PS_ASSERT_INT_POSITIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_LESS_THAN_OR_EQUAL(inner, size, NULL); + PS_ASSERT_INT_NONNEGATIVE(ringsOrder, NULL); + + int fibNum = 0; // Number of Fibonacci values + { + int fibIndex = 1, fibIndexMinus1 = 0; // Fibonnacci parameters + int radius = inner; + while (radius + fibIndex < size) { + radius++; + int fibNew = fibIndex + fibIndexMinus1; + fibIndexMinus1 = fibIndex; + fibIndex = fibNew; + radius += fibIndex; + fibNum++; + } + } + + int numInner = inner - 1; // Number of pixels in the inner region + int numOuter = fibNum; // Number of summed pixels in the outer region + + int numRings = numOuter + numInner; // Number of rings (not including the central pixel) + int numPoly = PM_SUBTRACTION_POLYTERMS(ringsOrder); // Number of polynomial variants of each ring + + int num = numRings * numPoly; // Total number of basis functions + + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(num, PM_SUBTRACTION_KERNEL_RINGS, + size, spatialOrder, mode); // The kernels + psStringAppend(&kernels->description, "RINGS(%d,%d,%d,%d)", size, inner, ringsOrder, spatialOrder); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "RINGS kernel: %d,%d,%d,%d --> %d elements", + size, inner, ringsOrder, spatialOrder, num); + + // Set the Gaussian kernel parameters + int fibIndex = 1, fibIndexMinus1 = 0; // Fibonnacci parameters + int radiusLast = 1; // Last radius + for (int i = 1, index = 0; i < numRings + 1; i++) { + float lower2; // Lower limit of radius^2 + float upper2; // Upper limit of radius^2 + if (i <= inner) { + // A ring every pixel width + float radius = i; + lower2 = PS_SQR(radius - 0.5); + upper2 = PS_SQR(radius + 0.5); + radiusLast = i; + } else { + // Rings Fibonacci distributed (2, 3, 5...) + int fibNew = fibIndex + fibIndexMinus1; + fibIndexMinus1 = fibIndex; + fibIndex = fibNew; + + float radiusLower = radiusLast + 1; + radiusLast = radiusLower + fibIndex; + float radiusUpper = radiusLast; + + lower2 = PS_SQR(radiusLower - 0.5); + upper2 = PS_SQR(radiusUpper + 0.5); + } + + psTrace("psModules.imcombine", 8, "Radius limits: %f --> %f\n", sqrtf(lower2), sqrtf(upper2)); + + // Iterate over (u,v) order + for (int uOrder = 0; uOrder <= (i == 0 ? 0 : ringsOrder); uOrder++) { + for (int vOrder = 0; vOrder <= (i == 0 ? 0 : ringsOrder - uOrder); vOrder++, index++) { + + psArray *data = psArrayAlloc(3); // Container for data + psVector *uCoords = data->data[0] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // u coords + psVector *vCoords = data->data[1] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_S32); // v coords + psVector *poly = data->data[2] = psVectorAllocEmpty(RINGS_BUFFER, PS_TYPE_F32); // Polynomial + + if (i == 0) { + // Central pixel is easy + uCoords->data.S32[0] = vCoords->data.S32[0] = 0; + poly->data.F32[0] = 1.0; + uCoords->n = vCoords->n = poly->n = 1; + radiusLast = 0; + } else { + int j = 0; // Index for data + double norm = 0.0; // Normalisation + for (int v = -size; v <= size; v++) { + int v2 = PS_SQR(v); // Square of v + float vPoly = powf(v/(float)size, vOrder); // Value of v^vOrder + + for (int u = -size; u <= size; u++) { + int u2 = PS_SQR(u); // Square of u + int distance2 = u2 + v2; // Distance from the centre + if (distance2 > lower2 && distance2 < upper2) { + float uPoly = powf(u/(float)size, uOrder); // Value of u^uOrder + + float polyVal = uPoly * vPoly; // Value of polynomial + if (polyVal != 0) { // No point adding it otherwise + uCoords->data.S32[j] = u; + vCoords->data.S32[j] = v; + poly->data.F32[j] = polyVal; + norm += polyVal; + + psVectorExtend(uCoords, RINGS_BUFFER, 1); + psVectorExtend(vCoords, RINGS_BUFFER, 1); + psVectorExtend(poly, RINGS_BUFFER, 1); + psTrace("psModules.imcombine", 9, "u = %d, v = %d, poly = %f\n", + u, v, poly->data.F32[j]); + j++; + } + } + } + } + // Normalise kernel component to unit sum + if (uOrder % 2 == 0 && vOrder % 2 == 0) { + psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); + // Add subtraction of 0,0 component to preserve photometric scaling + uCoords->data.S32[j] = 0; + vCoords->data.S32[j] = 0; + poly->data.F32[j] = -1.0; + psVectorExtend(uCoords, RINGS_BUFFER, 1); + psVectorExtend(vCoords, RINGS_BUFFER, 1); + psVectorExtend(poly, RINGS_BUFFER, 1); + } else { + norm = powf(size, uOrder) * powf(size, vOrder); + psBinaryOp(poly, poly, "*", psScalarAlloc(1.0 / norm, PS_TYPE_F32)); + } + } + + psTrace("psModules.imcombine", 8, "%ld pixels in kernel\n", uCoords->n); + + kernels->preCalc->data[index] = data; + kernels->u->data.S32[index] = uOrder; + kernels->v->data.S32[index] = vOrder; + + psTrace("psModules.imcombine", 7, "Kernel %d: %d %d %d\n", index, + i, uOrder, vOrder); + } + } + } + + return kernels; +} + +pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, int size, int spatialOrder, + const psVector *fwhms, const psVector *orders, int inner, + int binning, int ringsOrder, pmSubtractionMode mode) +{ + switch (type) { + case PM_SUBTRACTION_KERNEL_POIS: + return pmSubtractionKernelsPOIS(size, spatialOrder, mode); + case PM_SUBTRACTION_KERNEL_ISIS: + return pmSubtractionKernelsISIS(size, spatialOrder, fwhms, orders, mode); + case PM_SUBTRACTION_KERNEL_SPAM: + return pmSubtractionKernelsSPAM(size, spatialOrder, inner, binning, mode); + case PM_SUBTRACTION_KERNEL_FRIES: + return pmSubtractionKernelsFRIES(size, spatialOrder, inner, mode); + case PM_SUBTRACTION_KERNEL_GUNK: + return pmSubtractionKernelsGUNK(size, spatialOrder, fwhms, orders, inner, mode); + case PM_SUBTRACTION_KERNEL_RINGS: + return pmSubtractionKernelsRINGS(size, spatialOrder, inner, ringsOrder, mode); + default: + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unknown kernel type: %x", type); + return NULL; + } +} + + +pmSubtractionKernelsType pmSubtractionKernelsTypeFromString(const char *type) +{ + if (strcasecmp(type, "POIS") == 0) { + return PM_SUBTRACTION_KERNEL_POIS; + } + if (strcasecmp(type, "ISIS") == 0) { + return PM_SUBTRACTION_KERNEL_ISIS; + } + if (strcasecmp(type, "SPAM") == 0) { + return PM_SUBTRACTION_KERNEL_SPAM; + } + if (strcasecmp(type, "FRIES") == 0) { + return PM_SUBTRACTION_KERNEL_FRIES; + } + if (strcasecmp(type, "GUNK") == 0) { + return PM_SUBTRACTION_KERNEL_GUNK; + } + if (strcasecmp(type, "RINGS") == 0) { + return PM_SUBTRACTION_KERNEL_RINGS; + } + + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unrecognised kernel type: %s", type); + return PM_SUBTRACTION_KERNEL_NONE; +} Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionKernels.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionKernels.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionKernels.h (revision 16477) @@ -0,0 +1,185 @@ +#ifndef PM_SUBTRACTION_KERNELS_H +#define PM_SUBTRACTION_KERNELS_H + +#include + +/// Type of subtraction kernel +typedef enum { + PM_SUBTRACTION_KERNEL_NONE, ///< Nothing --- an error + PM_SUBTRACTION_KERNEL_POIS, ///< Pan-STARRS Optimal Image Subtraction --- delta functions + PM_SUBTRACTION_KERNEL_ISIS, ///< Traditional kernel --- gaussians modified by polynomials + PM_SUBTRACTION_KERNEL_SPAM, ///< Summed Pixels for Advanced Matching --- summed delta functions + PM_SUBTRACTION_KERNEL_FRIES, ///< Fibonacci Radius Increases Excellence of Subtraction + PM_SUBTRACTION_KERNEL_GUNK, ///< Grid United with Normal Kernel --- POIS and ISIS hybrid + PM_SUBTRACTION_KERNEL_RINGS, ///< Rings Instead of the Normal Gaussian Subtraction +} pmSubtractionKernelsType; + +/// Modes --- specifies which image to convolve +typedef enum { + PM_SUBTRACTION_MODE_ERR, // Error in the mode + PM_SUBTRACTION_MODE_TARGET, // Convolve image 1 to match target PSF + PM_SUBTRACTION_MODE_1, // Convolve image 1 + PM_SUBTRACTION_MODE_2, // Convolve image 2 + PM_SUBTRACTION_MODE_UNSURE, // Not sure yet which image to convolve so try to satisfy both + PM_SUBTRACTION_MODE_DUAL, // Dual convolution +} pmSubtractionMode; + +/// Kernels specification +typedef struct { + pmSubtractionKernelsType type; ///< Type of kernels --- allowing the use of multiple kernels + psString description; ///< Description of the kernel parameters + long num; ///< Number of kernel components (not including the spatial ones) + psVector *u, *v; ///< Offset (for POIS) or polynomial order (for ISIS) + psVector *widths; ///< Gaussian FWHMs (ISIS) + psVector *uStop, *vStop; ///< Width of kernel element (SPAM,FRIES only) + psArray *preCalc; ///< Array of images containing pre-calculated kernel (for ISIS) + int size; ///< The half-size of the kernel + int inner; ///< The size of an inner region + int spatialOrder; ///< The spatial order of the kernels + int bgOrder; ///< The order for the background fitting + pmSubtractionMode mode; ///< Mode for subtraction + int numCols, numRows; ///< Size of image (for normalisation), or zero to use image provided + psVector *solution1, *solution2; ///< Solution for the PSF matching +} pmSubtractionKernels; + +// Assertion to check pmSubtractionKernels +#define PM_ASSERT_SUBTRACTION_KERNELS_NON_NULL(KERNELS, RETURNVALUE) { \ + PS_ASSERT_PTR_NON_NULL(KERNELS, RETURNVALUE); \ + PS_ASSERT_STRING_NON_EMPTY((KERNELS)->description, RETURNVALUE); \ + PS_ASSERT_INT_POSITIVE((KERNELS)->num, RETURNVALUE); \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->u, RETURNVALUE); \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->v, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->u, PS_TYPE_S32, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->v, PS_TYPE_S32, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->u, (KERNELS)->num, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->v, (KERNELS)->num, RETURNVALUE); \ + if ((KERNELS)->type == PM_SUBTRACTION_KERNEL_ISIS) { \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->widths, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->widths, PS_TYPE_F32, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->widths, (KERNELS)->num, RETURNVALUE); \ + } \ + if ((KERNELS)->uStop || (KERNELS)->vStop) { \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->uStop, RETURNVALUE); \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->vStop, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->uStop, PS_TYPE_S32, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->vStop, PS_TYPE_S32, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->uStop, (KERNELS)->num, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->vStop, (KERNELS)->num, RETURNVALUE); \ + } \ + if ((KERNELS)->preCalc) { \ + PS_ASSERT_ARRAY_NON_NULL((KERNELS)->preCalc, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((KERNELS)->preCalc, (KERNELS)->num, RETURNVALUE); \ + } \ + PS_ASSERT_INT_NONNEGATIVE((KERNELS)->size, RETURNVALUE); \ + PS_ASSERT_INT_NONNEGATIVE((KERNELS)->inner, RETURNVALUE); \ + PS_ASSERT_INT_NONNEGATIVE((KERNELS)->spatialOrder, RETURNVALUE); \ + PS_ASSERT_INT_NONNEGATIVE((KERNELS)->bgOrder, RETURNVALUE); \ +} + +// Assertion to check that the solution is attached +#define PM_ASSERT_SUBTRACTION_KERNELS_SOLUTION(KERNELS, RETURNVALUE) { \ + PS_ASSERT_VECTOR_NON_NULL((KERNELS)->solution1, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->solution1, PS_TYPE_F64, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->solution1, \ + (KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder) + 1 + \ + PM_SUBTRACTION_POLYTERMS((KERNELS)->bgOrder), \ + RETURNVALUE); \ + if (kernels->mode == PM_SUBTRACTION_MODE_DUAL) { \ + PS_ASSERT_VECTOR_NON_NULL(kernels->solution2, RETURNVALUE); \ + PS_ASSERT_VECTOR_TYPE((KERNELS)->solution2, PS_TYPE_F64, RETURNVALUE); \ + PS_ASSERT_VECTOR_SIZE((KERNELS)->solution2, \ + (KERNELS)->num * PM_SUBTRACTION_POLYTERMS((KERNELS)->spatialOrder), \ + RETURNVALUE); \ + } \ +} + +/// Generate a delta-function grid for subtraction kernels (like the POIS kernel) +bool p_pmSubtractionKernelsAddGrid(pmSubtractionKernels *kernels, ///< The subtraction kernels to append to + int start, ///< Index at which to start appending + int size ///< Half-size of the grid + ); + +/// General allocator for pmSubtractionKernels +/// +/// Unlike the functions for the specific kernel type, this function does not set up the basis functions, but +/// merely allocates space for their storage. +pmSubtractionKernels *pmSubtractionKernelsAlloc(int numBasisFunctions, ///< Number of basis functions + pmSubtractionKernelsType type, ///< Kernel type + int size, ///< Half-size of kernel + int spatialOrder, ///< Order of spatial variations + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Generate POIS kernels +pmSubtractionKernels *pmSubtractionKernelsPOIS(int size, ///< Half-size of the kernel (in both dims) + int spatialOrder, ///< Order of spatial variations + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Generate ISIS kernels without the flux scaling built in +pmSubtractionKernels *p_pmSubtractionKernelsRawISIS(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + const psVector *fwhms, ///< Gaussian FWHMs + const psVector *orders, ///< Polynomial order of gaussians + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Generate ISIS kernels +pmSubtractionKernels *pmSubtractionKernelsISIS(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + const psVector *fwhms, ///< Gaussian FWHMs + const psVector *orders, ///< Polynomial order of gaussians + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Generate SPAM kernels +pmSubtractionKernels *pmSubtractionKernelsSPAM(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + int inner, ///< Inner radius to preserve unbinned + int binning, ///< Kernel binning factor + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Generate FRIES kernels +pmSubtractionKernels *pmSubtractionKernelsFRIES(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + int inner, ///< Inner radius to preserve unbinned + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Generate GUNK kernels +pmSubtractionKernels *pmSubtractionKernelsGUNK(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + const psVector *fwhms, ///< Gaussian FWHMs + const psVector *orders, ///< Polynomial order of gaussians + int inner, ///< Inner radius containing grid of delta functions + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Generate RINGS kernels +pmSubtractionKernels *pmSubtractionKernelsRINGS(int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + int inner, ///< Inner radius to preserve unbinned + int ringsOrder, ///< Polynomial order + pmSubtractionMode mode ///< Mode for subtraction + ); + + +/// Generate a kernel of a specified type +pmSubtractionKernels *pmSubtractionKernelsGenerate(pmSubtractionKernelsType type, ///< Kernel type + int size, ///< Half-size of the kernel + int spatialOrder, ///< Order of spatial variations + const psVector *fwhms, ///< Gaussian FWHMs + const psVector *orders, ///< Polynomial order of gaussians + int inner, ///< Inner radius to preserve unbinned + int binning, ///< Kernel binning factor + int ringsOrder, ///< Polynomial order for RINGS + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Return the appropriate type from a string +pmSubtractionKernelsType pmSubtractionKernelsTypeFromString(const char *string // String name for kernel type + ); + + +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMask.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMask.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMask.c (revision 16477) @@ -0,0 +1,240 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#include "pmSubtraction.h" + +#include "pmSubtractionMask.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private (file-static) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Mark a pixel as blank in the image, mask and weight +static inline void markBlank(psImage *image, // Image to mark as blank + psImage *mask, // Mask to mark as blank (or NULL) + psImage *weight, // Weight map to mark as blank (or NULL) + int x, int y, // Coordinates to mark blank + psMaskType blank // Blank mask value + ) +{ + image->data.F32[y][x] = NAN; + if (mask) { + mask->data.PS_TYPE_MASK_DATA[y][x] |= blank; + } + if (weight) { + weight->data.F32[y][x] = NAN; + } + return; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +psImage *pmSubtractionMask(const psImage *mask1, const psImage *mask2, psMaskType maskVal, + int size, int footprint, float badFrac, bool useFFT) +{ + PS_ASSERT_IMAGE_NON_NULL(mask1, NULL); + PS_ASSERT_IMAGE_TYPE(mask1, PS_TYPE_MASK, NULL); + if (mask2) { + PS_ASSERT_IMAGE_NON_NULL(mask2, NULL); + PS_ASSERT_IMAGE_TYPE(mask2, PS_TYPE_MASK, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(mask2, mask1, NULL); + } + PS_ASSERT_INT_NONNEGATIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); + if (isfinite(badFrac)) { + PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL); + PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(badFrac, 1.0, NULL); + } + + int numCols = mask1->numCols, numRows = mask1->numRows; // Size of the images + + // Dereference inputs for convenience + psMaskType **data1 = mask1->data.PS_TYPE_MASK_DATA; + psMaskType **data2 = NULL; + if (mask2) { + data2 = mask2->data.PS_TYPE_MASK_DATA; + } + + // First, a pass through to determine the fraction of bad pixels + if (isfinite(badFrac) && badFrac != 1.0) { + int numBad = 0; // Number of bad pixels + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (data1[y][x] & maskVal) { + numBad++; + continue; + } + if (data2 && data2[y][x] & maskVal) { + numBad++; + } + } + } + if (numBad > badFrac * numCols * numRows) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Fraction of bad pixels (%d/%d=%f) exceeds limit (%f)\n", + numBad, numCols * numRows, (float)numBad/(float)(numCols * numRows), badFrac); + return NULL; + } + } + + // Worried about the masks for bad pixels and bad stamps colliding, so make our own mask + psImage *mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); // The global mask + psImageInit(mask, 0); + psMaskType **maskData = mask->data.PS_TYPE_MASK_DATA; // Dereference for convenience + + // Block out a border around the edge of the image + + // Bottom stripe + for (int y = 0; y < PS_MIN(size + footprint, numRows); y++) { + for (int x = 0; x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + // Either side + for (int y = PS_MIN(size + footprint, numRows); y < numRows - size - footprint; y++) { + for (int x = 0; x < PS_MIN(size + footprint, numCols); x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + for (int x = PS_MAX(numCols - size - footprint, 0); x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + // Top stripe + for (int y = PS_MAX(numRows - size - footprint, 0); y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BORDER; + } + } + + // XXX Could do something smarter here --- we will get images that are predominantly masked (where the + // skycell isn't overlapped by a large fraction by the observation), so that convolving around every bad + // pixel is wasting time. As a first cut, I've put in a check on the fraction of bad pixels, but we could + // imagine looking for the edge of big regions and convolving just at the edge. As a second cut, allow + // use of FFT convolution. + + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (data1[y][x] & maskVal) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_1; + } + if (data2 && data2[y][x] & maskVal) { + maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_2; + } + } + } + + // Block out the entire stamp footprint around bad input pixels. + + // We want to block out with the CONVOLVE mask anything that would be bad if we convolved with a bad + // reference pixel (within 'size'). Then we want to block out with the FOOTPRINT mask everything within a + // footprint's distance of those (within 'footprint'). + + if (useFFT) { + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_BAD_1 | PM_SUBTRACTION_MASK_BAD_2, + PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1, + -size, size, -size, size, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_BAD_2, + PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2, + -size, size, -size, size, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_1, + PM_SUBTRACTION_MASK_FOOTPRINT_1, + -footprint, footprint, -footprint, footprint, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskFFT(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_2, + PM_SUBTRACTION_MASK_FOOTPRINT_2, + -footprint, footprint, -footprint, footprint, 0.5)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + } else { + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_BAD_1, + PM_SUBTRACTION_MASK_CONVOLVE_1 | PM_SUBTRACTION_MASK_FOOTPRINT_1, + -size, size, -size, size)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_BAD_2, + PM_SUBTRACTION_MASK_CONVOLVE_2 | PM_SUBTRACTION_MASK_FOOTPRINT_2, + -size, size, -size, size)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_1, + PM_SUBTRACTION_MASK_FOOTPRINT_1, + -footprint, footprint, -footprint, footprint)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 1."); + psFree(mask); + return NULL; + } + if (!psImageConvolveMaskDirect(mask, mask, PM_SUBTRACTION_MASK_CONVOLVE_2, + PM_SUBTRACTION_MASK_FOOTPRINT_2, + -footprint, footprint, -footprint, footprint)) { + psError(PS_ERR_UNKNOWN, false, "Unable to reconvolve bad pixels from mask 2."); + psFree(mask); + return NULL; + } + } + + return mask; +} + + +bool pmSubtractionBorder(psImage *image, psImage *weight, psImage *mask, + int size, psMaskType blank) +{ + PS_ASSERT_IMAGE_NON_NULL(image, false); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, false); + if (mask) { + PS_ASSERT_IMAGE_NON_NULL(mask, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(mask, image, false); + PS_ASSERT_IMAGE_TYPE(mask, PS_TYPE_MASK, false); + } + if (weight) { + PS_ASSERT_IMAGE_NON_NULL(weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image, false); + PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); + } + + int numCols = image->numCols, numRows = image->numRows; // Image dimensions + + for (int y = size; y < numRows - size; y++) { + for (int x = 0; x < size; x++) { + markBlank(image, mask, weight, x, y, blank); + } + for (int x = numCols - size; x < numCols; x++) { + markBlank(image, mask, weight, x, y, blank); + } + } + for (int y = 0; y < size; y++) { + for (int x = 0; x < numCols; x++) { + markBlank(image, mask, weight, x, y, blank); + } + } + for (int y = numRows - size; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + markBlank(image, mask, weight, x, y, blank); + } + } + + return true; +} Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMask.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMask.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMask.h (revision 16477) @@ -0,0 +1,25 @@ +#ifndef PM_SUBTRACTION_MASK_H +#define PM_SUBTRACTION_MASK_H + +#include + +/// Generate a mask for use in the subtraction process +psImage *pmSubtractionMask(const psImage *refMask, ///< Mask for the reference image (will be convolved) + const psImage *inMask, ///< Mask for the input image, or NULL + psMaskType maskVal, ///< Value to mask out + int size, ///< Half-size of the kernel (pmSubtractionKernels.size) + int footprint, ///< Half-size of the kernel footprint + float badFrac, ///< Maximum fraction of bad input pixels to accept + bool useFFT ///< Use FFT to do convolution? + ); + +/// Mark the non-convolved part of the image as blank +bool pmSubtractionBorder(psImage *image,///< Image + psImage *weight, ///< Weight map (or NULL) + psImage *mask, ///< Mask (or NULL) + int size, ///< Kernel half-size + psMaskType blank ///< Mask value for blank regions + ); + + +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMatch.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMatch.c (revision 16477) @@ -0,0 +1,654 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmSubtractionParams.h" +#include "pmSubtractionKernels.h" +#include "pmSubtractionStamps.h" +#include "pmSubtractionEquation.h" +#include "pmSubtraction.h" +#include "pmSubtractionMask.h" +#include "pmSubtractionMatch.h" + +#define KERNEL_MOSAIC 2 // Half-number of kernel instances in the mosaic image + +static bool useFFT = true; // Do convolutions using FFT + + +#define TESTING +//#define TESTING_MEMORY + +// Output memory usage information +static void memCheck(const char *where) +{ +#ifdef TESTING_MEMORY + psMemBlock **leaks = NULL; + int numLeaks = psMemCheckLeaks(0, &leaks, NULL, true); + size_t largestSize = 0; + psMemId largest = 0; + size_t totalSize = 0; + for (int i = 0; i < numLeaks; i++) { + psMemBlock *mb = leaks[i]; + totalSize += mb->userMemorySize; + if (mb->userMemorySize > largestSize) { + largestSize = mb->userMemorySize; + largest = mb->id; + } + } + psFree(leaks); + fprintf(stderr, "%s:\n", where); + fprintf(stderr, " Memory in use: %zd\n", totalSize); + fprintf(stderr, " Largest block: %ld\n", largest); + fprintf(stderr, " sbrk(): %zd\n", (size_t)sbrk(0)); +#endif + return; +} + + +static bool getStamps(pmSubtractionStampList **stamps, // Stamps to read + const pmReadout *ro1, // Readout 1 + const pmReadout *ro2, // Readout 2 + const psImage *subMask, // Mask for subtraction, or NULL + psImage *weight, // Weight map + const psRegion *region, // Region of interest, or NULL + float threshold, // Threshold for stamp finding + float stampSpacing, // Spacing between stamps + int size, // Kernel half-size + int footprint, // Convolution footprint for stamps + pmSubtractionMode mode // Mode for subtraction + ) +{ + psTrace("psModules.imcombine", 3, "Finding stamps...\n"); + *stamps = pmSubtractionStampsFind(*stamps, ro1->image, subMask, region, threshold, footprint, + stampSpacing, mode); + if (!*stamps) { + psError(PS_ERR_UNKNOWN, false, "Unable to find stamps."); + return false; + } + + memCheck(" find stamps"); + + psTrace("psModules.imcombine", 3, "Extracting stamps...\n"); + if (!pmSubtractionStampsExtract(*stamps, ro1->image, ro2 ? ro2->image : NULL, weight, size)) { + psError(PS_ERR_UNKNOWN, false, "Unable to extract stamps."); + return false; + } + + memCheck(" extract stamps"); + + return true; +} + + + + + +bool pmSubtractionMatch(pmReadout *conv1, pmReadout *conv2, const pmReadout *ro1, const pmReadout *ro2, + int footprint, float regionSize, float stampSpacing, float threshold, + const psArray *sources, const char *stampsName, + pmSubtractionKernelsType type, int size, int spatialOrder, + const psVector *isisWidths, const psVector *isisOrders, + int inner, int ringsOrder, int binning, bool optimum, const psVector *optFWHMs, + int optOrder, float optThreshold, int iter, float rej, psMaskType maskBad, + psMaskType maskBlank, float badFrac, pmSubtractionMode mode) +{ + PS_ASSERT_PTR_NON_NULL(conv1, false); + if (mode == PM_SUBTRACTION_MODE_DUAL) { + PS_ASSERT_PTR_NON_NULL(conv2, false); + } + PS_ASSERT_PTR_NON_NULL(ro1, false); + PS_ASSERT_IMAGE_NON_NULL(ro1->image, false); + PS_ASSERT_IMAGE_TYPE(ro1->image, PS_TYPE_F32, false); + if (ro1->mask) { + PS_ASSERT_IMAGE_NON_NULL(ro1->mask, false); + PS_ASSERT_IMAGE_TYPE(ro1->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->mask, ro1->image, false); + } + if (ro1->weight) { + PS_ASSERT_IMAGE_NON_NULL(ro1->weight, false); + PS_ASSERT_IMAGE_TYPE(ro1->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->weight, ro1->image, false); + } + if (ro2) { + PS_ASSERT_IMAGE_NON_NULL(ro2->image, false); + PS_ASSERT_IMAGE_TYPE(ro2->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->image, ro1->image, false); + if (ro2->mask) { + PS_ASSERT_IMAGE_NON_NULL(ro2->mask, false); + PS_ASSERT_IMAGE_TYPE(ro2->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->mask, ro1->image, false); + } + if (ro2->weight) { + PS_ASSERT_IMAGE_NON_NULL(ro2->weight, false); + PS_ASSERT_IMAGE_TYPE(ro2->weight, PS_TYPE_F32, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(ro2->weight, ro1->image, false); + } + } else if (!stampsName && !sources) { + psError(PS_ERR_UNEXPECTED_NULL, true, + "A list of sources is required for convolving to a target PSF."); + return false; + } + PS_ASSERT_INT_NONNEGATIVE(footprint, false); + // regionSize can be just about anything (except maybe negative, but it can be NAN) + PS_ASSERT_FLOAT_LARGER_THAN(stampSpacing, 0.0, false); + // Don't care what threshold is + if (sources) { + PS_ASSERT_ARRAY_NON_NULL(sources, false); + } + // stampsName may be anything + // We'll check kernel type when we allocate the kernels + PS_ASSERT_INT_POSITIVE(size, false); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, false); + if (isisWidths || isisOrders) { + PS_ASSERT_VECTOR_NON_NULL(isisWidths, false); + PS_ASSERT_VECTOR_TYPE(isisWidths, PS_TYPE_F32, false); + PS_ASSERT_VECTOR_NON_NULL(isisOrders, false); + PS_ASSERT_VECTOR_TYPE(isisOrders, PS_TYPE_S32, false); + PS_ASSERT_VECTORS_SIZE_EQUAL(isisWidths, isisOrders, false); + } + PS_ASSERT_INT_NONNEGATIVE(inner, false); + PS_ASSERT_INT_NONNEGATIVE(ringsOrder, false); + PS_ASSERT_INT_POSITIVE(binning, false); + if (optimum) { + PS_ASSERT_VECTOR_NON_NULL(optFWHMs, false); + PS_ASSERT_INT_NONNEGATIVE(optOrder, false); + PS_ASSERT_FLOAT_LARGER_THAN(optThreshold, 0.0, false); + PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(optThreshold, 1.0, false); + } + PS_ASSERT_INT_POSITIVE(iter, false); + PS_ASSERT_FLOAT_LARGER_THAN(rej, 0.0, false); + // Don't care about maskBad + // Don't care about maskBlank + if (isfinite(badFrac)) { + PS_ASSERT_FLOAT_LARGER_THAN(badFrac, 0.0, NULL); + PS_ASSERT_FLOAT_LESS_THAN_OR_EQUAL(badFrac, 1.0, NULL); + } + + // If the stamp footprint is smaller than the kernel size, then we won't get much signal in the outer + // parts of the kernel, which can result in bad matching artifacts. + if (footprint < size) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Stamp footprint (%d) should be larger than or equal to the kernel size (%d)", + footprint, size); + return false; + } + + // Reset the output readout, just in case + if (conv1->image) { + psFree(conv1->image); + conv1->image = NULL; + } + if (conv1->mask) { + psFree(conv1->mask); + conv1->mask = NULL; + } + if (conv1->weight) { + psFree(conv1->weight); + conv1->weight = NULL; + } + + // Where does our weight map come from? + psImage *weight = NULL; // Weight image to use + if (ro1->weight && ro2 && ro2->weight) { + weight = (psImage*)psBinaryOp(NULL, ro1->weight, "+", ro2->weight); + } else if (ro1->weight) { + weight = psMemIncrRefCounter(ro1->weight); + } else if (ro2) { + if (ro2->weight) { + weight = psMemIncrRefCounter(ro2->weight); + } else { + weight = (psImage*)psBinaryOp(NULL, ro1->image, "+", ro2->image); + } + } else { + weight = psMemIncrRefCounter(ro1->image); + } + + // Putting important variable declarations here, since they are freed after a "goto" if there is an error. + psImage *subMask = NULL; // Mask for subtraction + psRegion *region = NULL; // Iso-kernel region + psString regionString = NULL; // String for region + pmSubtractionStampList *stamps = NULL; // Stamps for matching PSF + pmSubtractionKernels *kernels = NULL; // Kernel basis functions + int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Image dimensions + + memCheck("start"); + + subMask = pmSubtractionMask(ro1->mask, ro2 ? ro2->mask : NULL, maskBad, size, footprint, + badFrac, useFFT); + if (!subMask) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate subtraction mask."); + psFree(weight); + return false; + } + + memCheck("mask"); + + // Get region of interest + int xRegions = 1, yRegions = 1; // Number of iso-kernel regions + float xRegionSize = 0, yRegionSize = 0; // Size of iso-kernel regions + if (isfinite(regionSize) && regionSize != 0.0) { + xRegions = numCols / regionSize + 1; + yRegions = numRows / regionSize + 1; + xRegionSize = (float)numCols / (float)xRegions; + yRegionSize = (float)numRows / (float)yRegions; + region = psRegionAlloc(NAN, NAN, NAN, NAN); + } + + // Iterate over iso-kernel regions + for (int j = 0; j < yRegions; j++) { + for (int i = 0; i < xRegions; i++) { + psTrace("psModules.imcombine", 1, "Subtracting region %d of %d...\n", + j * xRegions + i + 1, xRegions * yRegions); + if (region) { + *region = psRegionSet((int)(i * xRegionSize), (int)((i + 1) * xRegionSize), + (int)(j * yRegionSize), (int)((j + 1) * yRegionSize)); + psFree(regionString); + regionString = psRegionToString(*region); + psTrace("psModules.imcombine", 3, "Iso-kernel region: %s out of %d,%d\n", + regionString, numCols, numRows); + } + + if (sources) { + stamps = pmSubtractionStampsSetFromSources(sources, subMask, region, footprint, + stampSpacing, mode); + } else if (stampsName && strlen(stampsName) > 0) { + stamps = pmSubtractionStampsSetFromFile(stampsName, ro1->image, subMask, region, footprint, + stampSpacing, mode); + } + + // We get the stamps here; we will also attempt to get stamps at the first iteration, but it + // doesn't matter. + if (!getStamps(&stamps, ro1, ro2, subMask, weight, NULL, threshold, stampSpacing, + size, footprint, mode)) { + goto MATCH_ERROR; + } + + if (mode == PM_SUBTRACTION_MODE_UNSURE || mode == PM_SUBTRACTION_MODE_TARGET) { + pmSubtractionMode newMode = pmSubtractionOrder(stamps, footprint); // Subtraction mode + switch (newMode) { + case PM_SUBTRACTION_MODE_1: + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 1 to match image 2."); + break; + case PM_SUBTRACTION_MODE_2: + if (mode == PM_SUBTRACTION_MODE_TARGET) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Input PSF is larger than target PSF --- can't match image."); + goto MATCH_ERROR; + } + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Convolving image 2 to match image 1."); + break; + default: + psError(PS_ERR_UNKNOWN, false, "Unable to determine subtraction order."); + goto MATCH_ERROR; + } + mode = newMode; + } + + // Define kernel basis functions + if (optimum && (type == PM_SUBTRACTION_KERNEL_ISIS || type == PM_SUBTRACTION_KERNEL_GUNK)) { + kernels = pmSubtractionKernelsOptimumISIS(type, size, inner, spatialOrder, optFWHMs, optOrder, + stamps, footprint, optThreshold, mode); + if (!kernels) { + psErrorClear(); + psWarning("Unable to derive optimum ISIS kernel --- switching to default."); + } + } + if (kernels == NULL) { + // Not an ISIS/GUNK kernel, or the optimum kernel search failed + kernels = pmSubtractionKernelsGenerate(type, size, spatialOrder, isisWidths, isisOrders, + inner, binning, ringsOrder, mode); + } + + // Add analysis metadata + { + psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, + PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); + if (conv2) { + psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, + PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Subtraction kernels", kernels); + } + psMetadataAddS32(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, + PS_META_DUPLICATE_OK, "Subtraction kernels", mode); + if (conv2) { + psMetadataAddS32(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_MODE, + PS_META_DUPLICATE_OK, "Subtraction kernels", mode); + } + psRegion *subRegion; + if (region) { + subRegion = psMemIncrRefCounter(region); + } else { + subRegion = psRegionAlloc(0, numCols, 0, numRows); + } + psMetadataAddPtr(conv1->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, + PS_DATA_REGION | PS_META_DUPLICATE_OK, + "Region over which subtraction was performed", subRegion); + if (conv2) { + psMetadataAddPtr(conv2->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, + PS_DATA_REGION | PS_META_DUPLICATE_OK, + "Region over which subtraction was performed", subRegion); + } + psFree(subRegion); + } + + memCheck("kernels"); + + int numRejected = -1; // Number of rejected stamps in each iteration + for (int k = 0; k < iter && numRejected != 0; k++) { + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Iteration %d.", k); + + if (!getStamps(&stamps, ro1, ro2, subMask, weight, region, threshold, stampSpacing, + size, footprint, mode)) { + goto MATCH_ERROR; + } + + psTrace("psModules.imcombine", 3, "Calculating equation...\n"); + if (!pmSubtractionCalculateEquation(stamps, kernels)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + + memCheck(" calculate equation"); + + psTrace("psModules.imcombine", 3, "Solving equation...\n"); + + if (!pmSubtractionSolveEquation(kernels, stamps)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + + memCheck(" solve equation"); + + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + goto MATCH_ERROR; + } + + memCheck(" calculate deviations"); + + psTrace("psModules.imcombine", 3, "Rejecting stamps...\n"); + numRejected = pmSubtractionRejectStamps(stamps, deviations, subMask, rej, footprint); + if (numRejected < 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to reject stamps."); + psFree(deviations); + goto MATCH_ERROR; + } + psFree(deviations); + + memCheck(" reject stamps"); + } + + if (numRejected > 0) { + psTrace("psModules.imcombine", 3, "Solving equation...\n"); + if (!pmSubtractionSolveEquation(kernels, stamps)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + goto MATCH_ERROR; + } + psVector *deviations = pmSubtractionCalculateDeviations(stamps, kernels); // Stamp deviations + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + goto MATCH_ERROR; + } + pmSubtractionRejectStamps(stamps, deviations, subMask, NAN, footprint); + psFree(deviations); + } + psFree(stamps); + stamps = NULL; + + memCheck("solution"); + + { + psTrace("psModules.imcombine", 2, "Generating diagnostics...\n"); + // Generate image with convolution kernels + int fullSize = 2 * size + 1 + 1; // Full size of kernel + psImage *convKernels = psImageAlloc(5 * fullSize - 1, 5 * fullSize - 1, PS_TYPE_F32); + psImageInit(convKernels, NAN); + for (int j = -KERNEL_MOSAIC; j <= KERNEL_MOSAIC; j++) { + for (int i = -KERNEL_MOSAIC; i <= KERNEL_MOSAIC; i++) { + psImage *kernel = pmSubtractionKernelImage(kernels, (float)i / (float)KERNEL_MOSAIC, + (float)j / (float)KERNEL_MOSAIC, + false); // Image of the kernel + if (!kernel) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate kernel image."); + psFree(convKernels); + goto MATCH_ERROR; + } + + if (psImageOverlaySection(convKernels, kernel, (i + KERNEL_MOSAIC) * fullSize, + (j + KERNEL_MOSAIC) * fullSize, "=") == 0) { + psError(PS_ERR_UNKNOWN, false, "Unable to overlay kernel image."); + psFree(kernel); + psFree(convKernels); + goto MATCH_ERROR; + } + psFree(kernel); + } + } + + psString comment = NULL; // Comment for metadata + psStringAppend(&comment, "Subtraction kernel for region %s", regionString); + psMetadataAddImage(conv1->analysis, PS_LIST_TAIL, "SUBTRACTION.KERNEL.IMAGE", + PS_META_DUPLICATE_OK, comment, convKernels); + psFree(comment); + psFree(convKernels); + } + +#if 0 + { + // Generate images of the kernel components + psMetadata *header = psMetadataAlloc(); // Header + for (int i = 0; i < solution->n; i++) { + psString name = NULL; // Header keyword + psStringAppend(&name, "SOLN%04d", i); + psMetadataAddF64(header, PS_LIST_TAIL, name, 0, NULL, solution->data.F64[i]); + psFree(name); + } + psArray *kernelImages = pmSubtractionKernelSolutions(solution, kernels, 0.0, 0.0); + psFits *kernelFile = psFitsOpen("kernels.fits", "w"); + (void)psFitsWriteImageCube(kernelFile, header, kernelImages, NULL); + psFitsClose(kernelFile); + psFree(kernelImages); + psFree(header); + } +#endif + + memCheck("diag outputs"); + + psTrace("psModules.imcombine", 2, "Convolving...\n"); + if (!pmSubtractionConvolve(conv1, conv2, ro1, ro2, subMask, maskBlank, region, kernels, useFFT)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve image."); + goto MATCH_ERROR; + } + psFree(kernels); + kernels = NULL; + + // There is data in the readout now + conv1->data_exists = true; + if (conv1->parent) { + conv1->parent->data_exists = true; + conv1->parent->parent->data_exists = true; + } + if (mode == PM_SUBTRACTION_MODE_DUAL) { + conv2->data_exists = true; + if (conv2->parent) { + conv2->parent->data_exists = true; + conv2->parent->parent->data_exists = true; + } + } + } + } + psFree(region); + region = NULL; + psFree(regionString); + regionString = NULL; + psFree(subMask); + subMask = NULL; + psFree(weight); + weight = NULL; + + if (!pmSubtractionBorder(conv1->image, conv1->weight, conv1->mask, size, maskBlank)) { + psError(PS_ERR_UNKNOWN, false, "Unable to set border of convolved image."); + goto MATCH_ERROR; + } + + memCheck("convolution"); + + +#ifdef TESTING + { + psFits *fits = psFitsOpen("convolved1.fits", "w"); + psFitsWriteImage(fits, NULL, conv1->image, 0, NULL); + psFitsClose(fits); + + if (mode == PM_SUBTRACTION_MODE_DUAL) { + psFits *fits = psFitsOpen("convolved2.fits", "w"); + psFitsWriteImage(fits, NULL, conv2->image, 0, NULL); + psFitsClose(fits); + } + } +#endif + + return true; + +MATCH_ERROR: + psFree(region); + psFree(regionString); + psFree(subMask); + psFree(kernels); + psFree(stamps); + psFree(weight); + return false; +} + +// Calculate the second order moments for an image +static float subtractionOrderMoment(const psKernel *kernel, // Image for which to measure moments + int radius // Maximum radius + ) +{ + assert(kernel && kernel->kernel); + + int xMin = PS_MAX(kernel->xMin, -radius), xMax = PS_MIN(kernel->xMax, radius); // Bounds in x + int yMin = PS_MAX(kernel->yMin, -radius), yMax = PS_MIN(kernel->yMax, radius); // Bounds in y + + float xCentroid = 0.0, yCentroid = 0.0; // Centroid (first moment) + float sum = 0.0; // Sum (zero-th moment) + for (int y = yMin; y <= yMax; y++) { + for (int x = xMin; x <= xMax; x++) { + xCentroid += kernel->kernel[y][x] * x; + yCentroid += kernel->kernel[y][x] * y; + sum += kernel->kernel[y][x]; + } + } + xCentroid /= sum; + yCentroid /= sum; + + float eta20 = 0.0, eta02 = 0.0; // Second moments + for (int y = yMin; y <= yMax; y++) { + float yDiff = y - yCentroid; + for (int x = xMin; x <= xMax; x++) { + float xDiff = x - xCentroid; + eta20 += PS_SQR(xDiff) * kernel->kernel[y][x]; + eta02 += PS_SQR(yDiff) * kernel->kernel[y][x]; + } + } + + // Normalise to calculate the scale-invariant + float sum2 = PS_SQR(sum); + eta20 /= sum2; + eta02 /= sum2; + // eta11 /= sum2; + + return eta20 + eta02; +} + +#if 0 +// Calculate the deviations for a particular subtraction order +static psVector *subtractionOrderDeviation(float *sumKernel, // Sum of the kernel + pmSubtractionStampList *stamps, // Stamps to convolve + const pmSubtractionKernels *kernels, // Kernel basis functions + int footprint, // Stamp footprint + pmSubtractionMode mode // Mode of subtraction + ) +{ + assert(stamps); + assert(footprint >= 0); + assert(mode == PM_SUBTRACTION_MODE_1 || mode == PM_SUBTRACTION_MODE_2); + + if (!pmSubtractionCalculateEquation(stamps, kernels, footprint, mode)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return NULL; + } + + psVector *solution = pmSubtractionSolveEquation(NULL, stamps); + if (!solution) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate least-squares equation."); + return NULL; + } + + if (sumKernel) { + float sum = 0.0; // Sum of the kernel + psImage *image = pmSubtractionKernelImage(solution, kernels, 0.0, 0.0); // Image of kernel + for (int y = 0; y < image->numRows; y++) { + for (int x = 0; x < image->numCols; x++) { + sum += image->data.F32[y][x]; + } + } + psFree(image); + *sumKernel = sum; + } + + psVector *deviations = pmSubtractionCalculateDeviations(stamps, solution, footprint, kernels, mode); + psFree(solution); + if (!deviations) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate deviations."); + return NULL; + } + + return deviations; +} +#endif + +pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, int radius) +{ + PS_ASSERT_INT_POSITIVE(radius, PM_SUBTRACTION_MODE_ERR); + + psVector *mask = psVectorAlloc(stamps->num, PS_TYPE_MASK); // Mask for stamps + psVector *moments = psVectorAlloc(stamps->num, PS_TYPE_F32); // Moments + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) { + mask->data.PS_TYPE_MASK_DATA[i] = 0xff; + continue; + } + mask->data.PS_TYPE_MASK_DATA[i] = 0; + moments->data.F32[i] = subtractionOrderMoment(stamp->image1, radius) / + subtractionOrderMoment(stamp->image2, radius); + } + + psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); + if (!psVectorStats(stats, moments, NULL, mask, 0xff)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate statistics for moments ratio."); + psFree(mask); + psFree(moments); + psFree(stats); + return PM_SUBTRACTION_MODE_ERR; + } + psFree(moments); + psFree(mask); + + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Median ratio of second moments: %lf", stats->robustMedian); + pmSubtractionMode mode = (stats->robustMedian <= 1.0 ? PM_SUBTRACTION_MODE_1 : PM_SUBTRACTION_MODE_2); + psFree(stats); + + return mode; +} + + + Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMatch.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMatch.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionMatch.h (revision 16477) @@ -0,0 +1,56 @@ +#ifndef PM_SUBTRACTION_MATCH_H +#define PM_SUBTRACTION_MATCH_H + +#include + +#include +#include +#include +#include + +#define PM_SUBTRACTION_ANALYSIS_KERNEL "SUBTRACTION.KERNEL" // Name of kernel in the analysis metadata +#define PM_SUBTRACTION_ANALYSIS_MODE "SUBTRACTION.MODE" // Name of subtraction mode in the analysis metadata +#define PM_SUBTRACTION_ANALYSIS_REGION "SUBTRACTION.REGION" // Name of subtraction region in the analysis MD + + +/// Match two images +bool pmSubtractionMatch(pmReadout *conv1, ///< Output convolved data for image 1 + pmReadout *conv2, ///< Output convolved data for image 2 + const pmReadout *ro1, ///< Image 1 + const pmReadout *ro2, ///< Image 2 + // Stamp parameters + int footprint, ///< Stamp half-size + float regionSize, ///< Typical size of iso-kernel regions + float stampSpacing, ///< Typical spacing between stamps + float threshold, ///< Threshold for stamps + const psArray *sources, ///< Sources for stamps + const char *stampsName, ///< Filename for stamps + // Kernel parameters + pmSubtractionKernelsType type, ///< Kernel type + int size, ///< Kernel half-size + int order, ///< Spatial polynomial order + const psVector *widths, ///< ISIS Gaussian widths + const psVector *orders, ///< ISIS Polynomial orders + int inner, ///< Inner radius for various kernel types + int ringsOrder, ///< RINGS polynomial order + int binning, ///< SPAM kernel binning + bool optimum, ///< Search for optimum ISIS kernel? + const psVector *optFWHMs, ///< FWHMs for optimum search + int optOrder, ///< Maximum order for optimum search + float optThreshold, ///< Threshold for optimum search (0..1) + // Operational parameters + int iter, ///< Rejection iterations + float rej, ///< Rejection threshold + psMaskType maskBad, ///< Value to mask + psMaskType maskBlank, ///< Mask for blank region + float badFrac, ///< Maximum fraction of bad input pixels to accept + pmSubtractionMode mode ///< Mode of subtraction + ); + +/// Determine which image to convolve +pmSubtractionMode pmSubtractionOrder(pmSubtractionStampList *stamps, ///< Stamps that have been extracted + int footprint ///< Stamp half-size + ); + + +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionParams.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionParams.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionParams.c (revision 16477) @@ -0,0 +1,510 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmSubtractionStamps.h" +#include "pmSubtraction.h" +#include "pmSubtractionParams.h" + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +#if 0 +// Convolve the reference stamp by the kernel +static psKernel *convolveStamp(const pmSubtractionStamp *stamp, // Stamp to be convolved + const psKernel *kernel, // Kernel by which to convolve + int footprint // Size of area to be convolved + ) +{ + psKernel *convolution = psKernelAlloc(-footprint, footprint, -footprint, footprint); // Result + psKernel *reference = stamp->reference; // Reference stamp, to be convolved + + // Range of kernel + int uMin = kernel->xMin; + int uMax = kernel->xMax; + int vMin = kernel->yMin; + int vMax = kernel->yMax; + + for (int y = -footprint; y <= footprint; y++) { + psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution + for (int x = -footprint; x <= footprint; x++, conv++) { + *conv = 0.0; + + int xStart = x + uMin; // Start index for convolution + for (int v = vMin; v <= vMax; v++) { + psF32 *ref = &reference->kernel[y + v][xStart]; // Dereference reference image + psF32 *krnl = &kernel->kernel[v][uMin]; // Dereference kernel in x + for (int u = uMin; u <= uMax; u++, ref++, krnl++) { + *conv += *ref * *krnl; + } + } + } + } + + return convolution; +} +#endif + +/// Select the appropriate convolution, given the kernel basis function and subtraction mode +static inline psKernel *selectConvolution(const pmSubtractionStamp *stamp, // Stamp + int kernelIndex, // Index for kernel component + pmSubtractionMode mode // Mode of subtraction + ) +{ + switch (mode) { + case PM_SUBTRACTION_MODE_1: + return stamp->convolutions1->data[kernelIndex]; + case PM_SUBTRACTION_MODE_2: + return stamp->convolutions2->data[kernelIndex]; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + return NULL; // Unreached +} + +// Accumulate cross-term sums for a stamp +static void accumulateCross(double *sumI, // Sum of I(x)/sigma(x)^2 + double *sumII, // Sum of I(x)^2/sigma(x)^2 + double *sumIC, // Sum of I(x)conv(x)/sigma(x)^2 + const pmSubtractionStamp *stamp, // Stamp with weight + const psKernel *target, // Target stamp + int kernelIndex, // Index for kernel component + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction + ) +{ + psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest + + for (int y = -footprint; y <= footprint; y++) { + psF32 *in = &target->kernel[y][-footprint]; // Dereference input + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution + for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { + double temp = *in / *wt; // Temporary product + *sumI += temp; + *sumII += *in * temp; + *sumIC += *conv * temp; + } + } + return; +} + +// Accumulate convolution sums for a stamp +static void accumulateConvolutions(double *sumC, // Sum of conv(x)/sigma(x)^2 + double *sumCC, // Sum of conv(x)^2/sigma(x)^2 + const pmSubtractionStamp *stamp, // Stamp with input and weight + int kernelIndex, // Index for kernel component + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction + ) +{ + psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest + + for (int y = -footprint; y <= footprint; y++) { + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution + for (int x = -footprint; x <= footprint; x++, wt++, conv++) { + double convNoise = *conv / *wt; // Temporary product + *sumC += convNoise; + *sumCC += *conv * convNoise; + } + } + return; +} + +static double accumulateChi2(const psKernel *target, // Target stamp + pmSubtractionStamp *stamp, // Stamp with weight + int kernelIndex, // Index for kernel component + double coeff, // Coefficient of convolution + double bg, // Background term + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction + ) +{ + double chi2 = 0.0; + psKernel *weight = stamp->weight; // Weight, sigma(x)^2 + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest + + for (int y = -footprint; y <= footprint; y++) { + psF32 *in = &target->kernel[y][-footprint]; // Dereference input + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution + for (int x = -footprint; x <= footprint; x++, in++, wt++, conv++) { + chi2 += PS_SQR(*in - bg - coeff * *conv) / *wt; + } + } + + return chi2; +} + +// Return the initial value of chi^2 +static double initialChi2(const psKernel *target, // Target stamp + const pmSubtractionStamp *stamp, // Stamp with weight + int footprint, // Size of convolution + pmSubtractionMode mode // Mode of subtraction + ) +{ + psKernel *weight = stamp->weight; // Weight map + psKernel *source; // Source stamp + switch (mode) { + case PM_SUBTRACTION_MODE_1: + source = stamp->image1; + break; + case PM_SUBTRACTION_MODE_2: + source = stamp->image2; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + + double chi2 = 0.0; // Chi^2 + for (int y = -footprint; y <= footprint; y++) { + psF32 *in = &target->kernel[y][-footprint]; // Dereference input + psF32 *wt = &weight->kernel[y][-footprint]; // Dereference weight + psF32 *ref = &source->kernel[y][-footprint]; // Derference reference + for (int x = -footprint; x <= footprint; x++, in++, wt++, ref++) { + float diff = *in - *ref; // Temporary value + chi2 += PS_SQR(diff) / *wt; + } + } + + return chi2; +} + +// Subtract a convolution from the input +static void subtractConvolution(psKernel *target, // Target stamp + const pmSubtractionStamp *stamp, // Stamp with weight + int kernelIndex, // Index for kernel component + float coeff, // Coefficient of subtraction + float bg, // Background term + int footprint, // Size of region of interest + pmSubtractionMode mode // Mode of subtraction + ) +{ + psKernel *convolution = selectConvolution(stamp, kernelIndex, mode); // Convolution of interest + for (int y = -footprint; y <= footprint; y++) { + psF32 *in = &target->kernel[y][-footprint]; // Dereference input + psF32 *conv = &convolution->kernel[y][-footprint]; // Dereference convolution + for (int x = -footprint; x <= footprint; x++, in++, conv++) { + *in -= *conv * coeff + bg; + } + } + + return; +} + + +pmSubtractionKernels *pmSubtractionKernelsOptimumISIS(pmSubtractionKernelsType type, int size, int inner, + int spatialOrder, const psVector *fwhms, int maxOrder, + const pmSubtractionStampList *stamps, int footprint, + float tolerance, pmSubtractionMode mode) +{ + if (type != PM_SUBTRACTION_KERNEL_ISIS && type != PM_SUBTRACTION_KERNEL_GUNK) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Invalid kernel type: %x\n", type); + return NULL; + } + PS_ASSERT_INT_NONNEGATIVE(size, NULL); + PS_ASSERT_INT_NONNEGATIVE(inner, NULL); + PS_ASSERT_INT_NONNEGATIVE(spatialOrder, NULL); + PS_ASSERT_VECTOR_NON_NULL(fwhms, NULL); + PS_ASSERT_VECTOR_TYPE(fwhms, PS_TYPE_F32, NULL); + PS_ASSERT_INT_NONNEGATIVE(maxOrder, NULL); + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, NULL); + PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); + PS_ASSERT_FLOAT_LARGER_THAN(tolerance, 0.0, NULL); + + // Generate the kernels to test + int numGaussians = fwhms->n; // Number of Gaussians + int numKernels = numGaussians * (maxOrder + 1) * (maxOrder + 2) / 2; // Number of kernel components + psString params = NULL; // Parameter, for description + for (int i = 0; i < numGaussians; i++) { + psStringAppend(¶ms, "%.2f,", fwhms->data.F32[i]); + } + params[strlen(params) - 1] = '\0'; + + psVector *orders = psVectorAlloc(numGaussians, PS_TYPE_S32); // Polynomial orders + psVectorInit(orders, maxOrder); + pmSubtractionKernels *kernels = p_pmSubtractionKernelsRawISIS(size, spatialOrder, + fwhms, orders, mode); // Kernels + psFree(orders); + psFree(kernels->description); + kernels->description = NULL; + psStringAppend(&kernels->description, "OptISIS(%d,(%s),%d)", size, params, spatialOrder); + psFree(params); + + // Need to save the stamp inputs --- we're changing the values! + int numStamps = stamps->num; // Number of stamps + psArray *targets = psArrayAlloc(numStamps); // Deep copies of the targets + psVector *badStamps = psVectorAlloc(numStamps, PS_TYPE_U8); // Mark the bad stamps + psVectorInit(badStamps, 0); + for (int i = 0; i < numStamps; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (stamp->status != PM_SUBTRACTION_STAMP_CALCULATE && stamp->status != PM_SUBTRACTION_STAMP_USED) { + badStamps->data.U8[i] = 0xff; + continue; + } + psKernel *target; // Target image of interest + switch (mode) { + case PM_SUBTRACTION_MODE_1: + target = stamp->image2; + break; + case PM_SUBTRACTION_MODE_2: + target = stamp->image1; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + psImage *copy = psImageCopy(NULL, target->image, PS_TYPE_F32); // Copy of the image + targets->data[i] = psKernelAllocFromImage(copy, size + footprint, size + footprint); + psFree(copy); // Drop reference + } + + // Generate the convolutions, accumulate sums, and measure initial chi^2 + double sum1 = 0.0; // sum of 1/sigma(x,y)^2 + psVector *sumC = psVectorAlloc(numKernels, PS_TYPE_F64); // sum of R(x)*k(u)/sigma(x)^2 + psVector *sumCC = psVectorAlloc(numKernels, PS_TYPE_F64); // sum of [R(x)*k(u)]^2/sigma(x)^2 + psVectorInit(sumC, 0.0); + psVectorInit(sumCC, 0.0); + double lastChi2 = 0.0; // Chi^2 from last iteration + int numPixels = 0; // Number of pixels contributing to chi^2 + for (int i = 0; i < numStamps; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (badStamps->data.U8[i]) { + continue; + } + if (!pmSubtractionConvolveStamp(stamp, kernels, footprint)) { + psError(PS_ERR_UNKNOWN, false, "Unable to convolve stamp %d.", i); + psFree(targets); + psFree(kernels); + psFree(badStamps); + return NULL; + } + + // This sum is invariant to the kernel + psKernel *weight = stamp->weight; // Weight map for stamp + for (int v = -footprint; v <= footprint; v++) { + psF32 *wt = &weight->kernel[v][-footprint]; // Dereference weight map + for (int u = -footprint; u <= footprint; u++, wt++) { + sum1 += 1.0 / *wt; + } + } + if (!isfinite(sum1)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, + "Sum of 1/sigma^2 is non-finite for stamp %d (%d,%d)\n", + i, (int)stamp->x, (int)stamp->y); + psFree(targets); + psFree(kernels); + psFree(badStamps); + return NULL; + } + + for (int j = 0; j < numKernels; j++) { + accumulateConvolutions(&sumC->data.F64[j], &sumCC->data.F64[j], stamp, j, footprint, mode); + } + + lastChi2 += initialChi2(targets->data[i], stamp, footprint, mode); + numPixels += PS_SQR(2 * footprint + 1); + } + lastChi2 /= numPixels; + + // Rank the kernel components + psVector *ranking = psVectorAlloc(numKernels, PS_TYPE_S32); // Ranking of the kernel components + psVectorInit(ranking, -1); + int cutIndex = -1; // Index at which to cut off kernels + for (int iter = 0; iter < numKernels; iter++) { + int bestIndex = -1; // Index of best kernel component + double bestChi2 = INFINITY; // Value of best chi^2 + double bestCoeff = 0; // Value of best coefficient + double bestBG = 0; // Value of best background + + for (int i = 0; i < numKernels; i++) { + if (ranking->data.S32[i] >= 0) { + continue; + } + + double sumI = 0.0; // sum of I(x)/sigma(x)^2 + double sumII = 0.0; // sum of I(x)^2/sigma(x)^2 + double sumIC = 0.0; // sum of I(x)C(x)/sigma(x)^2 + + for (int j = 0; j < numStamps; j++) { + if (badStamps->data.U8[j]) { + continue; + } + pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest + accumulateCross(&sumI, &sumII, &sumIC, stamp, targets->data[j], i, footprint, mode); + } + + double invDet = 1.0 / (sum1 * sumCC->data.F64[i] - PS_SQR(sumC->data.F64[i])); // Determinant^-1 + double coeff = invDet * (sum1 * sumIC - sumC->data.F64[i] * sumI); // Coefficient for kernel + double bg = invDet * (sumCC->data.F64[i] * sumI - sumC->data.F64[i] * sumIC); // Background + + double chi2 = 0.0; // Chi^2 + for (int j = 0; j < numStamps; j++) { + if (badStamps->data.U8[j]) { + continue; + } + pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest + chi2 += accumulateChi2(targets->data[j], stamp, i, coeff, bg, footprint, mode); + } + + if (chi2 < bestChi2) { + bestIndex = i; + bestCoeff = coeff; + bestChi2 = chi2; + bestBG = bg; + } + + psTrace("psModules.imcombine", 8, "%d: %lf %lf %lf %lf %lf %lf\n", i, sum1, sumI, sumII, sumIC, + sumC->data.F64[i], sumCC->data.F64[i]); + psTrace("psModules.imcombine", 6, "%d: %lf %lf %lf\n", i, coeff, bg, chi2); + } + bestChi2 /= numPixels; + + if (bestIndex == -1) { + psError(PS_ERR_UNKNOWN, false, "Unable to find best kernel component in round %d.", iter); + psFree(targets); + psFree(sumC); + psFree(sumCC); + psFree(ranking); + psFree(kernels); + psFree(badStamps); + return NULL; + } + + // And the winner is.... + ranking->data.S32[bestIndex] = iter; + // Remove its contribution, and don't include it in the future. + for (int j = 0; j < numStamps; j++) { + if (badStamps->data.U8[j]) { + continue; + } + pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest + subtractConvolution(targets->data[j], stamp, bestIndex, bestCoeff, bestBG, footprint, mode); + } + + double diff = lastChi2 - bestChi2; // Difference in chi^2 between iterations + + psTrace("psModules.imcombine", 3, "The winner of round %d is %d (%f,%d,%d): %lf (%lf) --> %lf\n", + iter, bestIndex, kernels->widths->data.F32[bestIndex], kernels->u->data.S32[bestIndex], + kernels->v->data.S32[bestIndex], bestCoeff, diff, bestChi2); + + if (fabsf(diff) < tolerance) { + cutIndex = iter; + break; + } + + lastChi2 = bestChi2; + } + psFree(targets); + psFree(sumC); + psFree(sumCC); + + if (cutIndex < 0) { + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Unable to converge to tolerance %g\n", tolerance); + psFree(ranking); + psFree(kernels); + psFree(badStamps); + return NULL; + } + + int newSize = cutIndex + 1; // Size of new kernel basis set + psTrace("psModules.imcombine", 2, "Accepting %d kernels.\n", newSize); + psVector *uNew = psVectorAlloc(newSize, PS_TYPE_S32); + psVector *vNew = psVectorAlloc(newSize, PS_TYPE_S32); + psVector *widthsNew = psVectorAlloc(newSize, PS_TYPE_F32); + psArray *preCalcNew = psArrayAlloc(newSize); + psArray *convNew = psArrayAlloc(numStamps); + for (int i = 0; i < numStamps; i++) { + if (badStamps->data.U8[i]) { + continue; + } + convNew->data[i] = psArrayAlloc(newSize); + } + + for (int i = 0; i < numKernels; i++) { + int rank = ranking->data.S32[i]; // This kernel component's ranking + if (rank >= 0 && rank < newSize) { + uNew->data.S32[rank] = kernels->u->data.S32[i]; + vNew->data.S32[rank] = kernels->v->data.S32[i]; + widthsNew->data.F32[rank] = kernels->widths->data.F32[i]; + preCalcNew->data[rank] = psMemIncrRefCounter(kernels->preCalc->data[i]); + + for (int j = 0; j < numStamps; j++) { + if (badStamps->data.U8[j]) { + continue; + } + pmSubtractionStamp *stamp = stamps->stamps->data[j]; // Stamp of interest + psArray *convolutions = convNew->data[j]; // Convolutions for this stamp + convolutions->data[rank] = psMemIncrRefCounter(selectConvolution(stamp, i, mode)); + } + } + } + psFree(kernels->u); + psFree(kernels->v); + psFree(kernels->widths); + psFree(kernels->preCalc); + kernels->u = uNew; + kernels->v = vNew; + kernels->widths = widthsNew; + kernels->preCalc = preCalcNew; + kernels->num = newSize; + + for (int i = 0; i < numStamps; i++) { + if (badStamps->data.U8[i]) { + continue; + } + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + switch (mode) { + case PM_SUBTRACTION_MODE_1: + stamp->convolutions1 = convNew->data[i]; + stamp->convolutions2 = NULL; + break; + case PM_SUBTRACTION_MODE_2: + stamp->convolutions1 = NULL; + stamp->convolutions2 = convNew->data[i]; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + } + + psFree(badStamps); + psFree(ranking); + + // Maintain photometric scaling + if (type == PM_SUBTRACTION_KERNEL_ISIS) { + psKernel *subtract = kernels->preCalc->data[0]; // Kernel to subtract from the rest + for (int i = 1; i < newSize; i++) { + if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) { + psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest + psBinaryOp(kernel->image, kernel->image, "-", subtract->image); + } + } + } else if (type == PM_SUBTRACTION_KERNEL_GUNK) { + psStringPrepend(&kernels->description, "GUNK="); + psStringAppend(&kernels->description, "+POIS(%d,%d)", inner, spatialOrder); + + for (int i = 0; i < newSize; i++) { + if (kernels->u->data.S32[i] % 2 == 0 && kernels->v->data.S32[i] % 2 == 0) { + psKernel *kernel = kernels->preCalc->data[i]; // Kernel of interest + kernel->kernel[0][0] -= 1.0; + } + } + + if (!p_pmSubtractionKernelsAddGrid(kernels, numGaussians, inner)) { + psAbort("Should never get here."); + } + + kernels->type = PM_SUBTRACTION_KERNEL_GUNK; + } + + return kernels; +} Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionParams.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionParams.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionParams.h (revision 16477) @@ -0,0 +1,21 @@ +#ifndef PM_SUBTRACTION_PARAMS_H +#define PM_SUBTRACTION_PARAMS_H + +#include +#include +#include + +/// Generate a set of optimum kernels for ISIS (or GUNK) +pmSubtractionKernels *pmSubtractionKernelsOptimumISIS(pmSubtractionKernelsType type, ///< Kernel type + int size, ///< Half-size of kernel + int inner, ///< Inner radius for GUNK + int spatialOrder, ///< Spatial polynomial order + const psVector *fwhms, ///< Gaussian FWHMs to try + int maxOrder, ///< Maximum polynomial order + const pmSubtractionStampList *stamps, ///< Stamps + int footprint, ///< Convolution footprint for stamps + float tolerance, ///< Maximum difference in chi^2 + pmSubtractionMode mode // Mode for subtraction + ); + +#endif Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionStamps.c =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionStamps.c (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionStamps.c (revision 16477) @@ -0,0 +1,606 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +// All these includes required to get stamps out of an array of pmSources +#include "pmMoments.h" +#include "pmPeaks.h" +#include "pmResiduals.h" +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmSource.h" + + +#include "pmSubtraction.h" +#include "pmSubtractionStamps.h" + +#define STAMP_LIST_BUFFER 20 // Number of stamps to add to list at a time + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Private (file-static) functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +// Free function for pmSubtractionStampList +static void subtractionStampListFree(pmSubtractionStampList *list // Stamp list to free + ) +{ + psFree(list->stamps); + psFree(list->regions); + psFree(list->x); + psFree(list->y); + psFree(list->flux); +} + +// Free function for pmSubtractionStamp +static void subtractionStampFree(pmSubtractionStamp *stamp // Stamp to free + ) +{ + psFree(stamp->image1); + psFree(stamp->image2); + psFree(stamp->weight); + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + + psFree(stamp->matrix1); + psFree(stamp->matrix2); + psFree(stamp->matrixX); + psFree(stamp->vector1); + psFree(stamp->vector2); + +} + +// Is this region OK? +static bool checkStampRegion(int x, int y, // Coordinates of stamp + const psRegion *region // Region of interest + ) +{ + if (!region) { + return true; + } + return (x < region->x0 || x > region->x1 || y < region->y0 || y > region->y1) ? + false : true; +} + +// Is this position unmasked? +static bool checkStampMask(int x, int y, // Coordinates of stamp + const psImage *mask, // Mask + pmSubtractionMode mode // Mode for subtraction + ) +{ + if (!mask) { + return true; + } + if (x < 0 || x >= mask->numCols || y < 0 || y >= mask->numRows) { + return false; + } + + psMaskType maskVal = PM_SUBTRACTION_MASK_BORDER | PM_SUBTRACTION_MASK_REJ; // Mask value + switch (mode) { + case PM_SUBTRACTION_MODE_1: + case PM_SUBTRACTION_MODE_TARGET: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1; + break; + case PM_SUBTRACTION_MODE_2: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_2; + break; + case PM_SUBTRACTION_MODE_UNSURE: + case PM_SUBTRACTION_MODE_DUAL: + maskVal |= PM_SUBTRACTION_MASK_FOOTPRINT_1 | PM_SUBTRACTION_MASK_FOOTPRINT_2; + break; + default: + psAbort("Unsupported subtraction mode: %x", mode); + } + + return (mask->data.PS_TYPE_MASK_DATA[y][x] & maskVal) ? false : true; +} + +////////////////////////////////////////////////////////////////////////////////////////////////////////////// +// Public functions +////////////////////////////////////////////////////////////////////////////////////////////////////////////// + +pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, int numRows, const psRegion *region, + int footprint, float spacing) +{ + pmSubtractionStampList *list = psAlloc(sizeof(pmSubtractionStampList)); // Stamp list to return + psMemSetDeallocator(list, (psFreeFunc)subtractionStampListFree); + + // Get region in which to find stamps: [xMin:xMax,yMin:yMax] + int xMin = 0, xMax = numCols, yMin = 0, yMax = numRows; + if (region) { + xMin = PS_MAX(region->x0, xMin); + xMax = PS_MIN(region->x1, xMax); + yMin = PS_MAX(region->y0, yMin); + yMax = PS_MIN(region->y1, yMax); + } + int xSize = xMax - xMin, ySize = yMax - yMin; // Size of region of interest + int xStamps = (float)xSize / spacing + 1, yStamps = (float)ySize / spacing + 1; // Number of stamps + + list->num = xStamps * yStamps; + list->stamps = psArrayAlloc(list->num); + list->regions = psArrayAlloc(list->num); + + for (int y = 0, index = 0; y < yStamps; y++) { + int yStart = yMin + y * ((float)ySize / (float)(yStamps)); // Subregion starts here + int yStop = yMin + (y + 1) * ((float)ySize / (float)(yStamps)) - 1; // Subregion stops here + assert(yStart >= yMin && yStop < yMax); + + for (int x = 0; x < xStamps; x++, index++) { + int xStart = xMin + x * ((float)xSize / (float)(xStamps)); // Subregion starts here + int xStop = xMin + (x + 1) * ((float)xSize / (float)(xStamps)) - 1; // Subregion stops here + assert(xStart >= xMin && xStop < xMax); + + list->stamps->data[index] = pmSubtractionStampAlloc(); + psTrace("psModules.imcombine", 6, "Stamp region %d: [%d:%d,%d:%d]\n", + index, xStart, xStop, yStart, yStop); + list->regions->data[index] = psRegionAlloc(xStart, xStop, yStart, yStop); + } + } + + list->x = NULL; + list->y = NULL; + list->flux = NULL; + list->footprint = footprint; + + return list; +} + +pmSubtractionStamp *pmSubtractionStampAlloc(void) +{ + pmSubtractionStamp *stamp = psAlloc(sizeof(pmSubtractionStamp)); // Stamp to return + psMemSetDeallocator(stamp, (psFreeFunc)subtractionStampFree); + + stamp->x = NAN; + stamp->y = NAN; + stamp->flux = NAN; + stamp->xNorm = NAN; + stamp->yNorm = NAN; + stamp->status = PM_SUBTRACTION_STAMP_INIT; + + stamp->image1 = NULL; + stamp->image2 = NULL; + stamp->weight = NULL; + stamp->convolutions1 = NULL; + stamp->convolutions2 = NULL; + + stamp->matrix1 = NULL; + stamp->matrix2 = NULL; + stamp->matrixX = NULL; + stamp->vector1 = NULL; + stamp->vector2 = NULL; + + return stamp; +} + + +pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, const psImage *image, + const psImage *subMask, const psRegion *region, + float threshold, int footprint, float spacing, + pmSubtractionMode mode) +{ + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGE_TYPE(image, PS_TYPE_F32, NULL); + if (subMask) { + PS_ASSERT_IMAGE_NON_NULL(subMask, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL); + PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, NULL); + } + PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); + PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); + if (region) { + if (psRegionIsNaN(*region)) { + psString string = psRegionToString(*region); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) contains NAN values", string); + psFree(string); + return false; + } + if (region->x0 < 0 || region->x1 > image->numCols || + region->y0 < 0 || region->y1 > image->numRows) { + psString string = psRegionToString(*region); + psError(PS_ERR_BAD_PARAMETER_VALUE, true, "Input region (%s) does not fit in image (%dx%d)", + string, image->numCols, image->numRows); + psFree(string); + return false; + } + } + + int numRows = image->numRows, numCols = image->numCols; // Size of image + + if (!stamps) { + stamps = pmSubtractionStampListAlloc(numCols, numRows, region, footprint, spacing); + } + + int numStamps = stamps->num; // Number of stamp regions + int numFound = 0; // Number of stamps found + int numValid = 0; // Number of valid regions + for (int i = 0; i < numStamps; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + + // Only find a new stamp if we need to + if (stamp->status != PM_SUBTRACTION_STAMP_REJECTED && stamp->status != PM_SUBTRACTION_STAMP_INIT) { + continue; + } + numValid++; + + float xStamp = 0, yStamp = 0; // Coordinates of stamp + float fluxStamp = NAN; // Flux of stamp + bool goodStamp = false; // Found a good stamp? + + // A couple different ways of finding stamps: + if (stamps->x && stamps->y) { + // Get the next stamp from the list + psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Coordinate lists + psVector *fluxList = stamps->flux->data[i]; // List of stamp fluxes + + // Take stamp off the top of the (sorted) list + if (xList->n > 0) { + int index = xList->n - 1; // Index of new stamp + xStamp = xList->data.F32[index]; + yStamp = yList->data.F32[index]; + fluxStamp = fluxList->data.F32[index]; + + // Chop off the top of the list + xList->n = index; + yList->n = index; + fluxList->n = index; + + goodStamp = true; + } + } else { + // Use a simple method of automatically finding stamps --- take the highest pixel in the subregion + fluxStamp = threshold; + psRegion *subRegion = stamps->regions->data[i]; // Sub-region of interest + for (int y = subRegion->y0; y <= subRegion->y1; y++) { + for (int x = subRegion->x0; x <= subRegion->x1; x++) { + if (checkStampMask(x, y, subMask, mode) && image->data.F32[y][x] > fluxStamp) { + fluxStamp = image->data.F32[y][x]; + xStamp = x; + yStamp = y; + goodStamp = true; + } + } + } + } + + if (goodStamp) { + stamp->x = xStamp; + stamp->y = yStamp; + stamp->flux = fluxStamp; + + // Reset the postage stamps since we're making a new stamp + psFree(stamp->image1); + psFree(stamp->image2); + psFree(stamp->weight); + psFree(stamp->convolutions1); + psFree(stamp->convolutions2); + stamp->image1 = stamp->image2 = stamp->weight = NULL; + stamp->convolutions1 = stamp->convolutions2 = NULL; + + stamp->status = PM_SUBTRACTION_STAMP_FOUND; + numFound++; + psTrace("psModules.imcombine", 5, "Found stamp in subregion %d: %d,%d\n", + i, (int)stamp->x, (int)stamp->y); + } else { + stamp->status = PM_SUBTRACTION_STAMP_NONE; + } + } + + if (numValid > 0) { + psLogMsg("psModules.imcombine", PS_LOG_INFO, "Found %d stamps", numFound); + } + + return stamps; +} + + +pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, const psVector *y, const psVector *flux, + const psImage *image, const psImage *subMask, + const psRegion *region, int footprint, float spacing, + pmSubtractionMode mode) + +{ + PS_ASSERT_VECTOR_NON_NULL(x, NULL); + PS_ASSERT_VECTOR_TYPE(x, PS_TYPE_F32, NULL); + PS_ASSERT_VECTOR_NON_NULL(y, NULL); + PS_ASSERT_VECTOR_TYPE(y, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(y, x, NULL); + if (flux) { + PS_ASSERT_VECTOR_NON_NULL(flux, NULL); + PS_ASSERT_VECTOR_TYPE(flux, PS_TYPE_F32, NULL); + PS_ASSERT_VECTORS_SIZE_EQUAL(flux, x, NULL); + } else { + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + } + if (subMask) { + PS_ASSERT_IMAGE_NON_NULL(subMask, NULL); + PS_ASSERT_IMAGE_TYPE(subMask, PS_TYPE_MASK, NULL); + if (image) { + PS_ASSERT_IMAGE_NON_NULL(image, NULL); + PS_ASSERT_IMAGES_SIZE_EQUAL(image, subMask, NULL); + } + } + PS_ASSERT_FLOAT_LARGER_THAN(spacing, 0.0, NULL); + + int numStars = x->n; // Number of stars + pmSubtractionStampList *stamps = pmSubtractionStampListAlloc(subMask->numCols, subMask->numRows, + region, footprint, spacing); // Stamp list + int numStamps = stamps->num; // Number of stamps + + // Initialise the lists + stamps->x = psArrayAlloc(numStamps); + stamps->y = psArrayAlloc(numStamps); + stamps->flux = psArrayAlloc(numStamps); + for (int i = 0; i < numStamps; i++) { + stamps->x->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32); + stamps->y->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32); + stamps->flux->data[i] = psVectorAllocEmpty(STAMP_LIST_BUFFER, PS_TYPE_F32); + } + + // Put the stars into their appropriate subregions + for (int i = 0; i < numStars; i++) { + float xStamp = x->data.F32[i], yStamp = y->data.F32[i]; // Coordinates of stamp + int xPix = xStamp + 0.5, yPix = yStamp + 0.5; // Pixel coordinate of stamp + if (!checkStampRegion(xPix, yPix, region)) { + // It's not in the big region + psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because outside region", + xPix, yPix); + continue; + } + if (!checkStampMask(xPix, yPix, subMask, mode)) { + // Not a good stamp + psTrace("psModules.imcombine", 9, "Rejecting input stamp (%d,%d) because bad mask", + xPix, yPix); + continue; + } + + bool found = false; + for (int j = 0; j < numStamps && !found; j++) { + psRegion *subRegion = stamps->regions->data[j]; // Subregion of interest + if (checkStampRegion(xPix, yPix, subRegion)) { + psVector *xList = stamps->x->data[j], *yList = stamps->y->data[j]; // Pixel lists + psVector *fluxList = stamps->flux->data[j]; // Flux list + + int index = xList->n; // Index of new stamp candidate + + psVectorExtend(xList, STAMP_LIST_BUFFER, 1); + psVectorExtend(yList, STAMP_LIST_BUFFER, 1); + psVectorExtend(fluxList, STAMP_LIST_BUFFER, 1); + + xList->data.F32[index] = xStamp; + yList->data.F32[index] = yStamp; + + if (flux) { + fluxList->data.F32[index] = flux->data.F32[i]; + } else { + fluxList->data.F32[index] = image->data.F32[yPix][xPix]; + } + + found = true; + psTrace("psModules.imcombine", 9, "Putting input stamp (%d,%d) into subregion %d", + xPix, yPix, j); + } + } + + if (!found) { + psTrace("psModules.imcombine", 9, "Unable to find subregion for stamp (%d,%d)", + xPix, yPix); + } + } + + // Sort the list by flux, with the brightest last + for (int i = 0; i < numStamps; i++) { + psVector *xList = stamps->x->data[i], *yList = stamps->y->data[i]; // Pixel lists + psVector *fluxList = stamps->flux->data[i]; // Flux list + + psVector *indexes = psVectorSortIndex(NULL, fluxList); // Indices to sort flux + int num = indexes->n; // Number of candidate stamps in this subregion + + psVector *xSorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of x list + psVector *ySorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of y list + psVector *fluxSorted = psVectorAlloc(num, PS_TYPE_F32); // Sorted version of flux list + for (int j = 0; j < num; j++) { + int k = indexes->data.S32[j]; // Sorted index + xSorted->data.F32[j] = xList->data.F32[k]; + ySorted->data.F32[j] = yList->data.F32[k]; + fluxSorted->data.F32[j] = fluxList->data.F32[k]; + } + psFree(indexes); + + psFree(stamps->x->data[i]); + psFree(stamps->y->data[i]); + psFree(stamps->flux->data[i]); + + stamps->x->data[i] = xSorted; + stamps->y->data[i] = ySorted; + stamps->flux->data[i] = fluxSorted; + } + + + return stamps; +} + + +bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, psImage *image1, psImage *image2, + psImage *weight, int kernelSize) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + PS_ASSERT_IMAGE_NON_NULL(image1, false); + PS_ASSERT_IMAGE_TYPE(image1, PS_TYPE_F32, false); + if (image2) { + PS_ASSERT_IMAGE_NON_NULL(image2, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(image2, image1, false); + PS_ASSERT_IMAGE_TYPE(image2, PS_TYPE_F32, false); + } + PS_ASSERT_IMAGE_NON_NULL(weight, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(weight, image1, false); + PS_ASSERT_IMAGE_TYPE(weight, PS_TYPE_F32, false); + PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); + + int numCols = image1->numCols, numRows = image1->numRows; // Size of images + int size = kernelSize + stamps->footprint; // Size of postage stamps + + for (int i = 0; i < stamps->num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (!stamp || stamp->status != PM_SUBTRACTION_STAMP_FOUND) { + continue; + } + + if (isnan(stamp->xNorm)) { + stamp->xNorm = 2.0 * (stamp->x - (float)numCols/2.0) / (float)numCols; + } + if (isnan(stamp->yNorm)) { + stamp->yNorm = 2.0 * (stamp->y - (float)numRows/2.0) / (float)numRows; + } + + int x = stamp->x + 0.5, y = stamp->y + 0.5; // Stamp coordinates + if (x < size || x > numCols - size || y < size || y > numRows - size) { + psError(PS_ERR_UNKNOWN, false, "Stamp %d (%d,%d) is within the image border.\n", i, x, y); + return false; + } + + // Catch memory leaks --- these should have been freed and NULLed before + assert(stamp->image1 == NULL); + assert(stamp->image2 == NULL); + assert(stamp->weight == NULL); + + psRegion region = psRegionSet(x - size, x + size + 1, y - size, y + size + 1); // Region of interest + + psImage *sub1 = psImageSubset(image1, region); // Subimage with stamp + stamp->image1 = psKernelAllocFromImage(sub1, size, size); + psFree(sub1); // Drop reference + + if (image2) { + psImage *sub2 = psImageSubset(image2, region); // Subimage with stamp + stamp->image2 = psKernelAllocFromImage(sub2, size, size); + psFree(sub2); // Drop reference + } + + psImage *wtSub = psImageSubset(weight, region); // Subimage with stamp + stamp->weight = psKernelAllocFromImage(wtSub, size, size); + psFree(wtSub); // Drop reference + + stamp->status = PM_SUBTRACTION_STAMP_CALCULATE; + } + + return true; +} + +#if 0 +bool pmSubtractionStampsGenerate(pmSubtractionStampList *stamps, float fwhm, int kernelSize) +{ + PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(stamps, false); + PS_ASSERT_FLOAT_LARGER_THAN(fwhm, 0.0, false); + PS_ASSERT_INT_NONNEGATIVE(kernelSize, false); + + int size = kernelSize + stamps->footprint; // Size of postage stamps + int num = stamps->num; // Number of stamps + float sigma = fwhm / (2.0 * sqrtf(2.0 * log(2.0))); // Gaussian sigma + + for (int i = 0; i < num; i++) { + pmSubtractionStamp *stamp = stamps->stamps->data[i]; // Stamp of interest + if (!(stamp->status & PM_SUBTRACTION_STAMP_CALCULATE)) { + continue; + } + + float x = stamp->x, y = stamp->y; // Coordinates of stamp + float flux = stamp->flux; // Flux of star + if (!isfinite(flux)) { + psWarning("Unable to generate PSF for stamp %d --- bad flux.", i); + stamp->status = PM_SUBTRACTION_STAMP_REJECTED; + continue; + } + + float xStamp = x - (int)(x + 0.5); // x coordinate of star in stamp frame + float yStamp = y - (int)(y + 0.5); // y coordinate of star in stamp frame + + psFree(stamp->image2); + stamp->image2 = psKernelAlloc(-size, size, -size, size); + psKernel *target = stamp->image2; // Target stamp + + // Put in a Waussian, just for fun! + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + float z = (PS_SQR(u + xStamp) + PS_SQR(v + yStamp)) / (2.0 * PS_SQR(sigma)); + target->kernel[v][u] = flux / sigma * 0.5 * M_2_SQRTPI * M_SQRT1_2 / (1.0 + z + PS_SQR(z)); + } + } + + } + + return true; +} +#endif + +pmSubtractionStampList *pmSubtractionStampsSetFromSources(const psArray *sources, const psImage *subMask, + const psRegion *region, int footprint, + float spacing, pmSubtractionMode mode) +{ + PS_ASSERT_ARRAY_NON_NULL(sources, NULL); + // Let pmSubtractionStampsSet take care of the rest of the assertions + + int numSources = sources->n; // Number of stars + + psVector *x = psVectorAlloc(numSources, PS_TYPE_F32); // x coordinates + psVector *y = psVectorAlloc(numSources, PS_TYPE_F32); // y coordinates + psVector *flux = psVectorAlloc(numSources, PS_TYPE_F32); // Fluxes + + for (int i = 0; i < numSources; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (source->modelPSF) { + x->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + y->data.F32[i] = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + } else { + x->data.F32[i] = source->peak->xf; + y->data.F32[i] = source->peak->yf; + } + flux->data.F32[i] = powf(10.0, -0.4 * source->psfMag); + } + + pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, NULL, subMask, region, + footprint, spacing, mode); // Stamps to return + psFree(x); + psFree(y); + psFree(flux); + + if (!stamps) { + psError(PS_ERR_UNKNOWN, false, "Unable to set stamps from sources."); + } + + return stamps; +} + + +pmSubtractionStampList *pmSubtractionStampsSetFromFile(const char *filename, const psImage *image, + const psImage *subMask, const psRegion *region, + int footprint, float spacing, pmSubtractionMode mode) +{ + PS_ASSERT_STRING_NON_EMPTY(filename, NULL); + // Let pmSubtractionStampsSet take care of the rest of the assertions + + const char *format = (mode == PM_SUBTRACTION_MODE_TARGET ? "%f %f %f" : "%f %f"); // Format of file + psArray *data = psVectorsReadFromFile(filename, format); + if (!data) { + psError(PS_ERR_IO, false, "Unable to read stamps file %s", filename); + return NULL; + } + psVector *x = data->data[0], *y = data->data[1]; // Stamp positions + psVector *flux = (mode == PM_SUBTRACTION_MODE_TARGET ? data->data[2] : NULL); // Stamp fluxes + + // Correct for IRAF/FITS (unit-offset) positions to C (zero-offset) positions + psBinaryOp(x, x, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + psBinaryOp(y, y, "-", psScalarAlloc(1.0, PS_TYPE_F32)); + + pmSubtractionStampList *stamps = pmSubtractionStampsSet(x, y, flux, image, subMask, region, footprint, + spacing, mode); + psFree(data); + + return stamps; + +} Index: /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionStamps.h =================================================================== --- /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionStamps.h (revision 16477) +++ /tags/pap_root_080214/psModules/src/imcombine/pmSubtractionStamps.h (revision 16477) @@ -0,0 +1,126 @@ +#ifndef PM_SUBTRACTION_STAMPS_H +#define PM_SUBTRACTION_STAMPS_H + +#include + +#include "pmSubtractionKernels.h" + +/// Status of stamp +typedef enum { + PM_SUBTRACTION_STAMP_INIT, ///< Initial state + PM_SUBTRACTION_STAMP_FOUND, ///< Found a suitable source for this stamp + PM_SUBTRACTION_STAMP_CALCULATE, ///< Calculate matrix and vector values for this stamp + PM_SUBTRACTION_STAMP_USED, ///< Use this stamp + PM_SUBTRACTION_STAMP_REJECTED, ///< This stamp has been rejected + PM_SUBTRACTION_STAMP_NONE ///< No stamp in this region +} pmSubtractionStampStatus; + +/// A list of stamps +typedef struct { + long num; ///< Number of stamps + psArray *stamps; ///< The stamps + psArray *regions; ///< Regions for each stamp + psArray *x, *y; ///< Coordinates for possible stamps (or NULL) + psArray *flux; ///< Fluxes for possible stamps (or NULL) + int footprint; ///< Half-size of stamps +} pmSubtractionStampList; + +/// Allocate a list of stamps +pmSubtractionStampList *pmSubtractionStampListAlloc(int numCols, // Number of columns in image + int numRows, // Number of rows in image + const psRegion *region, // Region for stamps, or NULL + int footprint, // Half-size of stamps + float spacing // Rough average spacing between stamps + ); + +/// Assertion for stamp list to be valid +#define PM_ASSERT_SUBTRACTION_STAMP_LIST_NON_NULL(LIST, RETURNVALUE) { \ + PS_ASSERT_PTR_NON_NULL(LIST, RETURNVALUE); \ + PS_ASSERT_ARRAY_NON_NULL((LIST)->stamps, RETURNVALUE); \ + PS_ASSERT_ARRAY_NON_NULL((LIST)->regions, RETURNVALUE); \ + PS_ASSERT_INT_POSITIVE((LIST)->num, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((LIST)->stamps, (LIST)->num, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((LIST)->regions, (LIST)->num, RETURNVALUE); \ + PS_ASSERT_INT_NONNEGATIVE((LIST)->footprint, RETURNVALUE); \ + if ((LIST)->x || (LIST)->y || (LIST)->flux) { \ + PS_ASSERT_ARRAY_NON_NULL((LIST)->x, RETURNVALUE); \ + PS_ASSERT_ARRAY_NON_NULL((LIST)->y, RETURNVALUE); \ + PS_ASSERT_ARRAY_NON_NULL((LIST)->flux, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((LIST)->x, (LIST)->num, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((LIST)->y, (LIST)->num, RETURNVALUE); \ + PS_ASSERT_ARRAY_SIZE((LIST)->flux, (LIST)->num, RETURNVALUE); \ + } \ +} + +/// A stamp for image subtraction +typedef struct { + float x, y; ///< Position + float flux; ///< Flux + float xNorm, yNorm; ///< Normalised position + psKernel *image1; ///< Reference image postage stamp + psKernel *image2; ///< Input image postage stamp + psKernel *weight; ///< Weight image postage stamp, or NULL + psArray *convolutions1; ///< Convolutions of image 1 for each kernel component, or NULL + psArray *convolutions2; ///< Convolutions of image 2 for each kernel component, or NULL + psImage *matrix1, *matrix2; ///< Least-squares matrices for each image, or NULL + psImage *matrixX; ///< Cross-matrix (for mode DUAL), or NULL + psVector *vector1, *vector2; ///< Least-squares vectors for each image, or NULL + pmSubtractionStampStatus status; ///< Status of stamp +} pmSubtractionStamp; + +/// Allocate a stamp +pmSubtractionStamp *pmSubtractionStampAlloc(void); + +/// Find stamps on an image +pmSubtractionStampList *pmSubtractionStampsFind(pmSubtractionStampList *stamps, ///< Output stamps, or NULL + const psImage *image, ///< Image for which to find stamps + const psImage *mask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + float threshold, ///< Threshold for stamps in the image + int footprint, ///< Half-size for stamps + float spacing, ///< Rough spacing for stamps + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Set stamps based on a list of x,y +pmSubtractionStampList *pmSubtractionStampsSet(const psVector *x, ///< x coordinates for each stamp + const psVector *y, ///< y coordinates for each stamp + const psVector *flux, ///< Flux for each stamp, or NULL + const psImage *image, ///< Image for flux of stamp + const psImage *mask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + int footprint, ///< Half-size for stamps + float spacing, ///< Rough spacing for stamps + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Set stamps based on a list of sources +pmSubtractionStampList *pmSubtractionStampsSetFromSources( + const psArray *sources, ///< Sources for each stamp + const psImage *subMask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + int footprint, ///< Half-size for stamps + float spacing, ///< Rough spacing for stamps + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Set stamps based on values in a file +pmSubtractionStampList *pmSubtractionStampsSetFromFile( + const char *filename, ///< Filename of file containing x,y (or x,y,flux) on each line + const psImage *image, ///< Image for flux of stamp + const psImage *subMask, ///< Mask, or NULL + const psRegion *region, ///< Region to search, or NULL + int footprint, ///< Half-size for stamps + float spacing, ///< Rough spacing for stamps + pmSubtractionMode mode ///< Mode for subtraction + ); + +/// Extract stamps from the images +bool pmSubtractionStampsExtract(pmSubtractionStampList *stamps, ///< Stamps + psImage *image1, ///< Reference image + psImage *image2, ///< Input image (or NULL) + psImage *weight, ///< Weight (variance) map + int kernelSize ///< Kernel half-size + ); + +#endif