#! /usr/bin/env python
'''Libraries for moving objects, trajectories...'''

import utils
import cvutils

from math import sqrt
from numpy import median

try:
    from shapely.geometry import Polygon, Point as shapelyPoint
    from shapely.prepared import prep
    shapelyAvailable = True
except ImportError:
    print('Shapely library could not be loaded')
    shapelyAvailable = False

__metaclass__ = type

class Interval(object):
    '''Generic interval: a subset of real numbers (not iterable)'''
    def __init__(self, first=0, last=-1, revert = False):
        if revert and last<first:
            self.first=last
            self.last=first
        else:
            self.first=first
            self.last=last

    def __str__(self):
        return '[{0}, {1}]'.format(self.first, self.last)

    def __repr__(self):
        return self.__str__()

    def empty(self):
        return self.first > self.last

    def center(self):
        return (self.first+self.last)/2.

    def length(self):
        '''Returns the length of the interval'''
        return float(max(0,self.last-self.first))

    def equal(self, i2):
        return self.first==i2.first and self.last == i2.last

    def getList(self):
        return [self.first, self.last]

    def contains(self, instant):
        return (self.first<=instant and self.last>=instant)

    def inside(self, interval2):
        '''Indicates if the temporal interval of self is comprised in interval2'''
        return (self.first >= interval2.first) and (self.last <= interval2.last)

    @classmethod
    def union(cls, interval1, interval2):
        '''Smallest interval comprising self and interval2'''
        return cls(min(interval1.first, interval2.first), max(interval2.last, interval2.last))
        
    @classmethod
    def intersection(cls, interval1, interval2):
        '''Largest interval comprised in both self and interval2'''
        return cls(max(interval1.first, interval2.first), min(interval1.last, interval2.last))

    def distance(self, interval2):
        if not Interval.intersection(self, interval2).empty():
            return 0
        elif self.first > interval2.last:
            return self.first - interval2.last
        elif self.last < interval2.first:
            return interval2.first - self.last
        else:
            return None


def unionIntervals(intervals):
    'returns the smallest interval containing all intervals'
    inter = intervals[0]
    for i in intervals[1:]:
        inter = Interval.union(inter, i)
    return inter


class TimeInterval(Interval):
    '''Temporal interval: set of instants at fixed time step, between first and last, included
    
    For example: based on frame numbers (hence the modified length method)
    It may be modified directly by setting first and last'''

    def __init__(self, first=0, last=-1):
        super(TimeInterval, self).__init__(first, last, False)

    @staticmethod
    def fromInterval(inter):
        return TimeInterval(inter.first, inter.last)

    def __getitem__(self, i):
        if not self.empty():
            if isinstance(i, int):
                return self.first+i
            else:
                raise TypeError, "Invalid argument type."
            #elif isinstance( key, slice ):

    def __iter__(self):
        self.iterInstantNum = -1
        return self

    def next(self):
        if self.iterInstantNum >= self.length()-1:
            raise StopIteration
        else:
            self.iterInstantNum += 1
            return self[self.iterInstantNum]

    def length(self):
        '''Returns the length of the interval'''
        return float(max(0,self.last-self.first+1))

# class BoundingPolygon:
#     '''Class for a polygon bounding a set of points
#     with methods to create intersection, unions...
#     '''
# We will use the polygon class of Shapely

class STObject(object):
    '''Class for spatio-temporal object, i.e. with temporal and spatial existence 
    (time interval and bounding polygon for positions (e.g. rectangle)).

    It may not mean that the object is defined 
    for all time instants within the time interval'''

    def __init__(self, num = None, timeInterval = None, boundingPolygon = None):
        self.num = num
        self.timeInterval = timeInterval
        self.boundingPolygon = boundingPolygon

    def empty(self):
        return self.timeInterval.empty() or not self.boudingPolygon

    def getNum(self):
        return self.num

    def getFirstInstant(self):
        return self.timeInterval.first

    def getLastInstant(self):
        return self.timeInterval.last

    def getTimeInterval(self):
        return self.timeInterval

    def existsAtInstant(self, t):
        return self.timeInterval.contains(t)

    def commonTimeInterval(self, obj2):
        return TimeInterval.intersection(self.getTimeInterval(), obj2.getTimeInterval())

class Point(object):
    def __init__(self, x, y):
        self.x = x
        self.y = y

    def __str__(self):
        return '(%f,%f)'%(self.x,self.y)

    def __repr__(self):
        return self.__str__()

    def __add__(self, other):
        return Point(self.x+other.x, self.y+other.y)

    def __sub__(self, other):
        return Point(self.x-other.x, self.y-other.y)

    def __neg__(self):
        return Point(-self.x, -self.y)

    def multiply(self, alpha):
        return Point(self.x*alpha, self.y*alpha)

    def plot(self, options = 'o', **kwargs):
        from matplotlib.pylab import plot
        plot([self.x], [self.y], options, **kwargs)

    def norm2Squared(self):
        '''2-norm distance (Euclidean distance)'''
        return self.x**2+self.y**2

    def norm2(self):
        '''2-norm distance (Euclidean distance)'''
        return sqrt(self.norm2Squared())

    def norm1(self):
        return abs(self.x)+abs(self.y)
    
    def normMax(self):
        return max(abs(self.x),abs(self.y))

    def aslist(self):
        return [self.x, self.y]

    def astuple(self):
        return (self.x, self.y)

    def asint(self):
        return Point(int(self.x), int(self.y))

    if shapelyAvailable:
        def asShapely(self):
            return shapelyPoint(self.x, self.y)

    def project(self, homography):
        from numpy.core.multiarray import array
        projected = cvutils.projectArray(homography, array([[self.x], [self.y]]))
        return Point(projected[0], projected[1])

    def inPolygonNoShapely(self, polygon):
        '''Indicates if the point x, y is inside the polygon
        (array of Nx2 coordinates of the polygon vertices)

        taken from http://www.ariel.com.au/a/python-point-int-poly.html

        Use Polygon.contains if Shapely is installed'''

        n = polygon.shape[0];
        counter = 0;

        p1 = polygon[0,:];
        for i in range(n+1):
            p2 = polygon[i % n,:];
            if self.y > min(p1[1],p2[1]):
                if self.y <= max(p1[1],p2[1]):
                    if self.x <= max(p1[0],p2[0]):
                        if p1[1] != p2[1]:
                            xinters = (self.y-p1[1])*(p2[0]-p1[0])/(p2[1]-p1[1])+p1[0];
                        if p1[0] == p2[0] or self.x <= xinters:
                            counter+=1;
            p1=p2
        return (counter%2 == 1);

    @staticmethod
    def dot(p1, p2):
        'Scalar product'
        return p1.x*p2.x+p1.y*p2.y

    @staticmethod
    def cross(p1, p2):
        'Cross product'
        return p1.x*p2.y-p1.y*p2.x

    @staticmethod
    def cosine(p1, p2):
        return Point.dot(p1,p2)/(p1.norm2()*p2.norm2())

    @staticmethod
    def distanceNorm2(p1, p2):
        return (p1-p2).norm2()

    @staticmethod
    def plotAll(points, **kwargs):
        from matplotlib.pyplot import scatter
        scatter([p.x for p in points],[p.y for p in points], **kwargs)

    def similarOrientation(self, refDirection, cosineThreshold):
        'Indicates whether the cosine of the vector and refDirection is smaller than cosineThreshold'
        return Point.cosine(self, refDirection) >= cosineThreshold

    @staticmethod
    def timeToCollision(p1, p2, v1, v2, collisionThreshold):
        '''Computes exact time to collision with a distance threshold
        The unknown of the equation is the time to reach the intersection
        between the relative trajectory of one road user
        and the circle of radius collisionThreshold around the other road user'''
        from math import sqrt
        dv = v1-v2
        dp = p1-p2
        a = dv.norm2Squared()#(v1.x-v2.x)**2 + (v1.y-v2.y)**2
        b = 2*Point.dot(dv, dp)#2 * ((p1.x-p2.x) * (v1.x-v2.x) + (p1.y-p2.y) * (v1.y-v2.y))
        c = dp.norm2Squared() - collisionThreshold**2#(p1.x-p2.x)**2 + (p1.y-p2.y)**2 - collisionThreshold**2

        delta = b**2 - 4*a*c
        if delta >= 0:
            deltaRoot = sqrt(delta)
            ttc1 = (-b + deltaRoot)/(2*a)
            ttc2 = (-b - deltaRoot)/(2*a)
            if ttc1 >= 0 and ttc2 >= 0:
                ttc = min(ttc1,ttc2)
            elif ttc1 >= 0:
                ttc = ttc1
            elif ttc2 >= 0:
                ttc = ttc2
            else: # ttc1 < 0 and ttc2 < 0:
                ttc = None
        else:
            ttc = None
        return ttc


if shapelyAvailable:
    def pointsInPolygon(points, polygon):
        '''Optimized tests of a series of points within (Shapely) polygon '''
        prepared_polygon = prep(polygon)
        return filter(prepared_polygon.contains, points)

# Functions for coordinate transformation
# From Paul St-Aubin's PVA tools
def subsec_spline_dist(splines):
    ''' Prepare list of spline subsegments from a spline list. 
    
    Output:
    =======
    ss_spline_d[spline #][mode][station]
    
    where:
        mode=0: incremental distance
        mode=1: cumulative distance
        mode=2: cumulative distance with trailing distance
    '''

    from numpy import zeros
    ss_spline_d = []
    #Prepare subsegment distances
    for spline in range(len(splines)):
        ss_spline_d.append([[],[],[]])
        ss_spline_d[spline][0] = zeros(len(splines[spline])-1)  #Incremental distance
        ss_spline_d[spline][1] = zeros(len(splines[spline])-1)  #Cumulative distance
        ss_spline_d[spline][2] = zeros(len(splines[spline]))  #Cumulative distance with trailing distance
        for spline_p in range(len(splines[spline])):
            if spline_p > (len(splines[spline]) - 2):
                break
            ss_spline_d[spline][0][spline_p] = utils.pointDistanceL2(splines[spline][spline_p][0],splines[spline][spline_p][1],splines[spline][(spline_p+1)][0],splines[spline][(spline_p+1)][1])
            ss_spline_d[spline][1][spline_p] = sum(ss_spline_d[spline][0][0:spline_p])
            ss_spline_d[spline][2][spline_p] = ss_spline_d[spline][1][spline_p]#sum(ss_spline_d[spline][0][0:spline_p])
    
    ss_spline_d[spline][2][-1] = ss_spline_d[spline][2][-2] + ss_spline_d[spline][0][-1]

    return ss_spline_d

def ppldb2p(qx,qy, p0x,p0y, p1x,p1y):
    ''' Point-projection (Q) on line defined by 2 points (P0,P1). 
        http://cs.nyu.edu/~yap/classes/visual/03s/hw/h2/math.pdf
        '''
    if(p0x == p1x and p0y == p1y):
        return False,False
    try:
        #Approximate slope singularity by giving some slope roundoff; account for roundoff error
        if(round(p0x, 10) == round(p1x, 10)):
            p1x += 0.0000000001
        if(round(p0y, 10) == round(p1y, 10)):
            p1y += 0.0000000001            
        #make the calculation
        Y = (-(qx)*(p0y-p1y)-(qy*(p0y-p1y)**2)/(p0x-p1x)+p0x**2*(p0y-p1y)/(p0x-p1x)-p0x*p1x*(p0y-p1y)/(p0x-p1x)-p0y*(p0x-p1x))/(p1x-p0x-(p0y-p1y)**2/(p0x-p1x))
        X = (-Y*(p1y-p0y)+qx*(p1x-p0x)+qy*(p1y-p0y))/(p1x-p0x)
    except ZeroDivisionError:
        print('Error: Division by zero in ppldb2p. Please report this error with the full traceback:')
        print('qx={0}, qy={1}, p0x={2}, p0y={3}, p1x={4}, p1y={5}...'.format(qx, qy, p0x, p0y, p1x, p1y))
        import pdb; pdb.set_trace()  
    return X,Y

def getSYfromXY(qx, qy, splines, spline_assumption_threshold=0.5, mode=0):
    ''' Snap a point (coordinates qx, qy) to it's nearest subsegment of it's nearest spline (from the list splines).
        
        To force snapping to a single spline, pass the spline alone through splines (e.g. splines=[splines[splineNum]]).
        
        Output:
        =======
        [spline index, 
        subsegment leading point index, 
        snapped x coordinate, 
        snapped y coordinate, 
        subsegment distance, 
        spline distance,
        orthogonal point offset]
    '''
    
    #(buckle in, it gets ugly from here on out)
    ss_spline_d = subsec_spline_dist(splines)
    
    temp_dist_min = float('inf')
    #For each spline
    for spline in range(len(splines)):
        #For each spline point
        for spline_p in range(len(splines[spline])):
            if(spline_p > (len(splines[spline]) - 2)):
                break
            #Get point-intersection distance 
            X,Y = ppldb2p(qx,qy,splines[spline][spline_p][0],splines[spline][spline_p][1],splines[spline][spline_p+1][0],splines[spline][spline_p+1][1])
            if(X == False and Y == False):
                print('Error: Spline {0}, segment {1} has identical bounds and therefore is not a vector. Projection cannot continue.'.format(spline, spline_p))
                return [False,False,False,False,False,False,False]
            #Check to see if point is not contained by subspline
            if ss_spline_d[spline][0][spline_p]*1.05 < max(utils.pointDistanceL2(splines[spline][spline_p][0],splines[spline][spline_p][1],X,Y),utils.pointDistanceL2(splines[spline][spline_p+1][0],splines[spline][spline_p+1][1],X,Y)): continue          
            
            #Ok
            temp_dist = utils.pointDistanceL2(qx,qy,X,Y)
            if(temp_dist < temp_dist_min):
                temp_dist_min             = temp_dist
                snappedSpline             = spline
                snappedSplineLeadingPoint = spline_p
                snapped_x                 = X
                snapped_y                 = Y
            #Jump loop if significantly close
            if(temp_dist < spline_assumption_threshold): break

    try:
        #Get sub-segment distance
        subsegmentDistance = utils.pointDistanceL2(splines[snappedSpline][snappedSplineLeadingPoint][0],splines[snappedSpline][snappedSplineLeadingPoint][1],snapped_x,snapped_y)
        #Get total segment distance
        splineDistanceS = ss_spline_d[snappedSpline][1][snappedSplineLeadingPoint] + subsegmentDistance
        #Get offset distance    
        offsetY = utils.pointDistanceL2(qx,qy, snapped_x,snapped_y) 
        
        if(mode):        
            direction = getOrientation([splines[snappedSpline][snappedSplineLeadingPoint][0],splines[snappedSpline][snappedSplineLeadingPoint][1]] , [splines[snappedSpline][snappedSplineLeadingPoint+1][0],splines[snappedSpline][snappedSplineLeadingPoint+1][1]] , [qx,qy])       
            if(direction == 'left'): offsetY = -offsetY       
        
        return [snappedSpline, snappedSplineLeadingPoint, snapped_x, snapped_y, subsegmentDistance, splineDistanceS, offsetY]
    except:
        return [False,False,False,False,False,False,False]
    
    
def getXYfromSY(s, y, splineNum, splines, mode = 0):
    ''' Find X,Y coordinate from S,Y data. 
    if mode = 0 : return Snapped X,Y
    if mode !=0 : return Real X,Y
    ''' 
    
    #(buckle in, it gets ugly from here on out)
    ss_spline_d = subsec_spline_dist(splines)
    
    #Find subsegment
    snapped_x = None
    snapped_y = None
    for spline_ss_index in range(len(ss_spline_d[splineNum][1])):
        if(s < ss_spline_d[splineNum][1][spline_ss_index]):
            ss_value = s - ss_spline_d[splineNum][1][spline_ss_index-1]
            #Get normal vector and then snap
            vector_l_x = (splines[splineNum][spline_ss_index][0] - splines[splineNum][spline_ss_index-1][0])
            vector_l_y = (splines[splineNum][spline_ss_index][1] - splines[splineNum][spline_ss_index-1][1])
            magnitude  = sqrt(vector_l_x**2 + vector_l_y**2)
            n_vector_x = vector_l_x/magnitude
            n_vector_y = vector_l_y/magnitude
            snapped_x  = splines[splineNum][spline_ss_index-1][0] + ss_value*n_vector_x
            snapped_y  = splines[splineNum][spline_ss_index-1][1] + ss_value*n_vector_y

            #Real values (including orthogonal projection of y))
            real_x = snapped_x - y*n_vector_y 
            real_y = snapped_y + y*n_vector_x            
            break
    
    if mode == 0 or (not snapped_x):
        if(not snapped_x):
            snapped_x = splines[splineNum][-1][0]
            snapped_y = splines[splineNum][-1][1]                
        return [snapped_x,snapped_y]
    else:
        return [real_x,real_y]


class NormAngle(object):
    '''Alternate encoding of a point, by its norm and orientation'''

    def __init__(self, norm, angle):
        self.norm = norm
        self.angle = angle
    
    @staticmethod
    def fromPoint(p):
        from math import atan2
        norm = p.norm2()
        if norm > 0:
            angle = atan2(p.y, p.x)
        else:
            angle = 0.
        return NormAngle(norm, angle)

    def __add__(self, other):
        'a norm cannot become negative'
        return NormAngle(max(self.norm+other.norm, 0), self.angle+other.angle)

    def getPoint(self):
        from math import cos, sin
        return Point(self.norm*cos(self.angle), self.norm*sin(self.angle))


def predictPositionNoLimit(nTimeSteps, initialPosition, initialVelocity, initialAcceleration = Point(0,0)):
    '''Predicts the position in nTimeSteps at constant speed/acceleration'''
    return initialVelocity + initialAcceleration.multiply(nTimeSteps),initialPosition+initialVelocity.multiply(nTimeSteps) + initialAcceleration.multiply(nTimeSteps**2*0.5)

def predictPosition(position, speedOrientation, control, maxSpeed = None):
    '''Predicts the position (moving.Point) at the next time step with given control input (deltaSpeed, deltaTheta)
    speedOrientation is the other encoding of velocity, (speed, orientation)
    speedOrientation and control are NormAngle'''
    predictedSpeedTheta = speedOrientation+control
    if maxSpeed:
         predictedSpeedTheta.norm = min(predictedSpeedTheta.norm, maxSpeed)
    predictedPosition = position+predictedSpeedTheta.getPoint()
    return predictedPosition, predictedSpeedTheta


class FlowVector(object):
    '''Class to represent 4-D flow vectors,
    ie a position and a velocity'''
    def __init__(self, position, velocity):
        'position and velocity should be Point instances'
        self.position = position
        self.velocity = velocity

    def __add__(self, other):
        return FlowVector(self.position+other.position, self.velocity+other.velocity)

    def multiply(self, alpha):
        return FlowVector(self.position.multiply(alpha), self.velocity.multiply(alpha))

    def plot(self, options = '', **kwargs):
        from matplotlib.pylab import plot
        plot([self.position.x, self.position.x+self.velocity.x], [self.position.y, self.position.y+self.velocity.y], options, **kwargs)
        self.position.plot(options+'x', **kwargs)
    
    @staticmethod
    def similar(f1, f2, maxDistance2, maxDeltavelocity2):
        return (f1.position-f2.position).norm2Squared()<maxDistance2 and (f1.velocity-f2.velocity).norm2Squared()<maxDeltavelocity2

def intersection(p1, p2, p3, p4):
    ''' Intersection point (x,y) of lines formed by the vectors p1-p2 and p3-p4
        http://paulbourke.net/geometry/pointlineplane/'''
    dp12 = p2-p1
    dp34 = p4-p3
    #det = (p4.y-p3.y)*(p2.x-p1.x)-(p4.x-p3.x)*(p2.y-p1.y)
    det = dp34.y*dp12.x-dp34.x*dp12.y
    if det == 0:
        return None
    else:
        ua = (dp34.x*(p1.y-p3.y)-dp34.y*(p1.x-p3.x))/det
        return p1+dp12.multiply(ua)

# def intersection(p1, p2, dp1, dp2):
#     '''Returns the intersection point between the two lines 
#     defined by the respective vectors (dp) and origin points (p)'''
#     from numpy import matrix
#     from numpy.linalg import linalg
#     A = matrix([[dp1.y, -dp1.x],
#                 [dp2.y, -dp2.x]])
#     B = matrix([[dp1.y*p1.x-dp1.x*p1.y],
#                 [dp2.y*p2.x-dp2.x*p2.y]])
    
#     if linalg.det(A) == 0:
#         return None
#     else:
#         intersection = linalg.solve(A,B)
#         return Point(intersection[0,0], intersection[1,0])

def segmentIntersection(p1, p2, p3, p4):
    '''Returns the intersecting point of the segments [p1, p2] and [p3, p4], None otherwise'''

    if (Interval.intersection(Interval(p1.x,p2.x,True), Interval(p3.x,p4.x,True)).empty()) or (Interval.intersection(Interval(p1.y,p2.y,True), Interval(p3.y,p4.y,True)).empty()):
        return None
    else:
        inter = intersection(p1, p2, p3, p4)
        if (inter != None 
            and utils.inBetween(p1.x, p2.x, inter.x)
            and utils.inBetween(p3.x, p4.x, inter.x)
            and utils.inBetween(p1.y, p2.y, inter.y)
            and utils.inBetween(p3.y, p4.y, inter.y)):
            return inter
        else:
            return None

class Trajectory(object):
    '''Class for trajectories: temporal sequence of positions

    The class is iterable'''

    def __init__(self, positions=None):
        if positions != None:
            self.positions = positions
        else:
            self.positions = [[],[]]

    @staticmethod
    def generate(p, v, nPoints):
        t = Trajectory()
        p0 = Point(p.x, p.y)
        t.addPosition(p0)
        for i in xrange(nPoints-1):
            p0 += v
            t.addPosition(p0)
        return t, Trajectory([[v.x]*nPoints, [v.y]*nPoints])

    @staticmethod
    def load(line1, line2):
        return Trajectory([[float(n) for n in line1.split(' ')],
                           [float(n) for n in line2.split(' ')]])

    @staticmethod
    def fromPointList(points):
        t = Trajectory()
        for p in points:
            t.addPosition(p)
        return t

    def __len__(self):
        return len(self.positions[0])

    def length(self):
        return self.__len__()

    def __getitem__(self, i):
        if isinstance(i, int):
            return Point(self.positions[0][i], self.positions[1][i])
        else:
            raise TypeError, "Invalid argument type."
            #elif isinstance( key, slice ):

    def __str__(self):
        return ' '.join([self.__getitem__(i).__str__() for i in xrange(self.length())])

    def __repr__(self):
        return self.__str__()


    def __iter__(self):
        self.iterInstantNum = 0
        return self

    def next(self):
        if self.iterInstantNum >= self.length():
            raise StopIteration
        else:
            self.iterInstantNum += 1
            return self[self.iterInstantNum-1]

    def setPositionXY(self, i, x, y):
        if i < self.__len__():
            self.positions[0][i] = x
            self.positions[1][i] = y

    def setPosition(self, i, p):
        self.setPositionXY(i, p.x, p.y)

    def addPositionXY(self, x, y):
        self.positions[0].append(x)
        self.positions[1].append(y)

    def addPosition(self, p):
        self.addPositionXY(p.x, p.y)

    def duplicateLastPosition(self):
        self.positions[0].append(self.positions[0][-1])
        self.positions[1].append(self.positions[1][-1])

    @staticmethod
    def _plot(positions, options = '', withOrigin = False, lastCoordinate = None, timeStep = 1, **kwargs):
        from matplotlib.pylab import plot
        if lastCoordinate == None:
            plot(positions[0][::timeStep], positions[1][::timeStep], options, **kwargs)
        elif 0 <= lastCoordinate <= len(positions[0]):
            plot(positions[0][:lastCoordinate:timeStep], positions[1][:lastCoordinate:timeStep], options, **kwargs)
        if withOrigin:
            plot([positions[0][0]], [positions[1][0]], 'ro', **kwargs)

    def project(self, homography):
        return Trajectory(cvutils.projectTrajectory(homography, self.positions))

    def plot(self, options = '', withOrigin = False, timeStep = 1, **kwargs):
        Trajectory._plot(self.positions, options, withOrigin, None, timeStep, **kwargs)

    def plotAt(self, lastCoordinate, options = '', withOrigin = False, timeStep = 1, **kwargs):
        Trajectory._plot(self.positions, options, withOrigin, lastCoordinate, timeStep, **kwargs)

    def plotOnWorldImage(self, nPixelsPerUnitDistance, imageHeight, options = '', withOrigin = False, timeStep = 1, **kwargs):
        from matplotlib.pylab import plot
        imgPositions = [[x*nPixelsPerUnitDistance for x in self.positions[0]],
                        [-x*nPixelsPerUnitDistance+imageHeight for x in self.positions[1]]]
        Trajectory._plot(imgPositions, options, withOrigin, timeStep, **kwargs)

    def getXCoordinates(self):
        return self.positions[0]

    def getYCoordinates(self):
        return self.positions[1]

    def asArray(self):
        from numpy.core.multiarray import array
        return array(self.positions)
    
    def xBounds(self):
        # look for function that does min and max in one pass
        return Interval(min(self.getXCoordinates()), max(self.getXCoordinates()))
    
    def yBounds(self):
        # look for function that does min and max in one pass
        return Interval(min(self.getYCoordinates()), max(self.getYCoordinates()))
    
    def add(self, traj2):
        '''Returns a new trajectory of the same length'''
        if self.length() != traj2.length():
            print 'Trajectories of different lengths'
            return None
        else:
            return Trajectory([[a+b for a,b in zip(self.getXCoordinates(),traj2.getXCoordinates())],
                               [a+b for a,b in zip(self.getYCoordinates(),traj2.getYCoordinates())]])

    def subtract(self, traj2):
        '''Returns a new trajectory of the same length'''
        if self.length() != traj2.length():
            print 'Trajectories of different lengths'
            return None
        else:
            return Trajectory([[a-b for a,b in zip(self.getXCoordinates(),traj2.getXCoordinates())],
                               [a-b for a,b in zip(self.getYCoordinates(),traj2.getYCoordinates())]])

    def differentiate(self, doubleLastPosition = False):
        diff = Trajectory()
        for i in xrange(1, self.length()):
            diff.addPosition(self[i]-self[i-1])
        if doubleLastPosition:
            diff.addPosition(diff[-1])
        return diff

    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)
        from numpy import hypot
        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 similarOrientation(self, refDirection, cosineThreshold, minProportion = 0.5):
        '''Indicates whether the minProportion (<=1.) (eg half) of the trajectory elements (vectors for velocity) 
        have a cosine with refDirection is smaller than cosineThreshold'''
        count = 0
        lengthThreshold = float(self.length())*minProportion
        for p in self:
            if p.similarOrientation(refDirection, cosineThreshold):
                count += 1
            if count > lengthThreshold:
                return True
        return False

    def wiggliness(self):
        return self.cumulatedDisplacement()/float(Point.distanceNorm2(self.__getitem__(0),self.__getitem__(self.length()-1)))

    def getIntersections(self, p1, p2):
        '''Returns a list of the indices at which the trajectory 
        intersects with the segment of extremities p1 and p2 
        the list is empty if there is no crossing'''
        indices = []

        for i in xrange(self.length()-1):
            q1=self.__getitem__(i)
            q2=self.__getitem__(i+1)
            p = utils.segmentIntersection(q1, q2, p1, p2)
            if p:
                if q1.x != q2.x:
                    ratio = (p.x-q1.x)/(q2.x-q1.x)
                elif q1.y != q2.y:
                    ratio = (p.y-q1.y)/(q2.y-q1.y)
                else:
                    ratio = 0
                indices.append(i+ratio)
        return indices

    def getTrajectoryInInterval(self, inter):
        if inter.first >=0 and inter.last<= self.length():
            return Trajectory([self.positions[0][inter.first:inter.last],
                               self.positions[1][inter.first:inter.last]])
        else:
            return None
    
    def getTrajectoryInPolygonNoShapely(self, polygon):
        '''Returns the trajectory built with the set of points inside the polygon
        (array of Nx2 coordinates of the polygon vertices)'''
        traj = Trajectory()
        for p in self:
            if p.inPolygonNoShapely(polygon):
                traj.addPosition(p)
        return traj

    if shapelyAvailable:
        def getTrajectoryInPolygon(self, polygon):
            '''Returns the trajectory built with the set of points inside the (shapely) polygon'''
            traj = Trajectory()
            points = [p.asShapely() for p in self]
            for p in pointsInPolygon(points, polygon):
                    traj.addPositionXY(p.x, p.y)
            return traj
    
    @staticmethod
    def lcss(t1, t2, lcss):
        return lcss.compute(t1, t2)

class CurvilinearTrajectory(Trajectory):
    '''Sub class of trajectory for trajectories with curvilinear coordinates and lane assignements
    longitudinal coordinate is stored as first coordinate (exterior name S)
    lateral coordiante is stored as second coordinate'''

    def __init__(self, S = None, Y = None, lanes = None):
        if S == None or Y == None or len(S) != len(Y):
            self.positions = [[],[]]
            if S != None and Y != None and len(S) != len(Y):
                print("S and Y coordinates of different lengths\nInitializing to empty lists")
        else:
            self.positions = [S,Y]
        if lanes == None or len(lanes) != self.length():
            self.lanes = []
        else:
            self.lanes = lanes
        
    def __getitem__(self,i): 
        if isinstance(i, int):
            return [self.positions[0][i], self.positions[1][i], self.lanes[i]]
        else:
            raise TypeError, "Invalid argument type."
            #elif isinstance( key, slice ):

    def getSCoordinates(self):
        return self.getXCoordinates()
    
    def getLanes(self):
        return self.lanes

    def addPositionSYL(self, s, y, lane):
        self.addPositionXY(s,y)
        self.lanes.append(lane)

    def addPosition(self, p):
        'Adds position in the point format for curvilinear of list with 3 values'
        self.addPositionSYL(p[0], p[1], p[2])

    def setPosition(self, i, s, y, lane):
        self.setPositionXY(i, s, y)
        if i < self.__len__():
            self.lanes[i] = lane

    def differentiate(self, doubleLastPosition = False):
        diff = CurvilinearTrajectory()
        p1 = self[0]
        for i in xrange(1, self.length()):
            p2 = self[i]
            diff.addPositionSYL(p2[0]-p1[0], p2[1]-p1[1], p1[2])
            p1=p2
        if doubleLastPosition:
            diff.addPosition(diff[-1])
        return diff

    def getIntersections(self, S1, lane = None):
        '''Returns a list of the indices at which the trajectory 
        goes past the curvilinear coordinate S1
        (in provided lane if lane != None)
        the list is empty if there is no crossing'''
        indices = []
        for i in xrange(self.length()-1):
            q1=self.__getitem__(i)
            q2=self.__getitem__(i+1)
            if q1[0] <= S1 < q2[0] and (lane == None or (self.lanes[i] == lane and self.lanes[i+1] == lane)):
                indices.append(i+(S1-q1[0])/(q2[0]-q1[0]))
        return indices

##################
# Moving Objects
##################

userTypeNames = ['unknown',
                 'car',
                 'pedestrian',
                 'motorcycle',
                 'bicycle',
                 'bus',
                 'truck']

userType2Num = utils.inverseEnumeration(userTypeNames)

class MovingObject(STObject):
    '''Class for moving objects: a spatio-temporal object 
    with a trajectory and a geometry (constant volume over time) and a usertype (e.g. road user) coded as a number (see 
    '''

    def __init__(self, num = None, timeInterval = None, positions = None, velocities = None, geometry = None, userType = userType2Num['unknown']):
        super(MovingObject, self).__init__(num, timeInterval)
        self.positions = positions
        self.velocities = velocities
        self.geometry = geometry
        self.userType = userType
        self.features = []
        # compute bounding polygon from trajectory
        
    @staticmethod
    def generate(p, v, timeInterval):
        positions, velocities = Trajectory.generate(p, v, int(timeInterval.length())) 
        return MovingObject(timeInterval = timeInterval, positions = positions, velocities = velocities)

    def getObjectInTimeInterval(self, inter):
        '''Returns a new object extracted from self,
        restricted to time interval inter'''
        intersection = TimeInterval.intersection(inter, self.getTimeInterval())
        if not intersection.empty():
            trajectoryInterval = TimeInterval(intersection.first-self.getFirstInstant(), intersection.last-self.getFirstInstant())
            obj = MovingObject(self.num, intersection, self.positions.getTrajectoryInInterval(trajectoryInterval), self.geometry, self.userType)
            if self.velocities:
                obj.velocities = self.velocities.getTrajectoryInInterval(trajectoryInterval)
            return obj
        else:
            print 'The object does not exist at '+str(inter)
            return None

    def length(self):
        return self.timeInterval.length()

    def getPositions(self):
        return self.positions

    def getVelocities(self):
        return self.velocities

    def getUserType(self):
        return self.userType

    def getCurvilinearPositions(self):
        if hasattr(self, 'curvilinearPositions'):
            return self.curvilinearPositions
        else:
            return None

    def setUserType(self, userType):
        self.userType = userType

    def setFeatures(self, features):
        self.features = [features[i] for i in self.featureNumbers]

    def getSpeeds(self):
        return self.getVelocities().norm()

    def getPositionAt(self, i):
        return self.positions[i]

    def getVelocityAt(self, i):
        return self.velocities[i]

    def getPositionAtInstant(self, i):
        return self.positions[i-self.getFirstInstant()]

    def getVelocityAtInstant(self, i):
        return self.velocities[i-self.getFirstInstant()]

    def getXCoordinates(self):
        return self.positions.getXCoordinates()
    
    def getYCoordinates(self):
        return self.positions.getYCoordinates()
    
    def plot(self, options = '', withOrigin = False, timeStep = 1, withFeatures = False, **kwargs):
        if withFeatures:
            for f in self.features:
                f.positions.plot('r', True, timeStep, **kwargs)
            self.positions.plot('bx-', True, timeStep, **kwargs)
        else:
            self.positions.plot(options, withOrigin, timeStep, **kwargs)

    def plotOnWorldImage(self, nPixelsPerUnitDistance, imageHeight, options = '', withOrigin = False, timeStep = 1, **kwargs):
        self.positions.plotOnWorldImage(nPixelsPerUnitDistance, imageHeight, options, withOrigin, timeStep, **kwargs)

    def play(self, videoFilename, homography = None):
        cvutils.displayTrajectories(videoFilename, [self], homography, self.getFirstInstant(), self.getLastInstant())

    def speedDiagnostics(self, framerate = 1., display = False):
        from numpy import std
        from scipy.stats import scoreatpercentile
        speeds = framerate*self.getSpeeds()
        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:
            from matplotlib.pyplot import figure, plot, axis
            figure(1)
            self.plot()
            axis('equal')
            figure(2)
            plot(list(self.getTimeInterval()), speeds)

    @staticmethod
    def distances(obj1, obj2, instant1, _instant2 = None):
        from scipy.spatial.distance import cdist
        if _instant2 == None:
            instant2 = instant1
        else:
            instant2 = _instant2
        positions1 = [f.getPositionAtInstant(instant1).astuple() for f in obj1.features if f.existsAtInstant(instant1)]
        positions2 = [f.getPositionAtInstant(instant2).astuple() for f in obj2.features if f.existsAtInstant(instant2)]
        return cdist(positions1, positions2, metric = 'euclidean')
        
    @staticmethod
    def minDistance(obj1, obj2, instant1, instant2 = None):
        return MovingObject.distances(obj1, obj2, instant1, instant2).min()

    @staticmethod
    def maxDistance(obj1, obj2, instant, instant2 = None):
        return MovingObject.distances(obj1, obj2, instant1, instant2).max()

    def maxSize(self):
        '''Returns the max distance between features
        at instant there are the most features'''
        if hasattr(self, 'features'):
            nFeatures = -1
            tMaxFeatures = 0
            for t in self.getTimeInterval():
                n = len([f for f in self.features if f.existsAtInstant(t)])
                if n > nFeatures:
                    nFeatures = n
                    tMaxFeatures = t
            return MovingObject.maxDistance(self, self, tMaxFeatures)
        else:
            print('Load features to compute a maximum size')
            return None
			
    def setRoutes(self, startRouteID, endRouteID):
        self.startRouteID = startRouteID
        self.endRouteID = endRouteID
           
    def getInstantsCrossingLane(self, p1, p2):
        '''Returns the instant(s)
        at which the object passes from one side of the segment to the other
        empty list if there is no crossing'''
        indices = self.positions.getIntersections(p1, p2)
        return [t+self.getFirstInstant() for t in indices]

    def predictPosition(self, instant, nTimeSteps, externalAcceleration = Point(0,0)):
        '''Predicts the position of object at instant+deltaT, 
        at constant speed'''
        return predictPositionNoLimit(nTimeSteps, self.getPositionAtInstant(instant), self.getVelocityAtInstant(instant), externalAcceleration)

    def transformToCurvilinear(objects, alignments, ln_mv_av_win=3, verbose=0):
        ''' Add, for every object position, the class 'moving.CurvilinearTrajectory()'
            (curvilinearPositions instance) which holds information about the
            curvilinear coordinates using alignment metadata.
            From Paul St-Aubin's PVA tools
            ======

            Input:
            ======
            objects      = list of objects supplied by Traffic Intelligence
            alignments   = a list of alignments, where each alignment is a list of
                           points, where each point is a list of coordinates, e.g.
                           alignments[alpha][1][x] is coordinate x for point number
                           1 of the spline that represents alignment alpha.
                           alignments can also be a compatible object that mimics a
                           3-dimensional list using the __getitem__() method as is
                           the case in the PVAT specification.
            ln_mv_av_win = moving average window (in points) in which to smooth
                           lane changes. As per tools_math.cat_mvgavg(), this term
                           is a search *radius* around the center of the window.

            Output: (objects, dropped_traj)
            =======
            objects        = modified list of objects
            dropped_traj   = list of objects or object segments that were
                             truncated. The format is identical to that of objects.
            '''
        if(len(objects) <= 0): return None, None
        if(verbose >= 2):
            print('    -------------')
            print('    Transforming coordinates...')

        lane_readjustments     = 0
        dropped_traj           = [] 
        original_object_length = len(objects)
        reset_objects          = 0

        #if(not isinstance(objects, list)):
        #    objects = [objects]
        #    reset_objects = 1
        #if(not isinstance(objects[0], TrafIntMoving.MovingObject)):
        #    return objects, dropped_traj

        #For each object
        for i in range(len(objects)):
            objects[i].curvilinearPositions = CurvilinearTrajectory([],[],[])

            #For each point
            for point in range(len(objects[i].getXCoordinates())):
                [align, alignPoint, snapped_x, snapped_y, subsegmentDistance, S, Y] = PvaTools.Geo.getSYfromXY(objects[i].getXCoordinates()[point],objects[i].getYCoordinates()[point],alignments)

                # Error handling
                if(align is False):
                    if(verbose >= 2): print('    Warning: trajectory '+str(i)+' at point '+str(point)+' has alignment errors (spline snapping)')
                    dropped_traj.append(objects[i])
                    objects[i] = None
                    break
                else: objects[i].curvilinearPositions.addPosition(S, Y, align)

            if(objects[i] == None): continue

            ## Go back through points and correct lane  
            #Run through objects looking for outlier point
            smoothed_lanes = PvaTools.Math.cat_mvgavg(objects[i].curvilinearPositions.getLanes(),window=ln_mv_av_win)
            ## Recalculate smoothed lanes
            if(objects[i].curvilinearPositions.getLanes() != smoothed_lanes):
                for point in range(len(objects[i].getXCoordinates())):
                    if(objects[i].curvilinearPositions.getLanes()[point] != smoothed_lanes[point]):
                        [align, alignPoint, snapped_x, snapped_y, subsegmentDistance, S, Y] = PvaTools.Geo.getSYfromXY(objects[i].getXCoordinates()[point],objects[i].getYCoordinates()[point],[alignments[smoothed_lanes[point]]])

                        # Error handling
                        if(align is False):
                            ## This can be triggered by tracking errors when the trajectory jumps around passed another alignment.
                            if(verbose >= 4): print('    Warning: trajectory '+str(i)+' at point '+str(point)+' has alignment errors during trajectory smoothing and will not be corrected.')
                        else:
                            objects[i].curvilinearPositions.setPosition(point, S, Y, align)

        #Resize objects
        if(len(dropped_traj) > 0):
            objects = filter(None, objects)    
            if(verbose >= 2): print('    Filtering report: Trajectories dropped: '+str(len(dropped_traj))) 
        if(verbose >= 2): print('    Filtering report: Lane observation corrections per object: '+str(lane_readjustments/original_object_length))   

        if(reset_objects and len(objects) > 0): return objects[0], dropped_traj
        else:                                   return objects, dropped_traj


    def computeSmoothTrajectory(self, minCommonIntervalLength):
        '''Computes the trajectory as the mean of all features
        if a feature exists, its position is 
        
        Warning work in progress
        TODO? not use the first/last 1-.. positions'''
        from numpy import array, median
        nFeatures = len(self.features)
        if nFeatures == 0:
            print('Empty object features\nCannot compute smooth trajectory')
        else:
            # compute the relative position vectors
            relativePositions = {} # relativePositions[(i,j)] is the position of j relative to i
            for i in xrange(nFeatures):
                for j in xrange(i):
                    fi = self.features[i]
                    fj = self.features[j]
                    inter = fi.commonTimeInterval(fj)
                    if inter.length() >= minCommonIntervalLength:
                        xi = array(fi.getXCoordinates()[inter.first-fi.getFirstInstant():int(fi.length())-(fi.getLastInstant()-inter.last)])
                        yi = array(fi.getYCoordinates()[inter.first-fi.getFirstInstant():int(fi.length())-(fi.getLastInstant()-inter.last)])
                        xj = array(fj.getXCoordinates()[inter.first-fj.getFirstInstant():int(fj.length())-(fj.getLastInstant()-inter.last)])
                        yj = array(fj.getYCoordinates()[inter.first-fj.getFirstInstant():int(fj.length())-(fj.getLastInstant()-inter.last)])
                        relativePositions[(i,j)] = Point(median(xj-xi), median(yj-yi))
                        relativePositions[(j,i)] = -relativePositions[(i,j)]

    ###
    # User Type Classification
    ###
    def classifyUserTypeSpeedMotorized(self, threshold, aggregationFunc = median, ignoreNInstantsAtEnds = 0):
        '''Classifies slow and fast road users
        slow: non-motorized -> pedestrians
        fast: motorized -> cars'''
        if ignoreNInstantsAtEnds > 0:
            speeds = self.getSpeeds()[ignoreNInstantsAtEnds:-ignoreNInstantsAtEnds]
        else:
            speeds = self.getSpeeds()
        if aggregationFunc(speeds) >= threshold:
            self.setUserType(userType2Num['car'])
        else:
            self.setUserType(userType2Num['pedestrian'])

    def classifyUserTypeSpeed(self, speedProbabilities, aggregationFunc = median):
        '''Classifies road user per road user type
        speedProbabilities are functions return P(speed|class)
        in a dictionary indexed by user type names
        Returns probabilities for each class

        for simple threshold classification, simply pass non-overlapping indicator functions (membership)
        e.g. def indic(x):
        if abs(x-mu) < sigma:
        return 1
        else:
        return x'''
        if not hasattr(self, aggregatedSpeed):
            self.aggregatedSpeed = aggregationFunc(self.getSpeeds())
        userTypeProbabilities = {userType2Num[userTypename]: speedProbabilities[userTypename](self.aggregatedSpeed) for userTypename in speedProbabilities}
        self.setUserType(utils.argmaxDict(userTypeProbabilities))
        return userTypeProbabilities

    def initClassifyUserTypeHoGSVM(self, aggregationFunc = median):
        '''Initializes the data structures for classification

        TODO? compute speed for longest feature?
        Skip beginning and end of feature for speed? Offer options instead of median'''
        self.aggregatedSpeed = aggregationFunc(self.getSpeeds())
        self.userTypes = {}

    def classifyUserTypeHoGSVMAtInstant(self, img, pedBikeCarSVM, instant, homography, width, height, bikeCarSVM = None, pedBikeSpeedTreshold = float('Inf'), bikeCarSpeedThreshold = float('Inf'), px = 0.2, py = 0.2, pixelThreshold = 800):
        '''Extract the image box around the object and 
        applies the SVM model on it'''
        from numpy import array
        croppedImg, yCropMin, yCropMax, xCropMin, xCropMax = imageBox(img, self, instant, homography, width, height, px, py, pixelThreshold)
        if len(croppedImg) > 0: # != []
            hog = array([cvutils.HOG(croppedImg)], dtype = np.float32)
            if self.aggregatedSpeed < pedBikeSpeedTreshold or bikeCarSVM == None:
                self.userTypes[instant] = int(pedBikeCarSVM.predict(hog))
            elif self.aggregatedSpeed < bikeCarSpeedTreshold:
                self.userTypes[instant] = int(bikeCarSVM.predict(hog))
            else:
                self.userTypes[instant] = userType2Num['car']
        else:
            self.userTypes[instant] = userType2Num['unknown']

    def classifyUserTypeHoGSVM(self, images, pedBikeCarSVM, homography, width, height, bikeCarSVM = None, pedBikeSpeedTreshold = float('Inf'), bikeCarSpeedThreshold = float('Inf'), speedProbabilities = None, aggregationFunc = median, px = 0.2, py = 0.2, pixelThreshold = 800):
        '''Agregates SVM detections in each image and returns probability
        (proportion of instants with classification in each category)

        iamges is a dictionary of images indexed by instant
        With default parameters, the general (ped-bike-car) classifier will be used
        TODO? consider all categories?'''
        if not hasattr(self, aggregatedSpeed) or not hasattr(self, userTypes):
            print('Initilize the data structures for classification by HoG-SVM')
            self.initClassifyUserTypeHoGSVM(aggregationFunc)

        if len(self.userTypes) != self.length(): # if classification has not been done previously
            for t in self.getTimeInterval():
                if t not in self.userTypes:
                    self.classifyUserTypeHoGSVMAtInstant(images[t], pedBikeCarSVM, t, homography, width, height, bikeCarSVM, pedBikeSpeedTreshold, bikeCarSpeedThreshold, px, py, pixelThreshold)
        # compute P(Speed|Class)
        if speedProbabilities == None: # equiprobable information from speed
            userTypeProbabilities = {userType2Num['car']: 1., userType2Num['pedestrian']: 1., userType2Num['bicycle']: 1.}
        else:
            userTypeProbabilities = {userType2Num[userTypename]: speedProbabilities[userTypename](self.aggregatedSpeed) for userTypename in speedProbabilities}
        # result is P(Class|Appearance) x P(Speed|Class)
        nInstantsUserType = {userType2Num[userTypename]: 0 for userTypename in userTypeProbabilities}# number of instants the object is classified as userTypename
        for t in self.userTypes:
            nInstantsUserType[self.userTypes[t]] += 1
        for userTypename in userTypeProbabilities:
            userTypeProbabilities[userTypename] *= nInstantsUserType[userTypename]
        # class is the user type that maximizes usertype probabilities
        self.setUserType(utils.argmaxDict(userTypeProbabilities))

    def classifyUserTypeArea(self, areas, homography):
        '''Classifies the object based on its location (projected to image space)
        areas is a dictionary of matrix of the size of the image space 
        for different road users possible locations, indexed by road user type names

        TODO: areas could be a wrapper object with a contains method that would work for polygons and images (with wrapper class)
        skip frames at beginning/end?'''
        print('not implemented/tested yet')
        if not hasattr(self, projectedPositions):
            if homography != None:
                self.projectedPositions = obj.positions.project(homography)
            else:
                self.projectedPositions = obj.positions
        possibleUserTypes = {userType: 0 for userType in range(len(userTypenames))}
        for p in self.projectedPositions:
            for userTypename in areas:
                if areas[userTypename][p.x, p.y] != 0:
                    possibleUserTypes[userType2Enum[userTypename]] += 1
        # what to do: threshold for most common type? self.setUserType()
        return possibleUserTypes

    @staticmethod
    def collisionCourseDotProduct(movingObject1, movingObject2, instant):
        'A positive result indicates that the road users are getting closer'
        deltap = movingObject1.getPositionAtInstant(instant)-movingObject2.getPositionAtInstant(instant)
        deltav = movingObject2.getVelocityAtInstant(instant)-movingObject1.getVelocityAtInstant(instant)
        return Point.dot(deltap, deltav)

    @staticmethod
    def collisionCourseCosine(movingObject1, movingObject2, instant):
        'A positive result indicates that the road users are getting closer'
        return Point.cosine(movingObject1.getPositionAtInstant(instant)-movingObject2.getPositionAtInstant(instant), #deltap
                            movingObject2.getVelocityAtInstant(instant)-movingObject1.getVelocityAtInstant(instant)) #deltav

def plotRoadUsers(objects, colors):
    '''Colors is a PlottingPropertyValues instance'''
    from matplotlib.pyplot import figure, axis
    figure()
    for obj in objects:
        obj.plot(colors.get(obj.userType))
    axis('equal')


if __name__ == "__main__":
    import doctest
    import unittest
    suite = doctest.DocFileSuite('tests/moving.txt')
    #suite = doctest.DocTestSuite()
    unittest.TextTestRunner().run(suite)
    #doctest.testmod()
    #doctest.testfile("example.txt")