Index: trunk/ppStack/src/ppStack.h =================================================================== --- trunk/ppStack/src/ppStack.h (revision 23462) +++ trunk/ppStack/src/ppStack.h (revision 23573) @@ -7,4 +7,6 @@ #include #include + +#include "ppStackOptions.h" // Mask values for inputs @@ -106,11 +108,6 @@ /// Convolve image to match specified seeing bool ppStackMatch(pmReadout *readout, // Readout to be convolved; replaced with output - psArray **regions, // Array of regions used in each PSF matching, returned - psArray **kernels, // Array of kernels used in each PSF matching, returned - float *chi2, // Chi^2 from the stamps - float *weighting, // Stack weighting (1/noise^2) - psArray *sources, // Array of sources - const pmPSF *psf, // Target PSF - psRandom *rng, // Random number generator + ppStackOptions *options, // Options for stacking + int index, // Index of image to match const pmConfig *config // Configuration ); @@ -120,8 +117,7 @@ /// /// Corrects the source PSF photometry to a common system. Return the sum of the exposure times. -float ppStackSourcesTransparency(const psArray *sourceLists, // Sources for each input - psVector *inputMask, // Indicates bad input - const pmFPAview *view, // View to readout - const pmConfig *config // Configuration +bool ppStackSourcesTransparency(ppStackOptions *options, // Stacking options + const pmFPAview *view, // View to readout + const pmConfig *config // Configuration ); Index: trunk/ppStack/src/ppStackCombineInitial.c =================================================================== --- trunk/ppStack/src/ppStackCombineInitial.c (revision 23462) +++ trunk/ppStack/src/ppStackCombineInitial.c (revision 23573) @@ -16,5 +16,10 @@ psAssert(config, "Require configuration"); - psTimerStart("PPSTACK_FINAL"); + if (!options->convolve) { + // No need to do initial combination when we haven't convolved + return true; + } + + psTimerStart("PPSTACK_INITIAL"); psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe @@ -69,5 +74,4 @@ } - // call: ppStackReadoutInitial(config, outRO, readouts, subRegions, subKernels) psThreadJob *job = psThreadJobAlloc("PPSTACK_INITIAL_COMBINE"); // Job to start psArrayAdd(job->args, 1, thread); @@ -128,6 +132,6 @@ if (options->stats) { - psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_FINAL", 0, - "Time to make final stack", psTimerMark("PPSTACK_FINAL")); + psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_INITIAL", 0, + "Time to make final stack", psTimerMark("PPSTACK_INITIAL")); } Index: trunk/ppStack/src/ppStackConvolve.c =================================================================== --- trunk/ppStack/src/ppStackConvolve.c (revision 23462) +++ trunk/ppStack/src/ppStackConvolve.c (revision 23573) @@ -30,6 +30,6 @@ int num = options->num; // Number of inputs options->cells = psArrayAlloc(num); // Cells for convolved images --- a handle for reading again - options->subKernels = psArrayAlloc(num); // Subtraction kernels --- required in the stacking - options->subRegions = psArrayAlloc(num); // Subtraction regions --- required in the stacking + options->kernels = psArrayAlloc(num); // PSF-matching kernels --- required in the stacking + options->regions = psArrayAlloc(num); // PSF-matching regions --- required in the stacking int numGood = 0; // Number of good frames options->numCols = 0; @@ -38,5 +38,5 @@ psVectorInit(options->matchChi2, NAN); options->weightings = psVectorAlloc(num, PS_TYPE_F32); // Combination weightings for images (1/noise^2) - psVectorInit(options->weightings, NAN); + psVectorInit(options->weightings, 0.0); options->covariances = psArrayAlloc(num); // Covariance matrices @@ -78,12 +78,8 @@ // Background subtraction, scaling and normalisation is performed automatically by the image matching - psArray *regions = NULL, *kernels = NULL; // Regions and kernels used in subtraction psTimerStart("PPSTACK_MATCH"); - - if (!ppStackMatch(readout, ®ions, &kernels, &options->matchChi2->data.F32[i], - &options->weightings->data.F32[i], options->sourceLists->data[i], - options->psf, rng, config)) { + if (!ppStackMatch(readout, options, i, config)) { psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i); - options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_MATCH; + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_MATCH; psErrorClear(); continue; @@ -92,45 +88,48 @@ if (options->stats) { - pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis, - PM_SUBTRACTION_ANALYSIS_KERNEL);// Conv kernel psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_MATCH", PS_META_DUPLICATE_OK, "Time to match PSF", psTimerMark("PPSTACK_MATCH")); - psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK, - "Mean deviation for stamps", kernels->mean); - psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK, - "RMS deviation for stamps", kernels->rms); - psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK, - "Number of stamps", kernels->numStamps); - float deconv = psMetadataLookupF32(NULL, readout->analysis, - PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Deconvolution fraction - psMetadataAddF32(options->stats, PS_LIST_TAIL, "KERNEL.DECONV", PS_META_DUPLICATE_OK, - "Deconvolution fraction for kernel", deconv); psMetadataAddF32(options->stats, PS_LIST_TAIL, "PPSTACK.WEIGHTING", PS_META_DUPLICATE_OK, "Weighting for image", options->weightings->data.F32[i]); + + if (options->convolve) { + // Pull parameters out of convolution kernel + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, readout->analysis, + PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.MEAN", PS_META_DUPLICATE_OK, + "Mean deviation for stamps", kernels->mean); + psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.RMS", PS_META_DUPLICATE_OK, + "RMS deviation for stamps", kernels->rms); + psMetadataAddF32(options->stats, PS_LIST_TAIL, "STAMP.NUM", PS_META_DUPLICATE_OK, + "Number of stamps", kernels->numStamps); + float deconv = psMetadataLookupF32(NULL, readout->analysis, + PM_SUBTRACTION_ANALYSIS_DECONV_MAX); + psMetadataAddF32(options->stats, PS_LIST_TAIL, "KERNEL.DECONV", PS_META_DUPLICATE_OK, + "Deconvolution fraction for kernel", deconv); + } } psLogMsg("ppStack", PS_LOG_INFO, "Time to match image %d: %f sec", i, psTimerClear("PPSTACK_MATCH")); - options->subRegions->data[i] = regions; - options->subKernels->data[i] = kernels; - // Write the temporary convolved files - pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image - assert(hdu); - ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config); - psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask - pmConfigMaskWriteHeader(config, maskHeader); - ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config); - psFree(maskHeader); - psImageCovarianceTransfer(readout->variance, readout->covariance); - ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config); + if (options->convolve) { + pmHDU *hdu = readout->parent->parent->parent->hdu; // HDU for convolved image + assert(hdu); + ppStackWriteImage(options->imageNames->data[i], hdu->header, readout->image, config); + psMetadata *maskHeader = psMetadataCopy(NULL, hdu->header); // Copy of header, for mask + pmConfigMaskWriteHeader(config, maskHeader); + ppStackWriteImage(options->maskNames->data[i], maskHeader, readout->mask, config); + psFree(maskHeader); + psImageCovarianceTransfer(readout->variance, readout->covariance); + ppStackWriteImage(options->varianceNames->data[i], hdu->header, readout->variance, config); #ifdef TESTING - { - psString name = NULL; - psStringAppend(&name, "covariance_%d.fits", i); - ppStackWriteImage(name, hdu->header, readout->covariance->image, config); - pmStackVisualPlotTestImage(readout->covariance->image, name); - psFree(name); - } + { + psString name = NULL; + psStringAppend(&name, "covariance_%d.fits", i); + ppStackWriteImage(name, hdu->header, readout->covariance->image, config); + pmStackVisualPlotTestImage(readout->covariance->image, name); + psFree(name); + } #endif + } pmCell *inCell = readout->parent; // Input cell @@ -152,4 +151,5 @@ psFree(rng); + psFree(options->norm); options->norm = NULL; psFree(options->sourceLists); options->sourceLists = NULL; psFree(options->psf); options->psf = NULL; Index: trunk/ppStack/src/ppStackMatch.c =================================================================== --- trunk/ppStack/src/ppStackMatch.c (revision 23462) +++ trunk/ppStack/src/ppStackMatch.c (revision 23573) @@ -163,12 +163,9 @@ -bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, float *chi2, float *weighting, - psArray *sources, const pmPSF *psf, psRandom *rng, const pmConfig *config) +bool ppStackMatch(pmReadout *readout, ppStackOptions *options, int index, const pmConfig *config) { assert(readout); - assert(regions && !*regions); - assert(kernels && !*kernels); + assert(options); assert(config); - *weighting = 0.0; psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe @@ -197,81 +194,77 @@ } - pmReadout *output = pmReadoutAlloc(NULL); // Output readout, for holding results temporarily - - static int numInput = -1; // Index of input file - numInput++; + // Match the PSF + if (options->convolve) { + pmReadout *conv = pmReadoutAlloc(NULL); // Conv readout, for holding results temporarily #ifdef TESTING - // Read previously produced kernel - if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) { - const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root - assert(outName); - // Read convolution kernel - psString filename = NULL; // Output filename - psStringAppend(&filename, "%s.%d.kernel", outName, numInput); - psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename - psFree(filename); - psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel - psFree(resolved); - if (!fits || !pmReadoutReadSubtractionKernels(output, fits)) { - psError(PS_ERR_IO, false, "Unable to read previously produced kernel"); + // This is a hack to use the temporary convolved images and kernel generated previously. + // This makes the 'matching' operation much faster, allowing debugging of the stack process easier. + // It implicitly assumes the output root name is the same between invocations. + + // Read previously produced kernel + if (psMetadataLookupBool(NULL, config->arguments, "PPSTACK.DEBUG.STACK")) { + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.CONV.KERNEL", index); + psAssert(file, "Require file"); + + pmFPAview *view = pmFPAviewAlloc(0); // View to readout of interest + view->chip = view->cell = view->readout = 0; + psString filename = pmFPAfileNameFromRule(filerule->rule, file, view); // Filename of interest + + // Read convolution kernel + psString resolved = pmConfigConvertFilename(filename, config, false, false); // Resolved filename + psFree(filename); + psFits *fits = psFitsOpen(resolved, "r"); // FITS file for subtraction kernel + psFree(resolved); + if (!fits || !pmReadoutReadSubtractionKernels(conv, fits)) { + psError(PS_ERR_IO, false, "Unable to read previously produced kernel"); + psFitsClose(fits); + return false; + } psFitsClose(fits); - return false; - } - psFitsClose(fits); - - // Add in variance factor - pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, output->analysis, - PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels - float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor - psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR"); - if (!isfinite(vf)) { - vf = 1.0; - } - if (isfinite(vfItem->data.F32)) { - vfItem->data.F32 *= vf; - } else { - vfItem->data.F32 = vf; - } - - // Read image, mask, variance - const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for image - const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for mask - const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for variance map - psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images - psStringAppend(&imageName, "%s.%d.%s", outName, numInput, tempImage); - psStringAppend(&maskName, "%s.%d.%s", outName, numInput, tempMask); - psStringAppend(&varianceName, "%s.%d.%s", outName, numInput, tempVariance); - - if (!readImage(&readout->image, imageName, config) || !readImage(&readout->mask, maskName, config) || - !readImage(&readout->variance, varianceName, config)) { - psError(PS_ERR_IO, false, "Unable to read previously produced image."); + + // Add in variance factor + pmSubtractionKernels *kernels = psMetadataLookupPtr(NULL, conv->analysis, + PM_SUBTRACTION_ANALYSIS_KERNEL); // Kernels + float vf = pmSubtractionVarianceFactor(kernels, 0.0, 0.0, false); // Variance factor + psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR"); + if (!isfinite(vf)) { + vf = 1.0; + } + if (isfinite(vfItem->data.F32)) { + vfItem->data.F32 *= vf; + } else { + vfItem->data.F32 = vf; + } + + if (!readImage(&readout->image, options->imageNames->data[index], config) || + !readImage(&readout->mask, options->maskNames->data[index], config) || + !readImage(&readout->variance, options->varianceNames->data[index], config)) { + psError(PS_ERR_IO, false, "Unable to read previously produced image."); + psFree(imageName); + psFree(maskName); + psFree(varianceName); + return false; + } psFree(imageName); psFree(maskName); psFree(varianceName); - return false; - } - psFree(imageName); - psFree(maskName); - psFree(varianceName); - - psRegion *region = psMetadataLookupPtr(NULL, output->analysis, - PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region - - pmSubtractionAnalysis(readout->analysis, kernels, region, - readout->image->numCols, readout->image->numRows); - - psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel - psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // New covariance matrix - psFree(readout->covariance); - readout->covariance = covar; - psFree(kernel); - - } else { -#endif - - // Normal operations here - if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) { - assert(psf); - assert(sources); + + psRegion *region = psMetadataLookupPtr(NULL, conv->analysis, + PM_SUBTRACTION_ANALYSIS_REGION); // Convolution region + + pmSubtractionAnalysis(readout->analysis, kernels, region, + readout->image->numCols, readout->image->numRows); + + psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel + psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix + psFree(readout->covariance); + readout->covariance = covar; + psFree(kernel); + } else { +#endif + + // Normal operations here + psAssert(options->psf, "Require target PSF"); + psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list"); int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order @@ -284,6 +277,6 @@ float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold float sysError = psMetadataLookupF32(NULL, ppsub, "SYS"); // Relative systematic error in kernel - pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString( - psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE")); // Kernel type + const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type + pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type psVector *widths = psMetadataLookupPtr(NULL, ppsub, "ISIS.WIDTHS"); // ISIS Gaussian widths psVector *orders = psMetadataLookupPtr(NULL, ppsub, "ISIS.ORDERS"); // ISIS Polynomial orders @@ -308,41 +301,17 @@ } -#if 0 - // Testing the normalisation of the fake image - { - pmReadout *test = pmReadoutAlloc(NULL); // Test readout - psArray *sources = psArrayAlloc(1); // Array of sources - pmSource *source = pmSourceAlloc(); // Source - sources->data[0] = source; - source->peak = pmPeakAlloc(500, 500, 10000, PM_PEAK_LONE); - source->psfMag = -13.0; - pmReadoutFakeFromSources(test, 1000, 1000, sources, NULL, NULL, psf, 0.1, 0, false, true); - float sum = 0.0; - for (int y = 0; y < test->image->numRows; y++) { - for (int x = 0; x < test->image->numCols; x++) { - sum += test->image->data.F32[y][x]; - } - } - fprintf(stderr, "Photometry: %f --> %f = -13.0 ???\n", sum, -2.5*log10(sum)); - - psFits *fits = psFitsOpen("testphot.fits", "w"); - psFitsWriteImage(fits, NULL, test->image, 0, NULL); - psFitsClose(fits); - exit(0); - } -#endif - pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF // For the sake of stamps, remove nearby sources - psArray *stampSources = stackSourcesFilter(sources, footprint); // Filtered list of sources + psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index], + footprint); // Filtered list of sources if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, - stampSources, NULL, NULL, psf, NAN, footprint + size, + stampSources, NULL, NULL, options->psf, NAN, footprint + size, false, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF."); psFree(fake); psFree(optWidths); - psFree(output); + psFree(conv); return false; } @@ -357,5 +326,5 @@ pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; - psStringAppend(&name, "fake_%03d.fits", numInput); + psStringAppend(&name, "fake_%03d.fits", index); pmStackVisualPlotTestImage(fake->image, name); psFits *fits = psFitsOpen(name, "w"); @@ -367,5 +336,5 @@ pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; - psStringAppend(&name, "real_%03d.fits", numInput); + psStringAppend(&name, "real_%03d.fits", index); pmStackVisualPlotTestImage(readout->image, name); psFits *fits = psFitsOpen(name, "w"); @@ -384,5 +353,5 @@ PM_SUBTRACTION_ANALYSIS_KERNEL); // Conv kernel if (kernel) { - if (!pmSubtractionMatchPrecalc(output, NULL, readout, fake, readout->analysis, + if (!pmSubtractionMatchPrecalc(conv, NULL, readout, fake, readout->analysis, stride, sysError, maskVal, maskBad, maskPoor, poorFrac, badFrac)) { @@ -391,9 +360,9 @@ psFree(optWidths); psFree(stampSources); - psFree(output); + psFree(conv); return false; } } else { - if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, stride, regionSize, spacing, + if (!pmSubtractionMatch(conv, NULL, readout, fake, footprint, stride, regionSize, spacing, threshold, stampSources, stampsName, type, size, order, widths, orders, inner, ringsOrder, binning, penalty, @@ -405,5 +374,5 @@ psFree(optWidths); psFree(stampSources); - psFree(output); + psFree(conv); return false; } @@ -414,9 +383,9 @@ pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; - psStringAppend(&name, "conv_%03d.fits", numInput); - pmStackVisualPlotTestImage(output->image, name); + psStringAppend(&name, "conv_%03d.fits", index); + pmStackVisualPlotTestImage(conv->image, name); psFits *fits = psFitsOpen(name, "w"); psFree(name); - psFitsWriteImage(fits, hdu->header, output->image, 0, NULL); + psFitsWriteImage(fits, hdu->header, conv->image, 0, NULL); psFitsClose(fits); } @@ -424,9 +393,9 @@ pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; - psStringAppend(&name, "diff_%03d.fits", numInput); + psStringAppend(&name, "diff_%03d.fits", index); pmStackVisualPlotTestImage(fake->image, name); psFits *fits = psFitsOpen(name, "w"); psFree(name); - psBinaryOp(fake->image, output->image, "-", fake->image); + psBinaryOp(fake->image, conv->image, "-", fake->image); psFitsWriteImage(fits, hdu->header, fake->image, 0, NULL); psFitsClose(fits); @@ -444,5 +413,5 @@ // Set the variance factor psMetadataItem *vfItem = psMetadataLookup(readout->parent->concepts, "CELL.VARFACTOR"); - float vf = psMetadataLookupF32(NULL, output->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1); + float vf = psMetadataLookupF32(NULL, conv->analysis, PM_SUBTRACTION_ANALYSIS_VARFACTOR_1); if (!isfinite(vf)) { vf = 1.0; @@ -459,113 +428,90 @@ psFree(readout->variance); psFree(readout->covariance); - readout->image = psMemIncrRefCounter(output->image); - readout->mask = psMemIncrRefCounter(output->mask); - readout->variance = psMemIncrRefCounter(output->variance); - readout->covariance = psImageCovarianceTruncate(output->covariance, COVAR_FRAC); - } else { - // Fake the convolution - psRegion *region = psRegionAlloc(0, readout->image->numCols - 1, 0, readout->image->numRows - 1); - psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_REGION, - PS_DATA_REGION | PS_META_DUPLICATE_OK, "Fake subtraction region", region); - psFree(region); - pmSubtractionKernels *kernels = pmSubtractionKernelsPOIS(FAKE_SIZE, 0, penalty, - PM_SUBTRACTION_MODE_1); - // Set solution to delta function - kernels->solution1 = psVectorAlloc(kernels->num + 2, PS_TYPE_F64); - psVectorInit(kernels->solution1, 0.0); - int normIndex = PM_SUBTRACTION_INDEX_NORM(kernels); // Index for normalisation - kernels->solution1->data.F64[normIndex] = 1.0; - psMetadataAddPtr(output->analysis, PS_LIST_TAIL, PM_SUBTRACTION_ANALYSIS_KERNEL, - PS_DATA_UNKNOWN | PS_META_DUPLICATE_OK, "Fake subtraction kernel", kernels); - psFree(kernels); - } - + readout->image = psMemIncrRefCounter(conv->image); + readout->mask = psMemIncrRefCounter(conv->mask); + readout->variance = psMemIncrRefCounter(conv->variance); + readout->covariance = psImageCovarianceTruncate(conv->covariance, COVAR_FRAC); #ifdef TESTING - // Write convolution kernel + } +#endif + + // Extract the regions and solutions used in the image matching + // This stops them from being freed when we iterate back up the FPA + psArray *regions = options->regions->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match regions { - const char *outName = psMetadataLookupStr(NULL, config->arguments, "OUTPUT"); // Output root - assert(outName); - - psString filename = NULL; // Output filename - psStringAppend(&filename, "%s.%d.kernel", outName, numInput); - psString resolved = pmConfigConvertFilename(filename, config, true, false); // Resolved filename - psFree(filename); - psFits *fits = psFitsOpen(resolved, "w"); // FITS file for subtraction kernel - psFree(resolved); - pmReadoutWriteSubtractionKernels(output, fits); - psFitsClose(fits); - } - } -#endif - - readout->analysis = psMetadataCopy(readout->analysis, output->analysis); - -// Extract the regions and solutions used in the image matching -// This stops them from being freed when we iterate back up the FPA - *regions = psArrayAllocEmpty(ARRAY_BUFFER); // Array of regions - { - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION); - psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator - psFree(regex); - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_DATA_REGION); - *regions = psArrayAdd(*regions, ARRAY_BUFFER, item->data.V); - } - psFree(iter); - } - *kernels = psArrayAllocEmpty(ARRAY_BUFFER); // Array of kernels - { - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); - psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator - psFree(regex); - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_DATA_UNKNOWN); - // Set the normalisation dimensions, since these will be otherwise unavailable when reading the - // images by scans. - pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction - kernel->numCols = readout->image->numCols; - kernel->numRows = readout->image->numRows; - - *kernels = psArrayAdd(*kernels, ARRAY_BUFFER, kernel); - } - psFree(iter); - } - assert((*regions)->n == (*kernels)->n); - - // Record chi^2 - { - *chi2 = 0.0; - int num = 0; // Number of measurements of chi^2 - psString regex = NULL; // Regular expression - psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); - psMetadataIterator *iter = psMetadataIteratorAlloc(output->analysis, PS_LIST_HEAD, regex); // Iterator - psFree(regex); - psMetadataItem *item = NULL;// Item from iteration - while ((item = psMetadataGetAndIncrement(iter))) { - assert(item->type == PS_DATA_UNKNOWN); - pmSubtractionKernels *kernels = item->data.V; // Convolution kernels - *chi2 += kernels->mean; - num++; - } - psFree(iter); - *chi2 /= psImageCovarianceFactor(readout->covariance) * num; - } - - // Reject image completely if the maximum deconvolution fraction exceeds the limit - float deconv = psMetadataLookupF32(NULL, output->analysis, - PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Maximum deconvolution fraction - if (deconv > deconvLimit) { - psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n", - deconv, deconvLimit); - psFree(output); - return NULL; + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_REGION); + psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_REGION); + regions = psArrayAdd(regions, ARRAY_BUFFER, item->data.V); + } + psFree(iter); + } + psArray *kernels = options->kernels->data[index] = psArrayAllocEmpty(ARRAY_BUFFER); // Match kernels + { + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_UNKNOWN); + // Set the normalisation dimensions, since these will be otherwise unavailable when reading + // the images by scans. + pmSubtractionKernels *kernel = item->data.V; // Kernel used in subtraction + kernel->numCols = readout->image->numCols; + kernel->numRows = readout->image->numRows; + + kernels = psArrayAdd(kernels, ARRAY_BUFFER, kernel); + } + psFree(iter); + } + psAssert((regions)->n == (kernels)->n, "Number of match regions and kernels should match"); + + // Record chi^2 + { + double sum = 0.0; // Sum of chi^2 + int num = 0; // Number of measurements of chi^2 + psString regex = NULL; // Regular expression + psStringAppend(®ex, "^%s$", PM_SUBTRACTION_ANALYSIS_KERNEL); + psMetadataIterator *iter = psMetadataIteratorAlloc(conv->analysis, PS_LIST_HEAD, regex); + psFree(regex); + psMetadataItem *item = NULL;// Item from iteration + while ((item = psMetadataGetAndIncrement(iter))) { + assert(item->type == PS_DATA_UNKNOWN); + pmSubtractionKernels *kernels = item->data.V; // Convolution kernels + sum += kernels->mean; + num++; + } + psFree(iter); + options->matchChi2->data.F32[index] = sum / (psImageCovarianceFactor(readout->covariance) * num); + } + + // Reject image completely if the maximum deconvolution fraction exceeds the limit + float deconv = psMetadataLookupF32(NULL, conv->analysis, + PM_SUBTRACTION_ANALYSIS_DECONV_MAX); // Max deconvolution fraction + if (deconv > deconvLimit) { + psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting\n", + deconv, deconvLimit); + psFree(conv); + return NULL; + } + + readout->analysis = psMetadataCopy(readout->analysis, conv->analysis); + + psFree(conv); + } else { + // Match the normalisation + float norm = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation + psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(norm, PS_TYPE_F32)); + psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(norm), PS_TYPE_F32)); } // Ensure the background value is zero psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics for background + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator if (!psImageBackground(bg, NULL, readout->image, readout->mask, maskVal | maskBad, rng)) { psWarning("Can't measure background for image."); @@ -581,70 +527,15 @@ if (!psImageBackground(bg, NULL, readout->variance, readout->mask, maskVal | maskBad, rng)) { psError(PS_ERR_UNKNOWN, false, "Can't measure mean variance for image."); - psFree(output); + psFree(rng); + psFree(bg); return false; } - *weighting = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * - psImageCovarianceFactor(readout->covariance)); + options->weightings->data.F32[index] = 1.0 / (psStatsGetValue(bg, PS_STAT_ROBUST_MEDIAN) * + psImageCovarianceFactor(readout->covariance)); psMetadataAddF32(readout->analysis, PS_LIST_TAIL, "PPSTACK.WEIGHTING", 0, - "Weighting by 1/noise^2 for stack", *weighting); - + "Weighting by 1/noise^2 for stack", options->weightings->data.F32[index]); + + psFree(rng); psFree(bg); - -#if 0 -#define RADIUS 10 // Radius of photometry -#define MIN_ERR 0.05 // Minimum photometric error, mag -#define MAX_MAG -13 // Maximum magnitude for source - - // Ensure the normalisation is correct - // XXX Ideally, would like to do proper PSF-fit photometry, but will settle for aperture photometry - int numSources = sources->n; // Number of sources - psVector *ratio = psVectorAlloc(numSources, PS_TYPE_F32); // Flux ratios for sources - psVector *ratioMask = psVectorAlloc(numSources, PS_TYPE_MASK); // Mask for flux ratios - psVectorInit(ratioMask, 0xFF); - psImage *image = readout->image; // Image of interest - psImage *mask = readout->mask; // Mask of interest - int numCols = image->numCols, numRows = image->numRows; // Size of image - for (int i = 0; i < numSources; i++) { - pmSource *source = sources->data[i]; // Source of interest - if (!source || source->mode & SOURCE_MASK || !isfinite(source->psfMag) || !isfinite(source->errMag) || - source->errMag > MIN_ERR || source->psfMag > MAX_MAG) { - continue; - } - - float xSrc, ySrc; // Source coordinates - coordsFromSource(&xSrc, &ySrc, source); - int xMin = PS_MAX(0, xSrc - RADIUS), xMax = PS_MIN(numCols - 1, xSrc + RADIUS); // Bounds in x - int yMin = PS_MAX(0, ySrc - RADIUS), yMax = PS_MIN(numRows - 1, ySrc + RADIUS); // Bounds in y - int numPix = 0; // Number of pixels - float sum = 0.0; // Sum of pixels - for (int y = yMin; y <= yMax; y++) { - for (int x = xMin; x <= xMax; x++) { - if (mask->data.PS_TYPE_MASK_DATA[y][x] & maskBad) { - continue; - } - sum += image->data.F32[y][x]; - numPix++; - } - } - if (sum >= 0 && numPix > 0) { - float mag = -2.5 * log10(sum * M_PI * PS_SQR(RADIUS) / numPix); // Instrumental magnitude - ratio->data.F32[i] = mag - source->psfMag; - ratioMask->data.PS_TYPE_MASK_DATA[i] = 0; - } - } - - psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // Statistics - if (!psVectorStats(stats, ratio, NULL, ratioMask, 0xFF)) { - psWarning("Unable to measure normalisation --- assuming correct."); - } else { - psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image by %f (+/- %f) mag\n", - stats->robustMedian, stats->robustStdev); - float norm = powf(10.0, -0.4 * stats->robustMedian); // Normalisation to apply - psBinaryOp(image, image, "*", psScalarAlloc(norm, PS_TYPE_F32)); - } - psFree(stats); - psFree(ratio); - psFree(ratioMask); -#endif #ifdef TESTING @@ -652,14 +543,12 @@ pmHDU *hdu = pmHDUFromCell(readout->parent); psString name = NULL; - psStringAppend(&name, "convolved_%03d.fits", numInput); - pmStackVisualPlotTestImage(output->image, name); + psStringAppend(&name, "convolved_%03d.fits", index); + pmStackVisualPlotTestImage(readout->image, name); psFits *fits = psFitsOpen(name, "w"); psFree(name); - psFitsWriteImage(fits, hdu->header, output->image, 0, NULL); + psFitsWriteImage(fits, hdu->header, readout->image, 0, NULL); psFitsClose(fits); } #endif - - psFree(output); return true; Index: trunk/ppStack/src/ppStackOptions.c =================================================================== --- trunk/ppStack/src/ppStackOptions.c (revision 23462) +++ trunk/ppStack/src/ppStackOptions.c (revision 23573) @@ -18,7 +18,8 @@ psFree(options->inputMask); psFree(options->sourceLists); + psFree(options->norm); psFree(options->cells); - psFree(options->subKernels); - psFree(options->subRegions); + psFree(options->kernels); + psFree(options->regions); psFree(options->matchChi2); psFree(options->weightings); @@ -35,4 +36,5 @@ psMemSetDeallocator(options, (psFreeFunc)stackOptionsFree); + options->convolve = true; options->stats = NULL; options->statsFile = NULL; @@ -44,7 +46,8 @@ options->inputMask = NULL; options->sourceLists = NULL; + options->norm = NULL; options->cells = NULL; - options->subKernels = NULL; - options->subRegions = NULL; + options->kernels = NULL; + options->regions = NULL; options->numCols = 0; options->numRows = 0; Index: trunk/ppStack/src/ppStackOptions.h =================================================================== --- trunk/ppStack/src/ppStackOptions.h (revision 23462) +++ trunk/ppStack/src/ppStackOptions.h (revision 23573) @@ -8,4 +8,5 @@ typedef struct { // Setup + bool convolve; // Convolve images? psMetadata *stats; // Statistics for output FILE *statsFile; // File to which to write statistics @@ -17,8 +18,9 @@ psVector *inputMask; // Mask for inputs psArray *sourceLists; // Individual lists of sources for matching + psVector *norm; // Normalisation for each image // Convolve psArray *cells; // Cells for convolved images --- a handle for reading again - psArray *subKernels; // Subtraction kernels --- required in the stacking - psArray *subRegions; // Subtraction regions --- required in the stacking + psArray *kernels; // PSF-matching kernels --- required in the stacking + psArray *regions; // PSF-matching regions --- required in the stacking int numCols, numRows; // Size of image psVector *matchChi2; // chi^2 for stamps from matching Index: trunk/ppStack/src/ppStackPrepare.c =================================================================== --- trunk/ppStack/src/ppStackPrepare.c (revision 23462) +++ trunk/ppStack/src/ppStackPrepare.c (revision 23573) @@ -127,22 +127,22 @@ psVectorInit(options->inputMask, 0); - if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) { - pmFPAfileActivate(config->files, false, NULL); - ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true); - pmFPAview *view = ppStackFilesIterateDown(config); - if (!view) { - return false; - } - - // Generate list of PSFs - psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, - "^PPSTACK.INPUT$"); // Iterator - psMetadataItem *fileItem; // Item from iteration - psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope - int numCols = 0, numRows = 0; // Size of image - int index = 0; // Index for file - while ((fileItem = psMetadataGetAndIncrement(fileIter))) { - assert(fileItem->type == PS_DATA_UNKNOWN); - pmFPAfile *inputFile = fileItem->data.V; // An input file + pmFPAfileActivate(config->files, false, NULL); + ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true); + pmFPAview *view = ppStackFilesIterateDown(config); + if (!view) { + return false; + } + + psMetadataIterator *fileIter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, "^PPSTACK.INPUT$"); + psMetadataItem *fileItem; // Item from iteration + psArray *psfs = psArrayAlloc(num); // PSFs for PSF envelope + int numCols = 0, numRows = 0; // Size of image + int index = 0; // Index for file + while ((fileItem = psMetadataGetAndIncrement(fileIter))) { + assert(fileItem->type == PS_DATA_UNKNOWN); + pmFPAfile *inputFile = fileItem->data.V; // An input file + + // Get list of PSFs, to determine target PSF + if (options->convolve) { pmChip *chip = pmFPAviewThisChip(view, inputFile->fpa); // The chip: holds the PSF pmPSF *psf = psMetadataLookupPtr(NULL, chip->analysis, "PSPHOT.PSF"); // PSF @@ -173,57 +173,51 @@ numRows = naxis2; } - - pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources - psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources - if (!sources) { - psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout."); - return NULL; - } - options->sourceLists->data[index] = psMemIncrRefCounter(sources); - - // Re-do photometry if we don't trust the source lists - if (psMetadataLookupBool(NULL, recipe, "PHOT")) { - psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num); - pmFPAfileActivate(config->files, false, NULL); - ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index); - pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File - pmFPAview *photView = ppStackFilesIterateDown(config); - if (!photView) { - psFree(view); - return false; - } - - pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest - - if (!ppStackInputPhotometer(ro, sources, config)) { - psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources"); - psFree(view); - psFree(photView); - return false; - } - + } + + + pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources + psArray *sources = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.SOURCES"); // Sources + if (!sources) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find sources in readout."); + return NULL; + } + options->sourceLists->data[index] = psMemIncrRefCounter(sources); + + // Re-do photometry if we don't trust the source lists + if (psMetadataLookupBool(NULL, recipe, "PHOT")) { + psTrace("ppStack", 2, "Photometering input %d of %d....\n", index, num); + pmFPAfileActivate(config->files, false, NULL); + ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index); + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File + pmFPAview *photView = ppStackFilesIterateDown(config); + if (!photView) { + psFree(view); + return false; + } + + pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest + + if (!ppStackInputPhotometer(ro, sources, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources"); + psFree(view); psFree(photView); - if (!ppStackFilesIterateUp(config)) { - psFree(view); - return false; - } - pmFPAfileActivate(config->files, false, NULL); - ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true); - } - - index++; - } - psFree(fileIter); - - // Zero point calibration - options->sumExposure = ppStackSourcesTransparency(options->sourceLists, options->inputMask, - view, config); - if (!isfinite(options->sumExposure) || options->sumExposure <= 0) { - psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences"); - psFree(view); - return false; - } - - // Generate target PSF + return false; + } + + psFree(photView); + if (!ppStackFilesIterateUp(config)) { + psFree(view); + return false; + } + pmFPAfileActivate(config->files, false, NULL); + ppStackFileActivation(config, PPSTACK_FILES_PREPARE, true); + } + + index++; + } + psFree(fileIter); + + // Generate target PSF + if (options->convolve) { options->psf = ppStackPSF(config, numCols, numRows, psfs, options->inputMask); psFree(psfs); @@ -240,10 +234,16 @@ "Target PSF", options->psf); outChip->data_exists = true; - + } + + // Zero point calibration + if (!ppStackSourcesTransparency(options, view, config)) { + psError(PS_ERR_UNKNOWN, false, "Unable to calculate transparency differences"); psFree(view); - - if (!ppStackFilesIterateUp(config)) { - return false; - } + return false; + } + psFree(view); + + if (!ppStackFilesIterateUp(config)) { + return false; } Index: trunk/ppStack/src/ppStackReject.c =================================================================== --- trunk/ppStack/src/ppStackReject.c (revision 23462) +++ trunk/ppStack/src/ppStackReject.c (revision 23573) @@ -14,4 +14,9 @@ psAssert(options, "Require options"); psAssert(config, "Require configuration"); + + if (!options->convolve) { + // No need to do complicated rejection when we haven't convolved + return true; + } int num = options->num; // Number of inputs @@ -86,6 +91,6 @@ psPixels *reject = pmStackReject(options->inspect->data[i], options->numCols, options->numRows, - threshold, poorFrac, stride, options->subRegions->data[i], - options->subKernels->data[i]); // Rejected pixels + threshold, poorFrac, stride, options->regions->data[i], + options->kernels->data[i]); // Rejected pixels #ifdef TESTING @@ -141,6 +146,6 @@ psFree(options->inspect); options->inspect = NULL; - psFree(options->subKernels); options->subKernels = NULL; - psFree(options->subRegions); options->subRegions = NULL; + psFree(options->kernels); options->kernels = NULL; + psFree(options->regions); options->regions = NULL; if (numRejected >= num - 1) { Index: trunk/ppStack/src/ppStackSetup.c =================================================================== --- trunk/ppStack/src/ppStackSetup.c (revision 23462) +++ trunk/ppStack/src/ppStackSetup.c (revision 23573) @@ -11,4 +11,34 @@ #include "ppStackLoop.h" +#define BUFFER 16 // Buffer for name array + +// Generate an array of input filenames +static psArray *stackNameArray(pmConfig *config, // Configuration + const char *name // Name of file + ) +{ + psAssert(config, "Require configuration"); + psAssert(config, "Require file name"); + + psString regex = NULL; // Regular expression for selecting file + psStringAppend(®ex, "^%s$", name); + psMetadataIterator *iter = psMetadataIteratorAlloc(config->files, PS_LIST_HEAD, regex); // Iterator + psFree(regex); + psMetadataItem *item; // Item from iteration + psArray *array = psArrayAlloc(BUFFER); // Array of file names + while ((item = psMetadataGetAndIncrement(iter))) { + psAssert(item->type == PS_DATA_UNKNOWN, "Should be this type"); + pmFPAfile *file = item->data.V; // An input file + + psString filename = psStringCopy(file->filename); // Filename of interest + psArrayAdd(array, array->n, filename); + psFree(filename); // Drop reference + } + psFree(iter); + + return array; +} + + bool ppStackSetup(ppStackOptions *options, pmConfig *config) { @@ -18,4 +48,10 @@ psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe psAssert(recipe, "We've thrown an error on this before."); + + options->convolve = psMetadataLookupBool(NULL, config->arguments, "CONVOLVE"); // Convolve images? + if (!psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) { + psWarning("No PSFs provided --- unable to convolve to common PSF."); + options->convolve = false; + } int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs @@ -36,43 +72,50 @@ } - const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images - if (!tempDir) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe"); - return false; + // Generate temporary names for convolved images + if (options->convolve) { + const char *tempDir = psMetadataLookupStr(NULL, recipe, "TEMP.DIR"); // Directory for temporary images + if (!tempDir) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to find TEMP.DIR in recipe"); + return false; + } + + psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments, + "OUTPUT")); // Name for temporary files + const char *tempName = basename(outputName); + if (!tempName) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files."); + psFree(outputName); + return false; + } + + const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for images + const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for masks + const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for var maps + if (!tempImage || !tempMask || !tempVariance) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe"); + psFree(outputName); + return false; + } + + options->imageNames = psArrayAlloc(num); + options->maskNames = psArrayAlloc(num); + options->varianceNames = psArrayAlloc(num); + for (int i = 0; i < num; i++) { + psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images + psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage); + psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask); + psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance); + psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName); + options->imageNames->data[i] = imageName; + options->maskNames->data[i] = maskName; + options->varianceNames->data[i] = varianceName; + } + psFree(outputName); + } else { + options->imageNames = stackNameArray(config, "PPSTACK.INPUT"); + options->maskNames = stackNameArray(config, "PPSTACK.INPUT.MASK"); + options->varianceNames = stackNameArray(config, "PPSTACK.INPUT.VARIANCE"); } - - psString outputName = psStringCopy(psMetadataLookupStr(NULL, config->arguments, - "OUTPUT")); // Name for temporary files - const char *tempName = basename(outputName); - if (!tempName) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Unable to construct basename for temporary files."); - psFree(outputName); - return false; - } - - const char *tempImage = psMetadataLookupStr(NULL, recipe, "TEMP.IMAGE"); // Suffix for temporary images - const char *tempMask = psMetadataLookupStr(NULL, recipe, "TEMP.MASK"); // Suffix for temporary masks - const char *tempVariance = psMetadataLookupStr(NULL, recipe, "TEMP.VARIANCE"); // Suffix for temp var maps - if (!tempImage || !tempMask || !tempVariance) { - psError(PS_ERR_BAD_PARAMETER_VALUE, false, - "Unable to find TEMP.IMAGE, TEMP.MASK and TEMP.VARIANCE in recipe"); - psFree(outputName); - return false; - } - - options->imageNames = psArrayAlloc(num); - options->maskNames = psArrayAlloc(num); - options->varianceNames = psArrayAlloc(num); - for (int i = 0; i < num; i++) { - psString imageName = NULL, maskName = NULL, varianceName = NULL; // Names for convolved images - psStringAppend(&imageName, "%s/%s.%d.%s", tempDir, tempName, i, tempImage); - psStringAppend(&maskName, "%s/%s.%d.%s", tempDir, tempName, i, tempMask); - psStringAppend(&varianceName, "%s/%s.%d.%s", tempDir, tempName, i, tempVariance); - psTrace("ppStack", 5, "Temporary files: %s %s %s\n", imageName, maskName, varianceName); - options->imageNames->data[i] = imageName; - options->maskNames->data[i] = maskName; - options->varianceNames->data[i] = varianceName; - } - psFree(outputName); if (!pmConfigMaskSetBits(NULL, NULL, config)) { Index: trunk/ppStack/src/ppStackSources.c =================================================================== --- trunk/ppStack/src/ppStackSources.c (revision 23462) +++ trunk/ppStack/src/ppStackSources.c (revision 23573) @@ -9,7 +9,9 @@ //#define TESTING // Enable debugging output +#ifdef TESTING // Size of fake image; set by hand because it's trouble to get it from other places #define FAKE_COLS 4861 #define FAKE_ROWS 4913 +#endif #ifdef TESTING @@ -53,13 +55,17 @@ -float ppStackSourcesTransparency(const psArray *sourceLists, psVector *inputMask, - const pmFPAview *view, const pmConfig *config) +bool ppStackSourcesTransparency(ppStackOptions *options, const pmFPAview *view, const pmConfig *config) { - PS_ASSERT_ARRAY_NON_NULL(sourceLists, NAN); - PS_ASSERT_VECTOR_NON_NULL(inputMask, NAN); - PS_ASSERT_VECTOR_TYPE(inputMask, PS_TYPE_U8, NAN); - PS_ASSERT_VECTOR_SIZE(inputMask, sourceLists->n, NAN); - PS_ASSERT_PTR_NON_NULL(view, NAN); - PS_ASSERT_PTR_NON_NULL(config, NAN); + PS_ASSERT_PTR_NON_NULL(options, false); + PS_ASSERT_PTR_NON_NULL(view, false); + PS_ASSERT_PTR_NON_NULL(config, false); + + psArray *sourceLists = options->sourceLists; // Source lists for each input + psVector *inputMask = options->inputMask; // Mask for inputs + + PS_ASSERT_ARRAY_NON_NULL(sourceLists, false); + PS_ASSERT_VECTOR_NON_NULL(inputMask, false); + PS_ASSERT_VECTOR_TYPE(inputMask, PS_TYPE_U8, false); + PS_ASSERT_VECTOR_SIZE(inputMask, sourceLists->n, false); #if defined(TESTING) && 0 @@ -100,5 +106,5 @@ if (!airmassZP) { psError(PS_ERR_UNKNOWN, false, "Unable to find ZP.AIRMASS in recipe."); - return NAN; + return false; } @@ -139,5 +145,5 @@ exptime, airmass, expFilter); psFree(zp); - return NAN; + return false; } @@ -149,10 +155,10 @@ "Unable to find airmass term (ZP.AIRMASS) for filter %s", filter); psFree(zp); - return NAN; + return false; } } else if (strcmp(filter, expFilter) != 0) { psError(PS_ERR_BAD_PARAMETER_VALUE, false, "Filters don't match: %s vs %s", filter, expFilter); psFree(zp); - return NAN; + return false; } @@ -160,4 +166,6 @@ sumExpTime += exptime; } + + options->sumExposure = sumExpTime; psArray *matches = pmSourceMatchSources(sourceLists, radius, true); // List of matches @@ -165,5 +173,5 @@ psError(PS_ERR_UNKNOWN, false, "Unable to match sources"); psFree(zp); - return NAN; + return false; } @@ -177,5 +185,5 @@ if (!trans) { psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies"); - return NAN; + return false; } @@ -186,5 +194,5 @@ for (int i = 0; i < trans->n; i++) { if (!isfinite(trans->data.F32[i])) { - inputMask->data.U8[i] = PPSTACK_MASK_CAL; + inputMask->data.U8[i] |= PPSTACK_MASK_CAL; } } @@ -223,4 +231,5 @@ // m_0 = m_1 + zp_1 + trans_1 - c1 * airmass_0 - 2.5log(t_0) // We don't need to know the magnitude zero point for the filter, since it cancels out + options->norm = psVectorAlloc(num, PS_TYPE_F32); for (int i = 0; i < num; i++) { if (!isfinite(trans->data.F32[i])) { @@ -229,4 +238,5 @@ psArray *sources = sourceLists->data[i]; // Sources of interest float magCorr = airmassTerm - 2.5*log10(sumExpTime) - zp->data.F32[i] - trans->data.F32[i]; + options->norm->data.F32[i] = magCorr; psLogMsg("ppStack", PS_LOG_INFO, "Applying magnitude correction to image %d: %f\n", i, magCorr); @@ -248,5 +258,5 @@ psError(PS_ERR_UNKNOWN, false, "Unable to match sources"); psFree(zp); - return NAN; + return false; } psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej, @@ -259,4 +269,4 @@ #endif - return sumExpTime; + return true; }