Index: /tags/ipp-20101206/psModules/src/concepts/pmConceptsAverage.c
===================================================================
--- /tags/ipp-20101206/psModules/src/concepts/pmConceptsAverage.c	(revision 30036)
+++ /tags/ipp-20101206/psModules/src/concepts/pmConceptsAverage.c	(revision 30037)
@@ -136,4 +136,51 @@
     return average;
 }
+float medianWithDropouts (psList *sources, char *name) {
+
+    bool status;
+
+    psListIterator *sourcesIter = psListIteratorAlloc(sources, PS_LIST_HEAD, false); // Iterator for sources
+    pmCell *cell = NULL;                // Source cell from iteration
+
+    psStats *stats = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
+    psVector *values = psVectorAlloc(sources->n, PS_TYPE_F32);
+    int nvalues = 0;
+    while ((cell = psListGetAndIncrement(sourcesIter))) {
+        if (!cell) {
+            continue;
+        }
+
+        float value = psMetadataLookupF32(&status, cell->concepts, name);
+        if (!status) continue;
+        if (!isfinite(value)) continue;
+
+        values->data.F32[nvalues++] = value;
+    }
+    psFree (sourcesIter);
+    if (!nvalues) {
+        psWarning("no valid values found for %s\n", name);
+        psFree(values);
+        psFree(stats);
+        return INFINITY;
+    }
+    if (!(values = psVectorRealloc(values, nvalues))) {
+        psWarning("failed to reallocate values vector for %s\n", name);
+        psFree(stats);
+        return INFINITY;
+    }
+    if (!psVectorStats(stats, values, NULL, NULL, 0)) {
+        psWarning("psVectorStats failed for %s\n", name);
+        psFree(values);
+        psFree(stats);
+        return INFINITY;
+    }
+
+    psF32 median = psStatsGetValue(stats, PS_STAT_SAMPLE_MEDIAN);
+
+    psFree(values);
+    psFree(stats);
+
+    return median;
+}
 
 // Set a variety of concepts in a cell by averaging over several
@@ -145,5 +192,4 @@
     PS_ASSERT_INT_POSITIVE(sources->n, false);
 
-    float saturation = INFINITY;        // Saturation level
     float bad        = -INFINITY;       // Bad level
     double time      = 0.0;             // Time of observation
@@ -158,4 +204,5 @@
     float exposure  = averageWithDropouts (sources, "CELL.EXPOSURE");
     float darktime  = averageWithDropouts (sources, "CELL.DARKTIME");
+    float saturation = medianWithDropouts(sources, "CELL.SATURATION");
 
     // other concepts are a bit more "special"
@@ -221,13 +268,4 @@
         }
 
-        float cellSaturation = psMetadataLookupF32(NULL, cell->concepts, "CELL.SATURATION");
-	if (cellSaturation > 10000) {
-	    // do not allow invalid values to polute this calculation
-	    // XXX really need to do this on the basis of the fraction of the cell that contributes..
-	    // if a cell is completely masked, it should not be included.
-	    if (cellSaturation < saturation) {
-		saturation = cellSaturation;
-	    }
-	}
         float cellBad = psMetadataLookupF32(NULL, cell->concepts, "CELL.BAD");
         if (cellBad > bad) {