Index: trunk/psModules/src/imcombine/Makefile.am =================================================================== --- trunk/psModules/src/imcombine/Makefile.am (revision 15813) +++ trunk/psModules/src/imcombine/Makefile.am (revision 15814) @@ -14,6 +14,6 @@ pmSubtractionMatch.c \ pmSubtractionParams.c \ - pmSubtractionStamps.c \ - pmPSFEnvelope.c + pmSubtractionStamps.c +# pmPSFEnvelope.c pkginclude_HEADERS = \ @@ -28,6 +28,6 @@ pmSubtractionMatch.h \ pmSubtractionParams.h \ - pmSubtractionStamps.h \ - pmPSFEnvelope.h + pmSubtractionStamps.h +# pmPSFEnvelope.h CLEANFILES = *~ Index: trunk/psModules/src/imcombine/pmPSFEnvelope.c =================================================================== --- trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 15813) +++ trunk/psModules/src/imcombine/pmPSFEnvelope.c (revision 15814) @@ -72,6 +72,6 @@ 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 + 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) @@ -82,15 +82,14 @@ 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 + 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 + 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->peak = pmPeakAlloc(xFake - dx, yFake - dy, PEAK_FLUX, PM_PEAK_LONE); fake->type = PM_SOURCE_TYPE_STAR; @@ -104,17 +103,19 @@ fakes->data[index] = fake; - xOffset->data.S32[index] = 0; // dx; - yOffset->data.S32[index] = 0; // dy; + xOffset->data.S32[index] = dx; + yOffset->data.S32[index] = dy; } } // Generate fake images with each PSF, and take the envelope - psImage *envelope = psImageAlloc(numCols, numRows, PS_TYPE_F32); // Image with envelope of PSFs + psImage *envelope = psImageAlloc(fakeSize, fakeSize, 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)) { + pmResiduals *resid = psf->residuals;// PSF residuals + psf->residuals = NULL; + if (!pmReadoutFakeFromSources(fakeRO, fakeSize, fakeSize, fakes, xOffset, yOffset, psf, + NAN, radius, true)) { psError(PS_ERR_UNKNOWN, false, "Unable to generate fake readout."); psFree(envelope); @@ -122,6 +123,8 @@ psFree(xOffset); psFree(fakes); + psf->residuals = resid; return NULL; } + psf->residuals = resid; // Need to renormalise sources so they all have the same peak. You would think they do have the same @@ -159,4 +162,6 @@ } + + #ifdef TESTING { @@ -183,16 +188,56 @@ #endif - // Reset the fake source pixels to the envelope pixels + // Put the fake sources onto a full-size image 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->image = psImageAlloc(numCols, numRows, PS_TYPE_F32); + psImageInit(readout->image, 0.0); + for (int i = 0; i < numFakes; i++) { + pmSource *source = fakes->data[i]; // Fake source + // Position of source on fake image + int xFake = source->peak->x + xOffset->data.S32[i]; + int yFake = source->peak->y + yOffset->data.S32[i]; + psRegion region = psRegionSet(xFake - radius, xFake + radius, + yFake - radius, yFake + radius); // PSF region + psImage *subImage = psImageSubset(envelope, region); // Subimage of fake PSF + + // Position of source on "real" image + int x0 = source->peak->x - radius; + int y0 = source->peak->y - radius; + + if (!psImageOverlaySection(readout->image, subImage, x0, y0, "=")) { + psError(PS_ERR_UNKNOWN, false, "Unable to overlay PSF"); + psFree(readout); + psFree(xOffset); + psFree(yOffset); + psFree(fakes); + return NULL; + } + } + readout->weight = (psImage*)psBinaryOp(NULL, readout->image, "+", psScalarAlloc(1.0, PS_TYPE_F32)); readout->mask = psImageAlloc(numCols, numRows, PS_TYPE_MASK); psImageInit(readout->mask, 0); +#ifdef TESTING + { + // Write out the envelope + psFits *fits = psFitsOpen("psf_field_full.fits", "w"); + psFitsWriteImage(fits, NULL, readout->image, 0, NULL); + psFitsClose(fits); + } +#endif + 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 + 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 + + psFree(source->pixels); + psFree(source->weight); + psFree(source->maskView); + psFree(source->maskObj); + source->pixels = NULL; + source->weight = NULL; + source->maskView = NULL; + source->maskObj = NULL; if (!pmSourceDefinePixels(source, readout, x, y, (float)radius)) { @@ -205,4 +250,20 @@ } +#if 0 + // Need to do some magic so that pmSourceMoments behaves --- the only offsets it understands are the + // col0,row0 in the psImage. + int oldCol0 = source->pixels->col0; + int oldRow0 = source->pixels->row0; + int col0 = oldCol0 - xOffset->data.S32[i]; + int row0 = oldRow0 - yOffset->data.S32[i]; + + P_PSIMAGE_SET_COL0(source->pixels, col0); + P_PSIMAGE_SET_ROW0(source->pixels, row0); + P_PSIMAGE_SET_COL0(source->weight, col0); + P_PSIMAGE_SET_ROW0(source->weight, row0); + P_PSIMAGE_SET_COL0(source->maskObj, col0); + P_PSIMAGE_SET_ROW0(source->maskObj, row0); +#endif + if (!pmSourceMoments(source, radius)) { psError(PS_ERR_UNKNOWN, false, "Unable to measure moments for source."); @@ -213,5 +274,17 @@ return NULL; } - } + +#if 0 + P_PSIMAGE_SET_COL0(source->pixels, oldCol0); + P_PSIMAGE_SET_ROW0(source->pixels, oldRow0); + P_PSIMAGE_SET_COL0(source->weight, oldCol0); + P_PSIMAGE_SET_ROW0(source->weight, oldRow0); + P_PSIMAGE_SET_COL0(source->maskObj, oldCol0); + P_PSIMAGE_SET_ROW0(source->maskObj, oldRow0); +#endif + } + + psMemCheckCorruption(stderr, true); + pmPSFOptions *options = pmPSFOptionsAlloc(); // Options for fitting a PSF @@ -224,5 +297,5 @@ options->radius = radius; - options->psfTrendMode = PM_TREND_POLY_ORD; + options->psfTrendMode = PM_TREND_MAP; options->psfTrendNx = 1; options->psfTrendNy = 1; @@ -261,5 +334,6 @@ pmReadout *generated = pmReadoutAlloc(NULL); // Generated image - pmReadoutFakeFromSources(generated, numCols, numRows, fakes, xOffset, yOffset, psf, NAN, radius); + pmReadoutFakeFromSources(generated, numCols, numRows, fakes, xOffset, yOffset, psf, NAN, radius, + false); { psFits *fits = psFitsOpen("psf_field_model.fits", "w"); Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 15813) +++ trunk/psModules/src/objects/pmSource.c (revision 15814) @@ -6,6 +6,6 @@ * @author EAM, IfA: significant modifications. * - * @version $Revision: 1.46 $ $Name: not supported by cvs2svn $ - * @date $Date: 2007-12-13 21:45:08 $ + * @version $Revision: 1.47 $ $Name: not supported by cvs2svn $ + * @date $Date: 2007-12-13 22:09:35 $ * * Copyright 2004 Maui High Performance Computing Center, University of Hawaii @@ -703,4 +703,6 @@ if (*vMsk) { vMsk++; + psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to mask: %d\n", + col, row, (int)*vMsk); continue; } @@ -713,4 +715,8 @@ // radius is just a function of (xDiff, yDiff) if (!VALID_RADIUS(xDiff, yDiff, R2)) { +#if 1 + psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to position: %f %f\n", + col, row, xDiff, yDiff); +#endif continue; } @@ -722,4 +728,8 @@ // XXX EAM : should this limit be user-defined? if (PS_SQR(pDiff) < wDiff) { +#if 1 + psTrace("psModules.objects", 10, "Ignoring pixel %d,%d due to insignificance: %f, %f\n", + col, row, pDiff, wDiff); +#endif continue; } @@ -746,5 +756,6 @@ // XXX EAM - the limit is a bit arbitrary. make it user defined? if ((numPixels < 0.75*R2) || (Sum <= 0)) { - psTrace ("psModules.objects", 3, "no valid pixels for source\n"); + psTrace ("psModules.objects", 3, "insufficient valid pixels (%d vs %d; %f) for source\n", + numPixels, (int)(0.75*R2), Sum); psTrace("psModules.objects", 5, "---- end (false) ----\n"); return (false); Index: trunk/psModules/src/psmodules.h =================================================================== --- trunk/psModules/src/psmodules.h (revision 15813) +++ trunk/psModules/src/psmodules.h (revision 15814) @@ -110,4 +110,5 @@ // The following headers are from random locations, here because they cross bounds #include +//#include #endif