Changeset 34041 for branches/meh_branches/ppstack_test/psModules/src/objects/pmSource.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)
• psModules (modified) (1 prop)
• psModules/src/objects/pmSource.c (modified) (13 diffs)

branches/meh_branches/ppstack_test/psModules/src/objects/pmSource.c

@@ -39 +39 @@
 #include "pmSourceExtendedPars.h"
 #include "pmSourceDiffStats.h"
+#include "pmSourcePhotometry.h"
 #include "pmSource.h"
 
@@ -50 +51 @@
     psFree(tmp->pixels);
     psFree(tmp->variance);
+# if (HAVE_MODEL_VAR)
+    psFree(tmp->modelVar);
+# endif
     psFree(tmp->maskObj);
     psFree(tmp->maskView);
@@ -76 +80 @@
     psFree (source->pixels);
     psFree (source->variance);
+# if (HAVE_MODEL_VAR)
+    psFree (source->modelVar);
+# endif
     psFree (source->maskObj);
     psFree (source->maskView);
@@ -83 +90 @@
     source->pixels = NULL;
     source->variance = NULL;
+# if (HAVE_MODEL_VAR)
+    source->modelVar = NULL;
+# endif
     source->maskObj = NULL;
     source->maskView = NULL;
@@ -112 +122 @@
     source->pixels = NULL;
     source->variance = NULL;
+# if (HAVE_MODEL_VAR)
+    source->modelVar = NULL;
+# endif
     source->maskObj = NULL;
     source->maskView = NULL;
@@ -159 +172 @@
     source->radialAper = NULL;
     source->parent = NULL;
+    source->tmpPtr = NULL;
     source->imageID = -1;
+    source->nFrames = 0;
 
     psMemSetDeallocator(source, (psFreeFunc) sourceFree);
@@ -197 +212 @@
     source->pixels   = in->pixels   ? psImageCopyView(NULL, in->pixels)   : NULL;
     source->variance = in->variance ? psImageCopyView(NULL, in->variance) : NULL;
+# if (HAVE_MODEL_VAR)
+    source->modelVar = NULL;
+# endif
     source->maskView = in->maskView ? psImageCopyView(NULL, in->maskView) : NULL;
 
@@ -1026 +1044 @@
     bool addNoise = mode & PM_MODEL_OP_NOISE;
 
+# if (HAVE_MODEL_VAR)
+    bool addModelVar = mode & PM_MODEL_OP_MODELVAR;
+    if (addModelVar) psAssert (source->modelVar, "programming error");
+# endif
+
     // require the use of pmModelAddWithOffset if we are adding noise (because the model size and norm are rescaled)
     if (!addNoise && source->modelFlux) {
@@ -1047 +1070 @@
         }
 
+# if (HAVE_MODEL_VAR)
+        psF32 **target = addModelVar ? source->modelVar->data.F32 : source->pixels->data.F32;
+# else
         psF32 **target = source->pixels->data.F32;
+# endif
 
         for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
@@ -1062 +1089 @@
             }
         }
+# if (HAVE_MODEL_VAR)
+        if (!addModelVar) {
+            if (add) {
+                source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+            } else {
+                source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+            }
+        }
+# else
         if (add) {
             source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
@@ -1067 +1103 @@
             source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
         }
+# endif
         return true;
     }
@@ -1074 +1111 @@
         target = source->variance;
     }
+# if (HAVE_MODEL_VAR)
+    if (addModelVar) {
+        target = source->modelVar;
+    }
+# endif
 
     if (add) {
         status = pmModelAddWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+# if (HAVE_MODEL_VAR)
+        if (!addNoise && !addModelVar) source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+# else
         source->tmpFlags &= ~PM_SOURCE_TMPF_SUBTRACTED;
+# endif
     } else {
         status = pmModelSubWithOffset (target, source->maskObj, model, PM_MODEL_OP_FULL, maskVal, dx, dy);
+# if (HAVE_MODEL_VAR)
+        if (!addNoise && !addModelVar) source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+# else
         source->tmpFlags |= PM_SOURCE_TMPF_SUBTRACTED;
+# endif
     }
 
@@ -1156 +1206 @@
     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;