Changeset 20448

Timestamp:

Oct 28, 2008, 2:01:22 PM (18 years ago)

Author:

eugene

Message:

use a hard-wired 2x2 bilinear interpolation

File:

• trunk/psModules/src/objects/pmModel.c (modified) (7 diffs)

trunk/psModules/src/objects/pmModel.c

@@ -6 +6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.23 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-10-22 02:11:08 $
+ *  @version $Revision: 1.24 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-10-29 00:01:22 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -193 +193 @@
     // use the true source position for the residual model
     // the PSF model has presumably already been set for this coordinate
-    float xPos    = params->data.F32[PM_PAR_XPOS];
-    float yPos    = params->data.F32[PM_PAR_YPOS];
+    float XoSave  = params->data.F32[PM_PAR_XPOS];
+    float YoSave  = params->data.F32[PM_PAR_YPOS];
     float IoSave  = params->data.F32[PM_PAR_I0];
     float skySave = params->data.F32[PM_PAR_SKY];
 
+    // the options allow us to modify various aspects of the model
     if (mode & PM_MODEL_OP_NORM) {
         params->data.F32[PM_PAR_I0] = 1.0;
@@ -208 +209 @@
         params->data.F32[PM_PAR_YPOS] = image->row0 + 0.5*image->numRows;
     }
+
+    // apply optional relative offset
+    params->data.F32[PM_PAR_XPOS] += dx;
+    params->data.F32[PM_PAR_YPOS] += dy;
 
     // use these values for this realization
@@ -216 +221 @@
     int xBin = 1;
     int yBin = 1;
-    float xResidCenter = 0.0;
-    float yResidCenter = 0.0;
-
-    psImage *myRo = NULL;
-    psImage *myRx = NULL;
-    psImage *myRy = NULL;
-    psImageInterpolation *Ro = NULL;
-    psImageInterpolation *Rx = NULL;
-    psImageInterpolation *Ry = NULL;
-    if (model->residuals && (mode & (PM_MODEL_OP_RES0 | PM_MODEL_OP_RES1))) {
-        // if the residual image and object image don't match,
-        // supply an appropriately overlapped residual image
-        psImage *inRo = model->residuals->Ro;
-        psImage *inRx = model->residuals->Rx;
-        psImage *inRy = model->residuals->Ry;
-        xBin = model->residuals->xBin;
-        yBin = model->residuals->yBin;
-        xResidCenter = model->residuals->xCenter;
-        yResidCenter = model->residuals->yCenter;
-        if ((image->numCols != inRo->numCols) ||
-            (image->numRows != inRo->numRows)) {
-            myRo = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32);
-            myRx = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32);
-            myRy = psImageAlloc (image->numCols, image->numRows, PS_TYPE_F32);
-            // Difference between input and desired centres
-            int xDiff = (int)(inRo->numCols / 2) - (xPos - image->col0);
-            int yDiff = (int)(inRo->numRows / 2) - (yPos - image->row0);
-            xResidCenter -= xDiff;
-            yResidCenter -= yDiff;
-            for (int iy = 0; iy < myRo->numRows; iy++) {
-                int jy = iy + yDiff;
-                if ((jy < 0) || (jy >= inRo->numRows)) {
-                    for (int ix = 0; ix < myRo->numCols; ix++) {
-                        myRo->data.F32[iy][ix] = 0.0;
-                        myRx->data.F32[iy][ix] = 0.0;
-                        myRy->data.F32[iy][ix] = 0.0;
-                    }
-                    continue;
-                }
-                for (int ix = 0; ix < myRo->numCols; ix++) {
-                    int jx = ix + xDiff;
-                    if ((jx < 0) || (jx >= inRo->numCols)) {
-                        myRo->data.F32[iy][ix] = 0.0;
-                        myRx->data.F32[iy][ix] = 0.0;
-                        myRy->data.F32[iy][ix] = 0.0;
-                    } else {
-                        myRo->data.F32[iy][ix] = inRo->data.F32[jy][jx];
-                        myRx->data.F32[iy][ix] = inRx->data.F32[jy][jx];
-                        myRy->data.F32[iy][ix] = inRy->data.F32[jy][jx];
-                    }
-                }
-            }
-        } else {
-            myRo = psMemIncrRefCounter (inRo);
-            myRx = psMemIncrRefCounter (inRx);
-            myRy = psMemIncrRefCounter (inRy);
-        }
-
-        Ro = psImageInterpolationAlloc(PS_INTERPOLATE_BILINEAR, myRo, NULL, mask, 0, 0.0, 0.0, 1, 0, 0.0, 0);
-        Rx = psImageInterpolationAlloc(PS_INTERPOLATE_BILINEAR, myRx, NULL, NULL, 0, 0.0, 0.0, 1, 0, 0.0, 0);
-        Ry = psImageInterpolationAlloc(PS_INTERPOLATE_BILINEAR, myRy, NULL, NULL, 0, 0.0, 0.0, 1, 0, 0.0, 0);
-
-    }
+    float DX = 0.0;
+    float DY = 0.0;
+    int NX = 0;
+    int NY = 0;
+
+    psF32 **Ro = NULL;
+    psF32 **Rx = NULL;
+    psF32 **Ry = NULL;
+    psU8  **Rm = NULL;
+
+    if (model->residuals) {
+        DX = xBin*(image->col0 - xCenter) + model->residuals->xCenter + 0.5;
+        DY = yBin*(image->row0 - yCenter) + model->residuals->yCenter + 0.5;
+        Ro = (model->residuals->Ro)   ? model->residuals->Ro->data.F32 : NULL;
+        Rx = (model->residuals->Rx)   ? model->residuals->Rx->data.F32 : NULL;
+        Ry = (model->residuals->Ry)   ? model->residuals->Ry->data.F32 : NULL;
+        Rm = (model->residuals->mask) ? model->residuals->mask->data.U8 : NULL;
+        if (Ro) {
+            NX = model->residuals->Ro->numCols;
+            NY = model->residuals->Ro->numRows;
+        }           
+    }
+
+    // XXX trying to improve the speed and threadability of this function.
+    // note: model->residuals is a view to the item on pmPSF.
 
     for (psS32 iy = 0; iy < image->numRows; iy++) {
@@ -285 +252 @@
                 continue;
 
-            // Convert i/j to image coord space:
             // XXX should we use using 0.5 pixel offset?
-            imageCol = ix + image->col0 + dx;
-            imageRow = iy + image->row0 + dy;
+            // Convert to coordinate in parent image, with offset (dx,dy)
+            imageCol = ix + image->col0;
+            imageRow = iy + image->row0;
 
             x->data.F32[0] = (float) imageCol;
@@ -301 +268 @@
 
             // get the contribution from the residual model
-            if (Ro) {
-                // fractional image position
-                float ox = xBin*(imageCol + 0.5 - xCenter) + xResidCenter;
-                float oy = yBin*(imageRow + 0.5 - yCenter) + yResidCenter;
-
-                psU8 mflux = 0;
-                if (mode & PM_MODEL_OP_RES0) {
-                    double Fo = 0.0;
-                    psImageInterpolate (&Fo, NULL, &mflux, ox, oy, Ro);
-                    if (!mflux && isfinite(Fo)) {
-                        pixelValue += Io*Fo;
-                    }
-                }
-                // skip Rx,Ry if Ro is masked
-                if (!mflux && (mode & PM_MODEL_OP_RES1)) {
-                    double Fx = 0.0;
-                    double Fy = 0.0;
-                    psImageInterpolate (&Fx, NULL, &mflux, ox, oy, Rx);
-                    psImageInterpolate (&Fy, NULL, &mflux, ox, oy, Ry);
-                    if (!mflux && isfinite(Fx) && isfinite(Fy)) {
-                        pixelValue += Io*(xPos*Fx + yPos*Fy);
-                    }
-                }
+            if (Ro && (mode & PM_MODEL_OP_RES0)) {
+                // residual image position
+                float ry = yBin*iy + DY;
+                float rx = xBin*ix + DX;
+
+                int rx0 = rx - 0.5;
+                int rx1 = rx + 0.5;
+                int ry0 = ry - 0.5;
+                int ry1 = ry + 0.5;
+
+                if (rx0 < 0) goto skip;
+                if (ry0 < 0) goto skip;
+                if (rx1 >= NX) goto skip;
+                if (ry1 >= NY) goto skip;
+
+                // these go from 0.0 to 1.0 between the centers of the pixels 
+                float fx = rx - 0.5 - rx0;
+                float Fx = 1.0 - fx;
+                float fy = ry - 0.5 - ry0;
+                float Fy = 1.0 - fy;
+
+                // check the residual image mask (if set). give up if any of the 4 pixels are masked.
+                if (Rm) {
+                    if (Rm[ry0][rx0]) goto skip;
+                    if (Rm[ry0][rx1]) goto skip;
+                    if (Rm[ry1][rx0]) goto skip;
+                    if (Rm[ry1][rx1]) goto skip;
+                }
+
+                // a possible further optimization if we re-use these values
+                // XXX allow for masked pixels, and add pixel weights
+                float V0 = (Ro[ry0][rx0]*Fx + Ro[ry0][rx1]*fx);
+                float V1 = (Ro[ry1][rx0]*Fx + Ro[ry1][rx1]*fx);
+                float Vo = V0*Fy + V1*fy;
+                if (!isfinite(Vo)) goto skip;
+
+                float Vx = 0.0;
+                float Vy = 0.0;
+
+                // skip Rx,Ry if Ro is masked
+                if (Rx && Ry && (mode & PM_MODEL_OP_RES1)) {
+                    V0 = (Rx[ry0][rx0]*Fx + Rx[ry0][rx1]*fx);
+                    V1 = (Rx[ry1][rx0]*Fx + Rx[ry1][rx1]*fx);
+                    Vx = V0*Fy + V1*fy;
+
+                    V0 = (Ry[ry0][rx0]*Fx + Ry[ry0][rx1]*fx);
+                    V1 = (Ry[ry1][rx0]*Fx + Ry[ry1][rx1]*fx);
+                    Vy = V0*Fy + V1*fy;
+                }
+                if (!isfinite(Vx)) goto skip;
+                if (!isfinite(Vy)) goto skip;
+
+                // 2D residual variations are set for the true source position
+                pixelValue += Io*(Vo + XoSave*Vx + XoSave*Vy);
             }
 
+        skip:
             // add or subtract the value
             if (add) {
@@ -336 +335 @@
 
     // restore original values
+    params->data.F32[PM_PAR_XPOS] = XoSave;
+    params->data.F32[PM_PAR_YPOS] = YoSave;
     params->data.F32[PM_PAR_I0]   = IoSave;
     params->data.F32[PM_PAR_SKY]  = skySave;
-    params->data.F32[PM_PAR_XPOS] = xPos;
-    params->data.F32[PM_PAR_YPOS] = yPos;
 
     psFree(x);
-    psFree(Ro);
-    psFree(Rx);
-    psFree(Ry);
-    psFree(myRo);
-    psFree(myRx);
-    psFree(myRy);
     psTrace("psModules.objects", 3, "---- %s(true) end ----\n", __func__);
     return(true);