+ self._commands = [
+ SurfaceContactCommand(self), ForceCommand(self),
+ ]
+
+ def default_settings(self):
+ return [
+ Setting(section=self.setting_section, help=self.__doc__),
+ Setting(section=self.setting_section,
+ option='surface contact point algorithm',
+ value='wtk',
+ help='Select the surface contact point algorithm. See the documentation for descriptions of available algorithms.')
+ ]
+
+
+class SurfaceContactCommand (Command):
+ """Automatically determine a block's surface contact point.
+
+ Uses the block's `z piezo (m)` and `deflection (m)` arrays.
+ Stores the contact parameters in `block.info`'s `surface distance
+ offset (m)` and `surface deflection offset (m)`. Model-specific
+ fitting information is stored in `surface detection parameters`.
+
+ The adjusted data columns `surface distance (m)` and `surface
+ adjusted deflection (m)` are also added to the block.
+ """
+ def __init__(self, plugin):
+ super(SurfaceContactCommand, self).__init__(
+ name='zero block surface contact point',
+ arguments=[
+ CurveArgument,
+ Argument(name='block', aliases=['set'], type='int', default=0,
+ help="""
+Data block for which the force should be calculated. For an
+approach/retract force curve, `0` selects the approaching curve and `1`
+selects the retracting curve.
+""".strip()),
+ ],
+ help=self.__doc__, plugin=plugin)
+
+ def _run(self, hooke, inqueue, outqueue, params):
+ data = params['curve'].data[int(params['block'])] # HACK, int() should be handled by ui
+ # HACK? rely on params['curve'] being bound to the local hooke
+ # playlist (i.e. not a copy, as you would get by passing a
+ # curve through the queue). Ugh. Stupid queues. As an
+ # alternative, we could pass lookup information through the
+ # queue...
+ new = Data((data.shape[0], data.shape[1]+2), dtype=data.dtype)
+ new.info = copy.deepcopy(data.info)
+ new[:,:-2] = data
+ new.info['columns'].extend(
+ ['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)
+ surface_offset = z_data[i]
+ new.info['surface distance offset (m)'] = surface_offset
+ new.info['surface deflection offset (m)'] = deflection_offset
+ new.info['surface detection parameters'] = ps
+ new[:,-2] = z_data - surface_offset
+ new[:,-1] = d_data - deflection_offset
+ 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'
+ % self.plugin.config['surface contact point algorithm'])
+ return fn(z_data, d_data, outqueue)
+
+ def _find_contact_point_fmms(self, 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.
+ """
+ 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):
+ """Algorithm by Massimo Sandal.
+
+ Notes
+ -----
+ WTK: At least the commits are by Massimo, and I see no notes
+ attributing the algorithm to anyone else.
+
+ Algorithm:
+
+ * ?
+ """
+ #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]
+ yret=raw_plot.vectors[1][1]
+
+ first_point=[xext[0], yext[0]]
+ last_point=[xext[-1], yext[-1]]
+
+ #regr=scipy.polyfit(first_point, last_point,1)[0:2]
+ diffx=abs(first_point[0]-last_point[0])
+ diffy=abs(first_point[1]-last_point[1])
+
+ #using polyfit results in numerical errors. good old algebra.
+ a=diffy/diffx
+ b=first_point[1]-(a*first_point[0])
+ baseline=scipy.polyval((a,b), xext)
+
+ ysub=[item-basitem for item,basitem in zip(yext,baseline)]
+
+ contact=ysub.index(min(ysub))
+
+ return xext,ysub,contact
+
+ #now, exploit a ClickedPoint instance to calculate index...
+ dummy=ClickedPoint()
+ dummy.absolute_coords=(x_intercept,y_intercept)
+ dummy.find_graph_coords(xret2,yret)
+
+ if debug:
+ return dummy.index, regr, regr_contact
+ else:
+ return dummy.index
+
+ def _find_contact_point_wtk(self, z_data, d_data, outqueue=None):
+ """Algorithm by W. Trevor King.
+
+ Notes
+ -----
+ Uses :func:`analyze_surf_pos_data_wtk` internally.
+ """
+ reverse = z_data[0] > z_data[-1]
+ if reverse == True: # approaching, contact region on the right
+ d_data = d_data[::-1]
+ s = SurfacePositionModel(d_data)
+ s.info['force zero non-contact slope'] = True
+ offset,contact_slope,surface_index,non_contact_slope = s.fit(
+ outqueue=outqueue)
+ info = {
+ 'offset': offset,
+ 'contact slope': contact_slope,
+ 'surface index': surface_index,
+ 'non-contact slope': non_contact_slope,
+ 'reversed': reverse,
+ }
+ deflection_offset = offset + contact_slope*surface_index,
+ if reverse == True:
+ surface_index = len(d_data)-1-surface_index
+ return (numpy.round(surface_index), deflection_offset, info)
+
+class ForceCommand (Command):
+ """Calculate a block's `deflection (N)` array.
+
+ Uses the block's `deflection (m)` array and `spring constant
+ (N/m)`.
+ """
+ def __init__(self, plugin):
+ super(ForceCommand, self).__init__(
+ name='add block force array',
+ arguments=[
+ CurveArgument,
+ Argument(name='block', aliases=['set'], type='int', default=0,
+ help="""
+Data block for which the force should be calculated. For an
+approach/retract force curve, `0` selects the approaching curve and `1`
+selects the retracting curve.
+""".strip()),
+ ],
+ help=self.__doc__, plugin=plugin)
+
+ def _run(self, hooke, inqueue, outqueue, params):
+ data = params['curve'].data[int(params['block'])] # HACK, int() should be handled by ui
+ # HACK? rely on params['curve'] being bound to the local hooke
+ # playlist (i.e. not a copy, as you would get by passing a
+ # curve through the queue). Ugh. Stupid queues. As an
+ # alternative, we could pass lookup information through the
+ # queue.
+ new = Data((data.shape[0], data.shape[1]+1), dtype=data.dtype)
+ new.info = copy.deepcopy(data.info)
+ new[:,:-1] = data
+ new.info['columns'].append('deflection (N)')
+ d_data = data[:,data.info['columns'].index('surface adjusted deflection (m)')]
+ new[:,-1] = d_data * data.info['spring constant (N/m)']
+ params['curve'].data[int(params['block'])] = new # HACK, int() should be handled by ui
+
projects / hooke.git / commitdiff