Changeset 42842

Timestamp:

May 9, 2025, 11:21:40 AM (14 months ago)

Author:

eugene

Message:

add psphotFindFeatures; enable forced-photometry between cameras; do NOT skip negative-flux sources in full-force photometry; make PSF residual component generation more robust: IRLS vs median; exclude pixels with poor values (scaled value > 1.5, relative to normalized PSF); exclude pixels with too few inputs (<5); allow PSF clump to use clipped stats or robust stats

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src/Makefile.am (modified) (1 diff)
• src/psphot.h (modified) (1 diff)
• src/psphotArguments.c (modified) (2 diffs)
• src/psphotCheckStarDistribution.c (modified) (1 diff)
• src/psphotChoosePSF.c (modified) (1 diff)
• src/psphotDefineFiles.c (modified) (7 diffs)
• src/psphotFindFeatures.c (copied) (copied from branches/eam_branches/ipp-20230313/psphot/src/psphotFindFeatures.c )
• src/psphotFullForceArguments.c (modified) (2 diffs)
• src/psphotFullForceReadout.c (modified) (1 diff)
• src/psphotFullForceSummary.c (modified) (1 diff)
• src/psphotFullForceSummaryReadout.c (modified) (1 diff)
• src/psphotMakeResiduals.c (modified) (16 diffs)
• src/psphotMaskReadout.c (modified) (1 diff)
• src/psphotOutput.c (modified) (1 diff)
• src/psphotRoughClass.c (modified) (1 diff)
• src/psphotSourceStats.c (modified) (8 diffs)
• src/psphotVisual.c (modified) (2 diffs)

trunk/psphot/src/Makefile.am

@@ -247 +247 @@
         psphotSetNFrames.c             \
         psphotSourceMemory.c           \
+        psphotFindFeatures.c           \
         psphotChipParams.c              \
         psphotGalaxyParams.c

trunk/psphot/src/psphot.h

@@ -592 +592 @@
 extern bool psphotINpsphotStack;
 
+bool psphotFindFeatures (pmReadout *readout, psArray *sources);
+
 
 #endif

trunk/psphot/src/psphotArguments.c

@@ -97 +97 @@
 }
 
+void psphotSaveConfig (pmConfig *config);
+
 pmConfig *psphotArguments(int argc, char **argv) {
 
@@ -244 +246 @@
     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);

trunk/psphot/src/psphotCheckStarDistribution.c

@@ -89 +89 @@
         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;

trunk/psphot/src/psphotChoosePSF.c

@@ -114 +114 @@
     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

trunk/psphot/src/psphotDefineFiles.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+pmConfig *pmConfigMakeTemp (pmConfig *config);
 
 bool psphotINpsphotStack = false;
@@ -143 +144 @@
     }
 
+    // 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 +161 @@
 
     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 +170 @@
         // 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 +177 @@
 
     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 +183 @@
         }
     }
+
+    psFree (srcconfig);
 
     // XXX add in example PSF image thumbnails
@@ -214 +226 @@
     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);
+}
+

trunk/psphot/src/psphotFullForceArguments.c

@@ -3 +3 @@
 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 +112 @@
     // output position is fixed
     psMetadataAddStr (config->arguments, PS_LIST_TAIL, "OUTPUT", 0, "", argv[1]);
+
+    psphotSaveConfig (config);
 
     psTrace("psphot", 1, "Done with psphotFullForceArguments...\n");

trunk/psphot/src/psphotFullForceReadout.c

@@ -121 +121 @@
       // 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)
     }
 

trunk/psphot/src/psphotFullForceSummary.c

@@ -322 +322 @@
     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
 

trunk/psphot/src/psphotFullForceSummaryReadout.c

@@ -45 +45 @@
         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

trunk/psphot/src/psphotMakeResiduals.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+pmConfig *psphotGetConfig ();
+bool _psphotSaveCube (char *basename, char *extname, psArray *cube);
 
 # define RESIDUAL_SOFTENING 0.005 
@@ -93 +95 @@
     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 +120 @@
         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 +146 @@
         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 +180 @@
     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 +214 @@
     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 +219 @@
     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 +233 @@
 
             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 +268 @@
                     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 +291 @@
 
             // 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 +313 @@
             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 +323 @@
                 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 +337 @@
                 }
 
-                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 +352 @@
                 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 +435 @@
     psFree (fluxStats);
     psFree (fluxStatsDef);
-    psFree (fluxClip);
-    psFree (fluxClipDef);
 
     if (resid != NULL && psTraceGetLevel("psphot") > 5) {
@@ -340 +444 @@
     }
 
+    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;
+}

trunk/psphot/src/psphotMaskReadout.c

@@ -58 +58 @@
             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;
+        }
+      }
     }
 

trunk/psphot/src/psphotOutput.c

@@ -1 +1 @@
 # 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

trunk/psphot/src/psphotRoughClass.c

@@ -198 +198 @@
         }
 
-        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);

trunk/psphot/src/psphotSourceStats.c

@@ -120 +120 @@
         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 +145 @@
 
     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 +507 @@
     }
 
+    // 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 +547 @@
     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 +585 @@
           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 +601 @@
 
         // 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 +620 @@
         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 +628 @@
         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)

trunk/psphot/src/psphotVisual.c

@@ -298 +298 @@
     bool status = false;
     pmFPAfile *file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND");
+    if (!file) return false;
 
     pmReadout *backgnd = READOUT_OR_INTERNAL(view, file);
@@ -394 +395 @@
 
     psArray *peaks = detections->peaks;
+    if (!peaks) return false;
 
     // note: this uses the Ohana allocation tools: