Changeset 27839 for branches/simmosaic_branches/psphot/src/psphotApResid.c

Timestamp:

May 3, 2010, 8:45:22 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simmosaic_branches

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotApResid.c (modified) (17 diffs)

branches/simmosaic_branches/psphot/src

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

branches/simmosaic_branches/psphot/src/psphotApResid.c

@@ -4 +4 @@
 // measure the aperture residual statistics and 2D variations
 
-bool psphotApResid (pmConfig *config, pmReadout *readout, psArray *sources, pmPSF *psf)
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotApResid (pmConfig *config, const pmFPAview *view)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    psAssert (recipe, "missing recipe?");
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+        if (!psphotApResidReadout (config, view, "PSPHOT.INPUT", i, recipe)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to measure aperture residual for PSPHOT.INPUT entry %d", i);
+            return false;
+        }
+    }
+    return true;
+}
+
+bool psphotApResidReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, psMetadata *recipe)
 {
     int Nfail = 0;
@@ -13 +40 @@
     pmSource *source;
 
-    PS_ASSERT_PTR_NON_NULL(config, false);
-    PS_ASSERT_PTR_NON_NULL(readout, false);
-    PS_ASSERT_PTR_NON_NULL(sources, false);
-    PS_ASSERT_PTR_NON_NULL(psf, false);
-
     psTimerStart ("psphot.apresid");
 
-    // select the appropriate recipe information
-    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
-    assert (recipe);
+    // find the currently selected readout
+    pmFPAfile *file = pmFPAfileSelectSingle(config->files, filename, index); // 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?");
+
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping ap resid");
+        return true;
+    }
+
+    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    psAssert (psf, "missing psf?");
 
     // determine the number of allowed threads
@@ -33 +72 @@
     if (!measureAptrend) {
         // save nan values since these were not calculated
-        psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   NAN);
-        psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   NAN);
-        psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
-        psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
-        psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
-        psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
+        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
+        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
+        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
+        psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
         return true;
     }
@@ -53 +90 @@
     maskVal |= markVal;
 
-    // S/N limit to perform full non-linear fits
+    // clipping for extreme outliers
+    // XXX this is not currently defined in the recipe
     float MAX_AP_OFFSET = psMetadataLookupF32 (&status, recipe, "MAX_AP_OFFSET");
 
-    // this is the smallest radius allowed: need to at least extend growth curve down to this...
+    // options for how the photometry is calculated
+    // XXX are these sensible?
     bool IGNORE_GROWTH = psMetadataLookupBool (&status, recipe, "IGNORE_GROWTH");
     bool INTERPOLATE_AP = psMetadataLookupBool (&status, recipe, "INTERPOLATE_AP");
@@ -100 +139 @@
             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
@@ -110 +150 @@
             }
             psFree(job);
-
-# if (0)
-                int nskip = 0;
-                int nfail = 0;
-
-                if (!psphotApResidMags_Unthreaded (&nskip, &nfail, sources, psf, photMode, maskVal)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                    return false;
-                }
-                Nskip += nskip;
-                Nfail += nfail;
-# endif
-
         }
 
@@ -138 +165 @@
             } else {
                 psScalar *scalar = NULL;
-                scalar = job->args->data[4];
+                scalar = job->args->data[5];
                 Nskip += scalar->data.S32;
-                scalar = job->args->data[5];
+                scalar = job->args->data[6];
                 Nfail += scalar->data.S32;
             }
@@ -150 +177 @@
 
     // gather the stats to assess the aperture residuals
-    psVector *mask    = psVectorAllocEmpty (300, PS_TYPE_VECTOR_MASK);
     psVector *mag     = psVectorAllocEmpty (300, PS_TYPE_F32);
     psVector *xPos    = psVectorAllocEmpty (300, PS_TYPE_F32);
@@ -158 +184 @@
     Npsf = 0;
 
+# ifdef DEBUG    
+    FILE *f = fopen ("apresid.dat", "w");
+    psAssert (f, "failed open");
+# endif
+
     for (int i = 0; i < sources->n; i++) {
         source = sources->data[i];
@@ -168 +199 @@
         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 (!isfinite(source->apMag) || !isfinite(source->psfMag)) {
             continue;
         }
+
+        // XXX make this user-configurable?
+        if (source->errMag > 0.01) continue;
 
         // aperture residual for this source
@@ -182 +220 @@
         }
 
-        mag->data.F32[Npsf]     = source->psfMag;
-        apResid->data.F32[Npsf] = dap;
-        xPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_XPOS];
-        yPos->data.F32[Npsf]    = model->params->data.F32[PM_PAR_YPOS];
-
-        mask->data.PS_TYPE_VECTOR_MASK_DATA[Npsf] = 0;
-
-        dMag->data.F32[Npsf] = model->dparams->data.F32[PM_PAR_I0] / model->params->data.F32[PM_PAR_I0];
-
-        psVectorExtend (mag,     100, 1);
-        psVectorExtend (mask,    100, 1);
-        psVectorExtend (xPos,    100, 1);
-        psVectorExtend (yPos,    100, 1);
-        psVectorExtend (dMag,    100, 1);
-        psVectorExtend (apResid, 100, 1);
+# 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]);
+# 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");
+
+        psVectorAppend (mag, source->psfMag);
+        psVectorAppend (dMag,source->errMag);
+        psVectorAppend (apResid, dap);
+        psVectorAppend (xPos, model->params->data.F32[PM_PAR_XPOS]);
+        psVectorAppend (yPos, model->params->data.F32[PM_PAR_YPOS]);
         Npsf ++;
     }
@@ -202 +245 @@
     psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "measure aperture residuals for %d objects (%d skipped, %d failed, %ld invalid)\n",
               Npsf, Nskip, Nfail, sources->n - Npsf - Nskip - Nfail);
+
+# ifdef DEBUG
+    fclose (f);
+# endif
 
     // XXX choose a better value here?
@@ -209 +256 @@
     }
 
-    // XXX deprecating the old code which allowed the ApResid to be fitted as a function of flux and r^2/flux
-    // XXX is this asymmetric clipping still needed?  this analysis should come after neighbors are subtracted...
-    // 3hi/1lo sigma clipping on the rflux vs metric fit
-    // systematic error information
-    float errorScale = 0.0;
-    float errorFloor = 0.0;
-
+    // 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
+
+    // set the max order (0 = constant) which the number of psf stars can support:
+    // rule of thumb: require 3 stars per 'cell' (order+1)^2
+    int MaxOrderForStars = 0;
+    if (Npsf >= 12) MaxOrderForStars = 1; // 4 cells
+    if (Npsf >= 27) MaxOrderForStars = 2; // 9 cells
+    if (Npsf >= 48) MaxOrderForStars = 3; // 16 cells
+    if (Npsf >  75) MaxOrderForStars = 4; // 25 cells
+
+    pmTrend2DMode mode = PM_TREND_MAP;
+    if (mode == PM_TREND_MAP) {
+        MaxOrderForStars ++;
+    } 
+    APTREND_ORDER_MAX = PS_MIN (APTREND_ORDER_MAX, MaxOrderForStars);
+
+    psFree (psf->ApTrend);
+    psf->ApTrend = NULL;
     float errorFloorMin = FLT_MAX;
-    int entryMin = -1;
-
-    // Fit out the dap vs mag trend, iterate over spatial scale until error Floor increases.
-    // Stop if Npsf / (Nx * Ny) < 3
+
+    // as we loop over orders, we need to refer to the initial selection, but we modify the
+    // option values to match the current guess: save the max values here:
+    int NX = readout->image->numCols;
+    int NY = readout->image->numRows;
     for (int i = 1; i <= APTREND_ORDER_MAX; i++) {
-
-        if (!psphotApResidTrend (readout, psf, Npsf, i, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag)) {
-            break;
-        }
-
-        // store the resulting errorFloor values and the scales, redo the best
+        
+        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");
+
+        // store the minimum errorFloor and best ApTrend to keep
         if (errorFloor < errorFloorMin) {
             errorFloorMin = errorFloor;
-            entryMin = i;
+            psFree (psf->ApTrend);
+            psf->ApTrend = psMemIncrRefCounter(apTrend);
         }
-    }
-    if (entryMin == -1) {
+        psFree (apTrend);
+    }
+    if (psf->ApTrend == NULL) {
         psWarning("Failed to find a valid aperture residual value");
         goto escape;
     }
 
-    // XXX catch error condition
-    psphotApResidTrend (readout, psf, Npsf, entryMin, &errorScale, &errorFloor, mask, xPos, yPos, apResid, dMag);
+    // 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;
+
+    psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
+    psf->dApResid = errorFloorMin;
+    psf->nApResid = Npsf;
+
+    // save results for later output
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid);
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
+
+    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
+
+    psFree (xPos);
+    psFree (yPos);
+    psFree (apResid);
+    psFree (mag);
+    psFree (dMag);
+
+    psphotVisualPlotApResid (sources, psf->ApResid, psf->dApResid);
+
+    return true;
+
+escape:
+    // save nan values since these were not calculated
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
+    psMetadataAdd (readout->analysis, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
+
+    psFree (xPos);
+    psFree (yPos);
+    psFree (apResid);
+    psFree (mag);
+    psFree (dMag);
+    return true;
+    // this is a quality error, not a programming error
+}
+
+pmTrend2D *psphotApResidTrend (float *apResidSysErr, pmReadout *readout, int Nx, int Ny, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
+
+    // the mask marks the values not used to calculate the ApTrend
+    psVector *mask = psVectorAlloc(xPos->n, PS_TYPE_VECTOR_MASK);
+    psVectorInit (mask, 0);
+
+    // XXX allow user to set this optionally?
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN | PS_STAT_SAMPLE_STDEV);
+
+    // measure Trend2D for the current spatial scale
+    pmTrend2D *apTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
+
+    // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate:
+    // XXX use this or not?  probably not, since this is the point of the systematic error analysis
+    psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32);
+    for (int i = 0; i < dMag->n; i++) {
+        dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.005);
+    }
+
+    // XXX test for errors here
+    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;
+    }
+    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:
+    }
 
     // construct the fitted values and the residuals
-    psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos);
+    psVector *apResidFit = pmTrend2DEvalVector (apTrend, mask, 0xff, xPos, yPos);
     psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
-    psVector *dMagSys = (psVector *) psBinaryOp (NULL, (void *) dMag, "*", (void *) psScalarAlloc(errorScale, PS_TYPE_F32));
-
-    if (psTraceGetLevel("psphot") >= 2) {
-        FILE *dumpFile = fopen ("apresid.dat", "w");
+
+    // measure systematic error
+    *apResidSysErr = psVectorSystematicError (apResidRes, dMag, 0.10);
+    if (!isfinite(*apResidSysErr)) {
+        psWarning("Failed to find systematic error for %d x %d map", Nx, Ny);
+        psFree (apTrend);
+        return NULL;
+    }
+
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid\n", Nx, Ny);
+    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *apResidSysErr);
+
+    if (psTraceGetLevel("psphot") >= 4) {
+        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++) {
-            fprintf (dumpFile, "%f %f  %f %f %f  %f %f  %x\n",
+            fprintf (dumpFile, "%f %f  %f %f  %f %f  %x\n",
                      xPos->data.F32[i], yPos->data.F32[i],
-                     mag->data.F32[i], dMag->data.F32[i], dMagSys->data.F32[i],
+                     dMag->data.F32[i], hypot(dMag->data.F32[i], *apResidSysErr),
                      apResid->data.F32[i], apResidRes->data.F32[i],
                      mask->data.PS_TYPE_VECTOR_MASK_DATA[i]);
@@ -258 +421 @@
     }
 
-    // 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;
-
-    psf->ApResid  = pmTrend2DEval (psf->ApTrend, xc, yc); // ap-fit at chip center
-    psf->dApResid = errorFloor;
-    psf->nApResid = Npsf;
-
-    // save results for later output
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   0.0);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   0.0);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   psf->ApResid);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", psf->dApResid);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", psf->growth->apLoss);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", psf->nApResid);
-
-    psLogMsg ("psphot.apresid", PS_LOG_DETAIL, "aperture residual: %f +/- %f\n", psf->ApResid, psf->dApResid);
-    psLogMsg ("psphot.apresid", PS_LOG_INFO, "measure full-frame aperture residuals for %d sources: %f sec\n", Npsf, psTimerMark ("psphot.apresid"));
-
-    psFree (mag);
     psFree (mask);
-    psFree (xPos);
-    psFree (yPos);
-
-    psFree (apResid);
-    psFree (apResidFit);
-    psFree (apResidRes);
-
-    psFree (dMagSys);
-    psFree (dMag);
-
-    psphotVisualPlotApResid (sources);
-
-    return true;
-
-escape:
-    // save nan values since these were not calculated
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYBIAS",  PS_DATA_F32 | PS_META_REPLACE, "aperture sky bias",   NAN);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "SKYSAT",   PS_DATA_F32 | PS_META_REPLACE, "aperture-determined saturation",   NAN);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "APMIFIT",  PS_DATA_F32 | PS_META_REPLACE, "aperture residual",   NAN);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "DAPMIFIT", PS_DATA_F32 | PS_META_REPLACE, "ap residual scatter", NAN);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "APLOSS",   PS_DATA_F32 | PS_META_REPLACE, "aperture loss (mag)", NAN);
-    psMetadataAdd (recipe, PS_LIST_TAIL, "NAPMIFIT", PS_DATA_S32 | PS_META_REPLACE, "number of apresid stars", 0);
-
-    psFree (mag);
-    psFree (mask);
-    psFree (xPos);
-    psFree (yPos);
-    psFree (apResid);
-    psFree (dMag);
-    return false;
-}
-
-/*
-  (aprMag' - fitMag) = rflux*skyBias + ApTrend(x,y)
-  (aprMag - rflux*skyBias) - fitMag = ApTrend(x,y)
-  (aprMag - rflux*skyBias) = fitMag + ApTrend(x,y)
-*/
-
-bool psphotMagErrorScale (float *errorScale, float *errorFloor, psVector *dMag, psVector *dap, psVector *mask, int nGroup) {
-
-    psStats *statsS = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    psStats *statsM = psStatsAlloc (PS_STAT_SAMPLE_MEAN);
-
-    // measure the trend in bins with 10 values each; if < 10 total, use them all
-    int nBin = PS_MAX (dMag->n / nGroup, 1);
-
-    // output vectors for ApResid trend
-    psVector *dSo = psVectorAlloc (nBin, PS_TYPE_F32);
-    psVector *dMo = psVectorAlloc (nBin, PS_TYPE_F32);
-    psVector *dRo = psVectorAlloc (nBin, PS_TYPE_F32);
-
-    // use dMag to group the dMag and dap vectors
-    psVector *index = psVectorSortIndex (NULL, dMag);
-
-    // subset vectors for dMag and dap values within the given range
-    psVector *dMSubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
-    psVector *dASubset = psVectorAllocEmpty (nGroup, PS_TYPE_F32);
-    psVector *mkSubset = psVectorAllocEmpty (nGroup, PS_TYPE_VECTOR_MASK);
-
-    int n = 0;
-    for (int i = 0; i < dMo->n; i++) {
-        int j;
-        for (j = 0; (j < nGroup) && (n < dMag->n); j++, n++) {
-            int N = index->data.U32[n];
-            dMSubset->data.F32[j] = dMag->data.F32[N];
-            dASubset->data.F32[j] = dap->data.F32[N];
-            mkSubset->data.PS_TYPE_VECTOR_MASK_DATA[j] = mask->data.PS_TYPE_VECTOR_MASK_DATA[N];
-        }
-        dMSubset->n = j;
-        dASubset->n = j;
-        mkSubset->n = j;
-
-        psStatsInit (statsS);
-        psStatsInit (statsM);
-
-        if (j > 2) {
-            if (!psVectorStats (statsS, dASubset, NULL, mkSubset, 0xff)) {
-                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-                return false;
-            }
-            if (!psVectorStats (statsM, dMSubset, NULL, mkSubset, 0xff)) {
-                psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-                return false;
-            }
-            dSo->data.F32[i] = statsS->robustStdev;
-            dMo->data.F32[i] = statsM->sampleMean;
-            dRo->data.F32[i] = statsS->robustStdev / statsM->sampleMean;
-        } else {
-            dSo->data.F32[i] = NAN;
-            dMo->data.F32[i] = NAN;
-            dRo->data.F32[i] = NAN;
-        }
-    }
-    psFree (dMSubset);
-    psFree (dASubset);
-    psFree (mkSubset);
-
-    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
-    if (!psVectorStats (stats, dRo, NULL, NULL, 0)) {
-        psError(PS_ERR_UNKNOWN, false, "failure to measure stats");
-        return false;
-    }
-
-    *errorScale = stats->sampleMedian;
-    for (int i = 0; i < dSo->n; i++) {
-        *errorFloor = dSo->data.F32[i];
-        if (isfinite(*errorFloor)) break;
-    }
-
-    psFree (stats);
-    psFree (index);
-
-    psFree (dRo);
-    psFree (dMo);
-    psFree (dSo);
-
-    psFree (statsS);
-    psFree (statsM);
-
-    return true;
-}
-
-bool psphotApResidTrend (pmReadout *readout, pmPSF *psf, int Npsf, int scale, float *errorScale, float *errorFloor, psVector *mask, psVector *xPos, psVector *yPos, psVector *apResid, psVector *dMag) {
-
-    int Nx, Ny;
-
-    if (readout->image->numCols > readout->image->numRows) {
-        Nx = scale;
-        float AR = readout->image->numRows / (float) readout->image->numCols;
-        Ny = (int) (Nx * AR + 0.5);
-        Ny = PS_MAX (1, Ny);
-    } else {
-        Ny = scale;
-        float AR = readout->image->numRows / (float) readout->image->numCols;
-        Nx = (int) (Ny * AR + 0.5);
-        Nx = PS_MAX (1, Nx);
-    }
-
-    // require at least 10 stars per spatial bin
-    if (Npsf < 10*Nx*Ny) {
-        return false;
-    }
-
-    // the mask marks the values not used to calculate the ApTrend
-    psVectorInit (mask, 0);
-
-    // XXX stats structure for use by ApTrend : make parameters user setable
-    psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
-    stats->min = 2.0;
-    stats->max = 3.0;
-
-    // measure Trend2D for the current spatial scale
-    psFree (psf->ApTrend);
-    psf->ApTrend = pmTrend2DAlloc (PM_TREND_MAP, readout->image, Nx, Ny, stats);
-
-    // XXX somewhat arbitrary: soften the errors so the few bright stars do not totally dominate:
-    psVector *dMagSoft = psVectorAlloc (dMag->n, PS_TYPE_F32);
-    for (int i = 0; i < dMag->n; i++) {
-        dMagSoft->data.F32[i] = hypot(dMag->data.F32[i], 0.01);
-    }
-
-    // XXX test for errors here
-    pmTrend2DFit (psf->ApTrend, mask, 0xff, xPos, yPos, apResid, dMagSoft);
-
-    // construct the fitted values and the residuals
-    psVector *apResidFit = pmTrend2DEvalVector (psf->ApTrend, mask, 0xff, xPos, yPos);
-    psVector *apResidRes = (psVector *) psBinaryOp (NULL, (void *) apResid, "-", (void *) apResidFit);
-
-    // measure systematic errorFloor & systematic / photon scale factor
-    // XXX this is a bit arbitrary, but it forces ~3 stars from the bright bin per spatial bin
-    int nGroup = PS_MAX (3*Nx*Ny, 10);
-    psphotMagErrorScale (errorScale, errorFloor, dMag, apResidRes, mask, nGroup);
-
-    psLogMsg ("psphot.apresid", PS_LOG_INFO, "result of %d x %d grid (%d stars per bin)\n", Nx, Ny, nGroup);
-    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic error / photon error: %f\n", *errorScale);
-    psLogMsg ("psphot.apresid", PS_LOG_INFO, "systematic scatter floor: %f\n", *errorFloor);
-
     psFree (stats);
     psFree (dMagSoft);
@@ -461 +427 @@
     psFree (apResidRes);
 
-    return true;
+    return apTrend;
 }
 
@@ -475 +441 @@
     pmSourcePhotometryMode photMode = PS_SCALAR_VALUE(job->args->data[2],S32);
     psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
 
     for (int i = 0; i < sources->n; i++) {
@@ -485 +452 @@
         if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
 
-        if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
-            Nskip ++;
-            psTrace ("psphot", 3, "skip : bad source mag");
-            continue;
+        // replace object in image
+        if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+            pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
         }
 
-        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;
+        // 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;
     }
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    scalar = job->args->data[4];
+    scalar = job->args->data[5];
     scalar->data.S32 = Nskip;
 
-    scalar = job->args->data[5];
+    scalar = job->args->data[6];
     scalar->data.S32 = Nfail;
 
     return true;
 }
-
-# if (0)
-bool psphotApResidMags_Unthreaded (int *nskip, int *nfail, psArray *sources, pmPSF *psf, pmSourcePhotometryMode photMode, psImageMaskType maskVal) {
-
-    int Nskip = 0;
-    int Nfail = 0;
-
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = (pmSource *) sources->data[i];
-
-        if (source->type != PM_SOURCE_TYPE_STAR) SKIPSTAR ("NOT STAR");
-        if (source->mode &  PM_SOURCE_MODE_SATSTAR) SKIPSTAR ("SATSTAR");
-        if (source->mode &  PM_SOURCE_MODE_BLEND) SKIPSTAR ("BLEND");
-        if (source->mode &  PM_SOURCE_MODE_FAIL) SKIPSTAR ("FAIL STAR");
-        if (source->mode &  PM_SOURCE_MODE_POOR) SKIPSTAR ("POOR STAR");
-
-        if (!pmSourceMagnitudes (source, psf, photMode, maskVal)) {
-            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;
-        }
-    }
-
-    // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    *nskip = Nskip;
-    *nfail = Nfail;
-
-    return true;
-}
-# endif

@@ -18 +18 @@
 psphotVersionDefinitions.h
 psphotMomentsStudy
+psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack