Changeset 28484 for branches/pap/psphot

Timestamp:

Jun 24, 2010, 2:59:09 PM (16 years ago)

Author:

Paul Price

Message:

Merging trunk in advance of reintegrating into trunk.

Location:

branches/pap

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src/psphotApResid.c (modified) (14 diffs)
• psphot/src/psphotBlendFit.c (modified) (7 diffs)
• psphot/src/psphotFitSourcesLinear.c (modified) (7 diffs)
• psphot/src/psphotFitSourcesLinearStack.c (modified) (5 diffs)
• psphot/src/psphotGuessModels.c (modified) (2 diffs)
• psphot/src/psphotImageLoop.c (modified) (3 diffs)
• psphot/src/psphotMagnitudes.c (modified) (5 diffs)
• psphot/src/psphotRadiusChecks.c (modified) (6 diffs)
• psphot/src/psphotReadout.c (modified) (4 diffs)
• psphot/src/psphotRoughClass.c (modified) (5 diffs)
• psphot/src/psphotSourceStats.c (modified) (10 diffs)
• psphot/src/psphotTest.c (modified) (5 diffs)

branches/pap/psphot/src/psphotApResid.c

@@ -22 +22 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (i == chisqNum) continue; // skip chisq image
-        if (!psphotApResidReadout (config, view, filerule, i, recipe)) {
+        if (i == chisqNum) continue; // skip chisq image
+        if (!psphotApResidReadout (config, view, filerule, i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for %s entry %d", filerule, i);
-            return false;
-        }
+            return false;
+        }
     }
     return true;
@@ -56 +56 @@
 
     if (!sources->n) {
-        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping ap resid");
-        return true;
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping ap resid");
+        return true;
     }
 
@@ -66 +66 @@
     int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
     if (!status) {
-        nThreads = 0;
+        nThreads = 0;
     }
 
@@ -128 +128 @@
     for (int i = 0; i < cellGroups->n; i++) {
 
-        psArray *cells = cellGroups->data[i];
-
-        for (int j = 0; j < cells->n; j++) {
-
-            // allocate a job -- if threads are not defined, this just runs the job
-            psThreadJob *job = psThreadJobAlloc ("PSPHOT_APRESID_MAGS");
-
-            psArrayAdd(job->args, 1, cells->data[j]); // sources
-            psArrayAdd(job->args, 1, psf);
-            PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
-            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
-
-            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
-            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nfail
-
-            if (!psThreadJobAddPending(job)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                psFree (job);
-                return false;
-            }
-            psFree(job);
-        }
-
-        // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-            return false;
-        }
-
-        // we have only supplied one type of job, so we can assume the types here
-        psThreadJob *job = NULL;
-        while ((job = psThreadJobGetDone()) != NULL) {
-            if (job->args->n < 1) {
-                fprintf (stderr, "error with job\n");
-            } else {
-                psScalar *scalar = NULL;
-                scalar = job->args->data[5];
-                Nskip += scalar->data.S32;
-                scalar = job->args->data[6];
-                Nfail += scalar->data.S32;
-            }
-            psFree(job);
-        }
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_APRESID_MAGS");
+
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, psf);
+            PS_ARRAY_ADD_SCALAR(job->args, photMode, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
+
+            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nskip
+            PS_ARRAY_ADD_SCALAR(job->args, 0,        PS_TYPE_S32); // this is used as a return value for Nfail
+
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                return false;
+            }
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                psScalar *scalar = NULL;
+                scalar = job->args->data[5];
+                Nskip += scalar->data.S32;
+                scalar = job->args->data[6];
+                Nfail += scalar->data.S32;
+            }
+            psFree(job);
+        }
     }
 
@@ -184 +182 @@
     Npsf = 0;
 
-# ifdef DEBUG    
+# ifdef DEBUG
     FILE *f = fopen ("apresid.dat", "w");
     psAssert (f, "failed open");
@@ -199 +197 @@
         if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
 
-        if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");
-        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
-        if (source->mode &  PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT");
-            
+        if (source->mode &  PM_SOURCE_MODE_EXT_LIMIT) SKIPSTAR ("EXTENDED");
+        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) SKIPSTAR ("COSMIC RAY");
+        if (source->mode &  PM_SOURCE_MODE_DEFECT) SKIPSTAR ("DEFECT");
+
         if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
             continue;
         }
 
-        // XXX make this user-configurable?
-        if (source->errMag > 0.01) continue;
+        // XXX make this user-configurable?
+        if (source->errMag > 0.01) continue;
 
         // aperture residual for this source
@@ -221 +219 @@
 
 # ifdef DEBUG
-        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
-                 source->peak->xf, source->peak->yf, 
-                 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS], 
-                 source->psfMag, source->apMag, source->errMag,
-                 source->modelPSF->params->data.F32[PM_PAR_I0], 
-                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY], 
-                 source->modelPSF->params->data.F32[PM_PAR_7]);
+        fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
+                 source->peak->xf, source->peak->yf,
+                 source->modelPSF->params->data.F32[PM_PAR_XPOS], source->modelPSF->params->data.F32[PM_PAR_YPOS],
+                 source->psfMag, source->apMag, source->errMag,
+                 source->modelPSF->params->data.F32[PM_PAR_I0],
+                 source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], source->modelPSF->params->data.F32[PM_PAR_SYY],
+                 source->modelPSF->params->data.F32[PM_PAR_7]);
 # endif
-        if (!isfinite(source->psfMag)) psAbort ("nan in psfMag");
-        if (!isfinite(source->errMag)) psAbort ("nan in errMag");
-        if (!isfinite(source->apMag)) psAbort ("nan in apMag");
-        if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos");
-        if (!isfinite(model->params->data.F32[PM_PAR_YPOS])) psAbort ("nan in yPos");
+        if (!isfinite(source->psfMag)) psAbort ("nan in psfMag");
+        if (!isfinite(source->errMag)) psAbort ("nan in errMag");
+        if (!isfinite(source->apMag)) psAbort ("nan in apMag");
+        if (!isfinite(model->params->data.F32[PM_PAR_XPOS])) psAbort ("nan in xPos");
+        if (!isfinite(model->params->data.F32[PM_PAR_YPOS])) psAbort ("nan in yPos");
 
         psVectorAppend (mag, source->psfMag);
@@ -253 +251 @@
     if (Npsf < APTREND_NSTAR_MIN) {
         psWarning("Only %d valid aperture residual sources (need %d), giving up", Npsf, APTREND_NSTAR_MIN);
-        goto escape;
-    }
-
-    // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different 
-    // definition for 'order' (order_MAP = order_POLY + 1).  in addition, we have a 
+        goto escape;
+    }
+
+    // this is a bit tricky, because we have two cases (MAP vs POLY), and they have a different
+    // definition for 'order' (order_MAP = order_POLY + 1).  in addition, we have a
     // user-specified MAX order, which we should respect, regardless of the mode
 
@@ -270 +268 @@
     pmTrend2DMode mode = PM_TREND_MAP;
     if (mode == PM_TREND_MAP) {
-        MaxOrderForStars ++;
-    } 
+        MaxOrderForStars ++;
+    }
     APTREND_ORDER_MAX = PS_MIN (APTREND_ORDER_MAX, MaxOrderForStars);
 
@@ -283 +281 @@
     int NY = readout->image->numRows;
     for (int i = 1; i <= APTREND_ORDER_MAX; i++) {
-        
-        int Nx, Ny;
-        if (NX > NY) {
-            Nx = i;
-            Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5));
-        } else {
-            Ny = i;
-            Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5));
-        }
-
-        float errorFloor;
+
+        int Nx, Ny;
+        if (NX > NY) {
+            Nx = i;
+            Ny = PS_MAX (1, (int)(i * (NY / (float)(NX)) + 0.5));
+        } else {
+            Ny = i;
+            Nx = PS_MAX (1, (int)(i * (NX / (float)(NY)) + 0.5));
+        }
+
+        float errorFloor;
         pmTrend2D *apTrend = psphotApResidTrend (&errorFloor, readout, Nx, Ny, xPos, yPos, apResid, dMag);
-        if (!apTrend) {
-            continue;
-        }
-
-        // apply ApTrend results
-        // float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
-        // float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
-        // float ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
-        // if (!isfinite(ApResid)) psAbort("nan apresid @ center");
+        if (!apTrend) {
+            continue;
+        }
+
+        // apply ApTrend results
+        // float xc = 0.5*readout->image->numCols + readout->image->col0 + 0.5;
+        // float yc = 0.5*readout->image->numRows + readout->image->row0 + 0.5;
+        // float ApResid = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
+        // if (!isfinite(ApResid)) psAbort("nan apresid @ center");
 
         // store the minimum errorFloor and best ApTrend to keep
         if (errorFloor < errorFloorMin) {
             errorFloorMin = errorFloor;
-            psFree (psf->ApTrend);
-            psf->ApTrend = psMemIncrRefCounter(apTrend);
-        }
-        psFree (apTrend);
+            psFree (psf->ApTrend);
+            psf->ApTrend = psMemIncrRefCounter(apTrend);
+        }
+        psFree (apTrend);
     }
     if (psf->ApTrend == NULL) {
         psWarning("Failed to find a valid aperture residual value");
-        goto escape;
+        goto escape;
     }
 
@@ -383 +381 @@
     if (!pmTrend2DFit (apTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft)) {
         psWarning("Failed to fit trend for %d x %d map", Nx, Ny);
-        psFree (apTrend);
-        return NULL;
+        psFree (apTrend);
+        return NULL;
     }
     if (apTrend->mode == PM_TREND_MAP) {
-        // p_psImagePrint (2, apTrend->map->map, "ApTrend Before"); // XXX TEST:
-        psImageMapRepair (apTrend->map->map);
-        // p_psImagePrint (2, apTrend->map->map, "ApTrend After"); // XXX TEST:
+        // p_psImagePrint (2, apTrend->map->map, "ApTrend Before"); // XXX TEST:
+        psImageMapRepair (apTrend->map->map);
+        // p_psImagePrint (2, apTrend->map->map, "ApTrend After"); // XXX TEST:
     }
 
@@ -400 +398 @@
     if (!isfinite(*apResidSysErr)) {
         psWarning("Failed to find systematic error for %d x %d map", Nx, Ny);
-        psFree (apTrend);
-        return NULL;
+        psFree (apTrend);
+        return NULL;
     }
 
@@ -408 +406 @@
 
     if (psTraceGetLevel("psphot") >= 4) {
-        char filename[64];
-        snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny);
+        char filename[64];
+        snprintf (filename, 64, "apresid.%dx%d.dat", Nx, Ny);
         FILE *dumpFile = fopen (filename, "w");
         for (int i = 0; i < xPos->n; i++) {
@@ -457 +455 @@
         }
 
-        // clear the mask bit and set the circular mask pixels
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
-
-        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal);
-
-        // clear the mask bit 
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        // clear the mask bit and set the circular mask pixels
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
+
+        bool status = pmSourceMagnitudes (source, psf, photMode, maskVal);
+
+        // clear the mask bit
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
 
         // re-subtract the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
-        if (!status) {
-            Nskip ++;
-            psTrace ("psphot", 3, "skip : bad source mag");
-            continue;
-        }
-    
-        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
-            Nfail ++;
-            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
-            continue;
-        }
-        source->mode |= PM_SOURCE_MODE_AP_MAGS;
+        if (!status) {
+            Nskip ++;
+            psTrace ("psphot", 3, "skip : bad source mag");
+            continue;
+        }
+
+        if (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
+            Nfail ++;
+            psTrace ("psphot", 3, "fail : nan mags : %f %f", source->apMag, source->psfMag);
+            continue;
+        }
+        source->mode |= PM_SOURCE_MODE_AP_MAGS;
     }
 

branches/pap/psphot/src/psphotBlendFit.c

@@ -15 +15 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (!psphotBlendFitReadout (config, view, filerule, i, recipe)) {
+        if (!psphotBlendFitReadout (config, view, filerule, i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (non-linear) for %s entry %d", filerule, i);
-            return false;
-        }
+            return false;
+        }
     }
     return true;
@@ -48 +48 @@
 
     if (!sources->n) {
-        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend fit");
-        return true;
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend fit");
+        return true;
     }
 
@@ -87 +87 @@
     sources = psArraySort (sources, pmSourceSortBySN);
     if (!sources->n) {
-        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
-        return true;
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
+        return true;
     }
 
@@ -103 +103 @@
         for (int j = 0; j < cells->n; j++) {
 
-            // allocate a job -- if threads are not defined, this just runs the job
-            psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
-            psArray *newSources = psArrayAllocEmpty(16);
-
-            psArrayAdd(job->args, 1, readout);
-            psArrayAdd(job->args, 1, recipe);
-            psArrayAdd(job->args, 1, cells->data[j]); // sources
-            psArrayAdd(job->args, 1, psf);
-            psArrayAdd(job->args, 1, newSources); // return for new sources
-            psFree (newSources);
-
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
-            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
-
-            if (!psThreadJobAddPending(job)) {
-                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                psFree (job);
-                return NULL;
-            }
-            psFree(job);
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_BLEND_FIT");
+            psArray *newSources = psArrayAllocEmpty(16);
+
+            psArrayAdd(job->args, 1, readout);
+            psArrayAdd(job->args, 1, recipe);
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, psf);
+            psArrayAdd(job->args, 1, newSources); // return for new sources
+            psFree (newSources);
+
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfit
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Npsf
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Next
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nfail
+
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                return NULL;
+            }
 
 # if (0)
@@ -152 +150 @@
         }
 
-        // wait for the threads to finish and manage results
-        if (!psThreadPoolWait (false)) {
-            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-            return NULL;
-        }
-
-        // we have only supplied one type of job, so we can assume the types here
-        psThreadJob *job = NULL;
-        while ((job = psThreadJobGetDone()) != NULL) {
-            if (job->args->n < 1) {
-                fprintf (stderr, "error with job\n");
-            } else {
-                psScalar *scalar = NULL;
-                scalar = job->args->data[5];
-                Nfit += scalar->data.S32;
-                scalar = job->args->data[6];
-                Npsf += scalar->data.S32;
-                scalar = job->args->data[7];
-                Next += scalar->data.S32;
-                scalar = job->args->data[8];
-                Nfail += scalar->data.S32;
-
-                // add these back onto sources
-                psArray *newSources = job->args->data[4];
-                for (int j = 0; j < newSources->n; j++) {
-                    psArrayAdd (sources, 16, newSources->data[j]);
-                }
-            }
-            psFree(job);
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            return NULL;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                psScalar *scalar = NULL;
+                scalar = job->args->data[5];
+                Nfit += scalar->data.S32;
+                scalar = job->args->data[6];
+                Npsf += scalar->data.S32;
+                scalar = job->args->data[7];
+                Next += scalar->data.S32;
+                scalar = job->args->data[8];
+                Nfail += scalar->data.S32;
+
+                // add these back onto sources
+                psArray *newSources = job->args->data[4];
+                for (int j = 0; j < newSources->n; j++) {
+                    psArrayAdd (sources, 16, newSources->data[j]);
+                }
+            }
+            psFree(job);
             }
     }
@@ -278 +276 @@
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
                 Next ++;
-                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
+                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
                 continue;
             }
@@ -286 +284 @@
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
                 Npsf ++;
-                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
+                source->mode |= PM_SOURCE_MODE_NONLINEAR_FIT;
                 continue;
             }

branches/pap/psphot/src/psphotFitSourcesLinear.c

@@ -27 +27 @@
     for (int i = 0; i < num; i++) {
 
-        // find the currently selected readout
-        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
-        psAssert (file, "missing file?");
-
-        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
-        psAssert (readout, "missing readout?");
-
-        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
-        psAssert (detections, "missing detections?");
-
-        psArray *sources = detections->allSources;
-        psAssert (sources, "missing sources?");
-
-        pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
-        psAssert (psf, "missing psf?");
-
-        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) {
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
+        psAssert (file, "missing file?");
+
+        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+        psAssert (readout, "missing readout?");
+
+        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+        psAssert (detections, "missing detections?");
+
+        psArray *sources = detections->allSources;
+        psAssert (sources, "missing sources?");
+
+        pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+        psAssert (psf, "missing psf?");
+
+        if (!psphotFitSourcesLinearReadout (recipe, readout, sources, psf, final)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for %s entry %d", filerule, i);
-            return false;
-        }
+            return false;
+        }
     }
     return true;
@@ -59 +59 @@
 
     if (!sources->n) {
-        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit");
-        return true;
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit");
+        return true;
     }
 
@@ -150 +150 @@
 
     if (fitSources->n == 0) {
+        psFree(fitSources);
         return true;
     }
@@ -170 +171 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
+        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
         psSparseMatrixElement (sparse, i, i, f);
 
         // the formal error depends on the weighting scheme
         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
-            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor);
+            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor, maskVal);
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
@@ -183 +184 @@
 
         // find the image x model value
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
+        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
         psSparseVectorElement (sparse, i, f);
 
@@ -189 +190 @@
         switch (SKY_FIT_ORDER) {
           case 1:
-            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
+            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 1, f);
-            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
+            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 2, f);
 
           case 0:
-            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
+            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
             psSparseBorderElementB (border, i, 0, f);
             break;
@@ -215 +216 @@
 
             // got an overlap; calculate cross-product and add to output array
-            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
+            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor, maskVal);
             psSparseMatrixElement (sparse, j, i, f);
         }

branches/pap/psphot/src/psphotFitSourcesLinearStack.c

@@ -43 +43 @@
     for (int i = 0; i < objects->n; i++) {
         pmPhotObj *object = objects->data[i];
-        if (!object) continue;
-        if (!object->sources) continue;
+        if (!object) continue;
+        if (!object->sources) continue;
 
-        // XXX check an element of the group to see if we should use it
-        // if (!object->flags & PM_PHOT_OBJ_BAD) continue;
+        // XXX check an element of the group to see if we should use it
+        // if (!object->flags & PM_PHOT_OBJ_BAD) continue;
 
-        for (int j = 0; j < object->sources->n; j++) {
-          pmSource *source = object->sources->data[j];
-          if (!source) continue;
+        for (int j = 0; j < object->sources->n; j++) {
+          pmSource *source = object->sources->data[j];
+          if (!source) continue;
 
-          // turn this bit off and turn it on again if we keep this source
-          source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
+          // turn this bit off and turn it on again if we keep this source
+          source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
 
-          // generate model for sources without, or skip if we can't
-          if (!source->modelFlux) {
+          // generate model for sources without, or skip if we can't
+          if (!source->modelFlux) {
             if (!pmSourceCacheModel (source, maskVal)) continue;
-          }
+          }
 
-          source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
-          psArrayAdd (fitSources, 100, source);
-        }
+          source->mode |= PM_SOURCE_MODE_LINEAR_FIT;
+          psArrayAdd (fitSources, 100, source);
+        }
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n);
@@ -85 +85 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
+        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal);
         psSparseMatrixElement (sparse, i, i, f);
 
         // the formal error depends on the weighting scheme
         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
-            float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR);
+            float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR, maskVal);
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
@@ -97 +97 @@
 
         // find the image x model value
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
+        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal);
         psSparseVectorElement (sparse, i, f);
 
@@ -104 +104 @@
             pmSource *SRCj = fitSources->data[j];
 
-            // we only need to generate dot terms for source on the same image
-            if (SRCj->imageID != SRCi->imageID) { continue; }
+            // we only need to generate dot terms for source on the same image
+            if (SRCj->imageID != SRCi->imageID) { continue; }
 
             // skip over disjoint source images, break after last possible overlap
@@ -114 +114 @@
 
             // got an overlap; calculate cross-product and add to output array
-            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
+            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR, maskVal);
             psSparseMatrixElement (sparse, j, i, f);
         }

branches/pap/psphot/src/psphotGuessModels.c

@@ -21 +21 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (i == chisqNum) continue; // skip chisq image
+        if (i == chisqNum) continue; // skip chisq image
         if (!psphotGuessModelsReadout (config, view, filerule, i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on to guess models for %s entry %d", filerule, i);
@@ -106 +106 @@
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                psFree (job);
                 return false;
             }
-            psFree(job);
         }
 

branches/pap/psphot/src/psphotImageLoop.c

@@ -46 +46 @@
         if (!psphotMosaicChip(config, view, "PSPHOT.INPUT", "PSPHOT.LOAD")) ESCAPE ("Unable to mosaic chip.");
 
+        // Read WCS if easy.
+        // XXX Since we're mosaicking cells, we ignore the case where the WCS is defined for a cell.
+        {
+            pmChip *inChip = pmFPAviewThisChip(view, input->fpa); // Mosaicked chip
+            pmHDU *hduLow = pmHDUGetLowest(input->fpa, inChip, NULL);
+            if (hduLow && !pmAstromReadWCS(input->fpa, inChip, hduLow->header, 1.0)) {
+                psWarning("Unable to read WCS astrometry from header.");
+                psErrorClear();
+                pmHDU *hduHigh = pmHDUGetHighest(input->fpa, inChip, NULL);
+                if (hduHigh && hduHigh != hduLow &&
+                    !pmAstromReadWCS(input->fpa, chip, hduHigh->header, 1.0)) {
+                    psWarning("Unable to read WCS astrometry from primary header.");
+                    psErrorClear();
+                }
+            }
+        }
+
         // try to load other supporting data (PSF, SRC, etc).
         // do not re-load the following three files
@@ -67 +84 @@
 
                 // Update the header
-                pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
-                if (hdu && hdu != lastHDU) {
-                    psphotVersionHeaderFull(hdu->header);
-                    lastHDU = hdu;
+                {
+                    pmHDU *hdu = pmHDUGetHighest(input->fpa, chip, cell);
+                    if (hdu && hdu != lastHDU) {
+                        psphotVersionHeaderFull(hdu->header);
+                        lastHDU = hdu;
+                    }
                 }
 
-                // if an external mask is supplied, ensure that NAN pixels are also masked
-                if (readout->mask) {
-                    psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
-                    if (!pmReadoutMaskNonfinite(readout, maskSat)) {
-                        psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");
-                        psFree(view);
-                        return false;
-                    }
-                }
+                // if an external mask is supplied, ensure that NAN pixels are also masked
+                if (readout->mask) {
+                    psImageMaskType maskSat = pmConfigMaskGet("SAT", config); // Mask value for saturated pixels
+                    if (!pmReadoutMaskNonfinite(readout, maskSat)) {
+                        psError(psErrorCodeLast(), false, "Unable to mask non-finite pixels.");
+                        psFree(view);
+                        return false;
+                    }
+                }
 
                 // run the actual photometry analysis on this chip/cell/readout
@@ -91 +110 @@
             }
 
-            // drop all versions of the internal files
-            status = true;
-            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
-            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");
-            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");
-            if (!status) {
-                psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files");
-                psFree (view);
-                return false;
-            }
+            // drop all versions of the internal files
+            status = true;
+            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL");
+            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKMDL.STDEV");
+            status &= pmFPAfileDropInternal (config->files, "PSPHOT.BACKGND");
+            if (!status) {
+                psError(PSPHOT_ERR_PROG, false, "trouble dropping internal files");
+                psFree (view);
+                return false;
+            }
         }
         // save output which is saved at the chip level

branches/pap/psphot/src/psphotMagnitudes.c

@@ -18 +18 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (i == chisqNum) continue; // skip chisq image
-
-        // find the currently selected readout
-        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
-        psAssert (file, "missing file?");
-
-        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
-        psAssert (readout, "missing readout?");
-
-        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
-        psAssert (detections, "missing detections?");
-
-        psArray *sources = detections->allSources;
-        psAssert (sources, "missing sources?");
-
-        pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
-        psAssert (psf, "missing psf?");
-
-        if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) {
+        if (i == chisqNum) continue; // skip chisq image
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, i); // File of interest
+        psAssert (file, "missing file?");
+
+        pmReadout *readout = pmFPAviewThisReadout(view, file->fpa);
+        psAssert (readout, "missing readout?");
+
+        pmDetections *detections = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.DETECTIONS");
+        psAssert (detections, "missing detections?");
+
+        psArray *sources = detections->allSources;
+        psAssert (sources, "missing sources?");
+
+        pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+        psAssert (psf, "missing psf?");
+
+        if (!psphotMagnitudesReadout (config, recipe, view, readout, sources, psf)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
-            return false;
-        }
+            return false;
+        }
     }
     return true;
@@ -50 +50 @@
 
     if (!sources->n) {
-        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source magnitudes");
-        return true;
-    }
-        
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source magnitudes");
+        return true;
+    }
+
     psTimerStart ("psphot.mags");
 
@@ -122 +122 @@
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                psFree (job);
                 return false;
             }
-            psFree(job);
 
 # if (0)
@@ -184 +182 @@
         }
 
-        // clear the mask bit and set the circular mask pixels
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
+        // clear the mask bit and set the circular mask pixels
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, source->apRadius, "OR", markVal);
 
         status = pmSourceMagnitudes (source, psf, photMode, maskVal); // maskVal includes markVal
         if (status && isfinite(source->apMag)) Nap ++;
 
-        // clear the mask bit 
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+        // clear the mask bit
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
 
         // re-subtract the object, leave local sky
@@ -273 +271 @@
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                psFree (job);
                 return false;
             }
-            psFree(job);
-
         }
 

branches/pap/psphot/src/psphotRadiusChecks.c

@@ -4 +4 @@
 static float PSF_FIT_NSIGMA;
 static float PSF_FIT_PADDING;
-static float PSF_APERTURE = 0;  // radius to use in PSF aperture mags
-static float PSF_FIT_RADIUS = 0;        // radius to use in fitting (ignored if <= 0,
-                                        // and a per-object radius is calculated)
+static float PSF_APERTURE = 0;  // radius to use in PSF aperture mags
+static float PSF_FIT_RADIUS = 0;        // radius to use in fitting (ignored if <= 0,
+                                        // and a per-object radius is calculated)
 
 bool psphotInitRadiusPSF(const psMetadata *recipe, const psMetadata *analysis, const pmModelType type) {
@@ -17 +17 @@
     PSF_FIT_RADIUS =  psMetadataLookupF32(&status, analysis, "PSF_FIT_RADIUS");
     if (!status) {
-        PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
+        PSF_FIT_RADIUS = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS");
     }
 
     PSF_APERTURE =  psMetadataLookupF32(&status, analysis, "PSF_APERTURE");
     if (!status) {
-        PSF_APERTURE =  psMetadataLookupF32(&status, recipe, "PSF_APERTURE");
+        PSF_APERTURE =  psMetadataLookupF32(&status, recipe, "PSF_APERTURE");
+    }
+
+    // The PSF_FIT_RADIUS and PSF_APERTURE may not be set if the PSF was loaded and not chosen
+
+    if (PSF_FIT_RADIUS == 0.0) {
+        float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA");
+        if (!status) {
+            gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
+        }
+        float fitScale = psMetadataLookupF32(&status, recipe, "PSF_FIT_RADIUS_SCALE");
+        PSF_FIT_RADIUS = (int)(fitScale*gaussSigma);
+    }
+
+    if (PSF_APERTURE == 0.0) {
+        float gaussSigma = psMetadataLookupF32(&status, analysis, "MOMENTS_GAUSS_SIGMA");
+        if (!status) {
+            gaussSigma = psMetadataLookupF32(&status, recipe, "MOMENTS_GAUSS_SIGMA");
+        }
+        float apScale = psMetadataLookupF32(&status, recipe, "PSF_APERTURE_SCALE");
+        PSF_APERTURE = (int)(apScale*gaussSigma);
     }
 
@@ -38 +58 @@
     // set the fit radius based on the object flux limit and the model
     float radiusFit = PSF_FIT_RADIUS;
-    if (radiusFit <= 0) {               // use fixed radius
-        if (moments == NULL) {
-            radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);
-        } else {
-            radiusFit = model->modelRadius(model->params, 1.0);
-        }
-        model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
+    if (radiusFit <= 0) {               // use fixed radius
+        if (moments == NULL) {
+            radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);
+        } else {
+            radiusFit = model->modelRadius(model->params, 1.0);
+        }
+        model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
     } else {
-        model->fitRadius = radiusFit;
+        model->fitRadius = radiusFit;
     }
     if (isnan(model->fitRadius)) psAbort("error in radius");
-        
+
     if (source->mode & PM_SOURCE_MODE_SATSTAR) {
-        model->fitRadius *= 2;
+        model->fitRadius *= 2;
     }
 
@@ -73 +93 @@
     // set the fit radius based on the object flux limit and the model
     float radiusFit = PSF_FIT_RADIUS;
-    if (radiusFit <= 0) {               // use fixed radius
-        if (moments == NULL) {
-            radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);
-        } else {
-            radiusFit = model->modelRadius(model->params, 1.0);
-        }
-        model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
+    if (radiusFit <= 0) {               // use fixed radius
+        if (moments == NULL) {
+            radiusFit = model->modelRadius(model->params, PSF_FIT_NSIGMA*moments->dSky);
+        } else {
+            radiusFit = model->modelRadius(model->params, 1.0);
+        }
+        model->fitRadius = (RADIUS_TYPE)(radiusFit + PSF_FIT_PADDING);
     } else {
-        model->fitRadius = radiusFit;
+        model->fitRadius = radiusFit;
     }
     if (isnan(model->fitRadius)) psAbort("error in radius");
@@ -89 +109 @@
 
     if (source->mode &  PM_SOURCE_MODE_SATSTAR) {
-        model->fitRadius *= 2;
+        model->fitRadius *= 2;
     }
 
@@ -134 +154 @@
     float radius = 0.0;
     for (int j = 0; j < footprint->spans->n; j++) {
-        pmSpan *span = footprint->spans->data[j];
-
-        float dY  = span->y  - peak->yf;
-        float dX0 = span->x0 - peak->xf;
-        float dX1 = span->x1 - peak->xf;
-
-        radius = PS_MAX (radius, hypot(dY, dX0));
-        radius = PS_MAX (radius, hypot(dY, dX1));
+        pmSpan *span = footprint->spans->data[j];
+
+        float dY  = span->y  - peak->yf;
+        float dX0 = span->x0 - peak->xf;
+        float dX1 = span->x1 - peak->xf;
+
+        radius = PS_MAX (radius, hypot(dY, dX0));
+        radius = PS_MAX (radius, hypot(dY, dX1));
     }
 

branches/pap/psphot/src/psphotReadout.c

@@ -45 +45 @@
         return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
     }
+
     if (!psphotSubtractBackground (config, view, "PSPHOT.INPUT")) {
         return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
@@ -93 +94 @@
         return psphotReadoutCleanup (config, view, "PSPHOT.INPUT");
     }
+
     // if we were not supplied a PSF model, determine the IQ stats here (detections->newSources)
     if (!psphotImageQuality (config, view, "PSPHOT.INPUT")) { // pass 1
@@ -208 +210 @@
 
     // calculate source magnitudes
-    psphotMagnitudes(config, view, "PSPHOT.INPUT"); // pass 1 (detections->allSources)
-
+    if (!psphotMagnitudes(config, view, "PSPHOT.INPUT")) { // pass 1 (detections->allSources)
+      psErrorStackPrint(stderr, "Unable to do magnitudes.");
+        psErrorClear();
+    }
     if (!psphotEfficiency(config, view, "PSPHOT.INPUT")) { // pass 1
         psErrorStackPrint(stderr, "Unable to determine detection efficiencies from fake sources");
@@ -219 +223 @@
 
     // replace background in residual image
-    psphotSkyReplace (config, view, "PSPHOT.INPUT"); // pass 1
-
+    if (!psphotSkyReplace (config, view, "PSPHOT.INPUT")) { // pass 1
+      psErrorStackPrint(stderr, "Unable to replace sky");
+      psErrorClear();
+
+/*       psLogMsg("psphot", 3, "failed on psphotSkyReplace"); */
+/*       return(psphotReadoutCleanup(config, view, "PSPHOT.INPUT")); */
+    }
     // drop the references to the image pixels held by each source
-    psphotSourceFreePixels (config, view, "PSPHOT.INPUT"); // pass 1
-
+    if (!psphotSourceFreePixels (config, view, "PSPHOT.INPUT")) { // pass 1
+      psErrorStackPrint(stderr, "Unable to free source pixels");
+      psErrorClear();
+
+/*       psLogMsg ("psphot", 3, "failed on psphotSourceFreePixels"); */
+/*       return(psphotReadoutCleanup(config, view, "PSPHOT.INPUT")); */
+    }
     // create the exported-metadata and free local data
     return psphotReadoutCleanup(config, view, "PSPHOT.INPUT");

branches/pap/psphot/src/psphotRoughClass.c

@@ -25 +25 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (i == chisqNum) continue; // skip chisq image
-        if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) {
+        if (i == chisqNum) continue; // skip chisq image
+        if (!psphotRoughClassReadout (config, view, filerule, i, recipe)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on rough classification for %s entry %d", filerule, i);
-            return false;
-        }
+            return false;
+        }
     }
     return true;
@@ -50 +50 @@
     bool havePSF = false;
     if (psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF")) {
-        havePSF = true;
+        havePSF = true;
     }
 
@@ -60 +60 @@
 
     if (!sources->n) {
-        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification");
-        return true;
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping rough classification");
+        return true;
     }
 
@@ -78 +78 @@
                 psLogMsg ("psphot", 4, "Failed to determine rough classification for region %f,%f - %f,%f\n",
                          region->x0, region->y0, region->x1, region->y1);
+
+                // If in doubt, it's a PSF
+                for (int i = 0; i < sources->n; i++) {
+                    pmSource *source = sources->data[i]; // Source of interest
+                    if (!source || !source->peak) {
+                        continue;
+                    }
+                    if (source->peak->x <  region->x0) continue;
+                    if (source->peak->x >= region->x1) continue;
+                    if (source->peak->y <  region->y0) continue;
+                    if (source->peak->y >= region->y1) continue;
+                    source->type = PM_SOURCE_TYPE_STAR;
+                }
                 psFree (region);
                 continue;
@@ -124 +137 @@
         // XXX why not save the psfClump as a PTR?
 
-        float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
-        float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");
+        float PSF_SN_LIM = psMetadataLookupF32(&status, recipe, "PSF_SN_LIM"); psAssert (status, "missing PSF_SN_LIM");
+        float MOMENTS_AR_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_AR_MAX"); psAssert (status, "missing MOMENTS_AR_MAX");
 
-        float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE");
-        if (!status) {
-            PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");
-            psAssert (status, "missing PSF_CLUMP_GRID_SCALE");
-        }
-        float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX");
-        if (!status) {
-            MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX");
-            psAssert (status, "missing MOMENTS_SX_MAX");
-        }
-        float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX");
-        if (!status) {
-            MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX");
-            psAssert (status, "missing MOMENTS_SY_MAX");
-        }
+        float PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, analysis, "PSF_CLUMP_GRID_SCALE");
+        if (!status) {
+            PSF_CLUMP_GRID_SCALE = psMetadataLookupF32(&status, recipe, "PSF_CLUMP_GRID_SCALE");
+            psAssert (status, "missing PSF_CLUMP_GRID_SCALE");
+        }
+        float MOMENTS_SX_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SX_MAX");
+        if (!status) {
+            MOMENTS_SX_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SX_MAX");
+            psAssert (status, "missing MOMENTS_SX_MAX");
+        }
+        float MOMENTS_SY_MAX = psMetadataLookupF32(&status, analysis, "MOMENTS_SY_MAX");
+        if (!status) {
+            MOMENTS_SY_MAX = psMetadataLookupF32(&status, recipe, "MOMENTS_SY_MAX");
+            psAssert (status, "missing MOMENTS_SY_MAX");
+        }
 
         psfClump = pmSourcePSFClump (NULL, region, sources, PSF_SN_LIM, PSF_CLUMP_GRID_SCALE, MOMENTS_SX_MAX, MOMENTS_SY_MAX, MOMENTS_AR_MAX);

branches/pap/psphot/src/psphotSourceStats.c

@@ -3 +3 @@
 // convert detections to sources and measure their basic properties (moments, local sky, sky
 // variance) Note: this function only generates sources for the new peaks (peak->assigned).
-// The new sources are added to any existing sources on detections->newSources.  The sources 
+// The new sources are added to any existing sources on detections->newSources.  The sources
 // on detections->allSources are ignored.
 bool psphotSourceStats (pmConfig *config, const pmFPAview *view, const char *filerule, bool setWindow)
@@ -22 +22 @@
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        // if (i == chisqNum) continue; // skip chisq image
+        // if (i == chisqNum) continue; // skip chisq image
         if (!psphotSourceStatsReadout (config, view, filerule, i, recipe, setWindow)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to find initial detections for %s entry %d", filerule, i);
@@ -91 +91 @@
     // generate the array of sources, define the associated pixel
     if (!detections->newSources) {
-        detections->newSources = psArrayAllocEmpty (peaks->n);
+        detections->newSources = psArrayAllocEmpty (peaks->n);
     }
     sources = detections->newSources;
@@ -107 +107 @@
         // create a new source
         pmSource *source = pmSourceAlloc();
-        source->imageID = index;
+        source->imageID = index;
 
         // add the peak
@@ -184 +184 @@
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
-                psFree (job);
                 psFree(detections->newSources);
                 return false;
             }
-            psFree(job);
         }
 
@@ -194 +192 @@
         if (!psThreadPoolWait (false)) {
             psError(PS_ERR_UNKNOWN, false, "Failure in thread job PSPHOT_SOURCE_STATS");
-            psFree(detections->newSources);
+            psFree(detections->newSources);
             return false;
         }
@@ -315 +313 @@
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to launch thread job PSPHOT_SOURCE_STATS");
-                psFree (job);
                 return NULL;
             }
-            psFree(job);
         }
 
@@ -380 +376 @@
         pmSource *source = sources->data[i];
 
-        if (source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue;
-        source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
+        if (source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED) continue;
+        source->tmpFlags |= PM_SOURCE_TMPF_MOMENTS_MEASURED;
 
         // skip faint sources for moments measurement
@@ -493 +489 @@
         psLogMsg ("psphot", 3, "radius %.1f, nStars: %d, nSigma: %5.2f, X,  Y: %f, %f (%f, %f)\n", sigma[i], psfClump.nStars, psfClump.nSigma, psfClump.X, psfClump.Y, sqrt(psfClump.X) / sigma[i], sqrt(psfClump.Y) / sigma[i]);
 
+#if 0
+        // Modifying clump parameters without restoring!
         psMetadataAddS32 (analysis, PS_LIST_TAIL, "PSF.CLUMP.NREGIONS",  PS_META_REPLACE, "psf clump regions", 1);
         psMetadata *regionMD = psMetadataLookupPtr (&status, analysis, "PSF.CLUMP.REGION.000");
@@ -504 +502 @@
         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DX", PS_META_REPLACE, "psf clump center", psfClump.dX);
         psMetadataAddF32 (regionMD, PS_LIST_TAIL, "PSF.CLUMP.DY", PS_META_REPLACE, "psf clump center", psfClump.dY);
-
         psphotVisualPlotMoments (recipe, analysis, sources);
+#endif
 
         Sout[i] = sqrt(0.5*(psfClump.X + psfClump.Y)) / sigma[i];

branches/pap/psphot/src/psphotTest.c

@@ -25 +25 @@
 
 bool FillImage_Threaded (psThreadJob *job) {
-    
+
     psImage *image = job->args->data[0];
     int xs = PS_SCALAR_VALUE(job->args->data[1],S32);
@@ -36 +36 @@
     psRegion region = psRegionSet (xs, xs + dx, ys, ys + dy);
     for (int i = 0; i < 100; i++) {
-        psImage *subset = psImageSubset (image, region);
-        psImageInit (subset, value + i);
-        psFree (subset);
+        psImage *subset = psImageSubset (image, region);
+        psImageInit (subset, value + i);
+        psFree (subset);
     }
     return true;
@@ -46 +46 @@
 
     if (argc != 3) {
-        fprintf (stderr, "USAGE: psphotTest (output.fits) (nThreads)\n");
-        exit (2);
+        fprintf (stderr, "USAGE: psphotTest (output.fits) (nThreads)\n");
+        exit (2);
     }
 
@@ -62 +62 @@
 
     for (int ix = 0; ix < 1000; ix += 100) {
-        for (int iy = 0; iy < 1000; iy += 100) {
+        for (int iy = 0; iy < 1000; iy += 100) {
 
-            // allocate a job -- if threads are not defined, this just runs the job
-            psThreadJob *job = psThreadJobAlloc ("FILL_IMAGE");
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("FILL_IMAGE");
 
-            psArrayAdd(job->args, 1, image);
-            PS_ARRAY_ADD_SCALAR(job->args, ix, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, iy, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, ix + iy, PS_TYPE_S32);
+            psArrayAdd(job->args, 1, image);
+            PS_ARRAY_ADD_SCALAR(job->args, ix, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, iy, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, 100, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, ix + iy, PS_TYPE_S32);
 
-            // FillImage (image, ix, iy, 100, 100, ix + iy);
+            // FillImage (image, ix, iy, 100, 100, ix + iy);
 
-            if (!psThreadJobAddPending(job)) {
-                fprintf (stderr, "failure to run FillImage(1)");
-                psFree (job);
-                exit (1);
-            }
-            psFree(job);
-        }
+            if (!psThreadJobAddPending(job)) {
+                fprintf (stderr, "failure to run FillImage(1)");
+                exit (1);
+            }
+        }
     }
 
@@ -88 +86 @@
     // wait for the threads to finish and manage results
     if (!psThreadPoolWait (true)) {
-        fprintf (stderr, "failure to run FillImage (2)");
-        exit (1);
+        fprintf (stderr, "failure to run FillImage (2)");
+        exit (1);
     }