Index: trunk/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 15762) +++ trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 15762) @@ -0,0 +1,284 @@ +#ifdef HAVE_CONFIG_H +#include +#endif + + +#include +#include +#include + +#include "pmHDU.h" +#include "pmFPA.h" +#include "pmReadoutFake.h" + +#include "pmMoments.h" +#include "pmResiduals.h" +#include "pmGrowthCurve.h" +#include "pmTrend2D.h" +#include "pmPSF.h" +#include "pmModel.h" +#include "pmModelClass.h" +#include "pmModelUtils.h" +#include "pmPeaks.h" +#include "pmSource.h" +#include "pmSourceUtils.h" +#include "pmSourceFitModel.h" +#include "pmPSFtry.h" + + +#include "pmPSFEnvelope.h" + + + +#define TESTING // Enable test output +#define PEAK_FLUX 1.0e3 // peak flux for each source +#define WEIGHT_VAL 1.0e0 // Weighting for image +#define PSF_STATS PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV // Statistics options for measuring PSF + + + +// XXX To do: +// +// * PSF representation doesn't normalise (x,y) positions in polynomials: need to find out how I can have a +// star at (x1,y1) in the image that is "really" supposed to be at (x2,y2). +// +// * PSF variation when only a portion of the image is present (e.g., the edge of an FPA overlapping a +// skycell) may mean a disastrously weird PSF in the missing regions. To counter this, get a region of +// validity for each PSF (perhaps from an associated mask, or have the user work it out), and taper the PSF +// when outside this region (perhaps multiply the peak flux by a Gaussian whose arguments are the distance +// from the valid region, and a width that the user supplies). + + +// We deliberately do not include the calculation of and storing of residuals (data - model) for the PSF +// model, because (1) there is no code in psModules to do this, and we're not going to implement it here; and +// (2) this is intended to generate "nice" or "ideal" PSFs to feed into pmSubtraction (PSF matching code), so +// any residuals will hopefully be dealt with by that. + + +pmPSF *pmPSFEnvelope(int numCols, int numRows, // Size of original image + const psArray *inputs, // Input PSF models + int instances, // Number of instances per dimension + int radius, // Radius of each PSF + const char *modelName // Name of PSF model to use + ) +{ + PS_ASSERT_INT_POSITIVE(numCols, NULL); + PS_ASSERT_INT_POSITIVE(numRows, NULL); + PS_ASSERT_ARRAY_NON_NULL(inputs, NULL); + PS_ASSERT_INT_POSITIVE(instances, NULL); + PS_ASSERT_INT_POSITIVE(radius, NULL); + PS_ASSERT_STRING_NON_EMPTY(modelName, NULL); + + float xOrigSpacing = (float)(numCols - 2 * radius) / (float)(instances - 1); // Spacing between instances + float yOrigSpacing = (float)(numRows - 2 * radius) / (float)(instances - 1); // Spacing between instances + // int fakeSpacing = 2 * radius + 1; // Spacing between instances (x and y) in the fake image + // int fakeSize = instances * fakeSpacing; // Size of fake image + + // Generate list of fake sources (instances of the PSF) + int numFakes = PS_SQR(instances); // Number of fake sources + psArray *fakes = psArrayAlloc(numFakes); // Fake sources + psVector *xOffset = psVectorAlloc(numFakes, PS_TYPE_S32); // X offset from fake position to image + psVector *yOffset = psVectorAlloc(numFakes, PS_TYPE_S32); // Y offset from fake position to image + for (int j = 0, index = 0; j < instances; j++) { + float yOrig = j * yOrigSpacing + radius; // Source position in original image + // float yFake = j * fakeSpacing + radius; // Position in fake image + // int dy = yFake - yOrig; // Difference between fake and original position + + for (int i = 0; i < instances; i++, index++) { + float xOrig = i * xOrigSpacing + radius; // Source position in original image + // float xFake = i * fakeSpacing + radius; // Position in fake image + // int dx = xFake - xOrig; // Difference between fake and original position + + pmSource *fake = pmSourceAlloc(); // Fake source + // fake->peak = pmPeakAlloc(xFake - dx, yFake - dy, PEAK_FLUX, PM_PEAK_LONE); + fake->peak = pmPeakAlloc(xOrig, yOrig, PEAK_FLUX, PM_PEAK_LONE); + fake->type = PM_SOURCE_TYPE_STAR; + +#ifdef TESTING + // Required to generate model image + fake->psfMag = -2.5 * log10(PEAK_FLUX); +#endif + + psTrace("psModules.imcombine", 5, "Source %d: %.2f,%.2f\n", + index, xOrig, yOrig); + + fakes->data[index] = fake; + xOffset->data.S32[index] = 0; // dx; + yOffset->data.S32[index] = 0; // dy; + } + } + + // Generate fake images with each PSF, and take the envelope + psImage *envelope = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Image with envelope of PSFs + psImageInit(envelope, 0.0); + pmReadout *fakeRO = pmReadoutAlloc(NULL); // Fake readout + for (int i = 0; i < inputs->n; i++) { + pmPSF *psf = inputs->data[i]; // PSF of interest + if (!pmReadoutFakeFromSources(fakeRO, numCols, numRows, fakes, xOffset, yOffset, psf, + NAN, radius)) { + psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout."); + psFree(envelope); + psFree(yOffset); + psFree(xOffset); + psFree(fakes); + return NULL; + } + + // Need to renormalise sources so they all have the same peak. You would think they do have the same + // peak already, but it seems that the residual map messes things up by adding extra flux + for (int j = 0; j < numFakes; j++) { + pmSource *source = fakes->data[j]; // Fake source + float x = source->peak->xf + xOffset->data.S32[j]; // x coordinate of source + float y = source->peak->yf + yOffset->data.S32[j]; // y coordinate of source + +#if 0 + psImageInterpolateOptions *interp = psImageInterpolateOptionsAlloc(PS_INTERPOLATE_BILINEAR, + fakeRO->image, NULL, NULL, 0, + 0.0, 0.0, 1, 0, 0.0); + double flux; // Flux of peak + if (!psImageInterpolate(&flux, NULL, NULL, x, y, interp)) { + psError(PS_ERR_UNKNOWN, false, "Unable to interpolate source image."); + psFree(envelope); + psFree(yOffset); + psFree(xOffset); + psFree(fakes); + return NULL; + } +#else + double flux = fakeRO->image->data.F32[(int)y][(int)x]; +#endif + + float norm = PEAK_FLUX / flux; // Normalisation for source + psRegion region = psRegionSet(x - radius, x + radius, y - radius, y + radius); // PSF region + psImage *subImage = psImageSubset(fakeRO->image, region); // Subimage of fake PSF + psImage *subEnv = psImageSubset(envelope, region); // Subimage of envelope + psBinaryOp(subImage, subImage, "*", psScalarAlloc(norm, PS_TYPE_F32)); + psBinaryOp(subEnv, subEnv, "MAX", subImage); + psFree(subImage); + psFree(subEnv); + } + +#ifdef TESTING + { + // Write out the PSF field + psString name = NULL; + psStringAppend(&name, "psf_field_%03d.fits", i); + psFits *fits = psFitsOpen(name, "w"); + psFitsWriteImage(fits, NULL, fakeRO->image, 0, NULL); + psFitsClose(fits); + psFree(name); + } +#endif + + } + psFree(fakeRO); + +#ifdef TESTING + { + // Write out the envelope + psFits *fits = psFitsOpen("psf_field_envelope.fits", "w"); + psFitsWriteImage(fits, NULL, envelope, 0, NULL); + psFitsClose(fits); + } +#endif + + // Reset the fake source pixels to the envelope pixels + pmReadout *readout = pmReadoutAlloc(NULL); // Readout to contain envelope pixels + readout->image = envelope; + readout->weight = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(readout->weight, WEIGHT_VAL); + readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); + psImageInit(readout->mask, 0); + + for (int i = 0; i < numFakes; i++) { + pmSource *source = fakes->data[i]; // Fake source + float x = source->peak->xf + xOffset->data.S32[i]; // x coordinates of source + float y = source->peak->yf + yOffset->data.S32[i]; // y coordinates of source + + if (!pmSourceDefinePixels(source, readout, x, y, (float)radius)) { + psError(PS_ERR_UNKNOWN, false, "Unable to define pixels for source."); + psFree(readout); + psFree(yOffset); + psFree(xOffset); + psFree(fakes); + return NULL; + } + + if (!pmSourceMoments(source, radius, xOffset->data.S32[i], yOffset->data.S32[i])) { + psError(PS_ERR_UNKNOWN, false, "Unable to measure moments for source."); + psFree(readout); + psFree(yOffset); + psFree(xOffset); + psFree(fakes); + return NULL; + } + } + + pmPSFOptions *options = pmPSFOptionsAlloc(); // Options for fitting a PSF + // Don't assume Poisson errors --- our fake system isn't Poisson. + options->poissonErrorsPhotLMM = false; + options->poissonErrorsPhotLin = false; + options->poissonErrorsParams = false; + options->stats = psStatsAlloc(PSF_STATS); + + options->radius = radius; + + options->psfTrendMode = PM_TREND_POLY_ORD; + options->psfTrendNx = 1; + options->psfTrendNy = 1; + + options->psfFieldNx = numCols; + options->psfFieldNy = numRows; + options->psfFieldXo = 0; + options->psfFieldYo = 0; + + pmSourceFitModelInit (15, 0.01, WEIGHT_VAL, options->poissonErrorsPhotLMM); + + pmPSFtry *try = pmPSFtryModel(fakes, modelName, options, 0, 0xff); + psFree(options); + if (!try) { + psError(PS_ERR_UNKNOWN, false, "Unable to fit PSF model to PSF envelope."); + psFree(readout); + psFree(yOffset); + psFree(xOffset); + psFree(fakes); + return NULL; + } + + pmPSF *psf = psMemIncrRefCounter(try->psf); // Output PSF + psFree(try); + +#ifdef TESTING + { + // Need to translate peak flux --> integrated flux + pmModel *fakeModel = pmModelFromPSFforXY(psf, (float)numCols / 2.0, (float)numRows / 2.0, + 1.0); // Fake model, with central intensity of 1.0 + float flux0 = fakeModel->modelFlux(fakeModel->params); // Flux for central intensity of 1.0 + for (int i = 0; i < numFakes; i++) { + pmSource *source = fakes->data[i]; // Fake source + source->psfMag -= 2.5 * log10(flux0); + } + + pmReadout *generated = pmReadoutAlloc(NULL); // Generated image + pmReadoutFakeFromSources(generated, numCols, numRows, fakes, xOffset, yOffset, psf, NAN, radius); + { + psFits *fits = psFitsOpen("psf_field_model.fits", "w"); + psFitsWriteImage(fits, NULL, generated->image, 0, NULL); + psFitsClose(fits); + } + psBinaryOp(generated->image, generated->image, "-", envelope); + { + psFits *fits = psFitsOpen("psf_field_resid.fits", "w"); + psFitsWriteImage(fits, NULL, generated->image, 0, NULL); + psFitsClose(fits); + } + psFree(generated); + } +#endif + + psFree(yOffset); + psFree(xOffset); + psFree(fakes); + psFree(readout); + return psf; +} Index: trunk/psModules/src/imcombine/pmPSFEnvelope.h =================================================================== --- trunk/psModules/src/imcombine/pmPSFEnvelope.h (revision 15762) +++ trunk/psModules/src/imcombine/pmPSFEnvelope.h (revision 15762) @@ -0,0 +1,21 @@ +#ifndef PM_PSF_ENVELOPE_H +#define PM_PSF_ENVELOPE_H + +#include +#include +#include +#include +#include +#include + +/// Generate a PSF which is the envelope of an array of PSFs +/// +/// Generates multiple instances of the PSFs (distributed over an image) +pmPSF *pmPSFEnvelope(int numCols, int numRows, // Size of original image + const psArray *inputs, // Input PSF models + int instances, // Number of instances per dimension + int radius, // Radius of each PSF + const char *modelName // Name of PSF model to use + ); + +#endif Index: trunk/psModules/src/imcombine/pmSubtractionEquation.c =================================================================== --- trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 15761) +++ trunk/psModules/src/imcombine/pmSubtractionEquation.c (revision 15762) @@ -704,4 +704,5 @@ double diff = NAN; // Difference between iterations psVector *vector = NULL; // RHS for equations + int iter = 0; // Iteration number do { if (!vector) { @@ -765,9 +766,16 @@ // Get difference produced by iteration diff = 0.0; +#if 0 psVector *sol2 = kernels->solution2; // Solution of second equation for (int i = 0; i < numParams2; i++) { - diff += fabs(sol2->data.F64[i] - lastSol2->data.F64[i]) / sol2->data.F64[i]; + diff += fabs(sol2->data.F64[i] - lastSol2->data.F64[i]); } lastSol2 = psVectorCopy(lastSol2, sol2, PS_TYPE_F64); + psTrace("psModules.imcombine", 4, "Solution iteration %d: %lf\n", iter, diff); +#endif + + psVectorInit(kernels->solution2, 0.0); + + iter++; } while (diff > tolerance); @@ -877,8 +885,4 @@ } -#ifdef TESTING - psImageInit(residual->image, 0.0); -#endif - for (int j = 0; j < numKernels; j++) { psKernel *conv1 = convolutions1->data[j]; // Convolution of first image @@ -890,5 +894,5 @@ for (int y = - footprint; y <= footprint; y++) { for (int x = - footprint; x <= footprint; x++) { - residual->kernel[y][x] += conv1->kernel[y][x] * coefficient1 - + residual->kernel[y][x] += conv1->kernel[y][x] * coefficient1 + conv2->kernel[y][x] * coefficient2; }