Source Code for Module pyvision.surveillance.optical_flow

  1  # PyVision License 
  2  # 
  3  # Copyright (c) 2006-2008 David S. Bolme 
  4  # All rights reserved. 
  5  # 
  6  # Redistribution and use in source and binary forms, with or without 
  7  # modification, are permitted provided that the following conditions 
  8  # are met: 
  9  #  
 10  # 1. Redistributions of source code must retain the above copyright 
 11  # notice, this list of conditions and the following disclaimer. 
 12  #  
 13  # 2. Redistributions in binary form must reproduce the above copyright 
 14  # notice, this list of conditions and the following disclaimer in the 
 15  # documentation and/or other materials provided with the distribution. 
 16  #  
 17  # 3. Neither name of copyright holders nor the names of its contributors 
 18  # may be used to endorse or promote products derived from this software 
 19  # without specific prior written permission. 
 20  #  
 21  #  
 22  # THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 
 23  # ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 
 24  # LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR 
 25  # A PARTICULAR PURPOSE ARE DISCLAIMED.  IN NO EVENT SHALL THE REGENTS OR 
 26  # CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
 27  # EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
 28  # PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 
 29  # PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 
 30  # LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
 31  # NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
 32  # SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 
 33   
 34  ''' 
 35  Created on Dec 10, 2009 
 36   
 37  @author: bolme 
 38  ''' 
 39  import pyvision as pv 
 40  import cv 
 41  import numpy as np 
 42   
 43   
 44 -class OpticalFlow: 
 45      ''' 
 46      This class tracks the motion of the camera.  This can be used for things like 
 47      video stablization or to better understand camera motion across multiple 
 48      frames. 
 49       
 50      This uses the LK optical flow algorithm to track point correspondences from one 
 51      frame to another.   
 52       
 53      It can output the frame to frame camera motion as a homography. 
 54      ''' 
 55       
 56 -    def __init__(self,tile_size="ORIG"): 
 57          ''' 
 58          Create the camera tracker 
 59           
 60          @param tile_size: Frames will be resized to these dimensions before tracking. 
 61          ''' 
 62          self.frame_size = None 
 63          self.prev_frame = None 
 64          self.tile_size = tile_size 
 65          self.tracks = [] 
 66          self.bad_points = [] 
 67           
 68          # Divide the image in to gridXgrid regions 
 69          self.grid = 5 
 70           
 71          # Require a minimum number of points_b per grid 
 72          self.min_points = 5 
 73           
 74          self.n = 0 
 75           
 76          self.homography = None 
 77          self.homography_rev = None 
 78           
 79          self.prev_input = None 
 80           
 81 -    def update(self,frame): 
 82          ''' 
 83          This tracks the points_b to the next frame using the LK tracker. 
 84          Add more good points_b to track. 
 85           
 86          @param frame: update optical flow for the frame. 
 87          @type frame: pv.Image 
 88          ''' 
 89          if self.frame_size == None: 
 90              self.frame_size = frame.size 
 91               
 92          if self.prev_frame == None: 
 93              #Initialize the tracker 
 94               
 95              # Pick a frame size 
 96              if self.tile_size == "AUTO": 
 97                  # Rescale the image so that the largest demenion is between 256 and 512 
 98                  w,h = frame.size 
 99                  n = max(w,h) 
100                  while n > 512: 
101                      n = n/2 
102                   
103                  scale = n/float(max(w,h)) 
104                   
105                  w = int(round(scale*w)) 
106                  h = int(round(scale*h)) 
107                   
108                  self.tile_size = (w,h) 
109                   
110                   
111              if self.tile_size == "ORIG": 
112                  self.tile_size = frame.size 
113   
114              # allocate memory 
115              w,h = self.tile_size 
116              self.frame = cv.CreateImage((w,h), 8, 1) 
117              self.prev_frame = cv.CreateImage((w,h), 8, 1)  
118               
119              # Resize the frame 
120              cvim = frame.asOpenCVBW() 
121              cv.Resize(cvim,self.frame) 
122               
123              # Find good features to track 
124              self._selectTrackingPoints(self.frame) 
125               
126              # Init transform 
127               
128               
129          else: 
130              # Resize the frame 
131              cvim = frame.asOpenCVBW() 
132              cv.Resize(cvim,self.frame) 
133           
134           
135              # Track the points_b (optical flow) 
136              new_tracks = self._opticalFlow() 
137           
138               
139              # Compute the transform 
140              assert len(new_tracks) == len(self.tracks) 
141              n = len(new_tracks) 
142              src_points = cv.CreateMat(n,2,cv.CV_32F) 
143              dst_points = cv.CreateMat(n,2,cv.CV_32F) 
144       
145              # Copy data into matricies 
146              for i in range(len(new_tracks)): 
147                  src_points[i,0] = self.tracks[i].X() 
148                  src_points[i,1] = self.tracks[i].Y() 
149   
150                  dst_points[i,0] = new_tracks[i].X() 
151                  dst_points[i,1] = new_tracks[i].Y() 
152       
153              # Create homography matrix 
154              self.homography = cv.CreateMat(3,3,cv.CV_32F) 
155              self.homography_rev = cv.CreateMat(3,3,cv.CV_32F) 
156              mask = cv.CreateMat(1,n,cv.CV_8U) 
157              cv.FindHomography(dst_points,src_points,self.homography_rev,cv.CV_LMEDS,1.5,mask) 
158              cv.FindHomography(src_points,dst_points,self.homography,cv.CV_LMEDS,1.5,mask) 
159               
160              # Drop bad points_b 
161              self.tracks = [] 
162              self.bad_points = [] 
163              for i in range(n): 
164                  if mask[0,i]: 
165                      self.tracks.append(new_tracks[i]) 
166                  else: 
167                      self.bad_points.append(new_tracks[i]) 
168                       
169              self.n = len(self.tracks) 
170               
171              # Add new points_b 
172              self._selectTrackingPoints(self.frame) 
173               
174              self.prev_input.to_next = self.asHomography(forward = False) 
175              frame.to_prev = self.asHomography(forward = True) 
176               
177          self.prev_input = frame 
178   
179          cv.Copy(self.frame,self.prev_frame) 
180               
181 -    def _selectTrackingPoints(self,frame): 
182          ''' 
183          This uses the OpenCV get good features to track to initialize a set of tracking points_b. 
184          ''' 
185          quality = 0.01 
186          min_distance = 15 
187           
188          w,h = self.tile_size 
189          tw = w//self.grid 
190          th = h//self.grid 
191   
192          for i in range(self.grid): 
193              for j in range(self.grid): 
194                  ul = pv.Point(i*tw,j*th) 
195                  rect = pv.Rect(i*tw,j*th,tw,th) 
196                  count = 0 
197                  for pt in self.tracks: 
198                      if rect.containsPoint(pt): 
199                          count += 1 
200                       
201                  if count < self.min_points: 
202                      gray = cv.CreateImage ((tw,th), 8, 1) 
203                       
204                      faceim = cv.GetSubRect(frame, rect.asOpenCV()) 
205                      cv.Resize(faceim,gray) 
206   
207                      eig = cv.CreateImage ((tw,th), 32, 1) 
208                      temp = cv.CreateImage ((tw,th), 32, 1) 
209                   
210                      # search the good points_b 
211                      points_b = cv.GoodFeaturesToTrack (gray, eig, temp, 2*self.min_points, quality, min_distance, None, 3, 0, 0.04) 
212                       
213                      for pt in points_b: 
214                          self.tracks.append(ul+pv.Point(pt)) 
215                                           
216           
217 -    def _opticalFlow(self): 
218          ''' 
219          Compute the optical flow between frames using cv.CalcOpticalFlow 
220           
221          @returns: a list of tracks for the new image 
222          @rtype: list of pv.Point() 
223          ''' 
224          flags = 0 
225           
226          grey = self.frame 
227          prev_grey = self.prev_frame 
228       
229          pyramid = cv.CreateImage (cv.GetSize (grey), 8, 1) 
230          prev_pyramid = cv.CreateImage (cv.GetSize (grey), 8, 1) 
231       
232          cv_points = [] 
233          for each in self.tracks: 
234              cv_points.append((each.X(),each.Y())) 
235           
236          points_b, _, _,= cv.CalcOpticalFlowPyrLK ( 
237                      prev_grey,  
238                      grey,  
239                      prev_pyramid,  
240                      pyramid, 
241                      cv_points, 
242                      (5,5), 
243                      3,#pyr number 
244                      (cv.CV_TERMCRIT_ITER | cv.CV_TERMCRIT_EPS, 10, 0.01), 
245                      flags) 
246                   
247          result = [] 
248          for pt in points_b: 
249              result.append(pv.Point(pt)) 
250       
251          return result 
252   
253                   
254 -    def asHomography(self,forward=True): 
255          ''' 
256          Get the transformation as a homography. 
257           
258          @keyword forward: switch between the forward and reverse transform. 
259          @ktype forward: True|False 
260          ''' 
261          fw,fh = self.frame_size 
262          tw,th = self.frame_size 
263          if forward: 
264              if self.homography == None: 
265                  matrix = np.eye(3) 
266              else: 
267                  matrix = np.linalg.inv(pv.OpenCVToNumpy(self.homography)) 
268              matrix = np.dot(np.diag([tw/fw,th/fh,1.0]),matrix) 
269              perspective = pv.PerspectiveTransform(matrix,self.frame_size) 
270          else: 
271              if self.homography_rev == None: 
272                  matrix = np.eye(3) 
273              else: 
274                  matrix = np.linalg.inv(pv.OpenCVToNumpy(self.homography_rev)) 
275              matrix = np.dot(np.diag([tw/fw,th/fh,1.0]),matrix) 
276              perspective = pv.PerspectiveTransform(matrix,self.frame_size) 
277               
278          return perspective 
279   
280           
281 -    def annotateFrame(self,frame,color='white'): 
282          ''' 
283          Renders optical flow information to the frame. 
284           
285          @param frame: the frame that will be annotated 
286          @type frame: pv.Image 
287          @keyword type: 
288          @ktype type: "TRACKING"  
289           
290          ''' 
291          w,h = self.tile_size 
292          rect = pv.Rect(0,0,w,h) 
293          affine = pv.AffineFromRect(rect,frame.size) 
294          for pt in affine.transformPoints(self.tracks[:self.n]): 
295              frame.annotatePoint(pt,color=color) 
296   
297          for pt in affine.transformPoints(self.tracks[self.n:]): 
298              frame.annotatePoint(pt,color='red') 
299               
300          for pt in affine.transformPoints(self.bad_points): 
301              frame.annotatePoint(pt,color='gray') 
302               
303           
304          if self.homography != None: 
305              matrix = pv.OpenCVToNumpy(self.homography) 
306              perspective = pv.PerspectiveTransform(matrix,frame.size) 
307   
308              for pt in self.tracks[:self.n]: 
309                  # Transform the point multiple times to amplify the track  
310                  # for visualization 
311                  old = perspective.invertPoint(pt) 
312                  old = perspective.invertPoint(old) 
313                  old = perspective.invertPoint(old) 
314                  frame.annotateLine(pt,old,color=color) 
315