Index: /trunk/ppStack/src/ppStackCamera.c =================================================================== --- /trunk/ppStack/src/ppStackCamera.c (revision 17425) +++ /trunk/ppStack/src/ppStackCamera.c (revision 17426) @@ -72,4 +72,5 @@ { bool haveWeights = false; // Do we have weight maps? + bool havePSF = false; // Do we have PSFs? psMetadata *inputs = psMetadataLookupMetadata(NULL, config->arguments, "INPUTS"); // The inputs info @@ -169,23 +170,30 @@ // Add the psf file if (!psf || strlen(psf) == 0) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF %d", i); - return false; + if (havePSFs) { + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find PSF %d", i); + return false; + } } else { - psArray *psfFiles = psArrayAlloc(1); // Array of filenames for this FPA - psfFiles->data[0] = psMemIncrRefCounter(psf); - psMetadataAddArray(config->arguments, PS_LIST_TAIL, "PSF.FILENAMES", PS_META_REPLACE, - "Filenames of PSF files", psfFiles); - psFree(psfFiles); - - bool status; - pmFPAfile *psfFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPSTACK.INPUT.PSF", - "PSF.FILENAMES"); - if (!status) { - psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf); - return false; - } - if (psfFile->type != PM_FPA_FILE_PSF) { - psError(PS_ERR_IO, true, "PPSTACK.INPUT.PSF is not of type PSF"); - return false; + if (!havePSFs && i != 0) { + psWarning("PSF not provided for all inputs --- ignoring."); + } else { + psArray *psfFiles = psArrayAlloc(1); // Array of filenames for this FPA + psfFiles->data[0] = psMemIncrRefCounter(psf); + psMetadataAddArray(config->arguments, PS_LIST_TAIL, "PSF.FILENAMES", PS_META_REPLACE, + "Filenames of PSF files", psfFiles); + psFree(psfFiles); + + bool status; + pmFPAfile *psfFile = pmFPAfileBindFromArgs(&status, imageFile, config, "PPSTACK.INPUT.PSF", + "PSF.FILENAMES"); + if (!status) { + psError(PS_ERR_UNKNOWN, false, "Unable to define file from psf %d (%s)", i, psf); + return false; + } + if (psfFile->type != PM_FPA_FILE_PSF) { + psError(PS_ERR_IO, true, "PPSTACK.INPUT.PSF is not of type PSF"); + return false; + } + havePSFs = true; } } @@ -233,4 +241,5 @@ psMetadataAddS32(config->arguments, PS_LIST_TAIL, "INPUTS.NUM", 0, "Number of input files", i); + psMetadataAddBool(config->arguments, PS_LIST_TAIL, "HAVE.PSF", 0, "Have PSFs available?", havePSFs); // Output image Index: /trunk/ppStack/src/ppStackLoop.c =================================================================== --- /trunk/ppStack/src/ppStackLoop.c (revision 17425) +++ /trunk/ppStack/src/ppStackLoop.c (revision 17426) @@ -215,5 +215,5 @@ pmReadout *sources = NULL; // Readout with sources to use for PSF matching pmPSF *targetPSF = NULL; // Target PSF - { + if (psMetadataLookupBool(NULL, config->arguments, "HAVE.PSF")) { pmFPAfileActivate(config->files, false, NULL); fileActivation(config, prepareFiles, true); Index: /trunk/ppStack/src/ppStackMatch.c =================================================================== --- /trunk/ppStack/src/ppStackMatch.c (revision 17425) +++ /trunk/ppStack/src/ppStackMatch.c (revision 17426) @@ -13,5 +13,4 @@ //#define TESTING -//#define NO_CONVOLUTION bool ppStackMatch(pmReadout *readout, psArray **regions, psArray **kernels, @@ -22,5 +21,4 @@ assert(kernels && !*kernels); assert(sourcesRO); - assert(psf); assert(config); @@ -37,100 +35,99 @@ int renormWidth = psMetadataLookupS32(&mdok, config->arguments, "RENORM.WIDTH"); // Width for renormalisation box -#ifndef NO_CONVOLUTION - int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order - float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs - float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing - int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size - float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps - int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations - float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold - pmSubtractionKernelsType type = - pmSubtractionKernelsTypeFromString(psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE")); // Kernel type - psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths - psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders - int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius - int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order - int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel - psMaskType maskBlank = pmConfigMask(psMetadataLookupStr(NULL, recipe, "MASK.BLANK"), - config); // Mask for blank reg. - float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction - bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters? - float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search - float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search - float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search - float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search - int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search - - // These values are specified specifically for stacking - const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps - - psVector *optWidths = NULL; // Vector with FWHMs for optimum search - if (optimum) { - optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32); - } - - psArray *sources = NULL; // Sources in image - if (sourcesRO) { - sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES"); - } - - // Generate a fake image to match to - float maxMag = -INFINITY; // Maximum magnitude of sources - for (int i = 0; i < sources->n; i++) { - pmSource *source = sources->data[i]; // Source of interest - if (source->psfMag > maxMag && source->psfMag != MAG_IGNORE) { - maxMag = source->psfMag; - } - } - - pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF - - if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources, NULL, NULL, - psf, powf(10.0, -0.4 * maxMag), 0, false)) { - psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF."); + if (psMetadataLookupBool(&mdok, config->arguments, "HAVE.PSF")) { + assert(psf); + int order = psMetadataLookupS32(NULL, recipe, "SPATIAL.ORDER"); // Spatial polynomial order + float regionSize = psMetadataLookupF32(NULL, recipe, "REGION.SIZE"); // Size of iso-kernel regs + float spacing = psMetadataLookupF32(NULL, recipe, "STAMP.SPACING"); // Typical stamp spacing + int footprint = psMetadataLookupS32(NULL, recipe, "STAMP.FOOTPRINT"); // Stamp half-size + float threshold = psMetadataLookupF32(NULL, recipe, "STAMP.THRESHOLD"); // Threshold for stmps + int iter = psMetadataLookupS32(NULL, recipe, "ITER"); // Rejection iterations + float rej = psMetadataLookupF32(NULL, recipe, "REJ"); // Rejection threshold + pmSubtractionKernelsType type = + pmSubtractionKernelsTypeFromString(psMetadataLookupStr(NULL, recipe, "KERNEL.TYPE")); // Kernel type + psVector *widths = psMetadataLookupPtr(NULL, recipe, "ISIS.WIDTHS"); // ISIS Gaussian widths + psVector *orders = psMetadataLookupPtr(NULL, recipe, "ISIS.ORDERS"); // ISIS Polynomial orders + int inner = psMetadataLookupS32(NULL, recipe, "INNER"); // Inner radius + int ringsOrder = psMetadataLookupS32(NULL, recipe, "RINGS.ORDER"); // RINGS polynomial order + int binning = psMetadataLookupS32(NULL, recipe, "SPAM.BINNING"); // Binning for SPAM kernel + psMaskType maskBlank = pmConfigMask(psMetadataLookupStr(NULL, recipe, "MASK.BLANK"), + config); // Mask for blank reg. + float badFrac = psMetadataLookupF32(NULL, recipe, "BADFRAC"); // Maximum bad fraction + bool optimum = psMetadataLookupBool(&mdok, recipe, "OPTIMUM"); // Derive optimum parameters? + float optMin = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MIN"); // Minimum width for search + float optMax = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.MAX"); // Maximum width for search + float optStep = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.STEP"); // Step for search + float optThresh = psMetadataLookupF32(&mdok, recipe, "OPTIMUM.TOL"); // Tolerance for search + int optOrder = psMetadataLookupS32(&mdok, recipe, "OPTIMUM.ORDER"); // Order for search + + // These values are specified specifically for stacking + const char *stampsName = psMetadataLookupStr(NULL, config->arguments, "STAMPS"); // Filename for stamps + + psVector *optWidths = NULL; // Vector with FWHMs for optimum search + if (optimum) { + optWidths = psVectorCreate(optWidths, optMin, optMax, optStep, PS_TYPE_F32); + } + + psArray *sources = NULL; // Sources in image + if (sourcesRO) { + sources = psMetadataLookupPtr(&mdok, sourcesRO->analysis, "PSPHOT.SOURCES"); + } + + // Generate a fake image to match to + float maxMag = -INFINITY; // Maximum magnitude of sources + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; // Source of interest + if (source->psfMag > maxMag && source->psfMag != MAG_IGNORE) { + maxMag = source->psfMag; + } + } + + pmReadout *fake = pmReadoutAlloc(NULL); // Fake readout with target PSF + + if (!pmReadoutFakeFromSources(fake, readout->image->numCols, readout->image->numRows, sources, NULL, + NULL, psf, powf(10.0, -0.4 * maxMag), 0, false)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate fake image with target PSF."); + psFree(fake); + psFree(optWidths); + psFree(output); + return false; + } + +#ifdef TESTING + { + psFits *fits = psFitsOpen("fake.fits", "w"); + psFitsWriteImage(fits, NULL, fake->image, 0, NULL); + psFitsClose(fits); + } +#endif + + // Do the image matching + if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, regionSize, spacing, threshold, + sources, stampsName, type, size, order, widths, orders, inner, ringsOrder, + binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad, + maskBlank, badFrac, PM_SUBTRACTION_MODE_1)) { + psError(PS_ERR_UNKNOWN, false, "Unable to match images."); + psFree(fake); + psFree(optWidths); + psFree(output); + return false; + } psFree(fake); psFree(optWidths); - psFree(output); - return false; - } - -#ifdef TESTING - { - psFits *fits = psFitsOpen("fake.fits", "w"); - psFitsWriteImage(fits, NULL, fake->image, 0, NULL); - psFitsClose(fits); - } -#endif - - // Do the image matching - if (!pmSubtractionMatch(output, NULL, readout, fake, footprint, regionSize, spacing, threshold, - sources, stampsName, type, size, order, widths, orders, inner, ringsOrder, - binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad, - maskBlank, badFrac, PM_SUBTRACTION_MODE_1)) { - psError(PS_ERR_UNKNOWN, false, "Unable to match images."); - psFree(fake); - psFree(optWidths); - psFree(output); - return false; - } - psFree(fake); - psFree(optWidths); - - // Replace original images with convolved - psFree(readout->image); - psFree(readout->mask); - psFree(readout->weight); - readout->image = psMemIncrRefCounter(output->image); - readout->mask = psMemIncrRefCounter(output->mask); - readout->weight = psMemIncrRefCounter(output->weight); - -#else // NO_CONVOLUTION - // Fake the convolution - { + + // Replace original images with convolved + psFree(readout->image); + psFree(readout->mask); + psFree(readout->weight); + readout->image = psMemIncrRefCounter(output->image); + readout->mask = psMemIncrRefCounter(output->mask); + readout->weight = psMemIncrRefCounter(output->weight); + } 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(size, 0, PM_SUBTRACTION_MODE_1); + pmSubtractionKernels *kernels = pmSubtractionKernelsPOIS(0, 0, PM_SUBTRACTION_MODE_1); // Set solution to delta function kernels->solution1 = psVectorAlloc(kernels->num + 2, PS_TYPE_F64); @@ -142,5 +139,4 @@ psFree(kernels); } -#endif // Extract the regions and solutions used in the image matching