Changeset 27838 for branches/tap_branches/psphot/src/psphotFitSourcesLinear.c

Timestamp:

May 3, 2010, 8:41:49 AM (16 years ago)

Author:

eugene

Message:

updates from trunk

Location:

branches/tap_branches

Files:

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

branches/tap_branches/psphot/src

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack

branches/tap_branches/psphot/src/psphotFitSourcesLinear.c

@@ -12 +12 @@
 static bool SetBorderMatrixElements (psSparseBorder *border, pmReadout *readout, psArray *sources, bool constant_weights, int SKY_FIT_ORDER, psImageMaskType markVal);
 
-bool psphotFitSourcesLinear (pmReadout *readout, psArray *sources, const psMetadata *recipe, const pmPSF *psf, bool final) {
+// for now, let's store the detections on the readout->analysis for each readout
+bool psphotFitSourcesLinear (pmConfig *config, const pmFPAview *view, bool final)
+{
+    bool status = true;
+
+    // select the appropriate recipe information
+    psMetadata *recipe  = psMetadataLookupPtr (&status, config->recipes, PSPHOT_RECIPE);
+    assert (recipe);
+
+    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++) {
+
+        // find the currently selected readout
+        pmFPAfile *file = pmFPAfileSelectSingle(config->files, "PSPHOT.INPUT", 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 PSPHOT.INPUT entry %d", i);
+            return false;
+        }
+    }
+    return true;
+}
+
+bool psphotFitSourcesLinearReadout (psMetadata *recipe, pmReadout *readout, psArray *sources, pmPSF *psf, bool final) {
 
     bool status;
@@ -18 +57 @@
     float y;
     float f;
-    // float r;
+
+    if (!sources->n) {
+        psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping linear fit");
+        return true;
+    }
 
     psTimerStart ("psphot.linear");
@@ -59 +102 @@
     }
 
+    // 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++) {
@@ -120 +170 @@
 
         // 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 {
@@ -133 +183 @@
 
         // 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);
 
@@ -139 +189 @@
         switch (SKY_FIT_ORDER) {
           case 1:
-            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelWeight (SRCi, 1, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseBorderElementB (border, i, 1, f);
-            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelWeight (SRCi, 2, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseBorderElementB (border, i, 2, f);
 
           case 0:
-            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS);
+            f = pmSourceModelWeight (SRCi, 0, CONSTANT_PHOTOMETRIC_WEIGHTS, covarFactor);
             psSparseBorderElementB (border, i, 0, f);
             break;
@@ -165 +215 @@
 
             // 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);
         }
@@ -219 +269 @@
 
     // 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);
@@ -238 +289 @@
 
     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;
@@ -244 +300 @@
 
 // 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

@@ -19 +19 @@
 psphotMomentsStudy
 psphotPetrosianStudy
+psphotForced
+psphotMakePSF
+psphotStack