Index: trunk/psphot/src/psphotRadialApertures.c
===================================================================
--- trunk/psphot/src/psphotRadialApertures.c	(revision 32348)
+++ trunk/psphot/src/psphotRadialApertures.c	(revision 32633)
@@ -26,6 +26,12 @@
     int num = psphotFileruleCount(config, filerule);
 
+    // skip the chisq image (optionally?)
+    int chisqNum = psMetadataLookupS32 (&status, config->arguments, "PSPHOT.CHISQ.NUM");
+    if (!status) chisqNum = -1;
+
     // loop over the available readouts
     for (int i = 0; i < num; i++) {
+        if (i == chisqNum) continue; // skip chisq image
+
 	if (!psphotRadialAperturesReadout (config, view, filerule, i, recipe, entry)) {
             psError (PSPHOT_ERR_CONFIG, false, "failed on measure extended source aperture-like parameters for %s entry %d", filerule, i);
@@ -35,4 +41,11 @@
     return true;
 }
+
+// these values are used by all threads repeatedly (and are not modified)
+static psVector *aperRadii = NULL;
+static psVector *aperRadii2 = NULL;
+static float outerRadius = NAN;
+static float SN_LIM = NAN;
+static psImageMaskType maskVal = 0;
 
 // aperture-like measurements for extended sources
@@ -75,4 +88,25 @@
         nThreads = 0;
     }
+
+    // aperRadii stores the upper bounds of the annuli
+    // XXX keep the same name here as for the petrosian / elliptical apertures?
+    aperRadii = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
+    psAssert (aperRadii, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
+    psAssert (aperRadii->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
+
+    outerRadius = aperRadii->data.F32[aperRadii->n - 1];
+
+    // save the R^2 values as well for quicker comparison
+    aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
+    for (int i = 0; i < aperRadii->n; i++) {
+	aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
+    }
+
+    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
+    maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
+    assert (maskVal);
+
+    // S/N limit to perform full non-linear fits
+    SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
 
     // source analysis is done in S/N order (brightest first)
@@ -117,12 +151,9 @@
             psArrayAdd(job->args, 1, cells->data[j]); // sources
             psArrayAdd(job->args, 1, AnalysisRegion);
-            psArrayAdd(job->args, 1, recipe);
-
             PS_ARRAY_ADD_SCALAR(job->args, entry,  PS_TYPE_S32);
             PS_ARRAY_ADD_SCALAR(job->args, nEntry, PS_TYPE_S32);
-
             PS_ARRAY_ADD_SCALAR(job->args, 0, PS_TYPE_S32); // this is used as a return value for Nradial
 
-// set this to 0 to run without threading
+	    // set this to 0 to run without threading
 # if (1)	    
             if (!psThreadJobAddPending(job)) {
@@ -138,5 +169,5 @@
 	    }
 	    psScalar *scalar = NULL;
-	    scalar = job->args->data[6];
+	    scalar = job->args->data[5];
 	    Nradial += scalar->data.S32;
 	    psFree(job);
@@ -158,5 +189,5 @@
             } else {
 		psScalar *scalar = NULL;
-		scalar = job->args->data[6];
+		scalar = job->args->data[5];
 		Nradial += scalar->data.S32;
             }
@@ -166,4 +197,5 @@
     psFree (cellGroups);
     psFree(AnalysisRegion);
+    psFree (aperRadii2);
 
     psLogMsg ("psphot", PS_LOG_WARN, "radial source apertures: %f sec for %d objects\n", psTimerMark ("psphot.radial"), Nradial);
@@ -173,5 +205,4 @@
 bool psphotRadialApertures_Threaded (psThreadJob *job) {
 
-    bool status;
     int Nradial = 0;
 
@@ -180,23 +211,11 @@
     psArray *sources        = job->args->data[1];
     psRegion *region        = job->args->data[2];
-    psMetadata *recipe      = job->args->data[3];
-
-    int entry               = PS_SCALAR_VALUE(job->args->data[4],S32); // which psf-matched image are we working on? (0 == unmatched)
-    int nEntry              = PS_SCALAR_VALUE(job->args->data[5],S32); // total number of psf-matched images + 1 unmatched
-
-    // radMax stores the upper bounds of the annuli
-    // XXX keep the same name here as for the petrosian / elliptical apertures?
-    psVector *radMax = psMetadataLookupPtr (&status, recipe, "RADIAL.ANNULAR.BINS.UPPER");
-    psAssert (radMax, "annular bins (RADIAL.ANNULAR.BINS.UPPER) are not defined in the recipe");
-    psAssert (radMax->n, "no valid annular bins (RADIAL.ANNULAR.BINS.UPPER) are define");
-
-    // user-defined masks to test for good/bad pixels (build from recipe list if not yet set)
-    psImageMaskType maskVal = psMetadataLookupImageMask(&status, recipe, "MASK.PSPHOT"); // Mask value for bad pixels
-    assert (maskVal);
-
-    // S/N limit to perform full non-linear fits
-    float SN_LIM = psMetadataLookupF32 (&status, recipe, "RADIAL_APERTURES_SN_LIM");
-
-    float outerRadius = radMax->data.F32[radMax->n - 1];
+    int entry               = PS_SCALAR_VALUE(job->args->data[3],S32); // which psf-matched image are we working on? (0 == unmatched)
+    int nEntry              = PS_SCALAR_VALUE(job->args->data[4],S32); // total number of psf-matched images + 1 unmatched
+
+    // storage for the derived pixel values (these are passed into psphotRadialApertureSource)
+    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
+    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
 
     // choose the sources of interest
@@ -204,4 +223,9 @@
 
 	pmSource *source = sources->data[i];
+
+	// if we have checked the source validity on the basis of the object set, then 
+	// we either skip these tests below or we skip the source completely
+	if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_SKIP) continue;
+	if (source->tmpFlags & PM_SOURCE_TMPF_RADIAL_KEEP) goto keepSource;
 
 	// skip PSF-like and non-astronomical objects
@@ -222,4 +246,6 @@
 	if (source->peak->x > region->x1) continue;
 	if (source->peak->y > region->y1) continue;
+
+    keepSource:
 
 	// allocate pmSourceExtendedParameters, if not already defined
@@ -240,17 +266,7 @@
 	}
 
-	// we need to change the view for the radial aperture analysis, but we want to recover exactly
-	// the original view; the following elements get destroyed by pmSourceRedefinePixels so save them:
-	psImage *oldMaskObj   = psMemIncrRefCounter(source->maskObj);
-	psImage *oldModelFlux = psMemIncrRefCounter(source->modelFlux);
-	psImage *oldPSFimage  = psMemIncrRefCounter(source->psfImage);
-	psRegion oldRegion    = source->region;
-
 	Nradial ++;
 
-	// force source image to be a bit larger...
-	pmSourceRedefinePixels (source, readout, source->peak->xf, source->peak->yf, outerRadius + 2);
-
-	if (!psphotRadialApertureSource (source, recipe, maskVal, radMax, entry)) {
+	if (!psphotRadialApertureSource (source, readout, entry, pixRadius2, pixFlux, pixVar)) {
 	    psTrace ("psphot", 5, "failed to extract radial profile for source at %7.1f, %7.1f", source->moments->Mx, source->moments->My);
 	} else {
@@ -258,20 +274,19 @@
 	}
 
-	pmSourceRedefinePixelsByRegion (source, readout, oldRegion);
-	psFree(source->maskObj);   source->maskObj   = oldMaskObj;
-	psFree(source->modelFlux); source->modelFlux = oldModelFlux;
-	psFree(source->psfImage);  source->psfImage  = oldPSFimage;
-	
 	// re-subtract the object, leave local sky
 	pmSourceSub (source, PM_MODEL_OP_FULL, maskVal);
     }
-    psScalar *scalar = job->args->data[6];
+    psScalar *scalar = job->args->data[5];
     scalar->data.S32 = Nradial;
+
+    psFree (pixRadius2);
+    psFree (pixFlux);
+    psFree (pixVar);
 
     return true;
 }
 
-bool psphotRadialApertureSource (pmSource *source, psMetadata *recipe, psImageMaskType maskVal, const psVector *aperRadii, int entry) {
-
+bool psphotRadialApertureSource (pmSource *source, pmReadout *readout, int entry, psVector *pixRadius2, psVector *pixFlux, psVector *pixVar) {
+					    
     // if we are a child source, save the results to the parent source radial aperture array
     psArray *radialAperSet = source->radialAper;
@@ -285,59 +300,68 @@
     radialAperSet->data[entry] = radialAper;
 
-    // storage for the derived pixel values
-    psVector *pixRadius2 = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *pixFlux    = psVectorAllocEmpty(100, PS_TYPE_F32);
-    psVector *pixVar     = psVectorAllocEmpty(100, PS_TYPE_F32);
+    // find the largest aperture of interest (use only apertures with inner radii <=
+    // source->skyRadius)
+    int lastAp = aperRadii->n;
+    for (int i = 0; i < aperRadii->n; i++) {
+	if (aperRadii->data.F32[i] < source->skyRadius) continue;
+	lastAp = i + 1;
+	break;
+    }
 
     // outer-most radius for initial truncation
-    float Rmax  = aperRadii->data.F32[aperRadii->n - 1];
+    float Rmax  = aperRadii->data.F32[lastAp - 1];
     float Rmax2 = PS_SQR(Rmax);
 
-    // store the R^2 values for the apertures
-    psVector *aperRadii2 = psVectorAlloc(aperRadii->n, PS_TYPE_F32);
-    for (int i = 0; i < aperRadii->n; i++) {
-	aperRadii2->data.F32[i] = PS_SQR(aperRadii->data.F32[i]);
-    }
+    // in this function, the operatins are relative to the full image (readout->image, etc)
 
     float xCM = NAN, yCM = NAN;
     if (pmSourcePositionUseMoments(source)) {
-	xCM = source->moments->Mx - 0.5 - source->pixels->col0; // coord of peak in subimage
-	yCM = source->moments->My - 0.5 - source->pixels->row0; // coord of peak in subimage
+	xCM = source->moments->Mx; // index coord of peak in readout
+	yCM = source->moments->My; // index coord of peak in readout
     } else {
-	xCM = source->peak->xf - 0.5 - source->pixels->col0; // coord of peak in subimage
-	yCM = source->peak->yf - 0.5 - source->pixels->row0; // coord of peak in subimage
-    }
+	xCM = source->peak->xf; // index coord of peak in readout
+	yCM = source->peak->yf; // index coord of peak in readout
+    }
+
+    int Nx = readout->image->numCols;
+    int Ny = readout->image->numRows;
+
+    pixRadius2->n = 0;
+    pixFlux->n = 0;
+    pixVar->n = 0;
 
     // one pass through the pixels to select the valid pixels and calculate R^2
-    for (int iy = 0; iy < source->pixels->numRows; iy++) {
-
-	float yDiff = iy - yCM;
-	if (fabs(yDiff) > Rmax) continue;
-
-	float *vPix = source->pixels->data.F32[iy];
-	float *vWgt = source->variance->data.F32[iy];
-	psImageMaskType  *vMsk = (source->maskObj == NULL) ? NULL : source->maskObj->data.PS_TYPE_IMAGE_MASK_DATA[iy];
-
-	for (int ix = 0; ix < source->pixels->numCols; ix++, vPix++, vWgt++) {
-
-	    if (vMsk) {
-		if (*vMsk & maskVal) {
-		    vMsk++;
-		    continue;
-		}
-		vMsk++;
-	    }
-	    if (isnan(*vPix)) continue;
-
-	    float xDiff = ix - xCM;
-	    if (fabs(xDiff) > Rmax) continue;
+    for (int iy = -Rmax; iy < Rmax + 1; iy++) {
+
+	float yDiff = iy + 0.5 + yCM;  // y-coordinate at this offse
+	int yPix = (int) yDiff;
+
+	if (yPix < 0) continue;
+	if (yPix > Ny - 1) continue;
+	if (fabs(iy) > Rmax) continue;
+
+	float *vPix = readout->image->data.F32[yPix];
+	float *vWgt = readout->variance->data.F32[yPix];
+	psImageMaskType  *vMsk = readout->mask->data.PS_TYPE_IMAGE_MASK_DATA[yPix];
+
+	for (int ix = -Rmax; ix < Rmax + 1; ix++) {
+
+	    float xDiff = ix + 0.5 + xCM;  // x-coordinate at this offse
+	    int xPix = (int) xDiff;
+	    
+	    if (xPix < 0) continue;
+	    if (xPix > Nx - 1) continue;
+	    if (fabs(ix) > Rmax) continue;
+	    
+	    if (vMsk[xPix] & maskVal) continue;
+	    if (isnan(vPix[xPix])) continue;
 
 	    // radius is just a function of (xDiff, yDiff)
-	    float r2  = PS_SQR(xDiff) + PS_SQR(yDiff);
+	    float r2  = PS_SQR(ix) + PS_SQR(iy);
 	    if (r2 > Rmax2) continue;
 
 	    psVectorAppend(pixRadius2, r2);
-	    psVectorAppend(pixFlux, *vPix);
-	    psVectorAppend(pixVar, *vWgt);
+	    psVectorAppend(pixFlux, vPix[xPix]);
+	    psVectorAppend(pixVar, vWgt[xPix]);
 	}
     }
@@ -348,8 +372,15 @@
     psVector *fill    = psVectorAlloc(aperRadii->n, PS_TYPE_F32); // surface brightness of radial bin
 
+    // init the apertures of interest to 0.0, the rest go to NAN
     psVectorInit (flux,    0.0);
     psVectorInit (fluxStd, 0.0);
     psVectorInit (fluxErr, 0.0);
     psVectorInit (fill,    0.0);
+    for (int i = lastAp; i < flux->n; i++) {
+	flux->data.F32[i] = NAN;
+	fluxStd->data.F32[i] = NAN;
+	fluxErr->data.F32[i] = NAN;
+	fill->data.F32[i] = NAN;
+    }
 
     float *rPix2 = pixRadius2->data.F32;
@@ -358,6 +389,8 @@
 	int j = 0;
 	float *aRad2 = aperRadii2->data.F32;
-	for (; (*aRad2 < *rPix2) && (j < aperRadii2->n); j++, aRad2++);
-	for (; j < aperRadii2->n; j++, aRad2++) {
+	for (; (*aRad2 < *rPix2) && (j < lastAp); j++, aRad2++);
+
+	// XXX I can speed this up by only saving this single aperture
+	for (; j < lastAp; j++, aRad2++) {
 	    flux->data.F32[j]    += pixFlux->data.F32[i];
 	    fluxStd->data.F32[j] += PS_SQR(pixFlux->data.F32[i]);
@@ -371,7 +404,7 @@
        2) the fractional fill factor (count of valid pixels / effective area of the aperture
        3) the error on the flux within that aperture
-     */
-
-    for (int i = 0; i < flux->n; i++) {
+    */
+
+    for (int i = 0; i < lastAp; i++) {
 	// calculate the total flux for bin 'nOut'
 	float Area = M_PI*aperRadii2->data.F32[i];
@@ -383,5 +416,5 @@
 	// XXX report the total flux or the mask-corrected flux?
 	// flux->data.F32[i]    = SBmean * Area;
-	// fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / nPix;
+	// fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]) * Area / Pinx;
 
 	fluxErr->data.F32[i] = sqrt(fluxErr->data.F32[i]);
@@ -393,16 +426,21 @@
     }
     
+# if (1)
     radialAper->flux = flux;
     radialAper->fluxStdev = fluxStd;
     radialAper->fluxErr = fluxErr;
     radialAper->fill = fill;
-
-    psFree (aperRadii2);
-    psFree (pixRadius2);
-    psFree (pixFlux);
-    psFree (pixVar);
+# else
+    // XXX TEST
+    psFree(flux);
+    psFree(fluxStd);
+    psFree(fluxErr);
+    psFree(fill);
+# endif
 
     return true;
 }
+
+/*** below is a test to use a sort to speed this up, not very successfully ***/
 
 static int nCalls = 0;
@@ -495,18 +533,18 @@
 // *** pmSourceRadialProfileSortPair is a utility function for sorting a pair of vectors
 # define COMPARE_VECT(A,B) (radius->data.F32[A] < radius->data.F32[B])
-# define SWAP_VECT(TYPE,A,B) { \
-  float tmp; \
-  if (A != B) { \
-    tmp = radius->data.F32[A]; \
-    radius->data.F32[A] = radius->data.F32[B];	\
-    radius->data.F32[B] = tmp; \
-    tmp = pixFlux->data.F32[A]; \
-    pixFlux->data.F32[A] = pixFlux->data.F32[B]; \
-    pixFlux->data.F32[B] = tmp; \
-    tmp = pixVar->data.F32[A]; \
-    pixVar->data.F32[A] = pixVar->data.F32[B]; \
-    pixVar->data.F32[B] = tmp; \
-  } \
-}
+# define SWAP_VECT(TYPE,A,B) {					\
+	float tmp;						\
+	if (A != B) {						\
+	    tmp = radius->data.F32[A];				\
+	    radius->data.F32[A] = radius->data.F32[B];		\
+	    radius->data.F32[B] = tmp;				\
+	    tmp = pixFlux->data.F32[A];				\
+	    pixFlux->data.F32[A] = pixFlux->data.F32[B];	\
+	    pixFlux->data.F32[B] = tmp;				\
+	    tmp = pixVar->data.F32[A];				\
+	    pixVar->data.F32[A] = pixVar->data.F32[B];		\
+	    pixVar->data.F32[B] = tmp;				\
+	}							\
+    }
 
 bool psphotRadialAperturesSortFlux (psVector *radius, psVector *pixFlux, psVector *pixVar) {