Merge pull request #42 from stuwilkins/add_std_filter

ENH: Added filter to stddev

csxtools/fastccd/phocount.py

@@ -1,7 +1,7 @@
 from ..ext import phocount as ph
 
 
-def photon_count(data, thresh, mean_filter, nsum=3, nan=False):
+def photon_count(data, thresh, mean_filter, std_filter, nsum=3, nan=False):
     """Do single photon counting on CCD image
 
     This routine does single photon counting by cluster analysis. The image
@@ -19,6 +19,9 @@ def photon_count(data, thresh, mean_filter, nsum=3, nan=False):
     mean_filter : tuple
         Filter only the values of the mean which are within the limits of
         the tuple of the form (min, max)
+    std_filter : tuple
+        Filter only the values of the standard deviation  which are within
+        the limits of the tuple of the form (min, max)
     nsum : int
         The number of pixels to use to calculate the energy deposited by the
         photon. This should be 0 < nsum <= 9.
@@ -33,4 +36,4 @@ def photon_count(data, thresh, mean_filter, nsum=3, nan=False):
         photon hit. The second array is the standard deviation for the
         integrated intensity on each photon hit.
     """
-    return ph.count(data, thresh, mean_filter, nsum, nan)
+    return ph.count(data, thresh, mean_filter, std_filter, nsum, nan)

src/phocount.c

@@ -44,7 +44,8 @@
 
 int count(data_t *in, data_t *out, data_t *stddev, 
           int ndims, index_t *dims, 
-          data_t *thresh, data_t *mean_filter, int sum_max, int nan){
+          data_t *thresh, data_t *sum_filter, data_t *std_filter,
+          int sum_max, int nan){
   index_t nimages = dims[0];
   index_t M = dims[ndims-1];
   index_t N = dims[ndims-2];
@@ -64,103 +65,112 @@ int count(data_t *in, data_t *out, data_t *stddev,
   }   
 
   index_t i;
-#pragma omp parallel for shared(in, out, stddev) private(i)
-  for(i=0;i<nimages;i++){
-    // Find the start pointers of the image
-    data_t *inp = in + (i*imsize) - 1;
-    data_t *outp = out + (i*imsize) - 1;
-    data_t *stddevp = stddev + (i*imsize) - 1;
-
-    data_t pixel[9];
-
-    // Clear out the parts of the output array we don't use
-
-    index_t j,k;
-    for(j=0;j<(M+1);j++){
-      inp++;
-      outp++;
-      stddevp++;
-      *outp = nodata;
-      *stddevp = nodata;
-    }
-
-    // Now start the search
-    for(j=1;j<(N-1);j++){
-      for(k=1;k<(M-1);k++){
+#pragma omp parallel shared(in, out, stddev) 
+  {
+#pragma omp for
+    for(i=0;i<nimages;i++){
+      // Find the start pointers of the image
+      data_t *inp = in + (i*imsize) - 1;
+      data_t *outp = out + (i*imsize) - 1;
+      data_t *stddevp = stddev + (i*imsize) - 1;
+      data_t pixel[9];
+
+      // Clear out the parts of the output array we don't use
+
+      index_t j, k;
+      for(j=0;j<(M+1);j++){
         inp++;
         outp++;
         stddevp++;
-
         *outp = nodata;
         *stddevp = nodata;
+      }
 
-        if((*inp < thresh[0]) || (*inp >= thresh[1])){
-          continue;
-        }
+      // Now start the search
+      for(j=1;j<(N-1);j++){
+        for(k=1;k<(M-1);k++){
+          inp++;
+          outp++;
+          stddevp++;
+
+          *outp = nodata;
+          *stddevp = nodata;
 
-        // The pixel is above thresh
-        // Now get the surrounding 9 pixels. 
-        pixel[0] = *inp;
-        pixel[1] = *(inp - M - 1);
-        pixel[2] = *(inp - M);
-        pixel[3] = *(inp - M + 1);
-        pixel[4] = *(inp - 1);
-        pixel[5] = *(inp + 1);
-        pixel[6] = *(inp + M - 1);
-        pixel[7] = *(inp + M);
-        pixel[8] = *(inp + M + 1);
-
-        // Is this the brightest pixel?
+          if((*inp < thresh[0]) || (*inp >= thresh[1])){
+            continue;
+          }
+
+          // The pixel is above thresh
+          // Now get the surrounding 9 pixels. 
+
+          pixel[0] = *inp;
+          pixel[1] = *(inp - M - 1);
+          pixel[2] = *(inp - M);
+          pixel[3] = *(inp - M + 1);
+          pixel[4] = *(inp - 1);
+          pixel[5] = *(inp + 1);
+          pixel[6] = *(inp + M - 1);
+          pixel[7] = *(inp + M);
+          pixel[8] = *(inp + M + 1);
+
+          // Is this the brightest pixel?
+
+          int n;
+          int flag = 0;
+          for(n=1;n<9;n++){
+            if(pixel[n] > pixel[0]){
+              flag = 1;
+              break;
+            }
+          }
+
+          if(flag){
+            continue;
+          }
 
-        int n;
-        int flag = 0;
-        for(n=1;n<9;n++){
-          if(pixel[n] > pixel[0]){
-            flag = 1;
-            break;
+          // Sort the array
+          sort(pixel, 9);
+
+          data_t sum = 0;
+          data_t scnd_moment = 0;
+
+          for(n=0;n<sum_max;n++){
+            sum += pixel[n];
+            scnd_moment += pixel[n] * pixel[n];
           }
-        }
 
-        if(flag){
-          continue;
-        }
-
-        // Sort the array
-        sort(pixel, 9);
-
-        data_t sum = 0;
-        data_t scnd_moment = 0;
-        for(n=0;n<sum_max;n++){
-          sum += pixel[n];
-          scnd_moment += pixel[n] * pixel[n];
-        }
+          if((sum < sum_filter[0]) || (sum >= sum_filter[1])){
+            continue;
+          }
 
-        if((sum < mean_filter[0]) || (sum >= mean_filter[1])){
-          continue;
-        }
+          data_t std = pow((scnd_moment - (sum*sum) / sum_max) / sum_max, 0.5);
+
+          if((std < std_filter[0]) || (std >= std_filter[1])){
+            continue;
+          }
+
+          *stddevp = std;
+          *outp = sum;
 
-        data_t var = (scnd_moment - (sum*sum) / sum_max) / sum_max;
-        *stddevp = pow(var, 0.5);
-        *outp = sum;
+        } // for(k)
 
-      } // for(k)
+        for(k=0;k<2;k++){
+          outp++;
+          stddevp++;
+          inp++;
+          *stddevp = nodata;
+          *outp = nodata;
+        }
+      } // for(j)
 
-      for(k=0;k<2;k++){
+      for(j=0;j<(M+1);j++){
         outp++;
         stddevp++;
-        inp++;
-        *stddevp = nodata;
         *outp = nodata;
+        *stddevp = nodata;
       }
-    } // for(j)
-
-    for(j=0;j<(M+1);j++){
-      outp++;
-      stddevp++;
-      *outp = nodata;
-      *stddevp = nodata;
-    }
-  } // for(nimages)
+    } // for(nimages)
+  } // pragma omp 
 
   return 0;
 }

src/phocount.h

@@ -44,7 +44,8 @@ typedef float data_t;
 
 int count(data_t *in, data_t *out, data_t *stddev, 
           int ndims, index_t *dims, 
-          data_t *thresh, data_t *mean_filter, int sum_max, int nan);
+          data_t *thresh, data_t *sum_filter, data_t *std_filter,
+          int sum_max, int nan);
 void sort(data_t *array, int n);
 
 #endif

src/phocountmodule.c

@@ -52,12 +52,13 @@ static PyObject* phocount_count(PyObject *self, PyObject *args){
   PyArrayObject *out = NULL;
   npy_intp *dims;
   int ndims;
-  float thresh[2], mean_filter[2];
+  float thresh[2], sum_filter[2], std_filter[2];
   int sum_max;
   int nan = 0;
 
-  if(!PyArg_ParseTuple(args, "O(ff)(ff)i|p", &_input, &thresh[0], &thresh[1],
-                                             &mean_filter[0], &mean_filter[1], 
+  if(!PyArg_ParseTuple(args, "O(ff)(ff)(ff)i|p", &_input, &thresh[0], &thresh[1],
+                                             &sum_filter[0], &sum_filter[1], 
+                                             &std_filter[0], &std_filter[1], 
                                              &sum_max, &nan)){
     return NULL;
   }
@@ -92,7 +93,8 @@ static PyObject* phocount_count(PyObject *self, PyObject *args){
   // Ok now we don't touch Python Object ... Release the GIL
   Py_BEGIN_ALLOW_THREADS
 
-  count(input_p, out_p, stddev_p, ndims, dims, thresh, mean_filter, sum_max, nan);
+  count(input_p, out_p, stddev_p, ndims, dims, thresh, 
+        sum_filter, std_filter, sum_max, nan);
 
   Py_END_ALLOW_THREADS
 

tests/test_fastccd.py

@@ -35,7 +35,8 @@ def test_photon_count():
                               dtype=np.float32)[:nsum])
 
     op = photon_count(np.array([x, x, x], dtype=np.float32),
-                      thresh=(5, 13), mean_filter=(10, 30), nsum=nsum)
+                      thresh=(5, 13), mean_filter=(10, 30),
+                      std_filter=(0, 100), nsum=nsum)
 
     assert_array_equal(op[0], np.array([y, y, y]))
     assert_array_almost_equal(op[1], np.array([z, z, z]), decimal=6)