Changeset 35768 for trunk/psModules/src/objects/pmPSFtryFitPSF.c

Timestamp:

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

Author:

eugene

Message:

deprecate KiiOpen,KiiClose (now KapaOpen,etc); major rework of psEllipse translations : use common functions pmModelAxesToParams and pmModelParamsToAxes ; use new convergence method in pmPCM_MinimizeChisq; add convergence crerition options to psMinimization; threaded versions of pmPSFtryFitEXT and pmPSFtryFitPSF

Location:

trunk/psModules

Files:

• . (modified) (1 prop)
• src/objects/pmPSFtryFitPSF.c (modified) (5 diffs)

trunk/psModules/src/objects/pmPSFtryFitPSF.c

@@ -43 +43 @@
 #include "pmSourceVisual.h"
 
+static int Npsf = 0;
+
 // stage 3: Refit with fixed shape parameters.  This function uses the LMM fitting, but could
 // be re-written to use the simultaneous linear fitting (see psphotFitSourcesLinear.c)
 bool pmPSFtryFitPSF (pmPSFtry *psfTry, pmPSFOptions *options, psImageMaskType maskVal, psImageMaskType markVal) {
 
-    bool status;
-
     psTimerStart ("psf.fit");
 
@@ -57 +57 @@
     maskVal |= markVal;
 
-    // DEBUG code: save the PSF model fit data in detail
-# ifdef DEBUG
-    char filename[64];
-    snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy);
-    FILE *f = fopen (filename, "w");
-    psAssert (f, "failed open");
-# endif
-
-    int Npsf = 0;
+    Npsf = 0;
     for (int i = 0; i < psfTry->sources->n; i++) {
 
@@ -77 +69 @@
         }
 
-        // set shape for this model based on PSF
+        // set shape for this model based on PSF
         psFree (source->modelPSF);
-        source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
-        if (source->modelPSF == NULL) {
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-        // PSF fit and aperture mags use different radii
-        source->modelPSF->fitRadius = options->fitRadius;
-        source->apRadius = options->apRadius;
-
-        // set object mask to define valid pixels for PSF model fit
-        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
-
-        // fit the PSF model to the source
-        status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal);
-
-        // skip poor fits
-        if (!status) {
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-
-        // set object mask to define valid pixels for APERTURE magnitude
-        if (options->fitRadius != options->apRadius) {
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-            psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal);
+        source->modelPSF = pmModelFromPSF (source->modelEXT, psfTry->psf);
+        if (source->modelPSF == NULL) {
+            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_MODEL;
+            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : bad PSF fit\n", i, source->peak->x, source->peak->y);
+            return false;
         }
 
-        // This function calculates the psf and aperture magnitudes
-        status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag
-        if (!status || isnan(source->apMag)) {
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-            psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
-            psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
-            continue;
-        }
-
-        // clear object mask to define valid pixels
-        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
-
-        psfTry->fitMag->data.F32[i] = source->psfMag;
-        psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
-        psfTry->metricErr->data.F32[i] = source->psfMagErr;
-
-        // XXX this did not work: modifies shape of psf too much
-        // psfTry->metric->data.F32[i] = -2.5*log10(source->moments->Sum) - source->psfMag;
-
+        // do some actual work on this source
+        psThreadJob *job = psThreadJobAlloc ("PSF_TRY_FIT_PSF");
+        psArrayAdd(job->args, 1, source);
+        psArrayAdd(job->args, 1, psfTry);
+        psArrayAdd(job->args, 1, options);
+        
+        PS_ARRAY_ADD_SCALAR(job->args, i,        PS_TYPE_S32);
+
+        PS_ARRAY_ADD_SCALAR(job->args, maskVal,  PS_TYPE_IMAGE_MASK);
+        PS_ARRAY_ADD_SCALAR(job->args, markVal,  PS_TYPE_IMAGE_MASK);
+
+# if (1)
+        if (!psThreadJobAddPending(job)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
+            return false;
+        }
+# else
+        if (!pmPSFtryFitPSF_Threaded(job)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to create psf model.");
+            return false;
+        }
+        psFree(job);
+# endif
+    }
+
+    // wait for the threads to finish and manage results
+    if (!psThreadPoolWait (false, true)) {
+        psError(PS_ERR_UNKNOWN, false, "failure to model psf");
+        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);
+    }
+    psfTry->psf->nPSFstars = Npsf;
+
+    // DEBUG code: save the PSF model fit data in detail
 # ifdef DEBUG
+
+    char filename[64];
+    snprintf (filename, 64, "psffit.%dx%d.dat", psfTry->psf->trendNx, psfTry->psf->trendNy);
+    FILE *f = fopen (filename, "w");
+    psAssert (f, "failed open");
+
+    for (int i = 0; i < psfTry->sources->n; i++) {
+
+        // skip masked sources
+        if (psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] & PSFTRY_MASK_ALL) continue;
+
+        pmSource *source = psfTry->sources->data[i];
+
         fprintf (f, "%6.1f %6.1f : %6.1f %6.1f : %8.3f %8.3f %8.3f : %f : %f %f %f : %f\n",
                  source->peak->xf, source->peak->yf, 
@@ -136 +140 @@
                  source->modelPSF->params->data.F32[PM_PAR_SXX], source->modelPSF->params->data.F32[PM_PAR_SXY], 
                  source->modelPSF->params->data.F32[PM_PAR_SYY], source->modelPSF->params->data.F32[PM_PAR_7]);
-# endif
-
-        psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
-        Npsf ++;
-    }
-    psfTry->psf->nPSFstars = Npsf;
-
-# ifdef DEBUG
+    }
     fclose (f);
 # endif
@@ -159 +156 @@
     return true;
 }
+
+bool pmPSFtryFitPSF_Threaded (psThreadJob *job) {
+
+    pmSource *source =      job->args->data[0];
+    pmPSFtry *psfTry =      job->args->data[1];
+    pmPSFOptions *options = job->args->data[2];
+
+    int i = PS_SCALAR_VALUE(job->args->data[3], S32);
+
+    psImageMaskType maskVal = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType markVal = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
+
+    // PSF fit and aperture mags use different radii
+    source->modelPSF->fitRadius = options->fitRadius;
+    source->apRadius            = options->apRadius;
+
+    // set object mask to define valid pixels for PSF model fit
+    psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->fitRadius, "OR", markVal);
+
+    // fit the PSF model to the source
+    bool status = pmSourceFitModel (source, source->modelPSF, options->fitOptions, maskVal);
+
+    // skip poor fits
+    if (!status) {
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+        psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_PSF_FAIL;
+        psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : failed PSF fit\n", i, source->peak->x, source->peak->y);
+        return true;
+    }
+
+    // set object mask to define valid pixels for APERTURE magnitude
+    if (options->fitRadius != options->apRadius) {
+        psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+        psImageKeepCircle (source->maskObj, source->peak->x, source->peak->y, options->apRadius, "OR", markVal);
+    }
+
+    // This function calculates the psf and aperture magnitudes
+    status = pmSourceMagnitudes (source, psfTry->psf, PM_SOURCE_PHOT_INTERP, maskVal, markVal, options->apRadius); // raw PSF mag, AP mag
+    psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal)); // clear the circular mask
+
+    if (!status || isnan(source->apMag)) {
+        psfTry->mask->data.PS_TYPE_VECTOR_MASK_DATA[i] = PSFTRY_MASK_BAD_PHOT;
+        psTrace ("psModules.objects", 4, "dropping %d (%d,%d) : poor photometry\n", i, source->peak->x, source->peak->y);
+        return true;
+    }
+
+    psfTry->fitMag->data.F32[i] = source->psfMag;
+    psfTry->metric->data.F32[i] = source->apMag - source->psfMag;
+    psfTry->metricErr->data.F32[i] = source->psfMagErr;
+
+    psTrace ("psModules.object", 6, "keeping source %d (%d) of %ld\n", i, Npsf, psfTry->sources->n);
+    Npsf ++;
+    return true;
+}