Changeset 12941

Timestamp:

Apr 20, 2007, 2:02:35 PM (19 years ago)

Author:

magnier

Message:

added residual subtraction code to pmModelAdd

File:

• branches/eam_02_branch/psModules/src/objects/pmModel.c (modified) (9 diffs)

branches/eam_02_branch/psModules/src/objects/pmModel.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.9 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2006-10-31 19:38:44 $
+ *  @version $Revision: 1.9.6.1 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-21 00:02:35 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -21 +21 @@
 #include <string.h>
 #include <pslib.h>
+#include "pmResiduals.h"
 #include "pmModel.h"
 
@@ -49 +50 @@
 {
     psTrace("psModules.objects", 3, "---- %s() begin ----\n", __func__);
+
     pmModel *tmp = (pmModel *) psAlloc(sizeof(pmModel));
+    psMemSetDeallocator(tmp, (psFreeFunc) modelFree);
 
     tmp->type = type;
@@ -57 +60 @@
     tmp->radiusFit = 0;
     tmp->status = PM_MODEL_UNTRIED;
+    tmp->residuals = NULL;              // XXX should the model own this memory?
 
     psS32 Nparams = pmModelParameterCount(type);
@@ -72 +76 @@
     }
 
-    psMemSetDeallocator(tmp, (psFreeFunc) modelFree);
     psTrace("psModules.objects", 3, "---- %s() end ----\n", __func__);
     return(tmp);
@@ -88 +91 @@
     new->status = model->status;
     new->radiusFit = model->radiusFit;
+
+    // XXX note that model->residuals is just a reference
+    new->residuals = model->residuals;
 
     for (int i = 0; i < new->params->n; i++) {
@@ -147 +153 @@
     psF32 skyValue = params->data.F32[0];
     psF32 pixelValue;
-
-    for (psS32 i = 0; i < image->numRows; i++) {
-        for (psS32 j = 0; j < image->numCols; j++) {
-            if ((mask != NULL) && mask->data.U8[i][j])
+    
+    float xCenter = model->params->data.F32[PM_PAR_XPOS];
+    float yCenter = model->params->data.F32[PM_PAR_YPOS];
+    float Io = model->params->data.F32[PM_PAR_I0];
+
+    int xBin = 1;
+    int yBin = 1;
+    float xResidCenter = 0.0;
+    float yResidCenter = 0.0;
+
+    psImageInterpolateOptions *Ro = NULL;
+    psImageInterpolateOptions *Rx = NULL;
+    psImageInterpolateOptions *Ry = NULL;
+    if (model->residuals) {
+        Ro = psImageInterpolateOptionsAlloc(
+            PS_INTERPOLATE_BILINEAR,
+            model->residuals->Ro, NULL, NULL, 0, 0.0, 0.0, 1, 0, 0.0);
+        Rx = psImageInterpolateOptionsAlloc(
+            PS_INTERPOLATE_BILINEAR,
+            model->residuals->Rx, NULL, NULL, 0, 0.0, 0.0, 1, 0, 0.0);
+        Ry = psImageInterpolateOptionsAlloc(
+            PS_INTERPOLATE_BILINEAR,
+            model->residuals->Ry, NULL, NULL, 0, 0.0, 0.0, 1, 0, 0.0);
+
+        xBin = model->residuals->xBin;
+        yBin = model->residuals->yBin;
+        xResidCenter = model->residuals->xCenter;
+        yResidCenter = model->residuals->yCenter;
+    }
+
+    for (psS32 iy = 0; iy < image->numRows; iy++) {
+        for (psS32 ix = 0; ix < image->numCols; ix++) {
+            if ((mask != NULL) && mask->data.U8[iy][ix])
                 continue;
 
@@ -156 +191 @@
             // 'center' option changes meaning of i,j:
             if (center) {
-                imageCol = j - 0.5*image->numCols + model->params->data.F32[2];
-                imageRow = i - 0.5*image->numRows + model->params->data.F32[3];
+                imageCol = ix - 0.5*image->numCols + xCenter;
+                imageRow = iy - 0.5*image->numRows + yCenter;
             } else {
-                imageCol = j + image->col0;
-                imageRow = i + image->row0;
+                imageCol = ix + image->col0;
+                imageRow = iy + image->row0;
             }
 
@@ -173 +208 @@
             }
 
+            // get the contribution from the residual model
+            // XXX for a test, do this for all sources and all pixels
+            if (Ro) {
+                // fractional image position
+                // this is wrong for the 'center' case
+                float ox = xBin*(ix + 0.5 + image->col0 - xCenter) + xResidCenter;
+                float oy = yBin*(iy + 0.5 + image->row0 - yCenter) + yResidCenter;
+
+                psU8 mflux = 0;
+                double Fo = 0.0;
+                double Fx = 0.0;
+                double Fy = 0.0;
+                psImageInterpolate (&Fo, NULL, &mflux, ox, oy, Ro);
+                psImageInterpolate (&Fx, NULL, &mflux, ox, oy, Rx);
+                psImageInterpolate (&Fy, NULL, &mflux, ox, oy, Ry);
+
+                if (!mflux && isfinite(Fo) && isfinite(Fx) && isfinite(Fy)) {
+                    double flux = Fo + xCenter*Fx + yCenter*Fy;
+                    pixelValue += Io*flux;
+                }
+            }
 
             // add or subtract the value
-            if (add
-               ) {
-                image->data.F32[i][j] += pixelValue;
-            }
-            else {
-                image->data.F32[i][j] -= pixelValue;
+            if (add) {
+                image->data.F32[iy][ix] += pixelValue;
+            } else {
+                image->data.F32[iy][ix] -= pixelValue;
             }
         }
     }
     psFree(x);
+    psFree(Ro);
+    psFree(Rx);
+    psFree(Ry);
     psTrace("psModules.objects", 3, "---- %s(true) end ----\n", __func__);
     return(true);