Changeset 28692 for branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c

Timestamp:

Jul 20, 2010, 7:21:53 PM (16 years ago)

Author:

eugene

Message:

update PCM fitting (not quite ready)

File:

• branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c (modified) (1 diff)

branches/eam_branches/ipp-20100621/psModules/src/objects/pmPCMdata.c

@@ -1 @@
-static void pmPCMDataFree (pmPCMData *pcm) {
+/* @file  pmPCMdata.c
+ * structures and functions to support PSF-convolved model fitting
+ *
+ * @author EAM, IfA
+ *
+ * @version $Revision: 1.29 $
+ * @date $Date: 2009-02-16 22:30:50 $
+ * Copyright 2010 Institute for Astronomy, University of Hawaii
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include <stdio.h>
+#include <math.h>
+#include <string.h>
+#include <strings.h>
+#include <pslib.h>
+#include "pmHDU.h"
+#include "pmFPA.h"
+#include "pmFPAMaskWeight.h"
+
+#include "pmTrend2D.h"
+#include "pmResiduals.h"
+#include "pmGrowthCurve.h"
+#include "pmSpan.h"
+#include "pmFootprintSpans.h"
+#include "pmFootprint.h"
+#include "pmPeaks.h"
+#include "pmMoments.h"
+#include "pmModelFuncs.h"
+#include "pmModel.h"
+#include "pmModelUtils.h"
+#include "pmModelClass.h"
+#include "pmSourceMasks.h"
+#include "pmSourceExtendedPars.h"
+#include "pmSourceDiffStats.h"
+#include "pmSource.h"
+#include "pmSourceFitModel.h"
+#include "pmPCMdata.h"
+
+static void pmPCMdataFree (pmPCMdata *pcm) {
 
     if (pcm == NULL) return;
 
-    psFree (pcm->model);
+    psFree (pcm->modelFlux);
+    psFree (pcm->modelConvFlux);
+    psFree (pcm->dmodelsFlux);
+    psFree (pcm->dmodelsConvFlux);
+
     psFree (pcm->modelConv);
-    psFree (pcm->dmodels);
-    psFree (pcm->dmodelsConv);
+    psFree (pcm->psf);
+    psFree (pcm->constraint);
     return;
 }
 
-pmPCMData *pmPCMDataAlloc (
+pmPCMdata *pmPCMdataAlloc (
     const psVector *params,
     const psVector *paramMask,
     pmSource *source) {
 
-    pmPCMData *pcm = (pmPCMData *) psAlloc(sizeof(pmPCMData));
-    psMemSetDeallocator(pcm, (psFreeFunc) pmPCMDataFree);
+    pmPCMdata *pcm = (pmPCMdata *) psAlloc(sizeof(pmPCMdata));
+    psMemSetDeallocator(pcm, (psFreeFunc) pmPCMdataFree);
 
     // Allocate storage images for raw model and derivative images
-    pcm->model = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
-    pcm->dmodels = psArrayAlloc (params->n);
+    pcm->modelFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
+    pcm->dmodelsFlux = psArrayAlloc (params->n);
     for (psS32 n = 0; n < params->n; n++) {
-        pcm->dmodels->data[n] = NULL;
+        pcm->dmodelsFlux->data[n] = NULL;
         if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
-        pcm->dmodels->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
+        pcm->dmodelsFlux->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
     }
 
     // Allocate storage images for convolved model and derivative images
-    pcm->modelConv = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
-    pcm->dmodelsConv = psArrayAlloc (params->n);
+    pcm->modelConvFlux = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
+    pcm->dmodelsConvFlux = psArrayAlloc (params->n);
     for (psS32 n = 0; n < params->n; n++) {
-        pcm->dmodelsConv->data[n] = NULL;
+        pcm->dmodelsConvFlux->data[n] = NULL;
         if ((paramMask != NULL) && (paramMask->data.PS_TYPE_VECTOR_MASK_DATA[n])) { continue; }
-        pcm->dmodelsConv->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
-    }
+        pcm->dmodelsConvFlux->data[n] = psImageCopy (NULL, source->pixels, PS_TYPE_F32);
+    }
+
+    pcm->modelConv = NULL;
+    pcm->psf = NULL;
+    pcm->constraint = NULL;
+    pcm->nDOF = 0;
 
     return pcm;
 }
 
-psImage *pmPCMDataSaveImage (pmPCMData *pcm) {
-
-    psImage *model = psImageCopy (NULL, pcm->modelConv, PS_TYPE_F32);
+psImage *pmPCMdataSaveImage (pmPCMdata *pcm) {
+
+    psImage *model = psImageCopy (NULL, pcm->modelConvFlux, PS_TYPE_F32);
 
     return model;
 }
 
+psKernel *pmPCMkernelFromPSF (pmSource *source, int nPix) {
+
+    assert (source);
+    assert (source->psfImage); // XXX build if needed?
+
+    int x0 = source->peak->xf - source->psfImage->col0;
+    int y0 = source->peak->yf - source->psfImage->row0;
+
+    // need to decide on the size: dynamically? statically?
+    psKernel *psf = psKernelAlloc (-nPix, +nPix, -nPix, +nPix);
+
+    // XXX we should just re-construct a PSF at this location 
+    // psModelAdd (psf->image, NULL, source->modelPSF, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM | PM_MODEL_OP_CENTER);
+  
+    // if the realized PSF for this object does not cover the full kernel, give up for now
+    if (x0 + psf->xMin < 0) goto escape;
+    if (x0 + psf->xMax >= source->psfImage->numCols) goto escape;
+    if (y0 + psf->yMin < 0) goto escape;
+    if (y0 + psf->yMax >= source->psfImage->numRows) goto escape;
+
+    double sum = 0.0;
+    for (int j = psf->yMin; j <= psf->yMax; j++) {
+        for (int i = psf->xMin; i <= psf->xMax; i++) {
+            double value = source->psfImage->data.F32[y0 + j][x0 + i];
+            psf->kernel[j][i] = value;
+            sum += value;
+        }
+    }
+    assert (sum > 0.0);
+
+    // psf must be normalized (integral = 1.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;
+        }
+    }
+
+    return psf;
+
+escape:
+    psFree (psf);
+    return NULL;
+}
+
+pmPCMdata *pmPCMinit(pmSource *source, pmSourceFitOptions *fitOptions, pmModelType modelType, psImageMaskType maskVal, float psfSize) {
+
+    // make sure we savep a cached copy of the psf flux
+    pmSourceCachePSF (source, maskVal);
+
+    // convert the cached cached psf model for this source to a psKernel
+    psKernel *psf = pmPCMkernelFromPSF (source, psfSize);
+    if (!psf) return NULL;
+
+# 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
+
+    // allocate the model
+    pmModel *modelConv = pmModelAlloc(modelType);
+    if (!modelConv) {
+        psFree (psf);
+        return NULL;
+    }
+
+    // count the number of unmasked pixels:
+    int nPix = 0;
+    for (psS32 i = 0; i < source->pixels->numRows; i++) {
+        for (psS32 j = 0; j < source->pixels->numCols; j++) {
+            // XXX are we doing the right thing with the mask?
+            // skip masked points
+            if (source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[i][j]) {
+                continue;
+            }
+            // skip zero-variance points
+            if (source->variance->data.F32[i][j] == 0) {
+                continue;
+            }
+            // skip nan value points
+            if (!isfinite(source->pixels->data.F32[i][j])) {
+                continue;
+            }
+            nPix ++;
+        }
+    }    
+
+    psVector *params  = modelConv->params;
+
+    // create the minimization constraints
+    psMinConstraint *constraint = psMinConstraintAlloc();
+    constraint->paramMask = psVectorAlloc (params->n, PS_TYPE_VECTOR_MASK);
+    constraint->checkLimits = modelConv->modelLimits;
+
+    // set parameter mask based on fitting mode
+    int nParams = 0;
+    switch (fitOptions->mode) {
+      case PM_SOURCE_FIT_NORM:
+        // NORM-only model fits only source normalization (Io)
+        nParams = 1;
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+        break;
+      case PM_SOURCE_FIT_PSF:
+        // PSF model only fits x,y,Io
+        nParams = 3;
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_XPOS] = 0;
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_YPOS] = 0;
+        break;
+      case PM_SOURCE_FIT_EXT:
+        // EXT model fits all params (except sky)
+        nParams = params->n - 1;
+        psVectorInit (constraint->paramMask, 0);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
+        break;
+      case PM_SOURCE_FIT_INDEX:
+        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
+        psVectorInit (constraint->paramMask, 1);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_I0] = 0;
+        if (params->n == 7) {
+            nParams = 1;
+        } else {
+            nParams = 2;
+            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 0;
+        }
+        break;
+      case PM_SOURCE_FIT_NO_INDEX:
+        // PSF model only fits Io, index (PAR7) -- only Io for models with < 8 params
+        psVectorInit (constraint->paramMask, 0);
+        constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_SKY] = 1;
+        if (params->n == 7) {
+            nParams = params->n - 1;
+        } else {
+            nParams = params->n - 2;
+            constraint->paramMask->data.PS_TYPE_VECTOR_MASK_DATA[PM_PAR_7] = 1;
+        }
+        break;
+      default:
+        psAbort("invalid fitting mode");
+    }
+
+    // generate PCM data storage structure
+    pmPCMdata *pcm = pmPCMdataAlloc (params, constraint->paramMask, source);
+
+    pcm->modelConv = modelConv;
+    pcm->psf = psf;
+    pcm->constraint = constraint;
+
+    if (nPix <  nParams + 1) {
+        psTrace ("psModules.objects", 4, "insufficient valid pixels\n");
+        pcm->modelConv->flags |= PM_MODEL_STATUS_BADARGS;
+        abort ();
+        // XXX This should not be an abort!!
+    }
+    pcm->nPix = nPix;
+    pcm->nDOF = nPix - nParams - 1;
+
+    return pcm;
+}