changeset 616:0791b3b55b8f

Merge
author MohamedGomaa
date Wed, 10 Dec 2014 14:18:30 -0500
parents 0954aaf28231 (diff) 34111db7cecb (current diff)
children e7f6ca76b7db
files python/moving.py
diffstat 9 files changed, 362 insertions(+), 160 deletions(-) [+]
line wrap: on
line diff
--- a/python/events.py	Wed Dec 03 18:56:49 2014 +0000
+++ b/python/events.py	Wed Dec 10 14:18:30 2014 -0500
@@ -8,12 +8,39 @@
 import multiprocessing
 import itertools
 
-import sys
 import moving, prediction, indicators, utils
-sys.path.append("D:/behaviourAnalysis/libs")
-import trajLearning
 __metaclass__ = type
 
+def findRoute(prototypes,objects,i,j,noiseEntryNums,noiseExitNums,minSimilarity= 0.3, spatialThreshold=1.0, delta=180):
+	if i[0] not in noiseEntryNums: 
+		prototypesRoutes= [ x for x in sorted(prototypes.keys()) if i[0]==x[0]]
+	elif i[1] not in noiseExitNums:
+		prototypesRoutes=[ x for x in sorted(prototypes.keys()) if i[1]==x[1]]
+	else:
+		prototypesRoutes=[x for x in sorted(prototypes.keys())]
+	routeSim={}
+	lcss = utils.LCSS(similarityFunc=lambda x,y: (distanceForLCSS(x,y) <= spatialThreshold),delta=delta)
+	for y in prototypesRoutes: 
+		if y in prototypes.keys():
+			prototypesIDs=prototypes[y]
+			similarity=[]
+			for x in prototypesIDs:
+				s=lcss.computeNormalized(objects[j].positions, objects[x].positions)
+				similarity.append(s)
+			routeSim[y]=max(similarity)
+	route=max(routeSim, key=routeSim.get)
+	if routeSim[route]>=minSimilarity:
+		return route
+	else:
+		return i
+
+def getRoute(obj,prototypes,objects,noiseEntryNums,noiseExitNums,useDestination=True):
+    route=(obj.startRouteID,obj.endRouteID)
+    if useDestination:
+        if route not in prototypes.keys():
+            route= findRoute(prototypes,objects,route,obj.getNum(),noiseEntryNums,noiseExitNums)
+    return route
+
 class Interaction(moving.STObject):
     '''Class for an interaction between two road users 
     or a road user and an obstacle
@@ -126,17 +153,20 @@
                 minDistance[instant] = moving.MovingObject.minDistance(self.roadUser1, self.roadUser2, instant)
             self.addIndicator(indicators.SeverityIndicator(Interaction.indicatorNames[3], minDistance))
 
-    def computeCrossingsCollisions(self, predictionParameters, collisionDistanceThreshold, timeHorizon, computeCZ = False, debug = False, timeInterval = None, nProcesses = 1):
+    def computeCrossingsCollisions(self, predictionParameters, collisionDistanceThreshold, timeHorizon, computeCZ = False, debug = False, timeInterval = None, nProcesses = 1,usePrototypes=True,route1= (-1,-1),route2=(-1,-1),prototypes={},secondStepPrototypes={},nMatching={},objects=[],noiseEntryNums=[],noiseExitNums=[],minSimilarity=0.1,mostMatched=None,useDestination=True,useSpeedPrototype=True,acceptPartialLength=30, step=1):
         '''Computes all crossing and collision points at each common instant for two road users. '''
         self.collisionPoints={}
         self.crossingZones={}
         TTCs = {}
+        if usePrototypes:
+            route1= getRoute(self.roadUser1,prototypes,objects,noiseEntryNums,noiseExitNums,useDestination)
+            route2= getRoute(self.roadUser2,prototypes,objects,noiseEntryNums,noiseExitNums,useDestination)
 
         if timeInterval:
             commonTimeInterval = timeInterval
         else:
             commonTimeInterval = self.timeInterval
-        self.collisionPoints, self.crossingZones = prediction.computeCrossingsCollisions(predictionParameters, self.roadUser1, self.roadUser2, collisionDistanceThreshold, timeHorizon, computeCZ, debug, commonTimeInterval, nProcesses)
+        self.collisionPoints, self.crossingZones = prediction.computeCrossingsCollisions(predictionParameters, self.roadUser1, self.roadUser2, collisionDistanceThreshold, timeHorizon, computeCZ, debug, commonTimeInterval, nProcesses,usePrototypes,route1,route2,prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination,useSpeedPrototype,acceptPartialLength, step)
         for i, cp in self.collisionPoints.iteritems():
             TTCs[i] = prediction.SafetyPoint.computeExpectedIndicator(cp)
         # add probability of collision, and probability of successful evasive action
@@ -144,70 +174,8 @@
         
         if computeCZ:
             pPETs = {}
-            for i in list(commonTimeInterval)[:-1]:
-                if len(self.crossingZones[i]) > 0:
-                    pPETs[i] = prediction.SafetyPoint.computeExpectedIndicator(self.crossingZones[i])
-            self.addIndicator(indicators.SeverityIndicator(Interaction.indicatorNames[9], pPETs))
-
-    def computeCrosssingCollisionsPrototypeAtInstant(self, instant,route1,route2,predictionParameters, collisionDistanceThreshold, timeHorizon, prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity=0.1,mostMatched=None, computeCZ = False, debug = False,useDestination=True,useSpeedPrototype=True):
-        inter1=moving.Interval(self.roadUser1.timeInterval.first,instant)
-        inter2=moving.Interval(self.roadUser2.timeInterval.first,instant)
-        partialObjPositions1= self.roadUser1.getObjectInTimeInterval(inter1).positions
-        partialObjPositions2= self.roadUser2.getObjectInTimeInterval(inter2).positions
-        if useSpeedPrototype:
-            prototypeTrajectories1=trajLearning.findPrototypesSpeed(prototypes,secondStepPrototypes,nMatching,objects,route1,partialObjPositions1,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination)
-            prototypeTrajectories2=trajLearning.findPrototypesSpeed(prototypes,secondStepPrototypes,nMatching,objects,route2,partialObjPositions2,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination)
-        else:
-            prototypeTrajectories1=trajLearning.findPrototypes(prototypes,nMatching,objects,route1,partialObjPositions1,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched)
-            prototypeTrajectories2=trajLearning.findPrototypes(prototypes,nMatching,objects,route2,partialObjPositions2,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched)
-        if prototypeTrajectories1=={}:
-            print self.roadUser1.num, 'is abnormal at instant', str(instant)
-            return [],[]
-        elif prototypeTrajectories2=={}:
-            print self.roadUser2.num, 'is abnormal at instant', str(instant)
-            return [],[]
-        else:
-            currentInstant,collisionPoints, crossingZones = predictionParameters.computeCrossingsCollisionsAtInstant(instant, self.roadUser1, self.roadUser2, collisionDistanceThreshold, timeHorizon, computeCZ, debug,prototypeTrajectories1,prototypeTrajectories2)
-            return collisionPoints,crossingZones
-
-    def computeCrossingsCollisionsPrototype2stages(self, predictionParameters, collisionDistanceThreshold, timeHorizon, prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,step=10,minSimilarity=0.1,mostMatched=None, computeCZ = False, debug = False, timeInterval = None,acceptPartialLength=30,useDestination=True,useSpeedPrototype=True):
-        '''Computes all crossing and collision points at each common instant for two road users. '''
-        self.collisionPoints={}
-        self.crossingZones={}
-        TTCs = {}
-        route1=(self.roadUser1.startRouteID,self.roadUser1.endRouteID)
-        route2=(self.roadUser2.startRouteID,self.roadUser2.endRouteID)
-        if useDestination:
-            if route1 not in prototypes.keys():
-                route1= trajLearning.findRoute(prototypes,objects,route1,self.roadUser1.num,noiseEntryNums,noiseExitNums)
-            if route2 not in prototypes.keys():
-                route2= trajLearning.findRoute(prototypes,objects,route2,self.roadUser2.num,noiseEntryNums,noiseExitNums)
-				
-        if timeInterval:
-            commonTimeInterval = timeInterval
-        else:
-            commonTimeInterval = self.timeInterval
-        reCompute=False
-        for i in xrange(commonTimeInterval.first,commonTimeInterval.last,step): # incremental calculation of CP,CZ to save time
-            if  i-self.roadUser1.timeInterval.first >= acceptPartialLength and i-self.roadUser2.timeInterval.first >= acceptPartialLength:
-				self.collisionPoints[i], self.crossingZones[i] = self.computeCrosssingCollisionsPrototypeAtInstant(i,route1,route2,predictionParameters, collisionDistanceThreshold, timeHorizon, prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched, computeCZ, debug,useDestination,useSpeedPrototype)
-				if len(self.collisionPoints[i]) >0 or len(self.crossingZones[i])>0:
-					reCompute=True
-					break
-        if reCompute:
-            for i in list(commonTimeInterval)[:-1]: # do not look at the 1 last position/velocities, often with errors
-				if  i-self.roadUser1.timeInterval.first >= acceptPartialLength and i-self.roadUser2.timeInterval.first >= acceptPartialLength:
-					self.collisionPoints[i], self.crossingZones[i] = self.computeCrosssingCollisionsPrototypeAtInstant(i,route1,route2,predictionParameters, collisionDistanceThreshold, timeHorizon, prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched, computeCZ, debug,useDestination,useSpeedPrototype)		
-					if len(self.collisionPoints[i]) > 0:
-						TTCs[i] = prediction.SafetyPoint.computeExpectedIndicator(self.collisionPoints[i])
-        # add probability of collision, and probability of successful evasive action
-        self.addIndicator(indicators.SeverityIndicator(Interaction.indicatorNames[7], TTCs))
-        
-        if computeCZ:
-            pPETs = {}
-            for i in list(commonTimeInterval)[:-1]:
-                if i in self.crossingZones.keys() and len(self.crossingZones[i]) > 0:
-                    pPETs[i] = prediction.SafetyPoint.computeExpectedIndicator(self.crossingZones[i])
+            for i, cz in self.crossingZones.iteritems():
+                pPETs[i] = prediction.SafetyPoint.computeExpectedIndicator(cz)
             self.addIndicator(indicators.SeverityIndicator(Interaction.indicatorNames[9], pPETs))
 			
     def addVideoFilename(self,videoFilename):
--- a/python/moving.py	Wed Dec 03 18:56:49 2014 +0000
+++ b/python/moving.py	Wed Dec 10 14:18:30 2014 -0500
@@ -924,7 +924,8 @@
 
 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 
+    with a trajectory and a geometry (constant volume over time) 
+    and a usertype (e.g. road user) coded as a number (see userTypeNames)
     '''
 
     def __init__(self, num = None, timeInterval = None, positions = None, velocities = None, geometry = None, userType = userType2Num['unknown']):
@@ -1279,6 +1280,24 @@
         return Point.cosine(movingObject1.getPositionAtInstant(instant)-movingObject2.getPositionAtInstant(instant), #deltap
                             movingObject2.getVelocityAtInstant(instant)-movingObject1.getVelocityAtInstant(instant)) #deltav
 
+
+##################
+# Annotations
+##################
+
+class BBAnnotation(MovingObject):
+    '''Class for : a series of ground truth annotations using bounding boxes
+    Its center is the center of the containing shape
+    '''
+
+    def __init__(self, num = None, timeInterval = None, topPositions = None, bottomPositions = None, userType = userType2Num['unknown']):
+        super(BBAnnotation, self).__init__(num, timeInterval, Trajectory(), userType = userType)
+        self.topPositions = topPositions.getPositions()
+        self.bottomPositions = bottomPositions.getPositions()
+        for i in xrange(int(topPositions.length())):
+            self.positions.addPosition((topPositions.getPositionAt(i) + bottomPositions.getPositionAt(i)).multiply(0.5))
+        
+
 def plotRoadUsers(objects, colors):
     '''Colors is a PlottingPropertyValues instance'''
     from matplotlib.pyplot import figure, axis
--- a/python/objectsmoothing.py	Wed Dec 03 18:56:49 2014 +0000
+++ b/python/objectsmoothing.py	Wed Dec 10 14:18:30 2014 -0500
@@ -13,13 +13,6 @@
 			dist[f]= moving.Point.distanceNorm2(feat.getPositionAtInstant(t-1),f.getPositionAtInstant(t))
 	return min(dist, key=dist.get) # = utils.argmaxDict(dist)
 	
-def FeatureList(obj,minLengthParam=0.7):
-	featureList=[]
-	for i in obj.features:
-		if i.length>= minLengthParam*obj.length():
-			featureList.append(i.num)
-	return featureList
-	
 def getFeatures(obj,features,featureID):
 	#longestFeature = utils.argmaxDict({f:f.length() for i,f in enumerate(obj.features)})
 	t1,t3 = features[featureID].getFirstInstant(), features[featureID].getLastInstant()
@@ -80,7 +73,7 @@
 		p1= object.getPositionAtInstant(i)
 		p2= object.getPositionAtInstant(i+1)
 		velocities[i]=p2-p1		
-	velocities[object.timeInterval.last]= velocities[object.timeInterval.last-1]  # duplicate last point
+	velocities[object.getLastInstant()]= velocities[object.getLastInstant()-1]  # duplicate last point
 	if smoothing:
 		velX= [velocities[y].aslist()[0] for y in sorted(velocities.keys())]
 		velY= [velocities[y].aslist()[1] for y in sorted(velocities.keys())]
@@ -121,7 +114,7 @@
 		t+= jerk[i]* jerk[i]
 	return t
 	
-def getObject(obj,features,featureID,newNum,smoothing=False,halfWidth=3,create=False):
+def smoothObjectTrajectory(obj,features,featureID,newNum,smoothing=False,halfWidth=3,create=False):
 	results=[]	
 	bearing={}
 	if create:
@@ -131,7 +124,7 @@
 	for t in feature.getTimeInterval():
 		p1= feature.getPositionAtInstant(t)
 		p2= obj.getPositionAtInstant(t)
-		if t!=feature.timeInterval.last:
+		if t!=feature.getLastInstant():
 			p3= feature.getPositionAtInstant(t+1)
 		else:
 			p1= feature.getPositionAtInstant(t-1)
@@ -174,24 +167,22 @@
 		d1= translated[halfWidth]- feature.positions[halfWidth]
 		d2= translated[-halfWidth-1]- feature.positions[-halfWidth-1]
 		for i in xrange(halfWidth):
-			p1.x=feature.positions.__getitem__(i).x+d1.x
-			p2.x= feature.positions.__getitem__(-i-1).x+d2.x
-			p1.y=feature.positions.__getitem__(i).y+d1.y
-			p2.y= feature.positions.__getitem__(-i-1).y+d2.y
+			p1= feature.getPositionAt(i)+d1
+			p2= feature.getPositionAt(-i-1)+d2
 			translated.setPosition(i,p1)
 			translated.setPosition(-i-1,p2)
 		
-	newObj= moving.MovingObject(newNum,timeInterval=feature.timeInterval,positions=translated)
+	newObj= moving.MovingObject(newNum,timeInterval=feature.getTimeInterval(),positions=translated)
 	return newObj
 	
 def smoothObjectTrajectory(obj,features,newNum,minLengthParam=0.7,smoothing=False,plotResults=True,halfWidth=3,computeVelocities=True,optimize=True,create=False):
-	featureList=FeatureList(obj,minLengthParam=minLengthParam)
+	featureList=[i.num for i in obj.features if i.length() >= minLengthParam*obj.length()]
 	if featureList==[]:
 		featureList.append(longestFeature(obj))
 		create=True
 	objs=[]
 	for featureID in featureList:
-		objTMP=getObject(obj,features,featureID,newNum,smoothing=smoothing,halfWidth=halfWidth,create=create)
+		objTMP=smoothObjectTrajectory(obj,features,featureID,newNum,smoothing=smoothing,halfWidth=halfWidth,create=create)
 		objs.append(objTMP)
 	newTranslated = moving.Trajectory()
 	newInterval=[]
--- a/python/poly-utils.py	Wed Dec 03 18:56:49 2014 +0000
+++ b/python/poly-utils.py	Wed Dec 10 14:18:30 2014 -0500
@@ -6,26 +6,18 @@
 import numpy as np
 
 __metaclass__ = type
+from indicators import SeverityIndicator
 
-# inputs variables	
-#dirname= 'G:/0-phdstart/Code/trial/indicatorsNew/'
-#extension= '-indicatorsNew.csv'
-#indicatorsNames= {1:'Distance',2:'Cosine',3:'collision Course Angle',4:'Velocity Cosine',5:'Velocity Angle',6:'Speed Differential',7:'Collision Probability',8:'Severity Index',9:'TTC'}
-''' min Distance case'''
-dirname= 'G:/0-phdstart/Code/trial/minDistanceIndicator/'
-extension= '-minDistanceInd.csv'
-indicatorsNames= {1:'minDistance'}
 
-def loadNewInteractions(videoFilename,interactionType, roaduserNum1,roaduserNum2, selectedIndicators=[]):
+def loadNewInteractions(videoFilename,interactionType,dirname, extension, indicatorsNames, roaduserNum1,roaduserNum2, selectedIndicators=[]):
     '''Loads interactions from the POLY traffic event format'''
     from events import Interaction 
-    from indicators import SeverityIndicator
-    #filename= dirname + videoFilename + extension
-    filename= dirname + interactionType+ '-' + videoFilename + extension # case of min distance todo: change the saving format to be matched with all outputs
+    filename= dirname + videoFilename + extension
+    #filename= dirname + interactionType+ '-' + videoFilename + extension # case of min distance todo: change the saving format to be matched with all outputs
     file = utils.openCheck(filename)
     if (not file):
         return []
-    interactions = []
+    #interactions = []
     interactionNum = 0
     data= np.loadtxt(filename)
     indicatorFrameNums= data[:,0]
@@ -43,7 +35,92 @@
             values[t] = [data[i,index] for index in selectedIndicators]
         inter.addIndicator(SeverityIndicator('selectedIndicators', values))	
 		
-    interactions.append(inter)
+    #interactions.append(inter)
     file.close()
-    return interactions
+    #return interactions
+    return inter
+
+# Plotting results
+
+frameRate = 15.
+
+# To run in directory that contains the directories that contain the results (Miss-xx and Incident-xx)
+#dirname = '/home/nicolas/Research/Data/kentucky-db/'
+
+interactingRoadUsers = {'Miss/0404052336': [(0,3)] # 0,2 and 1 vs 3
+                        #,
+                        #'Incident/0306022035': [(1,3)]
+                        #,
+                        #'Miss/0208030956': [(4,5),(5,7)]
+                        }
+
+
+def getIndicatorName(filename, withUnit = False):
+    if withUnit:
+        unit = ' (s)'
+    else:
+        unit = ''
+    if 'collision-point' in filename:
+        return 'TTC'+unit
+    elif 'crossing' in filename:
+        return 'pPET'+unit
+    elif 'probability' in filename:
+        return 'P(UEA)'
+
+def getMethodName(fileprefix):
+    if fileprefix == 'constant-velocity':
+        return 'Con. Vel.'
+    elif fileprefix == 'normal-adaptation':
+        return 'Norm. Ad.'
+    elif fileprefix == 'point-set':
+        return 'Pos. Set'
+    elif fileprefix == 'evasive-action':
+        return 'Ev. Act.'
+    elif fileprefix == 'point-set-evasive-action':
+        return 'Pos. Set'
+
+indicator2TimeIdx = {'TTC':2,'pPET':2, 'P(UEA)':3}
 
+def getDataAtInstant(data, i):
+    return data[data[:,2] == i]
+
+def getPointsAtInstant(data, i):
+    return getDataAtInstant(i)[3:5]
+
+def getIndicator(data, roadUserNumbers, indicatorName):
+    if data.ndim ==1:
+        data.shape = (1,data.shape[0])
+
+    # find the order for the roadUserNumbers
+    uniqueObj1 = np.unique(data[:,0])
+    uniqueObj2 = np.unique(data[:,1])
+    found = False
+    if roadUserNumbers[0] in uniqueObj1 and roadUserNumbers[1] in uniqueObj2:
+        objNum1 = roadUserNumbers[0]
+        objNum2 = roadUserNumbers[1]
+        found = True
+    if roadUserNumbers[1] in uniqueObj1 and roadUserNumbers[0] in uniqueObj2:
+        objNum1 = roadUserNumbers[1]
+        objNum2 = roadUserNumbers[0]
+        found = True
+
+    # get subset of data for road user numbers
+    if found:
+        roadUserData = data[np.logical_and(data[:,0] == objNum1, data[:,1] == objNum2),:]
+        if roadUserData.size > 0:
+            time = np.unique(roadUserData[:,indicator2TimeIdx[indicatorName]])
+            values = {}
+            if indicatorName == 'P(UEA)':
+                tmp = roadUserData[:,4]
+                for k,v in zip(time, tmp):
+                    values[k]=v
+                return SeverityIndicator(indicatorName, values, mostSevereIsMax = False, maxValue = 1.), roadUserData
+            else:
+                for i in xrange(time[0],time[-1]+1):
+                    try:
+						tmp = getDataAtInstant(roadUserData, i)
+						values[i] = np.sum(tmp[:,5]*tmp[:,6])/np.sum(tmp[:,5])/frameRate
+                    except IOError:
+                    	values[i] = np.inf
+                return SeverityIndicator(indicatorName, values, mostSevereIsMax = False), roadUserData
+    return None, None
\ No newline at end of file
--- a/python/prediction.py	Wed Dec 03 18:56:49 2014 +0000
+++ b/python/prediction.py	Wed Dec 10 14:18:30 2014 -0500
@@ -5,6 +5,7 @@
 import math
 import random
 import numpy as np
+from utils import LCSS
 
 class PredictedTrajectory:
     '''Class for predicted trajectories with lazy evaluation
@@ -166,14 +167,94 @@
     savefig('predicted-trajectories-t-{0}.png'.format(currentInstant))
     close()
 
-def computeCrossingsCollisionsAtInstant(predictionParams,currentInstant, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ = False, debug = False,prototypeTrajectories1=None,prototypeTrajectories2=None):  
+def calculateProbability(nMatching,similarity,objects):
+	sumFrequencies=sum([nMatching[p] for p in similarity.keys()])
+	prototypeProbability={}
+	for i in similarity.keys():
+		prototypeProbability[i]= similarity[i] * float(nMatching[i])/sumFrequencies
+	sumProbabilities= sum([prototypeProbability[p] for p in prototypeProbability.keys()])
+	probabilities={}
+	for i in prototypeProbability.keys():
+		probabilities[objects[i]]= float(prototypeProbability[i])/sumProbabilities
+	return probabilities
+
+def findPrototypes(prototypes,nMatching,objects,route,partialObjPositions,noiseEntryNums,noiseExitNums,minSimilarity=0.1,mostMatched=None,spatialThreshold=1.0, delta=180):
+	''' behaviour prediction first step'''
+	if route[0] not in noiseEntryNums: 
+		prototypesRoutes= [ x for x in sorted(prototypes.keys()) if route[0]==x[0]]
+	elif route[1] not in noiseExitNums:
+		prototypesRoutes=[ x for x in sorted(prototypes.keys()) if route[1]==x[1]]
+	else:
+		prototypesRoutes=[x for x in sorted(prototypes.keys())]
+	lcss = LCSS(similarityFunc=lambda x,y: (distanceForLCSS(x,y) <= spatialThreshold),delta=delta)
+	similarity={}
+	for y in prototypesRoutes: 
+		if y in prototypes.keys():
+			prototypesIDs=prototypes[y]			
+			for x in prototypesIDs:
+				s=lcss.computeNormalized(partialObjPositions, objects[x].positions)
+				if s >= minSimilarity:
+					similarity[x]=s
+	
+	if mostMatched==None:
+		probabilities= calculateProbability(nMatching,similarity,objects)		
+		return probabilities
+	else:
+		mostMatchedValues=sorted(similarity.values(),reverse=True)[:mostMatched]
+		keys=[k for k in similarity.keys() if similarity[k] in mostMatchedValues]
+		newSimilarity={}
+		for i in keys:
+			newSimilarity[i]=similarity[i]
+		probabilities= calculateProbability(nMatching,newSimilarity,objects)		
+		return probabilities		
+		
+def findPrototypesSpeed(prototypes,secondStepPrototypes,nMatching,objects,route,partialObjPositions,noiseEntryNums,noiseExitNums,minSimilarity=0.1,mostMatched=None,useDestination=True,spatialThreshold=1.0, delta=180):
+	if useDestination:
+		prototypesRoutes=[route]
+	else:
+		if route[0] not in noiseEntryNums: 
+			prototypesRoutes= [ x for x in sorted(prototypes.keys()) if route[0]==x[0]]
+		elif route[1] not in noiseExitNums:
+			prototypesRoutes=[ x for x in sorted(prototypes.keys()) if route[1]==x[1]]
+		else:
+			prototypesRoutes=[x for x in sorted(prototypes.keys())]
+	lcss = LCSS(similarityFunc=lambda x,y: (distanceForLCSS(x,y) <= spatialThreshold),delta=delta)
+	similarity={}
+	for y in prototypesRoutes: 
+		if y in prototypes.keys():
+			prototypesIDs=prototypes[y]	
+			for x in prototypesIDs:
+				s=lcss.computeNormalized(partialObjPositions, objects[x].positions)
+				if s >= minSimilarity:
+					similarity[x]=s
+	
+	newSimilarity={}
+	for i in similarity.keys():
+		if i in secondStepPrototypes.keys():
+			for j in secondStepPrototypes[i]:
+				newSimilarity[j]=similarity[i]
+	probabilities= calculateProbability(nMatching,newSimilarity,objects)		
+	return probabilities
+	
+def getPrototypeTrajectory(obj,route,currentInstant,prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity=0.1,mostMatched=None,useDestination=True,useSpeedPrototype=True):
+    partialInterval=moving.Interval(obj.getFirstInstant(),currentInstant)
+    partialObjPositions= obj.getObjectInTimeInterval(partialInterval).positions	
+    if useSpeedPrototype:
+        prototypeTrajectories=findPrototypesSpeed(prototypes,secondStepPrototypes,nMatching,objects,route,partialObjPositions,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination)
+    else:
+        prototypeTrajectories=findPrototypes(prototypes,nMatching,objects,route,partialObjPositions,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched)
+    return 	prototypeTrajectories
+
+def computeCrossingsCollisionsAtInstant(predictionParams,currentInstant, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ = False, debug = False,usePrototypes=True,route1= (-1,-1),route2=(-1,-1),prototypes={},secondStepPrototypes={},nMatching={},objects=[],noiseEntryNums=[],noiseExitNums=[],minSimilarity=0.1,mostMatched=None,useDestination=True,useSpeedPrototype=True):  
     '''returns the lists of collision points and crossing zones'''
-    if prototypeTrajectories1==None: 
+    if usePrototypes:
+		prototypeTrajectories1=getPrototypeTrajectory(obj1,route1,currentInstant,prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination,useSpeedPrototype)
+		prototypeTrajectories2= getPrototypeTrajectory(obj2,route2,currentInstant,prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination,useSpeedPrototype)
+		predictedTrajectories1 = predictionParams.generatePredictedTrajectories(obj1, currentInstant,prototypeTrajectories1)
+		predictedTrajectories2 = predictionParams.generatePredictedTrajectories(obj2, currentInstant,prototypeTrajectories2) 	
+    else:
         predictedTrajectories1 = predictionParams.generatePredictedTrajectories(obj1, currentInstant)
         predictedTrajectories2 = predictionParams.generatePredictedTrajectories(obj2, currentInstant)		
-    else:
-        predictedTrajectories1 = predictionParams.generatePredictedTrajectories(obj1, currentInstant,prototypeTrajectories1)
-        predictedTrajectories2 = predictionParams.generatePredictedTrajectories(obj2, currentInstant,prototypeTrajectories2)
 
     collisionPoints = []
     crossingZones = []
@@ -205,7 +286,7 @@
         savePredictedTrajectoriesFigure(currentInstant, obj1, obj2, predictedTrajectories1, predictedTrajectories2, timeHorizon)
     return currentInstant,collisionPoints, crossingZones
 
-def computeCrossingsCollisions(predictionParams, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ = False, debug = False, timeInterval = None,nProcesses = 1,prototypeTrajectories1=None,prototypeTrajectories2=None):
+def computeCrossingsCollisions(predictionParams, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ = False, debug = False, timeInterval = None,nProcesses = 1,usePrototypes=True,route1= (-1,-1),route2=(-1,-1),prototypes={},secondStepPrototypes={},nMatching={},objects=[],noiseEntryNums=[],noiseExitNums=[],minSimilarity=0.1,mostMatched=None,useDestination=True,useSpeedPrototype=True,acceptPartialLength=30, step=1):
     '''Computes all crossing and collision points at each common instant for two road users. '''
     collisionPoints={}
     crossingZones={}
@@ -214,16 +295,35 @@
     else:
         commonTimeInterval = obj1.commonTimeInterval(obj2)
     if nProcesses == 1:
-        for i in list(commonTimeInterval)[:-1]: # do not look at the 1 last position/velocities, often with errors
-            i, cp, cz = computeCrossingsCollisionsAtInstant(predictionParams, i, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ, debug,prototypeTrajectories1,prototypeTrajectories2)
-            if len(cp) != 0:
-                collisionPoints[i] = cp
-            if len(cz) != 0:
-                 crossingZones[i] = cz
+        if usePrototypes:
+            firstInstant= next( (x for x in xrange(commonTimeInterval.first,commonTimeInterval.last) if x-obj1.getFirstInstant() >= acceptPartialLength and x-obj2.getFirstInstant() >= acceptPartialLength), commonTimeInterval.last)
+            commonTimeIntervalList1= list(xrange(firstInstant,commonTimeInterval.last-1)) # do not look at the 1 last position/velocities, often with errors
+            commonTimeIntervalList2= list(xrange(firstInstant,commonTimeInterval.last-1,step)) # do not look at the 1 last position/velocities, often with errors
+            for i in commonTimeIntervalList2: 
+                i, cp, cz = computeCrossingsCollisionsAtInstant(predictionParams, i, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ, debug,usePrototypes,route1,route2,prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination,useSpeedPrototype)
+                if len(cp) != 0:
+                    collisionPoints[i] = cp
+                if len(cz) != 0:
+                    crossingZones[i] = cz
+            if collisionPoints!={} or crossingZones!={}:
+                for i in commonTimeIntervalList1:
+                    if i not in commonTimeIntervalList2:
+                        i, cp, cz = computeCrossingsCollisionsAtInstant(predictionParams, i, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ, debug,usePrototypes,route1,route2,prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination,useSpeedPrototype)
+                        if len(cp) != 0:
+                            collisionPoints[i] = cp
+                        if len(cz) != 0:
+                            crossingZones[i] = cz						
+        else:
+            for i in list(commonTimeInterval)[:-1]: # do not look at the 1 last position/velocities, often with errors
+                i, cp, cz = computeCrossingsCollisionsAtInstant(predictionParams, i, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ, debug,usePrototypes,route1,route2,prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination,useSpeedPrototype)
+                if len(cp) != 0:
+                    collisionPoints[i] = cp
+                if len(cz) != 0:
+                    crossingZones[i] = cz
     else:
         from multiprocessing import Pool
         pool = Pool(processes = nProcesses)
-        jobs = [pool.apply_async(computeCrossingsCollisionsAtInstant, args = (predictionParams, i, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ, debug,prototypeTrajectories1,prototypeTrajectories2)) for i in list(commonTimeInterval)[:-1]]
+        jobs = [pool.apply_async(computeCrossingsCollisionsAtInstant, args = (predictionParams, i, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ, debug,usePrototypes,route1,route2,prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination,useSpeedPrototype)) for i in list(commonTimeInterval)[:-1]]
         #results = [j.get() for j in jobs]
         #results.sort()
         for j in jobs:
@@ -247,11 +347,11 @@
     def generatePredictedTrajectories(self, obj, instant):
         return []
 
-    def computeCrossingsCollisionsAtInstant(self, currentInstant, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ = False, debug = False,prototypeTrajectories1=None,prototypeTrajectories2=None):
-        return computeCrossingsCollisionsAtInstant(self, currentInstant, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ, debug,prototypeTrajectories1,prototypeTrajectories2)
+    def computeCrossingsCollisionsAtInstant(self, currentInstant, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ = False, debug = False,usePrototypes=True,route1= (-1,-1),route2=(-1,-1),prototypes={},secondStepPrototypes={},nMatching={},objects=[],noiseEntryNums=[],noiseExitNums=[],minSimilarity=0.1,mostMatched=None,useDestination=True,useSpeedPrototype=True):
+        return computeCrossingsCollisionsAtInstant(self, currentInstant, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ, debug,usePrototypes,route1,route2,prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination,useSpeedPrototype)
 
-    def computeCrossingsCollisions(self, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ = False, debug = False, timeInterval = None, nProcesses = 1,prototypeTrajectories1=None,prototypeTrajectories2=None):
-        return computeCrossingsCollisions(self, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ, debug, timeInterval, nProcesses,prototypeTrajectories1,prototypeTrajectories2)
+    def computeCrossingsCollisions(self, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ = False, debug = False, timeInterval = None, nProcesses = 1,usePrototypes=True,route1= (-1,-1),route2=(-1,-1),prototypes={},secondStepPrototypes={},nMatching={},objects=[],noiseEntryNums=[],noiseExitNums=[],minSimilarity=0.1,mostMatched=None,useDestination=True,useSpeedPrototype=True,acceptPartialLength=30, step=1):
+        return computeCrossingsCollisions(self, obj1, obj2, collisionDistanceThreshold, timeHorizon, computeCZ, debug, timeInterval, nProcesses,usePrototypes,route1,route2,prototypes,secondStepPrototypes,nMatching,objects,noiseEntryNums,noiseExitNums,minSimilarity,mostMatched,useDestination,useSpeedPrototype,acceptPartialLength, step)
 
     def computeCollisionProbability(self, obj1, obj2, collisionDistanceThreshold, timeHorizon, debug = False, timeInterval = None):
         '''Computes only collision probabilities
--- a/python/storage.py	Wed Dec 03 18:56:49 2014 +0000
+++ b/python/storage.py	Wed Dec 10 14:18:30 2014 -0500
@@ -34,6 +34,7 @@
     except sqlite3.OperationalError as error:
         printDBError(error)
 
+# TODO: add test if database connection is open
 # IO to sqlite
 def writeTrajectoriesToSqlite(objects, outputFilename, trajectoryType, objectNumbers = -1):
     """
@@ -293,20 +294,21 @@
     if trajectoryType == 'feature':
         statementBeginning = 'where trajectory_id '
     elif trajectoryType == 'object':
-        statementBeginning =  'and OF.object_id '
+        statementBeginning = 'and OF.object_id '
+    elif trajectoryType == 'bbtop' or 'bbbottom':
+        statementBeginning = 'where object_id '
     else:
         print('no trajectory type was chosen')
 
-    if type(objectNumbers) == int:
-        if objectNumbers == -1:
-            query = ''
-        else:
-            query = statementBeginning+'between 0 and {0} '.format(objectNumbers)
+    if objectNumbers is None:
+        query = ''
+    elif type(objectNumbers) == int:
+        query = statementBeginning+'between 0 and {0} '.format(objectNumbers)
     elif type(objectNumbers) == list:
         query = statementBeginning+'in ('+', '.join([str(n) for n in objectNumbers])+') '
     return query
 
-def loadTrajectoriesFromTable(connection, tableName, trajectoryType, objectNumbers = -1):
+def loadTrajectoriesFromTable(connection, tableName, trajectoryType, objectNumbers = None):
     '''Loads trajectories (in the general sense) from the given table
     can be positions or velocities
 
@@ -314,14 +316,21 @@
     cursor = connection.cursor()
 
     try:
+        idQuery = getTrajectoryIdQuery(objectNumbers, trajectoryType)
         if trajectoryType == 'feature':
-            trajectoryIdQuery = getTrajectoryIdQuery(objectNumbers, trajectoryType)
-            queryStatement = 'SELECT * from '+tableName+' '+trajectoryIdQuery+'order by trajectory_id, frame_number'
+            queryStatement = 'SELECT * from '+tableName+' '+idQuery+'ORDER BY trajectory_id, frame_number'
             cursor.execute(queryStatement)
             logging.debug(queryStatement)
         elif trajectoryType == 'object':
-            objectIdQuery = getTrajectoryIdQuery(objectNumbers, trajectoryType)
-            queryStatement = 'SELECT OF.object_id, P.frame_number, avg(P.x_coordinate), avg(P.y_coordinate) from '+tableName+' P, objects_features OF where P.trajectory_id = OF.trajectory_id '+objectIdQuery+'group by OF.object_id, P.frame_number order by OF.object_id, P.frame_number'
+            queryStatement = 'SELECT OF.object_id, P.frame_number, avg(P.x_coordinate), avg(P.y_coordinate) from '+tableName+' P, objects_features OF where P.trajectory_id = OF.trajectory_id '+idQuery+'group by OF.object_id, P.frame_number ORDER BY OF.object_id, P.frame_number'
+            cursor.execute(queryStatement)
+            logging.debug(queryStatement)
+        elif trajectoryType in ['bbtop', 'bbbottom']:
+            if trajectoryType == 'bbtop':
+                corner = 'top_left'
+            elif trajectoryType == 'bbbottom':
+                corner = 'bottom_right'
+            queryStatement = 'SELECT object_id, frame_number, x_'+corner+', y_'+corner+' FROM '+tableName+' '+trajectoryIdQuery+'ORDER BY object_id, frame_number'
             cursor.execute(queryStatement)
             logging.debug(queryStatement)
         else:
@@ -336,21 +345,36 @@
     for row in cursor:
         if row[0] != objId:
             objId = row[0]
-            if obj:
+            if obj != None and obj.length() == obj.positions.length():
                 objects.append(obj)
+            elif obj != None:
+                print('Object {} is missing {} positions'.format(obj.getNum(), int(obj.length())-obj.positions.length()))
             obj = moving.MovingObject(row[0], timeInterval = moving.TimeInterval(row[1], row[1]), positions = moving.Trajectory([[row[2]],[row[3]]]))
         else:
             obj.timeInterval.last = row[1]
             obj.positions.addPositionXY(row[2],row[3])
 
-    if obj:
+    if obj != None and obj.length() == obj.positions.length():
         objects.append(obj)
+    elif obj != None:
+        print('Object {} is missing {} positions'.format(obj.getNum(), int(obj.length())-obj.positions.length()))
 
     return objects
 
-def loadTrajectoriesFromSqlite(filename, trajectoryType, objectNumbers = -1):
-    '''Loads nObjects or the indices in objectNumbers from the database'''
-    connection = sqlite3.connect(filename) # add test if it open
+def loadUserTypesFromTable(cursor, trajectoryType, objectNumbers):
+    objectIdQuery = getTrajectoryIdQuery(objectNumbers, trajectoryType)
+    if objectIdQuery == '':
+        cursor.execute('SELECT object_id, road_user_type from objects')
+    else:
+        cursor.execute('SELECT object_id, road_user_type from objects where '+objectIdQuery[7:])
+    userTypes = {}
+    for row in cursor:
+        userTypes[row[0]] = row[1]
+    return userTypes
+
+def loadTrajectoriesFromSqlite(filename, trajectoryType, objectNumbers = None):
+    '''Loads the first objectNumbers objects or the indices in objectNumbers from the database'''
+    connection = sqlite3.connect(filename)
 
     objects = loadTrajectoriesFromTable(connection, 'positions', trajectoryType, objectNumbers)
     objectVelocities = loadTrajectoriesFromTable(connection, 'velocities', trajectoryType, objectNumbers)
@@ -384,26 +408,40 @@
                 obj.featureNumbers = featureNumbers[obj.getNum()]
 
             # load userType
-            if objectIdQuery == '':
-                cursor.execute('SELECT object_id, road_user_type from objects')
-            else:
-                cursor.execute('SELECT object_id, road_user_type from objects where '+objectIdQuery[7:])
-            userTypes = {}
-            for row in cursor:
-                userTypes[row[0]] = row[1]
-            
+            userTypes = loadUserTypesFromTable(cursor, trajectoryType, objectNumbers)
             for obj in objects:
                 obj.userType = userTypes[obj.getNum()]
              
         except sqlite3.OperationalError as error:
             printDBError(error)
-            return []
+            objects = []
 
     connection.close()
     return objects
 
-def removeFromSqlite(filename, dataType):
-    'Removes some tables in the filename depending on type of data'
+def loadGroundTruthFromSqlite(filename, gtType, gtNumbers = None):
+    'Loads bounding box annotations (ground truth) from an SQLite '
+    connection = sqlite3.connect(filename)
+    gt = []
+
+    if gtType == 'bb':
+        topCorners = loadTrajectoriesFromTable(connection, 'bounding_boxes', 'bbtop', gtNumbers)
+        bottomCorners = loadTrajectoriesFromTable(connection, 'bounding_boxes', 'bbbottom', gtNumbers)
+        userTypes = loadUserTypesFromTable(connection.cursor(), 'object', gtNumbers) # string format is same as object
+        
+        for t, b in zip(topCorners, bottomCorners):
+            num = t.getNum()
+            if t.getNum() == b.getNum():
+                annotation = moving.BBAnnotation(num, t.getTimeInterval(), t, b, userTypes[num])
+                gt.append(annotation)
+    else:
+        print ('Unknown type of annotation {}'.format(gtType))
+
+    connection.close()
+    return gt
+
+def deleteFromSqlite(filename, dataType):
+    'Deletes (drops) some tables in the filename depending on type of data'
     import os
     if os.path.isfile(filename):
         connection = sqlite3.connect(filename)
@@ -525,7 +563,7 @@
     connection.commit()
     connection.close()
 
-def loadBoundingBoxTable(filename):
+def loadBoundingBoxTableForDisplay(filename):
     connection = sqlite3.connect(filename)
     cursor = connection.cursor()
     boundingBoxes = {} # list of bounding boxes for each instant
@@ -534,9 +572,7 @@
         result = [row for row in cursor]
         if len(result) > 0:
             cursor.execute('SELECT * FROM bounding_boxes')
-            #objId = -1
             for row in cursor:
-                #if row[0] != objId:
                 boundingBoxes.setdefault(row[1], []).append([moving.Point(row[2], row[3]), moving.Point(row[4], row[5])])
     except sqlite3.OperationalError as error:
         printDBError(error)
@@ -544,6 +580,19 @@
     connection.close()
     return boundingBoxes
 
+def loadBoundingBoxTable(filename):
+    connection = sqlite3.connect(filename)
+    cursor = connection.cursor()
+    boundingBoxes = []
+    
+    try:
+        pass
+    except sqlite3.OperationalError as error:
+        printDBError(error)
+        return boundingBoxes
+    connection.close()
+    return boundingBoxes
+
 
 #########################
 # txt files
--- a/python/utils.py	Wed Dec 03 18:56:49 2014 +0000
+++ b/python/utils.py	Wed Dec 10 14:18:30 2014 -0500
@@ -61,15 +61,13 @@
     def nSamples(self):
         return sum(self.counts)
 
-def cumulativeDensityFunction(sample):
+def cumulativeDensityFunction(sample, normalized = False):
     '''Returns the cumulative density function of the sample of a random variable'''
-    from numpy.core.multiarray import array
-    from numpy.lib.function_base import unique
-    from numpy.core.fromnumeric import sum
-    a = array(sample)
-    a.sort()
-    xaxis = unique(a)
-    counts = [sum(a <= x) for x in xaxis]
+    from numpy import arange, cumsum
+    xaxis = sorted(sample)
+    counts = arange(1,len(sample)+1) # dtype = float
+    if normalized:
+        counts /= float(len(sample))
     return xaxis, counts
 
 class EmpiricalDiscreteDistribution(EmpiricalDistribution):
--- a/scripts/delete-tables.py	Wed Dec 03 18:56:49 2014 +0000
+++ b/scripts/delete-tables.py	Wed Dec 10 14:18:30 2014 -0500
@@ -5,10 +5,10 @@
 import utils
 import storage
 
-parser = argparse.ArgumentParser(description='The program deletes the tables in the database before saving new results (for objects, tables object_features and objects are dropped; for interactions, the tables interactions and indicators are dropped')
+parser = argparse.ArgumentParser(description='The program deletes (drops) the tables in the database before saving new results (for objects, tables object_features and objects are dropped; for interactions, the tables interactions and indicators are dropped')
 #parser.add_argument('configFilename', help = 'name of the configuration file')
 parser.add_argument('-d', dest = 'databaseFilename', help = 'name of the Sqlite database', required = True)
 parser.add_argument('-t', dest = 'dataType', help = 'type of the data to remove', required = True, choices = ['object','interaction', 'bb'])
 args = parser.parse_args()
 
-storage.removeFromSqlite(args.databaseFilename, args.dataType)
+storage.deleteFromSqlite(args.databaseFilename, args.dataType)
--- a/scripts/display-trajectories.py	Wed Dec 03 18:56:49 2014 +0000
+++ b/scripts/display-trajectories.py	Wed Dec 10 14:18:30 2014 -0500
@@ -63,5 +63,5 @@
 
 
 objects = storage.loadTrajectoriesFromSqlite(databaseFilename, args.trajectoryType)
-boundingBoxes = storage.loadBoundingBoxTable(databaseFilename)
+boundingBoxes = storage.loadBoundingBoxTableForDisplay(databaseFilename)
 cvutils.displayTrajectories(videoFilename, objects, boundingBoxes, homography, firstFrameNum, args.lastFrameNum, rescale = args.rescale, nFramesStep = args.nFramesStep, saveAllImages = args.saveAllImages, undistort = (undistort or args.undistort), intrinsicCameraMatrix = intrinsicCameraMatrix, distortionCoefficients = distortionCoefficients, undistortedImageMultiplication = undistortedImageMultiplication)