Index: /trunk/psModules/src/imcombine/Makefile.am =================================================================== --- /trunk/psModules/src/imcombine/Makefile.am (revision 13456) +++ /trunk/psModules/src/imcombine/Makefile.am (revision 13457) @@ -6,5 +6,5 @@ libpsmodulesimcombine_la_SOURCES = \ pmReadoutCombine.c \ - pmImageCombine.c \ + pmStack.c \ pmSubtraction.c \ pmSubtractionKernels.c \ @@ -13,5 +13,5 @@ pkginclude_HEADERS = \ pmReadoutCombine.h \ - pmImageCombine.h \ + pmStack.h \ pmSubtraction.h \ pmSubtractionKernels.h \ Index: /trunk/psModules/src/imcombine/pmStack.c =================================================================== --- /trunk/psModules/src/imcombine/pmStack.c (revision 13457) +++ /trunk/psModules/src/imcombine/pmStack.c (revision 13457) @@ -0,0 +1,912 @@ +/** @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.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-05-22 03:59:32 $ + * + * Copyright 2004-2007 Institute for Astronomy, University of Hawaii + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include + +#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 + psStats *stats; // Statistics +} combineBuffer; + +static void combineBufferFree(combineBuffer *buffer) +{ + psFree(buffer->pixels); + psFree(buffer->masks); + psFree(buffer->stats); + return; +} + +static 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->stats = psStatsAlloc(PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_STDEV); + + return buffer; +} + + +// Deallocator for the stack data +static void stackDataFree(pmStackData *data) +{ + psFree(data->detector); + psFree(data->sky); + psFree(data->pixels); + return; +} + + +// 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 + const psArray *inputs, // Stack data + const psVector *weights, // Weights for data, 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 + + if (buffer) { + psMemIncrRefCounter(buffer); + } else { + buffer = combineBufferAlloc(num); + } + + psVector *pixelData = buffer->pixels; // Values for the pixel of interest + psVector *pixelMasks = buffer->masks; // Masks for the pixel of interest + psStats *stats = buffer->stats; // Statistics + + // Extract the pixel and mask data + for (int i = 0; i < num; i++) { + pmStackData *data = inputs->data[i]; // Stack data of interest + psImage *image = data->sky->image; // Image of interest + psImage *mask = data->sky->mask; // Mask of interest + pixelData->data.F32[i] = image->data.F32[y][x]; + pixelMasks->data.PS_TYPE_MASK_DATA[i] = mask->data.PS_TYPE_MASK_DATA[y][x]; + } + + // Record the value derived with no clipping, because pixels rejected using the harsh clipping applied in + // the first pass might later be accepted. + if (!psVectorStats(stats, pixelData, pixelMasks, weights, maskVal)) { + // Can't do anything about an error except give it a NAN and mask + psErrorClear(); + image->data.F32[y][x] = NAN; + mask->data.PS_TYPE_MASK_DATA[y][x] = bad; + return false; + } + image->data.F32[y][x] = stats->sampleMean; + mask->data.PS_TYPE_MASK_DATA[y][x] = 0; + + 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++) { + numClipped = 0; + float limit = rej * stats->sampleStdev; // Rejection limit + for (int j = 0; j < num; j++) { + if (pixelMasks->data.PS_TYPE_MASK_DATA[j] & ignore) { + continue; + } + float diff = fabsf(pixelData->data.F32[j] - stats->sampleMean); // Absolute difference from mean + if (diff > limit) { + pmStackData *data = inputs->data[i]; // Stack data of interest + data->pixels = psPixelsAdd(data->pixels, PIXEL_LIST_BUFFER, x, y); + pixelMasks->data.PS_TYPE_MASK_DATA[j] |= bad; + numClipped++; + } + } + if (!psVectorStats(stats, pixelData, pixelMasks, weights, maskVal)) { + // Can't do anything about an error except give it a NAN and mask + psErrorClear(); + image->data.F32[y][x] = NAN; + mask->data.PS_TYPE_MASK_DATA[y][x] = bad; + return false; + } + } + + return true; +} + + +// Ensure the input array of pmStackData is valid, and get some details out of it +static bool validateInputData(bool *haveDetector, // Do we have the detector images? + bool *haveSky, // Do we have 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 = input->data[0]; // First image off the rank + PS_ASSERT_PTR_NON_NULL(data, false); + assert(psMemGetDeallocator(data) == (psFreeFunc)stackDataFree); // Ensure it's the right type + *haveDetector = (data->detector != NULL); + if (*haveDetector) { + PS_ASSERT_IMAGE_NON_NULL(data->detector->image, false); + PS_ASSERT_IMAGE_NON_NULL(data->detector->mask, false); + } + *haveSky = (data->sky != NULL); + if (*haveSky) { + PS_ASSERT_IMAGE_NON_NULL(data->sky->image, false); + PS_ASSERT_IMAGE_NON_NULL(data->sky->mask, false); + *numCols = data->sky->image->numCols; + *numRows = data->sky->image->numRows; + } + *havePixels = (data->pixels != NULL); + } + for (int i = 1; i < *num; i++) { + pmStackData *data = input->data[i]; // Stack data for this input + assert(psMemGetDeallocator(data) == (psFreeFunc)stackDataFree); // Ensure it's the right type + if ((*haveDetector && !data->detector) || (data->detector && !*haveDetector)) { + psError(PS_ERR_UNEXPECTED_NULL, true, + "The detector readout is specified in some but not all inputs."); + return false; + } + if ((*haveSky && !data->sky) || (data->sky && !*haveSky)) { + psError(PS_ERR_UNEXPECTED_NULL, true, "The sky cell 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; + } + if (*haveDetector) { + PS_ASSERT_IMAGE_NON_NULL(data->detector->image, false); + PS_ASSERT_IMAGE_NON_NULL(data->detector->mask, false); + PS_ASSERT_IMAGE_TYPE(data->detector->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_TYPE(data->detector->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->detector->image, data->detector->mask, false); + } + if (*haveSky) { + PS_ASSERT_IMAGE_NON_NULL(data->sky->image, false); + PS_ASSERT_IMAGE_NON_NULL(data->sky->mask, false); + PS_ASSERT_IMAGE_TYPE(data->sky->image, PS_TYPE_F32, false); + PS_ASSERT_IMAGE_TYPE(data->sky->mask, PS_TYPE_MASK, false); + PS_ASSERT_IMAGE_SIZE(data->sky->image, *numCols, *numRows, false); + PS_ASSERT_IMAGES_SIZE_EQUAL(data->sky->image, data->sky->mask, false); + } + } + + return true; +} + + + +// 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->sky->image; // Image to convolve + psImage *mask = data->sky->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, valuesMask, NULL, 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(values); + psFree(valuesMask); + + return true; +} + + + +// 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++) { + psArray *colMap = psArrayAlloc(numCols); // The pixel map columns + for (int x = 0; x < numCols; x++) { + colMap->data[x] = NULL; + } + map->data[y] = psArrayAlloc(numCols); + } + + for (int i = 0; i < input->n; i++) { + pmStackData *data = input->data[i]; + assert(data && 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 + 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 *sky, float seeing, float weight) +{ + pmStackData *data = psAlloc(sizeof(pmStackData)); // Stack data, to return + psMemSetDeallocator(data, (psFreeFunc)stackDataFree); + + data->detector = NULL; + + data->sky = sky; + data->pixels = NULL; + data->seeing = seeing; + 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 haveDetector; // Do we have the detector images? + bool haveSky; // Do we have the sky images? + bool havePixels; // Do we have lists of pixels? + int num; // Number of inputs + int numCols, numRows; // Size of (sky) images + if (!validateInputData(&haveDetector, &haveSky, &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 (!haveDetector && !haveSky) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Nothing to combine!"); + return false; + } + + if (!haveSky) { + // Need to generate the sky cell images + + // ... + +#if 0 + numCols = data->sky->image->numCols; + numRows = data->sky->image->numRows; +#endif + } + + // Pull out the weights + psVector *weights = psVectorAlloc(num, PS_TYPE_F32); + for (int i = 0; i < num; i++) { + pmStackData *data = input->data[i]; // Stack data for this input + weights->data.F32[i] = data->weight; + } + + combineBuffer *buffer = combineBufferAlloc(num); // Buffer for combination + + if (havePixels) { + // Only combine select pixels + + } 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; + } + + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + combinePixels(combinedImage, combinedMask, input, weights, x, y, + maskVal, bad, numIter, rej, buffer); + } + } + + } + + psFree(buffer); + + return true; +} + +bool pmStackReject(psArray *input, psMaskType maskVal, float extent, float threshold) +{ + bool haveDetector; // Do we have the detector images? + bool haveSky; // Do we have the sky images? + bool havePixels; // Do we have lists of pixels? + int num; // Number of inputs + int numCols, numRows; // Size of (sky) images + if (!validateInputData(&haveDetector, &haveSky, &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, seeing, "-", psScalarAlloc(PS_SQR(seeingMax), PS_TYPE_F32)); + (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 + + 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); + } + + } + + return true; +} + + + + + + + + + +#if 0 +/****************************************************************************** +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); +} +#endif Index: /trunk/psModules/src/imcombine/pmStack.h =================================================================== --- /trunk/psModules/src/imcombine/pmStack.h (revision 13457) +++ /trunk/psModules/src/imcombine/pmStack.h (revision 13457) @@ -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.1 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-05-22 03:59:32 $ + * 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 *detector; ///< Original (unwarped) readout from the detector + pmReadout *sky; ///< Warped readout (sky cell) + psPixels *pixels; ///< Pixels to inspect or reject + float seeing; ///< Seeing FWHM (pixels) + float weight; ///< Weight to apply +} pmStackData; + +/// Constructor +pmStackData *pmStackDataAlloc(pmReadout *sky, ///< Warped readout (sky cell) + float seeing, ///< Seeing FWHM (pixels) + 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 + +/// 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 + ); + +#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: /trunk/psModules/src/psmodules.h =================================================================== --- /trunk/psModules/src/psmodules.h (revision 13456) +++ /trunk/psModules/src/psmodules.h (revision 13457) @@ -67,5 +67,5 @@ // the following headers are from psModule:imcombine -#include +#include #include #include