Index: trunk/psModules/src/objects/pmSource.c =================================================================== --- trunk/psModules/src/objects/pmSource.c (revision 12949) +++ trunk/psModules/src/objects/pmSource.c (revision 13034) @@ -6,6 +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,5 +39,7 @@ 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,4 +49,5 @@ } +// free only the pixel data associated with this source void pmSourceFreePixels(pmSource *source) { @@ -55,9 +58,13 @@ 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,5 +83,7 @@ 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,24 +123,29 @@ // 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,8 +165,12 @@ 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,5 +201,6 @@ extend |= (mySource->pixels == NULL); extend |= (mySource->weight == NULL); - extend |= (mySource->mask == NULL); + extend |= (mySource->maskObj == NULL); + extend |= (mySource->maskView == NULL); // extend = true; @@ -190,12 +209,19 @@ 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,4 +233,6 @@ // 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,6 +486,6 @@ // 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,5 +600,5 @@ 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,5 +650,5 @@ 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,12 +747,124 @@ } -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: