Source Code for Module pyvision.point.PhaseCorrelation

  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  __author__ = "$Author: bolme $" 
 35  __version__ = "$Rev: 625 $" 
 36  __info__ = "$Id: PhaseCorrelation.py 625 2008-03-24 20:24:37Z bolme $" 
 37  __copyright__ = "Copyright 2006 David S Bolme" 
 38   
 39  import numpy as np 
 40  import pyvision as pv 
 41   
 42   
 43 -def PhaseCorrelation(tile1, tile2, phase_only=True, ilog=None): 
 44      ''' 
 45      Uses phase correlation to estimate the best integer displacement to align the images. 
 46      Also fits a quadradic to the correltaion surface to determine a sub pixel estimate of  
 47      displacement. 
 48       
 49      Returns four values as a tuple: max_corr, max_displacement, est_corr, est_displacement 
 50      max_corr         - maximum correlation value. 
 51      max_displacement - displacement needed to obtain the maximum correlation. (full pixel) 
 52      est_corr         - estimated corelation value if subpixel displacement is used. 
 53      est_displacement - estimated displacement (subpixel) 
 54       
 55      see http://en.wikipedia.org/wiki/Phase_correlation 
 56      ''' 
 57      if isinstance(tile1,pv.Image): 
 58          tile1 = tile1.asMatrix2D() 
 59      else: 
 60          tile1 = pv.Image(tile1).asMatrix2D() 
 61          raise TypeError("Please pass in a numpy array or a pyvision image.") 
 62   
 63      if isinstance(tile2,pv.Image): 
 64          tile2 = tile2.asMatrix2D() 
 65      else: 
 66          tile2 = pv.Image(tile2).asMatrix2D()     
 67       
 68      if tile1.shape != tile2.shape: 
 69          raise ValueError("Image tiles must have the same shape. [tile1 %s] != [tile2 %s]"%(tile1.shape,tile2.shape)) 
 70   
 71      # copy the data 
 72      tile1 = tile1.copy() 
 73      tile2 = tile2.copy() 
 74   
 75      # normalize the image tiles 
 76      tile1 = pv.meanUnit(tile1) 
 77      tile2 = pv.meanUnit(tile2) 
 78   
 79      # compute the fft 
 80      Ia = np.fft.fft2(tile1) 
 81      Ib = np.conjugate(np.fft.fft2(tile2)) 
 82   
 83      # build the normalized cross-power spectrum 
 84      ncs = Ia*Ib 
 85   
 86      if phase_only: 
 87          ncs = ncs/np.abs(ncs) 
 88   
 89      # build the power spectrum 
 90      pc = np.real(np.fft.ifft2(ncs)) 
 91       
 92      if ilog != None: 
 93          ilog.log(pv.Image(tile1),label="Tile1") 
 94          ilog.log(pv.Image(tile2),label="Tile2") 
 95          ilog.log(pv.Image(np.fft.fftshift(pc)),label="Correlation") 
 96   
 97      max_corr = pc.max() 
 98      max_elem = (pc == max_corr).nonzero() 
 99      max_elem = max_elem[0][0],max_elem[1][0] 
100   
101      max_point = list(max_elem) 
102      if max_elem[0]*2 > tile1.shape[0]: 
103          max_point[0] = -tile1.shape[0] + max_elem[0] 
104      if max_elem[1]*2 > tile1.shape[1]: 
105          max_point[1] = -tile1.shape[1] + max_elem[1] 
106   
107      est_corr, est_point = QuadradicEstimation(pc,max_point) 
108   
109      return max_corr, max_point, est_corr, est_point 
110   
111   
112   
113 -def QuadradicEstimation(pc, max_point): 
114      ''' 
115      Estimate the subpixel displacement based on fitting a quadradic surface 
116      to area surrounding the max point. 
117   
118      returns the estimated power spectrum at the subpixel point, and the displacement. 
119      ''' 
120      # fit a surface to the area surrounding the point 
121      # ax2 + by2 + cx + dy + e = z 
122      def create_row(x,y): 
123          return [x*x, y*y, x, y, 1] 
124   
125      A = [] 
126      b = [] 
127   
128      x = max_point[0] 
129      y = max_point[1] 
130      A.append(create_row(x,y)) 
131      b.append([pc[x,y]]) 
132   
133      x = max_point[0]+1 
134      y = max_point[1] 
135      A.append(create_row(x,y)) 
136      b.append([pc[x,y]]) 
137   
138      x = max_point[0]-1 
139      y = max_point[1] 
140      A.append(create_row(x,y)) 
141      b.append([pc[x,y]]) 
142   
143      x = max_point[0] 
144      y = max_point[1]+1 
145      A.append(create_row(x,y)) 
146      b.append([pc[x,y]]) 
147   
148      x = max_point[0] 
149      y = max_point[1]-1 
150      A.append(create_row(x,y)) 
151      b.append([pc[x,y]]) 
152   
153      A = np.array(A) 
154      b = np.array(b) 
155   
156      lss = np.linalg.lstsq(A,b) 
157      assert lss[2] == 5 #make sure the matrix was of full rank 
158      x = lss[0] 
159   
160      a = x[0,0] 
161      b = x[1,0] 
162      c = x[2,0] 
163      d = x[3,0] 
164      e = x[4,0] 
165   
166      # find the derivitive with respect to x and y solve for f`(x,y) = 0 to find the max 
167      # 2ax + c = 0 
168      # 2by + d = 0 
169      x = -c/(2*a) 
170      y = -d/(2*b) 
171   
172      # estimate the max at subpixel resolution 
173      emax = a*x*x + b*y*y + c*x + d*y + e 
174   
175      return emax, [x,y] 
176       
177  # TODO: Add UnitTests 
178  import unittest 
179   
180 -class _TestPhaseCorrelation(unittest.TestCase): 
181       
182 -    def test_Correlation(self): 
183          import os.path 
184   
185          filename_a = os.path.join(pv.__path__[0],'data','test','registration1a.jpg') 
186          filename_b = os.path.join(pv.__path__[0],'data','test','registration1b.jpg') 
187           
188          im_a = pv.Image(filename_a) 
189          im_b = pv.Image(filename_b) 
190   
191          out = PhaseCorrelation(im_a,im_b,phase_only=False) 
192           
193          self.assertAlmostEqual(out[0],0.87382160686468002,places=3) 
194          self.assertEqual(out[1][0],20) 
195          self.assertEqual(out[1][1],20) 
196          self.assertAlmostEqual(out[2],0.87388092414032315,places=3) 
197          self.assertAlmostEqual(out[3][0],19.978881341260816,places=3) 
198          self.assertAlmostEqual(out[3][1],19.986396178942329,places=3) 
199           
200 -    def test_PhaseCorrelation(self): 
201          import os.path 
202   
203          filename_a = os.path.join(pv.__path__[0],'data','test','registration1a.jpg') 
204          filename_b = os.path.join(pv.__path__[0],'data','test','registration1b.jpg') 
205           
206          im_a = pv.Image(filename_a) 
207          im_b = pv.Image(filename_b) 
208   
209          out = PhaseCorrelation(im_a,im_b,phase_only=False) 
210           
211          self.assertAlmostEqual(out[0],0.87382160686468002,places=3) 
212          self.assertEqual(out[1][0],20) 
213          self.assertEqual(out[1][1],20) 
214          self.assertAlmostEqual(out[2],0.87388092414032315,places=3) 
215          self.assertAlmostEqual(out[3][0],19.978881341260816,places=3) 
216          self.assertAlmostEqual(out[3][1],19.986396178942329,places=3) 
217