Changeset 33879

Timestamp:

May 16, 2012, 4:23:42 PM (14 years ago)

Author:

bills

Message:

Add pmSourceSmoothOp which adds or subtracts a smoothed version of a source's
model.
Remove instrumentation from psphotKronIterate.c

Location:

trunk

Files:

• psModules/src/objects/pmSource.c (modified) (3 diffs)
• psModules/src/objects/pmSource.h (modified) (2 diffs)
• psphot/src/psphotKronIterate.c (modified) (10 diffs)

trunk/psModules/src/objects/pmSource.c

@@ -39 +39 @@
 #include "pmSourceExtendedPars.h"
 #include "pmSourceDiffStats.h"
+#include "pmSourcePhotometry.h"
 #include "pmSource.h"
 
@@ -159 +160 @@
     source->radialAper = NULL;
     source->parent = NULL;
+    source->tmpPtr = NULL;
     source->imageID = -1;
     source->nFrames = 0;
@@ -1157 +1159 @@
     PAR[PM_PAR_SYY] = oldshape.sy;
     PAR[PM_PAR_SXY] = oldshape.sxy;
+
+    return true;
+}
+
+bool pmSourceSmoothOp (pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy)
+{
+//    assert (mode & PM_MODEL_OP_NOISE);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_PTR_NON_NULL(target, false);
+
+    if (add) {
+        psTrace ("psphot", 3, "adding smoothed object at %f,%f\n", source->peak->xf, source->peak->yf);
+    } else {
+        psTrace ("psphot", 3, "removing smooted object at %f,%f\n", source->peak->xf, source->peak->yf);
+    }
+
+    pmModel *model = pmSourceGetModel(NULL, source);
+    if (!model) {
+        return false;
+    }
+    pmSourceSmoothOpModel (model, source, mode, target, sigma, add, maskVal, dx, dy);
+
+    return true;
+}
+
+bool pmSourceSmoothOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy) 
+{
+    bool status;
+    psEllipseShape oldshape;
+    psEllipseShape newshape;
+    psEllipseAxes axes;
+
+    if (add) {
+        psTrace ("psphot", 4, "adding smoothed object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+    } else {
+        psTrace ("psphot", 4, "removing smoothed object at %f,%f\n", model->params->data.F32[PM_PAR_XPOS], model->params->data.F32[PM_PAR_YPOS]);
+    }
+
+    psF32 *PAR = model->params->data.F32;
+
+    // Isn't this hanging around somewhere?
+    float oldFlux = NAN;
+    if (!pmSourcePhotometryModel (NULL, &oldFlux, model)) return false;
+
+    // save original values
+    float oldI0  = PAR[PM_PAR_I0];
+    float oldsxx = PAR[PM_PAR_SXX];
+    float oldsyy = PAR[PM_PAR_SYY];
+    float oldsxy = PAR[PM_PAR_SXY];
+
+    // Since we are going to scale the flux correctly we need to get our
+    // factors of sqrt(2) right
+    oldshape.sx  = oldsxx / M_SQRT2;
+    oldshape.sy  = oldsyy / M_SQRT2;
+    oldshape.sxy = oldsxy;
+
+    // XXX can this be done more intelligently?
+    if (oldI0 == 0.0) return false;
+    if (!isfinite(oldI0)) return false;
+
+    // increase size and height of source
+    axes = psEllipseShapeToAxes (oldshape, 20.0);
+    axes.major = sqrt(PS_SQR(axes.major) + PS_SQR(sigma));
+    axes.minor = sqrt(PS_SQR(axes.minor) + PS_SQR(sigma));
+    newshape = psEllipseAxesToShape (axes);
+    PAR[PM_PAR_SXX] = newshape.sx * M_SQRT2;
+    PAR[PM_PAR_SYY] = newshape.sy * M_SQRT2;
+    PAR[PM_PAR_SXY] = newshape.sxy;
+
+    PAR[PM_PAR_I0]  = 1.0;
+
+    float newFlux;
+    if (!pmSourcePhotometryModel (NULL, &newFlux, model)) return false;
+
+    PAR[PM_PAR_I0]  = oldFlux / newFlux;
+
+    if (add) {
+        status = pmModelAddWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
+    } else {
+        status = pmModelSubWithOffset (target, source->maskObj, model, mode, maskVal, dx, dy);
+    }
+
+    // restore original values
+    PAR[PM_PAR_I0]  = oldI0;
+    PAR[PM_PAR_SXX] = oldsxx;
+    PAR[PM_PAR_SYY] = oldsyy;
+    PAR[PM_PAR_SXY] = oldsxy;
 
     return true;

trunk/psModules/src/objects/pmSource.h

@@ -117 +117 @@
     psArray *radialAper;                ///< radial flux in circular apertures
     pmSource *parent;                   ///< reference to the master source from which this is derived
+    psPtr *tmpPtr;                      ///< pointer that may be used to store data in a particular module. e.g. psphotKronIterate.
     int imageID;
     psU16 nFrames;
@@ -302 +303 @@
 bool pmSourceNoiseOp (pmSource *source, pmModelOpMode mode, float FACTOR, float SIZE, bool add, psImageMaskType maskVal, int dx, int dy);
 
+bool pmSourceSmoothOp (pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy);
+bool pmSourceSmoothOpModel (pmModel *model, pmSource *source, pmModelOpMode mode, psImage *target, float sigma, bool add, psImageMaskType maskVal, int dx, int dy);
+
 bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add, psImageMaskType maskVal, int dx, int dy);
 bool pmSourceCacheModel (pmSource *source, psImageMaskType maskVal);

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);
@@ -57 +54 @@
 bool psphotVisualRangeImage (int kapaFD, psImage *inImage, const char *name, int channel, float min, float max);
 
-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) {
 
@@ -165 +65 @@
     psTimerStart ("psphot.kron");
 
-#ifdef USE_SE_KR
-    if (!pass1) {
-        build_kr_list(sources);
-    }
-#endif
 
     // determine the number of allowed threads
@@ -237 +132 @@
 
     // 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++) {
 
@@ -247 +139 @@
         // (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
@@ -291 +173 @@
 
     }
-
-#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
@@ -377 +237 @@
     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);
@@ -465 +322 @@
             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, smoothedPixels);
@@ -606 +449 @@
             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;
@@ -621 +459 @@
         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?)