From: W. Trevor King Date: Wed, 2 Jun 2010 19:26:20 +0000 (-0400) Subject: Translate find_contact_point_fmms to new architecture. X-Git-Url: http://git.tremily.us/?p=hooke.git;a=commitdiff_plain;h=ca3b8810184e9ec49637d7e58bb700ee369e56b2 Translate find_contact_point_fmms to new architecture. Also * Add some docstrings, comments, and line wrappings. * Rename _find_contact_point* -> find_contact_point* so users will look at the docstrings. * Added curve to find_contact_point*() for driver-specific hacks and workarounds (specifically, the PicoForce trigger bug). --- diff --git a/hooke/plugin/vclamp.py b/hooke/plugin/vclamp.py index b83d5de..c57c11a 100644 --- a/hooke/plugin/vclamp.py +++ b/hooke/plugin/vclamp.py @@ -28,6 +28,7 @@ common velocity clamp analysis tasks. import copy import numpy +import scipy from ..command import Command, Argument, Failure, NullQueue from ..config import Setting @@ -100,7 +101,8 @@ class SurfacePositionModel (ModelFitter): if r2 >= 1: self._model_data[:r2] = p[0] + p[1] * numpy.arange(r2) if r2 < len(self._data)-1: - self._model_data[r2:] = p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2) + self._model_data[r2:] = \ + p[0] + p[1]*p[2] + p[3] * numpy.arange(len(self._data)-r2) return self._model_data def set_data(self, data, info=None): @@ -168,7 +170,7 @@ class SurfacePositionModel (ModelFitter): params.append(0.) # restore the non-contact slope parameter # check that the fit is reasonable, see the :meth:`model` docstring - # for parameter descriptions. + # for parameter descriptions. HACK: hardcoded cutoffs. if abs(params[3]*10) > abs(params[1]) : raise PoorFit('Slope in non-contact region, or no slope in contact') if params[2] < self.info['min position']+0.02*self.info['position range']: @@ -214,6 +216,14 @@ class SurfaceContactCommand (Command): The adjusted data columns `surface distance (m)` and `surface adjusted deflection (m)` are also added to the block. + + You can select the contact point algorithm with the creatively + named `surface contact point algorithm` configuration setting. + Currently available options are: + + * fmms (:meth:`find_contact_point_fmms`) + * ms (:meth:`find_contact_point_ms`) + * wtk (:meth:`find_contact_point_wtk`) """ def __init__(self, plugin): super(SurfaceContactCommand, self).__init__( @@ -243,8 +253,8 @@ selects the retracting curve. ['surface distance (m)', 'surface adjusted deflection (m)']) z_data = data[:,data.info['columns'].index('z piezo (m)')] d_data = data[:,data.info['columns'].index('deflection (m)')] - i,deflection_offset,ps = self._find_contact_point( - z_data, d_data, outqueue) + i,deflection_offset,ps = self.find_contact_point( + params['curve'], z_data, d_data, outqueue) surface_offset = z_data[i] new.info['surface distance offset (m)'] = surface_offset new.info['surface deflection offset (m)'] = deflection_offset @@ -254,88 +264,78 @@ selects the retracting curve. data = params['curve'].data[int(params['block'])] # HACK, int() should be handled by ui params['curve'].data[int(params['block'])] = new # HACK, int() should be handled by ui - def _find_contact_point(self, z_data, d_data, outqueue=None): - fn = getattr(self, '_find_contact_point_%s' + def find_contact_point(self, curve, z_data, d_data, outqueue=None): + """Railyard for the `find_contact_point_*` family. + + Uses the `surface contact point algorithm` configuration + setting to call the appropriate backend algorithm. + """ + fn = getattr(self, 'find_contact_point_%s' % self.plugin.config['surface contact point algorithm']) - return fn(z_data, d_data, outqueue) + return fn(curve, z_data, d_data, outqueue) - def _find_contact_point_fmms(self, z_data, d_data, outqueue=None): + def find_contact_point_fmms(self, curve, z_data, d_data, outqueue=None): """Algorithm by Francesco Musiani and Massimo Sandal. Notes ----- Algorithm: - * take care of the PicoForce trigger bug - exclude retraction portions with too high standard deviation - * fit the second half of the retraction curve to a line - * if the fit is not almost horizontal, take a smaller chunk and repeat - * otherwise, we have something horizontal - * so take the average of horizontal points and use it as a baseline - - Then, start from the rise of the retraction curve and look at - the first point below the baseline. + 0) Driver-specific workarounds, e.g. deal with the PicoForce + trigger bug by excluding retraction portions with excessive + deviation. + 1) Select the second half (non-contact side) of the retraction + curve. + 2) Fit the selection to a line. + 3) If the fit is not almost horizontal, halve the selection + and retrun to (2). + 4) Average the selection and use it as a baseline. + 5) Slide in from the start (contact side) of the retraction + curve, until you find a point with greater than baseline + deflection. That point is the contact point. """ - if not plot: - plot=self.plots[0] - - outplot=self.subtract_curves(1) - xret=outplot.vectors[1][0] - ydiff=outplot.vectors[1][1] - - xext=plot.vectors[0][0] - yext=plot.vectors[0][1] - xret2=plot.vectors[1][0] - yret=plot.vectors[1][1] - - #taking care of the picoforce trigger bug: we exclude portions of the curve that have too much - #standard deviation. yes, a lot of magic is here. - monster=True - monlength=len(xret)-int(len(xret)/20) - finalength=len(xret) - while monster: - monchunk=scipy.array(ydiff[monlength:finalength]) - if abs(np.std(monchunk)) < 2e-10: - monster=False - else: #move away from the monster - monlength-=int(len(xret)/50) - finalength-=int(len(xret)/50) - - - #take half of the thing - endlength=int(len(xret)/2) - - ok=False - - while not ok: - xchunk=yext[endlength:monlength] - ychunk=yext[endlength:monlength] - regr=scipy.stats.linregress(xchunk,ychunk)[0:2] - #we stop if we found an almost-horizontal fit or if we're going too short... - #FIXME: 0.1 and 6 here are "magic numbers" (although reasonable) - if (abs(regr[1]) > 0.1) and ( endlength < len(xret)-int(len(xret)/6) ) : - endlength+=10 - else: - ok=True - - - ymean=np.mean(ychunk) #baseline - - index=0 - point = ymean+1 - - #find the first point below the calculated baseline - while point > ymean: - try: - point=yret[index] - index+=1 - except IndexError: - #The algorithm didn't find anything below the baseline! It should NEVER happen - index=0 - return index - - return index - - def _find_contact_point_ms(self, z_data, d_data, outqueue=None): + if curve.info['filetype'] == 'picoforce': + # Take care of the picoforce trigger bug (TODO: example + # data file demonstrating the bug). We exclude portions + # of the curve that have too much standard deviation. + # Yes, a lot of magic is here. + check_start = len(d_data)-len(d_data)/20 + monster_start = len(d_data) + while True: + # look at the non-contact tail + non_monster = d_data[check_start:monster_start] + if non_monster.std() < 2e-10: # HACK: hardcoded cutoff + break + else: # move further away from the monster + check_start -= len(d_data)/50 + monster_start -= len(d_data)/50 + z_data = z_data[:monster_start] + d_data = d_data[:monster_start] + + # take half of the thing to start + selection_start = len(d_data)/2 + while True: + z_chunk = z_data[selection_start:] + d_chunk = d_data[selection_start:] + slope,intercept,r,two_tailed_prob,stderr_of_the_estimate = \ + scipy.stats.linregress(z_chunk, d_chunk) + # We stop if we found an almost-horizontal fit or if we're + # getting to small a selection. FIXME: 0.1 and 5./6 here + # are "magic numbers" (although reasonable) + if (abs(slope) < 0.1 # deflection (m) / surface (m) + or selection_start > 5./6*len(d_data)): + break + selection_start += 10 + + d_baseline = d_chunk.mean() + + # find the first point above the calculated baseline + i = 0 + while i < len(d_data) and d_data[i] < ymean: + i += 1 + return (i, d_baseline, {}) + + def find_contact_point_ms(self, curve, z_data, d_data, outqueue=None): """Algorithm by Massimo Sandal. Notes @@ -347,13 +347,6 @@ selects the retracting curve. * ? """ - #raw_plot=self.current.curve.default_plots()[0] - raw_plot=self.plots[0] - '''xext=self.plots[0].vectors[0][0] - yext=self.plots[0].vectors[0][1] - xret2=self.plots[0].vectors[1][0] - yret=self.plots[0].vectors[1][1] - ''' xext=raw_plot.vectors[0][0] yext=raw_plot.vectors[0][1] xret2=raw_plot.vectors[1][0] @@ -387,7 +380,7 @@ selects the retracting curve. else: return dummy.index - def _find_contact_point_wtk(self, z_data, d_data, outqueue=None): + def find_contact_point_wtk(self, curve, z_data, d_data, outqueue=None): """Algorithm by W. Trevor King. Notes