Changeset 33415 for branches/meh_branches/ppstack_test/psphot/src/psphotRadialApertures.c

Timestamp:

Mar 5, 2012, 5:19:48 PM (14 years ago)

Author:

mhuber

Message:

merging latest r33407 trunk changes to meh_branches/ppstack_test

Location:

branches/meh_branches/ppstack_test

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (2 props)
• psphot/src/psphotRadialApertures.c (modified) (16 diffs)

branches/meh_branches/ppstack_test/psphot/src

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest

branches/meh_branches/ppstack_test/psphot/src/psphotRadialApertures.c

@@ -3 +3 @@
 bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar);
 
+// this function measures the radial aperture fluxes for the set of readouts.  this function
+// may be called multiple times, presumably for different versions of PSF-matched or unmatched images.  
+
 // for now, let's store the detections on the readout->analysis for each readout
-bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule)
+bool psphotRadialApertures (pmConfig *config, const pmFPAview *view, const char *filerule, int entry)
 {
     bool status = true;
+
+    fprintf (stdout, "\n");
+    psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Radial Apertures ---");
 
     // select the appropriate recipe information
@@ -20 +26 @@
     int num = psphotFileruleCount(config, filerule);
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
-        if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe)) {
+        if (i == chisqNum) continue; // skip chisq image
+
+        if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe, entry)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i);
             return false;
@@ -30 +42 @@
 }
 
+// these values are used by all threads repeatedly (and are not modified)
+static psVector *aperRadii = NULL;
+static psVector *aperRadii2 = NULL;
+static float outerRadius = NAN;
+static float SN_LIM = NAN;
+static psImageMaskType maskVal = 0;
+
 // aperture-like measurements for extended sources
 // flux in simple, circular apertures
-bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe) {
+// 'entry' tells us which of the matched-PSF images we are working on (0 == unmatched image, also non-stack psphot)
+bool psphotRadialAperturesReadout (pmConfig *config, const pmFPAview *view, const char *filerule, int index, psMetadata *recipe, int entry) {
 
     bool status;
@@ -63 +83 @@
     }
 
-    // radMax stores the upper bounds of the annuli
+    // determine the number of allowed threads
+    int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
+    if (!status) {
+        nThreads = 0;
+    }
+
+    // aperRadii stores the upper bounds of the annuli
     // XXX keep the same name here as for the petrosian / elliptical apertures?
-    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
-    psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
-    psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
-    float outerRadius = radMax->data.F32[radMax->n - 1];
+    aperRadii = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    psAssert (aperRadii, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
+    psAssert (aperRadii->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
+
+    outerRadius = aperRadii->data.F32[aperRadii->n - 1];
+
+    // save the R^2 values as well for quicker comparison
+    aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
+    for (int i = 0; i < aperRadii->n; i++) {
+        aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
+    }
 
     // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
     assert (maskVal);
 
     // S/N limit to perform full non-linear fits
-    float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
+    SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
     // source analysis is done in S/N order (brightest first)
@@ -81 +114 @@
     sources = psArraySort (sources, pmSourceSortByFlux);
 
+    // XXX make this consistent with entry 0 == unmatched
+    int nEntry = 1;
+    psVector *fwhmValues = psMetadataLookupVector(&status, readout->analysis, "STACK.PSF.FWHM.VALUES");
+    if (fwhmValues) {
+        psAssert (entry < fwhmValues->n, "inconsistent matched-PSF entry");
+        nEntry = fwhmValues->n;
+    }
+    if (entry > 0) {
+        psLogMsg ("psphot", PS_LOG_DETAIL, "Radial Apertures for matched image %s : PSF FWHM = %f pixels\n", file->name, fwhmValues->data.F32[entry]);
+    } else {
+        psLogMsg ("psphot", PS_LOG_DETAIL, "Radial Apertures for unmatched image %s\n", file->name);
+    }
+
     // option to limit analysis to a specific region
     char *region = psMetadataLookupStr (&status, recipe, "ANALYSIS_REGION");
-    psRegion AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
-    if (psRegionIsNaN (AnalysisRegion)) psAbort("analysis region mis-defined");
+    psRegion *AnalysisRegion = psRegionAlloc(0,0,0,0);
+    *AnalysisRegion = psRegionForImage (readout->image, psRegionFromString (region));
+    if (psRegionIsNaN (*AnalysisRegion)) psAbort("analysis region mis-defined");
+
+    // choose Cx, Cy (see psphotThreadTools.c for overview of the concepts)
+    int Cx = 1, Cy = 1;
+    psphotChooseCellSizes (&Cx, &Cy, readout, nThreads);
+
+    psArray *cellGroups = psphotAssignSources (Cx, Cy, sources);
+
+    for (int i = 0; i < cellGroups->n; i++) {
+
+        psArray *cells = cellGroups->data[i];
+
+        for (int j = 0; j < cells->n; j++) {
+
+            // allocate a job -- if threads are not defined, this just runs the job
+            psThreadJob *job = psThreadJobAlloc ("PSPHOT_RADIAL_APERTURES");
+
+            psArrayAdd(job->args, 1, readout);
+            psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, AnalysisRegion);
+            PS_ARRAY_ADD_SCALAR(job->args, entry,  PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, nEntry, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial
+
+            // set this to 0 to run without threading
+# if (1)            
+            if (!psThreadJobAddPending(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree(AnalysisRegion);
+                return false;
+            } 
+# else
+            if (!psphotRadialApertures_Threaded(job)) {
+                psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                psFree(AnalysisRegion);
+                return false;
+            }
+            psScalar *scalar = NULL;
+            scalar = job->args->data[5];
+            Nradial += scalar->data.S32;
+            psFree(job);
+# endif
+        }
+
+        // wait for the threads to finish and manage results
+        if (!psThreadPoolWait (false, true)) {
+            psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+            psFree(AnalysisRegion);
+            return false;
+        }
+
+        // we have only supplied one type of job, so we can assume the types here
+        psThreadJob *job = NULL;
+        while ((job = psThreadJobGetDone()) != NULL) {
+            if (job->args->n < 1) {
+                fprintf (stderr, "error with job\n");
+            } else {
+                psScalar *scalar = NULL;
+                scalar = job->args->data[5];
+                Nradial += scalar->data.S32;
+            }
+            psFree(job);
+        }
+    }
+    psFree (cellGroups);
+    psFree(AnalysisRegion);
+    psFree (aperRadii2);
+
+    psLogMsg ("psphot", PS_LOG_WARN, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
+    return true;
+}
+ 
+bool psphotRadialApertures_Threaded (psThreadJob *job) {
+
+    int Nradial = 0;
+
+    // arguments: readout, sources, models, region, psfSize, maskVal, markVal
+    pmReadout *readout      = job->args->data[0];
+    psArray *sources        = job->args->data[1];
+    psRegion *region        = job->args->data[2];
+    int entry               = PS_SCALAR_VALUE(job->args->data[3],S32); // which psf-matched image are we working on? (0 == unmatched)
+    int nEntry              = PS_SCALAR_VALUE(job->args->data[4],S32); // total number of psf-matched images + 1 unmatched
+
+    // storage for the derived pixel values (these are passed into psphotRadialApertureSource)
+    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
 
     // choose the sources of interest
@@ -90 +223 @@
 
         pmSource *source = sources->data[i];
+
+        // if we have checked the source validity on the basis of the object set, then 
+        // we either skip these tests below or we skip the source completely
+        if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_SKIP) continue;
+        if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_KEEP) goto keepSource;
 
         // skip PSF-like and non-astronomical objects
@@ -104 +242 @@
 
         // limit selection by analysis region
-        if (source->peak->x < AnalysisRegion.x0) continue;
-        if (source->peak->y < AnalysisRegion.y0) continue;
-        if (source->peak->x > AnalysisRegion.x1) continue;
-        if (source->peak->y > AnalysisRegion.y1) continue;
+        if (source->peak->x < region->x0) continue;
+        if (source->peak->y < region->y0) continue;
+        if (source->peak->x > region->x1) continue;
+        if (source->peak->y > region->y1) continue;
+
+    keepSource:
 
         // allocate pmSourceExtendedParameters, if not already defined
-        if (!source->radialAper) {
-            source->radialAper = psArrayAlloc(1);
+        // XXX check that nPSFsizes is consistent with targets
+        if (source->parent) {
+            if (!source->parent->radialAper) {
+                source->parent->radialAper = psArrayAlloc(nEntry);
+            }
+        } else {
+            if (!source->radialAper) {
+                source->radialAper = psArrayAlloc(nEntry);
+            }
         }
 
@@ -119 +266 @@
         }
 
-        // we need to change the view for the radial aperture analysis, but we want to recover exactly
-        // the original view; the following elements get destroyed by pmSourceRedefinePixels so save them:
-        psImage *oldMaskObj   = psMemIncrRefCounter(source->maskObj);
-        psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux);
-        psImage *oldPSFimage  = psMemIncrRefCounter(source->psfImage);
-        psRegion oldRegion    = source->region;
-
         Nradial ++;
 
-        // force source image to be a bit larger...
-        pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
-
-        if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, 0)) {
+        if (!psphotRadialApertureSource (source, readout, entry, pixRadius2, pixFlux, pixVar)) {
             psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
         } else {
@@ -137 +274 @@
         }
 
-        pmSourceRedefinePixelsByRegion (source, readout, oldRegion);
-        psFree(source->maskObj);   source->maskObj   = oldMaskObj;
-        psFree(source->modelFlux); source->modelFlux = oldModelFlux;
-        psFree(source->psfImage);  source->psfImage  = oldPSFimage;
-        
         // re-subtract the object, leave local sky
         pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     }
-
-    psLogMsg ("psphot", PS_LOG_INFO, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
+    psScalar *scalar = job->args->data[5];
+    scalar->data.S32 = Nradial;
+
+    psFree (pixRadius2);
+    psFree (pixFlux);
+    psFree (pixVar);
+
     return true;
 }
 
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
-
+bool psphotRadialApertureSource (pmSource *source, pmReadout *readout, int entry, psVector *pixRadius2, psVector *pixFlux, psVector *pixVar) {
+                                            
     // if we are a child source, save the results to the parent source radial aperture array
     psArray *radialAperSet = source->radialAper;
@@ -163 +300 @@
     radialAperSet->data[entry] = radialAper;
 
-    // storage for the derived pixel values
-    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
+    // find the largest aperture of interest (use only apertures with inner radii <=
+    // source->skyRadius)
+    int lastAp = aperRadii->n;
+    for (int i = 0; i < aperRadii->n; i++) {
+        if (aperRadii->data.F32[i] < source->skyRadius) continue;
+        lastAp = i + 1;
+        break;
+    }
 
     // outer-most radius for initial truncation
-    float Rmax  = aperRadii->data.F32[aperRadii->n - 1];
+    float Rmax  = aperRadii->data.F32[lastAp - 1];
     float Rmax2 = PS_SQR(Rmax);
 
-    // store the R^2 values for the apertures
-    psVector *aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
-    for (int i = 0; i < aperRadii->n; i++) {
-        aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
-    }
+    // in this function, the operatins are relative to the full image (readout->image, etc)
 
     float xCM = NAN, yCM = NAN;
     if (pmSourcePositionUseMoments(source)) {
-        xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
-        yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+        xCM = source->moments->Mx; // index coord of peak in readout
+        yCM = source->moments->My; // index coord of peak in readout
     } else {
-        xCM = source->peak->xf - 0.5 - source->pixels->col0; // coord of peak in subimage
-        yCM = source->peak->yf - 0.5 - source->pixels->row0; // coord of peak in subimage
-    }
+        xCM = source->peak->xf; // index coord of peak in readout
+        yCM = source->peak->yf; // index coord of peak in readout
+    }
+
+    int Nx = readout->image->numCols;
+    int Ny = readout->image->numRows;
+
+    pixRadius2->n = 0;
+    pixFlux->n = 0;
+    pixVar->n = 0;
 
     // one pass through the pixels to select the valid pixels and calculate R^2
-    for (int iy = 0; iy < source->pixels->numRows; iy++) {
-
-        float yDiff = iy - yCM;
-        if (fabs(yDiff) > Rmax) continue;
-
-        float *vPix = source->pixels->data.F32[iy];
-        float *vWgt = source->variance->data.F32[iy];
-        psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy];
-
-        for (int ix = 0; ix < source->pixels->numCols; ix++, vPix++, vWgt++) {
-
-            if (vMsk) {
-                if (*vMsk & maskVal) {
-                    vMsk++;
-                    continue;
-                }
-                vMsk++;
-            }
-            if (isnan(*vPix)) continue;
-
-            float xDiff = ix - xCM;
-            if (fabs(xDiff) > Rmax) continue;
+    for (int iy = -Rmax; iy < Rmax + 1; iy++) {
+
+        float yDiff = iy + 0.5 + yCM;  // y-coordinate at this offse
+        int yPix = (int) yDiff;
+
+        if (yPix < 0) continue;
+        if (yPix > Ny - 1) continue;
+        if (fabs(iy) > Rmax) continue;
+
+        float *vPix = readout->image->data.F32[yPix];
+        float *vWgt = readout->variance->data.F32[yPix];
+        psImageMaskType  *vMsk = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[yPix];
+
+        for (int ix = -Rmax; ix < Rmax + 1; ix++) {
+
+            float xDiff = ix + 0.5 + xCM;  // x-coordinate at this offse
+            int xPix = (int) xDiff;
+            
+            if (xPix < 0) continue;
+            if (xPix > Nx - 1) continue;
+            if (fabs(ix) > Rmax) continue;
+            
+            if (vMsk[xPix] & maskVal) continue;
+            if (isnan(vPix[xPix])) continue;
 
             // radius is just a function of (xDiff, yDiff)
-            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            float r2  = PS_SQR(ix) + PS_SQR(iy);
             if (r2 > Rmax2) continue;
 
             psVectorAppend(pixRadius2, r2);
-            psVectorAppend(pixFlux, *vPix);
-            psVectorAppend(pixVar, *vWgt);
+            psVectorAppend(pixFlux, vPix[xPix]);
+            psVectorAppend(pixVar, vWgt[xPix]);
         }
     }
@@ -226 +372 @@
     psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
 
+    // init the apertures of interest to 0.0, the rest go to NAN
     psVectorInit (flux,    0.0);
     psVectorInit (fluxStd, 0.0);
     psVectorInit (fluxErr, 0.0);
     psVectorInit (fill,    0.0);
+    for (int i = lastAp; i < flux->n; i++) {
+        flux->data.F32[i] = NAN;
+        fluxStd->data.F32[i] = NAN;
+        fluxErr->data.F32[i] = NAN;
+        fill->data.F32[i] = NAN;
+    }
 
     float *rPix2 = pixRadius2->data.F32;
@@ -236 +389 @@
         int j = 0;
         float *aRad2 = aperRadii2->data.F32;
-        for (; (*aRad2 < *rPix2) && (j < aperRadii2->n); j++, aRad2++);
-        for (; j < aperRadii2->n; j++, aRad2++) {
+        for (; (*aRad2 < *rPix2) && (j < lastAp); j++, aRad2++);
+
+        // XXX I can speed this up by only saving this single aperture
+        for (; j < lastAp; j++, aRad2++) {
             flux->data.F32[j]    += pixFlux->data.F32[i];
             fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]);
@@ -249 +404 @@
        2) the fractional fill factor (count of valid pixels / effective area of the aperture
        3) the error on the flux within that aperture
-     */
-
-    for (int i = 0; i < flux->n; i++) {
+    */
+
+    for (int i = 0; i < lastAp; i++) {
         // calculate the total flux for bin 'nOut'
         float Area = M_PI*aperRadii2->data.F32[i];
@@ -261 +416 @@
         // XXX report the total flux or the mask-corrected flux?
         // flux->data.F32[i]    = SBmean * Area;
-        // fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
+        // fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / Pinx;
 
         fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
@@ -271 +426 @@
     }
     
+# if (1)
     radialAper->flux = flux;
     radialAper->fluxStdev = fluxStd;
     radialAper->fluxErr = fluxErr;
     radialAper->fill = fill;
-
-    psFree (aperRadii2);
-    psFree (pixRadius2);
-    psFree (pixFlux);
-    psFree (pixVar);
+# else
+    // XXX TEST
+    psFree(flux);
+    psFree(fluxStd);
+    psFree(fluxErr);
+    psFree(fill);
+# endif
 
     return true;
 }
+
+/*** below is a test to use a sort to speed this up, not very successfully ***/
 
 static int nCalls = 0;
@@ -373 +533 @@
 // *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
 # define COMPARE_VECT(A,B) (radius->data.F32[A] < radius->data.F32[B])
-# define SWAP_VECT(TYPE,A,B) { \
-  float tmp; \
-  if (A != B) { \
-    tmp = radius->data.F32[A]; \
-    radius->data.F32[A] = radius->data.F32[B];  \
-    radius->data.F32[B] = tmp; \
-    tmp = pixFlux->data.F32[A]; \
-    pixFlux->data.F32[A] = pixFlux->data.F32[B]; \
-    pixFlux->data.F32[B] = tmp; \
-    tmp = pixVar->data.F32[A]; \
-    pixVar->data.F32[A] = pixVar->data.F32[B]; \
-    pixVar->data.F32[B] = tmp; \
-  } \
-}
+# define SWAP_VECT(TYPE,A,B) {                                  \
+        float tmp;                                              \
+        if (A != B) {                                           \
+            tmp = radius->data.F32[A];                          \
+            radius->data.F32[A] = radius->data.F32[B];          \
+            radius->data.F32[B] = tmp;                          \
+            tmp = pixFlux->data.F32[A];                         \
+            pixFlux->data.F32[A] = pixFlux->data.F32[B];        \
+            pixFlux->data.F32[B] = tmp;                         \
+            tmp = pixVar->data.F32[A];                          \
+            pixVar->data.F32[A] = pixVar->data.F32[B];          \
+            pixVar->data.F32[B] = tmp;                          \
+        }                                                       \
+    }
 
 bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar) {

@@ -22 +22 @@
 psphotMakePSF
 psphotStack
+psphotModelTest