Changeset 35559

Timestamp:

May 9, 2013, 12:15:43 PM (13 years ago)

Author:

eugene

Message:

all EXT_FIT_ITER and EXT_FIT_MIN_TOL to be set independent of the PSF versions; make psphotModelTest more useful for EXT; include mask and mark in modelGuess functions (needed for trail angle guess); allow psphot to work even if variance is missing

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/psphotExtendedSourceFits.c (modified) (6 diffs)
• src/psphotFakeSources.c (modified) (1 diff)
• src/psphotMakeResiduals.c (modified) (1 diff)
• src/psphotModelTestArguments.c (modified) (2 diffs)
• src/psphotModelTestReadout.c (modified) (10 diffs)
• src/psphotSetThreads.c (modified) (1 diff)
• src/psphotSourceFits.c (modified) (4 diffs)
• src/psphotStackImageLoop.c (modified) (1 prop)

trunk/psphot/src/models/pmModel_STRAIL.c

@@ -475 +475 @@
 
 //fixed I think...no good way of guessing as far as I can tell
-bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
+bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
 {
     pmMoments *Smoments = source->moments;

trunk/psphot/src/models/pmModel_TEST1.c

@@ -123 +123 @@
 
 // make an initial guess for parameters
-bool PM_MODEL_GUESS (pmModel *model, pmSource *source)
+bool PM_MODEL_GUESS (pmModel *model, pmSource *source, psImageMaskType maskVal, psImageMaskType markVal)
 {
     pmMoments *moments = source->moments;

trunk/psphot/src/psphotExtendedSourceFits.c

@@ -90 +90 @@
     assert (markVal);
 
+    // source fitting parameters for extended source fits
+    int fitIter = psMetadataLookupS32(&status, recipe, "EXT_FIT_ITER"); // Max number of fit iterations
+    assert (status && fitIter > 0);
+
+    float fitMinTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MIN_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+        fitMinTol = psMetadataLookupF32 (&status, recipe, "PSF_FIT_TOL"); // Fit tolerance
+        if (!status || !isfinite(fitMinTol) || fitMinTol <= 0) {
+            psAbort("PSF_FIT_MIN_TOL (and PSF_FIT_TOL) not defined or positive");
+        }
+    }
+
+    float fitMaxTol = psMetadataLookupF32 (&status, recipe, "EXT_FIT_MAX_TOL"); // Fit tolerance
+    if (!status || !isfinite(fitMaxTol) || fitMaxTol <= 0) {
+        fitMaxTol = 1.0;
+    }
+
     // maskVal is used to test for rejected pixels, and must include markVal
     maskVal |= markVal;
@@ -181 +198 @@
 
             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);
@@ -293 +314 @@
     psRegion *region        = job->args->data[4];
     int psfSize             = PS_SCALAR_VALUE(job->args->data[5],S32);
-    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[6],PS_TYPE_IMAGE_MASK_DATA);
-    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[7],PS_TYPE_IMAGE_MASK_DATA);
+    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
@@ -301 +325 @@
     fitOptions->saveCovariance = true;  // XXX make this a user option?
     fitOptions->covarFactor    = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
-
-    // XXX for now, use the defaults for the rest:
-    // fitOptions->nIter         = fitIter;
-    // fitOptions->tol           = fitTol;
-    // fitOptions->poissonErrors = poisson;
-    // fitOptions->weight        = PS_SQR(skySig);
+    fitOptions->nIter          = fitIter;
+    fitOptions->minTol         = fitMinTol;
+    fitOptions->maxTol         = fitMaxTol;
 
     // choose the sources of interest
@@ -335 +356 @@
         // set the fit radius based on the first radial moment (also sets the mask pixels)
         psphotSetRadiusMomentsExact(&fitRadius, &windowRadius, readout, source, markVal); // NOTE : 6 allocs
+
+        // XXX WATCH OUT HERE!!
+        // fitRadius = 30;
 
         // UPDATE : we have changed the moments calculation.  There is now an iteration within 
@@ -498 +522 @@
 
     // change the value of a scalar on the array (wrap this and put it in psArray.h)
-    scalar = job->args->data[8];
+    scalar = job->args->data[11];
     scalar->data.S32 = Next;
 
-    scalar = job->args->data[9];
+    scalar = job->args->data[12];
     scalar->data.S32 = Nconvolve;
 
-    scalar = job->args->data[10];
+    scalar = job->args->data[13];
     scalar->data.S32 = NconvolvePass;
 
-    scalar = job->args->data[11];
+    scalar = job->args->data[14];
     scalar->data.S32 = Nplain;
 
-    scalar = job->args->data[12];
+    scalar = job->args->data[15];
     scalar->data.S32 = NplainPass;
 
-    scalar = job->args->data[13];
+    scalar = job->args->data[16];
     scalar->data.S32 = Nfaint;
 
-    scalar = job->args->data[14];
+    scalar = job->args->data[17];
     scalar->data.S32 = Nfail;
 

trunk/psphot/src/psphotFakeSources.c

@@ -24 +24 @@
 
         pmModelType modelType = pmModelClassGetType ("PS_MODEL_QGAUSS");
-        source->modelPSF = pmSourceModelGuess (source, modelType);
+        source->modelPSF = pmSourceModelGuess (source, modelType, 0, 0);
         sources->data[i] = source;
     }

trunk/psphot/src/psphotMakeResiduals.c

@@ -105 +105 @@
         if (model == NULL) continue;  // model must be defined
 
-        psImage *image  = psImageCopy (NULL, source->pixels,   PS_TYPE_F32);
-        psImage *variance = psImageCopy (NULL, source->variance,   PS_TYPE_F32);
-        psImage *mask   = psImageCopy (NULL, source->maskView, PS_TYPE_IMAGE_MASK);
+        psImage *image    = psImageCopy (NULL, source->pixels,     PS_TYPE_F32);
+        psImage *mask     = psImageCopy (NULL, source->maskView ? source->maskView : source->maskObj,  PS_TYPE_IMAGE_MASK);
+        psImage *variance = psImageCopy (NULL, source->variance ? source->variance : source->pixels,   PS_TYPE_F32);
         pmModelSub (image, mask, model, PM_MODEL_OP_FUNC, maskVal);
 

trunk/psphot/src/psphotModelTestArguments.c

@@ -151 +151 @@
     }
 
+    // chip selection is used to limit chips to be processed
+    if ((N = psArgumentGet (argc, argv, "-moments-radius"))) {
+        psArgumentRemove (N, &argc, argv);
+        psMetadataAddBool (config->arguments, PS_LIST_TAIL, "TEST_MOMENTS_RADIUS", PS_DATA_BOOL, "", true);
+    }
+
     // if these command-line options are supplied, load the file name lists into config->arguments
     // override any configuration-specified source for these files
@@ -160 +166 @@
     pmConfigFileSetsMD (config->arguments, &argc, argv, "SRC", "-src", "-srclist");
     pmConfigFileSetsMD (config->arguments, &argc, argv, "SRCTEXT", "-srctext", "-srctextlist");
+
+    pmConfigFileSetsMD (config->arguments, &argc, argv, "PSPHOT.PSF", "-psf",      "-psflist");
 
     if (argc == 1) {

trunk/psphot/src/psphotModelTestReadout.c

@@ -5 +5 @@
 
     bool status;
+
+    psTimerStart ("modelTest");
 
     pmModelClassSetLimits(PM_MODEL_LIMITS_LAX);
@@ -41 +43 @@
     }
 
+    // load the psf model, if suppled.  FWHM_MAJ,FWHM_MIN,etc are determined and saved on
+    // readout->analysis. NOTE: this function currently only loads from PSPHOT.PSF.LOAD
+    if (!psphotLoadPSF (config, view, filerule)) { // ??? need to supply 2 ?
+        psError (PSPHOT_ERR_UNKNOWN, false, "error loading psf model");
+        return psphotReadoutCleanup (config, view, filerule);
+    }
+
     float MIN_KRON_RADIUS = 5.0;
 
@@ -53 +62 @@
     // maskVal is used to test for rejected pixels, and must include markVal
     maskVal |= markVal;
+
+    // what fraction of the PSF is used? (radius in pixels : 2 -> 5x5 box)
+    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
+    assert (status);
 
     // find the currently selected readout
@@ -84 +97 @@
 
     // define the source of interest
-   float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
+    float xObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_X");
     float yObj     = psMetadataLookupF32 (&status, recipe, "TEST_FIT_Y");
     if (!isfinite(xObj) || !isfinite(yObj)) psAbort ("object position is not defined");
@@ -99 +112 @@
     source->type = PM_SOURCE_TYPE_EXTENDED;
 
+    bool TEST_FIT_CONVOLVED = psMetadataLookupBool (&status, recipe, "TEST_FIT_CONVOLVED");
+
     // find the local sky
     status = pmSourceLocalSky (source, PS_STAT_SAMPLE_MEDIAN, INNER, maskVal, markVal);
     if (!status) psAbort("pmSourceLocalSky error");
 
-    { 
+    bool testMomentsRadius = psMetadataLookupBool (&status, config->arguments, "TEST_MOMENTS_RADIUS");
+    if (testMomentsRadius) { 
         // XXX I want to test an iterative aperture for brighter sources
         float radius = mRADIUS;
@@ -109 +125 @@
 
             // get the source moments
-            status = pmSourceMoments (source, radius, 0.0, 0.0, MIN_KRON_RADIUS, maskVal);
+            status = pmSourceMoments (source, radius, 0.25*radius, 0.0, MIN_KRON_RADIUS, maskVal);
             if (!status) psAbort("psSourceMoments error");
 
@@ -121 +137 @@
 
     // get the source moments
-    status = pmSourceMoments (source, mRADIUS, 0.0, 0.0, MIN_KRON_RADIUS, maskVal);
+    status = pmSourceMoments (source, mRADIUS, 0.25*mRADIUS, 0.0, MIN_KRON_RADIUS, maskVal);
     if (!status) psAbort("psSourceMoments error");
 
@@ -135 +151 @@
     fprintf (stderr, "axes: %f @ (%f, %f)\n", axes.theta*180/M_PI, axes.major, axes.minor);
 
+    pmPSF *psf = psMetadataLookupPtr (&status, readout->analysis, "PSPHOT.PSF");
+    if (psf) {
+      // set PSF parameters for this model (apply 2D shape model to coordinates Xo, Yo)
+      source->peak->rawFlux = source->moments->Peak;
+      float Io = source->moments->Peak;
+      source->modelPSF = pmModelFromPSFforXY(psf, xObj, yObj, Io);
+    }
+
     // get the initial model parameter guess
-    pmModel *model = pmSourceModelGuess (source, modelType);
+    pmModel *model = pmSourceModelGuess (source, modelType, maskVal, markVal);
     source->modelEXT = model;
 
@@ -159 +183 @@
     fprintf (stderr, "peak: %f @ (%f, %f)\n", source->moments->Sum*area, (double)source->peak->x, (double)source->peak->y);
 
-    // 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);
+    if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
+        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);
+    }
 
     fprintf (stderr, "input parameters: \n");
@@ -190 +217 @@
     fitOptions->covarFactor   = psImageCovarianceFactorForAperture(readout->covariance, 10.0); // Covariance matrix
 
-    status = pmSourceFitModel (source, model, fitOptions, maskVal);
+    if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        fitOptions->mode = PM_SOURCE_FIT_NO_INDEX;
+    }
+    if (modelType == pmModelClassGetType("PS_MODEL_TRAIL")) {
+        fitOptions->mode = PM_SOURCE_FIT_TRAIL;
+    }
+
+    if (TEST_FIT_CONVOLVED) {
+      pmPCMdata *pcm = pmPCMinit (source, fitOptions, model, maskVal, psfSize);
+      if (modelType == pmModelClassGetType("PS_MODEL_SERSIC")) {
+        psphotSersicModelClassGuessPCM (pcm, source);
+      } else {
+        pmSourceModelGuessPCM (pcm, source, maskVal, markVal);
+      }
+      pmPCMupdate(pcm, source, fitOptions, model);
+      pmSourceFitPCM (pcm, source, fitOptions, maskVal, markVal, psfSize);
+      psFree (pcm);
+    } else {
+      status = pmSourceFitModel (source, model, fitOptions, maskVal);
+    }
 
     // measure the source mags

trunk/psphot/src/psphotSetThreads.c

@@ -40 +40 @@
     psFree(task);
 
-    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 15);
+    task = psThreadTaskAlloc("PSPHOT_EXTENDED_FIT", 18);
     task->function = &psphotExtendedSourceFits_Threaded;
     psThreadTaskAdd(task);

trunk/psphot/src/psphotSourceFits.c

@@ -457 +457 @@
 
     // use the source moments, etc to guess basic model parameters
-    pmModel *model = pmSourceModelGuess (source, modelType);
+    pmModel *model = 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);
@@ -625 +625 @@
         model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[i];
 
-        if (!model->modelGuess(model, source)) {
+        if (!model->modelGuess(model, source, maskVal, markVal)) {
             model->flags |= PM_MODEL_STATUS_BADARGS;
             return false;
@@ -648 +648 @@
     model->flags = PM_MODEL_STATUS_NONE; // do not attempt to handle failures here, let the next iteration deal with it
     model->params->data.F32[PM_PAR_7] = 0.5/indexGuess[iMin];
-    model->modelGuess(model, source);
+    model->modelGuess(model, source, maskVal, markVal);
 
     return true;
@@ -677 +677 @@
         model->params->data.F32[PM_PAR_7] = indexGuess[i];
         
-        if (!model->modelGuess(model, source)) {
+        if (!model->modelGuess(model, source, maskVal, markVal)) {
             model->flags |= PM_MODEL_STATUS_BADARGS;
             return false;