Understanding FFT- based algorithm to calculate image displacements with IDL programming language

1. Understanding FFT- based algorithm to calculate image displacements with IDL programming language By Cynthia Rodriguez University of Texas at San Antonio

2. Automatic image registration • Automatic image registration - is the process that aligns one image to a second image of the same scene • One image is known as the reference data, and the other image, is known as the input data (or sensed data), whichis matched relative to the reference data. • This is done so that the properties of an object being imaged can be addressable by the same coordinate pair from either one of the images.

3. FFT- based algorithm • Several different algorithms to solve automatic image registration • FFT stands for Fast Fourier Transform • Formula proven as a successful factor in solving automatic image registration

4. The Basic Idea Images I1 and I2 Differ by a shift (x0,y0) FFT based automatic registration relies on the Fourier shift theorem, which guarantees that the phase of a specifically defined “ratio” is equal to the phase difference between the images. Ex: Ratio = (x0,y0), [I2(x,y) = I1(x-x0, y-y0)] To find displacement between two images Compute the Ratio Ratio = F1 conj(F2 ) / | F1F2 | And Apply Inverse Fourier transform

5. The Basic Idea Images I1 and I2 Differ by a shift (x0,y0) X 0 1 2 Y 0 1 2 0 0 0 0 0 0 0.6 1 0.3 By applying the Inverse Fourier transform to the Ratio, we get an array of numbers returned that is zero everywhere except for a small area around a single point. By using the Max function we can find the maximum value The location of the Max value is exactly the displacement (x0,y0) that is needed to optimally register the images. Max value location (1, 2)

6. Pseudocode for shift_idl.pro • Declare and Initiate variables • ‘n’ representing columns • ‘m’ representing lines or rows • Create a byte array • Open and read image 1 file into im1 variable • Call function FFT for im1 • Repeat above steps for im2 • Calculate ratio F1 conj(F2 ) / | F1F2 | • Call function inverse FFT on ratio

7. Cont. Pseudocode for shift_idl.pro The Result from FFT • Function Inverse FFT returns an array, with approximate zeros everywhere except at the displacement that is used to register the images. • Call Max function to get the position of max value in array (the location of max is where the displacement (x0, y0) is) • To find x and y coordinate of shift, (x0, y0), use MOD and Divide operator • Last make a conditional statement to check for shifts to the west and north that are negative numbers. While shifts to the east and south are positive.

8. Rotating and Scaling Too! To find new value for M(x,y) • Convert abs(F(ε,n)) from rectangular coordinates (x, y) to log-polar coordinates (log (p), θ) to represent both scaling and rotation as shifts. • However computing (log (p), θ) for the rectangular grid coordinates leads to points that are not located exactly at points on the original grid. So interpolation is used to findthe correct value on the new grid. Bilinear interpolation is used in this program. X 0 1 2 Y 0 1 2 0 i 0 ij i 0.6 i 0.3 Bilinear Interpolation Get intensities surrounding around M(x,y) j then interpolate

9. Pseudocode for srs_idl.pro • Declare and initiate variables • Call function Loadimage ( opens and reads data into Im1 and Im2) • Call function Normalize ( gets grid values in array and makes the array a floating array and divides values by 255.0 to make values between 0 and 1) • Call procedure ShiftArr (makes values in array negative values, if divisible by 2 with a remainder equal to 1) • Call FFT (which is a built-in IDL function) does Fast Fourier transfrom on each image • Call abs (also built-in function) gets absolute values on each image • Call function Highpass (removes low frequency noise) to the absolute values.

10. Cont. Pseudocode for srs_idl.pro • Call function LogPolar on im1 and im2 (Transforms rectangular coordinates to log-polar coordinates.) • To carry this function out accurately, the polar plane must have the same number of rows as the rectangular plane. To achieve this, create a nested for loop to go through each index in the array and individually get each coordinate and transfer it to log polar coordinate. • After getting coordinates, Apply Bilinear interpolation which gets surrounding values of (x0, y0) then interpolate as M(x,y)=Mjk(1-t)(1-u) + Mj+1,k t(1-u) + Mj,k+1(1-t)u + Mj+1,k+1tu, t is a fractional part of x and u is a fractional part of y

11. Cont. Pseudocode for srs_idl.pro • Call FFT function on polar coordinates for im1 and im2 • Call function Ritio to both fft_polar1 and fft_polar2 F1 conj(F2 ) / | F1F2 | • Call FFT inverse function for Ratio • Apply absolute function to inverse Ratio • Find position of Max number in array of absolute values • With Max position find the values for scale and angle scale=b^(position MOD rows) angle=180.0 * (position/rows) /cols

12. Cont. Pseudocode for srs_idl.pro • Once the angle is found make a conditional statement to check to see if the angle is greater than or less than 90.0. • A positive angle means that the image is rotated to the east relative to the base image (clockwise), • A negative angle means that the image is rotated to the west (counterclockwise) • Assume angle computed is between [0,180]. • If angle is less than 90, the actual angle is the negative of the computed one. • If angle is greater than or equal to 90, then the actual rotation angle is 180 – computed angle. • With angle and scale, create new image3 with function ROT (constructs new image by doing reverse rotation and reverse scaling)

13. Cont. Pseudocode for srs_idl.pro Next Steps like shift_idl.pro • Call FFT function for image1 and image3 • Compute Ratio • Call inverse FFT on Ratio • Find Max number of array and get position • Get coordiantes of position • Check to see if coordinates are out of bounds and adjust accordingly for north and west shifts which are computed as negative.

14. Results • ENVI> .COMPILE "D:\idlproject\shift_idl.pro" • Compiled module: SHIFT_IDL. • ENVI> shift_idl • I = 7050 • this is max = ( 92546.2, -0.000154680) • Position = 50 20 • ENVI> .COMPILE "D:\idlproject\srs_idl.pro" • Compiled module: RITIO. Compiled module: LOGPOLAR. Compiled module: HIGHPASS. Compiled module: COMBINEFFT. Compiled module: LOADIMAGE. Compiled module: NORMALIZE. Compiled module: SHIFTARR. Compiled module: SRS_IDL. • ENVI> srs_idl • Loading images... • Doing FFT on two images... • Getting absolute value... • Transformming log_polar coordinates... • Doing FFT on polar images... • max absolute value position= 143217 • total number of elements 160000 • max absolute value = 10522.9 • computed scale = 1.28999 • computed angle = 18.9000 • this is IR • I = 7603 • this is max = ( 14344.0, 5.01551e-005) • Position = 3 19

15. Summary • The shift only program contains two forward FFTs and one inverse FFT • The scale, rotation, and shift program contains six forward FFTs and two inverse FFTs and requires more time and memory. • These two programs provide Automatic Registration THE END