Changeset 14348

Timestamp:

Jul 20, 2007, 2:05:10 PM (19 years ago)

Author:

eugene

Message:

various fixes to get PSFConv model working

Location:

trunk/psphot/src

Files:

• psphotModelTest.c (modified) (9 diffs)
• psphotPSFConvModel.c (modified) (3 diffs)

trunk/psphot/src/psphotModelTest.c

@@ -15 +15 @@
     if (!psMetadataLookupBool (&status, recipe, "TEST_FIT")) return false;
 
+    psTimerStart ("modelTest");
+
     // find the currently selected readout
     pmReadout  *readout = pmFPAfileThisReadout (config->files, view, "PSPHOT.INPUT");
@@ -54 +56 @@
     if (fitModeWord && !strcasecmp (fitModeWord, "DEFAULT")) fitMode = PM_SOURCE_FIT_EXT;
 
+    // construct the source structures
+    pmSource *source = pmSourceAlloc();
+    source->peak = pmPeakAlloc (xObj, yObj, 0, 0);
+    pmSourceDefinePixels (source, readout, xObj, yObj, OUTER);
+
     // in fitMode, psf sets the model type
     if (fitMode == PM_SOURCE_FIT_PSF) {
@@ -59 +66 @@
         if (!psf) psAbort("PSF_INPUT_FILE not supplied");
         modelType = psf->type;
+        source->type = PM_SOURCE_TYPE_STAR;
     }
     if (fitMode == PM_SOURCE_FIT_EXT) {
@@ -83 +91 @@
         modelType = pmModelSetType (modelName);
         if (modelType < 0) psAbort("unknown model %s", modelName);
+        source->type = PM_SOURCE_TYPE_EXTENDED;
     }
     if (fitMode == PM_SOURCE_FIT_PSF_X_EXT) {
@@ -111 +120 @@
         modelType = pmModelSetType (modelName);
         if (modelType < 0) psAbort("unknown model %s", modelName);
-    }
-
-    // construct the source structures
-    pmSource *source = pmSourceAlloc();
-    source->peak = pmPeakAlloc (xObj, yObj, 0, 0);
-    pmSourceDefinePixels (source, readout, xObj, yObj, OUTER);
+        source->type = PM_SOURCE_TYPE_EXTENDED;
+    }
 
     // find the local sky
@@ -140 +145 @@
     // get the initial model parameter guess
     pmModel *model = pmSourceModelGuess (source, modelType);
+    source->modelEXT = model;
 
     // if any parameters are defined by the user, take those values
@@ -162 +168 @@
     // for PSF fitting, set the shape parameters based on the PSF & source position
     if (fitMode == PM_SOURCE_FIT_PSF) {
-        pmModel *modelPSF = pmModelFromPSF (model, psf);
+        source->modelPSF = pmModelFromPSF (model, psf);
         psFree (model);
-        model = modelPSF;
+        model = source->modelPSF;
         params = model->params->data.F32;
     }
@@ -212 +218 @@
 
     // subtract object, leave local sky
-    pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
+    // pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL, maskVal);
+    pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
 
     fprintf (stderr, "output parameters: \n");
@@ -223 +230 @@
     psphotSaveImage (NULL, source->maskObj, "mask.fits");
 
+    psLogMsg ("psphot", PS_LOG_INFO, "model test : %f sec\n", psTimerMark ("modelTest"));
+
     exit (0);
 }

trunk/psphot/src/psphotPSFConvModel.c

@@ -12 +12 @@
 bool psphotPSFConvModel (pmSource *source, psMetadata *recipe, psMaskType maskVal) {
     
+    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);
@@ -17 +24 @@
     // 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)
-    psKernel *psf = psphotKernelFromPSF (source, 2);
+    // XXX for the moment, hard-wire the kernel to be 9x9 (4 pix radius)
+    psKernel *psf = psphotKernelFromPSF (source, psfSize);
+
+    // 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;
+# endif
 
     // generate copy of the model
@@ -69 +101 @@
     psTrace ("psphot", 4, "niter: %d, chisq: %f", myMin->iter, myMin->value);
 
+    // renormalize output model image
+    float Io = params->data.F32[PM_PAR_I0];
+    for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
+        for (int ix = 0; ix < source->modelFlux->numCols; ix++) {
+            source->modelFlux->data.F32[iy][ix] /= Io;
+        }
+    }
+
     // save the resulting chisq, nDOF, nIter
     modelConv->chisq = myMin->value;