Changeset 13034 for trunk/psModules/src/objects/pmSource.c

Timestamp:

Apr 25, 2007, 3:20:29 PM (19 years ago)

Author:

magnier

Message:

incorporating updates from eam_02_branch (cached models, pmPSF FITS IO)

File:

• trunk/psModules/src/objects/pmSource.c (modified) (14 diffs)

trunk/psModules/src/objects/pmSource.c

@@ -6 +6 @@
  *  @author EAM, IfA: significant modifications.
  *
- *  @version $Revision: 1.26 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2007-04-21 19:47:14 $
+ *  @version $Revision: 1.27 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2007-04-26 01:20:29 $
  *
  *  Copyright 2004 Maui High Performance Computing Center, University of Hawaii
@@ -39 +39 @@
     psFree(tmp->pixels);
     psFree(tmp->weight);
-    psFree(tmp->mask);
+    psFree(tmp->maskObj);
+    psFree(tmp->maskView);
+    psFree(tmp->modelFlux);
     psFree(tmp->moments);
     psFree(tmp->modelPSF);
@@ -47 +49 @@
 }
 
+// free only the pixel data associated with this source
 void pmSourceFreePixels(pmSource *source)
 {
@@ -55 +58 @@
     psFree (source->pixels);
     psFree (source->weight);
-    psFree (source->mask);
+    psFree (source->maskObj);
+    psFree (source->maskView);
+    psFree (source->modelFlux);
 
     source->pixels = NULL;
     source->weight = NULL;
-    source->mask = NULL;
+    source->maskObj = NULL;
+    source->maskView = NULL;
+    source->modelFlux = NULL;
     return;
 }
@@ -76 +83 @@
     source->pixels = NULL;
     source->weight = NULL;
-    source->mask = NULL;
+    source->maskObj = NULL;
+    source->maskView = NULL;
+    source->modelFlux = NULL;
     source->moments = NULL;
     source->blends = NULL;
@@ -114 +123 @@
     // the original values.
     pmSource *source = pmSourceAlloc ();
-    source->peak = psMemIncrRefCounter (in->peak);
-    source->pixels = psMemIncrRefCounter (in->pixels);
-    source->weight = psMemIncrRefCounter (in->weight);
-    source->mask = psMemIncrRefCounter (in->mask);
-    source->moments = psMemIncrRefCounter (in->moments);
-    source->blends = NULL;
-    source->modelPSF = NULL;
-    source->modelEXT = NULL;
+
+    // this is actually the same peak as the original, is this the correct way to handle this?
+    source->peak = pmPeakAlloc (in->peak->x, in->peak->y, in->peak->value, in->peak->type);
+    source->peak->xf = in->peak->xf;
+    source->peak->yf = in->peak->yf;
+    source->peak->flux = in->peak->flux;
+    source->peak->SN = in->peak->SN;
+
+    // copy the values in the moments structure
+    source->moments  =  pmMomentsAlloc();
+    *source->moments = *in->moments;
+
+    // These images are all views to the parent.  
+    // We want a new view, but pointing at the same pixels.
+    source->pixels   = psImageCopyView (NULL, in->pixels);
+    source->weight   = psImageCopyView (NULL, in->weight);
+    source->maskView = psImageCopyView (NULL, in->maskView);
+
+    // the maskObj is a unique mask array; create a new mask image
+    if (in->maskObj) {
+        source->maskObj = psImageCopy (NULL, in->maskObj, PS_TYPE_MASK);
+    }
+
     source->type = in->type;
     source->mode = in->mode;
-    psMemSetDeallocator(source, (psFreeFunc) sourceFree);
-
-    // default values are NAN
-    source->psfMag = NAN;
-    source->extMag = NAN;
-    source->errMag = NAN;
-    source->apMag  = NAN;
-    source->sky    = NAN;
-    source->skyErr = NAN;
-    source->pixWeight = NAN;
 
     return(source);
@@ -151 +165 @@
     srcRegion = psRegionForImage (readout->image, srcRegion);
 
-    mySource->pixels = psMemIncrRefCounter(psImageSubset(readout->image, srcRegion));
-    mySource->weight = psMemIncrRefCounter(psImageSubset(readout->weight, srcRegion));
-    mySource->mask   = psMemIncrRefCounter(psImageSubset(readout->mask,  srcRegion));
-    mySource->region = srcRegion;
+    // these images are subset images of the equivalent parents
+    mySource->pixels   = psMemIncrRefCounter(psImageSubset(readout->image, srcRegion));
+    mySource->weight   = psMemIncrRefCounter(psImageSubset(readout->weight, srcRegion));
+    mySource->maskView = psMemIncrRefCounter(psImageSubset(readout->mask,  srcRegion));
+    mySource->region   = srcRegion;
+
+    // the object mask is a copy, and used to define the source pixels
+    mySource->maskObj = psImageCopy (NULL, mySource->maskView, PS_TYPE_MASK);
 
     return true;
@@ -183 +201 @@
     extend |= (mySource->pixels == NULL);
     extend |= (mySource->weight == NULL);
-    extend |= (mySource->mask == NULL);
+    extend |= (mySource->maskObj == NULL);
+    extend |= (mySource->maskView == NULL);
 
     // extend = true;
@@ -190 +209 @@
         psFree (mySource->pixels);
         psFree (mySource->weight);
-        psFree (mySource->mask);
-
-        mySource->pixels = psMemIncrRefCounter(psImageSubset(readout->image,  newRegion));
-        mySource->weight = psMemIncrRefCounter(psImageSubset(readout->weight, newRegion));
-        mySource->mask   = psMemIncrRefCounter(psImageSubset(readout->mask,   newRegion));
-        mySource->region = newRegion;
-    }
-
+        psFree (mySource->maskView);
+
+        mySource->pixels   = psMemIncrRefCounter(psImageSubset(readout->image,  newRegion));
+        mySource->weight   = psMemIncrRefCounter(psImageSubset(readout->weight, newRegion));
+        mySource->maskView = psMemIncrRefCounter(psImageSubset(readout->mask,   newRegion));
+        mySource->region   = newRegion;
+
+        // re-copy the main mask pixels.  NOTE: the user will need to reset the object mask
+        // pixels (eg, with psImageKeepCircle)
+        mySource->maskObj = psImageCopy (mySource->maskObj, mySource->maskView, PS_TYPE_MASK);
+        
+        // drop the old modelFlux pixels and force the user to re-create
+        psFree (mySource->modelFlux);
+        mySource->modelFlux = NULL;
+    }
     return extend;
 }
@@ -207 +233 @@
 
 // XXX EAM include a S/N cutoff in selecting the sources?
+// XXX this function should probably accept the values, not a recipe. wrap with a 
+// psphot-specific function which applies the recipe values
 pmPSFClump pmSourcePSFClump(psArray *sources, psMetadata *recipe)
 {
@@ -458 +486 @@
         // XXX EAM : can we use the value of SATURATE if mask is NULL?
         inner = psRegionForSquare (source->peak->x, source->peak->y, 2);
-        inner = psRegionForImage (source->mask, inner);
-        int Nsatpix = psImageCountPixelMask (source->mask, inner, PM_MASK_SAT);
+        inner = psRegionForImage (source->maskView, inner);
+        int Nsatpix = psImageCountPixelMask (source->maskView, inner, PM_MASK_SAT);
 
         // saturated star (size consistent with PSF or larger)
@@ -572 +600 @@
     PS_ASSERT_PTR_NON_NULL(source->peak, false);
     PS_ASSERT_PTR_NON_NULL(source->pixels, false);
-    PS_ASSERT_PTR_NON_NULL(source->mask, false);
+    PS_ASSERT_PTR_NON_NULL(source->maskObj, false);
     PS_ASSERT_FLOAT_LARGER_THAN(radius, 0.0, false);
 
@@ -622 +650 @@
         psF32 *vPix = source->pixels->data.F32[row];
         psF32 *vWgt = source->weight->data.F32[row];
-        psU8  *vMsk = (source->mask == NULL) ? NULL : source->mask->data.U8[row];
+        psU8  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.U8[row];
 
         for (psS32 col = 0; col < source->pixels->numCols ; col++, vPix++, vWgt++) {
@@ -719 +747 @@
 }
 
-pmModel *pmSourceSelectModel (pmSource *source)
-{
+// construct a realization of the source model 
+bool pmSourceCacheModel (pmSource *source) {
+
+    // select appropriate model
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+        
+    // if we already have a cached image, re-use that memory
+    source->modelFlux = psImageCopy (source->modelFlux, source->pixels, PS_TYPE_F32);
+    psImageInit (source->modelFlux, 0.0);
+
+    // in some places (psphotEnsemble), we need a normalized version
+    // in others, we just want the model.  which is more commonly used?
+    // modelFlux always has unity normalization (I0 = 1.0)
+    pmModelAdd (source->modelFlux, source->maskObj, model, PM_MODEL_OP_FULL | PM_MODEL_OP_NORM);
+    return true;
+}
+
+// should we call pmSourceCacheModel if it does not exist?
+bool pmSourceOp (pmSource *source, pmModelOpMode mode, bool add) {
+
+    bool status;
+
+    if (add) {
+        psTrace ("psphot", 3, "replacing object at %f,%f\n", source->peak->xf, source->peak->yf);
+    } else {
+        psTrace ("psphot", 3, "removing object at %f,%f\n", source->peak->xf, source->peak->yf);
+    }
+
+    pmModel *model = pmSourceGetModel (NULL, source);
+    if (model == NULL) return false;  // model must be defined
+
+    if (source->modelFlux) {
+        // add in the pixels from the modelFlux image
+        int dX = source->modelFlux->col0 - source->pixels->col0;
+        int dY = source->modelFlux->row0 - source->pixels->row0;
+        assert (dX >= 0);
+        assert (dY >= 0);
+        assert (dX + source->modelFlux->numCols <= source->pixels->numCols);
+        assert (dY + source->modelFlux->numRows <= source->pixels->numRows);
+        
+        // modelFlux has unity normalization
+        float Io = model->params->data.F32[PM_PAR_I0];
+        if (mode & PM_MODEL_OP_NORM) {
+            Io = 1.0;
+        }
+
+        psU8 **mask = NULL;
+        if (source->maskObj) {
+            mask = source->maskObj->data.U8;
+        }
+
+        // XXX need to respect the source and model masks?
+        for (int iy = 0; iy < source->modelFlux->numRows; iy++) {
+            int oy = iy + dY;
+            for (int ix = 0; ix < source->modelFlux->numCols; ix++) {
+                int ox = ix + dX;
+                if (mask && mask[iy][ix]) continue;
+                float value = Io*source->modelFlux->data.F32[iy][ix];
+                if (add) {
+                    source->pixels->data.F32[oy][ox] += value;
+                } else {
+                    source->pixels->data.F32[oy][ox] -= value;
+                }
+            }
+        }
+        return true;
+    } 
+
+    if (add) {
+        status = pmModelAdd (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
+    } else {
+        status = pmModelSub (source->pixels, source->maskObj, model, PM_MODEL_OP_FULL);
+    }
+
+    return true;
+}
+
+bool pmSourceAdd (pmSource *source, pmModelOpMode mode) {
+    bool status = pmSourceOp (source, mode, true);
+    return status;
+}
+
+bool pmSourceSub (pmSource *source, pmModelOpMode mode) {
+    bool status = pmSourceOp (source, mode, false);
+    return status;
+}
+
+// given a source, which model is currently appropriate?
+// choose PSF or EXT based on source->type, but fall back on PSF
+// if the EXT model is NULL
+pmModel *pmSourceGetModel (bool *isPSF, const pmSource *source)
+{
+
+    pmModel *model;
+
+    if (isPSF) {
+      *isPSF = false;
+    }
+
     switch (source->type) {
     case PM_SOURCE_TYPE_STAR:
-        return source->modelPSF;
+        model = source->modelPSF;
+        if (model == NULL)
+            return NULL;
+        if (isPSF) {
+            *isPSF = true;
+        }
+        return model;
 
     case PM_SOURCE_TYPE_EXTENDED:
-        return source->modelEXT;
+        model = source->modelEXT;
+        if (!model && source->modelPSF) {
+            if (isPSF) {
+                *isPSF = true;
+            }
+            return source->modelPSF;
+        }
+        return (model);
+        break;
 
     default: