Index: trunk/psModules/src/objects/pmGrowthCurveGenerate.c
===================================================================
--- trunk/psModules/src/objects/pmGrowthCurveGenerate.c	(revision 20076)
+++ trunk/psModules/src/objects/pmGrowthCurveGenerate.c	(revision 20233)
@@ -5,6 +5,6 @@
  *  @author EAM, IfA
  *
- *  @version $Revision: 1.1 $ $Name: not supported by cvs2svn $
- *  @date $Date: 2008-10-13 01:56:42 $
+ *  @version $Revision: 1.2 $ $Name: not supported by cvs2svn $
+ *  @date $Date: 2008-10-17 22:58:41 $
  *
  *  Copyright 2004 Institute for Astronomy, University of Hawaii
@@ -39,4 +39,6 @@
 #include "pmErrorCodes.h"
 
+pmGrowthCurve *pmGrowthCurveForPosition (psImage *image, pmPSF *psf, bool ignore, psMaskType maskVal, psMaskType markVal, float xc, float yc);
+
 /*****************************************************************************/
 /* FUNCTION IMPLEMENTATION - PUBLIC                                          */
@@ -49,46 +51,103 @@
     PS_ASSERT_PTR_NON_NULL(readout->image, false);
 
-    // bool status;
-    float xc, yc, dx, dy;
+    // maskVal is used to test for rejected pixels, and must include markVal
+    maskVal |= markVal;
+
+    // XXX something of a hack: measure the growth curve at a number of points in the field and
+    // average them together
+
+    psArray *growths = psArrayAllocEmpty (100);
+
+    for (float ix = -0.4; ix <= +0.4; ix += 0.2) {
+	for (float iy = -0.4; iy <= +0.4; iy += 0.2) {
+
+	    // use the center of the center pixel of the image
+	    float xc = ix*readout->image->numCols + 0.5*readout->image->numCols + readout->image->col0 + 0.5;
+	    float yc = iy*readout->image->numRows + 0.5*readout->image->numRows + readout->image->row0 + 0.5;
+
+	    pmGrowthCurve *growth = pmGrowthCurveForPosition (readout->image, psf, ignore, maskVal, markVal, xc, yc);
+	    if (!growth) continue;
+
+	    psArrayAdd (growths, 100, growth);
+	    psFree (growth);
+	}
+    }
+    psAssert (growths->n, "cannot build growth curve (psf model is invalid everywhere)");
+
+    // just use a simple sample median to get the 'best' value from each growth curve...
+    psStats *stats = psStatsAlloc (PS_STAT_SAMPLE_MEDIAN);
+
+    psVector *values = psVectorAlloc (growths->n, PS_DATA_F32);
+
+    // median the values for the fitMags
+    values->n = 0;
+    for (int j = 0; j < growths->n; j++) {
+	pmGrowthCurve *growth = growths->data[j];
+	if (!isfinite(growth->fitMag)) continue;
+	psVectorAppend (values, growth->fitMag);
+    }
+    psVectorStats (stats, values, NULL, NULL, 0);
+    psf->growth->fitMag = stats->sampleMedian;
+
+    // loop over a range of source fluxes
+    // no need to interpolate since we have forced the object center
+    // to 0.5, 0.5 above
+    for (int i = 0; i < psf->growth->radius->n; i++) {
+
+	// median the values for each radial bin
+	values->n = 0;
+	for (int j = 0; j < growths->n; j++) {
+	    pmGrowthCurve *growth = growths->data[j];
+	    if (!isfinite(growth->apMag->data.F32[i])) continue;
+	    psVectorAppend (values, growth->apMag->data.F32[i]);
+	}
+	psVectorStats (stats, values, NULL, NULL, 0);
+	psf->growth->apMag->data.F32[i] = stats->sampleMedian;
+    }
+    psf->growth->apRef = psVectorInterpolate (psf->growth->radius, psf->growth->apMag, psf->growth->refRadius);
+    psf->growth->apLoss = psf->growth->fitMag - psf->growth->apRef;
+
+    psLogMsg ("psphot.growth", 4, "GrowthCurve : apLoss : %f (fitMag - apMag @ ref : %f - %f)\n", psf->growth->apLoss, psf->growth->fitMag, psf->growth->apRef);
+
+    psFree (growths);
+    psFree (stats);
+    psFree (values);
+
+    return true;
+}
+
+pmGrowthCurve *pmGrowthCurveForPosition (psImage *image, pmPSF *psf, bool ignore, psMaskType maskVal, psMaskType markVal, float xc, float yc) {
+
     float fitMag, apMag;
     float radius;
 
+    assert (psf->growth);
+
+    float minRadius = psf->growth->radius->data.F32[0];
+    pmGrowthCurve *growth = pmGrowthCurveAlloc (minRadius, psf->growth->maxRadius, psf->growth->refRadius);
+
+    float dx = growth->maxRadius + 1;
+    float dy = growth->maxRadius + 1;
+
     // create template model
     pmModel *modelRef = pmModelAlloc(psf->type);
 
-    // maskVal is used to test for rejected pixels, and must include markVal
-    maskVal |= markVal;
-
-    // XXXX A Serious hack: the psf might not be determined in the field center.
-    // for the moment, just try a few locations until it is defined!
-
-    pmModel *model = NULL;
-    for (int ix = -1; ix <= +1; ix ++) {
-	for (int iy = -1; iy <= +1; iy ++) {
-
-	    // use the center of the center pixel of the image
-	    xc = (0.5 + 0.3*ix)*readout->image->numCols + readout->image->col0 + 0.5;
-	    yc = (0.5 + 0.3*ix)*readout->image->numRows + readout->image->row0 + 0.5;
-	    dx = psf->growth->maxRadius + 1;
-	    dy = psf->growth->maxRadius + 1;
-
-	    // assign the x and y coords to the image center
-	    // create an object with center intensity of 1000
-	    modelRef->params->data.F32[PM_PAR_SKY] = 0;
-	    modelRef->params->data.F32[PM_PAR_I0] = 1000;
-	    modelRef->params->data.F32[PM_PAR_XPOS] = xc;
-	    modelRef->params->data.F32[PM_PAR_YPOS] = yc;
-
-	    // create modelPSF from this model
-	    model = pmModelFromPSF (modelRef, psf);
-	    if (model) goto got_model;
-	}
-    }
-    psAssert (model, "cannot build growth curve (psf model is invalid everywhere)");
-
-got_model:
+    // assign the x and y coords to the image center
+    // create an object with center intensity of 1000
+    modelRef->params->data.F32[PM_PAR_SKY] = 0;
+    modelRef->params->data.F32[PM_PAR_I0] = 1000;
+    modelRef->params->data.F32[PM_PAR_XPOS] = xc;
+    modelRef->params->data.F32[PM_PAR_YPOS] = yc;
+
+    // create modelPSF from this model
+    pmModel *model = pmModelFromPSF (modelRef, psf);
+    if (!model) {
+	psFree (growth);
+	return NULL;
+    }
+
     // measure the fitMag for this model
     pmSourcePhotometryModel (&fitMag, model);
-    psf->growth->fitMag = fitMag;
+    growth->fitMag = fitMag;
 
     // generate working image for this source
@@ -96,12 +155,12 @@
 
     // force region to stop at dimensions of image
-    region = psRegionForImage (readout->image, region);
+    region = psRegionForImage (image, region);
 
     // the view, image, and mask retain col0,row0
-    psImage *view = psImageSubset (readout->image, region);
-    psImage *image = psImageCopy (NULL, view, PS_TYPE_F32);
+    psImage *view = psImageSubset (image, region);
+    psImage *pixels = psImageCopy (NULL, view, PS_TYPE_F32);
     psImage *mask = psImageCopy (NULL, view, PS_TYPE_U8);
 
-    psImageInit (image, 0.0);
+    psImageInit (pixels, 0.0);
     psImageInit (mask, 0);
 
@@ -109,22 +168,22 @@
     // no need to mask the source here
     // XXX should we measure this for the analytical model only or the full model?
-    pmModelAdd (image, NULL, model, PM_MODEL_OP_FULL, maskVal);
-
-    // loop over a range of source fluxes
-    // no need to interpolate since we have forced the object center
-    // to 0.5, 0.5 above
-    bool completeGrowthCurve = true;            // do we have a complete curve of growth?
-    for (int i = 0; i < psf->growth->radius->n; i++) {
-
-        radius = psf->growth->radius->data.F32[i];
+    pmModelAdd (pixels, NULL, model, PM_MODEL_OP_FULL, maskVal);
+
+    // Loop over a range of radii.  No need to interpolate since we have forced the object
+    // center to 0.5, 0.5 above
+    for (int i = 0; i < growth->radius->n; i++) {
+
+        radius = growth->radius->data.F32[i];
 
         // mask the given aperture and measure the apMag
         psImageKeepCircle (mask, xc, yc, radius, "OR", markVal);
-        if (!pmSourcePhotometryAper (&apMag, model, image, mask, maskVal)) {
-            psError(PM_ERR_PHOTOM, false, "Measuring apMag for radius == %g", radius);
-            psErrorStackPrint (stderr, "failure to get growth curve");
-            psErrorClear();
-            completeGrowthCurve = false;
-            break;
+        if (!pmSourcePhotometryAper (&apMag, model, pixels, mask, maskVal)) {
+	    psFree (growth);
+	    psFree (view);
+	    psFree (pixels);
+	    psFree (mask);
+	    psFree (model);
+	    psFree (modelRef);
+	    return NULL;
         }
         psImageKeepCircle (mask, xc, yc, radius, "AND", PS_NOT_U8(markVal));
@@ -132,26 +191,18 @@
         // the 'ignore' mode is for testing
         if (ignore) {
-            psf->growth->apMag->data.F32[i] = fitMag;
+            growth->apMag->data.F32[i] = fitMag;
         } else {
-            psf->growth->apMag->data.F32[i] = apMag;
+            growth->apMag->data.F32[i] = apMag;
         }
     }
-
-    if (completeGrowthCurve) {
-        psf->growth->apRef = psVectorInterpolate (psf->growth->radius, psf->growth->apMag, psf->growth->refRadius);
-        psf->growth->apLoss = psf->growth->fitMag - psf->growth->apRef;
-    } else {
-        psf->growth->apRef = NAN;
-        psf->growth->apLoss = 0;
-    }
-
-    psLogMsg ("psphot.growth", 4, "GrowthCurve : apLoss : %f\n", psf->growth->apLoss);
-
+    
     psFree (view);
-    psFree (image);
+    psFree (pixels);
     psFree (mask);
     psFree (model);
     psFree (modelRef);
 
-    return completeGrowthCurve ? true : false;
+    // psLogMsg ("psModules", 4, "GrowthCurve for %f,%f\n", xc, yc);
+
+    return growth;
 }