hooke/plugin/flatfilts-rolf.py

   1 # Copyright (C) 2010 W. Trevor King <wking@drexel.edu>
   2 #
   3 # This file is part of Hooke.
   4 #
   5 # Hooke is free software: you can redistribute it and/or
   6 # modify it under the terms of the GNU Lesser General Public
   7 # License as published by the Free Software Foundation, either
   8 # version 3 of the License, or (at your option) any later version.
   9 #
  10 # Hooke is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU Lesser General Public License for more details.
  14 #
  15 # You should have received a copy of the GNU Lesser General Public
  16 # License along with Hooke.  If not, see
  17 # <http://www.gnu.org/licenses/>.
  18
  19 """Force spectroscopy curves filtering of flat curves
  20
  21 Other plugin dependencies:
  22 procplots.py (plot processing plugin)
  23 """
  24
  25 import xml.dom.minidom
  26
  27 import wx
  28 #import scipy
  29 import numpy
  30 from numpy import diff
  31 import os.path
  32
  33 #import pickle
  34
  35 import libpeakspot as lps
  36 #import curve as lhc
  37 import hookecurve as lhc
  38 import libhooke as lh
  39 import wxversion
  40 wxversion.select(lh.WX_GOOD)
  41
  42
  43 class flatfiltsCommands:
  44
  45     def do_flatfilt(self):
  46         '''
  47         FLATFILT
  48         (flatfilts.py)
  49         Filters out flat (featureless) curves of the current playlist,
  50         creating a playlist containing only the curves with potential
  51         features.
  52         ------------
  53         Syntax:
  54         flatfilt [min_npks min_deviation]
  55
  56         min_npks = minmum number of points over the deviation
  57         (default=4)
  58
  59         min_deviation = minimum signal/noise ratio
  60         (default=9)
  61
  62         If called without arguments, it uses default values, that
  63         should work most of the times.
  64         '''
  65         #TODO: should this be optional?
  66         medianfilter = 7
  67
  68         self.AppendToOutput('Processing playlist...')
  69         self.AppendToOutput('(Please wait)')
  70         features = []
  71         playlist = self.GetActivePlaylist()
  72         curves = playlist.curves
  73         curve_index = 0
  74         for curve in curves:
  75             curve_index += 1
  76             try:
  77                 notflat = self.has_features(curve)
  78                 feature_string = ''
  79                 if notflat != 1:
  80                     if notflat > 0:
  81                         feature_string = str(notflat) + ' features'
  82                     else:
  83                         feature_string = 'no features'
  84                 else:
  85                     feature_string = '1 feature'
  86                 output_string = ''.join(['Curve ', curve.name, '(', str(curve_index), '/', str(len(curves)), '): ', feature_string])
  87             except:
  88                 notflat = False
  89                 output_string = ''.join(['Curve ', curve.name, '(', str(curve_index), '/', str(len(curves)), '): cannot be filtered. Probably unable to retrieve force data from corrupt file.'])
  90             self.AppendToOutput(output_string)
  91             if notflat:
  92                 curve.features = notflat
  93                 features.append(curve_index - 1)
  94         if not features:
  95             self.AppendToOutput('Found nothing interesting. Check the playlist, could be a bug or criteria could be too stringent.')
  96         else:
  97             if len(features) < playlist.count:
  98                 self.AppendToOutput(''.join(['Found ', str(len(features)), ' potentially interesting curves.']))
  99                 self.AppendToOutput('Regenerating playlist...')
 100                 playlist_filtered = playlist.filter_curves(features)
 101                 self.AddPlaylist(playlist_filtered, name='flatfilt')
 102             else:
 103                 self.AppendToOutput('No curves filtered. Try different filtering criteria.')
 104
 105     def has_features(self, curve):
 106         '''
 107         decides if a curve is flat enough to be rejected from analysis: it sees if there
 108         are at least min_npks points that are higher than min_deviation times the absolute value
 109         of noise.
 110
 111         Algorithm original idea by Francesco Musiani, with my tweaks and corrections.
 112         '''
 113         medianfilter = 7
 114         mindeviation = self.GetIntFromConfig('flatfilts', 'flatfilt', 'min_deviation')
 115         minpeaks = self.GetIntFromConfig('flatfilts', 'flatfilt', 'min_npks')
 116         #medianfilter = self.GetIntFromConfig('flatfilt', 'median_filter')
 117         #mindeviation = self.GetIntFromConfig('convfilt', 'mindeviation')
 118         #minpeaks = self.GetIntFromConfig('convfilt', 'minpeaks')
 119
 120         retvalue = 0
 121
 122         #item.identify(self.drivers)
 123         #we assume the first is the plot with the force curve
 124         #do the median to better resolve features from noise
 125         flat_plot = self.plotmanip_median(curve.driver.default_plots()[0], curve, customvalue=medianfilter)
 126         flat_vects = flat_plot.vectors
 127         curve.driver.close_all()
 128         #needed to avoid *big* memory leaks!
 129         #del item.driver
 130         #del item
 131
 132         #absolute value of derivate
 133         yretdiff=diff(flat_vects[1][1])
 134         yretdiff=[abs(value) for value in yretdiff]
 135         #average of derivate values
 136         diffmean=numpy.mean(yretdiff)
 137         yretdiff.sort()
 138         yretdiff.reverse()
 139         c_pks=0
 140         for value in yretdiff:
 141             if value/diffmean > mindeviation:
 142                 c_pks += 1
 143             else:
 144                 break
 145
 146         if c_pks >= minpeaks:
 147             retvalue = c_pks
 148
 149         del flat_plot, flat_vects, yretdiff
 150
 151         return retvalue
 152
 153     ################################################################
 154     #-----CONVFILT-------------------------------------------------
 155     #-----Convolution-based peak recognition and filtering.
 156     #Requires the libpeakspot.py library
 157
 158     def has_peaks(self, plot, curve=None):
 159         '''
 160         Finds peak position in a force curve.
 161         FIXME: should be moved to libpeakspot.py
 162         '''
 163
 164         blindwindow = self.GetFloatFromConfig('flatfilts', 'convfilt', 'blindwindow')
 165         #need to convert the string that contains the list into a list
 166         convolution = eval(self.GetStringFromConfig('flatfilts', 'convfilt', 'convolution'))
 167         maxcut = self.GetFloatFromConfig('flatfilts', 'convfilt', 'maxcut')
 168         mindeviation = self.GetFloatFromConfig('flatfilts', 'convfilt', 'mindeviation')
 169         positive = self.GetBoolFromConfig('flatfilts', 'convfilt', 'positive')
 170         seedouble = self.GetIntFromConfig('flatfilts', 'convfilt', 'seedouble')
 171         stable = self.GetFloatFromConfig('flatfilts', 'convfilt', 'stable')
 172
 173         xret = plot.vectors[1][0]
 174         yret = plot.vectors[1][1]
 175         #Calculate convolution
 176         convoluted = lps.conv_dx(yret, convolution)
 177
 178         #surely cut everything before the contact point
 179         cut_index = self.find_contact_point(plot, curve)
 180         #cut even more, before the blind window
 181         start_x = xret[cut_index]
 182         blind_index = 0
 183         for value in xret[cut_index:]:
 184             if abs((value) - (start_x)) > blindwindow * (10 ** -9):
 185                 break
 186             blind_index += 1
 187         cut_index += blind_index
 188         #do the dirty convolution-peak finding stuff
 189         noise_level = lps.noise_absdev(convoluted[cut_index:], positive, maxcut, stable)
 190         above = lps.abovenoise(convoluted, noise_level, cut_index, mindeviation)
 191         peak_location, peak_size = lps.find_peaks(above, seedouble=seedouble)
 192         #take the maximum
 193         for i in range(len(peak_location)):
 194             peak = peak_location[i]
 195             maxpk = min(yret[peak - 10:peak + 10])
 196             index_maxpk = yret[peak - 10:peak + 10].index(maxpk) + (peak - 10)
 197             peak_location[i] = index_maxpk
 198
 199         return peak_location, peak_size
 200
 201     def exec_has_peaks(self, curve):
 202         '''
 203         encapsulates has_peaks for the purpose of correctly treating the curve objects in the convfilt loop,
 204         to avoid memory leaks
 205         '''
 206         #item.identify(self.drivers)
 207         #we assume the first is the plot with the force curve
 208         plot = curve.driver.default_plots()[0]
 209
 210         if self.HasPlotmanipulator('plotmanip_flatten'):
 211             #If flatten is present, use it for better recognition of peaks...
 212             #flatten = self._find_plotmanip('flatten') #extract flatten plot manipulator
 213             #plot = flatten(plot, item, customvalue=1)
 214             plot = self.plotmanip_flatten(plot, curve, customvalue=1)
 215
 216         peak_location, peak_size = self.has_peaks(plot, curve)
 217         #close all open files
 218         curve.driver.close_all()
 219         #needed to avoid *big* memory leaks!
 220         #del item.driver
 221         #del item
 222         return peak_location, peak_size
 223
 224     #------------------------
 225     #------commands----------
 226     #------------------------
 227     def do_peaks(self):
 228         '''
 229         PEAKS
 230         (flatfilts.py)
 231         Test command for convolution filter / test.
 232         ----
 233         Syntax: peaks [deviations]
 234         absolute deviation = number of times the convolution signal is above the noise absolute deviation.
 235         Default is 5.
 236         '''
 237
 238         #TODO: check if the following line gives us what we need
 239         curve = self.GetActiveCurve()
 240         defplots = curve.driver.default_plots()[0] #we need the raw, uncorrected plots
 241
 242         #if 'flatten' in self.config['plotmanips']:
 243         if self.HasPlotmanipulator('plotmanip_flatten'):
 244             #flatten=self._find_plotmanip('flatten') #extract flatten plot manipulator
 245             #defplots=flatten(defplots, self.current)
 246             defplots = self.plotmanip_flatten(defplots, curve, customvalue=0)
 247         else:
 248             self.AppendToOutput('The flatten plot manipulator is not loaded. Enabling it could give better results.')
 249
 250         peak_location, peak_size = self.has_peaks(defplots, curve)
 251         self.AppendToOutput('Found ' + str(len(peak_location)) + ' peaks.')
 252         self.AppendToOutput('peaks ' + curve.filename + ' ' + str(len(peak_location)))
 253         #to_dump = 'peaks ' + current_curve.filename + ' ' + str(len(peak_location))
 254         #self.outlet.push(to_dump)
 255
 256         #if no peaks, we have nothing to plot. exit.
 257         if peak_location:
 258             #otherwise, we plot the peak locations.
 259             xplotted_ret = curve.plots[0].vectors[1][0]
 260             yplotted_ret = curve.plots[0].vectors[1][1]
 261             xgood = [xplotted_ret[index] for index in peak_location]
 262             ygood = [yplotted_ret[index] for index in peak_location]
 263
 264             curve.plots[0].add_set(xgood, ygood)
 265             curve.plots[0].styles.append('scatter')
 266             curve.plots[0].colors.append('indigo')
 267
 268             #recplot = self._get_displayed_plot()
 269             #recplot.vectors.append([xgood,ygood])
 270             #if recplot.styles == []:
 271                 #recplot.styles = [None, None, 'scatter']
 272                 #recplot.colors = [None, None, None]
 273             #else:
 274                 #recplot.styles += ['scatter']
 275                 #recplot.colors += [None]
 276
 277             #self._send_plot([recplot])
 278             self.UpdatePlot()
 279
 280     def do_convfilt(self):
 281         '''
 282         CONVFILT
 283         (flatfilts.py)
 284         Filters out flat (featureless) curves of the current playlist,
 285         creating a playlist containing only the curves with potential
 286         features.
 287         ------------
 288         Syntax:
 289         convfilt [min_npks min_deviation]
 290
 291         min_npks = minmum number of peaks
 292         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
 293
 294         min_deviation = minimum signal/noise ratio *in the convolution*
 295         (to set the default, see convfilt.conf file; CONVCONF and SETCONF commands)
 296
 297         If called without arguments, it uses default values.
 298         '''
 299
 300         self.AppendToOutput('Processing playlist...')
 301         self.AppendToOutput('(Please wait)')
 302         minpeaks = self.GetIntFromConfig('flatfilts', 'convfilt', 'minpeaks')
 303         features = []
 304         playlist = self.GetActivePlaylist()
 305
 306         curves = self.GetActivePlaylist().curves
 307         curve_index = 0
 308         for curve in curves:
 309             curve_index += 1
 310             try:
 311                 peak_location, peak_size = self.exec_has_peaks(curve)
 312                 number_of_peaks = len(peak_location)
 313                 if number_of_peaks != 1:
 314                     if number_of_peaks > 0:
 315                         feature_string = str(number_of_peaks) + ' features'
 316                     else:
 317                         feature_string = 'no features'
 318                 else:
 319                     feature_string = '1 feature'
 320                 if number_of_peaks >= minpeaks:
 321                     feature_string += '+'
 322                 output_string = ''.join(['Curve ', curve.name, '(', str(curve_index), '/', str(len(curves)), '): ', feature_string])
 323             except:
 324                 peak_location = []
 325                 peak_size = []
 326                 output_string = ''.join(['Curve ', curve.name, '(', str(curve_index), '/', str(len(curves)), '): cannot be filtered. Probably unable to retrieve force data from corrupt file.'])
 327             self.AppendToOutput(output_string)
 328             if number_of_peaks >= minpeaks:
 329                 curve.peak_location = peak_location
 330                 curve.peak_size = peak_size
 331                 features.append(curve_index - 1)
 332
 333         #TODO: do we need this? Flattening might not be necessary/desired
 334         #Warn that no flattening had been done.
 335         if not self.HasPlotmanipulator('plotmanip_flatten'):
 336             self.AppendToOutput('Flatten manipulator was not found. Processing was done without flattening.')
 337             self.AppendToOutput('Try to enable it in the configuration file for better results.')
 338         if not features:
 339             self.AppendToOutput('Found nothing interesting. Check the playlist, could be a bug or criteria could be too stringent.')
 340         else:
 341             if len(features) < playlist.count:
 342                 self.AppendToOutput(''.join(['Found ', str(len(features)), ' potentially interesting curves.']))
 343                 self.AppendToOutput('Regenerating playlist...')
 344                 playlist_filtered = playlist.filter_curves(features)
 345                 self.AddPlaylist(playlist_filtered, name='convfilt')
 346             else:
 347                 self.AppendToOutput('No curves filtered. Try different filtering criteria.')
 348