Changeset 15000 for trunk/psphot/src/psphotChoosePSF.c

Timestamp:

Sep 24, 2007, 11:27:58 AM (19 years ago)

Author:

eugene

Message:

update from eam_branch_20070921

File:

• trunk/psphot/src/psphotChoosePSF.c (modified) (8 diffs)

trunk/psphot/src/psphotChoosePSF.c

@@ -11 +11 @@
     sources = psArraySort (sources, pmSourceSortBySN);
 
+    // structure to store user options defining the psf
+    pmPSFOptions *options = pmPSFOptionsAlloc ();
+
+    // load user options from the recipe. no need to check existence -- they are 
     // array to store candidate PSF stars
     int NSTARS = psMetadataLookupS32 (&status, recipe, "PSF_MAX_NSTARS");
-    PS_ASSERT (status, NULL);
+    assert (status);
 
     float PSF_MIN_DS = psMetadataLookupF32 (&status, recipe, "PSF_MIN_DS");
-    PS_ASSERT (status, NULL);
+    assert (status);
 
     // supply the measured sky variance for optional constant errors (non-poissonian)
     float SKY_SIG = psMetadataLookupF32 (&status, recipe, "SKY_SIG");
-    PS_ASSERT (status, NULL);
+    assert (status);
 
     // use poissonian errors or local-sky errors
-    bool POISSON_ERRORS = psMetadataLookupBool (&status, recipe, "POISSON_ERRORS");
-    PS_ASSERT (status, NULL);
+    options->poissonErrorsPhotLMM = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PHOT.LMM");
+    assert (status);
 
     // use poissonian errors or local-sky errors
-    bool PSF_PARAM_WEIGHTS = psMetadataLookupBool (&status, recipe, "PSF_PARAM_WEIGHTS");
-    PS_ASSERT (status, NULL);
+    options->poissonErrorsPhotLin = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PHOT.LIN");
+    assert (status);
+
+    // use poissonian errors or local-sky errors
+    options->poissonErrorsParams = psMetadataLookupBool (&status, recipe, "POISSON.ERRORS.PARAMS");
+    assert (status);
 
     // how to model the PSF variations across the field
-    psMetadata *md = psMetadataLookupMetadata (&status, recipe, "PSF.TREND.MASK");
-    PS_ASSERT (status, NULL);
+    options->psfTrendMode = pmTrend2DModeFromString (psMetadataLookupStr (&status, recipe, "PSF.TREND.MODE"));
+    assert (status);
+    
+    options->psfTrendNx = psMetadataLookupS32 (&status, recipe, "PSF.TREND.NX");
+    assert (status);
+
+    options->psfTrendNy = psMetadataLookupS32 (&status, recipe, "PSF.TREND.NY");
+    assert (status);
 
     // get the fixed PSF fit radius
     // XXX check that PSF_FIT_RADIUS < SKY_OUTER_RADIUS
-    float RADIUS = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
-    PS_ASSERT (status, NULL);
-
-    pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), POISSON_ERRORS);
-
-    psPolynomial2D *psfTrendMask = psPolynomial2DfromMetadata (md);
-    if (!psfTrendMask) {
-        psError(PSPHOT_ERR_PSF, true, "Unable to construct polynomial from PSF.TREND.MASK in the recipe");
-        return NULL;
-    }
+    options->radius = psMetadataLookupF32 (&status, recipe, "PSF_FIT_RADIUS");
+    assert (status);
+
+    options->stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+
+    // dimensions of the field for which the PSF is defined
+    options->psfFieldNx = readout->image->numCols;
+    options->psfFieldNy = readout->image->numRows;
+    options->psfFieldXo = readout->image->col0;
+    options->psfFieldYo = readout->image->row0;
+
+    pmSourceFitModelInit (15, 0.01, PS_SQR(SKY_SIG), options->poissonErrorsPhotLMM);
 
     psArray *stars = psArrayAllocEmpty (sources->n);
@@ -66 +82 @@
         psLogMsg ("psphot.choosePSF", PS_LOG_WARN, "Failed to find any PSF candidates");
         psFree (stars);
-        psFree (psfTrendMask);
         return NULL;
     }
@@ -91 +106 @@
         psMetadataItem *item = psListGetAndIncrement (iter);
         char *modelName = item->data.V;
-        models->data[i] = pmPSFtryModel (stars, modelName, RADIUS, POISSON_ERRORS, psfTrendMask, PSF_PARAM_WEIGHTS, maskVal, mark);
+        models->data[i] = pmPSFtryModel (stars, modelName, options, maskVal, mark);
     }
 
@@ -97 +112 @@
     psFree (list);
     psFree (stars);
-    psFree (psfTrendMask);
 
     // select the best of the models
@@ -128 +142 @@
     if (psTraceGetLevel("psphot") >= 5) {
         for (int i = PM_PAR_SXX; i < try->psf->params->n; i++) {
-            psPolynomial2D *poly = try->psf->params->data[i];
-            for (int nx = 0; nx <= poly->nX; nx++) {
-                for (int ny = 0; ny <= poly->nY; ny++) {
-                    if (poly->mask[nx][ny]) continue;
-                    fprintf (stderr, "%g x^%d y^%d\n", poly->coeff[nx][ny], nx, ny);
-                }
-            }
+            // XXX re-write the output or some other status info
         }
     }
@@ -236 +244 @@
             // use pmModelSub because modelFlux has not been generated
             assert (source->maskObj);
-            psImageKeepCircle (source->maskObj, x, y, RADIUS, "OR", PM_MASK_MARK);
+            psImageKeepCircle (source->maskObj, x, y, options->radius, "OR", PM_MASK_MARK);
             pmModelSub (source->pixels, source->maskObj, source->modelPSF, PM_MODEL_OP_FULL, maskVal);
-            psImageKeepCircle (source->maskObj, x, y, RADIUS, "AND", PS_NOT_U8(PM_MASK_MARK));
+            psImageKeepCircle (source->maskObj, x, y, options->radius, "AND", PS_NOT_U8(PM_MASK_MARK));
         }
 
@@ -283 +291 @@
         pmSourcesWritePSFs (try->sources, "psfstars.dat");
         pmSourcesWriteEXTs (try->sources, "extstars.dat", false);
-        psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
-        psMetadataConfigWrite (psfData, "psfmodel.dat");
-        psFree (psfData);
+        // XXX need alternative output function
+        // psMetadata *psfData = pmPSFtoMetadata (NULL, try->psf);
+        // psMetadataConfigWrite (psfData, "psfmodel.dat");
         psLogMsg ("psphot.choosePSF", PS_LOG_INFO, "wrote out psf-subtracted image, psf data, exiting\n");
         exit (0);
@@ -306 +314 @@
     psLogMsg ("psphot.pspsf", PS_LOG_INFO, "psf model %s, ApResid: %f +/- %f\n", modelName, psf->ApResid, psf->dApResid);
 
+    psFree (options);
     return (psf);
 }