Index: trunk/ppStack/src/Makefile.am =================================================================== --- trunk/ppStack/src/Makefile.am (revision 28252) +++ trunk/ppStack/src/Makefile.am (revision 28253) @@ -48,4 +48,5 @@ ppStackPhotometry.c \ ppStackFinish.c \ + ppStackTarget.c \ ppStackErrorCodes.c Index: trunk/ppStack/src/ppStack.h =================================================================== --- trunk/ppStack/src/ppStack.h (revision 28252) +++ trunk/ppStack/src/ppStack.h (revision 28253) @@ -26,4 +26,5 @@ typedef enum { PPSTACK_FILES_PREPARE, // Files for preparation + PPSTACK_FILES_TARGET, // Files for target generation PPSTACK_FILES_CONVOLVE, // Files for convolution PPSTACK_FILES_STACK, // Stack files @@ -115,6 +116,12 @@ void ppStackVersionPrint(void); +/// Generate target PSF image +psImage *ppStackTarget(ppStackOptions *options, // Options for stacking + pmConfig *config // Configuration + ); + /// Convolve image to match specified seeing bool ppStackMatch(pmReadout *readout, // Readout to be convolved; replaced with output + const psImage *target, // Target PSF image ppStackOptions *options, // Options for stacking int index, // Index of image to match Index: trunk/ppStack/src/ppStackConvolve.c =================================================================== --- trunk/ppStack/src/ppStackConvolve.c (revision 28252) +++ trunk/ppStack/src/ppStackConvolve.c (revision 28253) @@ -45,4 +45,13 @@ options->convCovars = psArrayAlloc(num); // Covariance matrices + psImage *target = NULL; // Target PSF image + if (options->convolve) { + target = ppStackTarget(options, config); + if (!target) { + psError(psErrorCodeLast(), false, "Unable to produce stack target image"); + return false; + } + } + psVector *renorms = psVectorAlloc(num, PS_TYPE_F32); // Renormalisation values for variances psVectorInit(renorms, NAN); @@ -71,4 +80,5 @@ psFree(fpaList); psFree(cellList); + psFree(target); return false; } @@ -87,4 +97,5 @@ psFree(fpaList); psFree(cellList); + psFree(target); return false; } @@ -93,5 +104,5 @@ psTimerStart("PPSTACK_MATCH"); options->origCovars->data[i] = psMemIncrRefCounter(readout->covariance); - if (!ppStackMatch(readout, options, i, config)) { + if (!ppStackMatch(readout, target, options, i, config)) { // XXX many things can cause a failure of ppStackMatch -- should some be handled differently? psErrorCode error = psErrorCodeLast(); // Error code @@ -102,4 +113,8 @@ case PPSTACK_ERR_IO: psError(error, false, "Unable to match image %d due to fatal error.", i); + psFree(rng); + psFree(fpaList); + psFree(cellList); + psFree(target); return false; // Non-fatal errors @@ -154,4 +169,5 @@ psFree(cellList); psFree(rng); + psFree(target); return false; } @@ -164,4 +180,5 @@ psFree(rng); psFree(maskHeader); + psFree(target); return false; } @@ -172,4 +189,5 @@ psFree(cellList); psFree(rng); + psFree(target); return false; } @@ -222,4 +240,5 @@ psFree(cellList); psFree(rng); + psFree(target); return false; } @@ -229,4 +248,5 @@ } psFree(rng); + psFree(target); psFree(options->sourceLists); options->sourceLists = NULL; Index: trunk/ppStack/src/ppStackFiles.c =================================================================== --- trunk/ppStack/src/ppStackFiles.c (revision 28252) +++ trunk/ppStack/src/ppStackFiles.c (revision 28253) @@ -16,8 +16,9 @@ static char *filesPrepare[] = { "PPSTACK.INPUT.PSF", "PPSTACK.INPUT.SOURCES", "PPSTACK.TARGET.PSF", NULL }; +/// Files required for generating convolution target +static char *filesTarget[] = { "PPSTACK.INPUT.VARIANCE", "PPSTACK.INPUT.MASK", NULL }; + /// Files required for the convolution static char *filesConvolve[] = { "PPSTACK.INPUT", "PPSTACK.INPUT.MASK", "PPSTACK.INPUT.VARIANCE", NULL }; - -// "PPSTACK.CONV.KERNEL", NULL }; /// Regular (convolved) stack files @@ -41,4 +42,5 @@ switch (list) { case PPSTACK_FILES_PREPARE: return filesPrepare; + case PPSTACK_FILES_TARGET: return filesTarget; case PPSTACK_FILES_CONVOLVE: return filesConvolve; case PPSTACK_FILES_STACK: return filesStack; Index: trunk/ppStack/src/ppStackMatch.c =================================================================== --- trunk/ppStack/src/ppStackMatch.c (revision 28252) +++ trunk/ppStack/src/ppStackMatch.c (revision 28253) @@ -13,7 +13,4 @@ #define MAG_IGNORE 50 // Ignore magnitudes fainter than this --- they're not real! #define FAKE_SIZE 1 // Size of fake convolution kernel -#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \ - PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources -#define NOISE_FRACTION 0.01 // Set minimum flux to this fraction of noise #define COVAR_FRAC 0.01 // Truncation fraction for covariance matrix @@ -49,77 +46,4 @@ #endif -// Get coordinates from a source -static void coordsFromSource(float *x, float *y, const pmSource *source) -{ - assert(x && y); - assert(source); - - if (source->modelPSF) { - *x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; - *y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; - } else { - *x = source->peak->xf; - *y = source->peak->yf; - } - return; -} - -static psArray *stackSourcesFilter(psArray *sources, // Source list to filter - int exclusion // Exclusion zone, pixels - ) -{ - psAssert(sources && sources->n > 0, "Require array of sources"); - if (exclusion <= 0) { - return psMemIncrRefCounter(sources); - } - - int num = sources->n; // Number of sources - psVector *x = psVectorAlloc(num, PS_TYPE_F32), *y = psVectorAlloc(num, PS_TYPE_F32); // Coordinates - int numGood = 0; // Number of good sources - for (int i = 0; i < num; i++) { - pmSource *source = sources->data[i]; // Source of interest - if (!source) { - continue; - } - coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], source); - numGood++; - } - x->n = y->n = numGood; - - psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); // kd tree - - psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources - psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of source - int numFiltered = 0; // Number of filtered sources - for (int i = 0; i < num; i++) { - pmSource *source = sources->data[i]; // Source of interest - if (!source) { - continue; - } - float xSource, ySource; // Coordinates of source - coordsFromSource(&xSource, &ySource, source); - - coords->data.F64[0] = xSource; - coords->data.F64[1] = ySource; - - long numWithin = psTreeWithin(tree, coords, exclusion); // Number within exclusion zone - psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone", - coords->data.F64[0], coords->data.F64[1], numWithin); - if (numWithin == 1) { - // Only itself inside the exclusion zone - filtered = psArrayAdd(filtered, filtered->n, source); - } else { - numFiltered++; - } - } - psFree(coords); - psFree(tree); - psFree(x); - psFree(y); - - psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood); - - return filtered; -} // Add background into the fake image @@ -202,5 +126,6 @@ -bool ppStackMatch(pmReadout *readout, ppStackOptions *options, int index, const pmConfig *config) +bool ppStackMatch(pmReadout *readout, const psImage *target, ppStackOptions *options, + int index, const pmConfig *config) { assert(readout); @@ -286,8 +211,6 @@ } else { #endif - // Normal operations here - psAssert(options->psf, "Require target PSF"); - psAssert(options->sourceLists && options->sourceLists->data[index], "Require source list"); + psAssert(target, "Require target PSF image"); int order = psMetadataLookupS32(NULL, ppsub, "SPATIAL.ORDER"); // Spatial polynomial order @@ -341,37 +264,6 @@ } - pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF - - psStats *bg = psStatsAlloc(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)) { - psError(PPSTACK_ERR_DATA, false, "Can't measure background for image."); - psFree(fake); - psFree(optWidths); - psFree(conv); - psFree(bg); - psFree(rng); - return false; - } - float minFlux = NOISE_FRACTION * bg->robustStdev; // Minimum flux level for fake image - psFree(rng); - psFree(bg); - - // For the sake of stamps, remove nearby sources - psArray *stampSources = stackSourcesFilter(options->sourceLists->data[index], - footprint); // Filtered list of sources - - bool oldThreads = pmReadoutFakeThreads(true); // Old threading state - if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, - stampSources, SOURCE_MASK, NULL, NULL, options->psf, - minFlux, footprint + size, false, true)) { - psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF."); - psFree(fake); - psFree(optWidths); - psFree(conv); - return false; - } - pmReadoutFakeThreads(oldThreads); - + pmReadout *fake = pmReadoutAlloc(NULL); + fake->image = psImageCopy(NULL, target, PS_TYPE_F32); fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK); @@ -420,5 +312,4 @@ psFree(fake); psFree(optWidths); - psFree(stampSources); psFree(conv); if (threads > 0) { @@ -436,5 +327,4 @@ psFree(fake); psFree(optWidths); - psFree(stampSources); psFree(conv); psFree(widthsCopy); @@ -446,5 +336,5 @@ if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing, - threshold, stampSources, stampsName, type, size, order, widthsCopy, + threshold, options->sources, stampsName, type, size, order, widthsCopy, orders, inner, ringsOrder, binning, penalty, optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, @@ -454,5 +344,4 @@ psFree(fake); psFree(optWidths); - psFree(stampSources); psFree(conv); psFree(widthsCopy); @@ -492,5 +381,4 @@ psFree(fake); psFree(optWidths); - psFree(stampSources); if (threads > 0) { Index: trunk/ppStack/src/ppStackTarget.c =================================================================== --- trunk/ppStack/src/ppStackTarget.c (revision 28253) +++ trunk/ppStack/src/ppStackTarget.c (revision 28253) @@ -0,0 +1,213 @@ +#include +#include +#include + +#include "ppStack.h" + +#define SOURCE_MASK (PM_SOURCE_MODE_FAIL | PM_SOURCE_MODE_DEFECT | PM_SOURCE_MODE_SATURATED | \ + PM_SOURCE_MODE_CR_LIMIT | PM_SOURCE_MODE_EXT_LIMIT) // Mask to apply to input sources + +// Get coordinates from a source +static void coordsFromSource(float *x, float *y, float *mag, const pmSource *source) +{ + assert(x && y); + assert(source); + + if (source->modelPSF) { + *x = source->modelPSF->params->data.F32[PM_PAR_XPOS]; + *y = source->modelPSF->params->data.F32[PM_PAR_YPOS]; + } else { + *x = source->peak->xf; + *y = source->peak->yf; + } + if (mag) { + *mag = source->psfMag; + } + return; +} + +// Filter a list of sources to exclude sources with near neighbours +static psArray *stackSourcesFilter(psArray *sources, // Source list to filter + int exclusion, // Exclusion zone, pixels + float minMagDiff // Minimum magnitude difference + ) +{ + psAssert(sources && sources->n > 0, "Require array of sources"); + if (exclusion <= 0) { + return psMemIncrRefCounter(sources); + } + + int num = sources->n; // Number of sources + psVector *x = psVectorAlloc(num, PS_TYPE_F32), *y = psVectorAlloc(num, PS_TYPE_F32); // Coordinates + psVector *mag = psVectorAlloc(num, PS_TYPE_F32); // Magnitudes + int numGood = 0; // Number of good sources + for (int i = 0; i < num; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source) { + continue; + } + coordsFromSource(&x->data.F32[numGood], &y->data.F32[numGood], &mag->data.F32[numGood], source); + numGood++; + } + x->n = y->n = mag->n = numGood; + + psTree *tree = psTreePlant(2, 2, PS_TREE_EUCLIDEAN, x, y); // kd tree + + psArray *filtered = psArrayAllocEmpty(numGood); // Filtered list of sources + psVector *coords = psVectorAlloc(2, PS_TYPE_F64); // Coordinates of source + int numFiltered = 0; // Number of filtered sources + for (int i = 0; i < num; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (!source) { + continue; + } + float xSource, ySource; // Coordinates of source + coordsFromSource(&xSource, &ySource, NULL, source); + + coords->data.F64[0] = xSource; + coords->data.F64[1] = ySource; + + psVector *indices = psTreeAllWithin(tree, coords, exclusion); // Number within exclusion zone + psTrace("ppStack", 9, "Source at %.0lf,%.0lf has %ld sources in exclusion zone", + coords->data.F64[0], coords->data.F64[1], numWithin); + if (indices->n == 1) { + // Only itself inside the exclusion zone + filtered = psArrayAdd(filtered, filtered->n, source); + } else { + float inMag = mag->data.F32[i]; // Input magnitude + bool filter = false; // Filter this source? + for (int j = 0; j < indices->n; j++) { + long index = indices->data.S64[j]; // Index of matching source + float compareMag = mag->data.F32[index]; // Magnitude of matching source + if (fabsf(inMag - compareMag) < minMagDiff) { + filter = true; + break; + } + } + if (!filter) { + filtered = psArrayAdd(filtered, filtered->n, source); + } else { + numFiltered++; + } + } + } + psFree(coords); + psFree(tree); + psFree(x); + psFree(y); + psFree(mag); + + psLogMsg("ppStack", PS_LOG_INFO, "Filtered out %d of %d sources", numFiltered, numGood); + + return filtered; +} + + +psImage *ppStackTarget(ppStackOptions *options, pmConfig *config) +{ + psAssert(options, "Require options"); + psAssert(config, "Require configuration"); + + psAssert(options->psf, "Require target PSF"); + psAssert(options->sourceLists, "Require source lists"); + + psMetadata *recipe = psMetadataLookupMetadata(NULL, config->recipes, PPSTACK_RECIPE); // ppStack recipe + psAssert(recipe, "We've thrown an error on this before."); + psMetadata *ppsub = psMetadataLookupMetadata(NULL, config->recipes, "PPSUB"); // PPSUB recipe + psAssert(recipe, "We've thrown an error on this before."); + + psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of input bits to mask + psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Input bits to mask + float targetFrac = psMetadataLookupF32(NULL, recipe, "TARGET.FRAC"); // Target min flux fraction of noise + + int num = options->num; // Number of inputs + int numCols = 0, numRows = 0; // Size of image + + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS); // Random number generator + float minVariance = INFINITY; // Minimum variance + for (int i = 0; i < num; i++) { + if (options->inputMask->data.U8[i]) { + continue; + } + psTrace("ppStack", 2, "Characterising image %d....\n", i); + pmFPAfileActivate(config->files, false, NULL); + ppStackFileActivationSingle(config, PPSTACK_FILES_TARGET, true, i); + + pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT.VARIANCE", i); // File to read + pmFPAview *view = ppStackFilesIterateDown(config); + if (!view) { + psFree(rng); + return NULL; + } + pmReadout *readout = pmFPAviewThisReadout(view, file->fpa); // Input readout + psFree(view); + + if (numCols == 0 && numRows == 0) { + numCols = readout->image->numCols; + numRows = readout->image->numRows; + } else if (numCols != readout->image->numCols || + numRows != readout->image->numRows) { + psError(PPSTACK_ERR_ARGUMENTS, true, "Sizes of input images don't match: %dx%d vs %dx%d", + readout->image->numCols, readout->image->numRows, numCols, numRows); + psFree(rng); + return NULL; + } + + psImage *variance = readout->variance, *mask = readout->mask; // Dereference images + for (int y = 0; y < numRows; y++) { + for (int x = 0; x < numCols; x++) { + if (!isfinite(variance->data.F32[y][x])) { + mask->data.PS_TYPE_IMAGE_MASK_DATA[y][x] |= maskVal; + } + } + } + + psStats *bg = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Statistics for background + if (!psImageBackground(bg, NULL, variance, mask, maskVal, rng)) { + psError(PPSTACK_ERR_DATA, false, "Can't measure background for image."); + psFree(rng); + return NULL; + } + + float norm = powf(10.0, -0.4 * options->norm->data.F32[i]); // Normalisation from stars + float meanVariance = bg->robustMedian * PS_SQR(norm); // Mean variance in normalised image + + if (meanVariance < minVariance) { + minVariance = meanVariance; + } + + if (!ppStackFilesIterateUp(config)) { + psFree(rng); + return NULL; + } + + ppStackMemDump("target"); + } + psFree(rng); + + float minFlux = targetFrac * sqrtf(minVariance); // Minimum flux for target image + + int footprint = psMetadataLookupS32(NULL, ppsub, "STAMP.FOOTPRINT"); // Stamp half-size + int size = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size + float minMagDiff = psMetadataLookupF32(NULL, recipe, "TARGET.MINMAG"); // Minimum magnitude difference + + // For the sake of stamps, remove nearby sources + psArray *stampSources = stackSourcesFilter(options->sources, footprint, minMagDiff); // Filtered list + + bool oldThreads = pmReadoutFakeThreads(true); // Old threading state + pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF + if (!pmReadoutFakeFromSources(fake, numCols, numRows, stampSources, SOURCE_MASK, NULL, NULL, options->psf, + minFlux, footprint + size, false, true)) { + psError(PPSTACK_ERR_DATA, false, "Unable to generate fake image with target PSF."); + psFree(fake); + return NULL; + } + pmReadoutFakeThreads(oldThreads); + + psFree(stampSources); + + psImage *target = psMemIncrRefCounter(fake->image); + psFree(fake); + + return target; +}