--- a/python/events.py	Tue Feb 17 02:21:31 2015 +0100
+++ b/python/events.py	Tue Feb 17 03:55:55 2015 +0100
@@ -200,13 +200,17 @@
         for i, cp in self.collisionPoints.iteritems():
             TTCs[i] = prediction.SafetyPoint.computeExpectedIndicator(cp)
         # add probability of collision, and probability of successful evasive action
-        self.addIndicator(indicators.SeverityIndicator(Interaction.indicatorNames[7], TTCs))
+        self.addIndicator(indicators.SeverityIndicator(Interaction.indicatorNames[7], TTCs, mostSevereIsMax=False))
         
         if computeCZ:
             pPETs = {}
             for i, cz in self.crossingZones.iteritems():
                 pPETs[i] = prediction.SafetyPoint.computeExpectedIndicator(cz)
-            self.addIndicator(indicators.SeverityIndicator(Interaction.indicatorNames[9], pPETs))
+            self.addIndicator(indicators.SeverityIndicator(Interaction.indicatorNames[9], pPETs, mostSevereIsMax=False))
+
+    def computePET(self, collisionDistanceThreshold):
+        # TODO add crossing zone
+        self.pet = moving.MovingObject.computePET(self.roadUser1, self.roadUser2, collisionDistanceThreshold)
 
     def addVideoFilename(self,videoFilename):
         self.videoFilename= videoFilename

--- a/python/indicators.py	Tue Feb 17 02:21:31 2015 +0100
+++ b/python/indicators.py	Tue Feb 17 03:55:55 2015 +0100
@@ -141,7 +141,7 @@
         TemporalIndicator.__init__(self, name, values, timeInterval, maxValue)
         self.mostSevereIsMax = mostSevereIsMax
 
-    def getMostSevereValue(self, minNInstants=1): # TODO use scoreatpercentile
+    def getMostSevereValue(self, minNInstants=1): # TODO use np.percentile
         from matplotlib.mlab import find
         from numpy.core.multiarray import array
         from numpy.core.fromnumeric import mean

--- a/python/moving.py	Tue Feb 17 02:21:31 2015 +0100
+++ b/python/moving.py	Tue Feb 17 03:55:55 2015 +0100
@@ -582,7 +582,7 @@
 def segmentLineIntersection(p1, p2, p3, p4):
     '''Indicates if the line going through p1 and p2 intersects inside p3, p4'''
     inter = intersection(p1, p2, p3, p4)
-    if inter != None and utils.inBetween(p3.x, p4.x, inter.x) and utils.inBetween(p3.y, p4.y, inter.y)):
+    if inter != None and utils.inBetween(p3.x, p4.x, inter.x) and utils.inBetween(p3.y, p4.y, inter.y):
         return inter
     else:
         return None
@@ -812,6 +812,7 @@
         intersects with the segment of extremities p1 and p2 
         the list is empty if there is no crossing'''
         indices = []
+        intersections = []
 
         for i in xrange(self.length()-1):
             q1=self.__getitem__(i)
@@ -825,15 +826,16 @@
                 else:
                     ratio = 0
                 indices.append(i+ratio)
+                intersections.append(p)
         return indices
 
     def getLineIntersections(self, p1, p2):
-         '''Returns a list of the indices at which the trajectory 
+        '''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 = []
         intersections = []
-
+        
         for i in xrange(self.length()-1):
             q1=self.__getitem__(i)
             q2=self.__getitem__(i+1)
@@ -1169,9 +1171,29 @@
         '''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)
+        indices, intersections = self.positions.getIntersections(p1, p2)
         return [t+self.getFirstInstant() for t in indices]
 
+    @staticmethod
+    def computePET(obj1, obj2, collisionDistanceThreshold):
+        '''Post-encroachment time based on distance threshold'''
+        #for i in xrange(int(obj1.length())-1):
+        #    for j in xrange(int(obj2.length())-1):
+        #        inter = segmentIntersection(obj1.getPositionAt(i), obj1.getPositionAt(i+1), obj2.getPositionAt(i), obj2.getPositionAt(i+1))
+        from scipy.spatial.distance import cdist
+        from numpy import zeros
+        positions1 = [p.astuple() for p in obj1.getPositions()]
+        positions2 = [p.astuple() for p in obj2.getPositions()]
+        pets = zeros((int(obj1.length()), int(obj2.length())))
+        for i,t1 in enumerate(obj1.getTimeInterval()):
+            for j,t2 in enumerate(obj2.getTimeInterval()):
+                pets[i,j] = abs(t1-t2)
+        distances = cdist(positions1, positions2, metric = 'euclidean')
+        if distances.min() <= collisionDistanceThreshold:
+            return pets[distances <= collisionDistanceThreshold].min()
+        else:
+            return None
+
     def predictPosition(self, instant, nTimeSteps, externalAcceleration = Point(0,0)):
         '''Predicts the position of object at instant+deltaT, 
         at constant speed'''

--- a/python/tests/moving.txt	Tue Feb 17 02:21:31 2015 +0100
+++ b/python/tests/moving.txt	Tue Feb 17 03:55:55 2015 +0100
@@ -53,9 +53,9 @@
 >>> Point.distanceNorm2(Point(3,4),Point(1,7))
 3.605551275463989
 
->>> Point(3,2).inPolygonNoShapely(np.array([[0,0],[1,0],[1,1],[0,1]]))
+>>> Point(3,2).inPolygon(np.array([[0,0],[1,0],[1,1],[0,1]]))
 False
->>> Point(3,2).inPolygonNoShapely(np.array([[0,0],[4,0],[4,3],[0,3]]))
+>>> Point(3,2).inPolygon(np.array([[0,0],[4,0],[4,3],[0,3]]))
 True
 
 >>> predictPositionNoLimit(10, Point(0,0), Point(1,1)) # doctest:+ELLIPSIS
@@ -84,10 +84,6 @@
 True
 >>> t1[1]
 (1.500000,3.500000)
->>> t1.getTrajectoryInPolygonNoShapely(np.array([[0,0],[4,0],[4,3],[0,3]]))
-(0.500000,0.500000)
->>> t1.getTrajectoryInPolygonNoShapely(np.array([[10,10],[14,10],[14,13],[10,13]])).length()
-0
 
 >>> t1.differentiate()
 (1.000000,3.000000) (1.000000,3.000000)
@@ -141,6 +137,14 @@
 >>> Point.midPoint(p1, p2)
 (0.500000,0.500000)
 
+>>> o1 = MovingObject.generate(Point(-5.,0.), Point(1.,0.), TimeInterval(0,10))
+>>> o2 = MovingObject.generate(Point(0.,-5.), Point(0.,1.), TimeInterval(0,10))
+>>> MovingObject.computePET(o1, o2, 0.1)
+0.0
+>>> o2 = MovingObject.generate(Point(0.,-5.), Point(0.,1.), TimeInterval(5,15))
+>>> MovingObject.computePET(o1, o2, 0.1)
+5.0
+
 >>> t = CurvilinearTrajectory(S = [1., 2., 3., 5.], Y = [0.5, 0.5, 0.6, 0.7], lanes = ['1']*4)
 >>> t.differentiate() # doctest:+ELLIPSIS
 [1.0, 0.0, '1'] [1.0, 0.099..., '1'] [2.0, 0.099..., '1']