Changeset 13631 for fact/tools/pyscripts/pyfact/drs_spikes.py
- Timestamp:
- 05/10/12 11:33:36 (13 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
fact/tools/pyscripts/pyfact/drs_spikes.py
r13458 r13631 89 89 return data.reshape(self.row, self.col) 90 90 91 class DRSSpikes_2D(object):92 """ remove spikes (single or double false readings) from DRS4 data93 Strategy:94 * filter the data, removing the signal, thus spike(s) are clearly visible95 * search single and double spikes96 * replace the spike by a value derived from the neighbors97 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 filter106 template_single: template of a single slice spike107 template_double: template of a two slice spike108 """109 110 self.threshold = threshold111 self.single_pattern = single_pattern * threshold112 self.double_pattern = double_pattern * threshold113 114 self.remove_signal = fir.RemoveSignal()115 116 self.user_action = user_action117 self.debug = debug118 119 def __call__(self, data):120 # shortcuts121 row, col = data.shape122 thr = self.threshold123 124 # these are: lists, which will contain positiones of spikes125 # lets see if this is feasible126 self.singles = []127 self.doubles = []128 singles = self.singles129 doubles = self.doubles130 131 # indi means indicator, i.e. a filter output, which indicates, where spikes132 # are positioned in the unfiltered data133 # indi is delayed w.r.t. data by 1134 self.indi = self.remove_signal(data)135 indi = self.indi136 137 # cand (candidates), is a tuple of two equal length np.arrays138 # each pair of array elements can be understood as coordinates, pointing out139 # 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 spike142 #143 # for double spikes, two neighboring slices fulfill the condition144 # which lead to something like:145 # cand =( array([ ... 3, 3, ... ]) , array([ ... 102, 103 ...]) )146 cand = np.where(indi[:,1:-2] > thr)147 self.cand = cand148 # in order to verify, that the candidate is really a single or double149 # 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 the151 # condition was fullfilled152 #153 # note: since indi was sliced in the np.where statement,154 # the resulting cand coordinates are reduced by 1 in the slice coordinate155 # this is actually what we want, since a spike has a distinctive low-high-low156 # pattern in the indicator. So we *want* the indicator slices to be shifted 1157 # to the left158 # and in addition, by pure chance, the coordinates in cand[1] point directly159 # 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_slices162 for i in range(4):163 cand_slices[i] = indi[ (cand[0], cand[1]+i )164 165 # search for single spikes166 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 a173 # single and like a double, ... nut with the current patters174 # 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 condition177 # only the first candidate will fulfill the 2nd condition178 # In case only the first part fulfilled the 1st conition, then179 # we are fine180 # In case only the second part triggered the first time, then181 # 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 spikes187 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 data198 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 data203 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 data208 209 210 211 91 def _test(): 212 92
Note:
See TracChangeset
for help on using the changeset viewer.