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

Timestamp:

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

Author:

Paul Price

Message:

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

File:

• trunk/stac/src/stacRejection.c (modified) (6 diffs)

trunk/stac/src/stacRejection.c

@@ -74 +74 @@
     psTrace("stac.rejection", 1, "Mapping rejection masks back to source....\n");
 
+    // Vectors for calculating the mean gradient
+    psVector *grads = psVectorAlloc(nImages, PS_TYPE_F32); // Gradient for each image
+    psVector *gradsMask = psVectorAlloc(nImages, PS_TYPE_U8); // Mask for gradient vector
+
     // Transform rejection masks back to source
     psArray *inputRej = psArrayAlloc(nImages);
@@ -107 +111 @@
         // calculate derivatives of the map, and use that as a buffer around the transformed position
         // in the input image.
+        psStats *median = psStatsAlloc(PS_STAT_SAMPLE_MEDIAN);
         for (int y = 0; y < nyInput; y++) {
             for (int x = 0; x < nxInput; x++) {
@@ -135 +140 @@
                                     (yPix >= 0) && (yPix <= ((psImage*)(inputs->data[j]))->numRows - 1)) {
                                     // Calculate the gradient
-                                    meanGrads += stacGradient(inputs->data[j], xPix, yPix);
+                                    grads->data.F32[j] = stacGradient(inputs->data[j], xPix, yPix);
+                                    gradsMask->data.U8[j] = 0;
                                     numGrads++;
+                                } else {
+                                    gradsMask->data.U8[j] = 1; // Mask this one
                                 }
+                            } else {
+                                gradsMask->data.U8[j] = 1; // Mask this one
                             }
                         }
                         if (numGrads > 0) {
-                            meanGrads /= (float)numGrads;
+                            (void)psVectorStats(median, grads, NULL, gradsMask, (psU8)1);
+                            meanGrads = median->sampleMedian;
                         } else {
-                            meanGrads = 0;
+                            meanGrads = 0.0;
                         }
 
 #ifdef TESTING
-                        gradient->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
-#endif
-
-                        if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
+                        //gradient->data.F32[y][x] = stacGradient(inputs->data[i], x, y) / meanGrads;
+                        gradient->data.F32[y][x] = meanGrads;
+#endif
+
+                        //if (stacGradient(inputs->data[i], x, y) < grad * meanGrads) {
+                        if (meanGrads > grad) {
                             mask->data.U8[y][x] = 1;
                             nBad++;
@@ -173 +186 @@
             }
         } // Iterating over pixels
+        psFree(median);
 
 #ifdef CRFLUX
@@ -194 +208 @@
         psFits *rejmapFile = psFitsAlloc(rejmapName);
         psFits *gradFile = psFitsAlloc(gradName);
-        if (!psFitsWriteImage(maskFile, NULL, mask, 0, NULL)) {
+        if (!psFitsWriteImage(maskFile, NULL, mask, 0)) {
             psErrorStackPrint(stderr, "Unable to write image: %s\n", maskName);
         }
         psTrace("stac", 1, "Mask image written to %s\n", maskName);
-        if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0, NULL)) {
+        if (!psFitsWriteImage(rejmapFile, NULL, rejmap, 0)) {
             psErrorStackPrint(stderr, "Unable to write image: %s\n", rejmapName);
         }
         psTrace("stac", 1, "Rejection map written to %s\n", rejmapName);
-        if (!psFitsWriteImage(gradFile, NULL, gradient, 0, NULL)) {
+        if (!psFitsWriteImage(gradFile, NULL, gradient, 0)) {
             psErrorStackPrint(stderr, "Unable to write image: %s\n", gradName);
         }
@@ -218 +232 @@
     }
 
+    psFree(grads);
+    psFree(gradsMask);
 
     psFree(inCoords);