Morphological Foot Model for Temperature Pattern Analysis Proposed for Diabetic Foot Disorders

Infrared thermography is a non-invasive and accessible tool that maps the surface temperature of a body. This technology is particularly useful for diabetic foot disorders, since it facilitates the identification of higher risk patients by frequent monitoring and therefore limits the incidence of disabling conditions. The aim of this work is to provide a methodology to explore the entire plantar aspects of both feet, based on infrared thermography, for the assessment of diabetic foot anomalies. A non-invasive methodology was established to identify areas of higher risk and track their progress via longitudinal monitoring. A standard morphological model was extracted from a group of healthy subjects, nine females and 13 males, by spatial image registration. This healthy foot model can be taken as a template for the assessment of temperature asymmetry, even in cases in which partial amputations or deformations are present. A pixel-wise comparison of the temperature patterns was carried out by Wilcoxon´s matched-pairs test using the corresponding template. For all the subjects, the left foot was compared to the contralateral foot, the right one, providing a map of statistically significant areas of variation, within the template, among the healthy subjects at different time points. In the female case, the main areas of variability were the boundaries of the feet, whereas for the male, in addition to this, substantial changes that exhibited a clear pattern were observed. A fast and simple monitoring tool is provided to be used for personalized medical diagnosis in patients affected by diabetic foot disorders.


Introduction
Abnormal plantar temperature in patients with diabetes may be an early sign indicating the appearance of foot disorders which include peripheral arterial disease, neuropathy, and infection among others [1,2]. These pathologies are often resistant to healing and require long-term medical care [3]. Preventive care by regular monitoring and correct assessment limits the incidence of disabling conditions and the lower-limb amputations associated with acute cases [4][5][6]. Therefore, improved medical care, by monitoring the patient's condition over time, is required to identify higher risk patients and prevent the significant impact caused by economic burden [7]. This demands fast, simple, autonomous, repeatable and precise screening protocols, avoiding error-prone and labor-intensive processes as well as the requirement of special training [6,7].
Infrared thermography is a non-invasive tool that maps the surface temperature of a body [8]. This marker-free, safe, accessible, contactless, portable, and easily repeatable technique has been widely applied in vascular diseases, particularly diabetic foot template. Consequently, the comparison of the contralateral foot can be readily performed on the same subject as well as the average healthy model extracted.

Image Acquisition
IR images were acquired with the affordable thermal camera model TE-Q1 Plus from Thermal Expert ™ (i3system Inc., Daejeon, Korea) which has been previously described and was previously calibrated [10]. RGB-D images were acquired with an Intel ® RealSense ™ D415 camera (Intel Corporation, Santa Clara, CA, USA). The resolutions of the RGB-D images were conditioned by the maximum value achievable by the depth sensor (1280 × 720 pixels), whereas for the IR images they are determined by the TE-Q1 Plus sensor (384 × 288 pixels). IR images were stored in Tag Image File Format (TIFF) and the RGB images in Portable Network Graphic (PNG) format. In both cases, 32 depth bits were used. These cameras were assembled together in a customized support, manufactured in a 3D printer, that kept the cameras horizontally aligned. RGB-D images were translated and cropped to match the coordinate system and size of the IR images [25,26]. Thus, each pixel was represented by five channels (R,G,B,D,IR).
An image acquisition campaign was carried out in May 2021 among coworkers at our facilities. The group of volunteers consisted of 22 healthy subjects, nine women and thirteen men, whose average ages were 36.4 ± 10.5 and 36.1 ± 7.4 years old, respectively. The average female European foot size was 37.5 ± 1.3, ranging between 36 and 40, whereas the male's was 42.8 ± 1.0 in a range between 41 and 45. No subject was part of the acquisition more than once and none of these volunteers presented partial amputations or deformations of the feet. Informed consent was obtained from all subjects involved in the image acquisition. In addition, acquired data were codified and anonymized to ensure subject confidentiality and data protection.
A non-constrained acquisition protocol was selected for the presented setup in which the working distance between the cameras and the subject's feet was set at approximately 80 cm. Furthermore, the acquisition was carried out in a room with controlled luminosity, humidity, and an average ambient temperature of 25 • C. Prior to image acquisition, the cameras were stabilized for approximately 15 min [10,31]. Four sets of images were acquired per subject. The first one, at T0, was captured as soon as the person took their shoes off and sat with their legs extended forward or laid down in a supine position with their feet off the ground. Thus, this takes the minimum amount of time since the person took their shoes off and image acquisition occurred at T0. The subsequent images were taken every five minutes up to 15 (T5, T10 and T15); meanwhile, the subject was at the same resting position keeping their feet off the ground. Therefore, a total of 88 images were available. It must be clarified that, prior to image acquisition, the subjects were not exercising or performing strenuous physical activity, bathing, or using hot or cold compresses for cleaning purposes.

Ground Truth Establishment
RGB images were labeled manually by three independent and unbiased researchers to extract the sole of the feet from the background. Each foot was segmented separately using the manual and the semi-automatic tools provided by the software application ITKSnap (version 3.8.0) [32].
Therefore, six masks were generated from each image acquired, corresponding to right (R) and left (L) feet delineated by each researcher, which facilitated the quantification of the inter-and intraresearcher variability. The masks were used as inputs for the Simultaneous Truth and Performance Level Estimation (STAPLE) algorithm [33]. The output masks, representing a probabilistic estimation of the true segmentation, were established as the ground truth segmentation.

Foot Model Extraction
A foot model was extracted for each gender and foot separately, employing all the images previously acquired. For simplicity, the procedure to extract the right foot model is described but an analogous procedure was employed for the left foot. The method consisted in selecting one of the images as a reference or fixed image, and subsequently registering, by an affine transformation, the remaining ones to that reference. A simplified scheme describing the workflow employed is presented in Figure 1. As mentioned in the previous section, the RGB images for each volunteer were segmented by three independent researchers. Then, the ground truth segmentation for each subject was obtained by applying the STAPLE algorithm (STAPLE #1), which will be employed as reference for a specific subject. For instance, as shown in Figure 1, the subject on the left was considered the reference to extract the final model. The other ground truth segmentations, corresponding to the remaining subjects, were spatially registered to the reference using the affine transformation from the SimpleElastix module (ELASTIX) [34] in Python 3 [35].
A registration procedure based on IR images is complicated since these images lack the spatial information required to apply the usual feature-based and intensity-based registration techniques. However, by employing a third modality, in this case the binary image corresponding to each mask (ground truth segmentation), registration can be successfully achieved, as previously reported [26,36]. This approach requires an image registration between modalities (RGB-D-IR) that is as precise as possible. Therefore, prior to the procedure to extract the foot model, any possible misalignment between the RGB-D images, or the associated ground truth segmentations, and the IR images were corrected by performing a 2D affine transformation. The registrations between masks were also based on an affine transformation [26]. When all the masks were in the same coordinate system, the STAPLE algorithm (STAPLE #2) was applied again in order to extract the right feet model associated with the subject used as reference. Since this approach may be dependant on the image selected as reference, the process was repeated as many times as there were available images. In this manner, all subjects were used as reference, and the one providing the highest overlap measures was selected as the final foot model.
Once the final right foot model was established as global reference, based on the overlap measures, the ground truth segmentations associated with each subject were subsequently registered to this model by another affine transformation. All these spatial transformations were also carried out using the SimpleElastix module [34] in Python 3 [35]. Since, as previously ensured, RGB-D-IR images were in the same coordinate system, the same affine transformation used for each segmentation mask was applied to the corresponding RGB-D and IR images using the Transformix module [34,37] in Python 3 [35]. That is, the file of parameters corresponding to each registration was employed in the Transformix module to bring the IR image to the reference model. Finally, the foot model was applied to the transformed IR image. A simplified scheme of the workflow is presented in Figure 2.

Evaluation Metrics
The STAPLE algorithm provides a probabilistic estimate of the true segmentation and a measure of the performance level represented by quantitative values for sensitivity and specificity [33].

Statistical Analysis
Statistical analysis was performed in RStudio v1.1.453 (RStudio, Inc., Boston, MA, USA) [41]. Wilcoxon's matched-pairs test [42,43] was used to determine intergroup differences at each time point as well as to assess longitudinal changes. When significance was detected, the probability value was indicated (p). The Hochberg and the false discovery rate adjustment methods [44,45] were employed as a correction for multiple testing. A significance level of 5% was considered in all tests.

Ground Truth Establishment
The overlap measures for the ground truth segmentations, which give an account of the inter-and intraresearcher variability, were averaged and are listed in Table 1. According to Wilcoxon's matched-pairs test, no intergroup differences, that is right versus left foot overlap measures, were found at any time point. An exception was found for specificity and FN at T5 (p < 0.05).
The fact that the workflow is based on the segmentation of the RGB images, instead of the IR ones, presents a major advantage. Many subjects present cold toes, in which tissue limits are difficult to define in the IR images. Therefore, the segmentation step is complicated, strongly dependent on the observations, and prone to errors.

Foot Model Extraction
The overlap measures corresponding to the foot models extracted, considering each subject as reference, were averaged and listed in Table 2. As explained in the previous section, although the model extraction was performed separately for each foot, the overlap measurements are displayed together for comparison purposes. The coefficients of variation (CV) for most of the overlap measures, estimated as the ratio between the standard deviation and the mean, were below 0.4% and 0.25% for the female and male models, respectively. An exception was observed for the incorrectly detected conditions, FP, and the incorrectly rejected conditions, FN, in which the CVs were up to 3% and 5% for the female and male models, respectively.
According to Wilcoxon's matched-pairs test, for the left foot in the female model, significant longitudinal differences (p < 0.05) were found for specificity and FN between T0 and T10. For the male model, significant longitudinal differences (p < 0.05) were found among the overlap measures for the left foot between the baseline, T0, and the subsequent time points, T5, T10, T15. However, for specificity and FN, significant changes were not detected. Regarding the right foot, no longitudinal changes were observed for the associated overlap measures for any of the models. For the female model, intergroup differences were found at all time points for specificity and FN, which were significantly higher for the right foot as compared to the left one (p < 0.01). Similarly, sensitivity and FP were significantly higher for the right foot but only at T10 and T15 (p < 0.01). For the male model, intergroup differences were found at all the time points for all the overlap measures (p < 0.05), which, as observed for the female model, presented higher values for the right foot. An exception was found at T15 in which specificity and FN presented comparable values for both feet.

Foot Model for All the Available Shoe Sizes
The model providing the highest overlap measures for each foot, usually DICE and IoU, was selected as the final model. Female and male models have been extracted with the aim to consider the range of sizes characteristic of each gender. The overlap measures for these models are listed in Table 3. The masks corresponding to the final models for the female and male groups of subjects are displayed in Figure 3a,b, respectively.  The discrepancies are not large between the different extracted models. For the left foot, according to Wilcoxon's test, no significant differences between the female and male models were found among the overlap measures. However, for the right foot, sensitivity and FN were significantly increased, whereas specificity and FP were decreased for the male model as compared to the female one (p < 0.01).

Foot Models for Variable Range of Shoe Sizes
The robustness of the final model was tested by creating additional models in which the available shoe sizes were grouped in alternative ranges. First, the most common shoe size was considered, being 37 and 43 for females and males, respectively. Then, a close range was established selecting one size above and below, that is, the female ranged from 36 to 38, whereas the male was between 42 and 44. The overlap measures are displayed in Table 4.
As can be seen, the performance of these models is quite similar to the final model extracted previously that employed all the available shoe sizes. In fact, according to Wilcoxon's test, no significant differences were detected among the overlap measures corresponding to these female models for any foot. Similarly, for the male models, no significant differences were found in the overlap measures for the left foot. However, regarding the right foot, DICE and IoU were significantly higher for the male model, extracted for a shoe size of 43, as compared to the close range model, extracted using show sizes between 42 and 44 (p < 0.05). Additionally, sensitivity was significantly higher for the model extracted with all shoe sizes in comparison to the one extracted using a close range (p < 0.01). Similarly, the model extracted for a single shoe size presented a higher sensitivity compared to the model extracted using a close range (p < 0.05). FP presented the opposite trend to the one observed for the sensitivity, in which values were higher for the close range model as compared to the single shoe size model (p < 0.05) and the model extracted using all shoe sizes (p < 0.01).

Partial Foot Amputations or Deformations
One of the most common limitations in the assessment of diabetic foot disorders is that the contralateral foot analysis is prevented by the existence of partial foot amputations or deformations. In these cases, there is not a straightforward method to perform a comparison if the same deformation is not present on both feet or simply because the severity of the amputation does not allow it [15]. Furthermore, despite the fact that these are common features exhibited by persons affected by diabetic foot disorders, these patients are usually excluded from the assessment studies [15,46].
The workflow presented here, in which a standard foot model is employed for image analysis, facilitates the assessment, to a certain extent, in case of partial amputations or deformations. The contralateral analysis of the missing areas cannot be performed. However, the remaining areas can be compared to the opposite foot of the same patient, via spatial registration to the standard model extracted for each foot, or directly to the model associated with a healthy subject.
In order to demonstrate the performance of the workflow described in Figure 2, the left foot segmentation of a random male subject was modified to simulate various degrees of partial foot amputations. Three configurations were implemented in which the extent of the amputation considered solely the big toe, several toes and up to half foot. Figure 4 illustrates the successful completion of the main steps within the workflow considering these partial foot amputations. Table 5 lists the overlap measures after registering the segmentation mask to the left foot model as displayed in Figure 4b,e,h. Table 5. Overlap measures between the model associated with the left foot and the corresponding segmentation simulating an amputation as shown in Figure 4. Values are expressed as percentages.   Regarding deformations, an analogous procedure to the one described above for the amputations was carried out. In this case, instead of simulations, a subject affected by a different degree of deformation in both feet was employed to demonstrate the performance of the workflow described in Figure 2. Figure 5 illustrates the output after the main steps in the workflow. Table 6 lists the overlap measures between the extracted model for each foot and the corresponding segmentation. It must be mentioned that this subject was not included in the workflow to extract the standard model (Figure 1), in which male subjects with European shoe sizes between 41 and 45 were considered. The acquisition protocol for this subject consisted of only two images, at T0 and T5, and for this reason, it was solely employed for demonstration purposes. In addition, the shoe size for this subject was 46, and therefore is out of the range considered to extract the standard model. That is, the performance of the method to accommodate to this particular size would have been superior if an extended size range would have been employed.  Figure 2 for the left foot whereas the lower one corresponds to the right foot. The first column, (a,d), displays the original IR images and the corresponding segmentation (pink mask). The second column, (b,e), shows the segmentation registered to and overlapped with the left and right models, respectively. The third column, (c,f), illustrates the final result in which the segmentation is applied to the original IR image for further processing. As observed, the standard foot model extracted provides a tool to handle partial foot amputations and deformations. However, the performance level has a strong dependence on the extension of these features and therefore is directly related to the degrees of freedom required to achieve the spacial registration between the model and the corresponding segmentation mask.

Statistical Analysis
Once the female and male models were established for each foot, the IR images were transformed following the workflow described in Figure 2. Figure 6a,c display the registered IR images for all the corresponding female and male subjects, for each foot, respectively. As can be seen, all feet are in the same coordinate system, although some mismatches can be observed near the toes as well as on the arch area. The variability observed may be due to small mismatches associated with the spatial registration as well as to the inherent heterogeneity of the group of subjects, since people vary in shape, gender, age, etc.
Subsequently, for further processing, the IR images corresponding to the right foot were flipped and registered to the left foot model. Figure 6b,d show all the overlapped images for the female and male subjects, respectively. A pixel-wise comparison was carried out to identify the foot areas more prone to variability within a group of healthy subjects. Since the feet of all the subjects were in the same spatial coordinates, all the left feet were compared to all the right feet, at each time point, using Wilcoxon's matched-pairs test. A correction for multiple comparisons was conducted using the false discovery rate. It must be noted that an analogous analysis was performed for each gender. Figure 7 shows the pixels in which changes are statistically significant (p < 0.05) at each time point and for each gender.
As can be observed, for the female subjects the largest pixels' variability is found at the boundaries. This may be caused by small mismatches in the several spatial registrations performed, previous to this step, within the workflow. In the male model, the pixels' variability is more heterogeneous and not so subtle. A part of the same variation observed for the female model within the boundaries, that could be reduced by the application of a morphological dilation, the extension of the areas prone to change is considerably larger.
A source of controversy in the assessment of diabetic foot disorders is the time required to achieve temperature normalization in a clinical setting prior to thermography acquisition. In some studies, normalization is expected after 10 min [15], whereas in others it is at 15 min [8]. As demonstrated, the difference between the plantar aspects of both feet at T0 does not suffice to pinpoint the areas in which an anomaly is found. By looking at Figure 7, for the female case, no substantial differences are subsequently found. However, for the male case, it can be seen that, just after 5 min, at T5, an entire new perspective can be appreciated. Furthermore, this trend is roughly maintained to T15. Therefore, at least a few minutes are required to achieve temperature normalization. This is an important issue to be considered by medical practitioners in a clinical setting but also when aiming at self-monitoring programs to be carried out at home by the patients [47]. Furthermore, in some studies, the mean temperature of the complete plantar surface is used for assessment [15] or even selected spots [14,47]. In addition to the contralateral pixel-wise comparison, the average temperature for each foot was calculated and the difference, left minus right foot, was recorded for each individual. Afterwards, the mean difference across subjects was calculated at each time point and displayed in Table 7. As can be noticed, the mean temperature difference is larger for the male group as compared to the female, although not statistical significance was detected according to Wilcoxon's test. Furthermore, the largest temperature difference for the female group of subjects was around 0.8 • C at T15, whereas for the male group, the largest difference was 1.6 • C at T0. In order to facilitate the visualization of the changes observed, the individual temperature differences at each time point are plotted in Figure 8. For some individuals in both groups, negative values were detected, this indicates that the right foot was hotter than the left one. Another source of controversy is the establishment of the threshold at which an acute situation is detected. As mentioned above, 2.2 • C is usually consider as the standard threshold at the clinical settings for DFU [7,17,18]. Furthermore, this threshold seems to be independent of the ambient temperature, the time of the day, and the walking activity of the subject [47]. However, in patients affected by diabetic foot infection, a plantar foot asymmetry of 1.35 • C suffices to seek urgent treatment [16]. As demonstrated here, some healthy subjects presented a mean temperature difference above this threshold. Therefore, the heterogeneity among subjects must be considered, and the establishment of a personalized threshold and treatment are mandatory. In fact, it is plausible that using the mean values solely to assess the foot asymmetry in a single measurement may not provide enough information. As previously reported, the above-threshold asymmetry must be confirmed longitudinally to improve its validity [47].
Once the model has been extracted, the methodology presented here would allow us to monitor the same patient longitudinally as long as required. Furthermore, the images acquired in time can be easily compared pixel-wise, facilitating the localization of areas, within the entire plantar aspects of both feet, which may present an anomaly. In this particular case, the comparison between left and right feet for the female group, was performed in 129.7 ± 14.1 s for 36 input images using a standard processor. For the male group, which consisted of 52 images, the analysis was performed in 127.4 ± 4.4 s. Each individual analysis, in which eight input images are used (T0, T5, T10, T15), was 125.5 ± 1.4 s. The spatial registration to the corresponding model depends on the number of images to be registered as well as the processor used. In the present case, using the same processor, an execution time of less than 2 s was estimated for the registration of the four images associated with a single subject. Thus, the methodology presented is executed fast and without substantial technical requirements. Furthermore, it is feasible to include it in the standard clinical settings as well as in a mobile application to be used by the patient at home.
The acquisition protocol associated with thermography, which is fast, simple and easily repeatable, is the major advantage as compared to a conventional thermometer. The temperature of the entire plantar aspects of both feet can be acquired in a single image by thermography. However, for practical reasons, when employing a thermometer, only a few specific locations are measured symmetrically in both feet. This is an extremely timeconsuming process. In cases in which the patient is required to measure the temperature several times a day [47], the assessment may be tedious for the patient as well as for the clinical practitioners. Furthermore, when considering the self-assessment performed by the patient at home, the procedure may be discouraging.
A major limitation of the present study is the low number of subjects included in addition to the absence of patients with diabetes. Since, across people and over time, there is an inherent heterogeneity in foot temperature, influenced by internal and external factors, only healthy subjects were considered in the present work. In this manner, confounding factors causing thermal asymmetry, such as inflammation, ulceration and/or infection, were reduced. Nevertheless, the age range of the healthy volunteers considered may not be representative of the profiles of patients affected with diabetic foot disorders as they are usually older populations. Furthermore, although within the volunteers only one presented deformations on the feet, none presented partial amputations. Therefore, these common features of the intended pathology were not considered for the extraction of the model. Despite these limitations, the methodology presented provides promising results regarding the detection of thermal asymmetry within the entire plantar aspects of both feet. Future work will be aimed at applying the presented methodology to patients with diabetes at risk of developing foot disorders as well as exploring alternative methods to assess the variation in the temperature patterns.

Conclusions
A non-invasive methodology, based on the thermography of the entire plantar aspects of both feet, was established to identify areas of higher variability and track their progress via longitudinal monitoring. Moreover, a standard model was extracted to provide a healthy foot morphological model, which can be taken as a template for further analysis of temperature patterns in diabetic foot disorders. Consequently, the comparison of the contralateral foot can be readily performed on the same subject as well as to the average healthy foot model. Most importantly, the extracted model can successfully handle the presence of deformities and partial amputations usually exhibited by patients affected by the above mentioned condition. Despite the heterogeneity observed in the contralateral analysis of healthy volunteers, the usefulness of the presented methodology has been demonstrated. Since a single measurement, solely using the mean values, may not provide enough information to assess foot asymmetry. Particularly, the methodology proposed provides a fast and simple monitoring tool to be used for personalized medical diagnosis. Funding: This research was funded by the IACTEC Technological Training program, grant number TF INNOVA 2016-2021. This work was completed while Abián Hernández was a beneficiary of a pre-doctoral grant given by the "Agencia Canaria de Investigacion, Innovacion y Sociedad de la Información (ACIISI)" of the "Consejería de Economía, Conocimiento y Empleo" of the "Gobierno de Canarias", which is partly financed by the European Social Fund (FSE) (POC 2014-2020, Eje 3 Tema Prioritario 74 (85%)).

Institutional Review Board Statement:
Not applicable as no clinical trial has been conducted and nor has a study involving patients.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study. Furthermore, acquired data were codified and anonymized to ensure subject confidentiality and data protection.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy restrictions.