Index: trunk/ppStack/src/ppStackMatch.c =================================================================== --- trunk/ppStack/src/ppStackMatch.c (revision 23379) +++ 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;