Changeset 17396 for trunk/psphot/src/psphotPSFConvModel.c

Timestamp:

Apr 8, 2008, 8:36:06 AM (18 years ago)

Author:

eugene

Message:

merging from eam_branch_20080324 : psphot work on extended source fitting and related I/O functions

File:

• trunk/psphot/src/psphotPSFConvModel.c (modified) (4 diffs)

trunk/psphot/src/psphotPSFConvModel.c

@@ -1 @@
-# include "psphot.h"
+# include "psphotInternal.h"
+# define USE_DELTA_PSF 0
 
 // save as static values so they may be set externally
 static psF32 PM_SOURCE_FIT_MODEL_NUM_ITERATIONS = 15;
 static psF32 PM_SOURCE_FIT_MODEL_TOLERANCE = 0.1;
-// static psF32 PM_SOURCE_FIT_MODEL_WEIGHT = 1.0;
-// static bool  PM_SOURCE_FIT_MODEL_PIX_WEIGHTS = true;
 
 // input source has both modelPSF and modelEXT.  on successful exit, we set the
 // modelConv to contain the fitted parameters, and the modelFlux to contain the 
 // convolved model image.
-bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
+pmModel *psphotPSFConvModel (pmReadout *readout, pmSource *source, pmModelType modelType, psMaskType maskVal, int psfSize) {
     
-    bool status;
-
-    int psfSize = psMetadataLookupS32 (&status, recipe, "PCM_BOX_SIZE");
-    if (!status) {
-        psfSize = 2;
-    }
-
     // make sure we save a cached copy of the psf flux
     pmSourceCachePSF (source, maskVal);
 
     // convert the cached cached psf model for this source to a psKernel
-    // XXX for the moment, hard-wire the kernel to be 5x5 (2 pix radius)
-    // XXX for the moment, hard-wire the kernel to be 9x9 (4 pix radius)
     psKernel *psf = psphotKernelFromPSF (source, psfSize);
+    if (!psf) return NULL;
 
-    // psf must be normalized (integral = 1.0)
-    double sum = 0.0;
-    for (int i = 0; i < psf->image->numRows; i++) {
-        for (int j = 0; j < psf->image->numCols; j++) {
-            sum += psf->image->data.F32[i][j];
-        }
-    }
-    assert (sum > 0.0);
-    for (int i = 0; i < psf->image->numRows; i++) {
-        for (int j = 0; j < psf->image->numCols; j++) {
-            psf->image->data.F32[i][j] /= sum;
-        }
-    }
-
-# if (0)
-    // XXX sanity check: convolve with delta function should behave like unconvolved version
-    for (int i = 0; i < psf->image->numRows; i++) {
-        for (int j = 0; j < psf->image->numCols; j++) {
-            psf->image->data.F32[i][j] = 0.0;
-        }
-    }
-    psf->image->data.F32[2][2] = 1.0;
+# if (USE_DELTA_PSF)
+    psImageInit (psf->image, 0.0);
+    psf->image->data.F32[(int)(0.5*psf->image->numRows)][(int)(0.5*psf->image->numCols)] = 1.0;
 # endif
 
@@ -54 +26 @@
     // XXX we could modify the parameter values or even the model 
     // here based on the observed seeing (some lookup table...)
-    pmModel *modelConv = pmModelCopy (source->modelEXT);
+
+    // use the source moments, etc to guess basic model parameters
+    pmModel *modelConv = pmSourceModelGuess (source, modelType);
+    if (!modelConv) {
+        psFree (psf);
+        return NULL;
+    }
+
+    // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
+    psEllipseShape psfShape;
+    psfShape.sx  = source->modelPSF->params->data.F32[PM_PAR_SXX] / M_SQRT2;
+    psfShape.sxy = source->modelPSF->params->data.F32[PM_PAR_SXY];
+    psfShape.sy  = source->modelPSF->params->data.F32[PM_PAR_SYY] / M_SQRT2;
+    psEllipseAxes psfAxes = psEllipseShapeToAxes (psfShape, 20.0);
+
+    // XXX test : modify the Io, SXX, SYY terms based on the psf SXX, SYY terms:
+    psEllipseShape extShape;
+    extShape.sx  = modelConv->params->data.F32[PM_PAR_SXX] / M_SQRT2;
+    extShape.sxy = modelConv->params->data.F32[PM_PAR_SXY];
+    extShape.sy  = modelConv->params->data.F32[PM_PAR_SYY] / M_SQRT2;
+    psEllipseAxes extAxes = psEllipseShapeToAxes (extShape, 20.0);
+
+    // decrease the initial guess ellipse by psf_minor axis:
+    psEllipseAxes extAxesMod;
+    extAxesMod.major = sqrt (PS_MAX (1.0, PS_SQR(extAxes.major) - PS_SQR(psfAxes.minor)));
+    extAxesMod.minor = sqrt (PS_MAX (1.0, PS_SQR(extAxes.minor) - PS_SQR(psfAxes.minor)));
+    extAxesMod.theta = extAxes.theta;
+
+    psEllipseShape extShapeMod = psEllipseAxesToShape (extAxesMod);
+    modelConv->params->data.F32[PM_PAR_SXX] = extShapeMod.sx * M_SQRT2;
+    modelConv->params->data.F32[PM_PAR_SXY] = extShapeMod.sxy;
+    modelConv->params->data.F32[PM_PAR_SYY] = extShapeMod.sy * M_SQRT2;
+
+    // increase the initial guess central intensity by 2pi r^2:
+    modelConv->params->data.F32[PM_PAR_I0] *= (1.0 + PS_SQR(psfAxes.minor) / PS_SQR(extAxesMod.minor));
+
     psVector *params  = modelConv->params;
     psVector *dparams = modelConv->dparams;
+
+    psphotCheckRadiusEXT (readout, source, modelConv);
 
     // create the minimization constraints
@@ -91 +100 @@
     psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
 
-    // renormalize output model image
+    // renormalize output model image (generated by fitting process)
     float Io = params->data.F32[PM_PAR_I0];
     for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
@@ -115 +124 @@
     onPic &= (params->data.F32[PM_PAR_YPOS] >= source->pixels->row0);
     onPic &= (params->data.F32[PM_PAR_YPOS] <  source->pixels->row0 + source->pixels->numRows);
-    if (!onPic) {
-        modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE;
-    }
+    if (!onPic) modelConv->flags |= PM_MODEL_STATUS_OFFIMAGE;
 
-    source->mode |= PM_SOURCE_MODE_FITTED;
-    source->modelConv = modelConv;
+    source->mode |= PM_SOURCE_MODE_FITTED; // XXX is this needed?
 
+    psFree(psf);
     psFree(myMin);
     psFree(covar);
     psFree(constraint);
 
-    bool retval = (onPic && fitStatus);
-    psTrace("psphot", 5, "---- %s(%d) end ----\n", __func__, retval);
-    return(retval);
+    return modelConv;
 }