Changeset 14778 for branches/eam_branch_20070830/psModules/src/objects/pmTrend2D.c

Timestamp:

Sep 7, 2007, 10:24:34 AM (19 years ago)

Author:

magnier

Message:

adding fitting, evaluation functions

File:

• branches/eam_branch_20070830/psModules/src/objects/pmTrend2D.c (modified) (6 diffs)

branches/eam_branch_20070830/psModules/src/objects/pmTrend2D.c

@@ -3 +3 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1.2.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-09-04 18:35:51 $
+ *  @version $Revision: 1.1.2.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-09-07 20:24:34 $
  *
  *  Copyright 2004 Institute for Astronomy, University of Hawaii
@@ -22 +22 @@
         return;
 
+    psFree (trend->stats);
     psFree (trend->poly);
     psFree (trend->map);
@@ -27 +28 @@
 }
 
-pmTrend2D *pmTrend2DAlloc (pmTrend2DMode mode, int nXout, int nYout, int nXtrend, int nYtrend, psStats *stats)
+pmTrend2D *pmTrend2DAlloc (pmTrend2DMode mode, psImage *image, int nXtrend, int nYtrend, psStats *stats)
 {
     pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
@@ -35 +36 @@
     trend->poly = NULL;
     trend->stats = psMemIncrRefCounter (stats);
+    trend->mode = mode;
         
     switch (mode) {
@@ -45 +47 @@
         break;
 
-      case PM_TREND_MAP:
-        // XXX *** I need the output image dimensions here!
+      case PM_TREND_MAP: {
+          // binning defines the map scale relationship
+          psImageBinning *binning = psImageBinningAlloc();
+          binning->nXfine = image->numCols;
+          binning->nYfine = image->numRows;
+          binning->nXruff = nXtrend;
+          binning->nYruff = nYtrend;
 
-        // function is defined over the range 0-1000, 0-1000
-        psImageBinning *binning = psImageBinningAlloc();
-        binning->nXfine = nXout;
-        binning->nYfine = nYout;
-        binning->nXruff = nXtrend;
-        binning->nYruff = nYtrend;
-
-        trend->map = psImageMapAlloc (NULL, binning, stats);
-        psFree (binning);
-        break;
+          trend->map = psImageMapAlloc (image, binning, stats);
+          psFree (binning);
+          break;
+      }
 
       default:
@@ -64 +65 @@
     return (trend);
 }
+
+pmTrend2D *pmTrend2DFieldAlloc (pmTrend2DMode mode, int nXfield, int nYfield, int nXtrend, int nYtrend, psStats *stats)
+{
+    pmTrend2D *trend = (pmTrend2D *) psAlloc(sizeof(pmTrend2D));
+    psMemSetDeallocator(trend, (psFreeFunc) pmTrend2DFree);
+
+    trend->map = NULL;
+    trend->poly = NULL;
+    trend->stats = psMemIncrRefCounter (stats);
+    trend->mode = mode;
+        
+    switch (mode) {
+      case PM_TREND_POLY_ORD:
+        trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_ORD, nXtrend, nYtrend);
+        break;
+
+      case PM_TREND_POLY_CHEB:
+        trend->poly = psPolynomial2DAlloc (PS_POLYNOMIAL_CHEB, nXtrend, nYtrend);
+        break;
+
+      case PM_TREND_MAP: {
+          // binning defines the map scale relationship
+          psImageBinning *binning = psImageBinningAlloc();
+          binning->nXfine = nXfield;
+          binning->nYfine = nYfield;
+          binning->nXruff = nXtrend;
+          binning->nYruff = nYtrend;
+
+          trend->map = psImageMapAlloc (NULL, binning, stats);
+          psFree (binning);
+          break;
+      }
+
+      default:
+        psAbort ("error");
+    }
+    return (trend);
+}
+
+bool pmTrend2DFit (pmTrend2D *trend, psVector *x, psVector *y, psVector *f, psVector *df) {
+
+    bool status;
+
+    psVector *mask = psVectorAlloc (x->n, PS_TYPE_MASK);
+
+    switch (trend->mode) {
+      case PM_TREND_POLY_ORD:
+      case PM_TREND_POLY_CHEB:
+        status = psVectorClipFitPolynomial2D (trend->poly, trend->stats, mask, 0xff, f, df, x, y);
+        // we can use the API here which adjusts the polynomial order based on the number
+        // of points in the image, and potentially based on the fractional range of the
+        // data?
+        break;
+
+      case PM_TREND_MAP:
+        // XXX supply fraction from trend elements
+        status = psImageMapGenerateScale (trend->map, x, y, f, df, 0.1);
+        break;
+
+      default:
+        psAbort ("error");
+    }
+    psFree (mask);
+    return status;
+}
+
+double pmTrend2DEval (pmTrend2D *trend, float x, float y) {
+
+    double result;
+
+    if (!trend) return 0.0;
+
+    switch (trend->mode) {
+      case PM_TREND_POLY_ORD:
+      case PM_TREND_POLY_CHEB:
+        result = psPolynomial2DEval (trend->poly, x, y);
+        break;
+
+      case PM_TREND_MAP:
+        result = psImageMapEval (trend->map, x, y);
+        break;
+
+      default:
+        psAbort ("error");
+    }
+    return result;
+}
+
+psVector *pmTrend2DEvalVector (pmTrend2D *trend, psVector *x, psVector *y) {
+
+    psVector *result;
+
+    switch (trend->mode) {
+      case PM_TREND_POLY_ORD:
+      case PM_TREND_POLY_CHEB:
+        result = psPolynomial2DEvalVector (trend->poly, x, y);
+        break;
+
+      case PM_TREND_MAP:
+        result = psImageMapEvalVector (trend->map, x, y);
+        break;
+
+      default:
+        psAbort ("error");
+    }
+    return result;
+}