Changeset 6439

Timestamp:

Feb 16, 2006, 4:48:25 PM (20 years ago)

Author:

magnier

Message:

added local sky variance, modified to work with the new psMinimizeLMChi2

Location:

branches/eam_rel9_p0/psModules/src/objects

Files:

• pmObjects.c (modified) (6 diffs)
• pmObjects.h (modified) (3 diffs)

branches/eam_rel9_p0/psModules/src/objects/pmObjects.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.5.4.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-08 07:16:49 $
+ *  @version $Revision: 1.5.4.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-17 02:48:25 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -786 +786 @@
         return(false);
     }
-    source->moments = pmMomentsAlloc();
+    if (source->moments == NULL) {
+        source->moments = pmMomentsAlloc();
+    }
     source->moments->Sky = (psF32) tmpF64;
+    psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
+    return (true);
+}
+
+// A complementary function to pmSourceLocalSky: calculate the local median variance
+bool pmSourceLocalSkyVariance(
+    pmSource *source,
+    psStatsOptions statsOptions,
+    psF32 Radius)
+{
+    psTrace(__func__, 3, "---- %s() begin ----\n", __func__);
+    PS_ASSERT_PTR_NON_NULL(source, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->weight, false);
+    PS_ASSERT_IMAGE_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->peak, false);
+    PS_ASSERT_INT_POSITIVE(Radius, false);
+    PS_ASSERT_INT_NONNEGATIVE(Radius, false);
+
+    psImage *image = source->weight;
+    psImage *mask  = source->mask;
+    pmPeak *peak  = source->peak;
+    psRegion srcRegion;
+
+    srcRegion = psRegionForSquare(peak->x, peak->y, Radius);
+    srcRegion = psRegionForImage(mask, srcRegion);
+
+    psImageMaskRegion(mask, srcRegion, "OR", PSPHOT_MASK_MARKED);
+    psStats *myStats = psStatsAlloc(statsOptions);
+    myStats = psImageStats(myStats, image, mask, 0xff);
+    psImageMaskRegion(mask, srcRegion, "AND", ~PSPHOT_MASK_MARKED);
+
+    psF64 tmpF64;
+    p_psGetStatValue(myStats, &tmpF64);
+    psFree(myStats);
+
+    if (isnan(tmpF64)) {
+        psTrace(__func__, 3, "---- %s(false) end ----\n", __func__);
+        return(false);
+    }
+    if (source->moments == NULL) {
+        source->moments = pmMomentsAlloc();
+    }
+    source->moments->dSky = (psF32) tmpF64;
     psTrace(__func__, 3, "---- %s(true) end ----\n", __func__);
     return (true);
@@ -1337 +1382 @@
     PS_ASSERT_PTR_NON_NULL(source->moments, false);
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
-    PS_FLOAT_COMPARE(0.0, radius, false);
+    PS_ASSERT_FLOAT_LARGER_THAN_OR_EQUAL(radius, 0.0, false);
 
     //
@@ -1671 +1716 @@
     yErr->n = nPix;
 
+    // XXX EAM : the new minimization API supplies the constraints as a struct
     psMinimization *myMin = psMinimizationAlloc(PM_SOURCE_FIT_MODEL_NUM_ITERATIONS,
                             PM_SOURCE_FIT_MODEL_TOLERANCE);
+    psMinConstrain *constrain = psMinConstrainAlloc();
 
     // PSF model only fits first 4 parameters, EXT model fits all
@@ -1684 +1731 @@
         }
     }
+    constrain->paramMask = paramMask;
 
     // Set the parameter range checks
     pmModelLimits modelLimits = pmModelLimits_GetFunction (model->type);
-    psVector *beta_lim = NULL;
-    psVector *params_min = NULL;
-    psVector *params_max = NULL;
-
-    // XXX EAM : in this implementation, I pass in the limits with the covar matrix.
-    //           in the SDRS, I define a new psMinimization which will take these in
-    psImage *covar = psImageAlloc (params->n, 3, PS_TYPE_F64);
-    modelLimits (&beta_lim, &params_min, &params_max);
-    for (int i = 0; i < params->n; i++) {
-        covar->data.F64[0][i] = beta_lim->data.F32[i];
-        covar->data.F64[1][i] = params_min->data.F32[i];
-        covar->data.F64[2][i] = params_max->data.F32[i];
-    }
+    modelLimits (&constrain->paramDelta, &constrain->paramMin, &constrain->paramMax);
+
+    psImage *covar = psImageAlloc (params->n, params->n, PS_TYPE_F64);
 
     psTrace (".pmObjects.pmSourceFitModel", 5, "fitting function\n");
-    fitStatus = psMinimizeLMChi2(myMin, covar, params, paramMask, x, y, yErr, modelFunc);
+
+    fitStatus = psMinimizeLMChi2(myMin, covar, params, constrain, x, y, yErr, modelFunc);
     for (int i = 0; i < dparams->n; i++) {
         if ((paramMask != NULL) && paramMask->data.U8[i])
@@ -1749 +1788 @@
     psFree(myMin);
     psFree(covar);
-    psFree(paramMask);
-    psFree(params_min);
-    psFree(params_max);
-    psFree(beta_lim);
+    psFree(constrain->paramMask);
+    psFree(constrain->paramMin);
+    psFree(constrain->paramMax);
+    psFree(constrain->paramDelta);
+    psFree(constrain);
 
     rc = (onPic && fitStatus);

branches/eam_rel9_p0/psModules/src/objects/pmObjects.h

@@ -10 +10 @@
  *  @author GLG, MHPCC
  *
- *  @version $Revision: 1.4.4.7 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-02-08 07:16:49 $
+ *  @version $Revision: 1.4.4.8 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2006-02-17 02:48:23 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -88 +88 @@
 typedef struct
 {
-    float x;    ///< X-coord of centroid.
-    float y;    ///< Y-coord of centroid.
+    float x;     ///< X-coord of centroid.
+    float y;     ///< Y-coord of centroid.
     float Sx;    ///< x-second moment.
     float Sy;    ///< y-second moment.
-    float Sxy;    ///< xy cross moment.
-    float Sum;    ///< Pixel sum above sky (background).
-    float Peak;    ///< Peak counts above sky.
-    float Sky;    ///< Sky level (background).
+    float Sxy;   ///< xy cross moment.
+    float Sum;   ///< Pixel sum above sky (background).
+    float Peak;  ///< Peak counts above sky.
+    float Sky;   ///< Sky level (background).
+    float dSky;  ///< local Sky variance
     float SN;    ///< approx signal-to-noise
-    int nPixels;   ///< Number of pixels used.
+    int nPixels; ///< Number of pixels used.
 }
 pmMoments;
@@ -360 +361 @@
 
 
+// A complementary function to pmSourceLocalSky: calculate the local sky variance
+bool pmSourceLocalSkyVariance(
+    pmSource *source,   ///< The input image (float)
+    psStatsOptions statsOptions, ///< The statistic used in calculating the background sky
+    float Radius   ///< The inner radius of the square annulus to exclude
+);
+
 /** pmSourceMoments()
  *