--- a/trafficintelligence/moving.py	Thu Mar 14 17:08:05 2019 -0400
+++ b/trafficintelligence/moving.py	Fri Mar 15 15:58:37 2019 -0400
@@ -229,6 +229,14 @@
     def __neg__(self):
         return Point(-self.x, -self.y)
 
+    def __mul__(self, alpha):
+        'Warning, returns a new Point'
+        return Point(self.x*alpha, self.y*alpha)
+
+    def divide(self, alpha):
+        'Warning, returns a new Point'
+        return Point(self.x/alpha, self.y/alpha)
+
     def __getitem__(self, i):
         if i == 0:
             return self.x
@@ -244,23 +252,18 @@
         else:
             return Point(-self.y, self.x)
 
+    def normalize(self):
+        return self.divide(self.norm2())
+
     def projectLocal(self, v, clockwise = True):
         'Projects point projected on v, v.orthogonal()'
-        e1 = v/v.norm2()
+        e1 = v.normalize()
         e2 = e1.orthogonal(clockwise)
         return Point(Point.dot(self, e1), Point.dot(self, e2))
 
     def rotate(self, theta):
         return Point(self.x*cos(theta)-self.y*sin(theta), self.x*sin(theta)+self.y*cos(theta))
         
-    def __mul__(self, alpha):
-        'Warning, returns a new Point'
-        return Point(self.x*alpha, self.y*alpha)
-
-    def divide(self, alpha):
-        'Warning, returns a new Point'
-        return Point(self.x/alpha, self.y/alpha)
-
     def plot(self, options = 'o', **kwargs):
         plot([self.x], [self.y], options, **kwargs)
 
@@ -353,9 +356,8 @@
         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 plotAll(points, options = '', **kwargs):
+        plot([p.x for p in points],[p.y for p in points], options, **kwargs)
 
     def similarOrientation(self, refDirection, cosineThreshold):
         'Indicates whether the cosine of the vector and refDirection is smaller than cosineThreshold'
@@ -398,6 +400,27 @@
     def agg(points, aggFunc = mean):
         return Point(aggFunc([p.x for p in points]), aggFunc([p.y for p in points]))
 
+    @staticmethod
+    def boundingRectangle(points, v):
+        '''Returns the bounding rectangle of the points, aligned on the vector v
+        A list of points is returned: front left, front right, rear right, rear left'''
+        e1 = v.normalize()
+        e2 = e1.orthogonal()
+        xCoords = []
+        yCoords = []
+        for p in points:
+            xCoords.append(Point.dot(e1, p))
+            yCoords.append(Point.dot(e2, p))
+        xmin = min(xCoords)
+        xmax = max(xCoords)
+        ymin = min(yCoords)
+        ymax = max(yCoords)
+        frontLeft = Point(xmax, ymin)
+        frontRight = Point(xmax, ymax)
+        rearLeft = Point(xmin, ymin)
+        rearRight = Point(xmin, ymax)
+        return [Point(Point.dot(e1, p), Point.dot(e2, p)) for p in [frontLeft, frontRight, rearRight, rearLeft]]
+    
 if shapelyAvailable:
     def pointsInPolygon(points, polygon):
         '''Optimized tests of a series of points within (Shapely) polygon (not prepared)'''
@@ -509,8 +532,8 @@
         d = s - alignment.getCumulativeDistance(i-1) # distance on subsegment
         #Get difference vector and then snap
         dv = alignment[i] - alignment[i-1]
-        magnitude  = dv.norm2()
-        normalizedV = dv.divide(magnitude)
+        #magnitude  = dv.norm2()
+        normalizedV = dv.normalize()
         #snapped = alignment[i-1] + normalizedV*d # snapped point coordinate along alignment
         # add offset finally
         orthoNormalizedV = normalizedV.orthogonal()
@@ -1761,6 +1784,18 @@
                         relativePositions[(i,j)] = Point(median(xj-xi), median(yj-yi))
                         relativePositions[(j,i)] = -relativePositions[(i,j)]
 
+    def computeBoundingPolygon(self, instant):
+        '''Returns a bounding box for the feature positions at instant
+        bounding box format is a list of points (4 in this case for a rectangle)
+
+        TODO add method argument if using different methods/shapes'''
+        if self.hasFeatures():
+            positions = [f.getPositionAtInstant(instant) for f in self.getFeatures() if f.existsAtInstant(instant)]
+            return Point.boundingRectangle(positions, self.getVelocityAtInstant(instant))
+        else:
+            print('Object {} has no features'.format(self.getNum()))
+            return None
+                        
     ###
     # User Type Classification
     ###

--- a/trafficintelligence/tests/moving.txt	Thu Mar 14 17:08:05 2019 -0400
+++ b/trafficintelligence/tests/moving.txt	Fri Mar 15 15:58:37 2019 -0400
@@ -75,6 +75,12 @@
 >>> Point.distanceNorm2(Point(3,4),Point(1,7))
 3.605551275463989
 
+>>> Point.boundingRectangle([Point(0,0), Point(1,0), Point(0,1), Point(1,1)], Point(1, 1))
+[(0.500000,1.500000), (1.500000,0.500000), (0.500000,-0.500000), (-0.500000,0.500000)]
+>>> Point.boundingRectangle([Point(0,0), Point(1,0), Point(0,1), Point(1,1)], Point(-1, -1))
+[(0.500000,-0.500000), (-0.500000,0.500000), (0.500000,1.500000), (1.500000,0.500000)]
+
+
 >>> Point(3,2).inPolygon(np.array([[0,0],[1,0],[1,1],[0,1]]))
 False
 >>> Point(3,2).inPolygon(np.array([[0,0],[4,0],[4,3],[0,3]]))