Index: /trunk/psModules/src/imcombine/pmSubtractionSimple.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionSimple.c (revision 35716) +++ /trunk/psModules/src/imcombine/pmSubtractionSimple.c (revision 35716) @@ -0,0 +1,239 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmHDUUtils.h" +#include "pmSubtractionTypes.h" +#include "pmSubtraction.h" +#include "pmSubtractionParams.h" +#include "pmSubtractionKernels.h" +#include "pmSubtractionStamps.h" +#include "pmSubtractionEquation.h" +#include "pmSubtractionAnalysis.h" +#include "pmSubtractionMask.h" +#include "pmSubtractionThreads.h" +#include "pmSubtractionVisual.h" +#include "pmErrorCodes.h" + +#include "pmSubtractionSimple.h" + +bool pmSubtractionSimpleMatch(pmReadout *conv1, + pmReadout *conv2, + const pmReadout *ro1, + const pmReadout *ro2, + const psArray *sources, + int size + ) { + // + // We've already validated the input values at this level + float sig2fwhm = 2.0 * sqrt(2.0 * log(2.0)); + float fwhm1,fwhm2; + float sigma1,sigma2,sigmaKern; + int convolution_direction = 0; + + // Allocate images, as this is usually done by subtractionMatchAlloc after this function is called. + psImageMaskType maskBad = 1; + int numCols = ro1->image->numCols; + int numRows = ro1->image->numRows; + + if (!conv1->image) { + conv1->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv1->image, NAN); + if (ro1->variance) { + if (!conv1->variance) { + conv1->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv1->variance, NAN); + } + if (!conv1->mask) { + conv1->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + } + psImageInit(conv1->mask, maskBad); + if (!conv2->image) { + conv2->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv2->image, NAN); + if (ro2->variance) { + if (!conv2->variance) { + conv2->variance = psImageAlloc(numCols, numRows, PS_TYPE_F32); + } + psImageInit(conv2->variance, NAN); + } + if (!conv2->mask) { + conv2->mask = psImageAlloc(numCols, numRows, PS_TYPE_IMAGE_MASK); + } + psImageInit(conv2->mask, maskBad); + + + + + psImage *image1 = ro1->image; + psImage *mask1 = ro1->mask; + psImage *var1 = ro1->variance; + + psImage *image2 = ro2->image; + psImage *mask2 = ro2->mask; + psImage *var2 = ro2->variance; + + psImage *imageC1 = conv1->image; + psImage *maskC1 = conv1->mask; + psImage *varC1 = conv1->variance; + + psImage *imageC2 = conv2->image; + psImage *maskC2 = conv2->mask; + psImage *varC2 = conv2->variance; + + // + // Determine Gaussian widths + pmSubtractionGetFWHMs(&fwhm1,&fwhm2); + sigma1 = fwhm1 / sig2fwhm; + sigma2 = fwhm2 / sig2fwhm; + if (sigma1 > sigma2) { + convolution_direction = 2; + sigmaKern = sqrt(PS_SQR(sigma1) - PS_SQR(sigma2)); + } + else if (sigma1 < sigma2) { + convolution_direction = 1; + sigmaKern = sqrt(PS_SQR(sigma2) - PS_SQR(sigma1)); + } + + // + // Determine Normalization scaling + psVector *kernelVec = pmSubtractionKernelSIMPLE(sigmaKern,0,size); // This is normalized to unity. + + float normalization = 1.0; + + for (int i = 0; i < 2 * size + 1; i++) { + kernelVec->data.F32[i] *= sqrt(normalization); // Needs to be this because we apply this vector twice. + } + + // + // Do convolutions + for (int y = 0; y < image1->numRows; y++) { + for (int x = 0; x < image1->numCols; x++) { + if (convolution_direction == 1) { + imageC1->data.F32[y][x] = 0.0; + maskC1->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; + varC1->data.F32[y][x] = 0.0; + + imageC2->data.F32[y][x] = image2->data.F32[y][x]; + maskC2->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = mask2->data.PS_TYPE_IMAGE_MASK_DATA[y][x]; + varC2->data.F32[y][x] = var2->data.F32[y][x]; + } + else if (convolution_direction == 2) { + imageC2->data.F32[y][x] = 0.0; + maskC2->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = 0; + varC2->data.F32[y][x] = 0.0; + + imageC1->data.F32[y][x] = image1->data.F32[y][x]; + maskC1->data.PS_TYPE_IMAGE_MASK_DATA[y][x] = mask1->data.PS_TYPE_IMAGE_MASK_DATA[y][x]; + varC1->data.F32[y][x] = var1->data.F32[y][x]; + } + for (int v = -size; v <= size; v++) { + for (int u = -size; u <= size; u++) { + if ((y-v >= 0)&&(x-u >= 0)&& + (y-v < numRows)&&(x-u < numCols)) { + if (convolution_direction == 1) { + imageC1->data.F32[y][x] += kernelVec->data.F32[v + size] * kernelVec->data.F32[u + size] * + image1->data.F32[y - v][x - u]; + } + else if (convolution_direction == 2) { + imageC2->data.F32[y][x] += kernelVec->data.F32[v + size] * kernelVec->data.F32[u + size] * + image2->data.F32[y - v][x - u]; + } + } + } + } + } + } + psFree(kernelVec); + + // Make a fake pmSubtractionKernels element so we can add it appropriately. + // I call it fake because we've successfully done everything at this point + // without having to define these things. + psVector *fwhms = psVectorAlloc(1,PS_TYPE_F32); + fwhms->data.F32[0] = sigmaKern * sig2fwhm; + psVector *orders = psVectorAlloc(1,PS_TYPE_S32); + orders->data.S32[0] = 0; + psRegion bounds; + bounds.x0 = 0; + bounds.y0 = 0; + bounds.x1 = numCols; + bounds.y1 = numRows; + pmSubtractionKernels *kernels = pmSubtractionKernelsAlloc(1,PM_SUBTRACTION_KERNEL_SIMPLE, + size,fwhms,orders,0,0.0,bounds, + convolution_direction); + pmSubtractionKernelsMakeDescription(kernels); // Need this defined + pmSubtractionKernelPreCalc *preCalc = pmSubtractionKernelPreCalcAlloc(PM_SUBTRACTION_KERNEL_SIMPLE,0,0,size,sigmaKern); + kernels->widths->data.F32[0] = sigmaKern * sig2fwhm; + kernels->u->data.S32[0] = 0; + kernels->v->data.S32[0] = 0; + if (kernels->preCalc->data[0]) { + psFree(kernels->preCalc->data[0]); + } + kernels->preCalc->data[0] = preCalc; + kernels->solution1 = psVectorAlloc(3,PS_TYPE_F64); + kernels->solution1->data.F32[0] = 1.0; + kernels->solution1->data.F32[1] = 0.0; + kernels->solution1->data.F32[2] = 0.0; +/* kernels->solution2 = psVectorAlloc(3,PS_TYPE_F64); */ +/* kernels->solution2->data.F32[0] = 1.0; */ + kernels->solution1err = psVectorAlloc(3,PS_TYPE_F32); + kernels->solution1err->data.F32[0] = 0.0; + kernels->solution1err->data.F32[1] = 0.0; + kernels->solution1err->data.F32[2] = 0.0; +/* kernels->solution2err = psVectorAlloc(3,PS_TYPE_F32); */ +/* kernels->solution2err->data.F32[0] = 0.0; */ +/* kernels->solution2err->data.F32[1] = 0.0; */ +/* kernels->solution2err->data.F32[2] = 0.0; */ + + + // Actually add it to the headers + + psMetadata *outAnalysis = psMetadataAlloc(); + psMetadata *outHeader = psMetadataAlloc(); + if (!pmSubtractionAnalysis(outAnalysis,outHeader,kernels,NULL,numCols,numRows)) { + psError(psErrorCodeLast(),false,"Unable to generate QA data"); + psFree(fwhms); + psFree(orders); + psFree(preCalc); + psFree(kernels); + } + psFree(fwhms); + psFree(orders); + psFree(preCalc); + psFree(kernels); + + if (conv1) { + conv1->analysis = psMetadataCopy(conv1->analysis, outAnalysis); + } + if (conv2) { + conv2->analysis = psMetadataCopy(conv2->analysis, outAnalysis); + } + + if (conv1 && conv1->parent) { + pmHDU *hdu = pmHDUFromCell(conv1->parent); + if (hdu) { + hdu->header = psMetadataCopy(hdu->header, outHeader); + } + } + if (conv2 && conv2->parent) { + pmHDU *hdu = pmHDUFromCell(conv2->parent); + if (hdu) { + hdu->header = psMetadataCopy(hdu->header, outHeader); + } + } + psFree(outAnalysis); + psFree(outHeader); + // Return + + return(true); +} Index: /trunk/psModules/src/imcombine/pmSubtractionSimple.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionSimple.h (revision 35716) +++ /trunk/psModules/src/imcombine/pmSubtractionSimple.h (revision 35716) @@ -0,0 +1,19 @@ +#ifndef PM_SUBTRACTION_SIMPLE_H +#define PM_SUBTRACTION_SIMPLE_H + +#include + +#include +#include +#include + +bool pmSubtractionSimpleMatch(pmReadout *conv1, + pmReadout *conv2, + const pmReadout *ro1, + const pmReadout *ro2, + const psArray *sources, + int size + ); + + +#endif