Changeset 34041 for branches/meh_branches/ppstack_test/psphot/src/psphotKronIterate.c

Timestamp:

Jun 19, 2012, 5:24:19 PM (14 years ago)

Author:

mhuber

Message:

merging latest r34040 trunk changes to branch

Location:

branches/meh_branches/ppstack_test

Files:

• . (modified) (1 prop)
• psphot (modified) (1 prop)
• psphot/src (modified) (1 prop)
• psphot/src/psphotKronIterate.c (modified) (21 diffs)

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

@@ -1 +1 @@
 # include "psphotInternal.h"
-# ifndef ROUND
-# define ROUND(X) ((int) ((X) + 0.5*SIGN(X)))
-# endif
-
-bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert);
-bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal);
+
+bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool oldWindow);
+bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal, bool applyWeight, psImage *smoothedPixels);
+
 
 bool psphotKronIterate (pmConfig *config, const pmFPAview *view, const char *filerule)
@@ -11 +9 @@
     bool status = true;
 
-    // return true;
-
     fprintf (stdout, "\n");
     psLogMsg ("psphot", PS_LOG_INFO, "--- psphot Kron Iterate ---");
@@ -46 +42 @@
         // psAssert (psf, "missing psf?");
 
-        if (!psphotKronIterateReadout (config, recipe, view, readout, sources, psf)) {
+        if (!psphotKronIterateReadout (config, recipe, view, filerule, readout, sources, psf, i)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed to measure magnitudes for %s entry %d", filerule, i);
             return false;
@@ -58 +54 @@
 bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
 
-bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, pmReadout *readout, psArray *sources, pmPSF *psf) {
+bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index) {
 
     bool status = false;
@@ -69 +65 @@
     psTimerStart ("psphot.kron");
 
+
     // determine the number of allowed threads
     int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
@@ -84 +81 @@
         MIN_KRON_RADIUS = 0.25*RADIUS;
     }
+
+    int KRON_ITERATIONS = psMetadataLookupS32 (&status, recipe, "KRON_ITERATIONS");
+    if (!status) {
+        KRON_ITERATIONS = 1;
+    }
+
+    bool KRON_APPLY_WEIGHT = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WEIGHT");
+    if (!status) {
+        KRON_APPLY_WEIGHT = true;
+    }
+
+    bool KRON_APPLY_WINDOW = psMetadataLookupBool (&status, recipe, "KRON_APPLY_WINDOW");
+    if (!status) {
+        KRON_APPLY_WINDOW = false;
+    }
+    bool KRON_SMOOTH = psMetadataLookupBool (&status, recipe, "KRON_SMOOTH");
+    if (!status) {
+        KRON_SMOOTH = false;
+    }
+    float KRON_SMOOTH_SIGMA = psMetadataLookupF32 (&status, recipe, "KRON_SMOOTH_SIGMA");
+    if (!status) {
+        KRON_SMOOTH_SIGMA = 1.7;
+    }
+    float KRON_SMOOTH_NSIGMA = psMetadataLookupS32 (&status, recipe, "KRON_SMOOTH_NSIGMA");
+    if (!status) {
+        KRON_SMOOTH_NSIGMA = 2;
+    }
+    /*
+     *  Parameter for calculating maximum integration radius based on source's surface
+     *  brightness
+     *  Given minimum surface brightness SBmin and a flux the maximum radius is found from
+     *
+     *  SBmin = source->flux / (pi * Rmax**2)
+     *  Rmax = sqrt (source->flux ) / sqrt (SBmin * pi)
+     *
+     *  Now what do we use for SBmin?
+     *
+     *  SBmin = ( some flux ) / (some area)
+     *  some flux ~ flux of Ns sigma PSF source
+     *  some area ~ K times area of a PSF
+     *  
+     * flux of Ns sigma source ~ Ns * SKY_STDEV * PSF_EFFECTIVE_AREA
+     * PSF_EFFECTIVE_AREA = 4 pi sigma_PSF^2
+     * (the 4 accounts for the fact that 1 sigma is not the total area, it is
+     * actually a larger region).
+     *
+     *   SBmin = Ns * SKY_STDEV * 4 * pi * sigma_PSF^2 / (K * pi * sigma_PSF^2) 
+     *         = Ns * SKY_STDEV * 4 / K
+     *
+     * We combine the two parameters Ns and K, and the constant 4 into a single recipe value
+     * KRON_SB_MIN_FACTOR with Ns = 5 and K = 2 the corresponding value 
+     * for KRON_SB_MIN_FACTOR is 5 * 4 / 2 = 10
+     */
+    float KRON_SB_MIN_FACTOR = psMetadataLookupF32 (&status, recipe, "KRON_SB_MIN_FACTOR");
+    if (!status) {
+        KRON_SB_MIN_FACTOR = 10;
+    }
+    float SKY_STDEV = psMetadataLookupF32 (&status, readout->analysis, "MSKY_DEV");
+    float KRON_SB_MIN_DIVISOR = sqrt ( M_PI * KRON_SB_MIN_FACTOR * SKY_STDEV );
 
     // bit-masks to test for good/bad pixels
@@ -114 +170 @@
         // set a window function for each source based on the moments
         // (this skips really bad sources (no peak, no moments, DEFECT)
-        psphotKronWindowSetSource (source, kronWindow, false, true);
+        psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW);
+    }
+
+    // We measure the Kron Radius on a smoothed copy of the readout image
+    psImage *smoothedImage = NULL;
+    if (KRON_SMOOTH) {
+        // Build the smoothed source image
+        // Replace the subtracted sources
+        psphotReplaceAllSourcesReadout(config, view, filerule, index, recipe, false);
+        // Copy the image and smooth
+        psTimerStart ("psphot.kron.smooth");
+        smoothedImage = psImageCopy(NULL, readout->image, PS_TYPE_F32);
+        psImageSmooth(smoothedImage, KRON_SMOOTH_SIGMA, KRON_SMOOTH_NSIGMA);
+        psLogMsg ("psphot.kron", PS_LOG_INFO, "smoothed image %f sec\n", psTimerMark ("psphot.kron.smooth"));
+
+        // remove the sources
+        psphotRemoveAllSourcesReadout( config, view, filerule, index, recipe, false );
+        // Now subtract smooth versions of the sources from the smoothed image
+        psTimerStart ("psphot.kron.smooth.sources");
+        for (int i=0; i< sources->n; i++) {
+            pmSource *source = sources->data[i];
+            // If source has been subtracted from the readout image subtract a "smoothed" version from the smoothedImage
+            if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
+                // cache copy of smoothedPixels in the source
+                // tmpPtr is for use by a single "module" and must be null otherwise
+                psAssert(source->tmpPtr == NULL, "source->tmpPtr is not null!");
+
+                psImage *smoothedPixels = psImageSubset(smoothedImage, source->region);
+                source->tmpPtr = (psPtr) smoothedPixels;
+                pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0);
+            }
+        }
+        psLogMsg ("psphot.kron", PS_LOG_INFO, "removed %ld smoothed sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth.sources"));
+
     }
 
@@ -136 +225 @@
             psArrayAdd(job->args, 1, kronWindow);
             psArrayAdd(job->args, 1, cells->data[j]); // sources
+            psArrayAdd(job->args, 1, smoothedImage);
             PS_ARRAY_ADD_SCALAR(job->args, markVal,            PS_TYPE_IMAGE_MASK);
             PS_ARRAY_ADD_SCALAR(job->args, maskVal,            PS_TYPE_IMAGE_MASK);
             PS_ARRAY_ADD_SCALAR(job->args, RADIUS,             PS_TYPE_F32);
             PS_ARRAY_ADD_SCALAR(job->args, MIN_KRON_RADIUS,    PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, KRON_ITERATIONS,    PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WEIGHT, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, (psS32) KRON_APPLY_WINDOW, PS_TYPE_S32);
+            PS_ARRAY_ADD_SCALAR(job->args, KRON_SMOOTH_SIGMA,  PS_TYPE_F32);
+            PS_ARRAY_ADD_SCALAR(job->args, KRON_SB_MIN_DIVISOR,PS_TYPE_F32);
 
 // set this to 0 to run without threading
@@ -150 +245 @@
             if (!psphotKronIterate_Threaded(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
-                psFree(AnalysisRegion);
+                // psFree(AnalysisRegion);
                 return false;
             }
@@ -174 +269 @@
     psFree (cellGroups);
     psFree (kronWindow);
+    if (KRON_SMOOTH) {
+        for (int i = 0; i < sources->n; i++) {
+            pmSource *source = sources->data[i];
+            psFree(source->tmpPtr);
+            source->tmpPtr = NULL;
+        }
+    }
+    psFree (smoothedImage);
 
     psLogMsg ("psphot.kron", PS_LOG_WARN, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
@@ -184 +287 @@
     psImage *kronWindow             = job->args->data[1];
     psArray *sources                = job->args->data[2];
-    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[3],PS_TYPE_IMAGE_MASK_DATA);
-    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
-    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[5],F32);
-    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
-
-    // XXX TEST : set iteration to 1
-    for (int j = 0; j < 1; j++) {
+    psImage *smoothedImage          = job->args->data[3];
+    psImageMaskType markVal         = PS_SCALAR_VALUE(job->args->data[4],PS_TYPE_IMAGE_MASK_DATA);
+    psImageMaskType maskVal         = PS_SCALAR_VALUE(job->args->data[5],PS_TYPE_IMAGE_MASK_DATA);
+    float RADIUS                    = PS_SCALAR_VALUE(job->args->data[6],F32);
+    float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[7],F32);
+    int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[8],S32);
+    bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[9],S32);
+    bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[10],S32);
+    float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[11],F32);
+    float KRON_SB_MIN_DIVISOR       = PS_SCALAR_VALUE(job->args->data[12],F32);
+
+    for (int j = 0; j < KRON_ITERATIONS; j++) {
         for (int i = 0; i < sources->n; i++) {
 
@@ -196 +304 @@
             if (!source->peak) continue; // XXX how can we have a peak-less source?
 
-            // allocate space for moments
+            // check status of this source's moments
             if (!source->moments) continue;
+            if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) continue;
+            if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) continue;
 
             // replace object in image
             bool reSubtract = false;
+            psImage *smoothedPixels = NULL;
             if (source->tmpFlags & PM_SOURCE_TMPF_SUBTRACTED) {
-                pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+                pmSourceAdd (source, PM_MODEL_OP_FULL, maskVal);
+                smoothedPixels = (psImage *) source->tmpPtr;
+                if (smoothedPixels) {
+                    pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, true, maskVal, 0, 0);
+                }
                 reSubtract = true;
             }
 
-            // use S/N to control max window size
-            // float kronSN = source->moments->KronFlux / source->moments->KronFluxErr;
-
-            // maxWindow -> 1.5*RADIUS for kronSN = 5.0, keeping S.B. constant (kronSN ~ flux)
-            // (kronSN / maxWindow^2) = (5.0 / (1.5 RADIUS)^2)
-            // maxWindow = 1.5 * RADIUS * sqrt(kronSN / 5.0)
-            // XXX float maxWindow = (isfinite(kronSN) && (kronSN > 5.0)) ? 1.5 * RADIUS * sqrt(kronSN / 5.0) : 1.5*RADIUS;
-
-            // iterate to the window radius
-            // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow);
-
-            // XXX TEST : use a window based on the radial profile numbers: max is skyRadius, min is RADIUS
-            // if we lack the skyRadius (eg MATCHED sources), go to the default value
-            float maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
+            // On first iteration set window radius to sky radius (if valid) on second iteration
+            // use a factor times the previous radial moment value up to a maximum value that
+            // depends on the surface brightness of the source
+            float maxWindow;
+            if (j == 0) {
+                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
+            } else {
+                if (KRON_SB_MIN_DIVISOR) {
+                    if (isfinite(source->moments->KronFlux) && (source->moments->KronFlux > 0)) {
+                        // Limit window radius based on surface brightness
+                        float Rmax = sqrt(source->moments->KronFlux) / KRON_SB_MIN_DIVISOR;
+
+                        if (source->moments->Mrf > 0) {
+                            maxWindow = PS_MIN(6.0*source->moments->Mrf, Rmax);
+                        } else {
+                            maxWindow = Rmax;
+                        }
+                    } else {
+                        maxWindow = RADIUS;
+                    }
+                } else {
+                    // old code
+                    maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS;
+                }
+            }
             float windowRadius = PS_MAX(RADIUS, maxWindow);
 
             // re-allocate image, weight, mask arrays for each peak with box big enough to fit BIG_RADIUS
-            pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
+            bool extend = pmSourceRedefinePixels (source, readout, source->peak->x, source->peak->y, windowRadius + 2);
             psAssert (source->pixels, "WTF?");
+            if (extend && smoothedPixels) {
+                psFree(source->tmpPtr);
+                smoothedPixels = psImageSubset(smoothedImage, source->region);
+                psAssert (smoothedPixels, "WTF?");
+                source->tmpPtr = (psPtr) smoothedPixels ;
+            }
+
 
             // clear the window function for this source based on the moments
-            psphotKronWindowSetSource (source, kronWindow, (j > 0), false);
-
-            // this function populates moments->Mrf,KronFlux,KronFluxErr
-            psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal);
-            psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
-
-            // set a window function for each source based on the moments
-            psphotKronWindowSetSource (source, kronWindow, true, true);
+            // Note: this function also applies cuts on the source and returns false if it
+            // does not meet the requirements for measuring the Kron Radius or Magnitude.
+            // Note: that it performs the cuts even if KRON_APPLY_WINDOW is false
+            if (psphotKronWindowSetSource (source, kronWindow, (j > 0), false, KRON_APPLY_WINDOW)) {
+
+                // this function populates moments->Mrf,KronFlux,KronFluxErr
+                psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, smoothedPixels);
+                psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
+
+                // set a window function for each source based on the moments
+                psphotKronWindowSetSource (source, kronWindow, true, true, KRON_APPLY_WINDOW);
+            }
 
             // if we subtracted it above, re-subtract the object, leave local sky
             if (reSubtract) {
-                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+                pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
+                if (smoothedPixels) {
+                    pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, false, maskVal, 0, 0);
+                }
             }
         }
     }
+
     return true;
 }
 
-bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal) {
+bool psphotKronWindowMag (pmSource *source, psImage *kronWindow, float radius, float minKronRadius, psImageMaskType maskVal,
+    bool applyWeight, psImage *smoothedPixels) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -253 +395 @@
 
     psF32 R2 = PS_SQR(radius);
-    float rsigma2 = 0.5 / PS_SQR(radius/2.0);
+    float rsigma2 =  applyWeight ? 0.5 / R2 : 0;
 
     // a note about coordinates: coordinates of objects throughout psphot refer to the primary
@@ -266 +408 @@
     // Xn  = SUM (x - xc)^n * (z - sky)
 
+
     psF32 RF = 0.0;
     psF32 RS = 0.0;
@@ -284 +427 @@
     int Ywo = source->pixels->row0;
 
-    psF32 **vPix = source->pixels->data.F32;
+    psF32 **vPix;
+    if (smoothedPixels) {
+        vPix = smoothedPixels->data.F32;
+    } else {
+        vPix = source->pixels->data.F32;
+    }
     psF32 **vWin = kronWindow->data.F32;
     psF32 **vWgt = source->variance->data.F32;
@@ -357 +505 @@
     float Var = 0.0;
     float Win = 0.0;
+
+    // set vPix to the source pixels (it may have been set to the
+    // smoothed image above)
+    vPix = source->pixels->data.F32;
+
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -395 +548 @@
 }
 
-bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert) {
+bool psphotKronWindowSetSource(pmSource *source, psImage *kronWindow, bool useKronRadius, bool insert, bool applyWindow) {
 
     if (!source) return false;
@@ -402 +555 @@
     if (source->type == PM_SOURCE_TYPE_DEFECT) return false;
     if (source->type == PM_SOURCE_TYPE_SATURATED) return false;
+    if (!(source->tmpFlags & PM_SOURCE_TMPF_MOMENTS_MEASURED)) return false;
+    if (source->mode & PM_SOURCE_MODE_MOMENTS_FAILURE) return false;
     psAssert(kronWindow, "need a window");
+
+    if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
+
+    if (!applyWindow) {
+        return true;
+    }
 
     // we have a source with moments Mx, My, Mxx, Mxy, Myy.  we just need to define a Gaussian that has 
@@ -413 +574 @@
     float Yo = source->moments->My;
 
-    float Mxx = source->moments->Mxx;
-    float Mxy = source->moments->Mxy;
-    float Myy = source->moments->Myy;
-
-    float Mmajor = 0.5*(Mxx + Myy) + 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
-    float Mminor = 0.5*(Mxx + Myy) - 0.5*sqrt(PS_SQR(Mxx - Myy) + 4.0*PS_SQR(Mxy));
-
-    // Mxx, Mxy, Myy define the elliptical shape, but Mrf defines the width 
-    float scale = PS_SQR(0.5 * source->moments->Mrf) / Mmajor;
-
-    float Sxx = scale * Mmajor * Mminor / Myy; // sigma_x^2
-    float Sxy = Mxy / (scale * Mmajor * Mminor);
-    float Syy = scale * Mmajor * Mminor / Mxx; // sigma_y^2
-
-    float Smajor = sqrt(Mmajor);
+    psEllipseMoments moments;
+    moments.x2 = source->moments->Mxx;
+    moments.y2 = source->moments->Myy;
+    moments.xy = source->moments->Mxy;
+
+    psEllipseAxes axes = psEllipseMomentsToAxes(moments, 20.0);
+    if (! (isfinite(axes.major) && isfinite(axes.minor) && isfinite(axes.theta)) ) {
+        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
+        // fprintf(dump, "invalid axes found id: %4d major: %f minor: %f theta: %f\n", source->id, axes.major, axes.minor, axes.theta);
+        return false;
+    }
+
+    // Why this factor of 0.5 ?
+    float scale = 0.5 * source->moments->Mrf / axes.major;
+    axes.major *= scale;
+    axes.minor *= scale;
+
+    psEllipseShape shape = psEllipseAxesToShape(axes);
+    if (! (isfinite(shape.sx) && isfinite(shape.sy) && isfinite(shape.sxy)) ) {
+        // Shall we log a proper warning? This happens often with matched sources (forced photometry)
+        // fprintf(dump, "invalid shape found id: %d sx: %f sy %f sxy: %f\n", source->id, shape.sx, shape.sy, shape.sxy);
+        return false;
+    }
+
+    float Smajor = axes.major;
 
     int minX = PS_MIN(PS_MAX(Xo - 5*Smajor, 0), Nx - 1);
@@ -434 +605 @@
     int maxY = PS_MIN(PS_MAX(Yo + 5*Smajor, 0), Ny - 1);
 
-    float rMxx = 0.5 / Sxx;
-    float rMyy = 0.5 / Syy;
-    
+    float rMxx = 0.5 / PS_SQR(shape.sx);
+    float rMyy = 0.5 / PS_SQR(shape.sy);
+    float Sxy = -1. * shape.sxy;    // factor of -1 is included to match the previous window function
+                                    // implementation. XXX: Is this correct?
+
     for (int iy = minY; iy < maxY; iy++) {
         for (int ix = minX; ix < maxX; ix++) {
@@ -447 +620 @@
             float f = insert ? 1.001 - exp(-z) : 1.0 / (1.001 - exp(-z));
 
-            kronWindow->data.F32[iy][ix] *= f;
+            kronWindow->data.F32[iy][ix] *= f;
         }
     }
+
     return true;
 }