changeset 877:d1ff6917d082

added savitzky golay filter for accelerations
author Nicolas Saunier <nicolas.saunier@polymtl.ca>
date Wed, 08 Mar 2017 17:46:28 -0500
parents c7e72d758049
children 8e8ec4ece66e
files python/moving.py
diffstat 1 files changed, 17 insertions(+), 17 deletions(-) [+]
line wrap: on
line diff
--- a/python/moving.py	Wed Mar 08 17:05:29 2017 -0500
+++ b/python/moving.py	Wed Mar 08 17:46:28 2017 -0500
@@ -9,6 +9,7 @@
 from matplotlib.pyplot import plot, text
 from scipy.stats import scoreatpercentile
 from scipy.spatial.distance import cdist
+from scipy.signal import savgol_filter
 
 try:
     from shapely.geometry import Polygon, Point as shapelyPoint
@@ -804,7 +805,7 @@
             diff.addPosition(diff[-1])
         return diff
 
-    def differentiateSG(self, window_length, polyorder, deriv=0, delta=1.0, axis=-1, mode='interp', cval=0.0, removeBothEnds = 2):
+    def differentiateSG(self, window_length, polyorder, deriv=0, delta=1.0, axis=-1, mode='interp', cval=0.0, nInstantsIgnoredAtEnds = 2):
         '''Differentiates the trajectory using the Savitsky Golay filter
 
         window_length : The length of the filter window (i.e. the number of coefficients). window_length must be a positive odd integer.
@@ -813,12 +814,12 @@
         delta : The spacing of the samples to which the filter will be applied. This is only used if deriv > 0. Default is 1.0.
         axis : The axis of the array x along which the filter is to be applied. Default is -1.
         mode : Must be mirror, constant, nearest, wrap or interp. This determines the type of extension to use for the padded signal to which the filter is applied. When mode is constant, the padding value is given by cval. See the Notes for more details on mirror, constant, wrap, and nearest. When the interp mode is selected (the default), no extension is used. Instead, a degree polyorder polynomial is fit to the last window_length values of the edges, and this polynomial is used to evaluate the last window_length // 2 output values.
-        cval : Value to fill past the edges of the input if mode is constant. Default is 0.0.'''
-        from scipy.signal import savgol_filter
+        cval : Value to fill past the edges of the input if mode is constant. Default is 0.0.
 
+        https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.savgol_filter.html#scipy.signal.savgol_filter'''
         if removeBothEnds >=1:
-            pos = [self.positions[0][removeBothEnds:-removeBothEnds],
-                   self.positions[1][removeBothEnds:-removeBothEnds]]
+            pos = [self.positions[0][nInstantsIgnoredAtEnds:-nInstantsIgnoredAtEnds],
+                   self.positions[1][nInstantsIgnoredAtEnds:-nInstantsIgnoredAtEnds]]
         else:
             pos = self.positions
         filtered = savgol_filter(pos, window_length, polyorder, deriv, delta, axis, mode, cval)
@@ -826,18 +827,8 @@
 
     def norm(self):
         '''Returns the list of the norms at each instant'''
-#        def add(x, y): return x+y
-#        sq = map(add, [x*x for x in self.positions[0]], [y*y for y in self.positions[1]])
-#        return sqrt(sq)
         return hypot(self.positions[0], self.positions[1])
 
-    # def cumulatedDisplacement(self):
-    #     'Returns the sum of the distances between each successive point'
-    #     displacement = 0
-    #     for i in xrange(self.length()-1):
-    #         displacement += Point.distanceNorm2(self.__getitem__(i),self.__getitem__(i+1))
-    #     return displacement
-
     def computeCumulativeDistances(self):
         '''Computes the distance from each point to the next and the cumulative distance up to the point
         Can be accessed through getDistance(idx) and getCumulativeDistance(idx)'''
@@ -1297,6 +1288,13 @@
         else:
             return speeds
 
+    def getAccelerations(self, window_length, polyorder, delta=1.0, axis=-1, mode='interp', cval=0.0, speeds = None, nInstantsIgnoredAtEnds = 0):
+        '''Returns the 1-D acceleration from the 1-D speeds
+        Caution about previously filtered data'''
+        if speeds is None:
+            speeds = self.getSpeeds(nInstantsIgnoredAtEnds)
+        return savgol_filter(speeds, window_length, polyorder, 1, delta, axis, mode, cval)
+        
     def getSpeedIndicator(self):
         from indicators import SeverityIndicator
         return SeverityIndicator('Speed', {t:self.getVelocityAtInstant(t).norm2() for t in self.getTimeInterval()})
@@ -1340,8 +1338,8 @@
     def play(self, videoFilename, homography = None, undistort = False, intrinsicCameraMatrix = None, distortionCoefficients = None, undistortedImageMultiplication = 1.):
         cvutils.displayTrajectories(videoFilename, [self], homography = homography, firstFrameNum = self.getFirstInstant(), lastFrameNumArg = self.getLastInstant(), undistort = undistort, intrinsicCameraMatrix = intrinsicCameraMatrix, distortionCoefficients = distortionCoefficients, undistortedImageMultiplication = undistortedImageMultiplication)
 
-    def speedDiagnostics(self, framerate = 1., display = False):
-        speeds = framerate*self.getSpeeds()
+    def speedDiagnostics(self, framerate = 1., display = False, nInstantsIgnoredAtEnds=0):
+        speeds = framerate*self.getSpeeds(nInstantsIgnoredAtEnds)
         coef = utils.linearRegression(range(len(speeds)), speeds)
         print('min/5th perc speed: {} / {}\nspeed diff: {}\nspeed stdev: {}\nregression: {}'.format(min(speeds), scoreatpercentile(speeds, 5), speeds[-2]-speeds[1], std(speeds), coef[0]))
         if display:
@@ -1351,6 +1349,8 @@
             axis('equal')
             figure(2)
             plot(list(self.getTimeInterval()), speeds)
+            figure(3)
+            plot(list(self.getTimeInterval()), self.getAccelerations(9, 3, speeds = speeds)) # arbitrary parameter
 
     @staticmethod
     def minMaxDistance(obj1, obj2):