Changeset 33953 for branches/eam_branches/ipp-20120405/psphot/src/psphotKronIterate.c

Timestamp:

May 29, 2012, 5:16:09 PM (14 years ago)

Author:

eugene

Message:

changes from trunk

Location:

branches/eam_branches/ipp-20120405/psphot

Files:

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

branches/eam_branches/ipp-20120405/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
@@ -87 +84 @@
     int KRON_ITERATIONS = psMetadataLookupS32 (&status, recipe, "KRON_ITERATIONS");
     if (!status) {
-        KRON_ITERATIONS = 2;
+        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;
     }
 
@@ -119 +138 @@
         // 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"));
+
     }
 
@@ -146 +198 @@
             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);
+            psArrayAdd(job->args, 1, smoothedImage);
 
 // set this to 0 to run without threading
-# if (0)
+# if (1)
             if (!psThreadJobAddPending(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
@@ -156 +212 @@
             if (!psphotKronIterate_Threaded(job)) {
                 psError(PS_ERR_UNKNOWN, false, "Unable to guess model.");
+                // psFree(AnalysisRegion);
                 return false;
             }
@@ -179 +236 @@
     psFree (cellGroups);
     psFree (kronWindow);
+    psFree (smoothedImage);
 
     psLogMsg ("psphot.kron", PS_LOG_WARN, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
@@ -194 +252 @@
     float MIN_KRON_RADIUS           = PS_SCALAR_VALUE(job->args->data[6],F32);
     int KRON_ITERATIONS             = PS_SCALAR_VALUE(job->args->data[7],S32);
-
-    // XXX TEST : set iteration to 1
+    bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[8],S32);
+    bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[9],S32);
+    float KRON_SMOOTH_SIGMA         = PS_SCALAR_VALUE(job->args->data[10],F32);
+    psImage *smoothedImage          = job->args->data[11];
+
+    // psImage *smoothedPixels = NULL;
     for (int j = 0; j < KRON_ITERATIONS; j++) {
         for (int i = 0; i < sources->n; i++) {
@@ -202 +264 @@
             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) {
+                    // psFree(source->tmpPtr);
+                    // smoothedPixels = psImageSubset(smoothedImage, source->region);
+                    pmSourceSmoothOp(source, PM_MODEL_OP_FUNC, smoothedPixels, KRON_SMOOTH_SIGMA, true, maskVal, 0, 0);
+                    // source->tmpPtr = smoothedPixels;
+                }
                 reSubtract = true;
             }
@@ -223 +295 @@
             // XXX float windowRadius = PS_MIN(PS_MAX(RADIUS, 4.0*source->moments->Mrf), maxWindow);
 
-            // Sextractor apparently takes something like the skyRadius window, measures the moments inside that window,
-            // then uses an elliptical outline based on the 2nd moments, but 6x the size
-
             // 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, windowRadius;
-            
-            if (j == 0) {
-                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
-                windowRadius = PS_MAX(RADIUS, maxWindow);
-            } else {
-                maxWindow = isfinite(source->moments->Mrf) ? 6.0*source->moments->Mrf : RADIUS;
-                windowRadius = PS_MAX(RADIUS, maxWindow);
-            }
-
-#ifdef notdef
-            fprintf(stderr, "Redefining pixels for source: %d %4d %4d new radius: %f\n", 
-                                        i, source->peak->x, source->peak->y, windowRadius+2);
-#endif
+            float maxWindow;
+            if (j == 0) {
+                maxWindow = isfinite(source->skyRadius) ? source->skyRadius : RADIUS;
+            } else {
+                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.
+            // Not that it performs that function 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);
+                    if (j + 1 == KRON_ITERATIONS) {
+                        // We're done with the smoothedPixels
+                        psFree(source->tmpPtr);
+                        source->tmpPtr = NULL;
+                        smoothedPixels = NULL;
+                    }
+                }
             }
         }
     }
+    // psFree(smoothedPixels);
     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);
@@ -273 +358 @@
 
     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
@@ -286 +371 @@
     // Xn  = SUM (x - xc)^n * (z - sky)
 
+
     psF32 RF = 0.0;
     psF32 RS = 0.0;
@@ -304 +390 @@
     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;
@@ -377 +468 @@
     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++) {
@@ -415 +511 @@
 }
 
-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;
@@ -422 +518 @@
     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 
@@ -433 +537 @@
     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);
@@ -454 +568 @@
     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++) {
@@ -467 +583 @@
             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;
 }