Changeset 13631


Ignore:
Timestamp:
05/10/12 11:33:36 (13 years ago)
Author:
neise
Message:
removed new class DRSSpikes_2D ... not needed
File:
1 edited

Legend:

Unmodified
Added
Removed
  • fact/tools/pyscripts/pyfact/drs_spikes.py

    r13458 r13631  
    8989        return data.reshape(self.row, self.col)
    9090
    91 class DRSSpikes_2D(object):
    92     """ remove spikes (single or double false readings) from DRS4 data
    93     Strategy:
    94     * filter the data, removing the signal, thus spike(s) are clearly visible
    95     * search single and double spikes
    96     * replace the spike by a value derived from the neighbors
    97    
    98     """
    99    
    100     def __init__(self, threshold=7.,
    101                  single_pattern=np.array( [-0.5, 1.0, -0.5]) ,
    102                  double_pattern=np.array([-1., 1., 1., -1.]),
    103                  user_action=lambda self, data: None,
    104                  debug = False):
    105         """ initialize spike filter
    106         template_single: template of a single slice spike
    107         template_double: template of a two slice spike
    108         """
    109 
    110         self.threshold = threshold
    111         self.single_pattern = single_pattern * threshold
    112         self.double_pattern = double_pattern * threshold
    113        
    114         self.remove_signal = fir.RemoveSignal()
    115 
    116         self.user_action = user_action
    117         self.debug = debug
    118 
    119     def __call__(self, data):
    120         # shortcuts
    121         row, col = data.shape
    122         thr = self.threshold
    123        
    124         # these are: lists, which will contain positiones of spikes
    125         # lets see if this is feasible
    126         self.singles = []
    127         self.doubles = []
    128         singles = self.singles
    129         doubles = self.doubles
    130 
    131         # indi means indicator, i.e. a filter output, which indicates, where spikes
    132         # are positioned in the unfiltered data
    133         # indi is delayed w.r.t. data by 1
    134         self.indi = self.remove_signal(data)
    135         indi = self.indi
    136 
    137         # cand (candidates), is a tuple of two equal length np.arrays
    138         # each pair of array elements can be understood as coordinates, pointing out
    139         # where the condition was fullfilled in indi.
    140         # e.g. in pixel = cand[0][0] around slice = cand[1][0]
    141         # there is probably a spike
    142         #
    143         # for double spikes, two neighboring slices fulfill the condition
    144         # which lead to something like:
    145         # cand =( array([ ... 3, 3, ... ]) , array([ ... 102, 103 ...]) )
    146         cand = np.where(indi[:,1:-2] > thr)
    147         self.cand = cand
    148         # in order to verify, that the candidate is really a single or double
    149         # spike, we compare the spike with a 3 or 4 slices pattern.
    150         # therefor we want to slice out 4 slices out of indi, where ever the
    151         # condition was fullfilled
    152         #
    153         # note: since indi was sliced in the np.where statement,
    154         # the resulting cand coordinates are reduced by 1 in the slice coordinate
    155         # this is actually what we want, since a spike has a distinctive low-high-low
    156         # pattern in the indicator. So we *want* the indicator slices to be shifted 1
    157         # to the left
    158         # and in addition, by pure chance, the coordinates in cand[1] point directly
    159         # to the real spike in data, since indi was delayed by one anyway.
    160         cand_slices = np.empty( (len(cand[0]), 4), dtype=np.int )
    161         self.cand_slices = cand_slices
    162         for i in range(4):
    163             cand_slices[i] = indi[ (cand[0], cand[1]+i )
    164            
    165         # search for single spikes
    166         sp = self.single_pattern * np.sign( self.single_pattern )
    167         for i, can in enumerate(cand_slices[:-1]):
    168             can *= np.sign(sp)
    169             if np.all( can > sp):
    170                 singles.append( (cand[0][i],can[1][i]) )
    171 
    172         # I guess in principle it is possible, that a candidate looks like a
    173         # single and like a double, ... nut with the current patters
    174         # but in case one changes the patterns ... then it might happen.
    175         # In addition the treatment of double spikes is maybed not smart:
    176         #   In case both parts of a double spike fulfill the 1st condition
    177         #       only the first candidate will fulfill the 2nd condition
    178         #   In case only the first part fulfilled the 1st conition, then
    179         #       we are fine
    180         #   In case only the second part triggered the first time, then
    181         #       we sliced out the wrong piece and it wouldn't fulfull the 2nd anyway.
    182         #
    183         #   This means, in case there are neighboring hits in cand,
    184         #   The 2nd neighbor will never fulfill the 2nd condition.
    185 
    186         # search for double spikes
    187         dp = self.double_pattern * np.sign( self.double_pattern )
    188         for i, can in enumerate( cand_slices ):
    189             can *= np.sign(dp)
    190             if np.all(can > dp):
    191                 doubles.append( (cand[0][i],can[1][i]) )
    192 
    193         self.user_action(self, data)
    194 
    195         data = self.remove_single_spikes(singles, data)
    196         data = self.remove_double_spikes(doubles, data)
    197         return data
    198 
    199     def remove_single_spikes(self, singles, data):
    200         for spike in singles:
    201             data[spike[0],spike[1]] = (data[spike[0],spike[1]-1] + data[spike[0],spike[1]+1]) / 2.
    202         return data
    203    
    204     def remove_double_spikes(self, doubles, data):
    205         for spike in doubles:
    206             data[spike[0],spike[1]:spike[1]+2] = (data[spike[0],spike[1]-1] + data[spike[0],spike[1]+2]) / 2.
    207         return data
    208 
    209 
    210 
    21191def _test():
    21292 
Note: See TracChangeset for help on using the changeset viewer.