6 Force spectroscopy files filtering of flat files.
9 procplots.py (plot processing plugin)
11 Copyright 2008 Massimo Sandal, Fabrizio Benedetti
12 with modifications by Dr. Rolf Schmidt (Concordia University, Canada)
14 This program is released under the GNU General Public License version 2.
17 import lib.libhooke as lh
19 wxversion.select(lh.WX_GOOD)
22 from numpy import diff, mean
24 import lib.peakspot as lps
27 class flatfiltsCommands:
29 def do_flatfilt(self):
33 Filters out flat (featureless) files of the current playlist,
34 creating a playlist containing only the files with potential
38 flatfilt [min_npks min_deviation]
40 min_npks = minmum number of points over the deviation
43 min_deviation = minimum signal/noise ratio
46 If called without arguments, it uses default values, that
47 should work most of the times.
50 self.AppendToOutput('Processing playlist...')
51 self.AppendToOutput('(Please wait)')
53 playlist = self.GetActivePlaylist()
54 files = playlist.files
56 for current_file in files:
59 current_file.identify(self.drivers)
60 notflat = self.has_features(copy.deepcopy(current_file))
64 feature_string = str(notflat) + ' features'
66 feature_string = 'no features'
68 feature_string = '1 feature'
69 output_string = ''.join(['Curve ', current_file.name, '(', str(file_index), '/', str(len(files)), '): ', feature_string])
72 output_string = ''.join(['Curve ', current_file.name, '(', str(file_index), '/', str(len(files)), '): cannot be filtered. Probably unable to retrieve force data from corrupt file.'])
73 self.AppendToOutput(output_string)
75 current_file.features = notflat
76 features.append(file_index - 1)
78 self.AppendToOutput('Found nothing interesting. Check the playlist, could be a bug or criteria could be too stringent.')
80 if len(features) < playlist.count:
81 self.AppendToOutput(''.join(['Found ', str(len(features)), ' potentially interesting files.']))
82 self.AppendToOutput('Regenerating playlist...')
83 playlist_filtered = playlist.filter_curves(features)
84 self.AddPlaylist(playlist_filtered, name='flatfilt')
86 self.AppendToOutput('No files filtered. Try different filtering criteria.')
88 def has_features(self, current_file):
90 decides if a curve is flat enough to be rejected from analysis: it sees if there
91 are at least min_npks points that are higher than min_deviation times the absolute value
94 Algorithm original idea by Francesco Musiani, with my tweaks and corrections.
96 #TODO: shoudl medianfilter be variable?
98 #medianfilter = self.GetIntFromConfig('flatfilts', 'flatfilt', 'median_filter')
99 mindeviation = self.GetIntFromConfig('flatfilts', 'flatfilt', 'min_deviation')
100 minpeaks = self.GetIntFromConfig('flatfilts', 'flatfilt', 'min_npks')
104 #we assume the first is the plot with the force curve
105 #do the median to better resolve features from noise
106 flat_curve = self.plotmanip_median(current_file.plot, current_file, customvalue=medianfilter)
108 #absolute value of derivate
109 yretdiff = diff(flat_curve.curves[lh.RETRACTION].y)
110 yretdiff = [abs(value) for value in yretdiff]
111 #average of derivate values
112 diffmean = mean(yretdiff)
116 for value in yretdiff:
117 if value / diffmean > mindeviation:
122 if c_pks >= minpeaks:
127 ################################################################
128 #-----CONVFILT-------------------------------------------------
129 #-----Convolution-based peak recognition and filtering.
130 #Requires the peakspot.py library
132 def has_peaks(self, plot=None, plugin=None):
134 Finds peak position in a force curve.
138 blindwindow = self.GetFloatFromConfig('flatfilts', 'convfilt', 'blindwindow')
139 #need to convert the string that contains the list into a list
140 convolution = eval(self.GetStringFromConfig('flatfilts', 'convfilt', 'convolution'))
141 maxcut = self.GetFloatFromConfig('flatfilts', 'convfilt', 'maxcut')
142 mindeviation = self.GetFloatFromConfig('flatfilts', 'convfilt', 'mindeviation')
143 positive = self.GetBoolFromConfig('flatfilts', 'convfilt', 'positive')
144 seedouble = self.GetIntFromConfig('flatfilts', 'convfilt', 'seedouble')
145 stable = self.GetFloatFromConfig('flatfilts', 'convfilt', 'stable')
147 blindwindow = self.GetFloatFromConfig(plugin.name, plugin.section, plugin.prefix + 'blindwindow')
148 #need to convert the string that contains the list into a list
149 convolution = eval(self.GetStringFromConfig(plugin.name, plugin.section, plugin.prefix + 'convolution'))
150 maxcut = self.GetFloatFromConfig(plugin.name, plugin.section, plugin.prefix + 'maxcut')
151 mindeviation = self.GetFloatFromConfig(plugin.name, plugin.section, plugin.prefix + 'mindeviation')
152 positive = self.GetBoolFromConfig(plugin.name, plugin.section, plugin.prefix + 'positive')
153 seedouble = self.GetIntFromConfig(plugin.name, plugin.section, plugin.prefix + 'seedouble')
154 stable = self.GetFloatFromConfig(plugin.name, plugin.section, plugin.prefix + 'stable')
157 plot = self.GetDisplayedPlotCorrected()
159 retraction = plot.curves[lh.RETRACTION]
160 #Calculate convolution
161 convoluted = lps.conv_dx(retraction.y, convolution)
163 #surely cut everything before the contact point
164 cut_index = self.find_contact_point(plot)
165 #cut even more, before the blind window
166 start_x = retraction.x[cut_index]
168 for value in retraction.x[cut_index:]:
169 if abs((value) - (start_x)) > blindwindow * (10 ** -9):
172 cut_index += blind_index
173 #do the dirty convolution-peak finding stuff
174 noise_level = lps.noise_absdev(convoluted[cut_index:], positive, maxcut, stable)
175 above = lps.abovenoise(convoluted, noise_level, cut_index, mindeviation)
176 peak_location, peak_size = lps.find_peaks(above, seedouble=seedouble)
178 for i in range(len(peak_location)):
179 peak = peak_location[i]
180 maxpk = min(retraction.y[peak - 10:peak + 10])
181 index_maxpk = retraction.y[peak - 10:peak + 10].index(maxpk) + (peak - 10)
182 peak_location[i] = index_maxpk
184 return peak_location, peak_size
186 def do_peaks(self, plugin=None, peak_location=None, peak_size=None):
188 Test command for convolution filter.
190 Syntax: peaks [deviations]
191 absolute deviation = number of times the convolution signal is above the noise absolute deviation.
196 color = self.GetColorFromConfig('flatfilts', 'peaks', 'color')
197 size = self.GetIntFromConfig('flatfilts', 'peaks', 'size')
199 color = self.GetColorFromConfig(plugin.name, plugin.section, plugin.prefix + 'color')
200 size = self.GetIntFromConfig(plugin.name, plugin.section, plugin.prefix + 'size')
202 plot = self.GetDisplayedPlotCorrected()
204 if peak_location is None and peak_size is None:
205 if not self.AppliesPlotmanipulator('flatten'):
206 self.AppendToOutput('The flatten plot manipulator is not loaded. Enabling it could give better results.')
208 peak_location, peak_size = self.has_peaks(plot)
209 if len(peak_location) != 1:
213 self.AppendToOutput('Found ' + str(len(peak_location)) + peak_str)
216 retraction = plot.curves[lh.RETRACTION]
218 peaks = lib.curve.Curve()
221 peaks.style = 'scatter'
222 peaks.title = 'Peaks'
223 peaks.x = [retraction.x[index] for index in peak_location]
224 peaks.y = [retraction.y[index] for index in peak_location]
226 plot.curves.append(peaks)
228 self.UpdatePlot(plot)
230 def do_convfilt(self):
232 Filters out flat (featureless) files of the current playlist,
233 creating a playlist containing only the files with potential
236 min_npks: minmum number of peaks
237 min_deviation: minimum signal/noise ratio *in the convolution*
240 self.AppendToOutput('Processing playlist...')
241 self.AppendToOutput('(Please wait)')
242 apply_plotmanipulators = self.GetStringFromConfig('flatfilts', 'convfilt', 'apply_plotmanipulators')
243 minpeaks = self.GetIntFromConfig('flatfilts', 'convfilt', 'minpeaks')
245 playlist = self.GetActivePlaylist()
247 files = self.GetActivePlaylist().files
249 for current_file in files:
253 current_file.identify(self.drivers)
254 if apply_plotmanipulators == 'all':
255 plot = self.ApplyPlotmanipulators(current_file.plot, current_file)
256 if apply_plotmanipulators == 'flatten':
257 plotmanipulator = self.GetPlotmanipulator('flatten')
258 plot = plotmanipulator.method(current_file.plot, current_file)
259 if apply_plotmanipulators == 'none':
260 plot = copy.deepcopy(current_file.plot)
262 peak_location, peak_size = self.has_peaks(plot)
263 number_of_peaks = len(peak_location)
264 if number_of_peaks != 1:
265 if number_of_peaks > 0:
266 feature_string = str(number_of_peaks) + ' features'
268 feature_string = 'no features'
270 feature_string = '1 feature'
271 if number_of_peaks >= minpeaks:
272 feature_string += '+'
273 output_string = ''.join(['Curve ', current_file.name, '(', str(file_index), '/', str(len(files)), '): ', feature_string])
277 output_string = ''.join(['Curve ', current_file.name, '(', str(file_index), '/', str(len(files)), '): cannot be filtered. Probably unable to retrieve force data from corrupt file.'])
278 self.AppendToOutput(output_string)
279 if number_of_peaks >= minpeaks:
280 current_file.peak_location = peak_location
281 current_file.peak_size = peak_size
282 features.append(file_index - 1)
284 #Warn that no flattening had been done.
285 if not self.HasPlotmanipulator('plotmanip_flatten'):
286 self.AppendToOutput('Flatten manipulator was not found. Processing was done without flattening.')
288 if not self.AppliesPlotmanipulator('flatten'):
289 self.AppendToOutput('Flatten manipulator was not applied.')
290 self.AppendToOutput('Try to enable the flatten plotmanipulator for better results.')
293 self.AppendToOutput('Found nothing interesting. Check the playlist, could be a bug or criteria could be too stringent.')
295 if len(features) < playlist.count:
296 self.AppendToOutput(''.join(['Found ', str(len(features)), ' potentially interesting files.']))
297 self.AppendToOutput('Regenerating playlist...')
298 playlist_filtered = playlist.filter_curves(features)
299 self.AddPlaylist(playlist_filtered, name='convfilt')
301 self.AppendToOutput('No files filtered. Try different filtering criteria.')