Index: trunk/stac/src/stacScales.c
===================================================================
--- trunk/stac/src/stacScales.c	(revision 3673)
+++ trunk/stac/src/stacScales.c	(revision 5743)
@@ -14,9 +14,9 @@
     assert(image->type.type == PS_TYPE_F32);
 
-// Will use robust median instead of sampling --- it's supposed to be fast.
 #if 1
     int size = image->numCols * image->numRows;	// Number of pixels in image
     int numSamples = size / sample;	// Number of samples in image
     psVector *values = psVectorAlloc(numSamples + 1, PS_TYPE_F32); // Vector containing sub-sample
+    psVector *mask = psVectorAlloc(numSamples + 1, PS_TYPE_U8);	// Mask for sample
 
     int offset = 0;			// Offset from start of the row
@@ -27,4 +27,9 @@
 	while (col < image->numCols) {
 	    values->data.F32[index] = image->data.F32[row][col];
+	    if (isnan(values->data.F32[index])) {
+		mask->data.U8[index] = 1;
+	    } else {
+		mask->data.U8[index] = 0;
+	    }
 	    col += sample;
 	    index++;
@@ -34,9 +39,11 @@
 
     psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
-    stats = psVectorStats(stats, values, NULL, NULL, 0);
+    stats = psVectorStats(stats, values, NULL, mask, 1);
     float median = stats->sampleMedian;
     psFree(stats);
     psFree(values);
+    psFree(mask);
 #else
+    // Will use robust median instead of sampling --- it's supposed to be fast.
     psStats *stats = psStatsAlloc(PS_STAT_ROBUST_MEDIAN); // Using a clipped mean because median is SLOW
     stats->clipSigma = 3.0;
@@ -162,12 +169,16 @@
 			}
 		    }
-		    psTrace("stac.scales", 9, "Star at %f,%f --> %f\n", coords->x, coords->y, sum);
 		    // Subtract background, renormalise to account for circular aperture
 		    if (numPix > 0 && sum > 0) {
 			sum -= offsets->data.F32[i] * (float)numPix;
 			photometry->data.F32[j] = sum * M_PI * aper2 / (float)numPix;
+			if (photometry->data.F32[j] > 0 && finite(photometry->data.F32[j])) {
+			    mask->data.U8[j] = 1;
+			    psTrace("stac.scales", 8, "Star at %f,%f --> %f\n", coords->x, coords->y, sum);
+			} else {
+			    mask->data.U8[j] = 0;
+			}
+		    } else {
 			mask->data.U8[j] = 0;
-		    } else {
-			mask->data.U8[j] = 1;
 		    }
 		}
@@ -182,18 +193,42 @@
 	psVector *ref = stars->data[0];	// The reference photometry
 	psVector *refMask = masks->data[0]; // The reference mask
+#if 0
 	psStats *stats = psStatsAlloc(PS_STAT_CLIPPED_MEAN); // Statistics
 	stats->clipSigma = 2.5;
 	stats->clipIter = 3;
+#else
+	psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN); // Statistics
+#endif
 	for (int i = 1; i < scales->n; i++) {
-	    psVector *compare = stars->data[i];	// The comparison photometry
-	    psVector *compareMask = masks->data[i]; // The comparison mask
-
-	    compare = (psVector*)psBinaryOp(compare, compare, "/", ref);
-	    compareMask = (psVector*)psBinaryOp(compareMask, compareMask, "+", refMask);
-
-	    stats = psVectorStats(stats, compare, NULL, compareMask, 3); // Use maskVal of 3 to catch 1 and 2
-
+	    psVector *input = stars->data[i];	// The comparison photometry
+	    psVector *inputMask = masks->data[i]; // The comparison mask
+
+	    psVector *compare = (psVector*)psBinaryOp(NULL, input, "/", ref);
+	    psVector *compareMask = (psVector*)psBinaryOp(NULL, inputMask, "*", refMask);
+	    (void)psBinaryOp(compareMask, psScalarAlloc(1, PS_TYPE_U8), "-", compareMask);
+
+#if 0
+	    psTrace("stac.scales", 9, "Getting scale for image %d...\n", i);
+	    for (int j = 0; j < compare->n; j++) {
+		if (compareMask->data.U8[j] & 1) {
+		    psTrace("stac.scales", 9, "Bad star %d: %f %f --> %f %d\n", j, input->data.F32[j],
+			    ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
+		} else {
+		    psTrace("stac.scales", 9, "Good star %d: %f %f --> %f %d\n", j, input->data.F32[j],
+			    ref->data.F32[j], compare->data.F32[j], compareMask->data.U8[j]);
+		}
+	    }
+#endif
+
+	    psVectorStats(stats, compare, NULL, compareMask, 1);
+
+#if 0
 	    scales->data.F32[i] = stats->clippedMean;
+#else
+	    scales->data.F32[i] = stats->sampleMedian;
+#endif
 	    psTrace("stac.scales", 5, "Scale for image %d is %f\n", i, scales->data.F32[i]);
+	    psFree(compare);
+	    psFree(compareMask);
 	}
 	psFree(stats);