title: Image Registration - Phase Correlation Algorithm using AForge.NET Framework (C#)
body: <p>I've been writing a code in C# to register the shift (x, y axis) between two images with this Phase Correlation algorithm, as shown in Wikipedia here: <a href="http://en.wikipedia.org/wiki/Phase_correlation" rel="nofollow">Phase Correlation</a>.</p>
  
  <p>To help me with that, I'm using AForge.NET Framework's ComplexImage class and its Fast Fourier Transform methods.</p>
  
  <p>Code is shown below:</p>
  
  <pre><code> Bitmap fixed = (Bitmap)System.Drawing.Image.FromFile("Lenna.png");
   Bitmap shifted = (Bitmap)System.Drawing.Image.FromFile("Lenna2.png");
  
   //apply greyscale
   fixed = Grayscale.CommonAlgorithms.BT709.Apply(fixed);
   shifted = Grayscale.CommonAlgorithms.BT709.Apply(shifted);
  
   ComplexImage fixedCplx = ComplexImage.FromBitmap(fixed);
   int height = fixedCplx.Data.GetLength(0);
   int width  = fixedCplx.Data.GetLength(1);
   ComplexImage shiftedCplx = ComplexImage.FromBitmap(shifted);
  
   fixedCplx.ForwardFourierTransform();
   shiftedCplx.ForwardFourierTransform();
  
   for (int y = 0; y &lt; height; y++)
   {
       for (int x = 0; x &lt; width; x++)
       {
           //Calculating elementwise complex conjugate of the shifted image 2d vector
           shiftedCplx.Data[y, x].Im = -1 * shiftedCplx.Data[y, x].Im;
           //Elementwise multiplication to obtain cross power spectrum
           shiftedCplx.Data[y, x] = Complex.Multiply(fixedCplx.Data[y, x], shiftedCplx.Data[y, x]);
           //Elementwise normalization
           shiftedCplx.Data[y, x] = Complex.Divide(shiftedCplx.Data[y, x], shiftedCplx.Data[y, x].Magnitude);
       }
   }
   shiftedCplx.BackwardFourierTransform();
  </code></pre>
  
  <p>Yet, the peak function that is supposed to appear when the inverse transform is calculated does not show itself... I have tested with the pictures on that Wikipedia link but all I get as a result is a white picture. Is there anything wrong with the code shown above?</p>
  
  <p>Some considerations:</p>
  
  <ul>
  <li>Perhaps I'm not understanding the math behind the algorithm. When taking the elementwise product [Ga].[Gb*], does this * symbol mean I need to conjugate each element of the matrix as I've been doing? Or is it the <a href="http://en.wikipedia.org/wiki/Conjugate_transpose" rel="nofollow">conjugate transpose</a> of the matrix?</li>
  <li>When normalizing [Ga].[Gb*]/|[Ga].[Gb*]|, is it correct to divide each matrix element (which are all complex numbers) by its magnitude? I've done this in my code and as expected, by the end of the algorithm each element in the <code>shiftedCplx.Data[]</code> array has magnitude = 1 (this means all relevant info at this point is their phase). </li>
  </ul>
  
  <p>It looks like the ´BackwardFourierTransform()` function isn't really working as intended, or somewhere along the way I've mangled the data.    </p>
  
  <p><strong>EDIT 1</strong></p>
  
  <p>Just edited the code as jaket pointed out in the comment.</p>
  
  <p>Also, one thing I just realized: When Fourier transforming an greyscale image, each value in the corresponding 2D vector has a range of 0 to 255. For me, it was implied that the inverse transform would yield some value between that range as well (even if data was normalized beforehand), but I discovered it is not so in this case. Most elements in the <code>shiftedCplx.Data[,]</code> vector would have values larger than 255 after being transformed back into the spatial domain, and of course, making a bitmap image out of it creates an image thats is mostly all white. Is there something I'm missing?</p>