Changeset 30624 for trunk/psphot/src/psphotRadialApertures.c

Timestamp:

Feb 13, 2011, 12:33:05 PM (15 years ago)

Author:

eugene

Message:

adding some of the metadata needed by PSPS to output headers; skip models with few valid pixels (eg, only the edge is showing); only do the linear fit on pixels within the fit radius; modify psf-match auto-scaling process; enable multiple target psfs for matched-psf aperture photometry; speed up analysis of the radial apertures

Location:

trunk/psphot

Files:

• . (modified) (1 prop)
• src/psphotRadialApertures.c (modified) (8 diffs)

trunk/psphot/src/psphotRadialApertures.c

@@ -102 +102 @@
         if (source->peak->x > AnalysisRegion.x1) continue;
         if (source->peak->y > AnalysisRegion.y1) continue;
+
+        // allocate pmSourceExtendedParameters, if not already defined
+        if (!source->radialAper) {
+            source->radialAper = psArrayAlloc(1);
+        }
 
         // replace object in image
@@ -116 +121 @@
         pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, 1.5*radius);
 
-        if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax)) {
+        if (!psphotRadialApertureSource (source, recipe, skynoise, maskVal, radMax, 0)) {
             psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
         } else {
@@ -130 +135 @@
 }
 
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax) {
-
-    // allocate pmSourceExtendedParameters, if not already defined
-    if (!source->radial) {
-        source->radial = pmSourceRadialAperturesAlloc ();
+bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
+
+    // if we are a child source, save the results to the parent source radial aperture array
+    psArray *radialAperSet = source->radialAper;
+    if (source->parent) {
+        radialAperSet = source->parent->radialAper;
+    }
+    psAssert(radialAperSet, "this should be defined before calling");
+    psAssert(radialAperSet->data[entry] == NULL, "why is this already defined?");
+
+    pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc ();
+    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);
+
+    // outer-most radius for initial truncation
+    float Rmax  = aperRadii->data.F32[aperRadii->n - 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]);
+    }
+
+    // center of the apertures
+    float xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
+    float yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+
+    // 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;
+
+            // radius is just a function of (xDiff, yDiff)
+            float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+            if (r2 > Rmax2) continue;
+
+            psVectorAppend(pixRadius2, r2);
+            psVectorAppend(pixFlux, *vPix);
+            psVectorAppend(pixVar, *vWgt);
+        }
+    }
+
+    psVector *flux    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
+    psVector *fluxErr = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
+    psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
+
+    psVectorInit (flux,    0.0);
+    psVectorInit (fluxErr, 0.0);
+    psVectorInit (fill,    0.0);
+
+    float *rPix2 = pixRadius2->data.F32;
+    for (int i = 0; i < pixRadius2->n; i++, rPix2++) {
+
+        float *aRad2 = aperRadii2->data.F32;
+        for (int j = 0; (*rPix2 < *aRad2) && (j < aperRadii2->n); j++, aRad2++) {
+            flux->data.F32[j]    += pixFlux->data.F32[i];
+            fluxErr->data.F32[j] += pixVar->data.F32[i];
+            fill->data.F32[j]    += 1.0;
+        }
+    }
+
+    for (int i = 0; i < flux->n; i++) {
+        // calculate the total flux for bin 'nOut'
+        float Area = M_PI*aperRadii2->data.F32[i];
+        fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
+        fill->data.F32[i] /= Area;
+        psTrace ("psphot", 5, "radial bins: %3d  %5.1f : %8.1f +/- %7.2f : %4.2f %6.1f\n", 
+                 i, aperRadii->data.F32[i], flux->data.F32[i], fluxErr->data.F32[i], fill->data.F32[i], Area);
     }
     
-    psVector *radius  = psVectorAllocEmpty(100, PS_TYPE_F32);
+    radialAper->flux = flux;
+    radialAper->fluxErr = fluxErr;
+    radialAper->fill = fill;
+
+    psFree (aperRadii2);
+    psFree (pixRadius2);
+    psFree (pixFlux);
+    psFree (pixVar);
+
+    return true;
+}
+
+static int nCalls = 0;
+static int nPass = 0;
+static int nPix = 0;
+
+bool psphotRadialApertureSource_With_Sort (pmSource *source, psMetadata *recipe, float skynoise, psImageMaskType maskVal, const psVector *radMax, int entry) {
+
+    psAssert(source->radialAper->data[entry] == NULL, "why is this already defined?");
+
+    pmSourceRadialApertures *radialAper = pmSourceRadialAperturesAlloc ();
+    source->radialAper->data[entry] = radialAper;
+
+    psVector *pixRadius  = psVectorAllocEmpty(100, PS_TYPE_F32);
     psVector *pixFlux = psVectorAllocEmpty(100, PS_TYPE_F32);
     psVector *pixVar  = psVectorAllocEmpty(100, PS_TYPE_F32);
@@ -144 +259 @@
         for (int ix = 0; ix < source->pixels->numCols; ix++) {
 
-            // 0.5 PIX: get radius as a function of pixel coord
+            // 0.5 PIX: get pixRadius as a function of pixel coord
             float x = ix + 0.5 - source->peak->xf + source->pixels->col0;
             float y = iy + 0.5 - source->peak->yf + source->pixels->row0;
@@ -150 +265 @@
             float r = hypot(x, y);
 
-            psVectorAppend(radius, r);
+            psVectorAppend(pixRadius, r);
             psVectorAppend(pixFlux, source->pixels->data.F32[iy][ix]);
             psVectorAppend(pixVar, source->variance->data.F32[iy][ix]);
-        }
-    }
-    psphotRadialAperturesSortFlux(radius, pixFlux, pixVar);
+            nPix ++;
+            // if (nPix % 10000 == 0) {fprintf (stderr, "?");}
+        }
+    }
+    psphotRadialAperturesSortFlux(pixRadius, pixFlux, pixVar);
 
     psVector *flux    = psVectorAllocEmpty(radMax->n, PS_TYPE_F32); // surface brightness of radial bin
@@ -174 +291 @@
 
     // XXX assume (or enforce) that the bins are contiguous and non-overlapping (Rmax[i] = Rmin[i+1])
-    for (int i = 0; !done && (i < radius->n); i++) {
-        if (radius->data.F32[i] > Rmax) {
+    for (int i = 0; !done && (i < pixRadius->n); i++) {
+        if (pixRadius->data.F32[i] > Rmax) {
             // calculate the total flux for bin 'nOut'
             float Area = M_PI*PS_SQR(Rmax);
@@ -185 +302 @@
                      nOut, radMax->data.F32[nOut], flux->data.F32[nOut], fluxErr->data.F32[nOut], fill->data.F32[nOut], Area);
 
+            nPass ++;
+            // if (nPass % 1000 == 0) {fprintf (stderr, "!");}
+
             nOut ++;
             if (nOut >= radMax->n) break;
@@ -195 +315 @@
     flux->n = fluxErr->n = fill->n = nOut;
     
-    psFree(source->radial->flux);
-    psFree(source->radial->fluxErr);
-    psFree(source->radial->fill);
-
-    source->radial->flux = flux;
-    source->radial->fluxErr = fluxErr;
-    source->radial->fill = fill;
-
-    psFree (radius);
+    radialAper->flux = flux;
+    radialAper->fluxErr = fluxErr;
+    radialAper->fill = fill;
+
+    psFree (pixRadius);
     psFree (pixFlux);
     psFree (pixVar);
 
+    nCalls ++;
+    // if (nCalls % 100 == 0) {fprintf (stderr, "*");}
     return true;
 }