{'fraction_non_alphanumeric': 0.030171980287639545, 'fraction_numerical': 0.007945288142411747, 'mean_word_length': 4.938682062512442, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Registration of plantar pressure images

In this work, five computational methodologies to register plantar pressure images are compared: (1) the first methodology is based on matching the external contours of the feet; (2) the second uses the phase correlation technique; (3) the third addresses the direct maximization of cross‐correlation using the Fourier transform; (4) the fourth minimizes the sum of squared differences using the Fourier transform; and (5) the fifth methodology iteratively optimizes an intensity (dis)similarity measure based on Powell's method. The accuracy and robustness of the five methodologies were assessed by using images from three common plantar pressure acquisition devices: a Footscan system, an EMED system, and a light reflection system. Using the residual error as a measure of accuracy, all methodologies revealed to be very accurate even in the presence of noise. The most accurate was the methodology based on the iterative optimization, when the mean squared error was minimized. It achieved a residual error inferior to 0.01 mm and 0.6 mm for non‐noisy and noisy images, respectively. On the other hand, the methodology based on image contour matching was the fastest, but its accuracy was the lowest. Copyright © 2011 John Wiley & Sons, Ltd.


INTRODUCTION
Plantar pressure distribution provides significant information for clinicians and researchers concerning the structure and function of the foot and the general mechanics of gait. It is, for example, extremely helpful in the diagnosis of foot complaints, development of footwear [1], and for gait analysis. Also, it may be used to compare the loads in the lower limb, either between injured and non-injured, or pre-traumatic and post-traumatic, or pre-operative and post-operative states. It enables comparisons between patients and control groups and provides detailed and specific information on each region of contact [2]. There are a number of different techniques to access the relevant pressure distribution, and, in most cases, the pedobarographic data can be converted into a discrete rectangular array. Therefore, by converting the plantar pressure at each sensor into pixel intensity, techniques of image processing and analysis can be used. Image registration is required by clinicians and researchers for lower limb comparisons, patient follow-up, identification of the main plantar pressure areas and foot classifications [3]. Also, plantar pressure image registration supports pixel-level statistics, which makes the extraction of biomechanically-relevant information more effective than the traditional regional techniques [4]. Several computational methodologies have been developed to carry out image registration, for example principal axes transformation [5], modal matching [6,7], principal axes combined with 590 F. P. M. OLIVEIRA AND J. M. R. S. TAVARES steepest descent gradient search [8], optimization with evolutionary algorithms [9], foot size and progression angle [10], contour matching [11], optimization of the cross-correlation (CC), and sum of squared differences (SSD) both computed in the frequency domain [12], phase correlation [12], and optimization of an image (dis)similarity measure using Powell's method [13].
The main aim of the present work is to compare the five latter methodologies cited earlier, which have revealed several interesting features, such as efficiency and robustness. For comparative purposes, plantar pressure images from three common pedobarographic devices, a Footscan system (RSscan, Olen, Belgium), an EMED system (Novel GmbH, Germany) and a light reflection based system, were used. These five methodologies are introduced in the next section. Afterward, the experimental results are presented and discussed along with the main advantages and disadvantages of each methodology.

Registration based on matching external contours
The registration methodology presented in [11] is based on matching the external contours of the feet represented in the input images. The main steps are as follows: (1) Extraction (i.e., segmentation) of the external contours presented in each input image.
(2) Associate a matching cost to each possible match among the contour points of both feet. This cost among the points is based on the curvature and distance to the corresponding center of pressure of the contours. (3) Search for the best global matching among the points of both contours, which is the one that has the minimum sum of the matching costs. This minimization process uses an assignment algorithm on the basis of dynamic programming and with a constraint to preserve the circular order of the contour points [14]. A matching example using this approach is shown in Figure 1. (4) Computation of the geometric transformation that best aligns the contour points previously matched, on the basis of the minimization of the distance among them by using least squares. (5) Alignment of the input images considering the geometric transformation obtained in the previous step.
The success of this methodology relies mainly on the similarity of the external contours. If the feet shapes are similar, the contours will also be, and consequently a good registration is expected. However, when the feet shapes are significantly dissimilar, the quality of the registration obtained cannot be guaranteed. In addition, this methodology only relies on the contour points, discarding all the Figure 1. External contour matching of two input feet: on the left, and by row, two plantar pressure images to be registered; in the middle, the extracted contours after a sampling process; and on the right, the match found. 591 information conveyed inside the contours; hence, variations in the pressure distribution occurring in these regions do not affect the registration obtained. The size of the images has a low influence on the processing time for this methodology, as it is mainly affected by the number of points of each contour, which are obviously less than the total number of pixels of the input images.

Registration using the phase correlation technique
The phase correlation (PC) technique is based on the shift property of the Fourier transform; that is, a shift of a function in the space domain is represented by a shift in its phase when the function is represented in the frequency domain.
If two input images to be registered are represented by the real functions f and g, and their Fourier transforms, F and G, respectively, and if g .x, y/ D f .x x 0 , y y 0 /, then, according to the shift property [15] and by computing the cross-power where G represents the complex conjugate of G.
By computing the inverse of the Fourier transform of the cross-power, a Dirac ı-distribution centered at .x 0 , y 0 / is attained. Therefore, the coordinates of the Dirac pulse indicate the optimal integer shift. The Fourier transform and the inverse Fourier transform can be efficiently computed using the fast Fourier transform and the inverse fast Fourier transform, respectively.
To increase the shift accuracy to a subpixel level, the neighborhood of the strongest peak is interpolated using a quadratic function, and then, the continuous coordinates associated with the maximum value are obtained. In this registration methodology, the interpolator is built using the 3 3 neighborhood, which is centered on the highest peak, and then least squares are used to compute the coefficients of the quadratic function.
The aforementioned procedure allows the determination of the shift, but requires an initial correction of possible scaling and rotation between the images. With the properties of the Fourier transform (see Appendix A), the scaling and rotation between the input plantar pressure images can be estimated from their spectrums. First, both spectrums are converted into the log-polar coordinate system. Then, the shift between both log-polar spectrums is determined using the phase correlation technique. And with this shift, the rotation and scaling of the spectrums can be estimated [12], and consequently, the rotation and scaling between the input images ( Figure 2).
Because of image spectrum symmetry, any two images that differ by a rotation of only rad, will have the same spectrum. Therefore, if there is a rotation angleˇ, it needs to be tested if it is really equal toˇor toˇC instead. In this registration methodology, this test is performed when the optimal shift between the input images is determined, considering both images after an angle correction ofˇandˇC . Thus, the rotation angle is defined on the basis of the image (rotatedǒ rˇC / that gives the highest peak in the search for the optimal shift. When compared with the external contour matching methodology, this methodology has the advantage of being global, that is, all the image pixels are used to estimate the optimal geometric transformation. On the other hand, if the input images are corrupted by noise, for instance, background noise, the quality of registration can be jeopardized. The computational processing time required by this methodology depends mainly on the size of the images to be registered. Therefore, this methodology is computationally more expensive than the previous one ("Section 2.1") because it requires the computation of, at least, three 2D-fast Fourier transform and two 2D-inverse fast Fourier transform [12]. Besides, before converting the two input images into the frequency domain to correct their possible non-squareness and to avoid wrap-around effects in the frequency domain, the images need to be pre-pad with zeros.

Registration based on direct maximization of cross correlation
This registration methodology is based on the assumption that two images are best registered when the computed CC is maximized. So, if two input images to be registered are represented by the real functions f and g, and their CC in function of a shift a is as follows: then, by the convolution definition, where N g .x/ D g . x/ and represents the convolution. And from the convolution theorem, where F represents the Fourier transform, and k is a constant that depends on the specific Fourier transform normalization. Therefore, computing the inverse of the Fourier transform of the product in Equation (5), the correlation for all shifts can be obtained. Then, the coordinates of the point that have the highest values represent the desired optimal integer shift.

593
The shift accuracy of this registration methodology is increased by using the same procedure as in "Section 2.2". The optimal scaling and rotation is also estimated using the approach in "Section 2.2". This methodology has the same disadvantages as the phase correlation technique methodology ("Section 2.2"); however, it presents one advantage: the value of the CC or the normalization of that value is frequently used to measure the similarity between the input images.

Registration based on direct minimization of the sum of squared differences
The SSD is computed using the intensity of the pixels of both images to be registered. Hence, the lower the SSD is, the better the registration is.
If two input images to be registered are represented by the real functions f and g, and their SSD in function of a shift a: The last two terms in Equation (6) are constants and can be easily computed on a pointwise multiplication basis. The remaining term, the first one, can be transformed into a convolution and efficiently computed using the Fourier transform, adopting the same procedure as was used in the evaluation of the CC ("Section 2.3"). Then, the coordinates of the point that has the lowest value represents the desired optimal integer shift. In this registration methodology, the shift accuracy is increased by using a similar technique as the one used in "Section 2.2", but with the difference that here, the search is for the pixel with the lowest intensity. To obtain the scaling and rotation that minimize the SSD, a procedure similar to the one used in the phase correlation was employed.
This registration methodology has exactly the same advantages and disadvantages as the methodology based on the direct optimization of the CC ("Section 2.3"). Besides, comparing Equations (3) and (6), it is expected that the optimal geometric transformations obtained by this methodology and obtained by the CC-based methodology are identical.

Registration based on the iterative optimization of an image intensity (dis)similarity measure
This family of registration methodologies is based on the optimization of an image (dis)similarity measure, usually related to the intensities of the image pixels. Thus, the geometric transformation that optimizes the (dis)similarity measure adopted is used to register the input images. There are several multidimensional optimization algorithms that can be used to optimize the (dis)similarity measure adopted, and a great variety of (dis)similarity measures. The convergence of the optimization algorithms depends highly on its optimization strategy, and also on the smoothness and capture range of the (dis)similarity measure used. However, in most cases, the optimization algorithms only achieve the parameters of the geometric transformation that successfully register the input images if these images are not significantly misaligned. Hence, to overcome such a limitation, a pre-alignment is usually performed before the optimization process.
In the present work, the optimization solution presented in [13] was used. It is a solution based on a two-step approach: in the first step, a pre-registration is obtained; afterward, in the second step, the algorithm searches iteratively for the geometric transformation that optimizes the adopted (dis)similarity measure ( Figure 3). The multidimensional optimization algorithm used in the second step is based on Powell's method [16]. The geometric transformations allowed in the solution implemented are as follows: rigid (shift in x-axis and y-axis, and rotation), similarity (shift in x-axis and y-axis, rotation, and linear global scaling), affine, projective, and polynomial up to fourth degree.
Three different image (dis)similarity measures were experimentally used: the MSE, which is a normalization of the SSD, the mutual information (MI) [17,18]; and a dissimilarity measure based on the exclusive-or (XOR) between the input images after binarization [9].  Let I 0 and I 1 be two input images to be registered, with N M pixels, and bin .I 0 / and bin .I 1 / the binarized images of images I 0 and I 1 , respectively.
The MSE is given as follows: Consequently, the lower the MSE is, the better the input images are registered.
The XOR between the input images is computed as where j . j is the cardinal function and˚is the XOR operator. In the binarization process, the value 0 is attributed to all image pixels that have a pressure intensity inferior to the minimum threshold of the acquisition system used and 1 to all the remainder pixels. Therefore, this measure provides a measure of non-overlapped pixels; thus, the lower the XOR value is, the better the registration is. This image dissimilarity measure is only adequate for the registration of shapes represented in images without background noise. The MI is defined by where H.I k / is the Shannon's entropy of the pixels in image I k and is the joint entropy. For image registration purposes, higher MI values imply higher registration quality. Here, the MI was computed as in [19] and using 32 bins in all experiments carried out.
In comparison with the other four registration methodologies presented earlier, the main advantages of this methodology are the possibility to obtain a superior registration result by using the most suitable image (dis)similarity for the case under evaluation, and the high accuracy that is attained when the optimization algorithm converges to the optimal value.
Unlike the remaining four methodologies, the computational processing time demanded by this methodology neither depend on the foot size (as in the contour method) nor on the size of the input images. In fact, the required processing time depends greatly on the ability of the optimization algorithm to find the convergence path. Besides, there are two deterministic factors that influence the processing time: the computation of the image (dis)similarity measure used, and the image interpolation approach employed in the image resampling process. Other important factors on the required processing time are the quality of the initial registration because the lower the quality is, the higher is the time needed for the optimization process and the smoothness of the (dis)similarity measure used.

Dataset
In the experimental evaluation, a dataset of 36 plantar pressure images acquired from three common pedobarographic systems (12 images per system) was randomly built. Each of the three subsets used contained normal, low arched, and high arched feet. All the data were acquired with the subjects walking along a straight path at their normal speed.
The first subset contained peak pressure images (45 63 pixels) collected using a 0.5 m Footscan system, with a pixel resolution of 5.08 7.62 mm and a pressure sensitivity of 0.7 N/cm 2 . The second subset contained 12 peak pressure images (32 55 pixels) acquired using an EMED system, with a resolution of two sensors per cm 2 (approximately equivalent to 7.07 7.07 mm). The pressure sensitivity of this system is 5 kPa with a minimum threshold of 10 kPa. The third subset contained 12 peak pressure images (160 288 pixels) acquired using a pedobarographic light reflection system [20,21] with a resolution of approximately 1.8 1.8 mm. The color images acquired by the image camera of this last system were converted to gray scale images. The calibration of this device, that is, the relation image pixel intensity or pressure applied, was not addressed because it is outside the requirements of this work.
The images acquired using the Footscan system were vertically stretched by a factor of 1.5 to correct for non-square sensor array spacing. Image transformations were performed (here and throughout) using bilinear interpolation resampling [22].
Regarding the light reflection system, before the registration process, the images acquired were intensity rescaled to reduce the background level and noise ( Figure 4). The rescaling was made on the basis of the histogram of the image pixel intensity: the higher peaks of the histogram represent the larger regions of the input images, which are the background. The image pixel intensity was classified into 32 bins with the same width. Then, all image pixels with intensities lower than the intensity associated to the bin that follows the bin with the maximal intensity, were set to 0 (zero). Afterwards, the remainder image pixels were linearly rescaled between 0 and 255.

Accuracy assessment
Two control (i.e., known) transformations, a rigid geometric transformation (involving a shift and a rotation) and a similarity geometric transformation (composed of a global scaling, a shift, and a rotation), were used with the dataset of 36 real plantar pressure images to assess the registration accuracy: To simulate a real intra-subject misalignment, the following rigid geometric transformation was applied: a shift of -40 mm, 20 mm, and a rotation of -25 o . The inter-subject registration was simulated by applying a similarity geometric transformation to each image of the experimental dataset. The similarity transformation applied, which is in line with the usual values found in this domain, was made up of a scaling factor of 1/1.3, a rotation of 12 o and a shift of 15 mm, -50 mm. Then, the transformed images were registered with the original ones. The registration accuracy was measured by comparing the parameters of the known transformation applied and the ones estimated by each of the registration methodologies under evaluation. The residual error (RE), that is, the square root of the mean squared difference between the exact position expected for each pixel and the position obtained using the registration methodologies, was used as measure of the registration accuracy.
The robustness of the methodologies under comparison against spatial localization noise was also studied by adding Gaussian noise to the images that were misaligned by the control geometric transformations and the registration errors were analyzed.
Additionally, the differences between the RE values obtained by the registration methodologies under comparison were assessed using one-way ANOVA and Dunnett's T3 post-hoc comparisons. The statistical analyses used SPSS 16.0 (SPSS Inc., Chicago, IL, USA).   Table I presents the mean RE values and the mean processing times obtained in the registration experiments when the rigid control geometric transformation was used. The significance (p value) of the differences between the methodologies under comparison are also included in Table I. This table was organized in terms of the registration methodology and the pedobarographic systems. Table II presents similar experiments as in Table I, but after adding Gaussian noise (mean = 0 mm; SD = 2.5 mm in each axis) to the misaligned images. The RE values obtained using the similarity control geometric transformation are shown in Table III.  Table IV includes the p values from Dunnett's T3 post-hoc comparisons. These post-hoc comparisons compare group means on each of the nine subsets of image pairs defined in "Section 2.7", built over all combinations of the seven registration solutions. In Table IV, the lower triangular parts of each 7 7 matrices correspond to the measures indicated by the labels in the leftmost column, 597 and the upper triangular parts correspond to the measures indicated by the labels in the rightmost column. The meaning of the group differences can be inferred from the values shown in Tables I-III. For example, the minXOR-PC comparison of the subset of images pairs from the EMED system after a rigid control transformation yielded a p value of 0.023; from Table I one can realize that the minXOR (iterative minimization of the XOR) had a lower mean RE value than the PC (phase correlation technique based methodology).
It should be noted that, in Tables I-IV, the results of the fifth methodology, which iteratively optimizes an intensity (dis)similarity measure (MSE, MI or XOR), were obtained using the external contour matching methodology in the pre-registration step because this approach had revealed its suitability in our previous studies [13]. Also, the minimization of the XOR similarly measure was not addressed using the images from the light reflection system because these images have considerable background noise thwarting successful registration results.
In Figure 5, three registration results are presented, considering the registration after the misalignment obtained by applying the control geometric transformation, with and without the addition of Gaussian noise (mean = 0 mm, SD = 2.5 mm in each axis) to the misaligned images.
The methodologies based on the direct optimization of the CC and SSD, and the phase correlation methodology achieved good and identical results. The mean residual errors obtained using these three methodologies were always inferior to 1.4 mm for all experiments, even for the noise-corrupted images. Among these three algorithms, the CC and the SSD-based methodologies achieved the same results, and in most cases, better results than the methodology based on phase correlation.

DISCUSSION
In all registration experiments, the residual errors obtained were always smaller than the resolution of the pedobarographic devices used, which guarantees that all methodologies under comparison are suitable for clinical and research use.
The most accurate methodology was the one based on the iterative optimization because it achieved the lowest residual errors. This result was already expected because in the tests performed, the optimization process started after a pre-registration very close to the optimal one. We could register the input images by applying the iterative optimization methodology directly without a preregistration step. However, because that methodology was developed to be robust only against small misalignments, the convergence to the optimal solution was not guaranteed.
There was no statistical significant difference (p = 0.16) between the mean RE values obtained from the registration of the images acquired by the three pedobarographic systems. However, the RE from the light reflection system tended to be lower than the RE from the Footscan and EMED 599   systems (Table I). This finding is explained by the superior spatial resolution of the light reflection system. No significant differences were observed for the RE from the EMED system and the Footscan system,.
The methodology based on the iterative minimization of the MSE leads to an RE always lower than 0.6 mm, which is considerably smaller than the resolution of the pedobarographic devices used. Among all the registration experiments carried out, this methodology attained the lowest mean residual error (with statistical significance in most cases). This finding indicates that the minimization of the MSE is generally a better registration option than the maximization of the MI or minimization of the XOR.
The methodology based on matching of the external contour was the fastest. However, its overall accuracy was not as good as the other methodologies. The mean residual errors obtained using this methodology were greater than those of the other methodologies (with statistical significance in most cases).
In the comparison made, only rigid and similarity geometric transformations were used because the methodologies based on the contour, matching the phase correlation and direct optimization of the CC and SSD, are only suitable for these kinds of geometric transformations. Figure 5. Three examples of registration results. In each row, from the left to the right: image used as template, misaligned image, overlapped images before the registration, overlapped images after the registration, and difference between the registered images. On the first row, the original image was acquired by the light reflection system; on the second row, the original image was acquired by the EMED system; and on the third row, the original image was acquired by the Footscan system. On the first two rows, the intra-subject registration was simulated, and the template images were artificially distorted (i.e., corrupted), by adding Gaussian noise (mean = 0; SD = 2.5 mm in each axis). On the third row, the inter-subject registration was simulated. (For visualization enhancement, the images were colored and the pixels of the images from the light reflection system with intensity lower than 20 were set to zero to hide some remaining background noise.) As described in "Section 2.1", the accuracy of the registration methodology based on the matching of the external contours is influenced by the shape of the feet to be registered and the quality of the contours extracted from the images. So, for the registration of abnormal feet, the quality is expected to decrease. The remaining methodologies can be more efficient for the intra-subject registration of abnormal feet because the foot shape has a much less influence on the registration process. Moreover, as the plantar pressure images to be registered are from the same foot, the pressure distribution should be similar in both images and the registration is facilitated.
In real inter-subject registration (i.e., the registration of plantar pressure images from different feet), the accuracy of all the methodologies compared may be reduced because the successful overlapping of all the foot regions is difficult when just similarity geometric transformations (composed by linear scalings, shifts, and rotations) are used. This problem can be overcome by considering curved geometric transformations; however, the resultant deformation of the foot shape can make its use in further analysis impossible, for example, for footprint index calculations.
Even using a not up-to-date computer, the processing time required by the five computational methodologies to register the images acquired by the Footscan and EMED systems were always very low (much less than a second, Tables I & III). Regarding the images acquired by the light 602 F. P. M. OLIVEIRA AND J. M. R. S. TAVARES reflection system, obviously, the processing time required was higher than the ones required to register the images acquired by the other systems, but remained at an acceptable level for clinical and research purposes.
As a final conclusion, one can state that all the registration methodologies tested revealed high accuracy, speed, and robustness against image noise and arbitrary shifts, rotations and moderated linear scalings.

APPENDIX A
Consider two functions f and g from R 2 to R and their Fourier transforms F and G, respectively.
Scaling property: If g .x, y/ D f .ax, by/, then G .u, v/ D 1 jabj F u a , v b . Log-polar transformation property: Suppose that a rotation of amplitudeˇfollowed by a scaling of amplitude s .s > 0/ was applied to the real plane. For simplicity, consider that the rotation and scaling were applied around the origin point. Then, a point with rectangular coordinates .x, y/ is transformed into a point with rectangular coordinates .s x cosˇ s y sinˇ, s x sinˇC s y cosˇ/. Then, if the point .x, y/ has log-polar coordinates .log r, Â/, then the point .s x cosˇ s y sinˇ, s x sinˇC s y cosˇ/ has log-polar coordinates .log .s r/ , Â Cˇ/ D .log s C log r, Â Cˇ/.
Thus, a scaling and rotation in a rectangular coordinate system correspond to a shift in a log-polar coordinate system.