Changeset 5768 for branches/rel9/psModules/src/objects/pmPSFtry.c

Timestamp:

Dec 12, 2005, 11:52:24 AM (21 years ago)

Author:

desonia

Message:

merged post-release fixes from CVS main.

File:

• branches/rel9/psModules/src/objects/pmPSFtry.c (modified) (7 diffs)

branches/rel9/psModules/src/objects/pmPSFtry.c

@@ +1 @@
+/** @file  pmPSFtry.c
+ *
+ *  XXX: need description of file purpose
+ *
+ *  @author EAM, IfA
+ *
+ *  @version $Revision: 1.1.8.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2005-12-12 21:52:22 $
+ *
+ *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
+ *
+ */
+
 # include <pslib.h>
-# include "psLibUtils.h"
 # include "pmObjects.h"
-# include "psModulesUtils.h"
 # include "pmPSF.h"
 # include "pmPSFtry.h"
@@ -41 +52 @@
     type           = pmModelSetType (modelName);
     test->psf      = pmPSFAlloc (type);
-    test->sources  = psMemCopy(sources);
+    test->sources  = psMemIncrRefCounter(sources);
     test->modelFLT = psArrayAlloc (sources->n);
     test->modelPSF = psArrayAlloc (sources->n);
@@ -79 +90 @@
     int Npsf = 0;
 
-    pmPSFtry *try
-    = pmPSFtryAlloc (sources, modelName);
+    pmPSFtry *psfTry = pmPSFtryAlloc (sources, modelName);
 
     // stage 1:  fit an independent model (freeModel) to all sources
     psTimerStart ("fit");
-    for (int i = 0; i < try
-                ->sources->n; i++) {
-
-            pmSource *source = try
-                                   ->sources->data[i];
-            pmModel  *model  = pmSourceModelGuess (source, try
-                                                       ->psf->type);
-            x = source->peak->x;
-            y = source->peak->y;
-
-            // set temporary object mask and fit object
-            // fit model as FLT, not PSF
-            psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
-            status = pmSourceFitModel (source, model, false);
-            psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
-
-            // exclude the poor fits
-            if (!status) {
-                try
-                    ->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
-                psFree (model);
-                continue;
-            }
-            try
-                ->modelFLT->data[i] = model;
-            Nflt ++;
-        }
+    for (int i = 0; i < psfTry->sources->n; i++) {
+
+        pmSource *source = psfTry->sources->data[i];
+        pmModel  *model  = pmSourceModelGuess (source, psfTry->psf->type);
+        x = source->peak->x;
+        y = source->peak->y;
+
+        // set temporary object mask and fit object
+        // fit model as FLT, not PSF
+        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
+        status = pmSourceFitModel (source, model, false);
+        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
+
+        // exclude the poor fits
+        if (!status) {
+            psfTry->mask->data.U8[i] = PSFTRY_MASK_FLT_FAIL;
+            psFree (model);
+            continue;
+        }
+        psfTry->modelFLT->data[i] = model;
+        Nflt ++;
+    }
     psLogMsg ("psphot.psftry", 4, "fit flt:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
     psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (FLT)\n", Nflt, sources->n);
 
     // make this optional?
-    // DumpModelFits (try->modelFLT, "modelsFLT.dat");
+    // DumpModelFits (psfTry->modelFLT, "modelsFLT.dat");
 
     // stage 2: construct a psf (pmPSF) from this collection of model fits
-    pmPSFFromModels (try
-                     ->psf, try
-                         ->modelFLT, try
-                             ->mask);
+    pmPSFFromModels (psfTry->psf, psfTry->modelFLT, psfTry->mask);
 
     // stage 3: refit with fixed shape parameters
     psTimerStart ("fit");
-    for (int i = 0; i < try
-                ->sources->n; i++) {
-            // masked for: bad model fit, outlier in parameters
-            if (try
-                    ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
-
-            pmSource *source = try
-                                   ->sources->data[i];
-            pmModel  *modelFLT = try
-                                     ->modelFLT->data[i];
-
-            // set shape for this model based on PSF
-            pmModel *modelPSF = pmModelFromPSF (modelFLT, try
-                                                    ->psf);
-            x = source->peak->x;
-            y = source->peak->y;
-
-            psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
-            status = pmSourceFitModel (source, modelPSF, true);
-
-            // skip poor fits
-            if (!status) {
-                try
-                    ->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL;
-                psFree (modelPSF);
-                goto next_source;
-            }
-
-            // otherwise, save the resulting model
-            try
-                ->modelPSF->data[i] = modelPSF;
-
-            // XXX : use a different aperture radius from the fit radius?
-            // XXX : use a different estimator for the local sky?
-            // XXX : pass 'source' as input?
-            if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
-                try
-                    ->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
-                goto next_source;
-            }
-
-            try
-                ->metric->data.F64[i] = obsMag - fitMag;
-            try
-                ->fitMag->data.F64[i] = fitMag;
-            Npsf ++;
+    for (int i = 0; i < psfTry->sources->n; i++) {
+        // masked for: bad model fit, outlier in parameters
+        if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
+            continue;
+
+        pmSource *source = psfTry->sources->data[i];
+        pmModel  *modelFLT = psfTry->modelFLT->data[i];
+
+        // set shape for this model based on PSF
+        pmModel *modelPSF = pmModelFromPSF (modelFLT, psfTry->psf);
+        x = source->peak->x;
+        y = source->peak->y;
+
+        psImageKeepCircle (source->mask, x, y, RADIUS, "OR", PSPHOT_MASK_MARKED);
+        status = pmSourceFitModel (source, modelPSF, true);
+
+        // skip poor fits
+        if (!status) {
+            psfTry->mask->data.U8[i] = PSFTRY_MASK_PSF_FAIL;
+            psFree (modelPSF);
+            goto next_source;
+        }
+
+        // otherwise, save the resulting model
+        psfTry->modelPSF->data[i] = modelPSF;
+
+        // XXX : use a different aperture radius from the fit radius?
+        // XXX : use a different estimator for the local sky?
+        // XXX : pass 'source' as input?
+        if (!pmSourcePhotometry (&fitMag, &obsMag, modelPSF, source->pixels, source->mask)) {
+            psfTry->mask->data.U8[i] = PSFTRY_MASK_BAD_PHOT;
+            goto next_source;
+        }
+
+        psfTry->metric->data.F64[i] = obsMag - fitMag;
+        psfTry->fitMag->data.F64[i] = fitMag;
+        Npsf ++;
 
 next_source:
-            psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
-
-        }
+        psImageKeepCircle (source->mask, x, y, RADIUS, "AND", ~PSPHOT_MASK_MARKED);
+
+    }
     psLogMsg ("psphot.psftry", 4, "fit psf:   %f sec for %d sources\n", psTimerMark ("fit"), sources->n);
     psTrace ("psphot.psftry", 3, "keeping %d of %d PSF candidates (PSF)\n", Npsf, sources->n);
 
     // make this optional
-    // DumpModelFits (try->modelPSF, "modelsPSF.dat");
+    // DumpModelFits (psfTry->modelPSF, "modelsPSF.dat");
 
     // XXX this function wants aperture radius for pmSourcePhotometry
-    pmPSFtryMetric (try
-                    , RADIUS);
+    pmPSFtryMetric (psfTry, RADIUS);
     psLogMsg ("psphot.pspsf", 3, "try model %s, ap-fit: %f +/- %f, sky bias: %f\n",
-              modelName, try
-                  ->psf->ApResid, try
-                      ->psf->dApResid, try
-                          ->psf->skyBias);
-
-    return (try
-           );
+              modelName,
+              psfTry->psf->ApResid,
+              psfTry->psf->dApResid,
+              psfTry->psf->skyBias);
+
+    return (psfTry);
 }
 
 
-bool pmPSFtryMetric (pmPSFtry *try
-                     , float RADIUS)
+bool pmPSFtryMetric (pmPSFtry *psfTry, float RADIUS)
 {
 
@@ -203 +193 @@
     int   nKeep, nSkip;
 
-    // the measured (aperture - fit) magnitudes (dA == try->metric)
+    // the measured (aperture - fit) magnitudes (dA == psfTry->metric)
     //   depend on both the true ap-fit (dAo) and the bias in the sky measurement:
     //     dA = dAo + dsky/flux
@@ -213 +203 @@
 
     // rflux = ten(0.4*fitMag);
-    psVector *rflux = psVectorAlloc (try
-                                     ->sources->n, PS_TYPE_F64);
-    for (int i = 0; i < try
-                ->sources->n; i++) {
-            if (try
-                    ->mask->data.U8[i] & PSFTRY_MASK_ALL) continue;
-            rflux->data.F64[i] = pow(10.0, 0.4*try
-                                     ->fitMag->data.F64[i]);
-        }
+    psVector *rflux = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
+    for (int i = 0; i < psfTry->sources->n; i++) {
+        if (psfTry->mask->data.U8[i] & PSFTRY_MASK_ALL)
+            continue;
+        rflux->data.F64[i] = pow(10.0, 0.4*psfTry->fitMag->data.F64[i]);
+    }
 
     // find min and max of (1/flux):
     psStats *stats = psStatsAlloc (PS_STAT_MIN | PS_STAT_MAX);
-    psVectorStats (stats, rflux, NULL, try
-                       ->mask, PSFTRY_MASK_ALL);
+    psVectorStats (stats, rflux, NULL, psfTry->mask, PSFTRY_MASK_ALL);
 
     // build binned versions of rflux, metric
@@ -240 +226 @@
     for (int i = 0; i < daBin->n; i++) {
 
-        psVector *tmp = psVectorAlloc (try
-                                       ->sources->n, PS_TYPE_F64);
+        psVector *tmp = psVectorAlloc (psfTry->sources->n, PS_TYPE_F64);
         tmp->n = 0;
 
         // accumulate data within bin range
-        for (int j = 0; j < try
-                    ->sources->n; j++) {
-                // masked for: bad model fit, outlier in parameters
-                if (try
-                        ->mask->data.U8[j] & PSFTRY_MASK_ALL) continue;
-
-                // skip points with extreme dA values
-                if (fabs(try
-                         ->metric->data.F64[j]) > 0.5) continue;
-
-                // skip points outside of this bin
-                if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin)
-                    continue;
-                if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin)
-                    continue;
-
-                tmp->data.F64[tmp->n] = try
-                                            ->metric->data.F64[j];
-                tmp->n ++;
-            }
+        for (int j = 0; j < psfTry->sources->n; j++) {
+            // masked for: bad model fit, outlier in parameters
+            if (psfTry->mask->data.U8[j] & PSFTRY_MASK_ALL)
+                continue;
+
+            // skip points with extreme dA values
+            if (fabs(psfTry->metric->data.F64[j]) > 0.5)
+                continue;
+
+            // skip points outside of this bin
+            if (rflux->data.F64[j] < rfBin->data.F64[i] - 0.5*dBin)
+                continue;
+            if (rflux->data.F64[j] > rfBin->data.F64[i] + 0.5*dBin)
+                continue;
+
+            tmp->data.F64[tmp->n] = psfTry->metric->data.F64[j];
+            tmp->n ++;
+        }
 
         // is this a valid point?
@@ -319 +302 @@
     stats = psVectorStats (stats, daResid, NULL, maskB, 1);
 
-    try
-        ->psf->ApResid = poly->coeff[0];
-    try
-        ->psf->dApResid = stats->clippedStdev;
-    try
-        ->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
+    psfTry->psf->ApResid = poly->coeff[0];
+    psfTry->psf->dApResid = stats->clippedStdev;
+    psfTry->psf->skyBias = poly->coeff[1] / (M_PI * PS_SQR(RADIUS));
 
     psFree (rflux);