6 Force spectroscopy curves filtering of flat curves
7 Licensed under the GNU LGPL version 2
9 Other plugin dependencies:
10 procplots.py (plot processing plugin)
12 from libhooke import WX_GOOD
14 wxversion.select(WX_GOOD)
16 import xml.dom.minidom
21 from numpy import diff
25 import libpeakspot as lps
26 import libhookecurve as lhc
29 class flatfiltsCommands:
32 #configurate convfilt variables
33 convfilt_configurator=ConvfiltConfig()
35 #different OSes have different path conventions
36 if self.config['hookedir'][0]=='/':
37 slash='/' #a Unix or Unix-like system
39 slash='\\' #it's a drive letter, we assume it's Windows
41 self.convfilt_config=convfilt_configurator.load_config(self.config['hookedir']+slash+'convfilt.conf')
43 def do_flatfilt(self,args):
47 Filters out flat (featureless) curves of the current playlist,
48 creating a playlist containing only the curves with potential
52 flatfilt [min_npks min_deviation]
54 min_npks = minmum number of points over the deviation
57 min_deviation = minimum signal/noise ratio
60 If called without arguments, it uses default values, that
61 should work most of the times.
70 min_deviation=int(args[1])
74 print 'Processing playlist...'
78 for item in self.current_list:
82 notflat=self.has_features(item, median_filter, min_npks, min_deviation)
83 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': features are ',notflat
86 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.'
90 item.curve=None #empty the item object, to further avoid memory leak
91 notflat_list.append(item)
93 if len(notflat_list)==0:
94 print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent'
97 print 'Found ',len(notflat_list),' potentially interesting curves'
98 print 'Regenerating playlist...'
100 self.current_list=notflat_list
101 self.current=self.current_list[self.pointer]
104 def has_features(self,item,median_filter,min_npks,min_deviation):
106 decides if a curve is flat enough to be rejected from analysis: it sees if there
107 are at least min_npks points that are higher than min_deviation times the absolute value
110 Algorithm original idea by Francesco Musiani, with my tweaks and corrections.
114 item.identify(self.drivers)
115 #we assume the first is the plot with the force curve
116 #do the median to better resolve features from noise
117 flat_plot=self.plotmanip_median(item.curve.default_plots()[0], item, customvalue=median_filter)
118 flat_vects=flat_plot.vectors
119 item.curve.close_all()
120 #needed to avoid *big* memory leaks!
124 #absolute value of derivate
125 yretdiff=diff(flat_vects[1][1])
126 yretdiff=[abs(value) for value in yretdiff]
127 #average of derivate values
128 diffmean=numpy.mean(yretdiff)
132 for value in yretdiff:
133 if value/diffmean > min_deviation:
141 del flat_plot, flat_vects, yretdiff
145 ################################################################
146 #-----CONVFILT-------------------------------------------------
147 #-----Convolution-based peak recognition and filtering.
148 #Requires the libpeakspot.py library
150 def has_peaks(self, plot, abs_devs):
152 Finds peak position in a force curve.
153 FIXME: should be moved in libpeakspot.py
155 xret=plot.vectors[1][0]
156 yret=plot.vectors[1][1]
157 #Calculate convolution.
158 convoluted=lps.conv_dx(yret, self.convfilt_config['convolution'])
160 #surely cut everything before the contact point
161 cut_index=self.find_contact_point()
163 #cut even more, before the blind window
164 start_x=xret[cut_index]
166 for value in xret[cut_index:]:
167 if abs((value) - (start_x)) > self.convfilt_config['blindwindow']*(10**-9):
170 cut_index+=blind_index
172 #do the dirty convolution-peak finding stuff
173 noise_level=lps.noise_absdev(convoluted[cut_index:], self.convfilt_config['positive'], self.convfilt_config['maxcut'], self.convfilt_config['stable'])
174 above=lps.abovenoise(convoluted,noise_level,cut_index,abs_devs)
175 peak_location,peak_size=lps.find_peaks(above)
178 for i in range(len(peak_location)):
179 peak=peak_location[i]
180 maxpk=min(yret[peak-10:peak+10])
181 index_maxpk=yret[peak-10:peak+10].index(maxpk)+(peak-10)
182 peak_location[i]=index_maxpk
184 return peak_location,peak_size
187 def exec_has_peaks(self,item,abs_devs):
189 encapsulates has_peaks for the purpose of correctly treating the curve objects in the convfilt loop,
190 to avoid memory leaks
192 item.identify(self.drivers)
193 #we assume the first is the plot with the force curve
194 plot=item.curve.default_plots()[0]
196 if 'flatten' in self.config['plotmanips']:
197 #If flatten is present, use it for better recognition of peaks...
198 flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
199 plot=flatten(plot, item, customvalue=1)
201 peak_location,peak_size=self.has_peaks(plot,abs_devs)
202 #close all open files
203 item.curve.close_all()
204 #needed to avoid *big* memory leaks!
207 return peak_location, peak_size
209 #------------------------
210 #------commands----------
211 #------------------------
212 def do_peaks(self,args):
216 Test command for convolution filter / test.
218 Syntax: peaks [deviations]
219 absolute deviation = number of times the convolution signal is above the noise absolute deviation.
223 args=self.convfilt_config['mindeviation']
232 defplots=self.current.curve.default_plots()[0] #we need the raw, uncorrected plots
234 if 'flatten' in self.config['plotmanips']:
235 flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
236 defplots=flatten(defplots, self.current)
238 print 'You have the flatten plot manipulator not loaded. Enabling it could give you better results.'
240 peak_location,peak_size=self.has_peaks(defplots,abs_devs)
241 print 'Found '+str(len(peak_location))+' peaks.'
242 to_dump='peaks '+self.current.path+' '+str(len(peak_location))
243 self.outlet.push(to_dump)
246 #if no peaks, we have nothing to plot. exit.
247 if len(peak_location)==0:
250 #otherwise, we plot the peak locations.
251 xplotted_ret=self.plots[0].vectors[1][0]
252 yplotted_ret=self.plots[0].vectors[1][1]
253 xgood=[xplotted_ret[index] for index in peak_location]
254 ygood=[yplotted_ret[index] for index in peak_location]
256 recplot=self._get_displayed_plot()
257 recplot.vectors.append([xgood,ygood])
258 if recplot.styles==[]:
259 recplot.styles=[None,None,'scatter']
261 recplot.styles+=['scatter']
263 self._send_plot([recplot])
265 def do_convfilt(self,args):
269 Filters out flat (featureless) curves of the current playlist,
270 creating a playlist containing only the curves with potential
274 convfilt [min_npks min_deviation]
276 min_npks = minmum number of peaks
277 (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
279 min_deviation = minimum signal/noise ratio *in the convolution*
280 (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
282 If called without arguments, it uses default values.
285 min_npks=self.convfilt_config['minpeaks']
286 min_deviation=self.convfilt_config['mindeviation']
290 min_npks=int(args[0])
291 min_deviation=int(args[1])
295 print 'Processing playlist...'
296 print '(Please wait)'
301 for item in self.current_list:
305 peak_location,peak_size=self.exec_has_peaks(item,min_deviation)
306 if len(peak_location)>=min_npks:
310 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': found '+str(len(peak_location))+' peaks.'+isok
312 peak_location,peak_size=[],[]
313 print 'Curve',item.path, 'is',c,'of',len(self.current_list),': cannot be filtered. Probably unable to retrieve force data from corrupt file.'
315 if len(peak_location)>=min_npks:
316 item.peak_location=peak_location
317 item.peak_size=peak_size
318 item.curve=None #empty the item object, to further avoid memory leak
319 notflat_list.append(item)
321 #Warn that no flattening had been done.
322 if not ('flatten' in self.config['plotmanips']):
323 print 'Flatten manipulator was not found. Processing was done without flattening.'
324 print 'Try to enable it in your configuration file for better results.'
326 if len(notflat_list)==0:
327 print 'Found nothing interesting. Check your playlist, could be a bug or criteria could be too much stringent'
330 print 'Found ',len(notflat_list),' potentially interesting curves'
331 print 'Regenerating playlist...'
333 self.current_list=notflat_list
334 self.current=self.current_list[self.pointer]
337 def do_setconv(self,args):
341 Sets the convfilt configuration variables
343 Syntax: setconv variable value
346 #FIXME: a general "set dictionary" function has to be built
348 print self.convfilt_config
350 if not (args[0] in self.convfilt_config.keys()):
351 print 'This is not an internal convfilt variable!'
352 print 'Run "setconv" without arguments to see a list of defined variables.'
356 print self.convfilt_config[args[0]]
359 self.convfilt_config[args[0]]=eval(args[1])
360 except NameError: #we have a string argument
361 self.convfilt_config[args[0]]=args[1]
364 #########################
365 #HANDLING OF CONFIGURATION FILE
366 class ConvfiltConfig:
368 Handling of convfilt configuration file
370 Mostly based on the simple-yet-useful examples of the Python Library Reference
371 about xml.dom.minidom
373 FIXME: starting to look a mess, should require refactoring
380 def load_config(self, filename):
381 myconfig=file(filename)
382 #the following 3 lines are needed to strip newlines. otherwise, since newlines
383 #are XML elements too, the parser would read them (and re-save them, multiplying
385 #yes, I'm an XML n00b
386 the_file=myconfig.read()
387 the_file_lines=the_file.split('\n')
388 the_file=''.join(the_file_lines)
390 self.config_tree=xml.dom.minidom.parseString(the_file)
392 def getText(nodelist):
393 #take the text from a nodelist
394 #from Python Library Reference 13.7.2
396 for node in nodelist:
397 if node.nodeType == node.TEXT_NODE:
401 def handleConfig(config):
402 noiseabsdev_elements=config.getElementsByTagName("noise_absdev")
403 convfilt_elements=config.getElementsByTagName("convfilt")
404 handleAbsdev(noiseabsdev_elements)
405 handleConvfilt(convfilt_elements)
407 def handleAbsdev(noiseabsdev_elements):
408 for element in noiseabsdev_elements:
409 for attribute in element.attributes.keys():
410 self.config[attribute]=element.getAttribute(attribute)
412 def handleConvfilt(convfilt_elements):
413 for element in convfilt_elements:
414 for attribute in element.attributes.keys():
415 self.config[attribute]=element.getAttribute(attribute)
417 handleConfig(self.config_tree)
418 #making items in the dictionary machine-readable
419 for item in self.config.keys():
421 self.config[item]=eval(self.config[item])
422 except NameError: #if it's an unreadable string, keep it as a string