Index: /trunk/ippconfig/recipes/filerules-mef.mdc =================================================================== --- /trunk/ippconfig/recipes/filerules-mef.mdc (revision 27364) +++ /trunk/ippconfig/recipes/filerules-mef.mdc (revision 27365) @@ -111,4 +111,5 @@ PPSUB.REF.VARIANCE INPUT @FILES FPA VARIANCE PPSUB.REF.SOURCES INPUT @FILES FPA CMF +PPSUB.INPUT.KERNEL INPUT @FILES FPA SUBKERNEL ## files used by ppstack Index: /trunk/ippconfig/recipes/filerules-simple.mdc =================================================================== --- /trunk/ippconfig/recipes/filerules-simple.mdc (revision 27364) +++ /trunk/ippconfig/recipes/filerules-simple.mdc (revision 27365) @@ -74,4 +74,5 @@ PPSUB.REF.VARIANCE INPUT @FILES FPA VARIANCE PPSUB.REF.SOURCES INPUT @FILES FPA CMF +PPSUB.INPUT.KERNEL INPUT @FILES FPA SUBKERNEL ## files used by ppstack Index: /trunk/ippconfig/recipes/filerules-split.mdc =================================================================== --- /trunk/ippconfig/recipes/filerules-split.mdc (revision 27364) +++ /trunk/ippconfig/recipes/filerules-split.mdc (revision 27365) @@ -84,4 +84,5 @@ PPSUB.REF.VARIANCE INPUT @FILES FPA VARIANCE PPSUB.REF.SOURCES INPUT @FILES FPA CMF +PPSUB.INPUT.KERNEL INPUT @FILES FPA SUBKERNEL ## files used by ppstack Index: /trunk/ppStack/src/ppStackMatch.c =================================================================== --- /trunk/ppStack/src/ppStackMatch.c (revision 27364) +++ /trunk/ppStack/src/ppStackMatch.c (revision 27365) @@ -407,5 +407,5 @@ if (threads > 0) { - pmSubtractionThreadsInit(readout, fake); + pmSubtractionThreadsInit(); } @@ -423,5 +423,5 @@ psFree(conv); if (threads > 0) { - pmSubtractionThreadsFinalize(readout, fake); + pmSubtractionThreadsFinalize(); } return false; @@ -440,5 +440,5 @@ psFree(widthsCopy); if (threads > 0) { - pmSubtractionThreadsFinalize(readout, fake); + pmSubtractionThreadsFinalize(); } return false; @@ -458,5 +458,5 @@ psFree(widthsCopy); if (threads > 0) { - pmSubtractionThreadsFinalize(readout, fake); + pmSubtractionThreadsFinalize(); } return false; @@ -495,5 +495,5 @@ if (threads > 0) { - pmSubtractionThreadsFinalize(readout, fake); + pmSubtractionThreadsFinalize(); } Index: /trunk/ppSub/src/Makefile.am =================================================================== --- /trunk/ppSub/src/Makefile.am (revision 27364) +++ /trunk/ppSub/src/Makefile.am (revision 27365) @@ -1,3 +1,3 @@ -bin_PROGRAMS = ppSub ppSubKernel +bin_PROGRAMS = ppSub ppSubKernel ppSubConvolve if HAVE_SVNVERSION @@ -55,4 +55,12 @@ ppSubKernel.c +ppSubConvolve_CPPFLAGS = $(PSLIB_CFLAGS) $(PSMODULE_CFLAGS) $(PPSTATS_CFLAGS) $(PSPHOT_CFLAGS) $(PPSUB_CFLAGS) +ppSubConvolve_LDFLAGS = $(PSLIB_LIBS) $(PSMODULE_LIBS) $(PPSTATS_LIBS) $(PSPHOT_LIBS) $(PPSUB_LIBS) + +ppSubConvolve_SOURCES = \ + ppSubConvolve.c \ + ppSubExit.c \ + ppSubVersion.c + noinst_HEADERS = \ ppSub.h Index: /trunk/ppSub/src/ppSubConvolve.c =================================================================== --- /trunk/ppSub/src/ppSubConvolve.c (revision 27365) +++ /trunk/ppSub/src/ppSubConvolve.c (revision 27365) @@ -0,0 +1,281 @@ +#include +#include +#include + +#include "ppSub.h" + +// Print usage information and die +static void usage(const char *program, // Name of the program + psMetadata *arguments, // Command-line arguments + pmConfig *config // Configuration + ) +{ + fprintf(stderr, "\nPan-STARRS image convolution\n\n"); + fprintf(stderr, "Usage: %s OUTPUT_ROOT\n\n", + program); + fprintf(stderr, "\n"); + psArgumentHelp(arguments); + psFree(config); + pmConfigDone(); + psLibFinalize(); + exit(PS_EXIT_CONFIG_ERROR); +} + +// Put filename in it's own list +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[]) +{ + pmFPAview *view = pmFPAviewAlloc(0); // Pointer into FPA hierarchy + + pmConfig *config = pmConfigRead(&argc, argv, PPSUB_RECIPE); // Configuration + if (!config) { + psError(psErrorCodeLast(), false, "Unable to read configuration"); + goto die; + } + + bool reference = false; // Input is actually the reference image? + int threads = 0; // Number of threads + { + psMetadata *arguments = config->arguments; // Command-line arguments + psMetadataAddStr(arguments, PS_LIST_TAIL, "-image", 0, "Input image", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-mask", 0, "Input mask", NULL); + psMetadataAddStr(arguments, PS_LIST_TAIL, "-kernel", 0, "Convolution kernel", NULL); + psMetadataAddBool(arguments, PS_LIST_TAIL, "-reference", 0, "Input is actually reference?", false); + psMetadataAddS32(arguments, PS_LIST_TAIL, "-threads", 0, "Threads to use", 0); + + if (argc == 1 || !psArgumentParse(arguments, &argc, argv) || argc != 2) { + usage(argv[0], arguments, config); + } + + const char *inImage = psMetadataLookupStr(NULL, arguments, "-image"); // Input image + const char *inMask = psMetadataLookupStr(NULL, arguments, "-mask"); // Input mask + const char *inKernel = psMetadataLookupStr(NULL, arguments, "-kernel"); // Input kernel + if (!inImage || !inKernel) { + usage(argv[0], arguments, config); + } + + reference = psMetadataLookupBool(NULL, arguments, "-reference"); + threads = psMetadataLookupS32(NULL, arguments, "-threads"); + fileList("PPSUB.INPUT", inImage, "Input image", config); + fileList("PPSUB.INPUT.MASK", inMask, "Input mask", config); + fileList("PPSUB.INPUT.KERNEL", inKernel, "Input kernel", config); + psMetadataAddStr(config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "Name of the output image", argv[1]); + + bool status; // Status of file definition + pmFPAfile *image = pmFPAfileDefineFromArgs(&status, config, "PPSUB.INPUT", "PPSUB.INPUT"); + if (!image || !status || image->type != PM_FPA_FILE_IMAGE) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define input image file."); + goto die; + } + pmFPAfile *mask = pmFPAfileBindFromArgs(&status, image, config, "PPSUB.INPUT.MASK", + "PPSUB.INPUT.MASK"); + if (!mask || !status || mask->type != PM_FPA_FILE_MASK) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define input mask file."); + goto die; + } + + pmFPAfile *kernel = pmFPAfileBindFromArgs(&status, image, config, "PPSUB.INPUT.KERNEL", + "PPSUB.INPUT.KERNEL"); + if (!kernel || !status || kernel->type != PM_FPA_FILE_SUBKERNEL) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define kernel file."); + goto die; + } + + pmFPAfile *output = pmFPAfileDefineFromFile(config, image, 1, 1, "PPSUB.INPUT.CONV"); + if (!output || output->type != PM_FPA_FILE_IMAGE) { + psError(PPSUB_ERR_CONFIG, false, "Unable to define output file."); + goto die; + } + output->save = true; + } + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSUB_RECIPE); // Recipe for ppSim + if (!recipe) { + psError(PPSUB_ERR_CONFIG, false, "Unable to find recipe %s", PPSUB_RECIPE); + goto die; + } + + // Read everything + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + view->chip = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + view->cell = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_BEFORE)) { + psError(psErrorCodeLast(), false, "File checks failed."); + goto die; + } + + + pmReadout *inRO = pmFPAfileThisReadout(config->files, view, "PPSUB.INPUT"); // Input readout + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, inRO->analysis, + PM_SUBTRACTION_ANALYSIS_KERNEL); // Convolution kernel + pmCell *outCell = pmFPAfileThisCell(config->files, view, "PPSUB.INPUT.CONV"); // Output cell + pmReadout *outRO = pmReadoutAlloc(outCell); // Output readout + pmConceptsCopyFPA(outRO->parent->parent->parent, inRO->parent->parent->parent, true, true); + pmHDU *inHDU = pmHDUFromCell(inRO->parent); // Input header + pmHDU *outHDU = pmHDUFromCell(outCell); // Output header + outHDU->header = psMetadataCopy(outHDU->header, inHDU->header); + ppSubVersionHeader(outHDU->header); + + pmReadout *input, *ref, *inConv, *refConv; // Input and output readouts + if (reference) { + ref = inRO; + refConv = outRO; + input = NULL; + inConv = NULL; + } else { + input = inRO; + inConv = outRO; + ref = NULL; + refConv = NULL; + } + + bool convolve = false; // Do we need to convolve? + float norm = 1.0; // Normalisation to apply + switch (kernels->mode) { + case PM_SUBTRACTION_MODE_1: + if (!reference) { + convolve = true; + } + norm = 1.0 / kernels->solution1->data.F64[PM_SUBTRACTION_INDEX_NORM(kernels)]; + break; + case PM_SUBTRACTION_MODE_2: + if (reference) { + convolve = true; + } + // No normalisation; should already be normalised appropriately + break; + case PM_SUBTRACTION_MODE_DUAL: { + // We only have a single input, so hack together a single-mode kernel + int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index of normalisation term + norm = 1.0 / kernels->solution1->data.F64[normIndex]; + convolve = true; + if (reference) { + for (int i = 0; i < kernels->solution2->n; i++) { + kernels->solution1->data.F64[i] = -kernels->solution2->data.F64[i]; + } + kernels->solution1->data.F64[normIndex] = 1.0; + kernels->solution1->data.F64[PM_SUBTRACTION_INDEX_BG(kernels)] = 0.0; + kernels->mode = PM_SUBTRACTION_MODE_2; + } else { + kernels->mode = PM_SUBTRACTION_MODE_1; + } + psFree(kernels->solution2); + kernels->solution2 = NULL; + break; + } + default: + psAbort("Unrecognised kernel mode: %x", kernels->mode); + } + + if (convolve) { + int stride = psMetadataLookupS32(NULL, recipe, "STRIDE"); // Size of convolution patches + float kernelErr = psMetadataLookupF32(NULL, recipe, "KERNEL.ERR"); // Relative sys error in kernel + float covarFrac = psMetadataLookupF32(NULL, recipe, "COVAR.FRAC"); // Fraction for covar calculation + float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction + float poorFrac = psMetadataLookupF32(NULL, recipe, "POOR.FRACTION"); // Fraction for "poor" + + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs + psImageMaskType maskPoor = pmConfigMaskGet("CONV.POOR", config); // Bits to mask for poor pixels + psImageMaskType maskBad = pmConfigMaskGet("BLANK", config); // Bits to mask for bad pixels + + if (threads > 0) { + pmSubtractionThreadsInit(); + } + if (!pmSubtractionMatchPrecalc(inConv, refConv, input, ref, inRO->analysis, + stride, kernelErr, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac)) { + psError(psErrorCodeLast(), false, "Unable to convolve images."); + if (threads > 0) { + pmSubtractionThreadsFinalize(); + } + goto die; + } + if (threads > 0) { + pmSubtractionThreadsFinalize(); + } + } else { + // Subtract background, since that's what ppSub will do with it + psImageMaskType maskVal = pmConfigMaskGet("MASK.VALUE", config); // Bits to mask in inputs + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); + psStatsInit(bg); + if (!psImageBackground(bg, NULL, inRO->image, inRO->mask, maskVal, rng)) { + psError(PPSUB_ERR_DATA, false, "Unable to measure background statistics."); + psFree(bg); + psFree(rng); + goto die; + } + psFree(rng); + outRO->image = (psImage*)psBinaryOp(outRO->image, inRO->image, "-", + psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32)); + psFree(bg); + outRO->data_exists = outRO->parent->data_exists = outRO->parent->parent->data_exists = true; + } + + if (norm != 1.0) { + psBinaryOp(outRO->image, outRO->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); + } + + psFree(outRO); + + die: + { + psExit exitValue = ppSubExitCode(PS_EXIT_SUCCESS); // Exit value + + // Write everything + view->chip = view->cell = view->readout = 0; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + } + view->readout = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + } + view->cell = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + } + view->chip = -1; + if (!pmFPAfileIOChecks(config, view, PM_FPA_AFTER)) { + psError(psErrorCodeLast(), false, "File checks failed."); + exitValue = ppSubExitCode(exitValue); + } + psFree(view); + + psFree(config); + pmConceptsDone(); + pmConfigDone(); + psLibFinalize(); + exit(exitValue); + } + + psAbort("Should never reach here."); +} Index: /trunk/ppSub/src/ppSubMatchPSFs.c =================================================================== --- /trunk/ppSub/src/ppSubMatchPSFs.c (revision 27364) +++ /trunk/ppSub/src/ppSubMatchPSFs.c (revision 27365) @@ -331,5 +331,5 @@ int threads = psMetadataLookupS32(NULL, config->arguments, "-threads"); // Number of threads if (threads > 0) { - pmSubtractionThreadsInit(inRO, refRO); + pmSubtractionThreadsInit(); } @@ -380,5 +380,5 @@ psFree(optWidths); - pmSubtractionThreadsFinalize(inRO, refRO); + pmSubtractionThreadsFinalize(); if (!success) { Index: /trunk/psModules/src/imcombine/pmStackReject.c =================================================================== --- /trunk/psModules/src/imcombine/pmStackReject.c (revision 27364) +++ /trunk/psModules/src/imcombine/pmStackReject.c (revision 27365) @@ -103,5 +103,5 @@ if (!pmSubtractionThreaded()) { - pmSubtractionThreadsInit(NULL, NULL); + pmSubtractionThreadsInit(); } Index: /trunk/psModules/src/imcombine/pmSubtractionMask.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionMask.c (revision 27364) +++ /trunk/psModules/src/imcombine/pmSubtractionMask.c (revision 27365) @@ -42,11 +42,27 @@ pmSubtractionMode mode) { - PM_ASSERT_READOUT_NON_NULL(ro1, NULL); - PM_ASSERT_READOUT_IMAGE(ro1, NULL); - PM_ASSERT_READOUT_MASK(ro1, NULL); - PM_ASSERT_READOUT_NON_NULL(ro2, NULL); - PM_ASSERT_READOUT_IMAGE(ro2, NULL); - PM_ASSERT_READOUT_MASK(ro2, NULL); - PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, NULL); + int numCols = 0, numRows = 0; // Size of the images + if (ro1) { + PM_ASSERT_READOUT_NON_NULL(ro1, NULL); + PM_ASSERT_READOUT_IMAGE(ro1, NULL); + PM_ASSERT_READOUT_MASK(ro1, NULL); + numCols = ro1->image->numCols; + numRows = ro1->image->numRows; + } + if (ro2) { + PM_ASSERT_READOUT_NON_NULL(ro2, NULL); + PM_ASSERT_READOUT_IMAGE(ro2, NULL); + PM_ASSERT_READOUT_MASK(ro2, NULL); + numCols = ro2->image->numCols; + numRows = ro2->image->numRows; + } + if (ro1 && ro2) { + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, NULL); + } + if (!ro1 && !ro2) { + psError(PS_ERR_UNEXPECTED_NULL, true, "No image provided."); + return false; + } + psAssert(numCols > 0 && numRows > 0, "There should be an image provided"); PS_ASSERT_INT_NONNEGATIVE(size, NULL); PS_ASSERT_INT_NONNEGATIVE(footprint, NULL); @@ -56,10 +72,9 @@ } - int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Size of the images - // Dereference inputs for convenience - psF32 **imageData1 = ro1->image->data.F32, **imageData2 = ro2->image->data.F32; - psImageMaskType **maskData1 = ro1->mask->data.PS_TYPE_IMAGE_MASK_DATA, - **maskData2 = ro2->mask->data.PS_TYPE_IMAGE_MASK_DATA; + psF32 **imageData1 = ro1 ? ro1->image->data.F32 : NULL; + psF32 **imageData2 = ro2 ? ro2->image->data.F32 : NULL; + psImageMaskType **maskData1 = ro1 ? ro1->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL; + psImageMaskType **maskData2 = ro2 ? ro2->mask->data.PS_TYPE_IMAGE_MASK_DATA : NULL; // First, a pass through to determine the fraction of bad pixels @@ -69,9 +84,9 @@ for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - if ((maskData1[y][x] & maskVal) || !isfinite(imageData1[y][x])) { + if (ro1 && ((maskData1[y][x] & maskVal) || !isfinite(imageData1[y][x]))) { numBad++; continue; } - if ((maskData2[y][x] & maskVal) || !isfinite(imageData2[y][x])) { + if (ro2 && ((maskData2[y][x] & maskVal) || !isfinite(imageData2[y][x]))) { numBad++; continue; @@ -128,8 +143,8 @@ for (int y = 0; y < numRows; y++) { for (int x = 0; x < numCols; x++) { - if (maskData1[y][x] & maskVal) { + if (ro1 && maskData1[y][x] & maskVal) { maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_1; } - if (maskData2[y][x] & maskVal) { + if (ro2 && maskData2[y][x] & maskVal) { maskData[y][x] |= PM_SUBTRACTION_MASK_BAD_2; } @@ -144,12 +159,12 @@ // Pixels that will be bad (or poor) if we convolve with a bad reference pixel - if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_1, PM_SUBTRACTION_MASK_CONVOLVE_1, - -size, size, -size, size)) { + if (ro1 && !psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_1, PM_SUBTRACTION_MASK_CONVOLVE_1, + -size, size, -size, size)) { psError(psErrorCodeLast(), false, "Unable to convolve bad pixels from mask 1."); psFree(mask); return NULL; } - if (!psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_2, PM_SUBTRACTION_MASK_CONVOLVE_2, - -size, size, -size, size)) { + if (ro2 && !psImageConvolveMask(mask, mask, PM_SUBTRACTION_MASK_BAD_2, PM_SUBTRACTION_MASK_CONVOLVE_2, + -size, size, -size, size)) { psError(psErrorCodeLast(), false, "Unable to convolve bad pixels from mask 2."); psFree(mask); @@ -174,6 +189,6 @@ psAbort("Unsupported subtraction mode: %x", mode); } - if (!psImageConvolveMask(mask, mask, maskRej, PM_SUBTRACTION_MASK_REJ, - -footprint, footprint, -footprint, footprint)) { + if (ro1 && ro2 && !psImageConvolveMask(mask, mask, maskRej, PM_SUBTRACTION_MASK_REJ, + -footprint, footprint, -footprint, footprint)) { psError(psErrorCodeLast(), false, "Unable to convolve bad pixels."); psFree(mask); Index: /trunk/psModules/src/imcombine/pmSubtractionMatch.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 27364) +++ /trunk/psModules/src/imcombine/pmSubtractionMatch.c (revision 27365) @@ -132,4 +132,6 @@ if (subMode != PM_SUBTRACTION_MODE_2) { PM_ASSERT_READOUT_NON_NULL(conv1, false); + PM_ASSERT_READOUT_NON_NULL(ro1, false); + PM_ASSERT_READOUT_IMAGE(ro1, false); if (conv1->image) { psFree(conv1->image); @@ -147,4 +149,6 @@ if (subMode != PM_SUBTRACTION_MODE_1) { PM_ASSERT_READOUT_NON_NULL(conv2, false); + PM_ASSERT_READOUT_NON_NULL(ro2, false); + PM_ASSERT_READOUT_IMAGE(ro2, false); if (conv2->image) { psFree(conv2->image); @@ -161,9 +165,7 @@ } - PM_ASSERT_READOUT_NON_NULL(ro1, false); - PM_ASSERT_READOUT_NON_NULL(ro2, false); - PM_ASSERT_READOUT_IMAGE(ro1, false); - PM_ASSERT_READOUT_IMAGE(ro2, false); - PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false); + if (ro1 && ro2) { + PS_ASSERT_IMAGES_SIZE_EQUAL(ro1->image, ro2->image, false); + } PS_ASSERT_INT_NONNEGATIVE(stride, false); if (isfinite(normFrac)) { @@ -379,6 +381,45 @@ } + int numCols, numRows; // Size of image + if (ro1) { + numCols = ro1->image->numCols; + numRows = ro1->image->numRows; + } else if (ro2) { + numCols = ro2->image->numCols; + numRows = ro2->image->numRows; + } else { + psAbort("No input image provided."); + } + pmSubtractionMaskInvalid(ro1, maskVal); pmSubtractionMaskInvalid(ro2, maskVal); + + // General background subtraction, since this is done in pmSubtractionMatch + { + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for background + if (ro1) { + psStatsInit(bg); + if (!psImageBackground(bg, NULL, ro1->image, ro1->mask, maskVal, rng)) { + psError(PM_ERR_DATA, false, "Unable to measure background statistics."); + psFree(bg); + psFree(rng); + return false; + } + psBinaryOp(ro1->image, ro1->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32)); + } + if (ro2) { + psStatsInit(bg); + if (!psImageBackground(bg, NULL, ro2->image, ro2->mask, maskVal, rng)) { + psError(PM_ERR_DATA, false, "Unable to measure background statistics."); + psFree(bg); + psFree(rng); + return false; + } + psBinaryOp(ro2->image, ro2->image, "-", psScalarAlloc((float)bg->robustMedian, PS_TYPE_F32)); + } + psFree(bg); + psFree(rng); + } psRegion bounds = psRegionSet(NAN, NAN, NAN, NAN); // Bounds of valid pixels @@ -396,5 +437,4 @@ psMetadata *outHeader = psMetadataAlloc(); // Output header values - int numCols = ro1->image->numCols, numRows = ro1->image->numRows; // Size of image subtractionMatchAlloc(conv1, conv2, ro1, ro2, subMask, maskBad, mode, numCols, numRows); @@ -404,6 +444,5 @@ psRegion *region = regions->data[i]; // Region of interest - if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region, - ro1->image->numCols, ro1->image->numRows)) { + if (!pmSubtractionAnalysis(outAnalysis, outHeader, kernel, region, numCols, numRows)) { psError(psErrorCodeLast(), false, "Unable to generate QA data"); psFree(outAnalysis); @@ -454,4 +493,8 @@ return false; } + + // We need both inputs + PM_ASSERT_READOUT_NON_NULL(ro1, false); + PM_ASSERT_READOUT_NON_NULL(ro2, false); PS_ASSERT_INT_NONNEGATIVE(footprint, false); Index: /trunk/psModules/src/imcombine/pmSubtractionThreads.c =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 27364) +++ /trunk/psModules/src/imcombine/pmSubtractionThreads.c (revision 27365) @@ -17,38 +17,5 @@ } -// Initialise a mutex in each of the input -static void subtractionMutexInit(pmReadout *ro) -{ - if (!ro) { - return; - } - - // XXX if (ro->image) { - // XXX psMutexInit(ro->image); - // XXX } - // XXX if (ro->variance) { - // XXX psMutexInit(ro->variance); - // XXX } - - return; -} - -static void subtractionMutexDestroy(pmReadout *ro) -{ - if (!ro) { - return; - } - - // XXX if (ro->image) { - // XXX psMutexDestroy(ro->image); - // XXX } - // XXX if (ro->variance) { - // XXX psMutexDestroy(ro->variance); - // XXX } - - return; -} - -void pmSubtractionThreadsInit(pmReadout *in1, pmReadout *in2) +void pmSubtractionThreadsInit(void) { if (threaded) { @@ -57,7 +24,4 @@ threaded = true; - - subtractionMutexInit(in1); - subtractionMutexInit(in2); { @@ -86,5 +50,5 @@ -void pmSubtractionThreadsFinalize(pmReadout *in1, pmReadout *in2) +void pmSubtractionThreadsFinalize(void) { if (!threaded) { @@ -97,6 +61,4 @@ psThreadTaskRemove("PSMODULES_SUBTRACTION_CONVOLVE"); - subtractionMutexDestroy(in1); - subtractionMutexDestroy(in1); return; } Index: /trunk/psModules/src/imcombine/pmSubtractionThreads.h =================================================================== --- /trunk/psModules/src/imcombine/pmSubtractionThreads.h (revision 27364) +++ /trunk/psModules/src/imcombine/pmSubtractionThreads.h (revision 27365) @@ -7,13 +7,12 @@ /// Set up threading for image matching /// -/// Sets up thread tasks, and initialises mutexes in readouts -void pmSubtractionThreadsInit(pmReadout *in1, pmReadout *in2 // Input readouts - ); +/// Sets up thread tasks +void pmSubtractionThreadsInit(void); /// Take down threading for image matching /// -/// Destroys thread tasks, and initialises mutexes in readouts -void pmSubtractionThreadsFinalize(pmReadout *in1, pmReadout *in2); +/// Destroys thread tasks +void pmSubtractionThreadsFinalize(void); #endif