Index: /trunk/ppStack/src/ppStack.h =================================================================== --- /trunk/ppStack/src/ppStack.h (revision 26897) +++ /trunk/ppStack/src/ppStack.h (revision 26898) @@ -13,8 +13,9 @@ typedef enum { PPSTACK_MASK_CAL = 0x01, // Photometric calibration failed - PPSTACK_MASK_MATCH = 0x02, // PSF-matching failed - PPSTACK_MASK_CHI2 = 0x04, // Chi^2 too deviant - PPSTACK_MASK_REJECT = 0x08, // Rejection failed - PPSTACK_MASK_BAD = 0x10, // Bad image (too many pixels rejected) + PPSTACK_MASK_PSF = 0x02, // PSF measurement failed + PPSTACK_MASK_MATCH = 0x04, // PSF-matching failed + PPSTACK_MASK_CHI2 = 0x08, // Chi^2 too deviant + PPSTACK_MASK_REJECT = 0x10, // Rejection failed + PPSTACK_MASK_BAD = 0x20, // Bad image (too many pixels rejected) PPSTACK_MASK_ALL = 0xff // All errors } ppStackMask; Index: /trunk/ppStack/src/ppStackCamera.c =================================================================== --- /trunk/ppStack/src/ppStackCamera.c (revision 26897) +++ /trunk/ppStack/src/ppStackCamera.c (revision 26898) @@ -373,6 +373,6 @@ psMetadataLookupBool(&mdok, config->arguments, "-photometry"); // perform photometry if (doPhotom) { - // This file, PSPHOT.INPUT, is just used as a carrier; output files (eg, PSPHOT.RESID) are defined by - // psphotDefineFiles + // This pmFPAfile, PSPHOT.INPUT, is just used as a carrier; output files (eg, + // PSPHOT.RESID) are defined by psphotDefineFiles pmFPAfile *psphotInput = pmFPAfileDefineFromFPA(config, output->fpa, 1, 1, "PSPHOT.INPUT"); if (!psphotInput) { @@ -380,4 +380,6 @@ return false; } + // specify the number of psphot input images + psMetadataAddS32 (config->arguments, PS_LIST_TAIL, "PSPHOT.INPUT.NUM", PS_META_REPLACE, "number of inputs", 1); // Define associated psphot input/output files Index: /trunk/ppStack/src/ppStackConvolve.c =================================================================== --- /trunk/ppStack/src/ppStackConvolve.c (revision 26897) +++ /trunk/ppStack/src/ppStackConvolve.c (revision 26898) @@ -59,9 +59,9 @@ pmFPAfileActivate(config->files, false, NULL); ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, i); - if (options->convolve) { - // XXX PPSTACK.CONV.KERNEL not defined unless convolve - // pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL"); - pmFPAfileActivateSingle(config->files, true, "PPSTACK.CONV.KERNEL", i); // Activated file - } + if (options->convolve) { + // XXX PPSTACK.CONV.KERNEL not defined unless convolve + // pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL"); + pmFPAfileActivateSingle(config->files, true, "PPSTACK.CONV.KERNEL", i); // Activated file + } pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", i); // File of interest @@ -94,4 +94,5 @@ options->origCovars->data[i] = psMemIncrRefCounter(readout->covariance); if (!ppStackMatch(readout, options, i, config)) { + // XXX many things can cause a failure of ppStackMatch -- should some be handled differently? psErrorStackPrint(stderr, "Unable to match image %d --- ignoring.", i); options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] |= PPSTACK_MASK_MATCH; Index: /trunk/ppStack/src/ppStackFinish.c =================================================================== --- /trunk/ppStack/src/ppStackFinish.c (revision 26897) +++ /trunk/ppStack/src/ppStackFinish.c (revision 26898) @@ -83,5 +83,5 @@ fprintf(options->statsFile, "%s", statsMDC); } - psFree((void *)statsMDC); + psFree(statsMDC); fclose(options->statsFile); options->statsFile = NULL; pmConfigRunFilenameAddWrite(config, "STATS", psMetadataLookupStr(NULL, config->arguments, "STATS")); Index: /trunk/ppStack/src/ppStackMatch.c =================================================================== --- /trunk/ppStack/src/ppStackMatch.c (revision 26897) +++ /trunk/ppStack/src/ppStackMatch.c (revision 26898) @@ -144,5 +144,5 @@ psMetadataAddU8(psphotRecipe, PS_LIST_TAIL, "MASK.PSPHOT", PS_META_REPLACE, "user-defined mask", maskBad); - psImage *binned = psphotBackgroundModel(ro, config); // Binned background model + psImage *binned = psphotModelBackgroundReadoutNoFile(ro, config); // Binned background model psImageBinning *binning = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.BACKGROUND.BINNING"); // Binning for model @@ -274,9 +274,11 @@ PM_SUBTRACTION_ANALYSIS_KERNEL); - pmSubtractionAnalysis(readout->analysis, NULL, kernels, region, + pmSubtractionAnalysis(conv->analysis, NULL, kernels, region, readout->image->numCols, readout->image->numRows); psKernel *kernel = pmSubtractionKernel(kernels, 0.0, 0.0, false); // Convolution kernel + bool oldThreads = psImageCovarianceSetThreads(true); // Old thread setting psKernel *covar = psImageCovarianceCalculate(kernel, readout->covariance); // Covariance matrix + psImageCovarianceSetThreads(oldThreads); psFree(readout->covariance); readout->covariance = covar; @@ -298,5 +300,9 @@ float rej = psMetadataLookupF32(NULL, ppsub, "REJ"); // Rejection threshold float kernelError = psMetadataLookupF32(NULL, ppsub, "KERNEL.ERR"); // Relative systematic error in kernel + float normFrac = psMetadataLookupF32(NULL, ppsub, "NORM.FRAC"); // Fraction of window for normalisn windw float sysError = psMetadataLookupF32(NULL, ppsub, "SYS.ERR"); // Relative systematic error in images + float skyErr = psMetadataLookupF32(NULL, ppsub, "SKY.ERR"); // Additional error in sky + float covarFrac = psMetadataLookupF32(NULL, ppsub, "COVAR.FRAC"); // Fraction for covariance calculation + const char *typeStr = psMetadataLookupStr(NULL, ppsub, "KERNEL.TYPE"); // Kernel type pmSubtractionKernelsType type = pmSubtractionKernelsTypeFromString(typeStr); // Kernel type @@ -315,4 +321,16 @@ float poorFrac = psMetadataLookupF32(&mdok, ppsub, "POOR.FRACTION"); // Fraction for "poor" + bool scale = psMetadataLookupBool(NULL, ppsub, "SCALE"); // Scale kernel parameters? + float scaleRef = psMetadataLookupF32(NULL, ppsub, "SCALE.REF"); // Reference for scaling + float scaleMin = psMetadataLookupF32(NULL, ppsub, "SCALE.MIN"); // Minimum for scaling + float scaleMax = psMetadataLookupF32(NULL, ppsub, "SCALE.MAX"); // Maximum for scaling + if (!isfinite(scaleRef) || !isfinite(scaleMin) || !isfinite(scaleMax)) { + psError(PS_ERR_BAD_PARAMETER_VALUE, false, + "Scale parameters (SCALE.REF=%f, SCALE.MIN=%f, SCALE.MAX=%f) not set in PPSUB recipe.", + scaleRef, scaleMin, scaleMax); + return false; + } + + // These values are specified specifically for stacking const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS");// Stamps filename @@ -358,8 +376,10 @@ fake->mask = psImageCopy(NULL, readout->mask, PS_TYPE_IMAGE_MASK); +#if 1 // Add the background into the target image psImage *bgImage = stackBackgroundModel(readout, config); // Image of background psBinaryOp(fake->image, fake->image, "+", bgImage); psFree(bgImage); +#endif #ifdef TESTING @@ -395,5 +415,5 @@ if (kernel) { if (!pmSubtractionMatchPrecalc(NULL, conv, fake, readout, readout->analysis, - stride, kernelError, maskVal, maskBad, maskPoor, + stride, kernelError, covarFrac, maskVal, maskBad, maskPoor, poorFrac, badFrac)) { psError(PS_ERR_UNKNOWN, false, "Unable to convolve images."); @@ -408,10 +428,27 @@ } } else { + // Scale the input parameters + psVector *widthsCopy = psVectorCopy(NULL, widths, PS_TYPE_F32); // Copy of kernel widths + if (scale && !pmSubtractionParamsScale(&size, &footprint, widthsCopy, + options->inputSeeing->data.F32[index], + options->targetSeeing, scaleRef, scaleMin, scaleMax)) { + psError(PS_ERR_UNKNOWN, false, "Unable to scale kernel parameters"); + psFree(fake); + psFree(optWidths); + psFree(stampSources); + psFree(conv); + psFree(widthsCopy); + if (threads > 0) { + pmSubtractionThreadsFinalize(readout, fake); + } + return false; + } + if (!pmSubtractionMatch(NULL, conv, fake, readout, footprint, stride, regionSize, spacing, - threshold, stampSources, stampsName, type, size, order, widths, + threshold, stampSources, stampsName, type, size, order, widthsCopy, orders, inner, ringsOrder, binning, penalty, - optimum, optWidths, optOrder, optThresh, iter, rej, sysError, - kernelError, maskVal, maskBad, maskPoor, poorFrac, badFrac, - PM_SUBTRACTION_MODE_2)) { + optimum, optWidths, optOrder, optThresh, iter, rej, normFrac, + sysError, skyErr, kernelError, covarFrac, maskVal, maskBad, maskPoor, + poorFrac, badFrac, PM_SUBTRACTION_MODE_2)) { psError(PS_ERR_UNKNOWN, false, "Unable to match images."); psFree(fake); @@ -419,4 +456,5 @@ psFree(stampSources); psFree(conv); + psFree(widthsCopy); if (threads > 0) { pmSubtractionThreadsFinalize(readout, fake); @@ -424,4 +462,5 @@ return false; } + psFree(widthsCopy); } @@ -495,10 +534,5 @@ 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); } @@ -526,10 +560,35 @@ } + // Kernel normalisation + { + 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_NORM); + 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_TYPE_F32); + float norm = item->data.F32; // Normalisation + sum += norm; + num++; + } + psFree(iter); + float conv = sum/num; // Mean normalisation from convolution + float stars = powf(10.0, -0.4 * options->norm->data.F32[index]); // Normalisation from stars + float renorm = stars / conv; // Renormalisation to apply + psLogMsg("ppStack", PS_LOG_INFO, "Renormalising image %d by %f (kernel: %f, stars: %f)\n", + index, renorm, conv, stars); + psBinaryOp(readout->image, readout->image, "*", psScalarAlloc(renorm, PS_TYPE_F32)); + psBinaryOp(readout->variance, readout->variance, "*", psScalarAlloc(PS_SQR(renorm), PS_TYPE_F32)); + } + // 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); + psWarning("Maximum deconvolution fraction (%f) exceeds limit (%f) --- rejecting image %d\n", + deconv, deconvLimit, index); psFree(conv); return NULL; Index: /trunk/ppStack/src/ppStackOptions.c =================================================================== --- /trunk/ppStack/src/ppStackOptions.c (revision 26897) +++ /trunk/ppStack/src/ppStackOptions.c (revision 26898) @@ -20,7 +20,9 @@ psFree(options->convVariances); psFree(options->psf); + psFree(options->inputSeeing); psFree(options->inputMask); psFree(options->sourceLists); psFree(options->norm); + psFree(options->sources); psFree(options->cells); psFree(options->kernels); @@ -43,4 +45,5 @@ options->convolve = true; options->matchZPs = true; + options->photometry = false; options->stats = NULL; options->statsFile = NULL; @@ -54,7 +57,10 @@ options->num = 0; options->psf = NULL; + options->inputSeeing = NULL; + options->targetSeeing = NAN; options->inputMask = NULL; options->sourceLists = NULL; options->norm = NULL; + options->sources = NULL; options->cells = NULL; options->kernels = NULL; Index: /trunk/ppStack/src/ppStackOptions.h =================================================================== --- /trunk/ppStack/src/ppStackOptions.h (revision 26897) +++ /trunk/ppStack/src/ppStackOptions.h (revision 26898) @@ -10,4 +10,5 @@ bool convolve; // Convolve images? bool matchZPs; // Adjust relative fluxes based on transparency analysis? + bool photometry; // Perform photometry? psMetadata *stats; // Statistics for output FILE *statsFile; // File to which to write statistics @@ -18,8 +19,11 @@ // Prepare pmPSF *psf; // Target PSF + psVector *inputSeeing; // Input seeing FWHMs + float targetSeeing; // Target seeing FWHM float sumExposure; // Sum of exposure times psVector *inputMask; // Mask for inputs psArray *sourceLists; // Individual lists of sources for matching psVector *norm; // Normalisation for each image + psArray *sources; // Matched sources // Convolve psArray *cells; // Cells for convolved images --- a handle for reading again Index: /trunk/ppStack/src/ppStackPSF.c =================================================================== --- /trunk/ppStack/src/ppStackPSF.c (revision 26897) +++ /trunk/ppStack/src/ppStackPSF.c (revision 26898) @@ -14,4 +14,6 @@ const psArray *psfs, const psVector *inputMask) { + bool mdok = false; + #ifndef TESTING // Get the recipe values @@ -24,4 +26,11 @@ int psfOrder = psMetadataLookupS32(NULL, recipe, "PSF.ORDER"); // Spatial order for PSF + psString maskValStr = psMetadataLookupStr(&mdok, recipe, "MASK.VAL"); // Name of bits to mask going in + if (!mdok || !maskValStr) { + psError(PS_ERR_UNKNOWN, false, "Unable to find MASK.VAL in recipe"); + return false; + } + psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask + for (int i = 0; i < psfs->n; i++) { if (inputMask->data.U8[i]) { @@ -32,6 +41,5 @@ // Solve for the target PSF - pmPSF *psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, - psfOrder, psfOrder); + pmPSF *psf = pmPSFEnvelope(numCols, numRows, psfs, psfInstances, psfRadius, psfModel, psfOrder, psfOrder, maskVal); if (!psf) { psError(PS_ERR_UNKNOWN, false, "Unable to determine output PSF."); Index: /trunk/ppStack/src/ppStackPhotometry.c =================================================================== --- /trunk/ppStack/src/ppStackPhotometry.c (revision 26897) +++ /trunk/ppStack/src/ppStackPhotometry.c (revision 26898) @@ -19,6 +19,5 @@ psAssert(recipe, "We've thrown an error on this before."); - bool mdok; // Status of MD lookup - if (!psMetadataLookupBool(&mdok, recipe, "PHOTOMETRY")) { + if (!options->photometry) { // Nothing to do return true; @@ -60,6 +59,5 @@ "Bits to use for mark", markValue); - pmCell *sourcesCell = pmFPAfileThisCell(config->files, photView, "PPSTACK.OUTPUT"); - psArray *inSources = psMetadataLookupPtr(NULL, sourcesCell->analysis, "PSPHOT.SOURCES"); + psArray *inSources = options->sources; if (!inSources) { psError(PS_ERR_UNKNOWN, false, "Unable to find input sources"); @@ -92,9 +90,12 @@ if (options->stats) { pmReadout *photRO = pmFPAviewThisReadout(photView, photFile->fpa); // Readout with the sources - psArray *sources = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.SOURCES"); // Sources - psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0, - "Number of sources detected", sources ? sources->n : 0); - psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0, - "Time to do photometry", psTimerMark("PPSTACK_PHOT")); + pmDetections *detections = psMetadataLookupPtr(NULL, photRO->analysis, "PSPHOT.DETECTIONS"); // detections + if (detections) { + psAssert (detections->allSources, "missing sources?"); + psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", detections->allSources->n); + } else { + psMetadataAddS32(options->stats, PS_LIST_TAIL, "NUM_SOURCES", 0, "Number of sources detected", 0); + } + psMetadataAddF32(options->stats, PS_LIST_TAIL, "TIME_PHOT", 0, "Time to do photometry", psTimerMark("PPSTACK_PHOT")); } psFree(photView); Index: /trunk/ppStack/src/ppStackPrepare.c =================================================================== --- /trunk/ppStack/src/ppStackPrepare.c (revision 26897) +++ /trunk/ppStack/src/ppStackPrepare.c (revision 26898) @@ -176,16 +176,18 @@ - bool redoPhot = psMetadataLookupBool(NULL, recipe, "PHOT"); - psArray *sources = NULL; - - if (options->convolve || options->matchZPs || redoPhot) { - pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources - 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); - } + bool redoPhot = psMetadataLookupBool(NULL, recipe, "PHOT"); + + pmDetections *detections = NULL; + if (options->convolve || options->matchZPs || options->photometry || redoPhot) { + pmReadout *ro = pmFPAviewThisReadout(view, inputFile->fpa); // Readout with sources + detections = psMetadataLookupPtr(NULL, ro->analysis, "PSPHOT.DETECTIONS"); // Sources + if (!detections) { + psError(PS_ERR_UNEXPECTED_NULL, false, "Unable to find source detections in readout."); + return NULL; + } + psAssert (detections->allSources, "missing sources?"); + + options->sourceLists->data[index] = psMemIncrRefCounter(detections->allSources); + } // Re-do photometry if we don't trust the source lists @@ -194,7 +196,7 @@ pmFPAfileActivate(config->files, false, NULL); ppStackFileActivationSingle(config, PPSTACK_FILES_CONVOLVE, true, index); - if (options->convolve) { - pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL"); - } + if (options->convolve) { + pmFPAfileActivate(config->files, true, "PPSTACK.CONV.KERNEL"); + } pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PPSTACK.INPUT", index); // File pmFPAview *photView = ppStackFilesIterateDown(config); @@ -206,5 +208,5 @@ pmReadout *ro = pmFPAviewThisReadout(view, file->fpa); // Readout of interest - if (!ppStackInputPhotometer(ro, sources, config)) { + if (!ppStackInputPhotometer(ro, detections->allSources, config)) { psError(PS_ERR_UNKNOWN, false, "Unable to do photometry on input sources"); psFree(view); @@ -225,4 +227,47 @@ } psFree(fileIter); + + psString log = psStringCopy("Input seeing FWHMs:\n"); // Log message + bool havePSFs = false; // Do we have any PSFs? + options->inputSeeing = psVectorAlloc(num, PS_TYPE_F32); + psVectorInit(options->inputSeeing, NAN); + for (int i = 0; i < num; i++) { + pmPSF *psf = psfs->data[i]; // PSF for image + if (!psf) { + continue; + } + havePSFs = true; + + if (psf->trendNx > 0 && psf->trendNy > 0) { + double sumFWHM = 0.0; // FWHM for image + int numFWHM = 0; // Number of FWHM measurements + for (int y = 0; y < psf->trendNy; y++) { + float yPos = ((float)y + 0.5) / (float)psf->trendNy * numRows; + for (int x = 0; x < psf->trendNx; x++) { + float xPos = ((float)x + 0.5) / (float)psf->trendNx * numCols; + float fwhm = pmPSFtoFWHM(psf, xPos, yPos); // FWHM for image + if (isfinite(fwhm)) { + sumFWHM += fwhm; + numFWHM++; + } + } + } + if (numFWHM == 0) { + options->inputSeeing->data.F32[i] = NAN; + options->inputMask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PPSTACK_MASK_PSF; + psLogMsg("ppStack", PS_LOG_INFO, "Unable to measure PSF FWHM for image %d --- rejected.", i); + } else { + options->inputSeeing->data.F32[i] = sumFWHM / (float)numFWHM; + } + } else { + options->inputSeeing->data.F32[i] = pmPSFtoFWHM(psf, 0.5 * numCols, 0.5 * numRows); + } + + psStringAppend(&log, "Input %d: %f\n", i, options->inputSeeing->data.F32[i]); + } + if (havePSFs) { + psLogMsg("ppStack", PS_LOG_INFO, "%s", log); + } + psFree(log); // Generate target PSF @@ -237,4 +282,6 @@ psMetadataAddPtr(config->arguments, PS_LIST_TAIL, "PSF.TARGET", PS_DATA_UNKNOWN, "Target PSF for stack", options->psf); + options->targetSeeing = pmPSFtoFWHM(options->psf, 0.5 * numCols, 0.5 * numRows); // FWHM for target + psLogMsg("ppStack", PS_LOG_INFO, "Target seeing FWHM: %f\n", options->targetSeeing); pmChip *outChip = pmFPAfileThisChip(config->files, view, "PPSTACK.OUTPUT"); // Output chip Index: /trunk/ppStack/src/ppStackReadout.c =================================================================== --- /trunk/ppStack/src/ppStackReadout.c (revision 26897) +++ /trunk/ppStack/src/ppStackReadout.c (revision 26898) @@ -126,5 +126,5 @@ int kernelSize = psMetadataLookupS32(NULL, ppsub, "KERNEL.SIZE"); // Kernel half-size - psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in + psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits @@ -219,5 +219,5 @@ bool safe = psMetadataLookupBool(&mdok, recipe, "SAFE"); // Be safe when combining small numbers of pixels - psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.IN"); // Name of bits to mask going in + psString maskValStr = psMetadataLookupStr(NULL, recipe, "MASK.VAL"); // Name of bits to mask going in psImageMaskType maskVal = pmConfigMaskGet(maskValStr, config); // Bits to mask going in to pmSubtractionMatch psString maskSuspectStr = psMetadataLookupStr(NULL, recipe, "MASK.SUSPECT"); // Name of suspect mask bits Index: /trunk/ppStack/src/ppStackSetup.c =================================================================== --- /trunk/ppStack/src/ppStackSetup.c (revision 26897) +++ /trunk/ppStack/src/ppStackSetup.c (revision 26898) @@ -22,5 +22,8 @@ psAssert(recipe, "We've thrown an error on this before."); - options->matchZPs = psMetadataLookupBool(NULL, recipe, "MATCH.ZERO.POINTS"); // Adjust zero points based on tranparency analysis? + // XXX : switch to this name? options->matchZPs = psMetadataLookupBool(NULL, recipe, "MATCH.ZERO.POINTS"); // Adjust zero points based on tranparency analysis? + options->matchZPs = psMetadataLookupBool(NULL, recipe, "ZP"); // Adjust zero points? + + options->photometry = psMetadataLookupBool(NULL, recipe, "PHOTOMETRY"); // Perform photometry? options->convolve = psMetadataLookupBool(NULL, recipe, "CONVOLVE"); // Convolve images? Index: /trunk/ppStack/src/ppStackSources.c =================================================================== --- /trunk/ppStack/src/ppStackSources.c (revision 26897) +++ /trunk/ppStack/src/ppStackSources.c (revision 26898) @@ -61,13 +61,13 @@ PS_ASSERT_PTR_NON_NULL(config, false); - if (!options->matchZPs) { - int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs - options->norm = psVectorAlloc(num, PS_TYPE_F32); - psVectorInit (options->norm, 0.0); - - // XXX do I need to set this? - // options->sumExposure = sumExpTime; - - return true; + if (!options->matchZPs && !options->photometry) { + int num = psMetadataLookupS32(NULL, config->arguments, "INPUTS.NUM"); // Number of inputs + options->norm = psVectorAlloc(num, PS_TYPE_F32); + psVectorInit (options->norm, 0.0); + + // XXX do I need to set this? + // options->sumExposure = sumExpTime; + + return true; } @@ -191,11 +191,68 @@ #endif - psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1, - iter2, starRej2, starSys2, starLimit, - transIter, transRej, transThresh); // Transparencies for each image - if (!trans) { - psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies"); - return false; - } + if (options->matchZPs) { + psVector *trans = pmSourceMatchRelphot(matches, zp, tol, iter1, starRej1, starSys1, + iter2, starRej2, starSys2, starLimit, + transIter, transRej, transThresh); // Transparencies per image + if (!trans) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure transparencies"); + return false; + } + for (int i = 0; i < trans->n; i++) { + if (!isfinite(trans->data.F32[i])) { + inputMask->data.U8[i] |= PPSTACK_MASK_CAL; + } + } + // M = m + c0 + c1 * airmass - 2.5log(t) + transparency + // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0 + // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1 + // 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 + if (options->matchZPs) { + options->norm = psVectorAlloc(num, PS_TYPE_F32); + for (int i = 0; i < num; i++) { + if (!isfinite(trans->data.F32[i])) { + continue; + } + 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); + + for (int j = 0; j < sources->n; j++) { + pmSource *source = sources->data[j]; // Source of interest + if (!source) { + continue; + } + source->psfMag += magCorr; + } + } + } + +#ifdef TESTING + dumpMatches("source_mags.dat", num, matches, zp, trans); +#endif + psFree(trans); + +#ifdef TESTING + // Double check: all transparencies should be zero + { + psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches + if (!matches) { + psError(PS_ERR_UNKNOWN, false, "Unable to match sources"); + psFree(zp); + return false; + } + psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej, + transThresh, starRej, starSys); + for (int i = 0; i < num; i++) { + fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]); + } + psFree(trans); + } +#endif + } + #if 0 @@ -216,84 +273,28 @@ #endif -#ifdef TESTING - dumpMatches("source_mags.dat", num, matches, zp, trans); -#endif - - for (int i = 0; i < trans->n; i++) { - if (!isfinite(trans->data.F32[i])) { - inputMask->data.U8[i] |= PPSTACK_MASK_CAL; - } - } - - // Save best matches SOMEWHERE for future photometry - // XXX this is a really poor output location; clean up the pmFPAfiles used in ppStack - pmCell *sourcesCell = pmFPAfileThisCell(config->files, view, "PPSTACK.OUTPUT"); - psArray *sourcesBest = psArrayAllocEmpty(matches->n); - - // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible - int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required - for (int i = 0; i < matches->n; i++) { - pmSourceMatch *match = matches->data[i]; // Match of interest - if (match->num < minMatches) { - continue; - } - - // We need to grab a single instance of this source: just take the first available - int image = match->image->data.S32[0]; // Index of image - int index = match->index->data.S32[0]; // Index of source within image - psArray *sources = sourceLists->data[image]; // Sources for image - pmSource *source = sources->data[index]; // Source of interest - - psArrayAdd(sourcesBest, sourcesBest->n, source); - } - psMetadataAdd(sourcesCell->analysis, PS_LIST_TAIL, "PSPHOT.SOURCES", PS_DATA_ARRAY | PS_META_REPLACE, - "psphot sources", sourcesBest); - psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n); - psFree(sourcesBest); + if (options->photometry) { + // Save best matches for future photometry + psArray *sourcesBest = options->sources = psArrayAllocEmpty(matches->n); + + // XXX something of a hack: require at least 2 detections or the nominated fraction of the max possible + int minMatches = PS_MAX(2, fracMatch * num);// Minimum number of matches required + for (int i = 0; i < matches->n; i++) { + pmSourceMatch *match = matches->data[i]; // Match of interest + if (match->num < minMatches) { + continue; + } + + // We need to grab a single instance of this source: just take the first available + int image = match->image->data.S32[0]; // Index of image + int index = match->index->data.S32[0]; // Index of source within image + psArray *sources = sourceLists->data[image]; // Sources for image + pmSource *source = sources->data[index]; // Source of interest + + psArrayAdd(sourcesBest, sourcesBest->n, source); + } + psLogMsg("ppStack", PS_LOG_INFO, "Selected %ld sources for photometry analysis", sourcesBest->n); + } psFree(matches); - - // M = m + c0 + c1 * airmass - 2.5log(t) + transparency - // Want sources to have m corresponding to airmass = 1 and t = sumExpTime and transparency = 0 - // m_0 + c1 * airmass_0 + 2.5log(t_0) - trans_0 = m_1 + c1 * airmass_1 + 2.5log(t_1) - trans_1 - // 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])) { - continue; - } - 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); - - for (int j = 0; j < sources->n; j++) { - pmSource *source = sources->data[j]; // Source of interest - if (!source) { - continue; - } - source->psfMag += magCorr; - } - } - psFree(trans); - -#ifdef TESTING - // Double check: all transparencies should be zero - { - psArray *matches = pmSourceMatchSources(sourceLists, radius); // List of matches - if (!matches) { - psError(PS_ERR_UNKNOWN, false, "Unable to match sources"); - psFree(zp); - return false; - } - psVector *trans = pmSourceMatchRelphot(matches, zp, iter, tol, starLimit, transIter, transRej, - transThresh, starRej, starSys); - for (int i = 0; i < num; i++) { - fprintf(stderr, "Transparency of image %d: %f\n", i, trans->data.F32[i]); - } - psFree(trans); - } -#endif return true;