Changeset 14338 for trunk/psphot/src/psphotModelWithPSF.c

Timestamp:

Jul 20, 2007, 10:34:30 AM (19 years ago)

Author:

eugene

Message:

finished PSFConvModel code; added structure to carry convolved model data

File:

• trunk/psphot/src/psphotModelWithPSF.c (modified) (16 diffs)

trunk/psphot/src/psphotModelWithPSF.c

@@ -6 +6 @@
     psVector *params,
     psMinConstraint *constraint,
-    const pmSource *source,
+    pmSource *source,
     const psKernel *psf,
     psMinimizeLMChi2Func func)
 {
-    psTrace("psLib.math", 3, "---- begin ----\n");
+    psTrace("psphot", 3, "---- begin ----\n");
     PS_ASSERT_PTR_NON_NULL(min, false);
     PS_ASSERT_VECTOR_NON_NULL(params, false);
@@ -46 +46 @@
     psF32 dLinear = 0.0;
 
+    // generate PCM data storage structure
+    pmPCMData *pcm = pmPCMDataAlloc (params, paramMask, source);
+
     // calculate initial alpha and beta, set chisq (min->value)
-    min->value = psphotModelWithPSF_SetABX(alpha, beta, params, paramMask, source, psf, func);
+    min->value = psphotModelWithPSF_SetABX(alpha, beta, params, paramMask, pcm, source, psf, func);
     if (isnan(min->value)) {
         min->iter = min->maxIter;
@@ -53 +56 @@
     }
     // dump some useful info if trace is defined
-    if (psTraceGetLevel("psLib.math") >= 6) {
+    if (psTraceGetLevel("psphot") >= 6) {
         p_psImagePrint(psTraceGetDestination(), alpha, "alpha guess (0)");
         p_psVectorPrint(psTraceGetDestination(), beta, "beta guess (0)");
     }
-    if (psTraceGetLevel("psLib.math") >= 5) {
+    if (psTraceGetLevel("psphot") >= 5) {
         p_psVectorPrint(psTraceGetDestination(), params, "params guess (0)");
     }
@@ -63 +66 @@
     // iterate until the tolerance is reached, or give up
     while ((min->iter < min->maxIter) && ((min->lastDelta > min->tol) || !isfinite(min->lastDelta))) {
-        psTrace("psLib.math", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
-        psTrace("psLib.math", 5, "Last delta is %f.  Min->tol is %f.\n", min->lastDelta, min->tol);
+        psTrace("psphot", 5, "Iteration number %d.  (max iterations is %d).\n", min->iter, min->maxIter);
+        psTrace("psphot", 5, "Last delta is %f.  Min->tol is %f.\n", min->lastDelta, min->tol);
 
 
@@ -75 +78 @@
 
         // dump some useful info if trace is defined
-        if (psTraceGetLevel("psLib.math") >= 6) {
+        if (psTraceGetLevel("psphot") >= 6) {
             p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (1)");
             p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (1)");
         }
-        if (psTraceGetLevel("psLib.math") >= 5) {
+        if (psTraceGetLevel("psphot") >= 5) {
             p_psVectorPrint(psTraceGetDestination(), Params, "params guess (1)");
         }
 
         // calculate Chisq for new guess, update Alpha & Beta
-        Chisq = psphotModelWithPSF_SetABX(Alpha, Beta, Params, paramMask, source, psf, func);
+        Chisq = psphotModelWithPSF_SetABX(Alpha, Beta, Params, paramMask, pcm, source, psf, func);
         if (isnan(Chisq)) {
             min->iter ++;
@@ -97 +100 @@
         psF32 rho = (min->value - Chisq) / dLinear;
 
-        psTrace("psLib.math", 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value,
+        psTrace("psphot", 5, "last chisq: %f, new chisq %f, delta: %f, rho: %f\n", min->value,
                 Chisq, min->lastDelta, rho);
 
         // dump some useful info if trace is defined
-        if (psTraceGetLevel("psLib.math") >= 6) {
+        if (psTraceGetLevel("psphot") >= 6) {
             p_psImagePrint(psTraceGetDestination(), Alpha, "alpha guess (2)");
             p_psVectorPrint(psTraceGetDestination(), Beta, "beta guess (2)");
@@ -114 +117 @@
             params = psVectorCopy(params, Params, PS_TYPE_F32);
             lambda *= 0.1;
+
+            // save the new convolved model image
+            psFree (source->modelFlux);
+            source->modelFlux = pmPCMDataSaveImage(pcm);
         } else {
             lambda *= 10.0;
@@ -119 +126 @@
         min->iter++;
     }
-    psTrace("psLib.math", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
+    psTrace("psphot", 5, "chisq: %f, last delta: %f, Niter: %d\n", min->value, min->lastDelta, min->iter);
 
     // construct & return the covariance matrix (if requested)
     if (covar != NULL) {
         if (!psMinLM_GuessABP(covar, Beta, Params, alpha, beta, params, paramMask, NULL, 0.0, NULL)) {
-            psTrace ("psLib.math", 5, "failure to calculate covariance matrix\n");
+            psTrace ("psphot", 5, "failure to calculate covariance matrix\n");
         }
     }
@@ -134 +141 @@
     psFree(Beta);
     psFree(Params);
+    psFree(pcm);
 
     if (min->iter == min->maxIter) {
-        psTrace("psLib.math", 3, "---- end (false) ----\n");
+        psTrace("psphot", 3, "---- end (false) ----\n");
         return(false);
     }
 
-    psTrace("psLib.math", 3, "---- end (true) ----\n");
+    psTrace("psphot", 3, "---- end (true) ----\n");
     return(true);
 }
@@ -149 +157 @@
     const psVector *params,
     const psVector *paramMask,
+    pmPCMData *pcm,
     const pmSource *source,
     const psKernel *psf,
@@ -168 +177 @@
     }
 
+    // generate the model and derivative images for this parameter set
+
+    // storage for model derivatives
     psVector *deriv = psVectorAlloc(params->n, PS_TYPE_F32);
-
-    // generate the model and derivative images for this parameter set
-    psImage *model = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
-    psArray *dmodels = psArrayAlloc (params->n);
-    for (psS32 n = 0; n < params->n; n++) {
-      if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
-      psImage *dmodel = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
-      dmodels->data[n] = dmodel;
-    }
 
     // working vector to store local coordinate
@@ -206 +209 @@
             coord->data.F32[1] = (psF32) (i + source->pixels->row0);
 
-            model->data.F32[i][j] = func (deriv, params, coord);
+            pcm->model->data.F32[i][j] = func (deriv, params, coord);
 
             for (int n = 0; n < params->n; n++) {
               if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
-              psImage *dmodel = dmodels->data[n];
+              psImage *dmodel = pcm->dmodels->data[n];
               dmodel->data.F32[i][j] = deriv->data.F32[n];
             }
         }
     }
-    psFree (coord);
+    psFree(coord);
+    psFree(deriv);
+
+    psphotSaveImage (NULL, pcm->model, "model1.fits");
 
     // convolve model and dmodel arrays with PSF
-    psImage *modelConv = psImageConvolveDirect (model, psf);
-    psArray *dmodelsConv = psArrayAlloc (params->n);
-    for (int n = 0; n < params->n; n++) {
-      if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
-      psImage *dmodel = dmodels->data[n];
-      psImage *dmodelConv = psImageConvolveDirect (dmodel, psf);
-      dmodelsConv->data[n] = dmodelConv;
-    }
+    psImageConvolveDirect (pcm->modelConv, pcm->model, psf);
+    for (int n = 0; n < pcm->dmodels->n; n++) {
+      if (pcm->dmodels->data[n] == NULL) continue;
+      psImage *dmodel = pcm->dmodels->data[n];
+      psImage *dmodelConv = pcm->dmodelsConv->data[n];
+      psImageConvolveDirect (dmodelConv, dmodel, psf);
+    }
+
+    // XXX TEST : SAVE IMAGES
+    psphotSaveImage (NULL, pcm->model, "model.fits");
+    psphotSaveImage (NULL, pcm->modelConv, "modelConv.fits");
+    psphotSaveImage (NULL, source->pixels, "obj.fits");
+    psphotSaveImage (NULL, source->maskObj, "mask.fits");
+    psphotSaveImage (NULL, source->weight, "weight.fits");
+    exit (0);
 
     // zero alpha and beta for summing below
@@ -248 +261 @@
             }
 
-            float ymodel  = modelConv->data.F32[i][j];
+            float ymodel  = pcm->modelConv->data.F32[i][j];
             float yweight = 1.0 / source->weight->data.F32[i][j];
             float delta = ymodel - source->pixels->data.F32[i][j];
@@ -263 +276 @@
                 continue;
               }
-              psImage *dmodel = dmodelsConv->data[n1];
+              psImage *dmodel = pcm->dmodelsConv->data[n1];
               float weight = dmodel->data.F32[i][j] * yweight;
               for (psS32 n2 = 0; n2 <= n1; n2++) {
@@ -269 +282 @@
                   continue;
                 }
-                dmodel = dmodelsConv->data[n2];
+                dmodel = pcm->dmodelsConv->data[n2];
                 alpha->data.F32[n1][n2] += weight * dmodel->data.F32[i][j];
               }
@@ -294 +307 @@
     }
 
-    psFree (model);
-    psFree (dmodels);
-    psFree (modelConv);
-    psFree (dmodelsConv);
-    psFree(deriv);
-
     return(chisq);
 }
 
-
+static void pmPCMDataFree (pmPCMData *pcm) {
+
+    if (pcm == NULL) return;
+
+    psFree (pcm->model);
+    psFree (pcm->modelConv);
+    psFree (pcm->dmodels);
+    psFree (pcm->dmodelsConv);
+    return;
+}
+
+pmPCMData *pmPCMDataAlloc (
+    const psVector *params,
+    const psVector *paramMask,
+    pmSource *source) {
+
+    pmPCMData *pcm = (pmPCMData *) psAlloc(sizeof(pmPCMData));
+    psMemSetDeallocator(pcm, (psFreeFunc) pmPCMDataFree);
+
+    // Allocate storage images for raw model and derivative images
+    pcm->model = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
+    pcm->dmodels = psArrayAlloc (params->n);
+    for (psS32 n = 0; n < params->n; n++) {
+      if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
+      pcm->dmodels->data[n] = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
+    }
+
+    // Allocate storage images for convolved model and derivative images
+    pcm->modelConv = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
+    pcm->dmodelsConv = psArrayAlloc (params->n);
+    for (psS32 n = 0; n < params->n; n++) {
+      if ((paramMask != NULL) && (paramMask->data.U8[n])) { continue; }
+      pcm->dmodelsConv->data[n] = psImageAlloc (source->pixels->numCols, source->pixels->numRows, PS_TYPE_F32);
+    }
+
+    return pcm;
+}
+
+psImage *pmPCMDataSaveImage (pmPCMData *pcm) {
+
+    psImage *model = psImageCopy (NULL, pcm->model, PS_TYPE_F32);
+    return model;
+}
 
 /*