Index: /trunk/ppSub/src/ppSubArguments.c =================================================================== --- /trunk/ppSub/src/ppSubArguments.c (revision 18288) +++ /trunk/ppSub/src/ppSubArguments.c (revision 18289) @@ -186,4 +186,5 @@ psMetadataAddS32(arguments, PS_LIST_TAIL, "-order", 0, "Spatial polynomial order", -1); psMetadataAddStr(arguments, PS_LIST_TAIL, "-type", 0, "Kernel type (ISIS|POIS|SPAM|FRIES|GUNK|RINGS)", NULL); + psMetadataAddF32(arguments, PS_LIST_TAIL, "-penalty", 0, "Penalty for wideness", NAN); psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-widths", 0, "ISIS Gaussian FWHMs (comma-separated)", NULL); psMetadataAddStr(arguments, PS_LIST_TAIL, "-isis-orders", 0, "ISIS polynomial orders (comma-separated)", NULL); @@ -260,4 +261,5 @@ VALUE_ARG_RECIPE_FLOAT("-rej", "REJ", F32); VALUE_ARG_RECIPE_FLOAT("-badfrac", "BADFRAC", F32); + VALUE_ARG_RECIPE_FLOAT("-penalty", "PENALTY", F32); valueArgRecipeStr(arguments, recipe, "-mask-bad", "MASK.BAD", config->arguments); Index: /trunk/ppSub/src/ppSubCamera.c =================================================================== --- /trunk/ppSub/src/ppSubCamera.c (revision 18288) +++ /trunk/ppSub/src/ppSubCamera.c (revision 18289) @@ -82,5 +82,5 @@ // Output image - pmFPAfile *output = pmFPAfileDefineSkycell(config, NULL, "PPSUB.OUTPUT"); + pmFPAfile *output = pmFPAfileDefineFromFile(config, input, 1, 1, "PPSUB.OUTPUT"); if (!output) { psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT")); @@ -94,5 +94,5 @@ // Output mask - pmFPAfile *outMask = pmFPAfileDefineSkycell(config, output->fpa, "PPSUB.OUTPUT.MASK"); + pmFPAfile *outMask = pmFPAfileDefineOutput(config, output->fpa, "PPSUB.OUTPUT.MASK"); if (!outMask) { psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.MASK")); @@ -107,5 +107,5 @@ // Output weight if (inputWeight && refWeight) { - pmFPAfile *outWeight = pmFPAfileDefineSkycell(config, output->fpa, "PPSUB.OUTPUT.WEIGHT"); + pmFPAfile *outWeight = pmFPAfileDefineOutput(config, output->fpa, "PPSUB.OUTPUT.WEIGHT"); if (!outWeight) { psError(PS_ERR_IO, false, _("Unable to generate output file from PPSUB.OUTPUT.WEIGHT")); @@ -119,8 +119,10 @@ } +#if 0 if (!pmFPAAddSourceFromFormat(output->fpa, "Subtraction", output->format)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate output FPA."); return false; } +#endif pmFPAfile *sources = pmFPAfileDefineFromArgs(&status, config, "PPSUB.SOURCES", "PPSUB.SOURCES"); @@ -153,5 +155,5 @@ // Internal-ish file for getting the PSF from the matched addition - pmFPAfile *psf = pmFPAfileDefineSkycell(config, psphot->fpa, "PSPHOT.PSF.LOAD"); + pmFPAfile *psf = pmFPAfileDefineOutputFromFile(config, output, "PSPHOT.PSF.LOAD"); if (!psf) { psError(PS_ERR_IO, false, "Failed to build FPA from PSPHOT.PSF.LOAD"); Index: /trunk/ppSub/src/ppSubReadout.c =================================================================== --- /trunk/ppSub/src/ppSubReadout.c (revision 18288) +++ /trunk/ppSub/src/ppSubReadout.c (revision 18289) @@ -12,5 +12,5 @@ #define WCS_TOLERANCE 0.001 // Tolerance for WCS -//#define TESTING // For test output +#define TESTING // For test output @@ -73,4 +73,5 @@ int ringsOrder = psMetadataLookupS32(NULL, config->arguments, "RINGS.ORDER"); // RINGS polynomial order int binning = psMetadataLookupS32(NULL, config->arguments, "SPAM.BINNING"); // Binning for SPAM kernel + float penalty = psMetadataLookupF32(NULL, config->arguments, "PENALTY"); // Penalty for wideness psMaskType maskBad = pmConfigMask(psMetadataLookupStr(NULL, config->arguments, "MASK.BAD"), config); // Value to mask @@ -131,6 +132,6 @@ if (!pmSubtractionMatch(inConv, refConv, inRO, refRO, footprint, regionSize, spacing, threshold, sources, - stampsName, type, size, order, widths, orders, inner, ringsOrder, - binning, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad, + stampsName, type, size, order, widths, orders, inner, ringsOrder, binning, + penalty, optimum, optWidths, optOrder, optThresh, iter, rej, maskBad, maskBlank, badFrac, mode)) { psError(PS_ERR_UNKNOWN, false, "Unable to match images."); @@ -184,10 +185,31 @@ #endif - // "Subtract" the mask and weight map + // Subtraction is: minuend - subtrahend + pmReadout *minuend = inConv; + pmReadout *subtrahend = refConv; + + if (reverse) { + pmReadout *temp = subtrahend; + subtrahend = minuend; + minuend = temp; + } + +#ifdef TESTING + { + pmReadoutMaskApply(minuend, maskBlank); + psFits *fits = psFitsOpen("minuend.fits", "w"); + psFitsWriteImage(fits, NULL, minuend->image, 0, NULL); + psFitsClose(fits); + } + { + pmReadoutMaskApply(subtrahend, maskBlank); + psFits *fits = psFitsOpen("subtrahend.fits", "w"); + psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL); + psFitsClose(fits); + } +#endif + outRO->mask = (psImage*)psBinaryOp(outRO->mask, inConv->mask, "|", refConv->mask); - if (inConv->weight && refConv->weight) { - outRO->weight = (psImage*)psBinaryOp(outRO->weight, inConv->weight, "+", refConv->weight); - } - outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true; + outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true; // It'll be there soon // Photometry is to be performed in two stages: @@ -198,6 +220,8 @@ pmPSF *psf = NULL; // PSF for photometry if (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) { - // We use a summed image as the basis for the PSF: this will have the maximum S/N. - outRO->image = (psImage*)psBinaryOp(outRO->image, inConv->image, "+", refConv->image); + outRO->image = psImageCopy(outRO->image, minuend->image, PS_TYPE_F32); + if (minuend->weight) { + outRO->weight = psImageCopy(outRO->weight, minuend->weight, PS_TYPE_F32); + } pmFPAfile *photFile = psMetadataLookupPtr(NULL, config->files, "PSPHOT.INPUT"); @@ -244,33 +268,12 @@ } - // Do the subtraction - { - // Subtraction is: minuend - subtrahend - pmReadout *minuend = inConv; - pmReadout *subtrahend = refConv; - - if (reverse) { - pmReadout *temp = subtrahend; - subtrahend = minuend; - minuend = temp; - } - - outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image); - -#ifdef TESTING - { - psFits *fits = psFitsOpen("minuend.fits", "w"); - psFitsWriteImage(fits, NULL, minuend->image, 0, NULL); - psFitsClose(fits); - } - { - psFits *fits = psFitsOpen("subtrahend.fits", "w"); - psFitsWriteImage(fits, NULL, subtrahend->image, 0, NULL); - psFitsClose(fits); - } -#endif - - pmReadoutMaskApply(outRO, maskBlank); - } + // Do the actual subtraction + outRO->image = (psImage*)psBinaryOp(outRO->image, minuend->image, "-", subtrahend->image); + if (inConv->weight && refConv->weight) { + outRO->weight = (psImage*)psBinaryOp(outRO->weight, inConv->weight, "+", refConv->weight); + } + outRO->data_exists = outCell->data_exists = outCell->parent->data_exists = true; + + pmReadoutMaskApply(outRO, maskBlank); psFree(inConv); @@ -304,6 +307,4 @@ } - psTraceSetLevel("psphot", 6); - // Photometry stage 2: find and measure sources on the subtracted image if (psMetadataLookupBool(NULL, config->arguments, "PHOTOMETRY")) { @@ -314,5 +315,5 @@ pmReadout *psfRO = pmFPAfileThisReadout(config->files, view, "PSPHOT.PSF.LOAD"); if (!psfRO) { - psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find file PSPHOT.PSF.LOAD"); + psError(PS_ERR_UNEXPECTED_NULL, true, "Unable to find readout for file PSPHOT.PSF.LOAD"); return false; } Index: /trunk/ppSub/test/fake.c =================================================================== --- /trunk/ppSub/test/fake.c (revision 18289) +++ /trunk/ppSub/test/fake.c (revision 18289) @@ -0,0 +1,182 @@ +#include +#include + +// Quick and dirty program to generate cross-directed PSFs for testing ppSub. +// To build: gcc fake.c -o fake -g -std=gnu99 `pslib-config --cflags --libs` + + +// PSF for images +#define X_AXIS_1 2.0 // Length of x axis (FWHM) for image 1 +#define Y_AXIS_1 3.0 // Length of y axis (FWHM) for image 1 +#define X_AXIS_2 3.0 // Length of x axis (FWHM) for image 1 +#define Y_AXIS_2 2.0 // Length of y axis (FWHM) for image 1 + +// Image parameters +#define NUMCOLS 1024 // Number of columns +#define NUMROWS 1024 // Number of rows +#define IMAGE_SCALE 0.2 // Arcsec per pixel + +// Source parameters +#define SOURCES_LUM 0.4 // Luminosity function +#define SOURCES_MAG 14.5 // Magnitude for density +#define SOURCES_DENSITY 30000.0 // Source density at nominated magnitude (per deg^2) +#define SOURCES_ZP 25.0 // Magnitude ZP +#define SOURCES_GAUSSIAN 0 // Gaussian sources (boolean)? + +// Noise properties +#define NOISE_1 10.0 // Noise in image 1 +#define NOISE_2 10.0 // Noise in image 2 +#define NOISE_FRAC 0.1 // Go out to this fraction of the noise + +// Conversion factor by which to multiply sigma to get FWHM. Roughly 2.35 +#define SIGMA_FWHM (2.0*sqrt((double)2.0*(log((double)2.0)))) + +// Add a star to an image +static void addstar(float x, float y, // Coordinates of source + float flux, // Total (integrated) flux of source + float xSigma, float ySigma, // Gaussian width of source + float noise, // Lowest noise level to consider + psImage *image // Image to which to add source + ) +{ + float peak = flux / 2.0 / M_PI / sqrtf(PS_SQR(xSigma) + PS_SQR(ySigma)); // Peak flux + if (peak < noise) { + return; + } + +#if SOURCES_GAUSSIAN + // Gaussian + float radius = sqrtf(2.0 * (PS_SQR(xSigma) + PS_SQR(ySigma)) * logf(peak / noise)); // Radius for source +#else + // Basic Waussian + float radius; + if (peak / noise > 10.0) { + // Need to treat the bright limit differently because the wings are much more noticeable + radius = sqrtf(2.0 * (PS_SQR(xSigma) + PS_SQR(ySigma)) * peak / noise); + } else { + radius = sqrtf(2.0) * hypot(1.0 / xSigma, 1.0 / ySigma) * sqrtf(2.0 * logf(peak / noise)); + } +#endif + + // Bounds of source + int xMin = PS_MAX(floorf(x - radius), 0); + int xMax = PS_MIN(ceilf(x + radius), NUMCOLS - 1); + int yMin = PS_MAX(floorf(y - radius), 0); + int yMax = PS_MIN(ceilf(y + radius), NUMROWS - 1); + + for (int j = yMin; j <= yMax; j++) { + float dy = j - y; + for (int i = xMin; i <= xMax; i++) { + float dx = i - x; +#if SOURCES_GAUSSIAN + // Gaussian + float value = peak * expf(- PS_SQR(dx) / 2.0 / PS_SQR(xSigma) - + PS_SQR(dy) / 2.0 / PS_SQR(ySigma)); +#else + // Basic Waussian + float z = 0.5 * (PS_SQR(dx / xSigma) + PS_SQR(dy / ySigma)); + float value = peak / (1.0 + z + PS_SQR(z)); +#endif + + image->data.F32[j][i] += value; + } + } + + return; +} + + +int main(int argc, char *argv[]) +{ + psRandom *rng = psRandomAlloc(PS_RANDOM_TAUS, 0); // Random number generator + + // Images to generate + psImage *image1 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32); + psImage *image2 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32); + + psImageInit(image1, 0.0); + psImageInit(image2, 0.0); + + float magNorm = SOURCES_DENSITY * NUMCOLS * NUMROWS * PS_SQR(IMAGE_SCALE / 3600.0) / + powf(10.0, (SOURCES_LUM * SOURCES_MAG)); // Normalisation for magnitudes + + float faintMag = -2.5 * log10(PS_MIN(NOISE_1 * sqrtf(X_AXIS_1 * Y_AXIS_1), + NOISE_2 * sqrtf(X_AXIS_2 * Y_AXIS_2)) * + NOISE_FRAC) + SOURCES_ZP; // Faintest magnitude + printf("Faintest mag: %f\n", faintMag); + + long numStars = magNorm * powf(10.0, SOURCES_LUM * faintMag); // Number of stars + printf("%ld stars\n", numStars); + + // Gaussian widths, in sigma + float xSigma1 = X_AXIS_1 / SIGMA_FWHM; + float ySigma1 = Y_AXIS_1 / SIGMA_FWHM; + float xSigma2 = X_AXIS_2 / SIGMA_FWHM; + float ySigma2 = Y_AXIS_2 / SIGMA_FWHM; + + for (int i = 0; i < numStars; i++) { + float x = psRandomUniform(rng) * NUMCOLS; // x coordinate + float y = psRandomUniform(rng) * NUMROWS; // y coordinate + float mag = log10((i + 1.0) / magNorm) / SOURCES_LUM; // Magnitude of source + float flux = powf(10.0, -0.4 * (mag - SOURCES_ZP)); // Total flux of source + + addstar(x, y, flux, xSigma1, ySigma1, NOISE_1 * NOISE_FRAC, image1); + addstar(x, y, flux, xSigma2, ySigma2, NOISE_2 * NOISE_FRAC, image2); + } + +#if 1 + psImage *var1 = (psImage*)psBinaryOp(NULL, image1, "+", psScalarAlloc(PS_SQR(NOISE_1), PS_TYPE_F32)); + psImage *var2 = (psImage*)psBinaryOp(NULL, image2, "+", psScalarAlloc(PS_SQR(NOISE_2), PS_TYPE_F32)); +#else + psImage *var1 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32); + psImage *var2 = psImageAlloc(NUMCOLS, NUMROWS, PS_TYPE_F32); + psImageInit(var1, PS_SQR(NOISE_1)); + psImageInit(var2, PS_SQR(NOISE_2)); +#endif + + for (int y = 0; y < NUMROWS; y++) { + for (int x = 0; x < NUMCOLS; x++) { + image1->data.F32[y][x] += psRandomGaussian(rng) * sqrtf(var1->data.F32[y][x]); + image2->data.F32[y][x] += psRandomGaussian(rng) * sqrtf(var2->data.F32[y][x]); + } + } +#if 1 + var1 = (psImage*)psBinaryOp(var1, image1, "+", psScalarAlloc(PS_SQR(NOISE_1), PS_TYPE_F32)); + var2 = (psImage*)psBinaryOp(var2, image2, "+", psScalarAlloc(PS_SQR(NOISE_2), PS_TYPE_F32)); +#else + psImageInit(var1, PS_SQR(NOISE_1)); + psImageInit(var2, PS_SQR(NOISE_2)); +#endif + + { + psFits *fits = psFitsOpen("image1.fits", "w"); + psFitsWriteImage(fits, NULL, image1, 0, NULL); + psFitsClose(fits); + } + { + psFits *fits = psFitsOpen("image2.fits", "w"); + psFitsWriteImage(fits, NULL, image2, 0, NULL); + psFitsClose(fits); + } + + { + psFits *fits = psFitsOpen("variance1.fits", "w"); + psFitsWriteImage(fits, NULL, var1, 0, NULL); + psFitsClose(fits); + } + { + psFits *fits = psFitsOpen("variance2.fits", "w"); + psFitsWriteImage(fits, NULL, var2, 0, NULL); + psFitsClose(fits); + } + + psFree(image1); + psFree(image2); + + psFree(var1); + psFree(var2); + + psFree(rng); + + return PS_EXIT_SUCCESS; +}