Changeset 29060 for branches/sc_branches/trunkTest/psphot/src/psphotBlendFit.c

Timestamp:

Aug 26, 2010, 9:18:39 AM (16 years ago)

Author:

Serge CHASTEL

Message:

Merging trunk in branch

Location:

branches/sc_branches/trunkTest

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src/psphotBlendFit.c (modified) (12 diffs)

branches/sc_branches/trunkTest/psphot/src/psphotBlendFit.c

@@ -65 +65 @@
     assert (status && fitIter > 0);
 
-    float fitTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_TOL"); // Fit tolerance
-    assert (status && isfinite(fitTol) && fitTol > 0);
+    float fitMinTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MIN_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+        fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance
+        if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+            psAbort("PSF_FIT_MIN_TOL (and PSF_FIT_TOL) not defined or positive");
+        }
+    }
+
+    float fitMaxTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) {
+        fitMaxTol = 1.0;
+    }
 
     bool poisson = psMetadataLookupBool(&status, recipe, "POISSON.ERRORS.PHOT.LMM"); // Poisson errors?
     assert (status);
 
+    float maxChisqDOF = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_CHISQ"); // Fit tolerance
+
     float skySig = psMetadataLookupF32(&status, recipe, "SKY_SIG");
     assert (status && isfinite(skySig) && skySig > 0);
 
     // Define source fitting parameters for extended source fits
-    pmSourceFitModelInit(fitIter, fitTol, PS_SQR(skySig), poisson);
+    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
+    fitOptions->nIter         = fitIter;
+    fitOptions->minTol        = fitMinTol;
+    fitOptions->maxTol        = fitMaxTol;
+    fitOptions->maxChisqDOF   = maxChisqDOF;
+    fitOptions->poissonErrors = poisson;
+    fitOptions->weight        = PS_SQR(skySig);
+    fitOptions->mode          = PM_SOURCE_FIT_PSF;
 
     psphotInitLimitsPSF (recipe, readout);
-    psphotInitLimitsEXT (recipe);
-    psphotInitRadiusPSF (recipe, readout->analysis, psf->type);
+    psphotInitLimitsEXT (recipe, readout);
+    psphotInitRadiusPSF (recipe, readout);
 
     // starts the timer, sets up the array of fitSets
@@ -88 +107 @@
     if (!sources->n) {
         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping blend");
+        psFree (fitOptions);
         return true;
     }
@@ -112 +132 @@
             psArrayAdd(job->args, 1, psf);
             psArrayAdd(job->args, 1, newSources); // return for new sources
+            psArrayAdd(job->args, 1, fitOptions); // default fit options
             psFree (newSources);
 
@@ -121 +142 @@
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree (fitOptions);
                 return NULL;
             }
-
-# if (0)
-            {
-                int nfit = 0;
-                int npsf = 0;
-                int next = 0;
-                int nfail = 0;
-                psArray *newSources = psArrayAllocEmpty(16);
-
-                if (!psphotBlendFit_Unthreaded (&nfit, &npsf, &next, &nfail, readout, recipe, cells->data[j], psf, newSources)) {
-                    psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                    return NULL;
-                }
-                Nfit += nfit;
-                Npsf += npsf;
-                Next += next;
-                Nfail += nfail;
-
-                // add these back onto sources
-                for (int k = 0; k < newSources->n; k++) {
-                    psArrayAdd (sources, 16, newSources->data[k]);
-                }
-                psFree (newSources);
-            }
-# endif
         }
 
@@ -153 +150 @@
         if (!psThreadPoolWait (false)) {
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            psFree (fitOptions);
             return NULL;
         }
@@ -163 +161 @@
             } else {
                 psScalar *scalar = NULL;
-                scalar = job->args->data[5];
+                scalar = job->args->data[6];
                 Nfit += scalar->data.S32;
-                scalar = job->args->data[6];
+                scalar = job->args->data[7];
                 Npsf += scalar->data.S32;
-                scalar = job->args->data[7];
+                scalar = job->args->data[8];
                 Next += scalar->data.S32;
-                scalar = job->args->data[8];
+                scalar = job->args->data[9];
                 Nfail += scalar->data.S32;
 
@@ -186 +184 @@
       psphotSaveImage (NULL, readout->image,  "image.v2.fits");
     }
+    psFree (fitOptions);
 
     psLogMsg ("psphot.psphotBlendFit", PS_LOG_INFO, "fit models: %f sec for %d objects (%d psf, %d ext, %d failed, %ld skipped)\n", psTimerMark ("psphot.fit.nonlinear"), Nfit, Npsf, Next, Nfail, sources->n - Nfit);
@@ -204 +203 @@
     psScalar *scalar = NULL;
 
-    pmReadout *readout  = job->args->data[0];
-    psMetadata *recipe  = job->args->data[1];
-    psArray *sources    = job->args->data[2];
-    pmPSF *psf          = job->args->data[3];
-    psArray *newSources = job->args->data[4];
+    pmReadout *readout             = job->args->data[0];
+    psMetadata *recipe             = job->args->data[1];
+    psArray *sources               = job->args->data[2];
+    pmPSF *psf                     = job->args->data[3];
+    psArray *newSources            = job->args->data[4];
+    pmSourceFitOptions *fitOptions = job->args->data[5];
 
     // bit-masks to test for good/bad pixels
@@ -269 +269 @@
         Nfit ++;
 
+        if (0) {
+            psF32 Mxx = source->moments->Mxx;
+            psF32 Myy = source->moments->Myy;
+            fprintf (stderr, "1: Mxx: %f, Myy: %f\n", Mxx, Myy);
+        }
+
         // try fitting PSFs or extended sources depending on source->mode
         // these functions subtract the resulting fitted source
         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-            if (psphotFitBlob (readout, source, newSources, psf, maskVal, markVal)) {
+            if (psphotFitBlob (readout, source, newSources, psf, fitOptions, maskVal, markVal)) {
                 source->type = PM_SOURCE_TYPE_EXTENDED;
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is ext", source->peak->xf, source->peak->yf);
@@ -280 +286 @@
             }
         } else {
-            if (psphotFitBlend (readout, source, psf, maskVal, markVal)) {
+            if (psphotFitBlend (readout, source, psf, fitOptions, maskVal, markVal)) {
                 source->type = PM_SOURCE_TYPE_STAR;
                 psTrace ("psphot", 5, "source at %7.1f, %7.1f is psf", source->peak->xf, source->peak->yf);
@@ -289 +295 @@
         }
 
+        if (0) {
+            psF32 Mxx = source->moments->Mxx;
+            psF32 Myy = source->moments->Myy;
+            fprintf (stderr, "2: Mxx: %f, Myy: %f\n", Mxx, Myy);
+        }
+
         psTrace ("psphot", 5, "source at %7.1f, %7.1f failed", source->peak->xf, source->peak->yf);
         Nfail ++;
@@ -298 +310 @@
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    scalar = job->args->data[5];
+    scalar = job->args->data[6];
     scalar->data.S32 = Nfit;
 
-    scalar = job->args->data[6];
+    scalar = job->args->data[7];
     scalar->data.S32 = Npsf;
 
-    scalar = job->args->data[7];
+    scalar = job->args->data[8];
     scalar->data.S32 = Next;
 
-    scalar = job->args->data[8];
+    scalar = job->args->data[9];
     scalar->data.S32 = Nfail;