Changeset 13458 for fact/tools/pyscripts/pyfact/drs_spikes.py

Timestamp:

04/26/12 11:43:46 (13 years ago)

Author:

neise

Message:

...testing...

File:

• fact/tools/pyscripts/pyfact/drs_spikes.py (modified) (1 diff)

fact/tools/pyscripts/pyfact/drs_spikes.py

@@ -89 +89 @@
         return data.reshape(self.row, self.col)
 
+class DRSSpikes_2D(object):
+    """ remove spikes (single or double false readings) from DRS4 data
+    Strategy:
+    * filter the data, removing the signal, thus spike(s) are clearly visible
+    * search single and double spikes
+    * replace the spike by a value derived from the neighbors
+    
+    """
+    
+    def __init__(self, threshold=7., 
+                 single_pattern=np.array( [-0.5, 1.0, -0.5]) ,
+                 double_pattern=np.array([-1., 1., 1., -1.]), 
+                 user_action=lambda self, data: None,
+                 debug = False):
+        """ initialize spike filter 
+        template_single: template of a single slice spike
+        template_double: template of a two slice spike
+        """
+
+        self.threshold = threshold
+        self.single_pattern = single_pattern * threshold
+        self.double_pattern = double_pattern * threshold
+        
+        self.remove_signal = fir.RemoveSignal()
+
+        self.user_action = user_action
+        self.debug = debug
+
+    def __call__(self, data):
+        # shortcuts
+        row, col = data.shape
+        thr = self.threshold
+        
+        # these are: lists, which will contain positiones of spikes
+        # lets see if this is feasible
+        self.singles = []
+        self.doubles = []
+        singles = self.singles
+        doubles = self.doubles
+
+        # indi means indicator, i.e. a filter output, which indicates, where spikes
+        # are positioned in the unfiltered data
+        # indi is delayed w.r.t. data by 1
+        self.indi = self.remove_signal(data)
+        indi = self.indi
+
+        # cand (candidates), is a tuple of two equal length np.arrays
+        # each pair of array elements can be understood as coordinates, pointing out
+        # where the condition was fullfilled in indi. 
+        # e.g. in pixel = cand[0][0] around slice = cand[1][0]
+        # there is probably a spike
+        #
+        # for double spikes, two neighboring slices fulfill the condition
+        # which lead to something like: 
+        # cand =( array([ ... 3, 3, ... ]) , array([ ... 102, 103 ...]) )
+        cand = np.where(indi[:,1:-2] > thr)
+        self.cand = cand
+        # in order to verify, that the candidate is really a single or double
+        # spike, we compare the spike with a 3 or 4 slices pattern.
+        # therefor we want to slice out 4 slices out of indi, where ever the 
+        # condition was fullfilled
+        #
+        # note: since indi was sliced in the np.where statement, 
+        # the resulting cand coordinates are reduced by 1 in the slice coordinate
+        # this is actually what we want, since a spike has a distinctive low-high-low 
+        # pattern in the indicator. So we *want* the indicator slices to be shifted 1
+        # to the left
+        # and in addition, by pure chance, the coordinates in cand[1] point directly 
+        # to the real spike in data, since indi was delayed by one anyway.
+        cand_slices = np.empty( (len(cand[0]), 4), dtype=np.int )
+        self.cand_slices = cand_slices
+        for i in range(4):
+            cand_slices[i] = indi[ (cand[0], cand[1]+i )
+            
+        # search for single spikes
+        sp = self.single_pattern * np.sign( self.single_pattern )
+        for i, can in enumerate(cand_slices[:-1]):
+            can *= np.sign(sp)
+            if np.all( can > sp):
+                singles.append( (cand[0][i],can[1][i]) )
+
+        # I guess in principle it is possible, that a candidate looks like a 
+        # single and like a double, ... nut with the current patters
+        # but in case one changes the patterns ... then it might happen.
+        # In addition the treatment of double spikes is maybed not smart:
+        #   In case both parts of a double spike fulfill the 1st condition
+        #       only the first candidate will fulfill the 2nd condition
+        #   In case only the first part fulfilled the 1st conition, then
+        #       we are fine
+        #   In case only the second part triggered the first time, then
+        #       we sliced out the wrong piece and it wouldn't fulfull the 2nd anyway.
+        # 
+        #   This means, in case there are neighboring hits in cand,
+        #   The 2nd neighbor will never fulfill the 2nd condition.
+
+        # search for double spikes
+        dp = self.double_pattern * np.sign( self.double_pattern )
+        for i, can in enumerate( cand_slices ):
+            can *= np.sign(dp)
+            if np.all(can > dp):
+                doubles.append( (cand[0][i],can[1][i]) )
+
+        self.user_action(self, data)
+
+        data = self.remove_single_spikes(singles, data)
+        data = self.remove_double_spikes(doubles, data)
+        return data
+
+    def remove_single_spikes(self, singles, data):
+        for spike in singles:
+            data[spike[0],spike[1]] = (data[spike[0],spike[1]-1] + data[spike[0],spike[1]+1]) / 2.
+        return data
+    
+    def remove_double_spikes(self, doubles, data):
+        for spike in doubles:
+            data[spike[0],spike[1]:spike[1]+2] = (data[spike[0],spike[1]-1] + data[spike[0],spike[1]+2]) / 2.
+        return data
+
+
 
 def _test():