Index: trunk/psphot/src/psphotImageMedian.c
===================================================================
--- trunk/psphot/src/psphotImageMedian.c	(revision 13013)
+++ trunk/psphot/src/psphotImageMedian.c	(revision 13411)
@@ -2,11 +2,11 @@
 static int npass = 0;
 
-// we have 4 possibilities: (INTERNAL or I/O file) and (exists or not) 
+// we have 4 possibilities: (INTERNAL or I/O file) and (exists or not)
 // select model pixels (from output background model file, or create internal file)
 static pmReadout *get_model_readout(const char *name, // name of internal/external file
-				    const pmConfig *config, // configuration information
-				    pmFPAview *view,
-				    pmFPA *inFPA,
-				    const psImageBinning *binning) {
+                                    const pmConfig *config, // configuration information
+                                    const pmFPAview *view,
+                                    pmFPA *inFPA,
+                                    const psImageBinning *binning) {
     pmReadout *model = NULL;
 
@@ -14,23 +14,23 @@
     pmFPAfile *file = psMetadataLookupPtr(&status, config->files, name);
     if (file == NULL) {
-	// we are not using PSPHOT.BACKMDL as an I/O file: define an internal version
+        // we are not using PSPHOT.BACKMDL as an I/O file: define an internal version
         model = pmFPAfileDefineInternal (config->files, name, binning->nXruff, binning->nYruff, PS_TYPE_F32);
     } else {
-	if (file->mode == PM_FPA_MODE_INTERNAL) {
-	    // we are not using PSPHOT.BACKMDL as an I/O file: already defined above
-	    model = file->readout;
-	} else {
-	    // we are using PSPHOT.BACKMDL as an I/O file: select readout or create
-	    model = pmFPAviewThisReadout (view, file->fpa);
-	    if (model == NULL) {
-		// readout does not yet exist: create from input
-		// XXX we have an inconsistency in this calculation here and in pmFPACopy
-		// XXX use the psImageBinning functions to set the output image size
-		pmFPAfileCopyStructureView (file->fpa, inFPA, binning->nXbin, binning->nYbin, view);
-		model = pmFPAviewThisReadout (view, file->fpa);
-		PS_ASSERT (binning->nXruff == model->image->numCols, false);
-		PS_ASSERT (binning->nYruff == model->image->numRows, false);
-	    }
-	}
+        if (file->mode == PM_FPA_MODE_INTERNAL) {
+            // we are not using PSPHOT.BACKMDL as an I/O file: already defined above
+            model = file->readout;
+        } else {
+            // we are using PSPHOT.BACKMDL as an I/O file: select readout or create
+            model = pmFPAviewThisReadout (view, file->fpa);
+            if (model == NULL) {
+                // readout does not yet exist: create from input
+                // XXX we have an inconsistency in this calculation here and in pmFPACopy
+                // XXX use the psImageBinning functions to set the output image size
+                pmFPAfileCopyStructureView (file->fpa, inFPA, binning->nXbin, binning->nYbin, view);
+                model = pmFPAviewThisReadout (view, file->fpa);
+                PS_ASSERT (binning->nXruff == model->image->numCols, false);
+                PS_ASSERT (binning->nYruff == model->image->numRows, false);
+            }
+        }
     }
 
@@ -40,5 +40,5 @@
 // generate the median in NxN boxes, clipping heavily
 // linear interpolation to generate full-scale model
-bool psphotImageMedian (pmConfig *config, pmFPAview *view)
+bool psphotImageMedian (pmConfig *config, const pmFPAview *view)
 {
     bool status = true;
@@ -70,37 +70,37 @@
     char *statsName = psMetadataLookupStr (&status, recipe, "SKY_STAT");
     if (statsName == NULL) {
-	statsName = defaultStatsName;
+        statsName = defaultStatsName;
     }
     psStatsOptions statsOptionLocation = psStatsOptionFromString(statsName);
     if (!(statsOptionLocation & (PS_STAT_SAMPLE_MEAN |
-				 PS_STAT_SAMPLE_MEDIAN |
-				 PS_STAT_ROBUST_MEDIAN |
-				 PS_STAT_ROBUST_QUARTILE |
-				 PS_STAT_CLIPPED_MEAN |
-				 PS_STAT_FITTED_MEAN |
-				 PS_STAT_FITTED_MEAN_V2 |
-				 PS_STAT_FITTED_MEAN_V3))) {
-	statsOptionLocation = PS_STAT_FITTED_MEAN;
+                                 PS_STAT_SAMPLE_MEDIAN |
+                                 PS_STAT_ROBUST_MEDIAN |
+                                 PS_STAT_ROBUST_QUARTILE |
+                                 PS_STAT_CLIPPED_MEAN |
+                                 PS_STAT_FITTED_MEAN |
+                                 PS_STAT_FITTED_MEAN_V2 |
+                                 PS_STAT_FITTED_MEAN_V3))) {
+        statsOptionLocation = PS_STAT_FITTED_MEAN;
     }
 
     psStatsOptions statsOptionWidth = PS_STAT_NONE;
     if (statsOptionLocation & (PS_STAT_SAMPLE_MEAN | PS_STAT_SAMPLE_MEDIAN)) {
-	statsOptionWidth = PS_STAT_SAMPLE_STDEV;
+        statsOptionWidth = PS_STAT_SAMPLE_STDEV;
     } else if (statsOptionLocation & (PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_QUARTILE)) {
 #if 1
-	statsOptionWidth = PS_STAT_ROBUST_STDEV; // not set; => NaN
+        statsOptionWidth = PS_STAT_ROBUST_STDEV; // not set; => NaN
 #else
-	statsOptionWidth = PS_STAT_FITTED_STDEV;
+        statsOptionWidth = PS_STAT_FITTED_STDEV;
 #endif
     } else if (statsOptionLocation & PS_STAT_FITTED_MEAN) {
-	statsOptionWidth = PS_STAT_FITTED_STDEV;
+        statsOptionWidth = PS_STAT_FITTED_STDEV;
     } else if (statsOptionLocation & PS_STAT_CLIPPED_MEAN) {
-	statsOptionWidth = PS_STAT_CLIPPED_STDEV;
+        statsOptionWidth = PS_STAT_CLIPPED_STDEV;
     } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V2) {
-	statsOptionWidth = PS_STAT_FITTED_STDEV_V2;
+        statsOptionWidth = PS_STAT_FITTED_STDEV_V2;
     } else if (statsOptionLocation & PS_STAT_FITTED_MEAN_V3) {
-	statsOptionWidth = PS_STAT_FITTED_STDEV_V3;
+        statsOptionWidth = PS_STAT_FITTED_STDEV_V3;
     } else {
-	psAbort("Unable to estimate variance of selected statsOptionLocations 0x%x", statsOptionLocation);
+        psAbort("Unable to estimate variance of selected statsOptionLocations 0x%x", statsOptionLocation);
     }
 
@@ -110,6 +110,6 @@
     // set range for old-version of sky statistic
     if (statsOptionLocation & PS_STAT_ROBUST_QUARTILE) {
-	stats->min = 0.25;
-	stats->max = 0.75;
+        stats->min = 0.25;
+        stats->max = 0.75;
     }
 
@@ -130,8 +130,8 @@
     if (!status) {
         if ((stats->options & PS_STAT_FITTED_MEAN) || (stats->options & PS_STAT_FITTED_MEAN_V2)) {
-	    stats->clipSigma = 1.0;
-	} else {
-	    stats->clipSigma = 3.0;
-	}
+            stats->clipSigma = 1.0;
+        } else {
+            stats->clipSigma = 3.0;
+        }
     }
 
@@ -172,52 +172,52 @@
     for (int iy = 0; iy < model->image->numRows; iy++) {
         for (int ix = 0; ix < model->image->numCols; ix++) {
-	    
-	    // convert the ruff grid cell to the equivalent fine grid cell
-	    // XXX we need to watch out for row0,col0
-	    ruffRegion = psRegionSet (ix, ix + 2, iy, iy + 2);
-	    fineRegion = psImageBinningSetFineRegion (binning, ruffRegion);
+
+            // convert the ruff grid cell to the equivalent fine grid cell
+            // XXX we need to watch out for row0,col0
+            ruffRegion = psRegionSet (ix, ix + 2, iy, iy + 2);
+            fineRegion = psImageBinningSetFineRegion (binning, ruffRegion);
             fineRegion = psRegionForImage (image, fineRegion);
 
             psImage *subset  = psImageSubset (image, fineRegion);
-	    if (!subset->numCols || !subset->numRows) {
-		psFree (subset);
-		continue;
-	    }
+            if (!subset->numCols || !subset->numRows) {
+                psFree (subset);
+                continue;
+            }
             psImage *submask = psImageSubset (mask, fineRegion);
 
-	    // reset the default values
-	    *stats = *statsDefaults;
+            // reset the default values
+            *stats = *statsDefaults;
 
             // Use the selected background statistic for the first pass
-	    // If it fails, fall back on the "ROBUST_MEDIAN" version
-	    // If both fail, set the pixel to NAN and (below) interpolate
-	    // XXX psImageBackground will probably be renamed psImageStats
-	    // XXX don't bother trying if there are no valid pixels...
-
-	    if (psImageBackground(stats, subset, submask, 0xff, rng)) {
-		if (stats->options & PS_STAT_ROBUST_QUARTILE) {
-		    modelData[iy][ix] = stats->robustMedian;
-		} else {	
-		    modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation);
-		}
-		modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
-	    } else {
-		psStatsOptions currentOptions = stats->options;
-		stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
-		if (!psImageBackground(stats, subset, submask, 0xff, rng)) {
-		    psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
-			       "(%dx%d, (row0,col0) = (%d,%d)",
-			       subset->numRows, subset->numCols, subset->row0, subset->col0);
-		    modelData[iy][ix] = modelStdevData[iy][ix] = NAN;
-		} else {
-		    modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
-		    modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV);
-		}
-		// drop errors caused by psImageBackground failures
-		// XXX we probably should trap and exit on serious failures
-		psErrorClear(); 
-		stats->options = currentOptions;
-	    }
-	    modelData[iy][ix] += SKY_BIAS;
+            // If it fails, fall back on the "ROBUST_MEDIAN" version
+            // If both fail, set the pixel to NAN and (below) interpolate
+            // XXX psImageBackground will probably be renamed psImageStats
+            // XXX don't bother trying if there are no valid pixels...
+
+            if (psImageBackground(stats, subset, submask, 0xff, rng)) {
+                if (stats->options & PS_STAT_ROBUST_QUARTILE) {
+                    modelData[iy][ix] = stats->robustMedian;
+                } else {
+                    modelData[iy][ix] = psStatsGetValue(stats, statsOptionLocation);
+                }
+                modelStdevData[iy][ix] = psStatsGetValue(stats, statsOptionWidth);
+            } else {
+                psStatsOptions currentOptions = stats->options;
+                stats->options = PS_STAT_ROBUST_MEDIAN | PS_STAT_ROBUST_STDEV;
+                if (!psImageBackground(stats, subset, submask, 0xff, rng)) {
+                    psLogMsg ("psphot", PS_LOG_WARN, "Failed to estimate background using ROBUST_MEDIAN for "
+                               "(%dx%d, (row0,col0) = (%d,%d)",
+                               subset->numRows, subset->numCols, subset->row0, subset->col0);
+                    modelData[iy][ix] = modelStdevData[iy][ix] = NAN;
+                } else {
+                    modelData[iy][ix] = psStatsGetValue (stats, PS_STAT_ROBUST_MEDIAN);
+                    modelStdevData[iy][ix] = psStatsGetValue(stats, PS_STAT_ROBUST_STDEV);
+                }
+                // drop errors caused by psImageBackground failures
+                // XXX we probably should trap and exit on serious failures
+                psErrorClear();
+                stats->options = currentOptions;
+            }
+            modelData[iy][ix] += SKY_BIAS;
             psFree (subset);
             psFree (submask);
@@ -227,43 +227,43 @@
     // patch over bad regions (use average of 8 possible neighbor pixels)
     // XXX consider testing all pixels against the 8 neighbors and replacing outliers...
-    double Count = 0;			// number of good pixels
-    double Value = 0;			// sum of good pixel's value
-    double ValueStdev = 0;		// sum of good pixel's standard deviations
+    double Count = 0;                   // number of good pixels
+    double Value = 0;                   // sum of good pixel's value
+    double ValueStdev = 0;              // sum of good pixel's standard deviations
     for (int iy = 0; iy < model->image->numRows; iy++) {
         for (int ix = 0; ix < model->image->numCols; ix++) {
-	    if (!isnan(modelData[iy][ix])) {
-		Value += modelData[iy][ix];
-		ValueStdev += modelStdevData[iy][ix];
-		Count++;
-		continue;
-	    }
-
-	    double value = 0;
-	    double count = 0;
-	    for (int jy = iy - 1; jy <= iy + 1; jy++) {
-		if (jy <   0) continue;
-		if (jy >= model->image->numRows) continue;
-		for (int jx = ix - 1; jx <= ix + 1; jx++) {
-		    if (!jx && !jy) continue;
-		    if (jx   <   0) continue;
-		    if (jx   >= model->image->numCols) continue;
-		    value += modelData[jy][jx];
-		    count += 1.0;
-		}
-	    }
-	    if (count > 0) modelData[iy][ix] = value / count;
-	}
+            if (!isnan(modelData[iy][ix])) {
+                Value += modelData[iy][ix];
+                ValueStdev += modelStdevData[iy][ix];
+                Count++;
+                continue;
+            }
+
+            double value = 0;
+            double count = 0;
+            for (int jy = iy - 1; jy <= iy + 1; jy++) {
+                if (jy <   0) continue;
+                if (jy >= model->image->numRows) continue;
+                for (int jx = ix - 1; jx <= ix + 1; jx++) {
+                    if (!jx && !jy) continue;
+                    if (jx   <   0) continue;
+                    if (jx   >= model->image->numCols) continue;
+                    value += modelData[jy][jx];
+                    count += 1.0;
+                }
+            }
+            if (count > 0) modelData[iy][ix] = value / count;
+        }
     }
     assert (Count > 0);
     Value /= Count;
     ValueStdev /= Count;
-	    
+
     // patch over remaining bad regions (use global average)
     for (int iy = 0; iy < model->image->numRows; iy++) {
         for (int ix = 0; ix < model->image->numCols; ix++) {
-	    if (!isnan(modelData[iy][ix])) continue;
-	    modelData[iy][ix] = Value;
-	    modelStdevData[iy][ix] = ValueStdev;
-	}
+            if (!isnan(modelData[iy][ix])) continue;
+            modelData[iy][ix] = Value;
+            modelStdevData[iy][ix] = ValueStdev;
+        }
     }
 
@@ -273,25 +273,25 @@
     psMetadataAddF32(recipe, PS_LIST_TAIL, "SKY_STDEV", PS_META_REPLACE, "sky stdev", ValueStdev);
     psLogMsg ("psphot", PS_LOG_INFO, "image sky : mean %f stdev %f", Value, ValueStdev);
-    
+
     // measure image and background stats and save for later output
     psStats *statsBck = psStatsAlloc (PS_STAT_SAMPLE_MEAN |
-				      PS_STAT_SAMPLE_STDEV |
-				      PS_STAT_MIN |
-				      PS_STAT_MAX);
+                                      PS_STAT_SAMPLE_STDEV |
+                                      PS_STAT_MIN |
+                                      PS_STAT_MAX);
     psImageStats (statsBck, model->image, NULL, 0);
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MODEL_MEAN",
-		      PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
+                      PS_META_REPLACE, "sky model mean",          statsBck->sampleMean);
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MODEL_STDEV",
-		      PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
+                      PS_META_REPLACE, "sky model stdev",         statsBck->sampleStdev);
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MODEL_MAX",
-		      PS_META_REPLACE, "sky model maximum value", statsBck->max);
+                      PS_META_REPLACE, "sky model maximum value", statsBck->max);
     psMetadataAddF32 (recipe, PS_LIST_TAIL, "SKY_MODEL_MIN",
-		      PS_META_REPLACE, "sky model minimum value", statsBck->min);
+                      PS_META_REPLACE, "sky model minimum value", statsBck->min);
     psMetadataAddS32 (recipe, PS_LIST_TAIL, "SKY_MODEL_NX",
-		      PS_META_REPLACE, "sky model size (x)",      model->image->numCols);
+                      PS_META_REPLACE, "sky model size (x)",      model->image->numCols);
     psMetadataAddS32 (recipe, PS_LIST_TAIL, "SKY_MODEL_NY",
-		      PS_META_REPLACE, "sky model size (y)",      model->image->numRows);
-    psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f", 
-	      statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
+                      PS_META_REPLACE, "sky model size (y)",      model->image->numRows);
+    psLogMsg ("psphot", PS_LOG_INFO, "background sky : min %f mean %f max %f stdev %f",
+              statsBck->min, statsBck->sampleMean, statsBck->max, statsBck->sampleStdev);
     psFree (statsBck);
 
@@ -299,19 +299,19 @@
     file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKGND");
     if (file) {
-	// we are using PSPHOT.BACKGND as an I/O file: select readout or create
-	if (file->mode == PM_FPA_MODE_INTERNAL) {
-	    background = file->readout;
-	} else {
-	    background = pmFPAviewThisReadout (view, file->fpa);
-	}
-	if (background == NULL) {
-	    // readout does not yet exist: create from input
-	    pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view);
-	    background = pmFPAviewThisReadout (view, file->fpa);
-	    if ((image->numCols != background->image->numCols) || (image->numRows != background->image->numRows)) {
-		psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for background dimensions");
-		return false;
-	    }
-	}
+        // we are using PSPHOT.BACKGND as an I/O file: select readout or create
+        if (file->mode == PM_FPA_MODE_INTERNAL) {
+            background = file->readout;
+        } else {
+            background = pmFPAviewThisReadout (view, file->fpa);
+        }
+        if (background == NULL) {
+            // readout does not yet exist: create from input
+            pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view);
+            background = pmFPAviewThisReadout (view, file->fpa);
+            if ((image->numCols != background->image->numCols) || (image->numRows != background->image->numRows)) {
+                psError (PSPHOT_ERR_PROG, true, "inconsistent sizes for background dimensions");
+                return false;
+            }
+        }
     } else {
         background = pmFPAfileDefineInternal (config->files, "PSPHOT.BACKGND", image->numCols, image->numRows, PS_TYPE_F32);
@@ -321,8 +321,8 @@
     // linear interpolation to full-scale
     if (!psImageUnbin (background->image, model->image, binning)) {
-	psError (PSPHOT_ERR_PROG, true, "failed to build background image");
-	return false;
-    }
-	
+        psError (PSPHOT_ERR_PROG, true, "failed to build background image");
+        return false;
+    }
+
     psLogMsg ("psphot", PS_LOG_MINUTIA, "build resampled image: %f sec\n", psTimerMark ("psphot"));
 
@@ -330,23 +330,23 @@
     file = psMetadataLookupPtr (&status, config->files, "PSPHOT.BACKSUB");
     if (file) {
-	// we are using PSPHOT.BACKSUB as an I/O file: select readout or create
-	backSub = pmFPAviewThisReadout (view, file->fpa);
-	if (backSub == NULL) {
-	    // readout does not yet exist: create from input
-	    pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view);
-	    backSub = pmFPAviewThisReadout (view, file->fpa);
-	}
+        // we are using PSPHOT.BACKSUB as an I/O file: select readout or create
+        backSub = pmFPAviewThisReadout (view, file->fpa);
+        if (backSub == NULL) {
+            // readout does not yet exist: create from input
+            pmFPAfileCopyStructureView (file->fpa, inFPA, 1, 1, view);
+            backSub = pmFPAviewThisReadout (view, file->fpa);
+        }
     }
 
     if (psTraceGetLevel("psphot") > 5) {
-	char name[256];
-	sprintf (name, "image.%02d.fits", npass);
-	psphotSaveImage (NULL, image, name);
-	sprintf (name, "back.%02d.fits", npass);
-	psphotSaveImage (NULL, background->image, name);
-	sprintf (name, "mask.%02d.fits", npass);
-	psphotSaveImage (NULL, mask, name);
-	sprintf (name, "backmdl.%02d.fits", npass);
-	psphotSaveImage (NULL, model->image, name);
+        char name[256];
+        sprintf (name, "image.%02d.fits", npass);
+        psphotSaveImage (NULL, image, name);
+        sprintf (name, "back.%02d.fits", npass);
+        psphotSaveImage (NULL, background->image, name);
+        sprintf (name, "mask.%02d.fits", npass);
+        psphotSaveImage (NULL, mask, name);
+        sprintf (name, "backmdl.%02d.fits", npass);
+        psphotSaveImage (NULL, model->image, name);
     }
 
@@ -354,15 +354,15 @@
     for (int j = 0; j < image->numRows; j++) {
         for (int i = 0; i < image->numCols; i++) {
-	    image->data.F32[j][i] -= backData[j][i];
-	    if (backSub) {
-		backSub->image->data.F32[j][i] = image->data.F32[j][i];
-	    }
+            image->data.F32[j][i] -= backData[j][i];
+            if (backSub) {
+                backSub->image->data.F32[j][i] = image->data.F32[j][i];
+            }
         }
     }
 
     if (psTraceGetLevel("psphot") > 5) {
-	char name[256];
-	sprintf (name, "backsub.%02d.fits", npass);
-	psphotSaveImage (NULL, image, name);
+        char name[256];
+        sprintf (name, "backsub.%02d.fits", npass);
+        psphotSaveImage (NULL, image, name);
     }
     npass ++;