Mercurial Hosting > traffic-intelligence

comparison python/moving.py @ 877:d1ff6917d082

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
added savitzky golay filter for accelerations
author Nicolas Saunier <nicolas.saunier@polymtl.ca>
date Wed, 08 Mar 2017 17:46:28 -0500
parents eb2f8ce2b39d
children 000555430b28
comparison
equal deleted inserted replaced
876:c7e72d758049 877:d1ff6917d082
7 from math import sqrt, atan2, cos, sin 7 from math import sqrt, atan2, cos, sin
8 from numpy import median, array, zeros, hypot, NaN, std, floor, float32 8 from numpy import median, array, zeros, hypot, NaN, std, floor, float32
9 from matplotlib.pyplot import plot, text 9 from matplotlib.pyplot import plot, text
10 from scipy.stats import scoreatpercentile 10 from scipy.stats import scoreatpercentile
11 from scipy.spatial.distance import cdist 11 from scipy.spatial.distance import cdist
12 from scipy.signal import savgol_filter
12 13
13 try: 14 try:
14 from shapely.geometry import Polygon, Point as shapelyPoint 15 from shapely.geometry import Polygon, Point as shapelyPoint
15 from shapely.prepared import prep, PreparedGeometry 16 from shapely.prepared import prep, PreparedGeometry
16 shapelyAvailable = True 17 shapelyAvailable = True
802 diff.addPosition(self[i]-self[i-1]) 803 diff.addPosition(self[i]-self[i-1])
803 if doubleLastPosition: 804 if doubleLastPosition:
804 diff.addPosition(diff[-1]) 805 diff.addPosition(diff[-1])
805 return diff 806 return diff
806 807
807 def differentiateSG(self, window_length, polyorder, deriv=0, delta=1.0, axis=-1, mode='interp', cval=0.0, removeBothEnds = 2): 808 def differentiateSG(self, window_length, polyorder, deriv=0, delta=1.0, axis=-1, mode='interp', cval=0.0, nInstantsIgnoredAtEnds = 2):
808 '''Differentiates the trajectory using the Savitsky Golay filter 809 '''Differentiates the trajectory using the Savitsky Golay filter
809 810
810 window_length : The length of the filter window (i.e. the number of coefficients). window_length must be a positive odd integer. 811 window_length : The length of the filter window (i.e. the number of coefficients). window_length must be a positive odd integer.
811 polyorder : The order of the polynomial used to fit the samples. polyorder must be less than window_length. 812 polyorder : The order of the polynomial used to fit the samples. polyorder must be less than window_length.
812 deriv : The order of the derivative to compute. This must be a nonnegative integer. The default is 0, which means to filter the data without differentiating. 813 deriv : The order of the derivative to compute. This must be a nonnegative integer. The default is 0, which means to filter the data without differentiating.
813 delta : The spacing of the samples to which the filter will be applied. This is only used if deriv > 0. Default is 1.0. 814 delta : The spacing of the samples to which the filter will be applied. This is only used if deriv > 0. Default is 1.0.
814 axis : The axis of the array x along which the filter is to be applied. Default is -1. 815 axis : The axis of the array x along which the filter is to be applied. Default is -1.
815 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. 816 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.
816 cval : Value to fill past the edges of the input if mode is constant. Default is 0.0.''' 817 cval : Value to fill past the edges of the input if mode is constant. Default is 0.0.
817 from scipy.signal import savgol_filter 818
818 819 https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.savgol_filter.html#scipy.signal.savgol_filter'''
819 if removeBothEnds >=1: 820 if removeBothEnds >=1:
820 pos = [self.positions[0][removeBothEnds:-removeBothEnds], 821 pos = [self.positions[0][nInstantsIgnoredAtEnds:-nInstantsIgnoredAtEnds],
821 self.positions[1][removeBothEnds:-removeBothEnds]] 822 self.positions[1][nInstantsIgnoredAtEnds:-nInstantsIgnoredAtEnds]]
822 else: 823 else:
823 pos = self.positions 824 pos = self.positions
824 filtered = savgol_filter(pos, window_length, polyorder, deriv, delta, axis, mode, cval) 825 filtered = savgol_filter(pos, window_length, polyorder, deriv, delta, axis, mode, cval)
825 return Trajectory(filtered) 826 return Trajectory(filtered)
826 827
827 def norm(self): 828 def norm(self):
828 '''Returns the list of the norms at each instant''' 829 '''Returns the list of the norms at each instant'''
829 # def add(x, y): return x+y
830 # sq = map(add, [x*x for x in self.positions[0]], [y*y for y in self.positions[1]])
831 # return sqrt(sq)
832 return hypot(self.positions[0], self.positions[1]) 830 return hypot(self.positions[0], self.positions[1])
833
834 # def cumulatedDisplacement(self):
835 # 'Returns the sum of the distances between each successive point'
836 # displacement = 0
837 # for i in xrange(self.length()-1):
838 # displacement += Point.distanceNorm2(self.__getitem__(i),self.__getitem__(i+1))
839 # return displacement
840 831
841 def computeCumulativeDistances(self): 832 def computeCumulativeDistances(self):
842 '''Computes the distance from each point to the next and the cumulative distance up to the point 833 '''Computes the distance from each point to the next and the cumulative distance up to the point
843 Can be accessed through getDistance(idx) and getCumulativeDistance(idx)''' 834 Can be accessed through getDistance(idx) and getCumulativeDistance(idx)'''
844 self.distances = [] 835 self.distances = []
1295 n = min(nInstantsIgnoredAtEnds, int(floor(self.length()/2.))) 1286 n = min(nInstantsIgnoredAtEnds, int(floor(self.length()/2.)))
1296 return speeds[n:-n] 1287 return speeds[n:-n]
1297 else: 1288 else:
1298 return speeds 1289 return speeds
1299 1290
1291 def getAccelerations(self, window_length, polyorder, delta=1.0, axis=-1, mode='interp', cval=0.0, speeds = None, nInstantsIgnoredAtEnds = 0):
1292 '''Returns the 1-D acceleration from the 1-D speeds
1293 Caution about previously filtered data'''
1294 if speeds is None:
1295 speeds = self.getSpeeds(nInstantsIgnoredAtEnds)
1296 return savgol_filter(speeds, window_length, polyorder, 1, delta, axis, mode, cval)
1297
1300 def getSpeedIndicator(self): 1298 def getSpeedIndicator(self):
1301 from indicators import SeverityIndicator 1299 from indicators import SeverityIndicator
1302 return SeverityIndicator('Speed', {t:self.getVelocityAtInstant(t).norm2() for t in self.getTimeInterval()}) 1300 return SeverityIndicator('Speed', {t:self.getVelocityAtInstant(t).norm2() for t in self.getTimeInterval()})
1303 1301
1304 def getPositionAt(self, i): 1302 def getPositionAt(self, i):
1338 self.positions.plotOnWorldImage(nPixelsPerUnitDistance, options, withOrigin, timeStep, None, **kwargs) 1336 self.positions.plotOnWorldImage(nPixelsPerUnitDistance, options, withOrigin, timeStep, None, **kwargs)
1339 1337
1340 def play(self, videoFilename, homography = None, undistort = False, intrinsicCameraMatrix = None, distortionCoefficients = None, undistortedImageMultiplication = 1.): 1338 def play(self, videoFilename, homography = None, undistort = False, intrinsicCameraMatrix = None, distortionCoefficients = None, undistortedImageMultiplication = 1.):
1341 cvutils.displayTrajectories(videoFilename, [self], homography = homography, firstFrameNum = self.getFirstInstant(), lastFrameNumArg = self.getLastInstant(), undistort = undistort, intrinsicCameraMatrix = intrinsicCameraMatrix, distortionCoefficients = distortionCoefficients, undistortedImageMultiplication = undistortedImageMultiplication) 1339 cvutils.displayTrajectories(videoFilename, [self], homography = homography, firstFrameNum = self.getFirstInstant(), lastFrameNumArg = self.getLastInstant(), undistort = undistort, intrinsicCameraMatrix = intrinsicCameraMatrix, distortionCoefficients = distortionCoefficients, undistortedImageMultiplication = undistortedImageMultiplication)
1342 1340
1343 def speedDiagnostics(self, framerate = 1., display = False): 1341 def speedDiagnostics(self, framerate = 1., display = False, nInstantsIgnoredAtEnds=0):
1344 speeds = framerate*self.getSpeeds() 1342 speeds = framerate*self.getSpeeds(nInstantsIgnoredAtEnds)
1345 coef = utils.linearRegression(range(len(speeds)), speeds) 1343 coef = utils.linearRegression(range(len(speeds)), speeds)
1346 print('min/5th perc speed: {} / {}\nspeed diff: {}\nspeed stdev: {}\nregression: {}'.format(min(speeds), scoreatpercentile(speeds, 5), speeds[-2]-speeds[1], std(speeds), coef[0])) 1344 print('min/5th perc speed: {} / {}\nspeed diff: {}\nspeed stdev: {}\nregression: {}'.format(min(speeds), scoreatpercentile(speeds, 5), speeds[-2]-speeds[1], std(speeds), coef[0]))
1347 if display: 1345 if display:
1348 from matplotlib.pyplot import figure, axis 1346 from matplotlib.pyplot import figure, axis
1349 figure(1) 1347 figure(1)
1350 self.plot() 1348 self.plot()
1351 axis('equal') 1349 axis('equal')
1352 figure(2) 1350 figure(2)
1353 plot(list(self.getTimeInterval()), speeds) 1351 plot(list(self.getTimeInterval()), speeds)
1352 figure(3)
1353 plot(list(self.getTimeInterval()), self.getAccelerations(9, 3, speeds = speeds)) # arbitrary parameter
1354 1354
1355 @staticmethod 1355 @staticmethod
1356 def minMaxDistance(obj1, obj2): 1356 def minMaxDistance(obj1, obj2):
1357 '''Computes the min max distance used for feature grouping''' 1357 '''Computes the min max distance used for feature grouping'''
1358 commonTimeInterval = obj1.commonTimeInterval(obj2) 1358 commonTimeInterval = obj1.commonTimeInterval(obj2)