Changeset 5743 for trunk/stac/src/stacScales.c

Timestamp:

Dec 7, 2005, 2:57:14 PM (21 years ago)

Author:

Paul Price

Message:

Importing from PAP cvs tree again. This will now be the working tree.

File:

• trunk/stac/src/stacScales.c (modified) (5 diffs)

trunk/stac/src/stacScales.c

@@ -14 +14 @@
     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 +27 @@
         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 +39 @@
 
     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 +169 @@
                         }
                     }
-                    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 +193 @@
         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);