Registration of pedobarographic data sets in frequency domain

The assumption behind the cross-correlation-based method is that two images are optimally registered when their cross-correlation is maximized. If one considers two functions fand gand their cross-correlation: then, the cross-correlation can be given in function of a shift aby: This equation can be written as a convolution: where . From the convolution theorem, one can obtain: Thus, computing the inverse of the Fourier transform of the product of the last equation, the crosscorrelation can be obtained for all shifts. Then, the coordinates of the point which has the higher value represent the desired optimal integer shift. Francisco P. M. Oliveira1, Todd C. Pataky2, João Manuel R. S. Tavares1 1Faculdade de Engenharia da Universidade do Porto / Instituto de Engenharia Mecânica e Gestão Industrial -Portugal 2Shinshu University -Japan ([email protected], [email protected], t[email protected]) REGISTRATION OF PEDOBAROGRAPHIC DATA SETS IN FREQUENCY DOMAIN Plantar pressure distribution conveys significant information for clinicians and researchers about the structure and function of the foot, general mechanics of gait, and is a helpful means to evaluate patients with foot complaints. Image registration, i.e. the process of optimally aligning homologous structures represented in images, can be very useful for clinicians and researchers, since tasks such as the identification of the main plantar pressure areas and classification of the foot type can be done automatically.Image registration also allows that clinicians accurately compare the plantar pressure of a patient over the time or with an atlas.In addition, pedobarographic image registration supports pixel-level statistics, which makes possible the extraction of biomechanically-relevant information from plantar pressure images more effectively than traditional regional techniques. Several methods had been developed to register pedobarographic images;for instance, principal axes, modal matching, contours’ matching and based on the optimization of image similarity measure. In this work, three image registration algorithms are presented, all using the Fourier transform and its properties. One of the algorithms determines the optimal alignment by maximizing the crosscorrelation (CC), another one by minimizing the sum of squared differences (SSD) and the last one using the phase correlation technique. 1 INTRODUCTION 2.1 Fourier transform 4 CONCLUSIONS The experimental results validate the methodologies presented.The accuracy using a control geometric transformation was awfully good for all algorithms, even in the presence of noise. Comparing the MSE results obtained using unknown transformations with the MSE reported in other studies, one can conclude that the accuracy is slight inferior to the best accuracy reported for the same data set.However, the processing time is several times inferior. Using the PC-based methodology, the MSE value is not as good as using CC or SSD-based methodologies. This fact was already expected because MSE is anormalization of SSD. Based on the accuracy, processing time and robustness to noise, one can conclude that the methodologies proposed are suitable for implementation in near real-time applications in clinics and laboratories. REFERENCES Castro, E.; Morandi, C. (1987). Registration of translated and rotated images using finite Fourier transforms.IEEE Transactions on Pattern Analysis and Machine Intelligence 9(5):700-703. Oliveira, F.; Tavares, J.M.R.S.; Pataky, T. (2009). Rapid pedobarographic image registration based on contour curvature and optimization. Journal of Biomechanics 42(15):2620-2623. Oliveira, F.; Pataky, T.; Tavares, J.M.R.S. (2010). Registration of pedobarographic image data in the frequency domain. Computer Methods in Biomechanics and Biomedical Engineering, DOI: 10.1080/10255840903573020 (in press). Pataky, T.; Goulermas, J.; Crompton, R. (2008). A comparison of seven methods of within-subjects rigid-body pedobarographic image registration. Journal of Biomechanics 41:3085-3089. Reddy, B.; Chatterji, B. (1996). An FFT-based technique for translation, rotation, and scaleinvariant image registration.IEEE Transactions on Image Processing 5(8):1266-1271. Acknowledgements: The first author would like to thank his PhD grant from Fundação Calouste Gulbenkian in Portugal. This work was partially done under the scope of projects PTDC/EEA-CRO/103320/2008, UTAustin/MAT/0009/2008 and UTAustin/CA/0047/2008 supported by Fundação para a Ciência e a Tecnologia in Portugal. 2 METHODS Fourier transform (2D): ( ) ( ) ∫ =dxxgxfCCfg ( ) ( ) ( ) dxaxgxfaCCfg ∫−= ( ) ( ) xgxg −= { } { } { } gfgf FFF =∗ 3 RESULTS Fig. 1: Algorithm of the cross-correlation-based method and the associated data pipeline.(FFT is the fast Fourier transform and IFFT is its inverse.) MSE Time [ms] METHODS [N/cm2]2SD Before Registration 23.6 22.8 - CC-based and SSD-based 4.06 2.11 33 PC-based 4.17 2.19 39 REFERENCE RESULTS Principal axes*8.71 10.1 100 Contours-based (with pseudooptimization) (Oliveira et el., 2009) 4.52 2.32 53 Min(MSE)*# (Pataky et el., 2008) 3.98 2.09 9010 By visual inspection, all three algorithms achieved very good registration results, even in the presence of noise. In Figs. 2, 3 and 4, three examples of registration results can be observed. In each figure, from the left to the right, the template image, source image, and overlapped images before and after registration are presented.It should be noticed that, just to aid their analysis, the images presented were post-processed. Some experimental results are indicated in Tables 1 and 2. Each of these results is the average value obtained in the registration of 30 pairs of images. Additionally, the average computational processing time is also presented in Table 2. The values presented in Table 2were obtained using a rigid geometric transformation, as in the reference works considered. The experiments referred in Tables 1 and 2 were done considering an angle resolution of 360º/256 and 1.033 as the base of the logarithm for the conversion of the spectrums to log-polar coordinate system. Table 2: Comparison in terms of registration accuracy, considering the MSE as similarity measure. *Implemented in Matlab.#Using optimization models based on evolutionary algorithms. ( ) ( ) ( ) { }( ) agfdxxagxfa fg ∗=−= ∫ CC ( )( ) ( ) ( ) ( ) ∫ ∫ +∞ ∞− +∞ ∞− +− == dxdyeyxfvufyxf vyuxi π 2 ,, ˆ ,F Shift property: if then ( ) ( ) 0012 ,, yyxxfyxf −−= ( ) ( ) ( ) vufevuf vyuxi , ˆ , ˆ 1 2 200 +− = π Linear scaling property: if then ( ) ( ) byaxfyxf ,, 12 = ( )       =b v a u f ab vuf , ˆ 1 , ˆ12 Rotation property:if then ( ) ( ) θθθθ cossin,sincos,12 yxyxfyxf +−+= ( ) ( ) θθθθ cossin,sincos ˆ , ˆ 12 vuvufvuf +−+= Convolution theorem: where *represents the convolution operation. ( ) ( ) ( ) gfgf FFF =∗ 2.2.2 SSD-based method The assumption behind the SSD-based method is that two images are optimally aligned when their sum of squared differences is minimized. The SSD of two images can be directly computed for all shifts.Let one consider the SSD in function of ashift a: Hence, the first and second terms do not depend on the shift and the last term can be directly computed for all shifts using the convolution theorem. This equation allows to conclude that the optimal shift obtained by minimizing the SSD is equal to the optimal shift obtained by maximizing the CC. ( ) ( ) ( )( ) ( ) ( ) ( ) ( ) ( ) ( ) ∫ ∫ ∫ ∫ −− −−+ =−−= dxaxgxf dxaxgaxgdxxfxf dxaxgxfaSSDfg 2 2 Inverse of Fourier transform (2D): ( ) ( ) ( ) ∫ ∫ +∞ ∞− +∞ ∞− + =dudvevufyxf vyuxi π 2 , ˆ , 2.2 Optimal shift We have developed three methods to get the optimal integer shift:Cross-correlation-based, SSD-based and Phase-correlation-based. 2.2.1 Cross-correlation-based method 2.3 Optimal scaling and rotation Based on the scaling and rotation properties of Fourier transform, the optimal scaling and rotation can be obtained (Reddy and Chatterji, 1996). The necessary first step is to convert the image spectrums to log-polar coordinates systems.Hence, the scaling and rotation became ashift in the log-polar coordinate system.Afterward, the optimal scaling and rotation can be obtained by determining the optimal shift of the log-polar spectrum images using one of the three solutions previously presented to get the optimal shift.Fig. 1 displays a diagram of this approach (Oliveira et al., 2010). 2.2.3 Phase-correlation-based method This last method is based on the shift property of the Fourier transform. Let one consider that gis ashift version of fby vector (x0,y0). Then, computing the cross-power, one has: where represents the complex conjugate of . By computing the inverse of Fourier transform of the cross-power, a Dirac δ-distribution is obtained (Castro and Morandi, 1987). Therefore, based on the coordinates of Dirac pulse, the optimal integer shift is obtained. ( ) ( ) ( ) ( ) ( ) 00 2 * * , ˆ , ˆ , ˆ , ˆvyuxi e vugvuf vugvuf +− = π * ˆ g g ˆ 2.4 Subpixel accuracy To get the subpixel accuracy, one interpolate the CC, SSD or PC matrices in the neighborhood of the optimal integer shift by a parabolic function. Then, this function is used to estimate the optimal real shift. 2.5 Data set and accuracy assessment The methodologies were tested on a data set constituted by 30 pairs of plantar pressure images with dimensions of 45×63 pixels. Two kinds of experiments were done: In the first one, the accuracy was assessed by comparing the geometric transformation estimated by the registration algorithms with a control geometric transformation applied to the images.In the second experiment, the accuracy was assessed by comparing the mean squared error (MSE) obtained after registration with the MSE reported in previous studies for the same data set.In this last experiment, the involved geometric transformations are unknown. Angle Scale Tx Ty METHODS [º] SD SD [pixel] SD [pixel] SD Control values -26.50 1.10 0.00 -2.50 CC and SSDbased -26.51 0.050 1.10 0.002 0.00 0.011 -2.50 0.020 PC-based -26.53 0.027 1.10 0.001 0.00 0.005 -2.49 0.014 AFTER ADDING GAUSSIAN NOISE (mean = 0, σ = 0.5 N/cm2) CC and SSDbased -26.50 0.097 1.10 0.002 0.00 0.020 -2.50 0.024 PC-based -26.53 0.066 1.10 0.001 0.00 0.016 -2.49 0.026 Table 1: Comparison between the control geometric transformation applied to 30 pedobarographic images and the corresponding geometric transformation estimated by the algorithms presented.(SD - Standard deviation, Tx -Translation along x axis, Ty -Translation along y axis) 2.6 Implementation The algorithms were implemented in C++, using Microsoft Visual Studio 8 and tested on a notebook PC with an AMD Turion64 2.0 GHz microprocessor, 1.0 GB of RAM and running Microsoft Windows XP. Fig. 2: Registration of two images of the same foot. Fig. 3: Registration of two images of different feet. Fig. 4: Registration of two images of the same foot after noise adding.