Index: trunk/psphot/src/Makefile.am =================================================================== --- trunk/psphot/src/Makefile.am (revision 42841) +++ trunk/psphot/src/Makefile.am (revision 42842) @@ -247,4 +247,5 @@ psphotSetNFrames.c \ psphotSourceMemory.c \ + psphotFindFeatures.c \ psphotChipParams.c \ psphotGalaxyParams.c Index: trunk/psphot/src/psphot.h =================================================================== --- trunk/psphot/src/psphot.h (revision 42841) +++ trunk/psphot/src/psphot.h (revision 42842) @@ -592,4 +592,6 @@ extern bool psphotINpsphotStack; +bool psphotFindFeatures (pmReadout *readout, psArray *sources); + #endif Index: trunk/psphot/src/psphotArguments.c =================================================================== --- trunk/psphot/src/psphotArguments.c (revision 42841) +++ trunk/psphot/src/psphotArguments.c (revision 42842) @@ -97,4 +97,6 @@ } +void psphotSaveConfig (pmConfig *config); + pmConfig *psphotArguments(int argc, char **argv) { @@ -244,4 +246,7 @@ psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]); + // save a copy for deep retrieval + psphotSaveConfig (config); + psTrace("psphot", 1, "Done with psphotArguments...\n"); return (config); Index: trunk/psphot/src/psphotCheckStarDistribution.c =================================================================== --- trunk/psphot/src/psphotCheckStarDistribution.c (revision 42841) +++ trunk/psphot/src/psphotCheckStarDistribution.c (revision 42842) @@ -89,4 +89,8 @@ if (source->type != PM_SOURCE_TYPE_STAR) continue; + // XXX for a test, do not include faint sources + // EAM XXX : make this optional and default to all + if (source->moments->SN < 20) continue; + float x = source->peak->xf; float y = source->peak->yf; Index: trunk/psphot/src/psphotChoosePSF.c =================================================================== --- trunk/psphot/src/psphotChoosePSF.c (revision 42841) +++ trunk/psphot/src/psphotChoosePSF.c (revision 42842) @@ -114,4 +114,10 @@ options->psfTrendNy = psMetadataLookupS32 (&status, recipe, "PSF.TREND.NY"); assert (status); + + // fit the PSF star positions during PSF model analysis (default false) + options->fitPSFstarCoords = psMetadataLookupBool (&status, recipe, "PSF_FIT_STAR_COORDS"); + if (!status) { + options->fitPSFstarCoords = false; + } // get the fixed PSF fit radius Index: trunk/psphot/src/psphotDefineFiles.c =================================================================== --- trunk/psphot/src/psphotDefineFiles.c (revision 42841) +++ trunk/psphot/src/psphotDefineFiles.c (revision 42842) @@ -1,3 +1,4 @@ # include "psphotInternal.h" +pmConfig *pmConfigMakeTemp (pmConfig *config); bool psphotINpsphotStack = false; @@ -143,6 +144,15 @@ } + // Use a temporary config for the source files, keeping the main user, site, + // system, files, arguments entries. This allows the sources to be from a + // different camera than the input image. NOTE: there is no check that these two + // images match in terms of size, pixel scale, etc. That is up to the user. + // The temporary config structure has PSCAMERA entries removed from the SKYCELL rules + // to allow forced photometry between cameras. + + pmConfig *srcconfig = pmConfigMakeTemp(config); + if (psMetadataLookupPtr(NULL, config->arguments, "SRC")) { - if (!pmFPAfileDefineFromArgs (&status, config, "PSPHOT.INPUT.CMF", "SRC")) { + if (!pmFPAfileDefineFromArgs (&status, srcconfig, "PSPHOT.INPUT.CMF", "SRC")) { psError(PSPHOT_ERR_CONFIG, false, "Failed to find/build PSPHOT.INPUT.CMF"); return status; @@ -151,5 +161,5 @@ if (psMetadataLookupPtr(NULL, config->arguments, "FORCE")) { - if (!pmFPAfileDefineFromArgs (&status, config, "PSPHOT.INPUT.CFF", "FORCE")) { + if (!pmFPAfileDefineFromArgs (&status, srcconfig, "PSPHOT.INPUT.CFF", "FORCE")) { psError(PSPHOT_ERR_CONFIG, false, "Failed to find/build PSPHOT.INPUT.CFF"); return status; @@ -160,5 +170,5 @@ // XXX cannot use pmFPAfileDefineFromArgs: this is explicitly a FITS-based I/O function // supply the attach the - if (!psphotLoadSRCTEXT(input->fpa, config)) { + if (!psphotLoadSRCTEXT(input->fpa, srcconfig)) { psError(PSPHOT_ERR_CONFIG, false, "Failed to load PSPHOT.INPUT.TEXT"); return status; @@ -167,5 +177,5 @@ if (psMetadataLookupPtr(NULL, config->arguments, "PSPHOT.PSF")) { - pmFPAfileBindFromArgs(&status, input, config, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); + pmFPAfileBindFromArgs(&status, input, srcconfig, "PSPHOT.PSF.LOAD", "PSPHOT.PSF"); if (!status) { psError(PSPHOT_ERR_CONFIG, false, "Failed to find/build PSPHOT.PSF.LOAD"); @@ -173,4 +183,6 @@ } } + + psFree (srcconfig); // XXX add in example PSF image thumbnails @@ -214,2 +226,67 @@ return rule; } + +pmConfig *pmConfigMakeTemp (pmConfig *config) { + pmConfig *altconfig = pmConfigAlloc(); + + // these are NULL on pmConfigAlloc + altconfig->user = psMemIncrRefCounter(config->user); // inherit from primary camera + altconfig->site = psMemIncrRefCounter(config->site); // inherit from primary camera + + psFree (altconfig->system); + altconfig->system = psMetadataCopy(NULL, config->system); // container for camera-specific recipe values (to be dropped) + + psFree (altconfig->files); + altconfig->files = psMemIncrRefCounter(config->files); // inherit from primary camera + + psFree (altconfig->arguments); + altconfig->arguments = psMemIncrRefCounter(config->arguments); // inherit from primary camera + + psFree (altconfig->recipes); + altconfig->recipes = psMetadataCopy(NULL, config->recipes); // container for camera-specific recipe values (to be dropped) + + // remove PSCAMERA entries from *-SKYCELL cameras. This allows skycell images from + // one camera to match those of another camera. XXX Make this optional to force + // matching cameras if desired? + bool mdok = false; + psMetadata *cameras = psMetadataLookupMetadata(&mdok, altconfig->system, "CAMERAS"); + psAssert(cameras, "missing cameras in system config info"); + + // iterate over the cameras and find ones with names like _*-SKYCELLS + // psMetadataIterator *camerasIter = psMetadataIteratorAlloc(cameras, PS_LIST_HEAD, "^_.+-SKYCELL$"); + psMetadataIterator *camerasIter = psMetadataIteratorAlloc(cameras, PS_LIST_HEAD, NULL); + psAssert(camerasIter, "unable to generate iterator?"); + + psMetadataItem *camerasItem = NULL; // Item from the metadata + while ((camerasItem = psMetadataGetAndIncrement(camerasIter))) { + psAssert(camerasItem->type == PS_DATA_METADATA, "camerasItem has invalid type"); + + char *name = camerasItem->name; + psAssert (name, "NULL name for item"); + + int nameLen = strlen(name); + if (name[0] != '_') continue; + if (nameLen <= 9) continue; + char *p = &name[nameLen - 8]; + if (strcmp(p, "-SKYCELL")) continue; + + // remove PSCAMERA entries from *-SKYCELL cameras + psMetadata *formats = psMetadataLookupMetadata(&mdok, camerasItem->data.md, "FORMATS"); // List of formats + psAssert (formats, "missing FORMATS in camera.config "); + + psMetadata *format = psMetadataLookupMetadata(&mdok, formats, "SKYCELL"); // one true format for skycells + psAssert (format, "missing SKYCELL metaformat for -SKYCELL entry "); + + psMetadata *rule = psMetadataLookupMetadata(&mdok, format, "RULE"); // one true format for skycells + psAssert (rule, "missing RULE in SKYCELL metaformat for -SKYCELL entry "); + + psString pscamera = psMetadataLookupStr(&mdok, rule, "PSCAMERA"); // one true format for skycells + if (!pscamera) continue; // already removed, or never supplied + + psMetadataRemoveKey (rule, "PSCAMERA"); // allow any camera skycell to match + } + psFree (camerasIter); + + return (altconfig); +} + Index: trunk/psphot/src/psphotFindFeatures.c =================================================================== --- trunk/psphot/src/psphotFindFeatures.c (revision 42842) +++ trunk/psphot/src/psphotFindFeatures.c (revision 42842) @@ -0,0 +1,436 @@ +# include "psphotInternal.h" + +/* Use Hough transform to detect strong linear features (e.g., satellite trails) in the image. + Sources which are plausibly associated with the feature are marked in the MODE2 bits. These + sources are excluded from the analysis of the moments and selection of the PSF stars. + Weaknesses: + - the angular samples and the length bins are hard-wired at 1 degree x 10 pixels. + - the minimum significance in the Hough image and the minimum number of sources associated + with the feature are hard-wired. + - the distance of a source from the feature is only calculated in 1D. The code is not careful + about the center of the linear feature compared to the edges of the bin. + - the features detected in the Hough transform are used to select sources which are then + fitted to a line. The line is used to select sources again and re-fit. Only two steps + are rigidly defined without a clear iteration end criterion. + + */ + +# define DCOS(THETA) cos(THETA*PS_RAD_DEG) +# define DSIN(THETA) sin(THETA*PS_RAD_DEG) +# define HOUGH_MIN_SOURCES 10 +# define HOUGH_MIN_SIGMA 5.0 + +typedef struct { + psImage *fArea; + psImage *image; + psImage *sigma; + int NXpix; + int NYpix; + float Rdel; + float Rmin; + float Rmax; + int Nt; + int Nr; + int nKeep; +} psphotFindFeatures_HoughInfo; + +bool psphotFindFeatures_HoughInit (psphotFindFeatures_HoughInfo *myHough, int NXpix, int NYpix, float Rdel); +bool psphotFindFeatures_HoughImage (psphotFindFeatures_HoughInfo *myHough, psArray *sources); +bool psphotFindFeatures_HoughSigma (psphotFindFeatures_HoughInfo *myHough); +float psphotFindFeatures_LineLength (float radius, float theta, int *edges, int NXpix, int NYpix); +int psphotFindFeatures_PeaksOnTrail (psArray *sources, float theta, float radius, float Rdel); +bool psphotFindFeatures_HoughFree (psphotFindFeatures_HoughInfo *myHough); + +bool psphotSavePeaks (char *filename, psArray *sources) { + + FILE *ftest = fopen (filename, "w"); + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + pmPeak *peak = source->peak; + bool onTrail = source->mode2 & PM_SOURCE_MODE2_ON_LINE; + fprintf (ftest, "%d %d %f %d\n", peak->x, peak->y, peak->detValue, onTrail); + } + fclose (ftest); + return true; +} + +bool psphotSavePoints (char *filename, psVector *x, psVector *y) { + + FILE *ftest = fopen (filename, "w"); + for (int i = 0; i < x->n; i++) { + fprintf (ftest, "%f %f\n", x->data.F32[i], y->data.F32[i]); + } + fclose (ftest); + return true; +} + +// examine the x,y distribution of the sources and attempt to find sources which are due to +// specific finds of features. For example, a bunch of detections in a line could be a +// satellite and should not be used for PSF modeling + +bool psphotFindFeatures (pmReadout *readout, psArray *sources) { + + psphotFindFeatures_HoughInfo myHough; + + int NXpix = readout->image->numCols; + int NYpix = readout->image->numRows; + + // calculate the Hough transform for the source positions + psphotFindFeatures_HoughInit (&myHough, NXpix, NYpix, 10.0); + + // generate the hough transform for the collection of sources + // and find only a single trail at a time. remove that trail + // and try again. + + int nTrails = 0; + int nPeaksOnTrails = 0; + bool foundTrail = true; // true to ensure we try a first pass + while (foundTrail) { + psphotFindFeatures_HoughImage (&myHough, sources); + psphotFindFeatures_HoughSigma (&myHough); + + // select the peaks, sort by the peak amplitude + psArray *lines = pmPeaksInImage (myHough.sigma, HOUGH_MIN_SIGMA); + psArraySort (lines, pmPeaksSortByRawFluxDescend); + + // if we do not find a trail in this pass, we stop + foundTrail = false; + + psImage *image = myHough.image; + psImage *sigma = myHough.sigma; + + // choose the strongest line (> HOUGH_MIN_SIGMA) with at least HOUGH_MIN_SOURCES + // and mark sources associated with this line + for (int i = 0; (i < lines->n) && !foundTrail; i++) { + pmPeak *line = lines->data[i]; + + float hough = image->data.F32[line->y][line->x]; + if (hough < HOUGH_MIN_SOURCES) continue; + + float nsigma = sigma->data.F32[line->y][line->x]; + + float theta = line->x; + float radius = line->y*myHough.Rdel + myHough.Rmin; + psLogMsg ("psphot", PS_LOG_INFO, "trail %d : %f %f (%f : %f)\n", i, theta, radius, hough, nsigma); + + int nPeaksOnThisTrail = psphotFindFeatures_PeaksOnTrail (sources, theta, radius, myHough.Rdel); + + if (nPeaksOnThisTrail < HOUGH_MIN_SOURCES) continue; + + // this is only a valid trail if it has enough sources + nTrails ++; + foundTrail = true; + nPeaksOnTrails += nPeaksOnThisTrail; + } + psFree (lines); + } + psphotFindFeatures_HoughFree (&myHough); + psLogMsg ("psphot", PS_LOG_INFO, "found %d peaks in %d trails\n", nPeaksOnTrails, nTrails); + + return true; +} + +bool psphotFindFeatures_HoughFree (psphotFindFeatures_HoughInfo *myHough) { + + psFree (myHough->fArea); + psFree (myHough->image); + psFree (myHough->sigma); + return true; +} + +// define the basic parameters for all Hough transform images +// the transform image is not populated here, but the fractional area image is +bool psphotFindFeatures_HoughInit (psphotFindFeatures_HoughInfo *myHough, int NXpix, int NYpix, float Rdel) { + + // diagonal length (absolute max value of radius) + int NRpix = sqrt(PS_SQR(NXpix) + PS_SQR(NYpix)); + + // radial bins + myHough->Rmin = -1*NRpix; + myHough->Rmax = +1*NRpix; + myHough->Rdel = Rdel; + + myHough->Nt = 180; // angle range from 0 to 180 in steps of 1 degree + myHough->Nr = (myHough->Rmax - myHough->Rmin) / myHough->Rdel; + + myHough->NXpix = NXpix; + myHough->NYpix = NYpix; + + myHough->fArea = psImageAlloc (myHough->Nt, myHough->Nr, PS_TYPE_F32); + myHough->image = psImageAlloc (myHough->Nt, myHough->Nr, PS_TYPE_F32); + myHough->sigma = psImageAlloc (myHough->Nt, myHough->Nr, PS_TYPE_F32); + + float ImageArea = myHough->NXpix * myHough->NYpix; + for (int theta = 0; theta < myHough->Nt; theta ++) { + for (int rbin = 0; rbin < myHough->Nr; rbin ++) { + + // for the given value of Radius, Theta, what is the line length from edge to edge + int edges[4]; + float radius = rbin * myHough->Rdel + myHough->Rmin; + float LineLength = psphotFindFeatures_LineLength (radius, theta, edges, myHough->NXpix, myHough->NYpix); + float LineArea = LineLength * myHough->Rdel; + myHough->fArea->data.F32[rbin][theta] = LineArea / ImageArea; + } + } + return true; +} + +bool psphotFindFeatures_HoughImage (psphotFindFeatures_HoughInfo *myHough, psArray *sources) { + + if (!sources) return false; + + psImage *image = myHough->image; + psImageInit (image, 0.0); + + int nKeep = 0; + + // generate the Hough transform image + for (int i = 0; i < sources->n; i++) { + pmSource *source = sources->data[i]; + pmPeak *peak = source->peak; + + // skip bad peaks and record the number actually used + if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue; + nKeep ++; + + for (int theta = 0; theta < myHough->Nt; theta ++) { + + float r = peak->x * DCOS(theta) + peak->y * DSIN(theta); + int rbin = (r - myHough->Rmin) / myHough->Rdel; + + if (theta < 0) continue; + if (rbin < 0) continue; + + if (theta >= myHough->Nt) continue; + if (rbin >= myHough->Nr) continue; + + image->data.F32[rbin][theta] ++; + } + } + myHough->nKeep = nKeep; + return true; +} + +bool psphotFindFeatures_HoughSigma (psphotFindFeatures_HoughInfo *myHough) { + + // Given the fractional area image and nKeep sources + // calculate the expected mean number of sources for each line at + // radius R and angle theta. use that to determine significance + + psImage *sigma = myHough->sigma; + psImage *image = myHough->image; + psImage *fArea = myHough->fArea; + + for (int theta = 0; theta < myHough->Nt; theta ++) { + for (int rbin = 0; rbin < myHough->Nr; rbin ++) { + float mean = myHough->nKeep * fArea->data.F32[rbin][theta]; + float meanMin = PS_MAX(mean, 1); + float nsigma = (image->data.F32[rbin][theta] - meanMin) / sqrt(meanMin); + sigma->data.F32[rbin][theta] = nsigma; + } + } + return true; +} + +// mark all sources within myHough->Rdel pix of this line +// I should probably calculate the distance from the line in 2D perpendicular direction +// NOTE: I am marking the sources on this trail, but if it does not have enough +// matches, it is not counted as a real peak. + +int psphotFindFeatures_PeaksOnTrail (psArray *sources, float theta, float radius, float Rdel) { + + float cosTheta = DCOS(theta); + float sinTheta = DSIN(theta); + + // find sources in the line and fit to 1D polynomial + psVector *xpos = psVectorAllocEmpty (100, PS_TYPE_F32); + psVector *ypos = psVectorAllocEmpty (100, PS_TYPE_F32); + psPolynomial1D *poly = psPolynomial1DAlloc (PS_POLYNOMIAL_ORD, 1); + + // if the angle is shallow, use image (x,y), otherwise use image (y,x) + bool noFlip = !((theta > 45) && (theta < 135)); + + // select only sources close to the line + for (int j = 0; j < sources->n; j++) { + pmSource *source = sources->data[j]; + if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue; // ignore sources already found + + pmPeak *peak = source->peak; + if (noFlip) { + float yLine = (radius - peak->x*cosTheta) / sinTheta; + if (fabs(peak->y - yLine) > Rdel) continue; + psVectorAppend (xpos, peak->xf); + psVectorAppend (ypos, peak->yf); + } else { + float xLine = (radius - peak->y*sinTheta) / cosTheta; + if (fabs(peak->x - xLine) > Rdel) continue; + psVectorAppend (xpos, peak->yf); + psVectorAppend (ypos, peak->xf); + } + } + if (xpos->n < HOUGH_MIN_SOURCES) goto escape; + + // fit the line + if (!psVectorFitPolynomial1D (poly, NULL, 0, ypos, NULL, xpos)) { + psError(PS_ERR_UNKNOWN, true, "Failed to find a good trail fit"); + goto escape; + } + + // reset vector lengths + xpos->n = 0; + ypos->n = 0; + + // select sources close to the line and refit + // select only sources close to the line + for (int j = 0; j < sources->n; j++) { + pmSource *source = sources->data[j]; + if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue; // ignore sources already found + + pmPeak *peak = source->peak; + if (noFlip) { + float yLine = psPolynomial1DEval (poly, peak->x); + if (fabs(peak->y - yLine) > Rdel) continue; + psVectorAppend (xpos, peak->xf); + psVectorAppend (ypos, peak->yf); + } else { + float xLine = psPolynomial1DEval (poly, peak->y); + if (fabs(peak->x - xLine) > Rdel) continue; + psVectorAppend (xpos, peak->yf); + psVectorAppend (ypos, peak->xf); + } + } + if (xpos->n < HOUGH_MIN_SOURCES) goto escape; + + if (!psVectorFitPolynomial1D (poly, NULL, 0, ypos, NULL, xpos)) { + psError(PS_ERR_UNKNOWN, true, "Failed to find a good trail fit"); + goto escape; + } + + // mark sources close to the line + int nPeaksOnThisTrail = 0; + for (int j = 0; j < sources->n; j++) { + pmSource *source = sources->data[j]; + if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue; // ignore sources already found + + pmPeak *peak = source->peak; + if (noFlip) { + float yLine = psPolynomial1DEval (poly, peak->x); + if (fabs(peak->y - yLine) > Rdel) continue; + } else { + float xLine = psPolynomial1DEval (poly, peak->y); + if (fabs(peak->x - xLine) > Rdel) continue; + } + source->mode2 |= PM_SOURCE_MODE2_ON_LINE; + nPeaksOnThisTrail ++; + } + if (nPeaksOnThisTrail < HOUGH_MIN_SOURCES) goto escape; + + psFree (xpos); + psFree (ypos); + psFree (poly); + return nPeaksOnThisTrail; + +escape: + psFree (xpos); + psFree (ypos); + psFree (poly); + return false; +} + +enum {EDGE_LEFT, EDGE_BOTTOM, EDGE_RIGHT, EDGE_TOP}; + +void psphotFindFeatures_WhichEdges (float radius, float sinTheta, float cosTheta, int *edges, int NXpix, int NYpix) { + + edges[EDGE_TOP] = false; + edges[EDGE_LEFT] = false; + edges[EDGE_RIGHT] = false; + edges[EDGE_BOTTOM] = false; + + // r = 0 goes through (0,0) -> left edge + if ((radius >= 0) && (radius < NYpix*sinTheta)) { + edges[EDGE_LEFT] = true; + } + + // if (radius = NXpix*cosTheta) : goes through (Nx,0) -> bottom edge + if (cosTheta > 0.0) { + if ((radius > 0) && (radius <= NXpix*cosTheta)) { + edges[EDGE_BOTTOM] = true; + } + } else { + if ((radius < 0) && (radius >= NXpix*cosTheta)) { + edges[EDGE_BOTTOM] = true; + } + } + + // if (radius = NXpix*cosTheta + NYpix*sinTheta) : goes through (Nx,Ny) -> right edge + if ((radius > NXpix*cosTheta) && (radius <= NXpix*cosTheta + NYpix*sinTheta)) { + edges[EDGE_RIGHT] = true; + } + + // if (radius = NYpix*sinTheta) : goes through (0,Ny) -> top edge + if (cosTheta > 0.0) { + if ((radius >= NYpix*sinTheta) && (radius < NXpix*cosTheta + NYpix*sinTheta)) { + edges[EDGE_TOP] = true; + } + } else { + if ((radius <= NYpix*sinTheta) && (radius > NXpix*cosTheta + NYpix*sinTheta)) + edges[EDGE_TOP] = true; + } +} + +float psphotFindFeatures_LineLength (float radius, float theta, int *edges, int NXpix, int NYpix) { + + float cosTheta = DCOS(theta); + float sinTheta = DSIN(theta); + + psphotFindFeatures_WhichEdges (radius, sinTheta, cosTheta, edges, NXpix, NYpix); + + if (fabs(theta - 0.0) < 0.01) return ((radius >= 0) && (radius < +NXpix)) ? NYpix : 0.0; + if (fabs(theta - 180.0) < 0.01) return ((radius <= 0) && (radius > -NXpix)) ? NYpix : 0.0; + if (fabs(theta - 90.0) < 0.01) return ((radius >= 0) && (radius < +NYpix)) ? NXpix : 0.0; + + if (edges[EDGE_LEFT] && edges[EDGE_RIGHT]) { + float yLT = (radius - 0*cosTheta) / sinTheta; + float yRT = (radius - NXpix*cosTheta) / sinTheta; + float dY = yLT - yRT; + float lineLength = hypot(NXpix, dY); + return lineLength; + } + if (edges[EDGE_LEFT] && edges[EDGE_BOTTOM]) { + float yLT = (radius - 0*cosTheta) / sinTheta; + float xBT = (radius - 0*sinTheta) / cosTheta; + float lineLength = hypot(xBT, yLT); + return lineLength; + } + if (edges[EDGE_LEFT] && edges[EDGE_TOP]) { + float yLT = (radius - 0*cosTheta) / sinTheta; + float xTP = (radius - NYpix*sinTheta) / cosTheta; + float dY = NYpix - yLT; + float lineLength = hypot(xTP, dY); + return lineLength; + } + + if (edges[EDGE_BOTTOM] && edges[EDGE_TOP]) { + float xBT = (radius - 0*sinTheta) / cosTheta; + float xTP = (radius - NYpix*sinTheta) / cosTheta; + float dX = xTP - xBT; + float lineLength = hypot(dX, NYpix); + return lineLength; + } + if (edges[EDGE_BOTTOM] && edges[EDGE_RIGHT]) { + float xBT = (radius - 0*sinTheta) / cosTheta; + float yRT = (radius - NXpix*cosTheta) / sinTheta; + float dX = NXpix - xBT; + float lineLength = hypot(dX, yRT); + return lineLength; + } + if (edges[EDGE_TOP] && edges[EDGE_RIGHT]) { + float xTP = (radius - NYpix*sinTheta) / cosTheta; + float yRT = (radius - NXpix*cosTheta) / sinTheta; + float dX = NXpix - xTP; + float dY = NYpix - yRT; + float lineLength = hypot(dX, dY); + return lineLength; + } + return 0.0; +} Index: trunk/psphot/src/psphotFullForceArguments.c =================================================================== --- trunk/psphot/src/psphotFullForceArguments.c (revision 42841) +++ trunk/psphot/src/psphotFullForceArguments.c (revision 42842) @@ -3,4 +3,5 @@ static void usage(const char *program, psMetadata *arg, pmConfig *config, int exitCode); static void writeHelpInfo(const char* program, pmConfig* config, FILE* ofile); +void psphotSaveConfig (pmConfig *config); pmConfig *psphotFullForceArguments(int argc, char **argv) { @@ -111,4 +112,6 @@ // output position is fixed psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]); + + psphotSaveConfig (config); psTrace("psphot", 1, "Done with psphotFullForceArguments...\n"); Index: trunk/psphot/src/psphotFullForceReadout.c =================================================================== --- trunk/psphot/src/psphotFullForceReadout.c (revision 42841) +++ trunk/psphot/src/psphotFullForceReadout.c (revision 42842) @@ -121,5 +121,6 @@ // linear fit to include all sources (subtract again) // NOTE : apply to ALL sources (extended + psf) - psphotFitSourcesLinear (config, view, filerule, true, true); // pass 2 (detections->allSources) + // do NOT skip negative-flux sources (last arg = false) + psphotFitSourcesLinear (config, view, filerule, true, false); // pass 2 (detections->allSources) } Index: trunk/psphot/src/psphotFullForceSummary.c =================================================================== --- trunk/psphot/src/psphotFullForceSummary.c (revision 42841) +++ trunk/psphot/src/psphotFullForceSummary.c (revision 42842) @@ -322,7 +322,20 @@ int numCols = psMetadataLookupS32(&status, input->fpa->hdu->header, "IMNAXIS1"); int numRows = psMetadataLookupS32(&status, input->fpa->hdu->header, "IMNAXIS2"); + char *photcode = psMetadataLookupStr(&status, input->fpa->hdu->header, "PHOTCODE"); +// char *filter = psMetadataLookupStr(&status, input->fpa->hdu->header, "FPA.FILTER"); +// char *telescope = psMetadataLookupStr(&status, input->fpa->hdu->header, "FPA.TELESCOPE"); +// char *instrument = psMetadataLookupStr(&status, input->fpa->hdu->header, "FPA.INSTRUMENT"); +// char *detector = psMetadataLookupStr(&status, input->fpa->hdu->header, "FPA.DETECTOR"); + psMetadataAddS32(output->fpa->hdu->header, PS_LIST_TAIL, "IMNAXIS1", PS_META_REPLACE, "", numCols); psMetadataAddS32(output->fpa->hdu->header, PS_LIST_TAIL, "IMNAXIS2", PS_META_REPLACE, "", numRows); - + psMetadataAddStr(output->fpa->hdu->header, PS_LIST_TAIL, "PHOTCODE", PS_META_REPLACE, "", photcode); + + // XXX EAM 20240204 : these are added to the header psMetadata but do not make it out to the file header? + // psMetadataAddStr(output->fpa->hdu->header, PS_LIST_TAIL, "FPA.FILTER", PS_META_REPLACE, "", filter); + // psMetadataAddStr(output->fpa->hdu->header, PS_LIST_TAIL, "FPA.TELESCOPE", PS_META_REPLACE, "", telescope); + // psMetadataAddStr(output->fpa->hdu->header, PS_LIST_TAIL, "FPA.INSTRUMENT", PS_META_REPLACE, "", instrument); + // psMetadataAddStr(output->fpa->hdu->header, PS_LIST_TAIL, "FPA.DETECTOR", PS_META_REPLACE, "", detector); + // XXX: Also add psphot version information Index: trunk/psphot/src/psphotFullForceSummaryReadout.c =================================================================== --- trunk/psphot/src/psphotFullForceSummaryReadout.c (revision 42841) +++ trunk/psphot/src/psphotFullForceSummaryReadout.c (revision 42842) @@ -45,4 +45,7 @@ outputReadout = pmReadoutAlloc(outputCell); } + + // this requires the readout to exist -- and does not work + // psphotAddPhotcode (config, view, "PSPHOT.FULLFORCE.OUTPUT"); // Get the exposure parameters for the output from recipe and set them on the output Index: trunk/psphot/src/psphotMakeResiduals.c =================================================================== --- trunk/psphot/src/psphotMakeResiduals.c (revision 42841) +++ trunk/psphot/src/psphotMakeResiduals.c (revision 42842) @@ -1,3 +1,5 @@ # include "psphotInternal.h" +pmConfig *psphotGetConfig (); +bool _psphotSaveCube (char *basename, char *extname, psArray *cube); # define RESIDUAL_SOFTENING 0.005 @@ -93,4 +95,14 @@ psVector *yC = psVectorAllocEmpty (100, PS_TYPE_F32); +# define TESTRESID 0 +# if (TESTRESID) + psArray *sourceRawCube = psArrayAllocEmpty (1); + psArray *sourceSigCube = psArrayAllocEmpty (1); + psArray *sourceVarCube = psArrayAllocEmpty (1); + psArray *sourceMskCube = psArrayAllocEmpty (1); + + psImage *refImage = NULL; +# endif + // build (DATA - MODEL) [an image] for each psf star psArray *input = psArrayAllocEmpty (100); @@ -108,4 +120,24 @@ psImage *mask = psImageCopy (NULL, source->maskView ? source->maskView : source->maskObj, PS_TYPE_IMAGE_MASK); psImage *variance = psImageCopy (NULL, source->variance ? source->variance : source->pixels, PS_TYPE_F32); + +# if (TESTRESID) + psImage *tmpImage = psImageCopy(NULL, image, PS_TYPE_F32); + if (refImage == NULL) { + refImage = tmpImage; + psArrayAdd (sourceRawCube, 1, tmpImage); + } else { + if (tmpImage->numCols != refImage->numCols) { + fprintf (stderr, "mismatched image sizes (%d : %f, %f)\n", i, model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); + continue; + } + if (tmpImage->numRows != refImage->numRows) { + fprintf (stderr, "mismatched image sizes (%d : %f, %f)\n", i, model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]); + continue; + } + psArrayAdd (sourceRawCube, 1, tmpImage); + psFree (tmpImage); + } +# endif + pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal); @@ -114,4 +146,10 @@ psBinaryOp (image, image, "/", psScalarAlloc(Io, PS_TYPE_F32)); psBinaryOp (variance, variance, "/", psScalarAlloc(Io*Io, PS_TYPE_F32)); + +# if (TESTRESID) + psArrayAdd (sourceSigCube, 1, image); + psArrayAdd (sourceVarCube, 1, variance); + psArrayAdd (sourceMskCube, 1, mask); +# endif // we interpolate the image and variance - include the mask or not? @@ -142,4 +180,32 @@ psImageInit (resid->mask, 0); + +# if (TESTRESID) + psArray *inputSigCube = psArrayAllocEmpty (1); + psArray *inputVarCube = psArrayAllocEmpty (1); + psArray *inputMskCube = psArrayAllocEmpty (1); + + // below, we are generating the residual images with interpolation + // to save these in an output cube, I need to generate place-holder images here + for (int i = 0; i < input->n; i++) { + psImage *inputSig = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_F32); + psImage *inputVar = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_F32); + psImage *inputMsk = psImageAlloc (resid->Ro->numCols, resid->Ro->numRows, PS_TYPE_IMAGE_MASK); + + psImageInit (inputSig, NAN); + psImageInit (inputVar, NAN); + psImageInit (inputMsk, 0); + + psArrayAdd (inputSigCube, 1, inputSig); + psArrayAdd (inputVarCube, 1, inputVar); + psArrayAdd (inputMskCube, 1, inputMsk); + + // XXX free these here, right? + psFree (inputSig); + psFree (inputVar); + psFree (inputMsk); + } +# endif + // x(resid) = (x(image) - Xo)*xBin + xCenter @@ -148,7 +214,4 @@ psVector *fmasks = psVectorAlloc (input->n, PS_TYPE_VECTOR_MASK); - // statistic to use to determine baseline for clipping - psStats *fluxClip = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); - psStats *fluxClipDef = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV); // statistic to use to determine output flux // XXX make API to convert statOption for MEAN/MEDIAN in to corresponding STDEV? @@ -156,4 +219,8 @@ psStats *fluxStatsDef = psStatsAlloc (statOption | PS_STAT_SAMPLE_STDEV); + // use this for the IRLS fitting below (defines niter) + psStats *statsIRLS = psStatsAlloc(PS_STAT_CLIPPED_MEAN); statsIRLS->clipIter = 10; // max number of iterations + psPolynomial2D *polyIRLS = psPolynomial2DAlloc(PS_POLYNOMIAL_ORD, 1, 1); + // Use psF64 to minimize overflow problems? psImage *A = psImageAlloc(3, 3, PS_TYPE_F64); // Least-squares matrix @@ -166,4 +233,11 @@ int nGoodPixel = 0; // pixel is off the image + + // skip pixels outside of the requested radius + float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows)); + if (radius > radiusMax) { + resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; + continue; + } // build the vector of data values for this output pixel @@ -194,5 +268,15 @@ nGoodPixel ++; } + +# if (TESTRESID) + psImage *inputSig = inputSigCube->data[i]; inputSig->data.F32[oy][ox] = flux; + psImage *inputVar = inputVarCube->data[i]; inputVar->data.F32[oy][ox] = flux; + psImage *inputMsk = inputMskCube->data[i]; inputMsk->data.PS_TYPE_IMAGE_MASK_DATA[oy][ox] = fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]; +# endif } + + if (nGoodPixel < 0.05*input->n) { + fprintf (stderr, "warning: %d, %d : residual pixel with few input pixels: %d vs %d\n", ox, oy, nGoodPixel, (int) input->n); + } // skip pixels with insufficient data @@ -207,15 +291,21 @@ // measure the robust median to determine a baseline reference value - *fluxClip = *fluxClipDef; + psStats *fluxClip = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV); if (!psVectorStats (fluxClip, fluxes, NULL, fmasks, fmaskVal)) { psError(PSPHOT_ERR_CONFIG, false, "Error calculating residual stats"); + psFree (fluxClip); return false; } - if (isnan(fluxClip->robustMedian)) { + if (isnan(fluxClip->sampleMedian)) { resid->Ro->data.F32[oy][ox] = 0.0; resid->Rx->data.F32[oy][ox] = 0.0; resid->Ry->data.F32[oy][ox] = 0.0; resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask; + psFree (fluxClip); continue; + } + + if (fabs(fluxClip->sampleMedian) > 1.5) { + fprintf (stderr, "warning: %d, %d : residual pixel has funny initial median value: %f\n", ox, oy, fluxClip->sampleMedian); } @@ -223,5 +313,5 @@ int nKeep = 0; for (int i = 0; i < fluxes->n; i++) { - float delta = fluxes->data.F32[i] - fluxClip->robustMedian; + float delta = fluxes->data.F32[i] - fluxClip->sampleMedian; float sigma = sqrt (dfluxes->data.F32[i]); float swing = fabs(delta) / sigma; @@ -233,4 +323,9 @@ if (!fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) nKeep++; } + psFree (fluxClip); + + if (nKeep < 5) { + fprintf (stderr, "warning: %d, %d : residual pixel with few good pixels: %d vs %d\n", ox, oy, nKeep, (int) input->n); + } if (SPATIAL_ORDER == 0) { @@ -242,17 +337,9 @@ } - float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows)); - if (radius > radiusMax) { - resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; - continue; - } + resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0; resid->Ro->data.F32[oy][ox] = psStatsGetValue(fluxStats, statOption); - resid->Rx->data.F32[oy][ox] = resid->Ry->data.F32[oy][ox] = 0.0; - if (isnan(resid->Ro->data.F32[oy][ox])) { resid->Ro->data.F32[oy][ox] = 0.0; - resid->Rx->data.F32[oy][ox] = 0.0; - resid->Ry->data.F32[oy][ox] = 0.0; resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = badMask; continue; @@ -265,52 +352,71 @@ assert (SPATIAL_ORDER == 1); - float radius = hypot((ox - 0.5*resid->Ro->numCols), (oy - 0.5*resid->Ro->numRows)); - if (radius > radiusMax) { - resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; - continue; - } - - psImageInit(A, 0.0); - psVectorInit(B, 0.0); - for (int i = 0; i < fluxes->n; i++) { - if (fmasks->data.PS_TYPE_VECTOR_MASK_DATA[i]) continue; - B->data.F64[0] += fluxes->data.F32[i]/dfluxes->data.F32[i]; - B->data.F64[1] += fluxes->data.F32[i]*xC->data.F32[i]/dfluxes->data.F32[i]; - B->data.F64[2] += fluxes->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i]; - - A->data.F64[0][0] += 1.0/dfluxes->data.F32[i]; - A->data.F64[1][0] += xC->data.F32[i]/dfluxes->data.F32[i]; - A->data.F64[2][0] += yC->data.F32[i]/dfluxes->data.F32[i]; - - A->data.F64[1][1] += PS_SQR(xC->data.F32[i])/dfluxes->data.F32[i]; - A->data.F64[2][2] += PS_SQR(yC->data.F32[i])/dfluxes->data.F32[i]; - A->data.F64[1][2] += xC->data.F32[i]*yC->data.F32[i]/dfluxes->data.F32[i]; - } - - A->data.F64[0][1] = A->data.F64[1][0]; - A->data.F64[0][2] = A->data.F64[2][0]; - A->data.F64[2][1] = A->data.F64[1][2]; - - if (!psMatrixGJSolve(A, B)) { - resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; - psWarning("Singular matrix solving for (y,x) = (%d,%d)'s residuals, masking", oy, ox); - continue; - } - - resid->Ro->data.F32[oy][ox] = B->data.F64[0]; - resid->Rx->data.F32[oy][ox] = B->data.F64[1]; - resid->Ry->data.F32[oy][ox] = B->data.F64[2]; - - float dRo = sqrt(A->data.F32[0][0]); + // we have the following vectors to describe this pixel: + // fluxes, dfluxes, fmasks, xC, yC + + psVectorIRLSFitPolynomial2D (polyIRLS, statsIRLS, fmasks, fmaskVal, fluxes, dfluxes, xC, yC); + + resid->Ro->data.F32[oy][ox] = polyIRLS->coeff[0][0]; + resid->Rx->data.F32[oy][ox] = polyIRLS->coeff[1][0]; + resid->Ry->data.F32[oy][ox] = polyIRLS->coeff[0][1]; + + // is the error well-defined? + float dRo = polyIRLS->coeffErr[0][0]; + + if (fabs(resid->Ro->data.F32[oy][ox]) > 1.5) { + fprintf (stderr, "warning: %d, %d : residual pixel has funny fit value: %f\n", ox, oy, resid->Ro->data.F32[oy][ox]); + } if (fabs(resid->Ro->data.F32[oy][ox]) < pixelSN*dRo/sqrt(nKeep)) { + // fprintf (stderr, "mask : %d, %d : %f +/- %f (%f for %d) : %d\n", ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]); resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox] = 1; resid->Ro->data.F32[oy][ox] = 0.0; resid->Rx->data.F32[oy][ox] = 0.0; resid->Ry->data.F32[oy][ox] = 0.0; - } + } else { + // fprintf (stderr, "keep : %d, %d : %f +/- %f (%f for %d) : %d\n", ox, oy, resid->Ro->data.F32[oy][ox], dRo, dRo/sqrt(nKeep), nKeep, resid->mask->data.PM_TYPE_RESID_MASK_DATA[oy][ox]); + } } } } + +# if (TESTRESID) + + pmConfig *config = psphotGetConfig(); + char *output = psMetadataLookupStr (&status, config->arguments, "OUTPUT"); + + // XXX // save the star coordinate vectors + // XXX + // XXX psString filename = NULL; + // XXX psStringAppend (&filename, "%s.instar.dat", output); + + char filename[1024]; + sprintf (filename, "%s.instar.dat", output); + + FILE *fout = fopen (filename, "w"); + for (int i = 0; i < xC->n; i++) { + fprintf (fout, "%d %f %f\n", i, xC->data.F32[i], yC->data.F32[i]); + } + fclose (fout); + + _psphotSaveCube (output, "raw.raw.fits", sourceRawCube); + _psphotSaveCube (output, "sig.raw.fits", sourceSigCube); + _psphotSaveCube (output, "var.raw.fits", sourceVarCube); + _psphotSaveCube (output, "msk.raw.fits", sourceMskCube); + + _psphotSaveCube (output, "sig.fits", inputSigCube); + _psphotSaveCube (output, "var.fits", inputVarCube); + _psphotSaveCube (output, "msk.fits", inputMskCube); + + psFree (sourceRawCube); + psFree (sourceSigCube); + psFree (sourceVarCube); + psFree (sourceMskCube); + + psFree (inputSigCube); + psFree (inputVarCube); + psFree (inputMskCube); + +# endif psFree (A); @@ -329,6 +435,4 @@ psFree (fluxStats); psFree (fluxStatsDef); - psFree (fluxClip); - psFree (fluxClipDef); if (resid != NULL && psTraceGetLevel("psphot") > 5) { @@ -340,5 +444,24 @@ } + psFree (statsIRLS); + psFree (polyIRLS); + psf->residuals = resid; return (resid != NULL) ? true : false; } + +bool _psphotSaveCube (char *basename, char *extname, psArray *cube) { + + psString filename = NULL; + psStringAppend (&filename, "%s.%s", basename, extname); + + psFits *fits = psFitsOpen (filename, "w"); + if (!psFitsWriteImageCube (fits, NULL, cube, NULL)) { + fprintf (stderr, "failed to write the cube %s\n", filename); + } + psFitsClose (fits); + + psFree (filename); + + return true; +} Index: trunk/psphot/src/psphotMaskReadout.c =================================================================== --- trunk/psphot/src/psphotMaskReadout.c (revision 42841) +++ trunk/psphot/src/psphotMaskReadout.c (revision 42842) @@ -58,4 +58,26 @@ return false; } + } + + // XXX test to repair variance + if (true) { + psImage *im = readout->image; + psImage *wt = readout->variance; + float maxVar = 0; + // find the max value + for (int j = 0; j < im->numRows; j++) { + for (int i = 0; i < im->numCols; i++) { + if (!isfinite(wt->data.F32[j][i])) continue; + maxVar = PS_MAX (wt->data.F32[j][i], maxVar); + } + } + // set pixels with nan variance and non-nan image to maxVar + for (int j = 0; j < im->numRows; j++) { + for (int i = 0; i < im->numCols; i++) { + if (!isfinite(im->data.F32[j][i])) continue; + if (isfinite(wt->data.F32[j][i])) continue; + wt->data.F32[j][i] = maxVar; + } + } } Index: trunk/psphot/src/psphotOutput.c =================================================================== --- trunk/psphot/src/psphotOutput.c (revision 42841) +++ trunk/psphot/src/psphotOutput.c (revision 42842) @@ -1,3 +1,13 @@ # include "psphotInternal.h" + +pmConfig *staticConfig = NULL; + +void psphotSaveConfig (pmConfig *config) { + staticConfig = config; +} + +pmConfig *psphotGetConfig () { + return staticConfig; +} // convert filerule to filerule.NUM and look up in the config->arguments metadata Index: trunk/psphot/src/psphotRoughClass.c =================================================================== --- trunk/psphot/src/psphotRoughClass.c (revision 42841) +++ trunk/psphot/src/psphotRoughClass.c (revision 42842) @@ -198,5 +198,9 @@ } - psfClump = pmSourcePSFClump (NULL, region, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX); + // use clipped stats or robust stats to define the clump + bool PSF_CLUMP_USE_CLIPPED_STATS = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_USE_CLIPPED_STATS"); + if (!status) { PSF_CLUMP_USE_CLIPPED_STATS = true; } + + psfClump = pmSourcePSFClump (NULL, region, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX, PSF_CLUMP_USE_CLIPPED_STATS); psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.X", PS_META_REPLACE, "psf clump center", psfClump.X); Index: trunk/psphot/src/psphotSourceStats.c =================================================================== --- trunk/psphot/src/psphotSourceStats.c (revision 42841) +++ trunk/psphot/src/psphotSourceStats.c (revision 42842) @@ -120,4 +120,5 @@ source->moments = pmMomentsAlloc(); + // XXX can this ever be set at this point? (we just allocated the source and moments) if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) { fprintf (stderr, "moment failure\n"); @@ -144,4 +145,7 @@ if (setWindow) { + // identify sources on lines, drop from MomentsWindow analysis + psphotFindFeatures (readout, sources); + if (!psphotSetMomentsWindow(recipe, readout->analysis, sources, maskVal)) { psError(PS_ERR_UNEXPECTED_NULL, false, "Failed to determine Moments Window!"); @@ -503,4 +507,8 @@ } + // use clipped stats or robust stats to define the clump + bool PSF_CLUMP_USE_CLIPPED_STATS = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_USE_CLIPPED_STATS"); + if (!status) { PSF_CLUMP_USE_CLIPPED_STATS = true; } + // when we set the window, we are not attempting to measure spatial variations; we can use a somewhat higher S/N limit // since we are using all sources (true?) @@ -539,11 +547,23 @@ sources = psArraySort (sources, pmSourceSortByFlux); + // generate an array of the sources which we want to use for this analysis + // XXX move max source number to config + psArray *sourceSubset = psArrayAllocEmpty (100); + for (int i = 0; (i < sources->n) && (i < 400); i++) { + + pmSource *source = sources->data[i]; + + // skip faint sources for moments measurement + if (sqrt(source->peak->detValue) < MIN_SN) continue; + + // skip sources on lines + if (source->mode2 & PM_SOURCE_MODE2_ON_LINE) continue; + psArrayAdd (sourceSubset, 100, source); + } + // loop over radii: for (int i = 0; i < nsigma; i++) { - - // XXX move max source number to config - for (int j = 0; (j < sources->n) && (j < 400); j++) { - - pmSource *source = sources->data[j]; + for (int j = 0; j < sourceSubset->n; j++) { + pmSource *source = sourceSubset->data[j]; psAssert (source->moments, "force moments to exist"); source->moments->nPixels = 0; @@ -565,6 +585,6 @@ sprintf (name, "moments.v%d.dat", i); FILE *fout = fopen (name, "w"); - for (int j = 0; j < sources->n; j++) { - pmSource *source = sources->data[j]; + for (int j = 0; j < sourceSubset->n; j++) { + pmSource *source = sourceSubset->data[j]; psAssert (source->moments, "force moments to exist"); source->moments->nPixels = 0; @@ -581,8 +601,8 @@ // determine the PSF parameters from the source moment values - pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX); + pmPSFClump psfClump = pmSourcePSFClump (NULL, NULL, sourceSubset, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_SX_MIN, MOMENTS_SY_MIN, MOMENTS_AR_MAX, PSF_CLUMP_USE_CLIPPED_STATS); psLogMsg ("psphot", 3, "sigma guess (pix) %.1f, nStars: %d of %d in clump, nSigma: %5.2f, X, Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nTotal, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]); - Rmin[i] = pmSourceMinKronRadius(sources, PSF_SN_LIM); + Rmin[i] = pmSourceMinKronRadius(sourceSubset, PSF_SN_LIM); #if 0 @@ -600,5 +620,5 @@ psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY); if (pmVisualTestLevel("psphot.moments.full", 2)) { - psphotVisualPlotMoments (recipe, analysis, sources); + psphotVisualPlotMoments (recipe, analysis, sourceSubset); } #endif @@ -608,4 +628,6 @@ Sout[i] = (Nout[i] == 0) ? NAN : sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i]; } + + psFree (sourceSubset); // we are looking for sigma for which Sout = 0.65 (or some other value) Index: trunk/psphot/src/psphotVisual.c =================================================================== --- trunk/psphot/src/psphotVisual.c (revision 42841) +++ trunk/psphot/src/psphotVisual.c (revision 42842) @@ -298,4 +298,5 @@ bool status = false; pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND"); + if (!file) return false; pmReadout *backgnd = READOUT_OR_INTERNAL(view, file); @@ -394,4 +395,5 @@ psArray *peaks = detections->peaks; + if (!peaks) return false; // note: this uses the Ohana allocation tools: