Changeset 14907

Timestamp:

Sep 20, 2007, 9:17:16 AM (19 years ago)

Author:

magnier

Message:

added psImageMapFit psImageMapClipFit

Location:

branches/eam_branch_20070830/psLib/src/imageops

Files:

• psImageMap.h (modified) (2 diffs)
• psImageMapFit.c (modified) (6 diffs)

branches/eam_branch_20070830/psLib/src/imageops/psImageMap.h

@@ -7 +7 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-17 20:17:21 $
+ *  @version $Revision: 1.1.2.6 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-20 19:17:16 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -50 +50 @@
 
 // fit the image map to a set of points
-bool psImageMapFit (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df);
+bool psImageMapFit (psImageMap *map, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df);
+
+// fit the image map to a set of points
+bool psImageMapClipFit (psImageMap *map, psStats *stats, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df);
 
 /// @}

branches/eam_branch_20070830/psLib/src/imageops/psImageMapFit.c

@@ -7 +7 @@
  *  @author Eugene Magnier, IfA
  *
- *  @version $Revision: 1.1.2.4 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-17 20:19:36 $
+ *  @version $Revision: 1.1.2.5 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-20 19:17:16 $
  *
  *  Copyright 2007 Institute for Astronomy, University of Hawaii
@@ -20 +20 @@
 #include "psError.h"
 #include "psAbort.h"
+#include "psTrace.h"
 
 // XXX for testing
@@ -44 +45 @@
 
 // map defines the output image dimensions and scaling.
-bool psImageMapFit (psImageMap *map, psVector *x, psVector *y, psVector *f, psVector *df) {
+bool psImageMapFit (psImageMap *map, psVector *mask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df) {
 
     int I, J;
+
+    // XXX Add Asserts
 
     // dimensions of the output map image
     int Nx = map->binning->nXruff;
     int Ny = map->binning->nYruff;
+
+    // no spatial information, just calculate mean & stdev
+    if ((Nx == 1) && (Ny == 1)) {
+        psStats *stats = psStatsAlloc (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV);
+        
+        // XXX does ROBUST_MEDIAN work with weight?
+        psVectorStats (stats, f, NULL, mask, maskValue);
+
+        map->map->data.F32[0][0]   = stats->robustMedian;
+        map->error->data.F32[0][0] = stats->robustStdev;
+        psFree (stats);
+        return true;
+    }
+
+    if (Nx == 1) {
+        psAbort ("un-implemented edge case");
+        goto insert;
+    }
+    if (Ny == 1) {
+        psAbort ("un-implemented edge case");
+        goto insert;
+    }
 
     // set up the redirection table so we can use sA[-1][-1], etc
@@ -104 +129 @@
             // J = (n + jn) + nX*(m + jm)
             for (int i = 0; i < x->n; i++) {
+
+                if (mask && (mask->data.U8[i] & maskValue)) continue;
+
                 // base coordinate offset for this point (x,y) relative to this map element (n,m)
                 // float dx = x->data.F32[i] - psImageBinningGetFineX (map->binning, n + 0.5);
@@ -198 +226 @@
             sA[+1][ 0] = dx_rx_ry_ry + dx_rx_py_py;
             sA[+1][+1] = dx_rx_dy_ry;
-
-            // XXX need to add three additional edge cases:
-            if ((Nx == 1) && (Ny == 1)) {
-                psAbort ("un-implemented edge case");
-                goto insert;
-            }
-            if (Nx == 1) {
-                psAbort ("un-implemented edge case");
-                goto insert;
-            }
-            if (Ny == 1) {
-                psAbort ("un-implemented edge case");
-                goto insert;
-            }
 
         insert:
@@ -299 +313 @@
     return true;
 }
+
+// measure residuals on each pass and clip outliers based on stats
+bool psImageMapClipFit (psImageMap *map, psStats *stats, psVector *inMask, psMaskType maskValue, psVector *x, psVector *y, psVector *f, psVector *df) {
+
+    // XXX add in full PS_ASSERTS
+    assert (map);
+    assert (stats);
+    assert (x);
+    assert (y);
+    assert (f);
+    assert (df);
+
+    // the user supplies one of various stats option pairs,
+    // determine the desired mean and stdev STATS options:
+    psStatsOptions meanOption = stats->options & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN | PS_STAT_ROBUST_MEDIAN | PS_STAT_CLIPPED_MEAN | PS_STAT_FITTED_MEAN | PS_STAT_FITTED_MEAN_V2 | PS_STAT_FITTED_MEAN_V3 | PS_STAT_FITTED_MEAN_V4);
+    psStatsOptions stdevOption = stats->options & (PS_STAT_SAMPLE_STDEV | PS_STAT_ROBUST_STDEV | PS_STAT_CLIPPED_STDEV | PS_STAT_FITTED_STDEV | PS_STAT_FITTED_STDEV_V2 | PS_STAT_FITTED_STDEV_V3 | PS_STAT_FITTED_STDEV_V4);
+    if (!meanOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid mean stats option selected");
+        return false;
+    }
+    if (!stdevOption) {
+        psError(PS_ERR_UNKNOWN, true, "no valid stdev stats option selected");
+        return false;
+    }
+
+    // clipping range defined by min and max and/or clipSigma
+    psF32 minClipSigma;
+    psF32 maxClipSigma;
+    if (isfinite(stats->max)) {
+        maxClipSigma = fabs(stats->max);
+    } else {
+        maxClipSigma = fabs(stats->clipSigma);
+    }
+    if (isfinite(stats->min)) {
+        minClipSigma = fabs(stats->min);
+    } else {
+        minClipSigma = fabs(stats->clipSigma);
+    }
+
+    psVector *mask = inMask;
+    if (!inMask) {
+        mask = psVectorAlloc (x->n, PS_TYPE_U8);
+    }
+
+    // vector to store residuals
+    psVector *resid = psVectorAlloc(f->n, PS_TYPE_F32);
+
+    psTrace("psLib.imageops", 4, "stats->clipIter is %d\n", stats->clipIter);
+    psTrace("psLib.imageops", 4, "(minClipSigma, maxClipSigma) is (%.2f, %.2f)\n", minClipSigma, maxClipSigma);
+
+    for (psS32 N = 0; N < stats->clipIter; N++) {
+        psTrace("psLib.imageops", 6, "Loop iteration %d.  Calling psImageMapFit()\n", N);
+        psS32 Nkeep = 0;
+        if (!psImageMapFit(map, mask, maskValue, x, y, f, df)) {
+            psError(PS_ERR_UNKNOWN, false, "Could not fit image map.\n");
+            psFree(resid);
+            if (!inMask) psFree (mask);
+            return false;
+        }
+
+        psVector *fit = psImageMapEvalVector(map, x, y);
+        if (fit == NULL) {
+            psError(PS_ERR_UNKNOWN, false, "Failure in psImageMapEvalVector().\n");
+            psFree(resid);
+            if (!inMask) psFree (mask);
+            return false;
+        }
+        for (int i = 0 ; i < f->n ; i++) {
+            resid->data.F32[i] = (f->data.F32[i] - fit->data.F32[i]);
+        }
+
+        if (!psVectorStats(stats, resid, NULL, mask, maskValue)) {
+            psError(PS_ERR_UNKNOWN, false, "Failure to compute statistics on the resid vector.\n");
+            psFree(resid);
+            psFree(fit);
+            if (!inMask) psFree (mask);
+            return false;
+        }
+
+        double meanValue = psStatsGetValue (stats, meanOption);
+        double stdevValue = psStatsGetValue (stats, stdevOption);
+
+        psTrace("psLib.imageops", 5, "Mean is %f\n", meanValue);
+        psTrace("psLib.imageops", 5, "Stdev is %f\n", stdevValue);
+        psF32 minClipValue = -minClipSigma*stdevValue;
+        psF32 maxClipValue = +maxClipSigma*stdevValue;
+
+        // set mask if pts are not valid
+        // we are masking out any point which is out of range
+        // recovery is not allowed with this scheme
+        for (psS32 i = 0; i < resid->n; i++) {
+            // XXX this prevents recovery of previously masked values
+            if (mask->data.U8[i] & maskValue) {
+                continue;
+            }
+
+            if ((resid->data.F32[i] - meanValue > maxClipValue) || (resid->data.F32[i] - meanValue < minClipValue)) {
+                psTrace("psLib.imageops", 6, "Masking element %d  : %f vs %f : resid is %f\n", i, f->data.F32[i], fit->data.F32[i], resid->data.F32[i]);
+                mask->data.U8[i] |= 0x01;
+                continue;
+            }
+            Nkeep++;
+        }
+
+        // We should probably exit this loop if no new elements were masked since the fit won't
+        // change.
+        psTrace("psLib.imageops", 6, "keeping %d of %ld pts for fit\n", Nkeep, x->n);
+        stats->clippedNvalues = Nkeep;
+        psFree(fit);
+    }
+
+    // Free local temporary variables
+    psFree(resid);
+    if (!inMask) psFree (mask);
+    return true;
+}