Index: trunk/psModules/src/objects/pmModel.c =================================================================== --- trunk/psModules/src/objects/pmModel.c (revision 20307) +++ trunk/psModules/src/objects/pmModel.c (revision 20448) @@ -6,6 +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,9 +193,10 @@ // 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,4 +209,8 @@ 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,67 +221,29 @@ 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,8 +252,8 @@ 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,29 +268,61 @@ // 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,16 +335,10 @@ // 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);