Changeset 27840 for branches/simtest_nebulous_branches/psphot/src/psphotFitSourcesLinearStack.c

Timestamp:

May 3, 2010, 8:50:52 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/simtest_nebulous_branches

Files:

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

branches/simtest_nebulous_branches/psphot/src

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

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

branches/simtest_nebulous_branches/psphot/src/psphotFitSourcesLinearStack.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+// XXX need array of covar factors for each image
+// XXX define the 'good' / 'bad' flags?
 
-// 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
+# define COVAR_FACTOR 1.0
 
-// 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
-
-static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, psImageMaskType markVal);
-
-bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, psMetadata *recipe, pmPSF *psf, 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);
 
     psTimerStart ("psphot.linear");
@@ -33 +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);
 
-    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");
-    }
+    bool CONSTANT_PHOTOMETRIC_WEIGHTS = psMetadataLookupBool(&status, recipe, "CONSTANT_PHOTOMETRIC_WEIGHTS");
+    psAssert (status, "You must provide a value for the BOOL recipe CONSTANT_PHOTOMETRIC_WEIGHTS");
+
+    // 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
-    for (int i = 0; i < sources->n; i++) {
-        pmSource *source = sources->data[i];
+    for (int i = 0; i < objects->n; i++) {
+        pmPhotObj *object = objects->data[i];
+        if (!object) continue;
+        if (!object->sources) continue;
 
-        // turn this bit off and turn it on again if we pass this test
-        source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
+        // XXX check an element of the group to see if we should use it
+        // if (!object->flags & PM_PHOT_OBJ_BAD) continue;
 
-        // skip non-astronomical objects (very likely defects)
-        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+        for (int j = 0; j < object->sources->n; j++) {
+          pmSource *source = object->sources->data[j];
+          if (!source) continue;
 
-        // do not include CRs in the full ensemble fit
-        if (source->mode & PM_SOURCE_MODE_CR_LIMIT) continue;
+          // turn this bit off and turn it on again if we keep this source
+          source->mode &= ~PM_SOURCE_MODE_LINEAR_FIT;
 
-        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) {
+            if (!pmSourceCacheModel (source, maskVal)) continue;
+          }
 
-        // 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);
+    psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "built fitSources: %f sec (%ld objects)\n", psTimerMark ("psphot.linear"), objects->n);
 
     if (fitSources->n == 0) {
@@ -83 +76 @@
     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)
@@ -95 +86 @@
 
         // diagonal elements of the sparse matrix (auto-cross-product)
-        f = pmSourceModelDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        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);
+            float var = pmSourceModelDotModel (SRCi, SRCi, false, COVAR_FACTOR);
             errors->data.F32[i] = 1.0 / sqrt(var);
         } else {
@@ -106 +97 @@
         }
 
-
         // find the image x model value
-        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS);
+        f = pmSourceDataDotModel (SRCi, SRCi, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
         psSparseVectorElement (sparse, i, f);
-
-        f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
-        psSparseBorderElementB (border, i, 0, f);
 
         // loop over all other stars following this one
@@ -118 +105 @@
             pmSource *SRCj = fitSources->data[j];
 
+            // 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);
+            f = pmSourceModelDotModel (SRCi, SRCj, CONSTANT_PHOTOMETRIC_WEIGHTS, COVAR_FACTOR);
             psSparseMatrixElement (sparse, j, i, f);
         }
@@ -133 +123 @@
     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;
     constraint.paramMin   = 0.0;
@@ -142 +129 @@
 
     // 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);
-
-    // 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
@@ -164 +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
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
-        source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "sub models: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
 
     // 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, COVAR_FACTOR);
     }
     psLogMsg ("psphot.ensemble", PS_LOG_MINUTIA, "get chisqs: %f sec (%d elements)\n", psTimerMark ("psphot.linear"), sparse->Nelem);
@@ -186 +166 @@
     psFree (norm);
     psFree (errors);
-    psFree (border);
 
     psLogMsg ("psphot.ensemble", PS_LOG_INFO, "measure ensemble of PSFs: %f sec\n", psTimerMark ("psphot.linear"));
-
-    psphotVisualShowResidualImage (readout);
-    psphotVisualShowFlags (sources);
 
     return true;
 }
 
-// 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, psImageMaskType markVal) {
+// sort by X (ascending)
+int pmPhotObjSortByX (const void **a, const void **b)
+{
+    pmPhotObj *objA = *(pmPhotObj **)a;
+    pmPhotObj *objB = *(pmPhotObj **)b;
 
-    // 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);
+    psF32 fA = objA->x;
+    psF32 fB = objB->x;
 
-    // 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->radiusFit, "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;
-        }
-    }
-
-    // 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);
-
-    return true;
+    psF32 diff = fA - fB;
+    if (diff > FLT_EPSILON) return (+1);
+    if (diff < FLT_EPSILON) return (-1);
+    return (0);
 }

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