Changeset 35769 for trunk

Timestamp:

Jul 3, 2013, 2:43:13 PM (13 years ago)

Author:

eugene

Message:

thread psphotModelBackground, psphotAddOrSubNoise, psphotGuessModel; add recipe options to control selection of extended models fitted; iterate on the trail window a bit

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src (modified) (1 prop)
• src/models/pmModel_STRAIL.c (modified) (1 diff)
• src/models/pmModel_TEST1.c (modified) (1 diff)
• src/psphot.h (modified) (4 diffs)
• src/psphotAddNoise.c (modified) (6 diffs)
• src/psphotBlendFit.c (modified) (6 diffs)
• src/psphotChoosePSF.c (modified) (2 diffs)
• src/psphotExtendedSourceFits.c (modified) (18 diffs)
• src/psphotGuessModels.c (modified) (5 diffs)
• src/psphotKronIterate.c (modified) (1 diff)
• src/psphotModelBackground.c (modified) (9 diffs)
• src/psphotModelGroupInit.c (modified) (1 diff)
• src/psphotModelTestReadout.c (modified) (5 diffs)
• src/psphotOutput.c (modified) (1 diff)
• src/psphotRadialPlot.c (modified) (1 diff)
• src/psphotReplaceUnfit.c (modified) (2 diffs)
• src/psphotSetThreads.c (modified) (2 diffs)
• src/psphotSourceFits.c (modified) (8 diffs)
• src/psphotStackImageLoop.c (modified) (1 prop)
• src/psphotVisual.c (modified) (1 diff)
• test/tap_psphot_galaxies.pro (copied) (copied from branches/eam_branches/ipp-20130509/psphot/test/tap_psphot_galaxies.pro )
• test/tap_psphot_galaxygrid.pro (copied) (copied from branches/eam_branches/ipp-20130509/psphot/test/tap_psphot_galaxygrid.pro )

trunk/psphot/src/models/pmModel_STRAIL.c

@@ -563 +563 @@
     // convert to shape terms (SXX,SYY,SXY)
     // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;

trunk/psphot/src/models/pmModel_TEST1.c

@@ -264 +264 @@
     // convert to shape terms (SXX,SYY,SXY)
     // XXX user-defined value for limit?
-    if (!pmPSF_FitToModel (PAR, 0.1)) {
+    bool useReff = pmModelUseReff (model->type);
+    if (!pmPSF_FitToModel (PAR, 0.1, useReff)) {
         psError(PM_ERR_PSF, false, "Failed to fit object at (r,c) = (%.1f,%.1f)", Xo, Yo);
         return false;

trunk/psphot/src/psphot.h

@@ -67 +67 @@
 bool            psphotModelBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
 bool            psphotModelBackgroundReadoutFileIndex (pmConfig *config, const pmFPAview *view, const char *filerule, int index);
+bool            psphotModelBackground_Threaded (psThreadJob *job);
 
 bool            psphotMaskBackground (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -119 +120 @@
 bool            psphotAddOrSubNoise (pmConfig *config, const pmFPAview *view, const char *filerule, bool add);
 bool            psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add);
+bool            psphotAddOrSubNoise_Threaded (psThreadJob *job);
 
 bool            psphotExtendedSourceAnalysis (pmConfig *config, const pmFPAview *view, const char *filerule);
@@ -232 +234 @@
 bool            psphotFitBlob (pmReadout *readout, pmSource *source, psArray *sources, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
 bool            psphotFitPSF (pmReadout *readout, pmSource *source, pmPSF *psf, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
-pmModel        *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal);
+pmModel        *psphotFitEXT (pmModel *guessModel, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal);
 psArray        *psphotFitDBL (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, psImageMaskType maskVal, psImageMaskType markVal);
 
@@ -243 +245 @@
 bool            psphotFitInit (int nThreads);
 bool            psphotFitSummary (void);
+bool            psphotFitSummaryExtended (void);
+bool            psphotFitInitExtended (void);
 
 bool            psphotLoadPSF (pmConfig *config, const pmFPAview *view, const char *filerule);

trunk/psphot/src/psphotAddNoise.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+bool psphotAddOrSubNoise_Threaded (psThreadJob *job);
 bool psphotMaskSource(pmSource *source, bool add, psImageMaskType maskVal);
 
@@ -32 +33 @@
 }
 
-static int Nmasked = 0;
-
 // the return state indicates if any sources were actually replaced
 bool psphotAddOrSubNoiseReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, bool add) {
 
     bool status = false;
+
+    psTimerStart ("psphot.noise");
 
     // find the currently selected readout
@@ -55 +56 @@
     if (!sources) return false;
 
-    psTimerStart ("psphot.noise");
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
@@ -84 +89 @@
 
     psphotVisualShowImage (readout);
+
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_ADD_NOISE");
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
+            PS_ARRAY_ADD_SCALAR(job->args, FACTOR,   PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, SIZE,     PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, add,      PS_TYPE_U8);
+
+# if (1)
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise.");
+                return false;
+            }
+# else
+            if (!psphotAddOrSubNoise_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise.");
+                return false;
+            }
+            psFree(job);
+# endif
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false, true)) {
+            psFree(cellGroups);
+            psError(PS_ERR_UNKNOWN, false, "Unable to add/sub noise.");
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            // we have no returned data from this operation
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            }
+            psFree(job);
+        }
+    }
+
+    if (add) {
+        psLogMsg ("psphot.noise", PS_LOG_WARN, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
+    } else {
+        psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
+    }
+
+    psFree (cellGroups);
+
+    psphotVisualShowImage (readout);
+
+    return true;
+}
+
+bool psphotAddOrSubNoise_Threaded (psThreadJob *job) {
+
+    psArray *sources = job->args->data[0];
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[1], PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[2], PS_TYPE_IMAGE_MASK_DATA);
+    psF32 FACTOR = PS_SCALAR_VALUE(job->args->data[3], F32);
+    psF32 SIZE = PS_SCALAR_VALUE(job->args->data[4], F32);
+    psBool add  = PS_SCALAR_VALUE(job->args->data[5], U8);
 
     // loop over all source
@@ -104 +184 @@
         psphotMaskSource (source, add, markVal);
     }
-    if (add) {
-        psLogMsg ("psphot.noise", PS_LOG_WARN, "add noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
-    } else {
-        psLogMsg ("psphot.noise", PS_LOG_WARN, "sub noise for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.noise"));
-    }
-    fprintf (stderr, "masked %d objects\n", Nmasked);
-
-    psphotVisualShowImage (readout);
-
     return true;
 }
@@ -142 +213 @@
         }
     }
-    Nmasked ++;
-
-    return true;
-}
-
+    return true;
+}
+

trunk/psphot/src/psphotBlendFit.c

@@ -86 +86 @@
     if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) {
         fitMaxTol = 1.0;
+    }
+
+    bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        chisqConvergence = true;
+    }
+
+    int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        gainFactorMode = 0;
     }
 
@@ -119 +131 @@
     fitOptions->mode          = PM_SOURCE_FIT_PSF;
     fitOptions->covarFactor   = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+
+    fitOptions->gainFactorMode = gainFactorMode;
+    fitOptions->chisqConvergence = chisqConvergence;
 
     psphotInitLimitsPSF (recipe, readout);
@@ -211 +226 @@
 
     psLogMsg ("psphot.psphotBlendFit", PS_LOG_WARN, "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);
+    psphotFitSummary ();
 
     psphotVisualShowResidualImage (readout, false);
@@ -271 +287 @@
 
         // skip non-astronomical objects (very likely defects)
-        if (source->mode &  PM_SOURCE_MODE_BLEND) continue;
-        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) continue;
-        if (source->type == PM_SOURCE_TYPE_DEFECT) continue;
-        if (source->type == PM_SOURCE_TYPE_SATURATED) continue;
+        if (source->mode &  PM_SOURCE_MODE_BLEND) goto skip_blend;
+        if (source->mode &  PM_SOURCE_MODE_CR_LIMIT) goto skip_cr;
+        if (source->type == PM_SOURCE_TYPE_DEFECT) goto skip_defect;
+        if (source->type == PM_SOURCE_TYPE_SATURATED) goto skip_sat;
 
         // skip saturated stars modeled with a radial profile 
-        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) continue;
+        if (source->mode2 & PM_SOURCE_MODE2_SATSTAR_PROFILE) goto skip_sat;
 
         // skip DBL second sources (ie, added by psphotFitBlob
-        if (source->mode &  PM_SOURCE_MODE_PAIR) continue;
+        if (source->mode &  PM_SOURCE_MODE_PAIR) goto skip_blend;
 
         // do not include MOMENTS_FAILURES in the fit
-        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
+        if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) goto skip_generic;
 
         // limit selection to some SN limit
         if (source->mode & PM_SOURCE_MODE_EXT_LIMIT) {
-            if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) continue;
+            if (source->moments->KronFlux < FIT_SN_LIM * source->moments->KronFluxErr) goto skip_generic;
         } else {
-            if (sqrt(source->peak->detValue) < FIT_SN_LIM) continue;
+            if (sqrt(source->peak->detValue) < FIT_SN_LIM) goto skip_generic;
         }
         // exclude sources outside optional analysis region
-        if (source->peak->xf < AnalysisRegion.x0) continue;
-        if (source->peak->yf < AnalysisRegion.y0) continue;
-        if (source->peak->xf > AnalysisRegion.x1) continue;
-        if (source->peak->yf > AnalysisRegion.y1) continue;
+        if (source->peak->xf < AnalysisRegion.x0) goto skip_generic;
+        if (source->peak->yf < AnalysisRegion.y0) goto skip_generic;
+        if (source->peak->xf > AnalysisRegion.x1) goto skip_generic;
+        if (source->peak->yf > AnalysisRegion.y1) goto skip_generic;
 
         // if model is NULL, we don't have a starting guess
-        if (source->modelPSF == NULL) continue;
+        if (source->modelPSF == NULL) goto skip_generic;
 
         // skip sources which are insignificant flux?
@@ -307 +323 @@
                      source->modelPSF->params->data.F32[2],
                      source->modelPSF->params->data.F32[3]);
-            continue;
+            goto skip_generic;
         }
 
@@ -357 +373 @@
         pmSourceCacheModel (source, maskVal);
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+        continue;
+
+    skip_blend:
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf);
+        continue;
+    skip_cr:
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf);
+        continue;
+    skip_defect:
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf);
+        continue;
+    skip_sat:
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf);
+        continue;
+    skip_generic:
+        psTrace ("psphot", 5, "source at %7.1f, %7.1f skipped", source->peak->xf, source->peak->yf);
+        continue;
     }
 

trunk/psphot/src/psphotChoosePSF.c

@@ -158 +158 @@
     }
     float maxChisqDOF = psMetadataLookupF32 (&status, recipe, "PSF_FIT_MAX_CHISQ"); // Fit tolerance
+
+    bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        chisqConvergence = true;
+    }
+
+    int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        gainFactorMode = 0;
+    }
 
     // options which modify the behavior of the model fitting
@@ -169 +181 @@
     options->fitOptions->mode          = PM_SOURCE_FIT_PSF;
     options->fitOptions->covarFactor   = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
-
     
+    options->fitOptions->gainFactorMode   = gainFactorMode;
+    options->fitOptions->chisqConvergence = chisqConvergence;
+
     psArray *stars = psArrayAllocEmpty (sources->n);
 

trunk/psphot/src/psphotExtendedSourceFits.c

@@ -1 +1 @@
 # include "psphotInternal.h"
+bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius);
 
 // for now, let's store the detections on the readout->analysis for each readout
@@ -52 +53 @@
     psTimerStart ("psphot.extended");
 
+    psphotFitInitExtended();
+
     // find the currently selected readout
     pmFPAfile *file = pmFPAfileSelectSingle(config->files, filerule, index); // File of interest
@@ -70 +73 @@
     if (!sources->n) {
         psLogMsg ("psphot", PS_LOG_INFO, "no sources, skipping source size");
+        psphotSersicModelClassCleanup();
         return true;
     }
@@ -80 +84 @@
     // do not thread if we are trying to study the fitting process
     if (psTraceGetLevel ("psphot.psphotFitEXT") >= 6) {
-      nThreads = 0;
+        nThreads = 0;
     }
 
@@ -107 +111 @@
     }
 
+    bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        chisqConvergence = true;
+    }
+
+    int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance
+    if (!status) {
+        // default to the old method (chisqConvergence)
+        gainFactorMode = 0;
+    }
+
+    // Define source fitting parameters for extended source fits
+    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
+    fitOptions->mode           = PM_SOURCE_FIT_EXT;
+    fitOptions->saveCovariance = true;  // XXX make this a user option?
+    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
+    fitOptions->nIter          = fitIter;
+    fitOptions->minTol         = fitMinTol;
+    fitOptions->maxTol         = fitMaxTol;
+
+    fitOptions->gainFactorMode   = gainFactorMode;
+    fitOptions->chisqConvergence = chisqConvergence;
+
     // maskVal is used to test for rejected pixels, and must include markVal
     maskVal |= markVal;
 
     // select the collection of desired models
-    psMetadata *models = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS");
+    psMetadata *allModels = psMetadataLookupMetadata (&status, recipe, "EXTENDED_SOURCE_MODELS");
     if (!status) {
         psWarning ("extended source model fits requested but model model is missing (EXTENDED_SOURCE_MODELS)\n");
+        psFree (fitOptions);
+        psphotSersicModelClassCleanup();
         return true;
     }
+    if (allModels->list->n == 0) {
+        psWarning ("extended source model fits requested but no models are specified\n");
+        psFree (fitOptions);
+        psphotSersicModelClassCleanup();
+        return true;
+    }
+
+    psMetadata *models = NULL;
+    char *modelSelection = psMetadataLookupStr (&status, recipe, "EXTENDED_SOURCE_MODELS_SELECTION");
+    if (!modelSelection || !status) {
+        models = psMemIncrRefCounter(allModels);
+        goto got_models;
+    } 
+
+    if (!strcasecmp(modelSelection, "ALL")) {
+        models = psMemIncrRefCounter(allModels);
+        goto got_models;
+    }
+
+    if (!strcasecmp(modelSelection, "NONE")) {
+        psWarning ("extended source model fits requested but no models are selected (EXTENDED_SOURCE_MODELS_SELECTION = NONE)\n");
+        psFree (fitOptions);
+        psphotSersicModelClassCleanup();
+        return true;
+    }
+
+    psArray *selection = psStringSplitArray (modelSelection, ",", false);
+    if (selection->n == 0) {
+        psWarning ("extended source model fits requested but model selection string is empty (EXTENDED_SOURCE_MODELS_SELECTION)\n");
+        psFree (fitOptions);
+        psFree (selection);
+        psphotSersicModelClassCleanup();
+        return true;
+    }
+    models = psMetadataAlloc();
+    for (int i = 0; i < selection->n; i++) {
+        psMetadata *model = psMetadataLookupMetadata (&status, allModels, selection->data[i]);
+        if (!model) {
+            psWarning ("extended source model selection string includes an invalid name %s\n", (char *) selection->data[i]);
+            continue;
+        }
+        psMetadataAddMetadata (models, PS_LIST_TAIL, selection->data[i], PS_META_REPLACE, "extended model", model);
+    }
+    psFree (selection);
+
     if (models->list->n == 0) {
-        psWarning ("extended source model fits requested but no models are specified\n");
-        return true;
-    }
+        psWarning ("extended source model fits requested but no valid models in selection string (EXTENDED_SOURCE_MODELS_SELECTION = %s)\n", modelSelection);
+        psFree (fitOptions);
+        psFree (models);
+        psphotSersicModelClassCleanup();
+        return true;
+    }
+
+got_models:
 
     psphotInitRadiusEXT (recipe, readout);
@@ -128 +208 @@
     while ((item = psMetadataGetAndIncrement (iter)) != NULL) {
 
-      if (item->type != PS_DATA_METADATA) {
-        // XXX we could cull the bad entries or build a validated model folder
-        psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);
-      }
-
-      psMetadata *model = (psMetadata *) item->data.md;
-
-      // check on the model type
-      char *modelName = psMetadataLookupStr (&status, model, "MODEL");
-      int modelType = pmModelClassGetType (modelName);
-      if (modelType < 0) {
-        psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName);
-      }
-      psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType);
-
-      // check on the SNLIM, set a float value
-      char *SNword = psMetadataLookupStr (&status, model, "SNLIM");
-      if (!status) {
-        psAbort("SNLIM not defined for extended source model %s\n", item->name);
-      }
-      float SNlim = atof (SNword);
-      psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim);
-
-      // check on the PSF-Convolution status
-      char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED");
-      if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) {
-        psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name);
-      }
-      bool convolved = !strcasecmp (convolvedWord, "true");
-      psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved);
+        if (item->type != PS_DATA_METADATA) {
+            // XXX we could cull the bad entries or build a validated model folder
+            psAbort ("Invalid type for EXTENDED_SOURCE_MODEL entry %s, not a metadata folder", item->name);
+        }
+
+        psMetadata *model = (psMetadata *) item->data.md;
+
+        // check on the model type
+        char *modelName = psMetadataLookupStr (&status, model, "MODEL");
+        int modelType = pmModelClassGetType (modelName);
+        if (modelType < 0) {
+            psAbort ("Unknown model class for EXTENDED_SOURCE_MODEL entry %s: %s", item->name, modelName);
+        }
+        psMetadataAddS32 (model, PS_LIST_TAIL, "MODEL_TYPE", PS_META_REPLACE, "", modelType);
+
+        // check on the SNLIM, set a float value
+        char *SNword = psMetadataLookupStr (&status, model, "SNLIM");
+        if (!status) {
+            psAbort("SNLIM not defined for extended source model %s\n", item->name);
+        }
+        float SNlim = atof (SNword);
+        psMetadataAddF32 (model, PS_LIST_TAIL, "SNLIM_VALUE", PS_META_REPLACE, "", SNlim);
+
+        // check on the PSF-Convolution status
+        char *convolvedWord = psMetadataLookupStr (&status, model, "PSF_CONVOLVED");
+        if (!status || (strcasecmp (convolvedWord, "true") && strcasecmp (convolvedWord, "false"))) {
+            psAbort ("PSF_CONVOLVED entry invalid or missing for EXTENDED_SOURCE_MODEL entry %s", item->name);
+        }
+        bool convolved = !strcasecmp (convolvedWord, "true");
+        psMetadataAddBool (model, PS_LIST_TAIL, "PSF_CONVOLVED_VALUE", PS_META_REPLACE, "", convolved);
+
+        if (convolved) {
+            psLogMsg ("psphot", PS_LOG_INFO, "using convolved model class %s (%s)", modelName, item->name);
+        } else {
+            psLogMsg ("psphot", PS_LOG_INFO, "using simple    model class %s (%s)", modelName, item->name);
+        }
     }
     psFree (iter);
@@ -195 +281 @@
             psMetadataIterator *iter = psMetadataIteratorAlloc (models, PS_LIST_HEAD, NULL);
             psArrayAdd(job->args, 1, iter);
-            psArrayAdd(job->args, 1, AnalysisRegion); // XXX make a pointer
+            psArrayAdd(job->args, 1, AnalysisRegion);
+            psArrayAdd(job->args, 1, fitOptions);
 
             PS_ARRAY_ADD_SCALAR(job->args, psfSize, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, fitIter, PS_TYPE_S32);
-            PS_ARRAY_ADD_SCALAR(job->args, fitMinTol, PS_TYPE_F32);
-            PS_ARRAY_ADD_SCALAR(job->args, fitMaxTol, PS_TYPE_F32);
-
             PS_ARRAY_ADD_SCALAR(job->args, maskVal, PS_TYPE_IMAGE_MASK);
             PS_ARRAY_ADD_SCALAR(job->args, markVal, PS_TYPE_IMAGE_MASK);
@@ -218 +301 @@
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
                 psFree(AnalysisRegion);
+                psFree (fitOptions);
+                psFree (models);
+                psphotSersicModelClassCleanup();
                 return false;
             } 
@@ -224 +310 @@
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
                 psFree(AnalysisRegion);
+                psFree (fitOptions);
+                psFree (models);
+                psphotSersicModelClassCleanup();
                 return false;
             }
             psScalar *scalar = NULL;
-            scalar = job->args->data[8];
+            scalar = job->args->data[9];
             Next += scalar->data.S32;
-            scalar = job->args->data[9];
+            scalar = job->args->data[10];
             Nconvolve += scalar->data.S32;
-            scalar = job->args->data[10];
+            scalar = job->args->data[11];
             NconvolvePass += scalar->data.S32;
-            scalar = job->args->data[11];
+            scalar = job->args->data[12];
             Nplain += scalar->data.S32;
-            scalar = job->args->data[12];
+            scalar = job->args->data[13];
             NplainPass += scalar->data.S32;
-            scalar = job->args->data[13];
+            scalar = job->args->data[14];
             Nfaint += scalar->data.S32;
-            scalar = job->args->data[14];
+            scalar = job->args->data[15];
             Nfail += scalar->data.S32;
             psFree(job->args->data[3]); // iterator allocated above
@@ -250 +339 @@
             psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
             psFree(AnalysisRegion);
+            psFree (fitOptions);
+            psFree (models);
+            psphotSersicModelClassCleanup();
             return false;
         }
@@ -260 +352 @@
             } else {
                 psScalar *scalar = NULL;
-                scalar = job->args->data[8];
+                scalar = job->args->data[9];
                 Next += scalar->data.S32;
-                scalar = job->args->data[9];
+                scalar = job->args->data[10];
                 Nconvolve += scalar->data.S32;
-                scalar = job->args->data[10];
+                scalar = job->args->data[11];
                 NconvolvePass += scalar->data.S32;
-                scalar = job->args->data[11];
+                scalar = job->args->data[12];
                 Nplain += scalar->data.S32;
-                scalar = job->args->data[12];
+                scalar = job->args->data[13];
                 NplainPass += scalar->data.S32;
-                scalar = job->args->data[13];
+                scalar = job->args->data[14];
                 Nfaint += scalar->data.S32;
-                scalar = job->args->data[14];
+                scalar = job->args->data[15];
                 Nfail += scalar->data.S32;
                 psFree(job->args->data[3]); // metadata iterator allocated above
@@ -281 +373 @@
     psFree (cellGroups);
     psFree(AnalysisRegion);
+    psFree (fitOptions);
+    psFree (models);
 
     psphotSersicModelClassCleanup();
@@ -290 +384 @@
     psLogMsg ("psphot", PS_LOG_INFO, "  %d plain models (%d passed)\n", Nplain, NplainPass);
     psLogMsg ("psphot", PS_LOG_INFO, "  %d too faint to fit, %d failed\n", Nfaint, Nfail);
+
+    psphotFitSummaryExtended();
+
     return true;
 }
@@ -308 +405 @@
 
     // arguments: readout, sources, models, region, psfSize, maskVal, markVal
-    pmReadout *readout      = job->args->data[0];
-    psArray *sources        = job->args->data[1];
-    psMetadata *models      = job->args->data[2];
+    pmReadout *readout       = job->args->data[0];
+    psArray *sources         = job->args->data[1];
+    psMetadata *models       = job->args->data[2];
     psMetadataIterator *iter = job->args->data[3];
-    psRegion *region        = job->args->data[4];
-    int psfSize             = PS_SCALAR_VALUE(job->args->data[5],S32);
-    int fitIter             = PS_SCALAR_VALUE(job->args->data[6],S32);
-    float fitMinTol         = PS_SCALAR_VALUE(job->args->data[7],F32);
-    float fitMaxTol         = PS_SCALAR_VALUE(job->args->data[8],F32);
-    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[9],PS_TYPE_IMAGE_MASK_DATA);
-    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[10],PS_TYPE_IMAGE_MASK_DATA);
-
-    // Define source fitting parameters for extended source fits
-    pmSourceFitOptions *fitOptions = pmSourceFitOptionsAlloc();
-    fitOptions->mode           = PM_SOURCE_FIT_EXT;
-    fitOptions->saveCovariance = true;  // XXX make this a user option?
-    fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
-    fitOptions->nIter          = fitIter;
-    fitOptions->minTol         = fitMinTol;
-    fitOptions->maxTol         = fitMaxTol;
+    psRegion *region         = job->args->data[4];
+    pmSourceFitOptions *fitOptions = job->args->data[5];
+
+    int psfSize             = PS_SCALAR_VALUE(job->args->data[6],S32);
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[8],PS_TYPE_IMAGE_MASK_DATA);
+
+    // psTraceSetLevel ("psLib.math.psMinimizeLMChi2_Alt", 5);
 
     // choose the sources of interest
@@ -348 +437 @@
         if (source->peak->y > region->y1) continue;
 
+
+        // XXX for a test, just do the obvious trail
+        // XXX if (source->peak->xf < 1100) continue;
+        // XXX if (source->peak->xf > 1400) continue;
+        // XXX if (source->peak->yf >  245) continue;
+
         // replace object in image
         if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
@@ -384 +479 @@
           assert (status);
 
-          // limit selection to some SN limit
-          // assert (source->peak); // how can a source not have a peak?
           // limit selection to some SN limit
           bool skipSource = false;
@@ -424 +517 @@
               }
           } else {
-              psFree (source->modelFlux);
-              source->modelFlux = NULL;
-              modelFit = psphotFitEXT (readout, source, fitOptions, modelType, maskVal, markVal);
-              if (!modelFit) {
-                  psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
-                  Nfail ++;
-                  continue;
-              }
-              psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", 
-                       source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
-              Nplain ++;
-              if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
-                  NplainPass ++;
-                  source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
-              }
+              bool doneFits = false;
+              while (!doneFits) {
+                  psFree (source->modelFlux);
+                  source->modelFlux = NULL;
+                  modelFit = psphotFitEXT (modelFit, readout, source, fitOptions, modelType, maskVal, markVal);
+                  if (!modelFit) {
+                      psTrace ("psphot", 5, "failed to fit plain model for object at %f, %f", source->moments->Mx, source->moments->My);
+                      Nfail ++;
+                      doneFits = true;
+                      continue;
+                  }
+                  psTrace ("psphot", 4, "fit plain model for %f, %f : %s chisq = %f (npix: %d, niter: %d)\n", source->moments->Mx, source->moments->My, pmModelClassGetName (modelFit->type), modelFit->chisq, modelFit->nPix, modelFit->nIter);
+                  Nplain ++;
+                  if (!(modelFit->flags & (PM_MODEL_STATUS_BADARGS | PM_MODEL_STATUS_NONCONVERGE | PM_MODEL_STATUS_OFFIMAGE))) {
+                      NplainPass ++;
+                      source->mode |= PM_SOURCE_MODE_EXTENDED_FIT;
+                  }
+                  doneFits = true;
+
+                  if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
+                      // calculate the object end points:
+                      float *PAR = modelFit->params->data.F32;
+                      float Xmin = PAR[PM_PAR_XPOS] - 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]);
+                      float Ymin = PAR[PM_PAR_YPOS] - 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]);
+                      float Xmax = PAR[PM_PAR_XPOS] + 0.5*PAR[PM_PAR_LENGTH]*sin(PAR[PM_PAR_THETA]);
+                      float Ymax = PAR[PM_PAR_YPOS] + 0.5*PAR[PM_PAR_LENGTH]*cos(PAR[PM_PAR_THETA]);
+
+                      if (false && (source->peak->xf > 1100) && 
+                          (source->peak->xf < 1400) && 
+                          (source->peak->yf < 245)) {
+                          fprintf (stderr, "src vs fit : %d %d - %d %d | %f %f - %f %f\n", 
+                                   source->pixels->col0, source->pixels->row0, 
+                                   source->pixels->col0 + source->pixels->numCols, source->pixels->row0 + source->pixels->numRows, 
+                                   Xmin, Ymin, Xmax, Ymax);
+                      }
+                      if (PAR[PM_PAR_LENGTH] > 0.9*fitRadius) {
+                          doneFits = false;
+                          fprintf (stderr, "update window : %f %f : %f -> %f\n", source->peak->xf, source->peak->yf, fitRadius, 2*fitRadius);
+                          psphotSetWindowTrail (&fitRadius, &windowRadius, readout, source, markVal, fitRadius*2.0);
+                      }
+                  }
+              }
           }
-
-          // XXX deprecate?
-          // XXX pmSourceExtFitPars *extFitPars = pmSourceExtFitParsAlloc();
-          // XXX extFitPars->Mxx = source->moments->Mxx;
-          // XXX extFitPars->Mxy = source->moments->Mxy;
-          // XXX extFitPars->Myy = source->moments->Myy;
-          // XXX extFitPars->Mrf = source->moments->Mrf;
-          // XXX extFitPars->Mrh = source->moments->Mrh;
-          // XXX extFitPars->peakMag = (source->peak->rawFlux > 0) ? -2.5*log10(source->peak->rawFlux) : NAN;
-          
-          // save kron mag, but assign apMag & psfMag on output (not yet calculated)
-          // XXX extFitPars->krMag = -2.5*log10(source->moments->KronFlux);
-
-          // psArrayAdd (source->extFitPars, 4, extFitPars);
-          // psFree (extFitPars);
+          // XXX really need to do this in a cleaner way:
+          if (!modelFit) continue;
 
           // test for fit quality / result
@@ -519 +626 @@
         psTrace ("psphot", 5, "extended source model for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
     }
-    psFree (fitOptions);
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
+    scalar = job->args->data[9];
+    scalar->data.S32 = Next;
+
+    scalar = job->args->data[10];
+    scalar->data.S32 = Nconvolve;
+
     scalar = job->args->data[11];
-    scalar->data.S32 = Next;
+    scalar->data.S32 = NconvolvePass;
 
     scalar = job->args->data[12];
-    scalar->data.S32 = Nconvolve;
+    scalar->data.S32 = Nplain;
 
     scalar = job->args->data[13];
-    scalar->data.S32 = NconvolvePass;
+    scalar->data.S32 = NplainPass;
 
     scalar = job->args->data[14];
-    scalar->data.S32 = Nplain;
+    scalar->data.S32 = Nfaint;
 
     scalar = job->args->data[15];
-    scalar->data.S32 = NplainPass;
-
-    scalar = job->args->data[16];
-    scalar->data.S32 = Nfaint;
-
-    scalar = job->args->data[17];
     scalar->data.S32 = Nfail;
 
     return true;
 }
+
+# define PAD_WINDOW 3.0
+bool psphotSetWindowTrail (float *fitRadius, float *windowRadius, pmReadout *readout, pmSource *source, psImageMaskType markVal, float newRadius) {
+
+    psRegion newRegion;
+
+    psAssert (source, "source not defined??");
+    psAssert (source->moments, "moments not defined??");
+
+    *fitRadius = newRadius;
+    *windowRadius = *fitRadius + PAD_WINDOW;
+
+    // check to see if new region is completely contained within old region
+    newRegion = psRegionForSquare (source->peak->xf, source->peak->yf, *windowRadius);
+    newRegion = psRegionForImage (readout->image, newRegion);
+
+    // redefine the pixels to match
+    pmSourceRedefinePixelsByRegion (source, readout, newRegion);
+
+    // clear the circular mask
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); 
+
+    // set the mask to flag the excluded pixels
+    psImageKeepCircle (source->maskObj, source->peak->xf, source->peak->yf, *fitRadius, "OR", markVal);
+
+    return true;
+}

trunk/psphot/src/psphotGuessModels.c

@@ -29 +29 @@
 }
 
+int NpixTotal = 0;
+
 // construct an initial PSF model for each object (new sources only)
 bool psphotGuessModelsReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index) {
@@ -86 +88 @@
 
     psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    NpixTotal = 0;
 
     for (int i = 0; i < cellGroups->n; i++) {
@@ -103 +107 @@
             PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
 
+# if (1)
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
                 return false;
             }
-        }
+# else
+            if (!psphotGuessModel_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                return false;
+            }
+            psFree(job);
+# endif
+        }
+
 
         // wait for the threads to finish and manage results
@@ -143 +156 @@
     psFree (cellGroups);
 
-    psLogMsg ("psphot.models", PS_LOG_WARN, "built models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.models"));
+    psLogMsg ("psphot.models", PS_LOG_WARN, "built models for %ld objects: %f sec (%d pixels)\n", sources->n, psTimerMark ("psphot.models"), NpixTotal);
     return true;
 }
@@ -240 +253 @@
 
         pmSourceCacheModel (source, maskVal);  // ALLOC x14 (!)
+        NpixTotal += source->pixels->numCols * source->pixels->numRows;
     }
 

trunk/psphot/src/psphotKronIterate.c

@@ -166 +166 @@
     }
 
+    fprintf (stderr, "--- starting KRON ---\n");
+
     // We measure the Kron Radius on a smoothed copy of the readout image
     psImage *smoothedImage = NULL;
     if (KRON_SMOOTH) {
+        psTimerStart ("psphot.kron.smooth");
+
         // Build the smoothed source image
         // Replace the subtracted sources
         psphotReplaceAllSourcesReadout(config, view, filerule, index, recipe, false);
+        psLogMsg ("psphot.kron", PS_LOG_INFO, "replaced %ld raw sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth"));
+
+        // smoothedImage = psImageCopy(NULL, readout->image, PS_TYPE_F32);
+        // psImageSmooth(smoothedImage, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA);
+        
         // Copy the image and smooth
         psTimerStart ("psphot.kron.smooth");
-        smoothedImage = psImageCopy(NULL, readout->image, PS_TYPE_F32);
-        psImageSmooth(smoothedImage, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA);
+        bool oldThreads = psImageConvolveSetThreads(true); // Old value of threading in psImageConvolve
+        smoothedImage = psImageSmoothNoMask_Threaded (NULL, readout->image, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA, 0.25);
+        psImageConvolveSetThreads(oldThreads);
         psLogMsg ("psphot.kron", PS_LOG_INFO, "smoothed image %f sec\n", psTimerMark ("psphot.kron.smooth"));
 
         // remove the sources
+        psTimerStart ("psphot.kron.smooth");
         psphotRemoveAllSourcesReadout( config, view, filerule, index, recipe, false );
+        psLogMsg ("psphot.kron", PS_LOG_INFO, "removed %ld raw sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth"));
+
         // Now subtract smooth versions of the sources from the smoothed image
         psTimerStart ("psphot.kron.smooth.sources");

trunk/psphot/src/psphotModelBackground.c

@@ -28 +28 @@
 }
 
+// track background cell failures for log reporting purposes
+static int nFailures = 0;
 
 // Generate the background model
@@ -33 +35 @@
 // linear interpolation to generate full-scale model
 //
-// NOTE that the 'analysis' metedata pass in here is used to store the binning information.
+// NOTE that the 'analysis' metedata passed in here is used to store the binning information.
 // This may be the analysis for this readout, but it may be the analysis for the pmFPAfile
 // corresponding to the model.  Other information about the background model is saved on the
@@ -109 +111 @@
 
     const psStatsOptions statsOption = statsOptionLocation | statsOptionWidth;
-    psStats *stats = psStatsAlloc (statsOption);
+    psStats *statsDefaults = psStatsAlloc (statsOption);
 
     // set range for old-version of sky statistic
     if (statsOptionLocation & PS_STAT_ROBUST_QUARTILE) {
-        stats->min = 0.25;
-        stats->max = 0.75;
+        statsDefaults->min = 0.25;
+        statsDefaults->max = 0.75;
     }
 
     // set user-option for number of pixels per region
-    stats->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
+    statsDefaults->nSubsample = psMetadataLookupF32 (&status, recipe, "IMSTATS_NPIX");
     if (!status) {
-        stats->nSubsample = 1000;
+        statsDefaults->nSubsample = 1000;
     }
 
     // optionally set the binsize
-    stats->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS");
+    statsDefaults->binsize = psMetadataLookupF32 (&status, recipe, "SKY_HISTOGRAM_BINS");
     if (status) {
-        stats->options |= PS_STAT_USE_BINSIZE;
+        statsDefaults->options |= PS_STAT_USE_BINSIZE;
     }
 
     // optionally set the sigma clipping
-    stats->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA");
+    statsDefaults->clipSigma = psMetadataLookupF32 (&status, recipe, "SKY_CLIP_SIGMA");
     if (!status) {
-        if (stats->options & PS_STAT_FITTED_MEAN) {
-            stats->clipSigma = 1.0;
+        if (statsDefaults->options & PS_STAT_FITTED_MEAN) {
+            statsDefaults->clipSigma = 1.0;
         } else {
-            stats->clipSigma = 3.0;
+            statsDefaults->clipSigma = 3.0;
         }
     }
-
-    // stats is not initialized by psStats???  use this to save the input options
-    // XXX re-work this to use the new psStatsInit function
-    psStats *statsDefaults = psStatsAlloc (statsOption);
-    *statsDefaults = *stats;
 
     // we save the binning structure for use in psphotMagnitudes
@@ -154 +151 @@
 
     // XXXX we can thread this here by running blocks in parallel
-
-    int nFailures = 0;
     psImageBackgroundInit();
+    nFailures = 0; // reset for this pass
 
     // we have Nx * Ny model points, but we can use a window which is larger (or smaller) than
@@ -163 +159 @@
 
     // measure clipped median for subimages
-    psRegion ruffRegion = psRegionSet (0,0,0,0);
-    psRegion fineRegion = psRegionSet (0,0,0,0);
     for (int iy = 0; iy < model->numRows; iy++) {
         for (int ix = 0; ix < model->numCols; ix++) {
-
-            // convert the ruff grid cell to the equivalent fine grid cell
-            ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample, 
-                                      iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample);
-            fineRegion = psImageBinningSetFineRegion (binning, ruffRegion);
-            fineRegion = psRegionForImage (image, fineRegion);
-
-            psImage *subset  = psImageSubset (image, fineRegion);
-            if (!subset->numCols || !subset->numRows) {
-                psFree (subset);
-                continue;
+            
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_MODEL_BACKGROUND");
+
+            psArrayAdd(job->args, 1, image);
+            psArrayAdd(job->args, 1, mask);
+            psArrayAdd(job->args, 1, binning);
+            psArrayAdd(job->args, 1, rng);
+            psArrayAdd(job->args, 1, statsDefaults);
+
+            psArrayAdd(job->args, 1, modelData);
+            psArrayAdd(job->args, 1, modelStdevData);
+
+            PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+
+            PS_ARRAY_ADD_SCALAR(job->args, ix, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, iy, PS_TYPE_S32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, dXsample, PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, dYsample, PS_TYPE_F32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, statsOptionLocation, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, statsOptionWidth, PS_TYPE_S32);
+
+            PS_ARRAY_ADD_SCALAR(job->args, NAN, PS_TYPE_F32); // this is used as a return value for dQvalue
+
+# if (1)
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                return NULL;
             }
-            psImage *submask = psImageSubset (mask, fineRegion);
-
-            // reset the default values
-            *stats = *statsDefaults;
-
-            // Use the selected background statistic for the first pass
-            // If it fails, fall back on the "ROBUST_MEDIAN" version
-            // If both fail, set the pixel to NAN and (below) interpolate
-            // XXX psImageBackground will probably be renamed psImageStats
-            // XXX don't bother trying if there are no valid pixels...
-
-            psVector *sample = NULL;
-
-            // turn on stats tracing in desired cells
-            # if (0)
-            psMetadata *plots = psMetadataLookupPtr (&status, recipe, "DIAGNOSTIC.PLOTS");
-            assert (plots);
-
-            int xPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.X");
-            assert (status);
-            int yPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.Y");
-            assert (status);
-
-            bool gotX = (xPlot < 0) || (xPlot == ix);
-            bool gotY = (yPlot < 0) || (yPlot == iy);
-
-            if (gotX && gotY) {
-                (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6);
-                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
-            } else {
-                (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0);
-                (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
+# else
+            if (!psphotModelBackground_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Failure to model background.");
+                return NULL;
             }
-            # endif
-
-            if (psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
-                if (stats->options & PS_STAT_ROBUST_QUARTILE) {
-                    modelData[iy][ix] = stats->robustMedian;
-                } else {
-                    modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation);
-                }
-                modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
-                // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
-                psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
-
-                // supply sample to plotting routing
-                psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
-            } else {
-                psStatsOptions currentOptions = stats->options;
-                stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
-                if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
-                    if ((nFailures < 3) || (nFailures % 100 == 0)) {
-                        psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
-                                  "(%dx%d, (row0,col0) = (%d,%d)",
-                                  subset->numRows, subset->numCols, subset->row0, subset->col0);
-                    }
-                    nFailures ++;
-                    modelData[iy][ix] = modelStdevData[iy][ix] = NAN;
-                } else {
-                    modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
-                    modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV);
-
-                    // supply sample to plotting routing
-                    psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
-                }
-                // drop errors caused by psImageBackground failures
-                // XXX we probably should trap and exit on serious failures
-                psErrorClear();
-                stats->options = currentOptions;
-            }
-            psFree(sample);
-            modelData[iy][ix] += SKY_BIAS;
-            psFree (subset);
-            psFree (submask);
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                psScalar *scalar = job->args->data[14];
+                float dQvalue = scalar->data.F32;
+                psVectorAppend (dQ, dQvalue);
+            }
+            psFree(job);
+# endif
         }
+    }
+
+    // wait for the threads to finish and manage results
+    if (!psThreadPoolWait (false, true)) {
+        psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+        return NULL;
+    }
+
+    // we have only supplied one type of job, so we can assume the types here
+    psThreadJob *job = NULL;
+    while ((job = psThreadJobGetDone()) != NULL) {
+        if (job->args->n < 1) {
+            fprintf (stderr, "error with job\n");
+        } else {
+            psScalar *scalar = job->args->data[14];
+            float dQvalue = scalar->data.F32;
+            psVectorAppend (dQ, dQvalue);
+        }
+        psFree(job);
     }
 
@@ -302 +274 @@
     if (Count == 0) {
         psError (PSPHOT_ERR_DATA, true, "failed to build background image");
-        psFree(stats);
         psFree(statsDefaults);
         psFree(binning);
@@ -318 +289 @@
             modelData[iy][ix] = Value;
             modelStdevData[iy][ix] = ValueStdev;
+        }
+    }
+
+    // apply artificial offset if desired
+    for (int iy = 0; iy < model->numRows; iy++) {
+        for (int ix = 0; ix < model->numCols; ix++) {
+            modelData[iy][ix] += SKY_BIAS;
         }
     }
@@ -364 +342 @@
     psFree(dQ);
 
-    psFree(stats);
     psFree(statsDefaults);
     psFree(binning);
@@ -433 +410 @@
     return true;
 }
+
+bool psphotModelBackground_Threaded (psThreadJob *job) {
+
+    psImage *image          = job->args->data[0];
+    psImage *mask           = job->args->data[1];
+
+    psImageBinning *binning = job->args->data[2];
+
+    psRandom *rng           = job->args->data[3];
+    psStats *statsDefaults  = job->args->data[4];
+
+    psF32 **modelData       = job->args->data[5];
+    psF32 **modelStdevData  = job->args->data[6];
+
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
+
+    int ix                  = PS_SCALAR_VALUE(job->args->data[8],S32);
+    int iy                  = PS_SCALAR_VALUE(job->args->data[9],S32);
+
+    float dXsample          = PS_SCALAR_VALUE(job->args->data[10],F32);
+    float dYsample          = PS_SCALAR_VALUE(job->args->data[11],F32);
+
+    psStatsOptions statsOptionLocation = PS_SCALAR_VALUE(job->args->data[12],S32);
+    psStatsOptions statsOptionWidth    = PS_SCALAR_VALUE(job->args->data[13],S32);
+
+    // convert the ruff grid cell to the equivalent fine grid cell
+    psRegion ruffRegion = psRegionSet (ix + 0.5 - 0.5*dXsample, ix + 0.5 + 0.5*dXsample, iy + 0.5 - 0.5*dYsample, iy + 0.5 + 0.5*dYsample);
+    psRegion fineRegion = psImageBinningSetFineRegion (binning, ruffRegion);
+    fineRegion = psRegionForImage (image, fineRegion);
+
+    psImage *subset  = psImageSubset (image, fineRegion);
+    if (!subset->numCols || !subset->numRows) {
+        psFree (subset);
+        return false; // XXX do we / should we fail on this?
+    }
+    psImage *submask = psImageSubset (mask, fineRegion);
+
+    psStats *stats = psStatsAlloc (PS_STAT_NONE);
+    *stats = *statsDefaults;
+
+    // Use the selected background statistic for the first pass
+    // If it fails, fall back on the "ROBUST_MEDIAN" version
+    // If both fail, set the pixel to NAN and (later) interpolate
+
+    psVector *sample = NULL;
+    float dQvalue = NAN;
+
+    if (psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
+        if (stats->options & PS_STAT_ROBUST_QUARTILE) {
+            modelData[iy][ix] = stats->robustMedian;
+        } else {
+            modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation);
+        }
+        modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
+
+        // fprintf (stderr, "dQ : %f - %f - %f = %f\n", stats->robustLQ, stats->robustMedian, stats->robustUQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
+        // XXX this operation is not thread safe -- move out somehow...
+        // psVectorAppend (dQ, stats->robustUQ + stats->robustLQ - 2*stats->robustMedian);
+        dQvalue = stats->robustUQ + stats->robustLQ - 2*stats->robustMedian;
+        // return dQvalue to main thread
+
+        // supply sample to plotting routing
+        // only allow in a non-threaded context -- can plot more than one cell
+        // psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
+    } else {
+        // psStatsOptions currentOptions = stats->options;
+        stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
+        if (!psImageBackground(stats, &sample, subset, submask, maskVal, rng)) {
+            if ((nFailures < 3) || (nFailures % 100 == 0)) {
+                psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
+                          "(%dx%d, (row0,col0) = (%d,%d)",
+                          subset->numRows, subset->numCols, subset->row0, subset->col0);
+            }
+            nFailures ++; // static
+            modelData[iy][ix] = modelStdevData[iy][ix] = NAN;
+        } else {
+            modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
+            modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV);
+
+            // supply sample to plotting routing
+            // only allow in a non-threaded context -- can plot more than one cell
+            // psphotDiagnosticPlots (config, "IMAGE.BACKGROUND.CELL.HISTOGRAM", ix, iy, modelData[iy][ix], modelStdevData[iy][ix], sample);
+        }
+        // drop errors caused by psImageBackground failures
+        // NOTE : psStats raises errors when it cannot find a solution; this 
+        // probably should not be errors but exit stats info only, but c'est la code
+        // XXX we probably should trap and exit on serious failures
+        psErrorClear();
+        // stats->options = currentOptions; // this is not needed (set by init above)
+    }
+    psFree (stats);
+    psFree (sample);
+    psFree (subset);
+    psFree (submask);
+
+    // return the dQvalue to the calling thread
+    psScalar *scalar = job->args->data[14];
+    scalar->data.F32 = dQvalue;
+
+    return true;
+}
+
+// code for in-line plotting from above
+// turn on stats tracing in desired cells
+// XXX this code may need to be re-worked for a threaded context
+# if (0)
+    psMetadata *plots = psMetadataLookupPtr (&status, recipe, "DIAGNOSTIC.PLOTS");
+    assert (plots);
+
+    int xPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.X");
+    assert (status);
+    int yPlot = psMetadataLookupS32 (&status, plots, "IMAGE.BACKGROUND.CELL.HISTOGRAM.Y");
+    assert (status);
+
+    bool gotX = (xPlot < 0) || (xPlot == ix);
+    bool gotY = (yPlot < 0) || (yPlot == iy);
+
+    if (gotX && gotY) {
+        (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 6);
+        (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 6);
+    } else {
+        (void) psTraceSetLevel ("psLib.math.vectorFittedStats", 0);
+        (void) psTraceSetLevel ("psLib.math.vectorRobustStats", 0);
+    }
+# endif

trunk/psphot/src/psphotModelGroupInit.c

@@ -4 +4 @@
 // psModule/src/objects/models
 
-# include "models/pmModel_TEST1.c"
-# include "models/pmModel_STRAIL.c"
+// # include "models/pmModel_TEST1.c"
+// # include "models/pmModel_STRAIL.c"
 
 static pmModelClass userModels[] = {
-    {"PS_MODEL_TEST1", 7, pmModelFunc_TEST1,  pmModelFlux_TEST1,  pmModelRadius_TEST1,  pmModelLimits_TEST1,  pmModelGuess_TEST1, pmModelFromPSF_TEST1, pmModelParamsFromPSF_TEST1, pmModelFitStatus_TEST1, NULL},
-    {"PS_MODEL_STRAIL", 9, pmModelFunc_STRAIL,  pmModelFlux_STRAIL,  pmModelRadius_STRAIL,  pmModelLimits_STRAIL,  pmModelGuess_STRAIL, pmModelFromPSF_STRAIL, pmModelParamsFromPSF_STRAIL, pmModelFitStatus_STRAIL, NULL},
+    // {"PS_MODEL_TEST1", 7, pmModelFunc_TEST1,  pmModelFlux_TEST1,  pmModelRadius_TEST1,  pmModelLimits_TEST1,  pmModelGuess_TEST1, pmModelFromPSF_TEST1, pmModelParamsFromPSF_TEST1, pmModelFitStatus_TEST1, NULL},
+    // {"PS_MODEL_STRAIL", 9, pmModelFunc_STRAIL,  pmModelFlux_STRAIL,  pmModelRadius_STRAIL,  pmModelLimits_STRAIL,  pmModelGuess_STRAIL, pmModelFromPSF_STRAIL, pmModelParamsFromPSF_STRAIL, pmModelFitStatus_STRAIL, NULL},
 };
 

trunk/psphot/src/psphotModelTestReadout.c

@@ -161 +161 @@
     // get the initial model parameter guess
     pmModel *model = pmSourceModelGuess (source, modelType, maskVal, markVal);
+    if (!model) {
+      fprintf (stderr, "failed to generate model guess\n");
+      exit (2);
+    }
+
     source->modelEXT = model;
 
@@ -186 +191 @@
         fprintf (stderr, "guess: %f @ (%f, %f)\n", params[6]*180/M_PI, params[4], params[5]);
     } else {
-        // list model input shape
-        psEllipseShape shape;
-        shape.sx  = 1.4 / model->params->data.F32[4];
-        shape.sy  = 1.4 / model->params->data.F32[5];
-        shape.sxy = model->params->data.F32[6];
-        axes = psEllipseShapeToAxes (shape, 20.0);
-        fprintf (stderr, "guess: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
+        bool useReff = pmModelUseReff (modelType);
+        pmModelParamsToAxes (&axes, params[PM_PAR_SXX], params[PM_PAR_SXY], params[PM_PAR_SYY], useReff);
+        fprintf (stderr, "guess: %f @ (%f, %f) : %f\n", axes.theta*180/M_PI, axes.major, axes.minor, params[PM_PAR_SXY]);
     }
 
@@ -217 +218 @@
     fitOptions->covarFactor   = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
 
+    bool chisqConvergence = psMetadataLookupBool (&status, recipe, "LMM_FIT_CHISQ_CONVERGENCE"); // Fit tolerance
+    if (!status) chisqConvergence = true;
+    fitOptions->chisqConvergence = chisqConvergence;
+
+    int gainFactorMode = psMetadataLookupS32 (&status, recipe, "LMM_FIT_GAIN_FACTOR_MODE"); // Fit tolerance
+    if (!status) gainFactorMode = 0;
+    fitOptions->gainFactorMode = gainFactorMode;
+
     if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
         fitOptions->mode = PM_SOURCE_FIT_NO_INDEX;
@@ -231 +240 @@
         pmSourceModelGuessPCM (pcm, source, maskVal, markVal);
       }
+
+      // if we provide a test guess value we want to use that value!
+      for (int i = 0; i < nParams; i++) {
+          if (i == PM_PAR_XPOS) continue;
+          if (i == PM_PAR_YPOS) continue;
+
+          char name[32];
+          sprintf (name, "TEST_FIT_PAR%d", i);
+          float value = psMetadataLookupF32 (&status, recipe, name);
+          if (status && isfinite (value)) {
+              params[i] = value;
+          }
+      }
+
       pmPCMupdate(pcm, source, fitOptions, model);
       pmSourceFitPCM (pcm, source, fitOptions, maskVal, markVal, psfSize);
@@ -258 +281 @@
     }
 
+    if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
+        fprintf (stderr, "result: %f @ (%f, %f)\n", params[6]*180/M_PI, params[4], params[5]);
+    } else {
+        bool useReff = pmModelUseReff (modelType);
+        pmModelParamsToAxes (&axes, params[PM_PAR_SXX], params[PM_PAR_SXY], params[PM_PAR_SYY], useReff);
+        fprintf (stderr, "result: %f @ (%f, %f) : %f\n", axes.theta*180/M_PI, axes.major, axes.minor, params[PM_PAR_SXY]);
+    }
+
     // write out
     psphotSaveImage (NULL, source->pixels, "resid.fits");

trunk/psphot/src/psphotOutput.c

@@ -379 +379 @@
         psF32 *outPar = source->modelEXT->params->data.F32;
 
-        psEllipseShape shape;
-
-        shape.sx  = outPar[PM_PAR_SXX] / M_SQRT2;
-        shape.sy  = outPar[PM_PAR_SYY] / M_SQRT2;
-        shape.sxy = outPar[PM_PAR_SXY];
-
-        psEllipsePol pol = pmPSF_ModelToFit (outPar);
+        bool useReff = pmModelUseReff (source->modelEXT->type);
+
+        psEllipseAxes axes1;
+        pmModelParamsToAxes (&axes1, outPar[PM_PAR_SXX], outPar[PM_PAR_SXY], outPar[PM_PAR_SYY], useReff);
+
+        psEllipsePol pol = pmPSF_ModelToFit (outPar, useReff);
         inPar[PM_PAR_E0] = pol.e0;
         inPar[PM_PAR_E1] = pol.e1;
         inPar[PM_PAR_E2] = pol.e2;
-        pmPSF_FitToModel (inPar, 0.1);
-
-        psEllipseAxes axes1 = psEllipseShapeToAxes (shape, 20.0);
+        pmPSF_FitToModel (inPar, 0.1, useReff);
+
         psEllipseAxes axes2 = psEllipsePolToAxes(pol, 0.1);
 

trunk/psphot/src/psphotRadialPlot.c

@@ -15 +15 @@
     if (nCount > 00) {
         if (*kapa != 0) {
-            KiiClose (*kapa);
+            KapaClose (*kapa);
             *kapa = 0;
         }

trunk/psphot/src/psphotReplaceUnfit.c

@@ -77 +77 @@
     // maskVal is used to test for rejected pixels, and must include markVal
     maskVal |= markVal;
+
+    int NpixTotal = 0;
 
     for (int i = 0; i < sources->n; i++) {
@@ -99 +101 @@
 
         pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+        NpixTotal += source->pixels->numCols * source->pixels->numRows;
     }
 
     psphotVisualShowImage(readout);
-    psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec\n", sources->n, psTimerMark ("psphot.replace"));
+    psLogMsg ("psphot.replace", PS_LOG_INFO, "replaced models for %ld objects: %f sec (%d pixels)\n", sources->n, psTimerMark ("psphot.replace"), NpixTotal);
     return true;
 }

trunk/psphot/src/psphotSetThreads.c

@@ -5 +5 @@
     psThreadTask *task = NULL;
 
+    pmPSFThreads ();
+
+    task = psThreadTaskAlloc("PSPHOT_MODEL_BACKGROUND", 15);
+    task->function = &psphotModelBackground_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
     task = psThreadTaskAlloc("PSPHOT_GUESS_MODEL", 5);
     task->function = &psphotGuessModel_Threaded;
+    psThreadTaskAdd(task);
+    psFree(task);
+
+    task = psThreadTaskAlloc("PSPHOT_ADD_NOISE", 6);
+    task->function = &psphotAddOrSubNoise_Threaded;
     psThreadTaskAdd(task);
     psFree(task);
@@ -40 +52 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 18);
+    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 16);
     task->function = &psphotExtendedSourceFits_Threaded;
     psThreadTaskAdd(task);

trunk/psphot/src/psphotSourceFits.c

@@ -3 +3 @@
 // given a source with an existing modelPSF, attempt a full PSF fit, subtract if successful
 
+static int NfitBlend = 0;
+
 static int NfitPSF = 0;
-static int NfitBlend = 0;
+static int NfitIterPSF = 0;
+
 static int NfitDBL = 0;
+static int NfitIterDBL = 0;
+static int NfitPixDBL = 0;
+
 static int NfitEXT = 0;
+static int NfitIterEXT = 0;
+static int NfitPixEXT = 0;
+
 static int NfitPCM = 0;
+static int NfitIterPCM = 0;
+static int NfitPixPCM = 0;
 
 bool psphotFitInit (int nThreads) {
@@ -17 +28 @@
 bool psphotFitSummary () {
 
-    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf, %5d blend, %5d ext, %5d dbl : %6.2f sec\n",
-              NfitPSF, NfitBlend, NfitEXT, NfitDBL, psTimerMark ("psphot.fits"));
+    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d psf (%d iter), %5d blend: %5d ext (%d iter, %d pix), %5d dbl (%d iter, %d pix) : %6.2f sec\n",
+              NfitPSF, NfitIterPSF, NfitBlend, NfitEXT, NfitIterEXT, NfitPixEXT, NfitDBL, NfitIterDBL, NfitPixDBL, psTimerMark ("psphot.fits"));
+    return true;
+}
+
+bool psphotFitSummaryExtended () {
+
+    psLogMsg ("psphot.pspsf", PS_LOG_INFO, "fitted %5d ext (%d iter, %d pix), %5d pcm (%d iter, %d pix)\n",
+              NfitEXT, NfitIterEXT, NfitPixEXT, NfitPCM, NfitIterPCM, NfitPixPCM);
+    return true;
+}
+
+bool psphotFitInitExtended () {
+    NfitEXT = 0;
+    NfitIterEXT = 0;
+    NfitPixEXT = 0;
+    NfitPCM = 0;
+    NfitIterPCM = 0;
+    NfitPixPCM = 0;
     return true;
 }
@@ -169 +197 @@
     options.mode = PM_SOURCE_FIT_PSF;
     pmSourceFitModel (source, PSF, &options, maskVal);
+    NfitIterPSF += PSF->nIter;
 
     if (!isfinite(PSF->params->data.F32[PM_PAR_I0])) psAbort("nan in fit");
@@ -268 +297 @@
     { 
         // XXX need to handle failures better here
-        EXT = psphotFitEXT (readout, source, fitOptions, modelTypeEXT, maskVal, markVal);
+        EXT = psphotFitEXT (NULL, readout, source, fitOptions, modelTypeEXT, maskVal, markVal);
         if (!EXT) goto escape;
         if (!isfinite(EXT->params->data.F32[PM_PAR_I0])) goto escape;
@@ -440 +469 @@
     options.mode = PM_SOURCE_FIT_PSF;
     pmSourceFitSet (source, modelSet, &options, maskVal);
+    NfitIterDBL += DBL->nIter;
+    NfitPixDBL += DBL->nDOF;
+
     return (modelSet);
 }
 
-pmModel *psphotFitEXT (pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
+pmModel *psphotFitEXT (pmModel *guessModel, pmReadout *readout, pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, psImageMaskType markVal) {
 
     if ((source->moments->Mxx < 1e-3) || (source->moments->Myy < 1e-3)) {
@@ -457 +489 @@
 
     // use the source moments, etc to guess basic model parameters
-    pmModel *model = pmSourceModelGuess (source, modelType, maskVal, markVal);
+    pmModel *model = guessModel ? guessModel : pmSourceModelGuess (source, modelType, maskVal, markVal);
     if (!model) {
         psTrace ("psphot", 5, "failed to generate a model for source: moments: %f %f\n", source->moments->Mxx, source->moments->Myy);
@@ -507 +539 @@
 
     pmSourceFitModel (source, model, &options, maskVal);
+    NfitIterEXT += model->nIter;
+    NfitPixEXT += model->nDOF;
 
 # if (PS_TRACE_ON) 
@@ -587 +621 @@
     // psTraceSetLevel("psLib.math.psMinimizeLMChi2", 5);
     pmSourceFitPCM (pcm, source, &options, maskVal, markVal, psfSize);  // NOTE : 1687 allocs in here
+    NfitIterPCM += pcm->modelConv->nIter;
+    NfitPixPCM += pcm->modelConv->nDOF;
     if (TIMING) { t5 = psTimerMark ("psphotFitPCM"); }
 

trunk/psphot/src/psphotVisual.c

@@ -29 +29 @@
 bool psphotVisualClose(void)
 {
-    if(kapa1 != -1) KiiClose(kapa1);
-    if(kapa2 != -1) KiiClose(kapa2);
-    if(kapa3 != -1) KiiClose(kapa3);
+    if(kapa1 != -1) KapaClose(kapa1);
+    if(kapa2 != -1) KapaClose(kapa2);
+    if(kapa3 != -1) KapaClose(kapa3);
     return true;
 }