Changeset 27547 for trunk/psphot/src/psphotFitSourcesLinearStack.c

Timestamp:

Mar 31, 2010, 7:02:48 PM (16 years ago)

Author:

eugene

Message:

minor work on psphotStack

File:

• trunk/psphot/src/psphotFitSourcesLinearStack.c (modified) (13 diffs)

trunk/psphot/src/psphotFitSourcesLinearStack.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
-// fit flux (and optionally sky model) to all reasonable sources
-// with the linear fitting process.  sources must have an associated
-// model with selected pixels, and the fit radius must be defined
-
-// given the set of sources, each of which points to the pixels in the
-// science image, we construct a set of simulated sources with their own pixels.
-// these are used to determine the simultaneous linear fit of fluxes.
-// the analysis is performed wrt the simulated pixel values
-
-bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, bool final) {
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotFitSourcesLinearStack (pmConfig *config, const pmFPAview *view, bool final)
+{
+    bool status = true;
+
+    int num = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.INPUT.NUM");
+    psAssert (status, "programming error: must define PSPHOT.INPUT.NUM");
+
+    // loop over the available readouts
+    for (int i = 0; i < num; i++) {
+        if (!psphotFitSourcesLinearReadoutStack (config, view, "PSPHOT.INPUT", i, final)) {
+            psError (PSPHOT_ERR_CONFIG, false, "failed to fit sources (linear) for PSPHOT.INPUT entry %d", i);
+            return false;
+        }
+    }
+    return true;
+}
+
+bool psphotFitSourcesLinearReadout (pmConfig *config, const pmFPAview *view, const char *filename, int index, bool final) {
 
     bool status;
@@ -18 +27 @@
     // float r;
 
+    // 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 linear fit");
+        return true;
+    }
+
+    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    psAssert (sources, "missing psf?");
+
     psTimerStart ("psphot.linear");
 
@@ -36 +70 @@
     // storage array for fitSources
     psArray *fitSources = psArrayAllocEmpty (sources->n);
+
+    // option to limit analysis to a specific region
+    char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
+    psRegion AnalysisRegion = psRegionFromString (region);
+    AnalysisRegion = psRegionForImage (readout->image, AnalysisRegion);
+    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
 
     bool CONSTANT_PHOTOMETRIC_WEIGHTS =
@@ -42 +82 @@
         psAbort("You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
     }
+    int SKY_FIT_ORDER = psMetadataLookupS32(&status, recipe, "SKY_FIT_ORDER");
+    if (!status) {
+        SKY_FIT_ORDER = 0;
+    }
+    bool SKY_FIT_LINEAR = psMetadataLookupBool(&status, recipe, "SKY_FIT_LINEAR");
+    if (!status) {
+        SKY_FIT_LINEAR = false;
+    }
+
+    // XXX test: choose a larger-than expected radius:
+    float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor);
+
+    // XXX do not apply covarFactor for the moment...
+    // covarFactor = 1.0;
 
     // select the sources which will be used for the fitting analysis
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
-
-        // turn this bit off and turn it on again if we pass this test
-        source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
-
-        // skip non-astronomical objects (very likely defects)
-        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
-
-        // do not include CRs in the full ensemble fit
-        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
-
-        if (final) {
-            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) continue;
-        } else {
-            if (source->mode & PM_SOURCE_MODE_BLEND) continue;
-        }
-
-        // generate model for sources without, or skip if we can't
-        if (!source->modelFlux) {
+    for (int i = 0; i < objects->n; i++) {
+        pmPhotObj *object = objects->data[i];
+        if (!object) continue;
+        if (!object->sources) continue;
+
+        // 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;
+
+          // 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) {
             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"), sources->n);
@@ -81 +133 @@
     psVector *errors = psVectorAlloc (fitSources->n, PS_TYPE_F32);
 
-    // create the border matrix (includes the sparse matrix)
-    // for just sky: 1 row; for x,y terms: 3 rows
+    // create the sparse matrix
     psSparse *sparse = psSparseAlloc (fitSources->n, 100);
-    psSparseBorder *border = psSparseBorderAlloc (sparse, 1);
 
     // fill out the sparse matrix elements and border elements (B)
@@ -93 +143 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
         psSparseMatrixElement (sparse, i, i, f);
 
         // the formal error depends on the weighting scheme
         if (CONSTANT_PHOTOMETRIC_WEIGHTS) {
-            float var = pmSourceModelDotModel (SRCi, SRCi, false);
+            float var = pmSourceModelDotModel (SRCi, SRCi, false, covarFactor);
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
@@ -104 +154 @@
         }
 
-
         // find the image x model value
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
         psSparseVectorElement (sparse, i, f);
-
-        f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
-        psSparseBorderElementB (border, i, 0, f);
 
         // loop over all other stars following this one
         for (int j = i + 1; j < fitSources->n; j++) {
             pmSource *SRCj = fitSources->data[j];
+
+            // XXX I need to know if this source is on the same image as SRCi --
+            if (!sameImge) { continue; }
 
             // skip over disjoint source images, break after last possible overlap
@@ -123 +172 @@
 
             // got an overlap; calculate cross-product and add to output array
-            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseMatrixElement (sparse, j, i, f);
         }
@@ -130 +179 @@
     psSparseResort (sparse);
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
-
-    // set the sky, sky_x, sky_y components of border matrix
-    SetBorderMatrixElements (border, readout, fitSources, CONSTANT_PHOTOMETRIC_WEIGHTS, SKY_FIT_ORDER, markVal);
 
     psSparseConstraint constraint;
@@ -140 +186 @@
 
     // solve for normalization terms (need include local sky?)
-    psVector *norm = psSparseSolve (NULL, constraint, sparse, 5);
-
+    psVector *norm = NULL;
+    norm = psSparseSolve (NULL, constraint, sparse, 5);
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
 
@@ -170 +216 @@
 
     // measure chisq for each source
-    for (int i = 0; final && (i < fitSources->n); i++) {
+    // for (int i = 0; final && (i < fitSources->n); i++) {
+    for (int i = 0; i < fitSources->n; i++) {
         pmSource *source = fitSources->data[i];
         if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue;
         pmModel *model = pmSourceGetModel (NULL, source);
-        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal);
+        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor);
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
@@ -182 +229 @@
     psFree (fitSources);
     psFree (norm);
+    psFree (skyfit);
     psFree (errors);
     psFree (border);
@@ -188 +236 @@
 
     psphotVisualShowResidualImage (readout);
-    psphotVisualShowFlags (sources);
+    psphotVisualPlotChisq (sources);
+    // psphotVisualShowFlags (sources);
+
+    // We have to place this visualization here because the models are not realized until
+    // psphotGuessModels or fitted until psphotFitSourcesLinear.
+    psphotVisualShowPSFStars (recipe, psf, sources);
 
     return true;
 }
+
+// XXX do we need this?
+// XXX disallow the simultaneous sky fit and remove this code...
+
+// Calculate the weight terms for the sky fit component of the matrix.  This function operates
+// on the pixels which correspond to all of the sources of interest.  These elements fill in
+// the border matrix components in the sparse matrix equation.
+static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal) {
+
+    // generate the image-wide weight terms
+    // turn on MARK for all image pixels
+    psRegion fullArray = psRegionSet (0, 0, 0, 0);
+    fullArray = psRegionForImage (readout->mask, fullArray);
+    psImageMaskRegion (readout->mask, fullArray, "OR", markVal);
+
+    // turn off MARK for all object pixels
+    for (int i = 0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        pmModel *model = pmSourceGetModel (NULL, source);
+        if (model == NULL) continue;
+        float x = model->params->data.F32[PM_PAR_XPOS];
+        float y = model->params->data.F32[PM_PAR_YPOS];
+        psImageMaskCircle (source->maskView, x, y, model->fitRadius, "AND", PS_NOT_IMAGE_MASK(markVal));
+    }
+
+    // accumulate the image statistics from the masked regions
+    psF32 **image  = readout->image->data.F32;
+    psF32 **variance = readout->variance->data.F32;
+    psImageMaskType  **mask   = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA;
+
+    double w, x, y, x2, xy, y2, xc, yc, wt, f, fo, fx, fy;
+    w = x = y = x2 = xy = y2 = fo = fx = fy = 0;
+
+    int col0 = readout->image->col0;
+    int row0 = readout->image->row0;
+
+    for (int j = 0; j < readout->image->numRows; j++) {
+        for (int i = 0; i < readout->image->numCols; i++) {
+            if (mask[j][i]) continue;
+            if (constant_weights) {
+                wt = 1.0;
+            } else {
+                wt = variance[j][i];
+            }
+            f = image[j][i];
+            w   += 1/wt;
+            fo  += f/wt;
+            if (SKY_FIT_ORDER == 0) continue;
+
+            xc  = i + col0;
+            yc  = j + row0;
+            x  +=    xc/wt;
+            y  +=    yc/wt;
+            x2 += xc*xc/wt;
+            xy += xc*yc/wt;
+            y2 += yc*yc/wt;
+            fx +=  f*xc/wt;
+            fy +=  f*yc/wt;
+        }
+    }
+
+    // turn off MARK for all image pixels
+    psImageMaskRegion (readout->mask, fullArray, "AND", PS_NOT_IMAGE_MASK(markVal));
+
+    // set the Border T elements
+    psSparseBorderElementG (border, 0, fo);
+    psSparseBorderElementT (border, 0, 0, w);
+    if (SKY_FIT_ORDER > 0) {
+        psSparseBorderElementG (border, 0, fx);
+        psSparseBorderElementG (border, 0, fy);
+        psSparseBorderElementT (border, 1, 0, x);
+        psSparseBorderElementT (border, 2, 0, y);
+        psSparseBorderElementT (border, 0, 1, x);
+        psSparseBorderElementT (border, 1, 1, x2);
+        psSparseBorderElementT (border, 2, 1, xy);
+        psSparseBorderElementT (border, 0, 2, y);
+        psSparseBorderElementT (border, 1, 2, xy);
+        psSparseBorderElementT (border, 2, 2, y2);
+    }
+
+    return true;
+}