Changeset 4630

Timestamp:

Jul 27, 2005, 3:26:57 PM (21 years ago)

Author:

eugene

Message:

various tests and changes

Location:

trunk/psphot/src

Files:

• LocalSky.c (modified) (7 diffs)
• apply_psf_model.c (modified) (6 diffs)
• choose_psf_model.c (modified) (1 diff)
• fit_galaxies.c (modified) (4 diffs)
• mark_psf_source.c (modified) (5 diffs)
• mark_psf_sources.c (modified) (1 diff)
• onesource.c (modified) (2 diffs)
• psphot-utils.c (modified) (4 diffs)
• psphot.h (modified) (2 diffs)
• setup.c (modified) (1 diff)
• subtract_psf_source.c (modified) (2 diffs)

trunk/psphot/src/LocalSky.c

@@ -62 +62 @@
 bool pmSourceFitModel_EAM(psSource *source,
                           psModel *model,
-                          const bool PSF,
-                          float SOFT)
+                          const bool PSF)
 {
     PS_PTR_CHECK_NULL(source, false);
@@ -74 +73 @@
     psBool onPic     = true;
     psBool rc        = true;
-    psF32  Ro;
+    psF32  Ro, ymodel;
 
 
@@ -81 +80 @@
 
     psModelFunc modelFunc = psModelFunc_GetFunction (model->type);
+    psModelLimits modelLimits = psModelLimits_GetFunction (model->type);
 
     psVector *params = model->params;
@@ -86 +86 @@
     psVector *paramMask = NULL;
 
+    psVector *beta_lim = NULL;
+    psVector *params_min = NULL;
+    psVector *params_max = NULL;
+
     int nParams = PSF ? params->n - 4 : params->n;
-    psF32 Xo = params->data.F32[2];
-    psF32 Yo = params->data.F32[3];
+    psF32 So = params->data.F32[0];
 
     // find the number of valid pixels
@@ -122 +125 @@
                 y->data.F32[tmpCnt] = source->pixels->data.F32[i][j];
 
-                // XXX just for test purposes, use the raw radius.  a better
-                // choice is to ask what is the value of z and scale by that
-                Ro = hypot ((Xo - coord->data.F32[0]), (Yo - coord->data.F32[1]));
-
-                // XXX enhance the noise (doubled if R = 10 pix, etc. note the square)
-                yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * (1 + PS_SQR(Ro/SOFT)));
+                ymodel = modelFunc (NULL, model->params, coord);
+                
+                // this test enhances the noise based on deviation from the model flux
+                Ro = 1.0 + fabs (y->data.F32[tmpCnt] - ymodel) / sqrt(PS_SQR(ymodel - So) + PS_SQR(So));
+                yErr->data.F32[tmpCnt] = sqrt (source->noise->data.F32[i][j] * Ro);
                 // XXX EAM : note the wasted effort: we carry dY^2, take sqrt(), then 
                 //           the minimization function calculates sq()
@@ -147 +149 @@
         paramMask->data.U8[i] = 1;
       }
-    }       
-
-    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    }  
+
+    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
+    modelLimits (&beta_lim, &params_min, &params_max);
+    for (int i = 0; i < params->n; i++) {
+        covar->data.F64[0][i] = beta_lim->data.F32[i];
+        covar->data.F64[1][i] = params_min->data.F32[i];
+        covar->data.F64[2][i] = params_max->data.F32[i];
+    }
 
     psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
+    fitStatus = psMinimizeLMChi2_EAM(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
     for (int i = 0; i < dparams->n; i++) {
         if ((paramMask != NULL) && paramMask->data.U8[i]) continue;
@@ -326 +334 @@
 }
 
+bool pmModelFitStatus (psModel *model) {
+
+    bool status;
+
+    psModelFitStatusFunc statusFunc = psModelFitStatusFunc_GetFunction (model->type);
+    status = statusFunc (model);
+
+    return (status);
+}
+
+// XXX EAM this implementation of MinLM includes limits on params & dparams
+psBool psMinimizeLMChi2_EAM(psMinimization *min,
+                            psImage *covar,
+                            psVector *params,
+                            const psVector *paramMask,
+                            const psArray *x,
+                            const psVector *y,
+                            const psVector *yErr,
+                            psMinimizeLMChi2Func func)
+{
+    PS_PTR_CHECK_NULL(min, NULL);
+    PS_VECTOR_CHECK_NULL(params, NULL);
+    PS_VECTOR_CHECK_EMPTY(params, NULL);
+    PS_PTR_CHECK_NULL(x, NULL);
+    PS_VECTOR_CHECK_NULL(y, NULL);
+    PS_VECTOR_CHECK_EMPTY(y, NULL);
+    PS_VECTOR_CHECK_SIZE_EQUAL(x, y, NULL);
+    PS_PTR_CHECK_NULL(func, NULL);
+
+    // this function has test and current values for several things
+    // the current best value is in lower case
+    // the next guess value is in upper case
+
+    // allocate internal arrays (current vs Guess)
+    psImage *alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    psImage *Alpha   = psImageAlloc (params->n, params->n, PS_TYPE_F64);
+    psVector *beta   = psVectorAlloc (params->n, PS_TYPE_F64);
+    psVector *Beta   = psVectorAlloc (params->n, PS_TYPE_F64);
+    psVector *Params = psVectorAlloc (params->n, PS_TYPE_F32);
+    psVector *dy     = NULL;
+    psF64 Chisq = 0.0;
+    psF64 lambda = 0.001;
+
+    psVector *beta_lim = NULL;
+    psVector *param_min = NULL;
+    psVector *param_max = NULL;
+
+    // if we are provided a covar image, we expect to find these three vectors in first three rows
+    if (covar != NULL) {
+        beta_lim  = psVectorAlloc (params->n, PS_TYPE_F32);
+        param_min = psVectorAlloc (params->n, PS_TYPE_F32);
+        param_max = psVectorAlloc (params->n, PS_TYPE_F32);
+        for (int i = 0; i < params->n; i++) {
+            beta_lim->data.F32[i] = covar->data.F64[0][i];
+            param_min->data.F32[i] = covar->data.F64[1][i];
+            param_max->data.F32[i] = covar->data.F64[2][i];
+        }
+        psImageRecycle (covar, params->n, params->n, PS_TYPE_F64);
+    }
+        
+
+    // why is this needed here??? the initial guess on params is provided by the user
+    Params = psVectorCopy (Params, params, PS_TYPE_F32);
+
+    // the user provides the error or NULL.  we need to convert
+    // to appropriate weights
+    dy = psVectorAlloc (y->n, PS_TYPE_F32);
+    if (yErr != NULL) {
+        for (int i = 0; i < dy->n; i++) {
+            if (yErr->data.F32[i] == 0.0) {
+                dy->data.F32[i] = 1.0;
+                // mask this?  bad pixel, obviously...
+            } else {
+                dy->data.F32[i] = 1.0 / PS_SQR (yErr->data.F32[i]);
+            }
+        }
+    } else {
+        for (int i = 0; i < dy->n; i++) {
+            dy->data.F32[i] = 1.0;
+        }
+    }
+
+    // calculate initial alpha and beta, set chisq (min->value)
+    min->value = p_psMinLM_SetABX (alpha, beta, params, paramMask, x, y, dy, func);
+    if (isnan(min->value)) {
+        min->iter = min->maxIter;
+        return (false);
+    }
+    # ifndef PS_NO_TRACE
+    // dump some useful info if trace is defined
+    if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
+        p_psImagePrint  (psTraceGetDestination(), alpha, "alpha guess");
+        p_psVectorPrint (psTraceGetDestination(), beta, "beta guess");
+        p_psVectorPrint (psTraceGetDestination(), params, "params guess");
+    }
+    if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) {
+        p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
+    }
+    # endif /* PS_NO_TRACE */
+
+    // iterate until the tolerance is reached, or give up
+    while ((min->lastDelta > min->tol) && (min->iter < min->maxIter)) {
+
+        // set a new guess for Alpha, Beta, Params
+        p_psMinLM_GuessABP_EAM (Alpha, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, lambda);
+
+        // measure linear model prediction
+        psF64 dLinear = p_psMinLM_dLinear (Beta, beta, lambda);
+
+        # ifndef PS_NO_TRACE
+        // dump some useful info if trace is defined
+        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
+            p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
+            p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
+            p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
+        }
+        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") == 4) {
+            p_psVectorPrintRow (psTraceGetDestination(), Params, "params guess");
+        }
+        # endif /* PS_NO_TRACE */
+
+        // calculate Chisq for new guess, update Alpha & Beta
+        Chisq = p_psMinLM_SetABX (Alpha, Beta, Params, paramMask, x, y, dy, func);
+
+        // XXX EAM alternate convergence criterion:
+        // compare the delta (min->value - Chisq) with the
+        // expected delta from the linear model (dLinear)
+        // accept new guess (if improvement), or increase lambda
+        psF64 rho = (min->value - Chisq) / dLinear;
+
+        psTrace (".psLib.dataManip.psMinimizeLMChi2", 4, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value, Chisq, min->lastDelta, rho);
+        # ifndef PS_NO_TRACE
+        // dump some useful info if trace is defined
+        if (psTraceGetLevel (".psLib.dataManip.psMinimizeLMChi2") >= 5) {
+            p_psImagePrint  (psTraceGetDestination(), Alpha, "alpha guess");
+            p_psVectorPrint (psTraceGetDestination(), Beta, "beta guess");
+            p_psVectorPrint (psTraceGetDestination(), Params, "params guess");
+        }
+        # endif /* PS_NO_TRACE */
+
+        /* if (Chisq < min->value) {  */
+        if (rho > 0.0) {  
+            min->lastDelta = (min->value - Chisq) / (dy->n - params->n);
+            min->value = Chisq;
+            alpha  = psImageCopy (alpha, Alpha, PS_TYPE_F64);
+            beta   = psVectorCopy (beta, Beta, PS_TYPE_F64);
+            params = psVectorCopy (params, Params, PS_TYPE_F32);
+            lambda *= 0.1;
+        } else {
+            lambda *= 10.0;
+        }
+        min->iter ++;
+    }
+    psTrace (".psLib.dataManip.psMinimizeLMChi2", 3, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
+
+    // construct & return the covariance matrix (if requested)
+    if (covar != NULL) {
+      p_psMinLM_GuessABP_EAM (covar, Beta, Params, alpha, beta, params, paramMask, beta_lim, param_min, param_max, 0.0);
+    }      
+
+    // free the internal temporary data
+    psFree (alpha);
+    psFree (Alpha);
+    psFree (beta);
+    psFree (Beta);
+    psFree (Params);
+    psFree (dy);
+
+    if (min->iter == min->maxIter) {
+      return (false);
+    } 
+    return (true);
+}
+
+// XXX EAM : can we use static copies of LUv, LUm, A?
+psBool p_psMinLM_GuessABP_EAM (psImage  *Alpha,
+                               psVector *Beta,
+                               psVector *Params,
+                               const psImage  *alpha,
+                               const psVector *beta,
+                               const psVector *params,
+                               const psVector *paramMask,
+                               const psVector *beta_lim,
+                               const psVector *params_min,
+                               const psVector *params_max,
+                               psF64 lambda)
+{
+
+    # define USE_LU_DECOMP 1
+    # if (USE_LU_DECOMP)
+    psVector *LUv = NULL;
+    psImage  *LUm = NULL;
+    psImage  *A   = NULL;
+    psF32    det;
+
+    // LU decomposition version
+    psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using LUD version\n");
+
+    // set new guess values (creates matrix A)
+    A = psImageCopy (NULL, alpha, PS_TYPE_F64);
+    for (int j = 0; j < params->n; j++) {
+        if ((paramMask != NULL) && (paramMask->data.U8[j])) continue;
+        A->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
+    }
+
+    // solve A*beta = Beta (Alpha = 1/A)
+    // these operations do not modify the input values (creates LUm, LUv)
+    LUm   = psMatrixLUD (NULL, &LUv, A);
+    Beta  = psMatrixLUSolve (Beta, LUm, beta, LUv);
+    Alpha = psMatrixInvert (Alpha, A, &det);
+
+    # else
+    // gauss-jordan version
+    psTrace (".psLib.dataManip.psMinLM_GuessABP", 5, "using Gauss-J version");
+
+    // set new guess values (creates matrix A)
+    Beta = psVectorCopy (Beta, beta, PS_TYPE_F64);
+    Alpha = psImageCopy (Alpha, alpha, PS_TYPE_F64);
+    for (int j = 0; j < params->n; j++) {
+        if ((paramMask != NULL) && (paramMask->data.U8[j])) continue;
+        Alpha->data.F64[j][j] = alpha->data.F64[j][j] * (1.0 + lambda);
+    }
+
+    psGaussJordan (Alpha, Beta);
+    # endif
+
+    // apply Beta to get new Params values
+    for (int j = 0; j < params->n; j++) {
+        if ((paramMask != NULL) && (paramMask->data.U8[j])) continue;
+        // Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
+        // continue;
+        // compare Beta to beta limits
+        if (beta_lim != NULL) {
+            if (fabs(Beta->data.F64[j]) > fabs(beta_lim->data.F32[j])) {
+                Beta->data.F64[j] = (Beta->data.F64[j] > 0) ? fabs(beta_lim->data.F32[j]) : -fabs(beta_lim->data.F32[j]);
+            }
+        }
+        Params->data.F32[j] = params->data.F32[j] - Beta->data.F64[j];
+        // compare new params to param limits
+        if (params_max != NULL) {
+            Params->data.F32[j] = PS_MIN (Params->data.F32[j], params_max->data.F32[j]);
+        }
+        if (params_min != NULL) {
+            Params->data.F32[j] = PS_MAX (Params->data.F32[j], params_min->data.F32[j]);
+        }
+    }
+
+    # if (USE_LU_DECOMP)
+    psFree (A);
+    psFree (LUm);
+    psFree (LUv);
+    # endif
+
+    return true;
+}
+

trunk/psphot/src/apply_psf_model.c

@@ -22 +22 @@
     float shapeNsigma = psMetadataLookupF32 (&status, config, "PSF_SHAPE_NSIGMA");
     float OUTER       = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
+    float FIT_MIN_SN  = psMetadataLookupF32 (&status, config, "FIT_MIN_SN");
+    float FIT_MAX_CHI = psMetadataLookupF32 (&status, config, "FIT_MAX_CHI");
 
     // set the object surface-brightness limit for fitted pixels
@@ -34 +36 @@
 
         // skip non-astronomical objects (very likely defects)
-        if (source->type == PS_SOURCE_DEFECT) continue;
+        if (source->type == PS_SOURCE_DEFECT) continue; // XX should I try these anyway?
         if (source->type == PS_SOURCE_SATURATED) continue;
 
+        source->modelPSF = NULL;
+
         // use the source moments, etc to guess basic model parameters
-        psModel  *model  = pmSourceModelGuess (source, psf->type); 
+        psModel *model  = pmSourceModelGuess (source, psf->type); 
 
         // set PSF parameters for this model
@@ -44 +48 @@
         x = model->params->data.F32[2];
         y = model->params->data.F32[3];
-        // XXX I need to check if the model center has moved too much relative to the peak
 
         // set the fit radius based on the object flux limit and the model
@@ -50 +53 @@
         model->radius = modelRadius (model->params, FLUX_LIMIT) + FIT_PADDING;
         if (isnan(model->radius)) {
-          fprintf (stderr, "error in radius\n");
-          continue;
+          psAbort ("apply_psf_model", "error in radius");
         }
         
@@ -65 +67 @@
         psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f);
         if (!status || (model->params->data.F32[1] < 0)) {
-          // if the fit fails, we need to change the classification
           psLogMsg ("psphot", 3, "PSF fit failed for %f, %f (%d iterations)\n", x, y, model->nIter);
-          source->type = PS_SOURCE_OTHER;  // better choice?
+          source->type = PS_SOURCE_FAIL_FIT_PSF;  // better choice?
           continue;
         }
@@ -75 +76 @@
         Nfit ++;
 
-        mark_psf_source (source, shapeNsigma, SATURATE);
+        mark_psf_source (source, shapeNsigma, FIT_MIN_SN, FIT_MAX_CHI, SATURATE);
         if (subtract_psf_source (source)) {
           Nsub ++;

trunk/psphot/src/choose_psf_model.c

@@ -62 +62 @@
         if (test->mask->data.U8[i]) {
           source->type = PS_SOURCE_OTHER;
-          // XXX is this the right type to go to?
+          source->modelPSF = NULL;
         } else {
           source->modelPSF = test->modelPSF->data[i];

trunk/psphot/src/fit_galaxies.c

@@ -1 +1 @@
 # include "psphot.h"
-
-// fit selected galaxy model (GAUSS) to all bright objects of type GALAXY
 
 bool fit_galaxies (psImageData *imdata, psMetadata *config, psArray *sources, psStats *skyStats) 
 { 
-    bool  status;
+    bool  status, goodfit;
     float x;
     float y;
@@ -13 +11 @@
     int   Niter = 0;
 
-    float MOMENT_R = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
-    // float RADIUS   = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
-    float snFaint  = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
-    float OUTER    = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
-    float FIT_NSIGMA  = psMetadataLookupF32 (&status, config, "FIT_NSIGMA");
-    float FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING");
+    float  MOMENT_R    = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
+    float  snFaint     = psMetadataLookupF32 (&status, config, "FAINT_SN_LIM");
+    float  OUTER       = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
+    float  FIT_NSIGMA  = psMetadataLookupF32 (&status, config, "FIT_NSIGMA");
+    float  FIT_PADDING = psMetadataLookupF32 (&status, config, "FIT_PADDING");
+    char  *modelName   = psMetadataLookupPtr (&status, config, "GAL_MODEL");
 
     float FLUX_LIMIT  = FIT_NSIGMA * skyStats->sampleStdev;
 
-    psModelType   modelType   = psModelSetType ("PS_MODEL_SGAUSS");
+    psModelType   modelType   = psModelSetType (modelName);
     psModelRadius modelRadius = psModelRadius_GetFunction (modelType);
-
-    psTraceSetLevel (".psModules.pmSourceMoments", 5);
 
     psTimerStart ("psphot");
     for (int i = 0; i < sources->n; i++) {
         psSource *source = sources->data[i];
-        if (source->type != PS_SOURCE_GALAXY) continue;
-        if (source->moments->SN < snFaint) continue;
+
+        // sources which should not be fitted
+        // skip all valid stars
+        if (source->type == PS_SOURCE_PSFSTAR) continue;
+        if (source->type == PS_SOURCE_SATSTAR) continue;
+        if (source->type == PS_SOURCE_GOODSTAR) continue;
+        // skip all likely defects
+        if (source->type == PS_SOURCE_DEFECT) continue;
+        if (source->type == PS_SOURCE_SATURATED) continue;
+        // 
+        if (source->type == PS_SOURCE_FAINTSTAR) continue;
+        if (source->type == PS_SOURCE_POOR_FIT_PSF) continue;
+
+        // XXX when do we pick these up again?
+        if (source->moments->SN < snFaint) {
+          source->type = PS_SOURCE_FAINT_GALAXY;  // better choice?
+          continue;
+        }
 
         // recalculate the source moments using the galaxy radius (larger)
         status = pmSourceMoments_EAM (source, MOMENT_R);
         if (!status) {
-          fprintf (stderr, "invalid moments, skipping\n");
+          source->type = PS_SOURCE_DROP_GALAXY;  // better choice?
           continue;
         }
@@ -46 +58 @@
         y = model->params->data.F32[3];
 
-        // need a better model guess and a better radius choice
-        // when radius is not fixed, we will need to check if new radius fits on image
-
         // set the fit radius based on the object flux limit and the model
         // FLUX_LIMIT should be set based on local sky model (not global median)
-        // model->radius = 25.0;
         model->radius = modelRadius (model->params, FLUX_LIMIT) + FIT_PADDING;
-        if (isnan(model->radius)) {
-          fprintf (stderr, "error in radius\n");
-          continue;
-        }
+        if (isnan(model->radius)) psAbort ("fit_galaxies", "error in radius");
+
         if (model->radius > OUTER) {
           // allocate image, noise, mask arrays for each peak (square of radius OUTER)
@@ -64 +70 @@
         // fit as FLT, not PSF (skip poor fits)
         psImageKeepCircle (source->mask, x, y, model->radius, OR, 0x80);
-        status = pmSourceFitModel (source, model, false);
+        status = pmSourceFitModel_EAM (source, model, false);
         psImageKeepCircle (source->mask, x, y, model->radius, AND, 0x7f);
-        if (!status || (model->params->data.F32[1] < 0)) {
+        if (!status) {
           // if the fit fails, we need to change the classification
           psLogMsg ("psphot", 3, "GAL fit failed for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
-          source->type = PS_SOURCE_OTHER;  // better choice?
+          source->type = PS_SOURCE_FAIL_FIT_GAL;  // better choice?
+          source->modelFLT = model;
           Nfail ++;
           continue;
         }
 
+        goodfit = pmModelFitStatus (model);
+        if (!goodfit) {
+          // if the fit fails, we need to change the classification
+          psLogMsg ("psphot", 3, "GAL fit poor for %f, %f (%d iterations, %f radius)\n", x, y, model->nIter, model->radius);
+          source->type = PS_SOURCE_POOR_FIT_GAL;  // better choice?
+          source->modelFLT = model;
+          Nfail ++;
+          continue;
+        }
+
+        source->type = PS_SOURCE_GALAXY;
         source->modelFLT = model;
         Niter += model[0].nIter;

trunk/psphot/src/mark_psf_source.c

@@ -15 +15 @@
 // PS_SOURCE_BRIGHTSTAR 
 # define MIN_DS 0.01
-bool mark_psf_source (psSource *source, float shapeNsigma, float SATURATE)
+bool mark_psf_source (psSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE)
 { 
+    int keep;
     float dSX, dSY, SX, SY, SN;
-    float nSx, nSy;
+    float nSx, nSy, Chi;
 
     if (source->modelPSF == NULL) return (false);
@@ -26 +27 @@
     if (source->modelPSF->params->data.F32[1] >= SATURATE) {
         if (source->type == PS_SOURCE_PSFSTAR) {
-            psLogMsg ("psphot", 3, "PSFSTAR marked saturated\n");
+            psLogMsg ("psphot", 3, "PSFSTAR marked SATSTAR\n");
         }
         source->type = PS_SOURCE_SATSTAR;
@@ -32 +33 @@
     } 
     if (source->type == PS_SOURCE_SATSTAR) {
-        psLogMsg ("psphot", 4, "SATSTAR marked bright (fitted peak below saturation)\n");
-        source->type = PS_SOURCE_BRIGHTSTAR;
+        psLogMsg ("psphot", 4, "SATSTAR marked GOODSTAR (fitted peak below saturation)\n");
+        source->type = PS_SOURCE_GOODSTAR;
     }
 
@@ -41 +42 @@
     dSX = source->modelPSF->dparams->data.F32[4];
     dSY = source->modelPSF->dparams->data.F32[5];
+    Chi = source->modelPSF->chisq / source->modelPSF->nDOF;
 
     nSx = dSX / hypot (MIN_DS, 1 / (SX * SN));
@@ -48 +50 @@
     // dsx_o = hypot (1/(SX*SN), MIN_DSX)
 
-    // assign PS_SOURCE_BRIGHTSTAR to bright objects within PSF region of dparams[]
-    if ((fabs(nSx) < shapeNsigma) && (fabs(nSy) < shapeNsigma)) {
+    // assign PS_SOURCE_GOODSTAR to bright objects within PSF region of dparams[]
+    keep = TRUE;
+    keep &= (fabs(nSx) < shapeNsigma);
+    keep &= (fabs(nSy) < shapeNsigma);
+    keep &= (SN > minSN);
+    keep &= (Chi < maxChi);
+    if (keep) {
         if (source->type == PS_SOURCE_PSFSTAR) return (true);
-        source->type = PS_SOURCE_BRIGHTSTAR;
+        source->type = PS_SOURCE_GOODSTAR;
         return (true);
     }
     
     if (source->type == PS_SOURCE_PSFSTAR) {
-        psLogMsg ("psphot", 3, "PSFSTAR demoted based on dSx, dSy\n");
+        psLogMsg ("psphot", 3, "PSFSTAR demoted based on fit quality\n");
     }
 
+    // object appears to be small, suspected defect
+    if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) {
+        source->type = PS_SOURCE_DEFECT;
+        return (false);
+    }
+
+    // object appears to be large, suspected galaxy
     if ((nSx >= shapeNsigma) || (nSy >= shapeNsigma)) {
         source->type = PS_SOURCE_GALAXY;
         return (false);
     }
-    // replace DEFECT with COSMIC?
-    if ((nSx <= -shapeNsigma) || (nSy <= -shapeNsigma)) {
-        source->type = PS_SOURCE_DEFECT;
+
+    // object appears to be extremely faint: what is this?
+    if (SN < minSN) {
+        source->type = PS_SOURCE_FAINTSTAR;
         return (false);
     }
-    psTrace (".psphot.mark_psf_source", 2, "unexpected result: unmarked object\n");
+
+    // these are pooly fitted, probable stars near other stars?
+    source->type = PS_SOURCE_POOR_FIT_PSF;
     return (false);
 }       

trunk/psphot/src/mark_psf_sources.c

@@ -37 +37 @@
         nSy = dSY * SY * SN;
 
-        // assign PS_SOURCE_BRIGHTSTAR to bright objects within PSF region of dparams[]
+        // assign PS_SOURCE_GOODSTAR to bright objects within PSF region of dparams[]
         if ((fabs(nSx) < shapeNsigma) && (fabs(nSy) < shapeNsigma)) {
             if (SN > snFaint) {
-                source->type = PS_SOURCE_BRIGHTSTAR;
+                source->type = PS_SOURCE_GOODSTAR;
             } else {
                 source->type = PS_SOURCE_FAINTSTAR;

trunk/psphot/src/onesource.c

@@ -12 +12 @@
     float INNER  = psMetadataLookupF32 (&status, config, "INNER_RADIUS");
     float OUTER  = psMetadataLookupF32 (&status, config, "OUTER_RADIUS");
-    float MRAD   = psMetadataLookupF32 (&status, config, "PSF_MOMENTS_RADIUS");
+    float MRAD   = psMetadataLookupF32 (&status, config, "GAL_MOMENTS_RADIUS");
     float RADIUS = psMetadataLookupF32 (&status, config, "FIT_RADIUS");
     float SOFT   = psMetadataLookupF32 (&status, config, "SOFT_RADIUS");
@@ -21 +21 @@
 
     status = pmSourceLocalSky_EAM (source, PS_STAT_SAMPLE_MEDIAN, INNER);
-    status = pmSourceMoments (source, MRAD);
+    status = pmSourceMoments_EAM (source, MRAD);
     source->peak->counts = source->moments->Peak;
 
     psModel *model = pmSourceModelGuess (source, modelType); 
+    fprintf (stderr, "guess slope: %f\n", model->params->data.F32[7]);
 
     psImageKeepCircle (source->mask, x, y, RADIUS, OR, 0x80);
-    status = pmSourceFitModel_EAM (source, model, false, SOFT);
+    status = pmSourceFitModel_EAM (source, model, false);
 
     pmSourcePhotometry (&fitMag, &obsMag, model, source->pixels, source->mask);

trunk/psphot/src/psphot-utils.c

@@ -79 +79 @@
 {
 
+    double dP, flux;
     int i, j;
     FILE *f;
@@ -96 +97 @@
         model = (psModel  *) source->modelPSF;
         if (model == NULL) continue;
-        if (source->type == PS_SOURCE_GALAXY) continue;
-        if (source->type == PS_SOURCE_DEFECT) continue;
-        if (source->type == PS_SOURCE_SATURATED) continue;
+
+        // valid source types for this function
+        if (source->type == PS_SOURCE_SATSTAR) goto valid;
+        if (source->type == PS_SOURCE_PSFSTAR) goto valid;
+        if (source->type == PS_SOURCE_GOODSTAR) goto valid;
+        continue;
+
+    valid:
         params = model->params;
         dparams = model->dparams;
-        fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f  ", 
-                 params[0].data.F32[2], params[0].data.F32[3], 
-                 params[0].data.F32[0], -2.5*log10(params[0].data.F32[1]), (dparams[0].data.F32[1]/params[0].data.F32[1]));
+        
+        dP = 0;
+        dP += PS_SQR(dparams[0].data.F32[2]);
+        dP += PS_SQR(dparams[0].data.F32[3]);
+        dP = sqrt (dP);
+        
+        psModelFlux modelFluxFunc = psModelFlux_GetFunction (model->type);
+        flux = modelFluxFunc (params);
+
+        fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f %7.4f  ", 
+                 params[0].data.F32[2], 
+                 params[0].data.F32[3], 
+                 params[0].data.F32[0], 
+                 -2.5*log10(flux), 
+                 (dparams[0].data.F32[1]/params[0].data.F32[1]),
+                 dP);
         for (j = 0; j < model->params->n - 4; j++) {
             fprintf (f, "%9.6f ", params[0].data.F32[j+4]);
         }
+        fprintf (f, " : ");
         for (j = 0; j < model->params->n - 4; j++) {
             fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]);
         }
-        fprintf (f, ": %2d %7.4f %7.4f %7.2f (%d %d)\n", 
-                 source[0].type, log10(model[0].chisq), source[0].moments->SN, 
-                 model[0].radius, model[0].nDOF, model[0].nIter);
+        fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", 
+                 source[0].type, 
+                 log10(model[0].chisq/model[0].nDOF), 
+                 source[0].moments->SN, 
+                 model[0].radius,
+                 model[0].nDOF, 
+                 model[0].nIter);
     }
     fclose (f);
@@ -122 +146 @@
 {
 
+    double dP;
     int i, j;
     FILE *f;
@@ -139 +164 @@
         model = (psModel  *) source->modelFLT;
         if (model == NULL) continue;
+        if (source->type == PS_SOURCE_GALAXY) goto valid;
+        continue;
+
+    valid:
+        params = model->params;
+        dparams = model->dparams;
+
+        // XXX these are hardwired for SGAUSS : this should be pushed into the
+        // model functions as an abstract function 
+        dP = 0;
+        dP += PS_SQR(dparams[0].data.F32[4] / params[0].data.F32[4]);
+        dP += PS_SQR(dparams[0].data.F32[5] / params[0].data.F32[5]);
+        dP += PS_SQR(dparams[0].data.F32[7] / params[0].data.F32[7]);
+        dP = sqrt (dP);
+
+        fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f %7.4f  ", 
+                 params[0].data.F32[2], 
+                 params[0].data.F32[3], 
+                 params[0].data.F32[0], 
+                 -2.5*log10(params[0].data.F32[1]), 
+                 (dparams[0].data.F32[1]/params[0].data.F32[1]),
+                 dP);
+
+        for (j = 4; j < model->params->n; j++) {
+            fprintf (f, "%9.6f ", params[0].data.F32[j]);
+        }
+        fprintf (f, " : ");
+        for (j = 4; j < model->params->n; j++) {
+            fprintf (f, "%9.6f ", dparams[0].data.F32[j]);
+        }
+        fprintf (f, ": %2d %7.3f %7.3f %7.2f %4d %2d\n", 
+                 source[0].type, log10(model[0].chisq/model[0].nDOF), 
+                 source[0].moments->SN, 
+                 model[0].radius, 
+                 model[0].nDOF, 
+                 model[0].nIter);
+    }
+    fclose (f);
+    return true;
+}
+
+// dump the sources to an output file
+bool DumpModelNULL (psArray *sources, char *filename) 
+{
+
+    int i;
+    FILE *f;
+
+    f = fopen (filename, "w");
+    if (f == NULL) {
+        psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename);
+        return false;
+    }
+
+    psMoments *empty = pmMomentsAlloc ();
+
+    // write sources with models first 
+    for (i = 0; i < sources->n; i++) {
+        psSource *source = (psSource *) sources->data[i];
+
+        // skip these sources (in PSF or FLT)
+        if (source->type == PS_SOURCE_GALAXY) continue;
         if (source->type == PS_SOURCE_PSFSTAR) continue;
         if (source->type == PS_SOURCE_SATSTAR) continue;
-        if (source->type == PS_SOURCE_BRIGHTSTAR) continue;
-        if (source->type == PS_SOURCE_FAINTSTAR) continue;
-        if (source->type == PS_SOURCE_OTHER) continue;
-        if (source->type == PS_SOURCE_DEFECT) continue;
-        if (source->type == PS_SOURCE_SATURATED) continue;
-
-        params = model->params;
-        dparams = model->dparams;
-        fprintf (f, "%7.1f %7.1f  %5.1f %7.3f  %7.4f  ", 
-                 params[0].data.F32[2], params[0].data.F32[3], 
-                 params[0].data.F32[0], -2.5*log10(params[0].data.F32[1]), (dparams[0].data.F32[1]/params[0].data.F32[1]));
-        for (j = 0; j < model->params->n - 4; j++) {
-            fprintf (f, "%9.6f ", params[0].data.F32[j+4]);
-        }
-        for (j = 0; j < model->params->n - 4; j++) {
-            fprintf (f, "%9.6f ", dparams[0].data.F32[j+4]);
-        }
-        fprintf (f, ": %2d %7.4f %7.4f %7.2f (%d %d)\n", source[0].type, log10(model[0].chisq), source[0].moments->SN, model[0].radius, model[0].nDOF, model[0].nIter);
-    }
-    fclose (f);
-    return true;
-}
-
-// dump the sources to an output file
-bool DumpModelNULL (psArray *sources, char *filename) 
-{
-
-    int i;
-    FILE *f;
-
-    f = fopen (filename, "w");
-    if (f == NULL) {
-        psLogMsg ("DumpObjects", 3, "can't open output file for moments%s\n", filename);
-        return false;
-    }
-
-    // write sources with models first 
-    for (i = 0; i < sources->n; i++) {
-        psSource *source = (psSource *) sources->data[i];
-        if (source->modelFLT != NULL) continue;
-        if (source->modelPSF != NULL) continue;
-        if (source->moments == NULL) continue;
-        fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n", 
-                 source->peak->x, source->peak->y, source->peak->counts, 
-                 source->moments->x, source->moments->y, 
-                 source->moments->Sx, source->moments->Sy, 
-                 source->moments->Sum, source->moments->Peak, 
-                 source->moments->Sky, source->moments->SN, 
-                 source->moments->nPixels, source->type); 
+        if (source->type == PS_SOURCE_GOODSTAR) continue;
+
+        if (source->moments == NULL) {
+          fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n", 
+                   source->peak->x, source->peak->y, source->peak->counts, 
+                   empty->x, empty->y, 
+                   empty->Sx, empty->Sy, 
+                   empty->Sum, empty->Peak, 
+                   empty->Sky, empty->SN, 
+                   empty->nPixels, source->type); 
+        } else {
+          fprintf (f, "%5d %5d  %7.1f  %7.1f %7.1f  %6.3f %6.3f  %8.1f %7.1f %7.1f %7.1f  %4d %2d\n", 
+                   source->peak->x, source->peak->y, source->peak->counts, 
+                   source->moments->x, source->moments->y, 
+                   source->moments->Sx, source->moments->Sy, 
+                   source->moments->Sum, source->moments->Peak, 
+                   source->moments->Sky, source->moments->SN, 
+                   source->moments->nPixels, source->type); 
+        }
     }
     fclose (f);

trunk/psphot/src/psphot.h

@@ -49 +49 @@
 bool subtract_galaxies (psArray *sources, psMetadata *config);
 bool subtract_psf_source (psSource *source);
-bool mark_psf_source (psSource *source, float shapeNsigma, float SATURATE);
+bool mark_psf_source (psSource *source, float shapeNsigma, float minSN, float maxChi, float SATURATE);
 bool find_defects (psImage *zapmask, psArray *sources, psMetadata *config, pmPSF *psf);
 bool basic_classes (psArray *sources, psMetadata *config);
@@ -80 +80 @@
 bool pmSourceDefinePixels(psSource *mySource, const psImageData *imdata, psF32 x, psF32 y, psF32 Radius);
 bool pmSourceLocalSky_EAM (psSource *source, psStatsOptions statsOptions, psF32 Radius);
-bool pmSourceFitModel_EAM(psSource *source, psModel *model, const bool PSF, float SOFT);
+bool pmSourceFitModel_EAM(psSource *source, psModel *model, const bool PSF);
 bool pmSourceMoments_EAM(psSource *source, psF32 radius);
+bool pmModelFitStatus (psModel *model);
+psBool p_psMinLM_GuessABP_EAM (psImage  *Alpha, psVector *Beta, psVector *Params, const psImage  *alpha, const psVector *beta, const psVector *params, const psVector *paramMask, const psVector *beta_lim, const psVector *params_min, const psVector *params_max, psF64 lambda);
+psBool psMinimizeLMChi2_EAM(psMinimization *min, psImage *covar, psVector *params, const psVector *paramMask, const psArray *x, const psVector *y, const psVector *yErr, psMinimizeLMChi2Func func);
+psF64 p_psMinLM_dLinear (const psVector *Beta, const psVector *beta, psF64 lambda);
 
 // fitsource utilities

trunk/psphot/src/setup.c

@@ -75 +75 @@
     }
 
-    float XBORDER  = psMetadataLookupF32 (&status, config, "XBORDER");
-    float YBORDER  = psMetadataLookupF32 (&status, config, "YBORDER");
-    psRegion *keep = psRegionAlloc (XBORDER, -XBORDER, YBORDER, -YBORDER);
+    float XMIN  = psMetadataLookupF32 (&status, config, "XMIN");
+    float XMAX  = psMetadataLookupF32 (&status, config, "XMAX");
+    float YMIN  = psMetadataLookupF32 (&status, config, "YMIN");
+    float YMAX  = psMetadataLookupF32 (&status, config, "YMAX");
+    psRegion *keep = psRegionAlloc (XMIN, XMAX, YMIN, YMAX);
     keep           = psRegionForImage (keep, image, keep);
     psImageKeepRegion (mask, keep, OR, 0x01);

trunk/psphot/src/subtract_psf_source.c

@@ -7 +7 @@
 { 
   float sky;
+  psModel *model;
+  psImage *pixels;
 
-  // non-stellar sources are ignored
-  if (source->type == PS_SOURCE_OTHER) return (false);
-  if (source->type == PS_SOURCE_GALAXY) return (false);
-  if (source->type == PS_SOURCE_DEFECT) return (false);
-  if (source->type == PS_SOURCE_SATURATED) return (false);
+  // only subtract successful fits
+  if (source->type == PS_SOURCE_SATSTAR) goto valid;
+  if (source->type == PS_SOURCE_PSFSTAR) goto valid;
+  if (source->type == PS_SOURCE_GOODSTAR) goto valid;
+  return (false);
 
-  psModel *model = source->modelPSF;
+valid:
+  model = source->modelPSF;
   if (model == NULL) {
     psLogMsg ("psphot.subtract_psf_source", 2, "missing model for %f, %f\n", source->moments->x, source->moments->y);
@@ -20 +23 @@
   }         
 
-  psImage *pixels = source->pixels;
+  pixels = source->pixels;
 
   // subtract object, leave local sky