Changeset 33877 for trunk/psphot/src/psphotKronIterate.c

Timestamp:

May 16, 2012, 4:21:50 PM (14 years ago)

Author:

bills

Message:

Do Kron measurements on smoothed image. Version with instrumentation.

File:

• trunk/psphot/src/psphotKronIterate.c (modified) (23 diffs)

trunk/psphot/src/psphotKronIterate.c

@@ -1 +1 @@
 # include "psphotInternal.h"
 
+#define NO_USE_SE_KR   1
+#define NO_SAVE_IMAGES 1
+
 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, bool smooth);
+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)
@@ -41 +45 @@
         // 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;
@@ -53 +57 @@
 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 pass1 = true;
+// This is defined in pmSource.c
+extern psString psphot_outroot;
+
+#ifdef USE_SE_KR
+static int   nKR = 0;
+// input list indexed by list file's source->id
+static float kr_sextractor[10000];
+static int used[10000];
+
+int numMatched = 0;
+
+#define SE_KR_USED 0x800
+
+static void find_se_kr(pmSource *source) {
+    int sid = source->id;
+    if (isfinite(kr_sextractor[sid])) {
+        source->moments->Mrf = kr_sextractor[sid] / 2.5;
+        source->tmpFlags |= SE_KR_USED;
+        numMatched++;
+    }
+}
+
+static void getxy(pmSource *source, float *pX, float *pY) {
+    float x, y;
+    if (source->modelPSF) {
+        psF32 *PAR = source->modelPSF->params->data.F32;
+        x = PAR[PM_PAR_XPOS];
+        y = PAR[PM_PAR_YPOS];
+    } else if (pmSourcePositionUseMoments(source)) {
+        x = source->moments->Mx;
+        y = source->moments->My;
+    } else {
+        x = source->peak->x;
+        y = source->peak->y;
+    }
+    *pX = x;
+    *pY = y;
+}
+
+static void build_kr_list(psArray *sources) {
+    // write list of source->id X Y to a file
+    pid_t pid = getpid();
+    char filename[40];
+    sprintf(filename, "sources.%d.list", pid);
+    FILE *out = fopen(filename, "w");
+    for (int i=0; i < sources->n; i++) {
+        pmSource *source = sources->data[i];
+        float x, y;
+        getxy(source, &x, &y);
+        fprintf(out, "%6.3f %6.3f %5d\n", x, y, source->id);
+    }
+    fclose(out);
+    // run gcompare on that file
+
+    char matchedfilename[20];
+    sprintf(matchedfilename, "sources.%d.matched", pid);
+    char command[80];
+    sprintf(command, "rungcompare %s %s", filename, matchedfilename);
+    int rc = system(command);
+    if (rc) {
+        fprintf(stderr, "failed to read matched list %d %d\n", rc, rc >> 8);
+        return;
+    }
+
+    unlink(filename);
+
+    for (int i = 0; i < 10000; i++) {
+        used[i] = -1;
+        kr_sextractor[i] = NAN;
+    }
+    int sid;
+    int ipp_idet;
+    float kr;
+    // read the results to build the kr_sextractor
+    FILE *in = fopen(matchedfilename, "r");
+    while (fscanf(in, "%d %f %d", &sid, &kr, &ipp_idet) > 0) {
+        // take first match for each sid ...
+        if (!isfinite(kr_sextractor[sid])) { 
+            // .. and for each ipp_idet Since the lists are sorted in
+            // order of SN this makes it more likely to get the right match
+            if (used[ipp_idet] == -1) {
+                kr_sextractor[sid] = kr;
+                used[ipp_idet] = sid;
+                nKR++;
+            } else {
+                fprintf(stderr, "match for %d %d is already used for %d\n",
+                    sid, ipp_idet, used[ipp_idet]);
+            }
+        }
+    }
+    fclose(in);
+    unlink(matchedfilename);
+}
+
+#endif // USE_SE_KR
+
+bool psphotKronIterateReadout(pmConfig *config, psMetadata *recipe, const pmFPAview *view, const char * filerule, pmReadout *readout, psArray *sources, pmPSF *psf, int index) {
 
     bool status = false;
@@ -64 +165 @@
     psTimerStart ("psphot.kron");
 
+#ifdef USE_SE_KR
+    if (!pass1) {
+        build_kr_list(sources);
+    }
+#endif
+
     // determine the number of allowed threads
     int nThreads = psMetadataLookupS32(&status, config->arguments, "NTHREADS"); // Number of threads
@@ -97 +204 @@
     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;
     }
 
@@ -122 +237 @@
 
     // start with the currently known moments (Mxx, Myy, Mxy) and generate a window image
+#ifdef USE_SE_KR
+    int previous = numMatched;
+#endif
     for (int i = 0; i < sources->n; i++) {
 
@@ -129 +247 @@
         // (this skips really bad sources (no peak, no moments, DEFECT)
         psphotKronWindowSetSource (source, kronWindow, false, true, KRON_APPLY_WINDOW);
-    }
+#ifdef USE_SE_KR
+        if (!pass1) {
+            find_se_kr(source);
+        }
+#endif
+    }
+#ifdef USE_SE_KR
+    if (!pass1) {
+        fprintf(stdout, "Matched %d sources previous %d\n", numMatched, previous);
+    }
+#endif
+
+    // 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_WARN, "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_WARN, " removed %ld smoothed sources %f sec\n", sources->n, psTimerMark ("psphot.kron.smooth.sources"));
+
+    }
+
+#ifdef SAVE_IMAGES
+    psphot_outroot = psMetadataLookupStr(&status, config->arguments, "OUTPUT");
+    {
+        // Save the background subtracted image
+        psString fn = NULL;
+        psStringAppend(&fn, "%s.p%d.src.sub.fits", psphot_outroot, pass1 ? 1 : 2);
+        psphotSaveImage(0, readout->image, fn);
+        psFree(fn);
+        fn = NULL;
+        if (KRON_SMOOTH) {
+            psStringAppend(&fn, "%s.p%d.src.sub.sm.fits", psphot_outroot, pass1 ? 1 : 2);
+            psphotSaveImage(0, smoothedImage, fn);
+            psFree(fn);
+        }
+        if (KRON_APPLY_WINDOW) {
+            psStringAppend(&fn, "%s.p%d.window.fits", psphot_outroot, pass1 ? 1 : 2);
+            psphotSaveImage(0, kronWindow, fn);
+            psFree(fn);
+        }
+    }
+#endif
 
     // threaded measurement of the source magnitudes
@@ -157 +340 @@
             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, (psS32) KRON_SMOOTH, 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.");
@@ -192 +376 @@
     psFree (cellGroups);
     psFree (kronWindow);
+    psFree (smoothedImage);
+    if (pass1) {
+        pass1 = false;
+    }
 
     psLogMsg ("psphot.kron", PS_LOG_WARN, "measure masked kron magnitudes : %f sec for %ld objects\n", psTimerMark ("psphot.kron"), sources->n);
@@ -209 +397 @@
     bool KRON_APPLY_WEIGHT          = PS_SCALAR_VALUE(job->args->data[8],S32);
     bool KRON_APPLY_WINDOW          = PS_SCALAR_VALUE(job->args->data[9],S32);
-    bool KRON_SMOOTH                = PS_SCALAR_VALUE(job->args->data[10],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++) {
@@ -224 +414 @@
             // 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;
             }
@@ -251 +449 @@
 
             // 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
+            // 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)) {
 
+#ifdef SAVE_IMAGES
+#ifdef SAVE_SOURCE_IMAGES
+                if (j == 1 && !pass1) {
+                    char fn[80];
+                    if (smoothedPixels) {
+                        sprintf(fn, "%s.s.%d.p%d.sm.fits", psphot_outroot, source->id, pass1 ? 1 : 2);
+                        psphotSaveImage(0, smoothedPixels, fn);
+                    }
+                    sprintf(fn, "%s.s.%d.p%d.fits", psphot_outroot, source->id, pass1 ? 1 : 2);
+                    psphotSaveImage(0, source->pixels, fn);
+                }
+#endif
+#endif
+
                 // this function populates moments->Mrf,KronFlux,KronFluxErr
-                psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, KRON_SMOOTH);
+                psphotKronWindowMag (source, kronWindow, windowRadius, MIN_KRON_RADIUS, maskVal, KRON_APPLY_WEIGHT, smoothedPixels);
                 psImageMaskPixels (source->maskObj, "AND", PS_NOT_IMAGE_MASK(markVal));
 
@@ -267 +489 @@
             // 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 applyWeight, bool smooth) {
+    bool applyWeight, psImage *smoothedPixels) {
 
     PS_ASSERT_PTR_NON_NULL(source, false);
@@ -296 +528 @@
     // Xn  = SUM (x - xc)^n * (z - sky)
 
+
     psF32 RF = 0.0;
     psF32 RS = 0.0;
@@ -315 +548 @@
 
     psF32 **vPix;
-    psImage *smoothedImage = NULL;;
-    if (smooth) {
-        smoothedImage = psImageCopy(NULL, source->pixels, PS_TYPE_F32);
-        psImageSmooth(smoothedImage, 1.7, 2);
-        vPix = smoothedImage->data.F32;
+    if (smoothedPixels) {
+        vPix = smoothedPixels->data.F32;
     } else {
         vPix = source->pixels->data.F32;
@@ -376 +606 @@
             psF32 rs = 0.5 * (fDiff1 + fDiff2);
 
+#ifdef BE_SIMPLE_MINDED
+            rf = vPix[row][col] * sqrt(r2);
+            rs = vPix[row][col];
+#endif
+
             RF  += rf;
             RS  += rs;
@@ -386 +621 @@
         Mrf = MIN (radius, Mrf);
     }
+#ifdef USE_SE_KR
+    if (source->tmpFlags & SE_KR_USED) {
+        // we are using the sextratcor KR for this source set it back
+        Mrf = source->moments->Mrf;
+    }
+#endif
 
     // Calculate the Kron magnitude (make this block optional?)
@@ -399 +640 @@
     // smoothed image above)
     vPix = source->pixels->data.F32;
+
 
     for (psS32 row = 0; row < source->pixels->numRows ; row++) {
@@ -433 +675 @@
     source->moments->KronFlux    = Sum;
     source->moments->KronFluxErr = sqrt(Var);
-    if (smoothedImage) {
-        psFree(smoothedImage);
-    }
 
     return true;
@@ -451 +690 @@
     psAssert(kronWindow, "need a window");
 
-    if (source->moments->Mrf < 0) return false;
+    if (!isfinite(source->moments->Mrf) || source->moments->Mrf < 0 ) return false;
 
     if (!applyWindow) {
@@ -478 +717 @@
     }
 
-    float scale = fabs(0.5 * source->moments->Mrf / axes.major);
+    // Why this factor of 0.5 ?
+    float scale = 0.5 * source->moments->Mrf / axes.major;
     axes.major *= scale;
     axes.minor *= scale;
@@ -489 +729 @@
     }
 
-    float Sxx = PS_SQR(shape.sx);
-    float Syy = PS_SQR(shape.sy);
-    float Sxy = -1. * shape.sxy;
     float Smajor = axes.major;
 
@@ -499 +736 @@
     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++) {