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

Timestamp:

Apr 11, 2010, 5:08:29 PM (16 years ago)

Author:

eugene

Message:

updates to support psphotStack

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/psphotFitSourcesLinearStack.c (modified) (11 diffs)

trunk/psphot/src

@@ -21 +21 @@
 psphotForced
 psphotMakePSF
+psphotStack

@@ -21 +21 @@
 psphotForced
 psphotMakePSF
+psphotStack

trunk/psphot/src/psphotFitSourcesLinearStack.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+// XXX need array of covar factors for each image
+// XXX define the 'good' / 'bad' flags?
 
-// 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;
+# define COVAR_FACTOR 1.0
 
-    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 psphotFitSourcesLinearStack (pmConfig *config, psArray *objects, bool final) {
 
     bool status;
-    float x;
-    float y;
     float f;
-    // 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");
@@ -65 +27 @@
     maskVal |= markVal;
 
-    // source analysis is done in spatial order
-    sources = psArraySort (sources, pmSourceSortByY);
+    // analysis is done in spatial order (to speed up overlap search)
+    // sort by first element in each source list
+    objects = psArraySort (objects, pmPhotObjSortByX);
 
     // storage array for fitSources
-    psArray *fitSources = psArrayAllocEmpty (sources->n);
+    psArray *fitSources = psArrayAllocEmpty (objects->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 = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
+    psAssert (status, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
 
-    bool CONSTANT_PHOTOMETRIC_WEIGHTS =
-        psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
-    if (!status) {
-        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;
+    // XXX store a local static array of covar factors for each of the images (by image ID)
+    // float covarFactor = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+    // psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "covariance factor: %f\n", covarFactor);
 
     // select the sources which will be used for the fitting analysis
@@ -104 +47 @@
         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;
+        // 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++) {
@@ -123 +66 @@
         }
     }
-    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), sources->n);
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n);
 
     if (fitSources->n == 0) {
@@ -143 +86 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
+        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
         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, COVAR_FACTOR);
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
@@ -155 +98 @@
 
         // find the image x model value
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
+        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
         psSparseVectorElement (sparse, i, f);
 
@@ -162 +105 @@
             pmSource *SRCj = fitSources->data[j];
 
-            // XXX I need to know if this source is on the same image as SRCi --
-            if (!sameImge) { 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
-            if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) break;
-            if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue;
-            if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) continue;
-            if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue;
+            if (SRCj->pixels->row0 + SRCj->pixels->numRows < SRCi->pixels->row0) continue;  // source(i) is above source(j)
+            if (SRCi->pixels->row0 + SRCi->pixels->numRows < SRCj->pixels->row0) continue;  // source(i) is below source(j)
+            if (SRCj->pixels->col0 + SRCj->pixels->numCols < SRCi->pixels->col0) continue;  // source(i) is right of source(j)
+            if (SRCi->pixels->col0 + SRCi->pixels->numCols < SRCj->pixels->col0) break;     // source(i) is left of source(j) [no other source(j) can overlap source(i)]
 
             // 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, COVAR_FACTOR);
             psSparseMatrixElement (sparse, j, i, f);
         }
@@ -190 +133 @@
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "solve matrix: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
 
-    // XXXX **** philosophical question:
-    // we measure bright objects in three passes: 1) linear fit; 2) non-linear fit; 3) linear fit:
-    // should retain the chisq and errors from the intermediate non-linear fit?
-    // the non-linear fit provides better values for the position errors, and for
-    // extended sources, the shape errors
-
     // adjust I0 for fitSources and subtract
     for (int i = 0; i < fitSources->n; i++) {
@@ -208 +145 @@
         model->params->data.F32[PM_PAR_I0] = norm->data.F32[i];
         model->dparams->data.F32[PM_PAR_I0] = errors->data.F32[i];
-        // XXX is the value of 'errors' modified by the sky fit?
 
         // subtract object
@@ -221 +157 @@
         if (source->mode & PM_SOURCE_MODE_NONLINEAR_FIT) continue;
         pmModel *model = pmSourceGetModel (NULL, source);
-        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, covarFactor);
+        pmSourceChisq (model, source->pixels, source->maskObj, source->variance, maskVal, COVAR_FACTOR);
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
@@ -229 +165 @@
     psFree (fitSources);
     psFree (norm);
-    psFree (skyfit);
     psFree (errors);
-    psFree (border);
 
     psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
-
-    psphotVisualShowResidualImage (readout);
-    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...
+// sort by X (ascending)
+int pmPhotObjSortByX (const void **a, const void **b)
+{
+    pmPhotObj *objA = *(pmPhotObj **)a;
+    pmPhotObj *objB = *(pmPhotObj **)b;
 
-// 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) {
+    psF32 fA = objA->x;
+    psF32 fB = objB->x;
 
-    // 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;
+    psF32 diff = fA - fB;
+    if (diff > FLT_EPSILON) return (+1);
+    if (diff < FLT_EPSILON) return (-1);
+    return (0);
 }

@@ -21 +21 @@
 psphotForced
 psphotMakePSF
+psphotStack