- Research
- Open Access

# Quantitative assessment of pain-related thermal dysfunction through clinical digital infrared thermal imaging

- Christophe L Herry
^{1}Email author and - Monique Frize
^{1, 2}

**3**:19

https://doi.org/10.1186/1475-925X-3-19

© Herry and Frize; licensee BioMed Central Ltd. 2004

**Received: **07 April 2004

**Accepted: **28 June 2004

**Published: **28 June 2004

## Abstract

### Background

The skin temperature distribution of a healthy human body exhibits a contralateral symmetry. Some nociceptive and most neuropathic pain pathologies are associated with an alteration of the thermal distribution of the human body. Since the dissipation of heat through the skin occurs for the most part in the form of infrared radiation, infrared thermography is the method of choice to study the physiology of thermoregulation and the thermal dysfunction associated with pain. Assessing thermograms is a complex and subjective task that can be greatly facilitated by computerised techniques.

### Methods

This paper presents techniques for automated computerised assessment of thermal images of pain, in order to facilitate the physician's decision making. First, the thermal images are pre-processed to reduce the noise introduced during the initial acquisition and to extract the irrelevant background. Then, potential regions of interest are identified using fixed dermatomal subdivisions of the body, isothermal analysis and segmentation techniques. Finally, we assess the degree of asymmetry between contralateral regions of interest using statistical computations and distance measures between comparable regions.

### Results

The wavelet domain-based Poisson noise removal techniques compared favourably against Wiener and other wavelet-based denoising methods, when qualitative criteria were used. It was shown to improve slightly the subsequent analysis. The automated background removal technique based on thresholding and morphological operations was successful for both noisy and denoised images with a correct removal rate of 85% of the images in the database. The automation of the regions of interest (ROIs) delimitation process was achieved successfully for images with a good contralateral symmetry. Isothermal division complemented well the fixed ROIs division based on dermatomes, giving a more accurate map of potentially abnormal regions. The measure of distance between histograms of comparable ROIs allowed us to increase the sensitivity and specificity rate for the classification of 24 images of pain patients when compared to common statistical comparisons.

### Conclusions

We developed a complete set of automated techniques for the computerised assessment of thermal images to assess pain-related thermal dysfunction.

## Keywords

## Background

The skin temperature distribution of a healthy human body exhibits a contralateral symmetry [1]. Temperature distribution that shows asymmetrical patterns is usually a strong indicator of abnormality [2–4], but the converse is not always true since some pathological conditions may exhibit bilateral thermal dysfunction. In such cases other signs of abnormalities in the temperature distribution need to be found [5, 6]. Some nociceptive and most neuropathic pain pathologies are associated with an alteration of the thermal distribution of the human body in the form of hyperthermic or hypothermic regions [5]. Since the dissipation of heat through the skin occurs for the most part in the form of infrared radiation, infrared thermography is the method of choice to study the physiology of thermoregulation and the thermal dysfunction associated with pain. The early literature on medical thermography focused on qualitative interpretation of thermograms; this involved determining abnormal thermal variations of the skin by means of a visual assessment of pseudo coloured or grey-level thermograms with the help of isothermal displays, visual localisation of hot or cold spots, and visual detection of symmetry [7–12].

The task of decrypting thermograms and extracting useful and reliable information was complex, even for highly trained medical thermographers, since it relied upon the subjectivity of the human visual ability to distinguish between variations in intensity levels representing temperature distribution in thermograms. In addition, the use of pseudo-colours for mapping the temperatures of a thermogram was also criticised for its subjectivity due to the psychological effect of certain colours, which may skew the observer's performance [13]. As a result, thermographic research examined general quantification techniques for specific problems in order to reduce the subjectivity of the assessment of thermograms [14]. Many past and recent publications discuss thermal dysfunction associated with pain, however, to our knowledge none so far applied comprehensive computerised techniques to the assessment of thermal images of persons experiencing pain.

## Methods

### Objectives

The overall goal of this work was to automate as much as possible a computerised assessment of thermal images of pain in order to support clinicians' decision making. Our approach consists of several steps. First, the thermal images are pre-processed to reduce the noise introduced during the initial acquisition of the images and to extract irrelevant background. Then, potential regions of interest are identified in a semi-automated manner, using fixed dermatomal subdivisions of the body; they are also identified in an automated manner based on an isothermal analysis and segmentation techniques. Finally, we assess the degree of asymmetry between contralateral regions of interest using statistical computations and distance measures between comparable regions.

### Data collection

### Methodology

#### Pre-processing of thermal images

The characteristics of the thermal imaging system used in this study was assumed to be an almost ideal infrared system, where most of the noise types inherent to the design are lower than the photon noise generated by the background and signal itself. It is known that photon noise follows a Poisson distribution and is signal-dependent [16][17]. Common noise removal techniques based on a Gaussian approximation to Poisson noise do not work well due to the signal-dependent, spatially variant nature of photon noise. In order to remove the noise efficiently while retaining the relevant information contained in thermal images, we adapted a wavelet-based Poisson noise removal technique that was suggested by Nowak and Baraniuk [18]. They considered the successive images with increasing intensity levels from the photon counting process, until they obtain the final image at the end of the observation period. The authors measured the error between a wavelet-domain filtered version of each image with the next image, until the final image. The sum of the errors is called the PREdictive Sum of Squares (PRESS), which is dependent on the filter (or threshold) used. Minimising the PRESS criterion yields the optimal filter (or threshold). The authors also showed that the sub-images are not required in order to compute the optimal threshold, since the latter reduces to:

_{+}sets to zero the negative thresholds.

The noise variance is estimated by projecting the image onto the square of the wavelet and scaling functions, yielding what the authors called a Discrete Wavelet Squared Transform (DWST). The DWST may be implemented efficiently using recursions similar to that of the Fast Wavelet Transform. The choice of the wavelet bases is important to preserve the original features of the thermal images. In particular, they need to be well localised and have a certain degree of smoothness. Possible candidates among orthogonal wavelets are Daubechies wavelets and Battle-Lemarié spline orthogonal wavelets. They have been shown to provide acceptable results in medical image processing.

- 1.
The histogram of the image is computed and smoothed with a Gaussian function to obtain a desired number of modes.

- 2.
The extreme values of the resulting function are then computed as well as their corresponding grey-level.

- 3.
A first threshold is set to the strongest minimum, that is the last minimum remaining after successive smoothing. A second threshold is selected to be the first stable minimum after successive smoothing. The final threshold is then given by a weighted sum of the two thresholds.

- 4.
Isolated pixels wrongly labelled as background pixels are corrected. The approximate boundaries of the image are superposed to the thresholded images and any holes created are filled.

- 5.
The remaining regions of interest are then roughly segmented and the smaller regions are discarded as they likely correspond to background parts. The same procedure is performed on the background regions.

- 6.
Finally the remaining possible holes in the foreground regions are filled based on the size and shape of the regions.

#### Identification of region of interests (ROIs)

In order to compare the results that were derived from our population of normal patients with the values published in the literature [3, 4], we divided the body into multiple squares or rectangles that other authors have used, for example Goodman [3] or Uematsu [4].

The manual division by Uematsu followed the areas of the skin called dermatomes that are innervated by the major peripheral nerves. Being able to link the data from those regions of interest with the corresponding nervous system that feeds them is believed to be more meaningful for the assessment of the thermograms. Goodman used a grid of squares, for the back and front, and rectangles for the arms, legs and extremities. The grid was adapted to the bodily dimensions of each subject and could be moved around to cover any desired areas. The statistics were computed for various numbers of slices or frames of the unit square or rectangle and results from the analysis compared to that of Uematsu's. For comparison purpose, we devised a hybrid of the two divisions for the back, chest and legs, which consists in Uematsu's division with two additional regions in the back.

- 1.
The vertical symmetry line separating contralateral parts of the body is first sought by detecting the edges and averaging the middle distance between the outer boundaries of the body.

- 2.
For the right and left sides of the body, we defined the function corresponding to the distance between the outer boundaries and the symmetry line, thus obtaining two functions representing the left and right profiles of the body.

- 3.
The profile resulting from averaging both left and right profiles is smoothed with a Gaussian function and we look for extreme values of the resulting function.

- 4.
The first significant maximum is found to be representative of the extremity of the shoulder and the level corresponding to 80% of the maximum is used to find the coordinate of the extreme point for choosing the ROIs.

- 5.
The flank level is found through a similar method for the upper back images but looking for the first significant minimum of the pro-file function instead of the maximum. For the lower back images, a search for the minimum distance between left and right side, above the upper point of the buttock cleft, provides acceptable results.

- 6.
The buttock cleft is found through an analysis of the horizontal profile lines across the lower part of the images, when a middle edge line is present. The distance between the lower point and higher points is computed for each profile lines, yielding a function of this distance. The first minimum of this function is chosen to be approximately the upper extremity of the buttock cleft.

- 7.
Finally, the knee is easily detected by looking at the minimal distance between each edge lines for both legs.

Based on the previous landmarks, the positioning of the various boxes delimiting the ROIs defined by both Uematsu and Goodman was determined in a straightforward manner. A second method was implemented to determine more specific regions of interest than dermatomes. It is based on an isothermal analysis of the entire pre-processed image, to determine hotter and colder regions. The step for the computation of each isotherm is chosen according to the range of intensity values covered in the image to provide some adaptivity to the generation of isotherms. In practice, the step is set to one tenth of the total range of intensity values.

Also, whenever the sensitivity setting is available, we looked for isothermal regions at regular intervals of 0.5°C, for comparison with common isothermal analysis found in the literature. The same rule is applied for the hotter regions, starting from the maximum and subtracting successive increments of the step. Typically, only the first two to three steps are necessary to yield significant regions. In addition, smaller regions of less than 10 pixels were discarded by thresholding the areas of the regions returned by the previous isothermal computation. Statistical comparisons with regions smaller than the threshold area would have limited significance in most cases.

Moreover, whenever two regions of very different sizes (one large and one small) are found to be too close to each other, they are merged into one larger region. The positions and pixels in the remaining regions are kept for the statistical analysis.

Finally, in order to avoid artifact regions located around the contour of the body area of interest, the image is eroded along the boundaries using a morphological eroding operator. For all methods, a simple algorithm is used to compute the reflection of the regions of interest with respect to the vertical line of symmetry, allowing us to assess the asymmetry of the resulting intensity/temperature distribution.

#### Asymmetry analysis

We considered a comprehensive set of statistical features that may be able to capture more completely the nature of the distributions in the regions of interest. Those features include the classical moments (mean, variance, skewness and kurtosis), the extreme values, the area, the heat content defined by Montoro and Anbar as the area times the mean temperature [20] and finally the entropy of the region, which may be approximated as:

where *p*
_{
k
}is the number of pixels in each grey-level *k* ∈ [1, *L*], where L is the number of grey levels, which we chose to be 256.

The outcome for a particular region of interest usually relies on the computation of the difference of statistics between this region and its contralateral counterpart. Some authors also considered the statistics of the pixel-by-pixel difference between two comparable regions [21]. However, we felt that both these approaches were omitting important parameters of the temperature distribution of ROIs, since they looked at a few statistical features derived from assumptions on the nature of this distribution. As a result, we proposed a new statistical approach using distances between histograms of comparable ROIs in order to assess their degree of similitude as it is extensively used in image registration and pattern recognition. The following distance measures were considered: Manhattan or absolute distance, Euclidian distance, maximum distance, chi-squared distance, Jeffrey-divergence distance and Mallows distance [22, 23]. The first three distances are simply the first three moments of the classical Minkowski distance based on the metric of the same name. The general formula for the Minkowski distance is:

where H and G denote the histograms and *p* = 1 for the Manhattan distance, *p* = 2 for the Euclidian distance and *p* = ∞ for the maximum distance. The chi-square distance is given by:

The Jeffrey divergence is defined by:

The Mallows distance is:

where *H*
_{(i)}and *G*
_{(i)}are the sorted histogram values and *p* was chosen to be 2.

In addition, we also calculated the Kolmogorov-Smirnov statistic, which is defined as the maximal distance between two cumulative distributions:

where *H*
_{
cum
}and *G*
_{
cum
}are the cumulative histograms corresponding to *H* and *G* respectively. The Kolmogorov-Smirnov test derived from this statistic was performed for each pair of symmetric regions of interest, since it showed encouraging results in Vavilov *et al*.'s paper [24]. This test has the advantage of being non-parametric, that is, it does not assume that the population of pixels within each region follows a specific distribution.

## Results and discussion

### Performance of the denoising method

The Discrete Wavelet Squared Transform (DWST) was implemented in a relatively straightforward manner using the one-dimensional and two-dimensional wavelet transforms available in the Wavelab library from the department of Statistics of Stanford University [25]. In order to reduce the memory requirements of the direct algorithm, we chose to decompose each 128 × 128 thermal image into four 64 × 64 sub-images and remove the noise from each of these sub-images. In practice, it was found that artifacts affect the quality of the resulting image at the border of each sub-image.

A few orthogonal wavelet bases were investigated on thermal images from our database with various noise levels. Nowak and Baraniuk [18] suggested using a simple Haar basis or a Daubechies wavelet basis for nuclear medicine images. Our tests on thermal images showed that the best results quantitatively (based on the reduction in noise variance) and visually (presence of artifacts and ringing effect for instance) were achieved with a Battle-Lemarié spline wavelet basis of the 3rd order.

Furthermore, we varied the parameter a in equation 1 and found that α > 1 yielded better denoising results. In particular, 4.5 < α < 7.5 consistently produced the best trade off between reduction of noise and conservation of the image structure.

Since the computation of mean square errors between original images and denoised images was not readily available, the performance of the denoising was assessed in terms of percentage of edges and contours remaining after denoising, for a comparable noise level. We found that the resulting images from the Press-optimal wavelet-domain filter retained the contours and edges of the thermal images better than when a Fourier domain Wiener filter or a wavelet-domain filter (applied to the square root of the image) was used.

In order to determine the percentage of edges remaining after denoising we summed the total number of pixels returned by a Canny edge detector for increasing threshold level, up to a level where no edge can be detected. We assumed that the Canny edge detector was not affected significantly by the noise and therefore we normalised the number of pixels to that of the noisy image.

### Background extraction efficiency

Rate of incorrect background extraction. Rate of incorrect background extraction for images of the lower back, upper-back and legs. The first figure indicates the number of images whose background was incorrectly extracted. The second figure is the total number of images.

With Denoising | Without Denoising | |
---|---|---|

Lower Back images | 9/72 | 22/72 |

Upper Back images | 12/85 | 15/85 |

Legs | 22/205 | 30/205 |

### Automated selection of ROIs

The automated selection of ROIs using a compromise between Uematsu's technique and the more complete technique by Goodman was fairly successful. The determination of the limits for the ROIs based on the shape of the body and on its intensity distribution enabled us to approximately place the boxes corresponding to each ROI. A simple iterative procedure was performed on the boxes to make sure they do not contain any background. However, the difference in body shapes and the lack of symmetry of body areas in some thermal images made it difficult to adapt the ROIs to cover the best area.

Furthermore, we assumed that thermal images would be symmetrical with respect to a vertical line of symmetry; this was verified in most thermal images. However, other thermal images showed areas of the body bent toward one side of the border of the image, or rotated. As a result, even if the limits for the ROIs were correctly identified, the final ROIs did not cover the wanted areas. It may therefore be necessary to consider alternative symmetry curves and define a proper reflection method to get corresponding pixels from contralateral regions.

In addition, we found that division into rectangular boxes based on the dimension of the body may fail to include an abnormality if the latter happens to be in the middle of two boxes. This outlines the necessity for a complementary ROI delimitation method such as the isothermal analysis that works well for asymmetric images and targets potential abnormalities. It produced similar results in most cases: targeted ROIs over possible abnormalities in the thermal pattern distribution in the form of hot or cold spots.

ROIs covering both contralateral parts of the body and ROIs along the line of symmetry were flagged for later analysis, ensuring that the assessment does not fail to detect bilateral abnormalities. A strong gradient of temperature between a flagged region and its surrounding area indicates an abnormality according to our procedure.

We also compared the ROIs produced by the isothermal analysis of fully pre-processed images and of the noisy images with the extracted background. We found that noisy images tend to return a greater number of smaller ROIs that would otherwise be merged into larger ROIs in the case of denoised images. However, smaller regions may not necessarily reflect a true elevation or decrease of temperature but instead indicate spurious, noisy pixels. This justifies discarding regions smaller than a threshold, which was chosen to be 20 pixels for our study. In addition, statistical comparisons of very small samples may not be significant, which further justifies our decision.

### Some performance aspects of thermographic analysis

The performance of our quantitative analysis of thermal images in terms of complexity and CPU time was greatly affected by the denoising stage and by our choice of the PRESS-Optimal filter. For instance, the DWST theoretically required *O*(21^{3} × 128^{2}) computations per resolution level, but our implementation using smaller sub-images increased the complexity. The average CPU time to process a single image was about 4 minutes, using a Pentium III at 1.2 GHz with 512 MB of memory, MATLAB release 12.1 and Wavelab 802 [25]. In comparison, the rest of the analysis (background extraction, identification of ROIs, asymmetry analysis) required typically less than 3 seconds. Although it is difficult to assess the overall complexity of the thermographic analysis, it is certainly heavily weighted by the denoising step. The other denoising techniques considered require less CPU time than the PRESS-Optimal method, but do not perform as well. If we assume that the CPU time is a decisive factor in the analysis of large databases of thermal images, it would be desirable to improve our implementation of the PRESS-Optimal method or investigate other denoising methods that perform equally well.

### Thermal asymmetry of the control population

Normal Values for the back. Mean Temperature Difference plus or minus standard deviation for our control population of healthy volunteers.

ROIs | Number of Cases | Mean Temperature Difference ± Standard Deviation |
---|---|---|

Thoracic (medial) | 12 | 0.09 ± 0.09 |

Shoulder (posterior) | 12 | 0.12 ± 0.10 |

Thoracic (lateral) | 12 | 0.09 ± 0.09 |

Thigh (anterior) | 13 | 0.15 ± 0.11 |

Thigh (posterior) | 14 | 0.11 ± 0.09 |

Knee (anterior) | 13 | 0.12 ± 0.13 |

Knee (posterior) | 14 | 0.15 ± 0.12 |

Leg (anterior) | 13 | 0.15 ± 0.12 |

Leg (posterior) | 14 | 0.10 ± 0.06 |

All the mean temperature difference values are well within the limits for normal value specified by Ue-matsu and Goodman. Our values are in fact sig-nificantly lower, by the order of 0.1°C, which may partly be explained by the relatively small number of cases in our control population. This certainly had an effect on the standard deviation of the mean temperature differences. Another reason for these low values may be the size of the regions of interest, which for practical reasons was smaller than that of the previous authors.

Distance results for the back. Mean, standard deviation and threshold of distance measures for the back, for our control population. Mah., Eucl., Max., Chi, JD, Mal, KS stand for: Mahalanobis, Euclidian, Maximum, Chi-square, Jeffrey-Divergence, Mallows, Kolmogorov-Smirnov distances respectively. *CI* stands for confidence interval.

Distance results: noisy image | |||||||
---|---|---|---|---|---|---|---|

Mahn. dist. | Eucl. dist. | Max. dist. | Chi sq. dist. | JD dist. | Mallows dist. | KS dist. | |

Mean | 0.492 | 0.083 | 0.032 | 0.110 | 0.084 | 1.092 | 0.171 |

Std. dev. | 0.152 | 0.025 | 0.013 | 0.057 | 0.048 | 0.571 | 0.103 |

Threshold (99%CI) | 0.796 | 0.133 | 0.057 | 0.223 | 0.180 | 2.235 | 0.377 |

Distance results: denoised image | |||||||

Mahn. dist. | Eucl. dist. | Max. dist. | Chi sq. dist. | JD dist. | Mallows dist. | KS dist. | |

Mean | 0.565 | 0.105 | 0.042 | 0.142 | 0.108 | 1.427 | 0.205 |

Std. dev. | 0.180 | 0.039 | 0.020 | 0.070 | 0.060 | 0.798 | 0.111 |

Threshold (99%CI) | 0.924 | 0.184 | 0.082 | 0.282 | 0.229 | 3.022 | 0.427 |

Distance results for the legs. Mean, standard deviation and threshold of distance measures for the legs, for our control population. Mah., Eucl., Max., Chi, JD, Mal, KS stand for: Mahalanobis, Euclidian, Maximum, Chi-square, Jeffrey-Divergence, Mallows, Kolmogorov-Smirnov distances respectively. *CI* stands for confidence interval.

Distance results: noisy image | |||||||
---|---|---|---|---|---|---|---|

Mahn. dist. | Eucl. dist. | Max. dist. | Chi sq. dist. | JD dist. | Mallows dist. | KS dist. | |

Mean | 0.663 | 0.099 | 0.037 | 0.195 | 0.101 | 0.509 | 0.181 |

Std. dev. | 0.150 | 0.022 | 0.013 | 0.075 | 0.030 | 0.247 | 0.091 |

Threshold (99%CI) | 0.964 | 0.142 | 0.063 | 0.345 | 0.160 | 1.004 | 0.364 |

Distance results: denoised image | |||||||

Mahn. dist. | Eucl. dist. | Max. dist. | Chi sq. dist. | JD dist. | Mallows dist. | KS dist. | |

Mean | 0.694 | 0.106 | 0.038 | 0.206 | 0.109 | 0.528 | 0.192 |

Std. dev. | 0.166 | 0.025 | 0.014 | 0.083 | 0.035 | 0.279 | 0.100 |

Threshold (99%CI) | 1.027 | 0.155 | 0.066 | 0.373 | 0.178 | 1.086 | 0.392 |

Results from noisy and denoised images did not differ significantly and confidence intervals were overlapping. This also confirms that the denoising does not affect the statistical properties of thermal images and thus does not induce a loss of critical information, while we showed that it facilitated the previous analysis.

Using the threshold values from table 3 we assessed the efficiency of each distance measure to discriminate between normal and abnormal regions of interest. For our normal population, this was accomplished by looking at the distances generated from a complete isothermal analysis of thermal images of the back and legs. The distance with the least number of misclassifications may be seen as the more discriminative with respect to normal regions, that is the distance with the best specificity (Number of true negatives with respect to the number of false positives) or the least number of regions identified as positive or abnormal.

Specificity values for the control population. Specificity values of statistical distance measures considered, from a complete isothermal analysis of thermograms of the back and legs for our control population.

Specificity for thermograms of the back | |||||||
---|---|---|---|---|---|---|---|

Mahn. dist. | Eucl. dist. | Max. dist. | Chi sq. dist. | JD dist. | Mallows dist. | KS dist. | |

Specificity | 0.95 | 0.95 | 0.93 | 0.95 | 0.94 | 0.95 | 0.86 |

Specificity for thermograms of the legs | |||||||

Mahn. dist. | Eucl. dist. | Max. dist. | Chi sq. dist. | JD dist. | Mallows dist. | KS dist. | |

Specificity | 0.97 | 0.99 | 0.98 | 0.97 | 0.98 | 0.98 | 0.95 |

Finally, we noted that contrary to Vavilov et al. who found that the Kolmogorov-Smirnov statistics allowed the accurate classification of normal regions from abnormal and potentially pathological regions, our experiment suggested that other distance measures might be more suitable to this task. All the other distance measures have similar specificity values. In order to validate our approach and draw some conclusions on the performance of our novel statistical tool for the assessment of pain in thermal images, we first analysed a set of 72 thermal images of the upper back with the fixed division into rectangular ROIs. Then the usual statistics were computed as well as distance between distributions of comparable ROIs.

Out of 222 identified regions of interest (and their symmetric parts when available), 33 regions were classified as abnormal using the Manhattan distance, 38 for the Euclidian distance, 27 for the Maximum distance, 36 for the Chi-square distance, 30 for the Jeffrey-Divergence distance, 39 for the Mallows distance and 37 for the Kolmogorov-Smirnov distance. The number of abnormal regions produced by the comparison of mean, variance, skewness, kurtosis, maximum and minimum was only 6. The abnormal regions were not necessarily consistent across the distance measures. The abnormal regions produced by the Euclidian distance, the Chi-squared distance, the Kolmogorov-Smirnov distance and to a lesser extent those produced by the Manhattan and Jeffrey-divergence distances were generally equivalent. The abnormal regions produced by the three other methods were relatively different. It is interesting to note that the comparison of the mean, variance, skewness, kurtosis, maximum and minimum returned a very low number of abnormal regions as opposed to the distance measures. The threshold for the mean was based on the values of table 2. The other parameters were used to determine whether or not the comparison of the mean was meaningful.

Furthermore, we noted that the performance of the Mallows distance seemed to be dependent on the size of the sample, where the sample was less than a hundred pixels typically. This is especially problematic when used in conjunction with the isothermal analysis, which may return several small regions. The Jeffrey-Divergence suffered from a similar drawback and as a result we need to consider one of the other distances for very small samples. The Minkowski distances, the chi-square distance and the Kolmogorov-Smirnov distance performed equally well when tested on small samples provided by the isothermal analysis. The five distances correctly identified small abnormal regions for a series of 5 images with apparent abnormalities, while the Jeffrey-Divergence and Mallows distances missed 4 out of 5.

Sensitivity and specificity values for actual pain patients. Sensitivity and specificity values for thermograms of actual pain patients. Mean Δ*T* refers to the comparison based on the mean, variance, skewness, kurtosis, maximum and minimum.

Mahn. dist. | Eucl. dist. | Max. dist. | Chi sq. dist. | JD dist. | Mallows dist. | KS dist. | Mean Temp. diff. | |
---|---|---|---|---|---|---|---|---|

Sensitivity | 0.67 | 0.78 | 0.33 | 0.72 | 0.67 | 0.55 | 0.71 | 0.5 |

Specificity | 0.67 | 0.83 | 0.83 | 0.67 | 0.67 | 0.5 | 0.43 | 0.67 |

These results confirm that the sensitivity of the Euclidian, Chi-square and Kolmogorov-Smirnov distance seem to be significantly higher than that of other distances. Also, the specificity of the Kolmogorov-Smirnov distance is lowest, again confirming the results from the analysis of the control population. Based on the assessment of these 24 images, the best distance appeared to be the Euclidian distance with a fairly high sensitivity and specificity. The Manhattan distance, the Mallows distance and the basic statistics method had relatively low sensitivity. This suggests these methods would dismiss as normal a number of patients with abnormal patterns. Both Manhattan distance and Mallows distance did not perform as well as the results of table 5 suggests.

## Conclusions

This paper investigated the computerised assessment of pain through digital infrared thermal imaging. In particular, we looked at the overall digital processing of medical thermograms, from the preprocessing requirements to the decision-support system inputs. In an attempt to increase the objectivity and to facilitate the assessment of pain by medical specialists, we automated as much as possible the assessment of thermal images. The automation process included the pre-processing and identification of the regions of interest, a statistical analysis and a set of features that can be used by a physician to make a diagnosis.

The pre-processing of thermal images had been overlooked in the medical thermography literature, which focused on the statistical and asymmetry analysis. In order to facilitate the identification of the regions of interest, we proposed applying denois-ing techniques to thermal images. We modelled the noise by a Poisson distribution with an additive Gaussian white noise component easily removed by classical denoising techniques. Then, we identi-fied the need for a noise removal method that takes the Poisson nature of the noise (and of the signal) into account. This was achieved through a Poisson noise removal technique in the Wavelet domain called PRESS-optimal that was available and used for nuclear medicine images [18]. We showed that the resulting denoised image was more easily processed in the subsequent analysis, without any sig-nificant loss of statistical information.

The next problem was to remove the background in order to focus on the areas of interest, again facilitating the following identification of the regions of interest for the statistical comparison. We proposed a simple yet efficient histogram thresholding technique with additional morphological and logical steps, which performed well on both our databases of thermal images.

The identification of the regions of interest consisted of an automated division into a rectangular grid following Uematsu et al. and Goodman et al.'s papers [3, 4] for comparison with our control population of normal healthy subjects. This was done by identifying specific reference points on the body, through an analysis of contours and profiles across the image. In addition, an isothermal analysis was performed to determine specific and localised regions, especially the hotter and colder ones that may be indicative of abnormalities.

Finally, we proposed a new statistical tool for the classification of regions of interest as normal or abnormal, based on the computation of several distances between histograms of the ROIs. The results of the quantitative assessment of thermal images from our control population and from pain patients showed that the method based on the Euclidian distance seemed to outperform the other statistical methods considered. Overall, we covered many of the most important steps in the computerised assessment of thermal images and removed the user input, thus providing a highly automated assessment of pain in thermal images.

## Authors contributions

CH conceived the algorithms and performed the data analysis and participated in the study design and coordination. MF participated in the study design and coordination and provided the database of images and expertise in infrared imaging. All authors read and approved the final manuscript.

## Declarations

### Acknowledgements

The authors thank the National Science and Engineering Research Council of Canada (NSERC) for funding this study. The authors also thank The Terry Fox Foundation in Fredericton, NB for funding the 1984 study and Dr. G. Quartey, at Moncton Hospital, NB, for his expertise in pain assessment, patient base, and clinical results.

## Authors’ Affiliations

## References

- Houdas Y, Ring EFJ, Eds:
*Human Body Temperature: its Measurement and Regulation*New York, NY: Plenum Press 1982.Google Scholar - Aubry-Frize M, Quartey GRC, Evans H, LaPalme D:
**The Thermographic Detection of Pain.***In Proceedings of the 3rd Canadian Clinical Engineering Conference: September 1981; Saskatoon, Canada*1981, 82–83.Google Scholar - Goodman PH, Murphy MG, Siltanen GL, Kelley MP, Rucker L:
**Normal Temperature Asymmetry of the Back and Extremities by Computer-Assisted Infrared Imaging.***Thermology*1986,**1:**195–202.Google Scholar - Uematsu S, Edwin DH, Jankel WR, Kozikowski J, Trattner M:
**Quantification of Thermal Asymmetry, Part 1: Normal Values and Reproducibility.***Journal of Neurosurgery*1988,**69:**552–555.View ArticleGoogle Scholar - Hooshmand H, Hashmi M, Phillips EM:
**Infrared Thermal Imaging as a Tool in Pain Management – An 11 Year Study. Part I of II.***Thermology International*2001,**11**(2):53–65.Google Scholar - Hooshmand H, Hashmi M, Phillips EM:
**Infrared Thermal Imaging as a Tool in Pain Management – An 11 Year Study. Part II: Clinical Applications.***Thermology International*2001,**11**(2):117–129.Google Scholar - Ellis WV, Morris JM, Swartz AA: Screening Thermography of Chronic Back Pain Patients with Negative Neuromusculoskeletal Findings. Thermology 1989.,3(2):Google Scholar
- Goldberg GS: Infrared Imaging and Magnetic Resonance Imaging Correlated in 35 Cases. Thermology 1986.,1(4):Google Scholar
- Hubbard JE, Hoyt C:
**Pain Evaluation in 805 Studies by Infrared Imaging.***Thermology*1986,**1**(3):161–166.Google Scholar - Kim YS, Cho YE:
**Correlation of Pain Severity with Thermography.***In Proceedings of the 17th Annual International Conference of the IEEE Engineering in Medicine and Biology Society and 21st Canadian Medical and Biological Engineering Conference*1995,**2:**1699–1700.MathSciNetGoogle Scholar - LeRoy PL:
**Update in Interventional Digitalized Infra-Red thermal Imaging in Pain Management.***In Proceedings of the 17th Annual International Conference of the IEEE Engineering in Medicine and Biology Society and 21st Canadian Medical and Biological Engineering Conference*1995,**2:**1707–1708. 10.1109/IEMBS.1995.579902MathSciNetGoogle Scholar - Thomas D, Siahamis G, Marion M, Boyle C:
**Computerised Infrared Thermography and Isotropic Bone Scanning in Tennis Elbow.***Ann Rheum Dis*1992,**51:**103–107.View ArticleGoogle Scholar - Anbar M:
**Objective Assessment of Clinical Computerized Thermal Images.***In Medical Imaging V: Image Processing, SPIE-The International Society for Optical Engineers*1991,**1445:**479–484.Google Scholar - Collins AJ, Ring EFJ, Cosh JA, Bacon PA:
**Quantitation of Thermography in Arthritis using Multiisothermal Analysis.***Ann Rheum Dis*1974,**33:**113–115.View ArticleGoogle Scholar - Pochaczevsky R:
**Technical Guidelines, Edition 2.***Thermology*1986,**2**(2):108–112.Google Scholar - Wolfe WL, Ed:
*Handbook of military infrared technology*Washington DC: Office of Naval Research, Department of Navy 1965.Google Scholar - Keyes RJ, Ed:
*Optical and infrared detectors. No. 19 in Topics in applied physics*Springer-Verlag 1977.Google Scholar - Nowak RD, Baraniuk RG:
**Wavelet-Domain Filtering for Photon Imaging Systems.***IEEE Transactions on Image Processing*1999,**8**(5):666–678. 10.1109/83.760334View ArticleGoogle Scholar - Tsai DM:
**A Fast Thresholding Selection Procedure for Multimodal and Unimodal Histograms.***Pattern Recognition Letters*1995,**16**(6):653–666. 10.1016/0167-8655(95)80011-HView ArticleGoogle Scholar - Montoro J, Anbar M:
**New Modes of Data Handling in Computerized Thermography.***In Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society*1988, 845–847.View ArticleGoogle Scholar - Mabuchi K, Chinzei T, Fujimasa I, Haeno S, Abe Y, Yonezawa T:
**An Image-Processing Program for the Evaluation of Asymmetrical Thermal Distributions.***In Proceedings of the 19th Annual International Conference of the IEEE Engineering in Medicine and Biology Society; Chicago, IL*1997, 725–728.Google Scholar - Puzicha J, Buhmann JM, Rubner Y, Tomasi C:
**Empirical Evaluation of Dissimilarity Measures for Color and Texture.***In Proceedings the IEEE International Conference on Computer Vision (ICCV1999)*1999, 1165–1173.View ArticleGoogle Scholar - Levina E, Bickel P:
**The Earth Mover's Distance is the Mallows Distance: Some Insights from Statistics.***In Proceedings of the Eight International Conference on Computer Vision (ICCV2001)*2001,**2:**251–256. 10.1109/ICCV.2001.937632View ArticleGoogle Scholar - Vavilov VP, Vavilova EV, Popov DN:
**Statistical Analysis of the Human Body Temperature Asymmetry as the Basis for Detecting Pathologies by Means of IR Thermography.***In Thermosense XXIII**(Edited by: Rozlosnik AE, Dinwiddie RB).*SPIE–The International Society for Optical Engineers 2001,**3061:**482–491.View ArticleGoogle Scholar -
**Wavelab 802: Matlab wavelet toolbox from the Department of Statistics, Stanford University**[http://www-stat.stanford.edu/~wavelab/]

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article: verbatim copying and redistribution of this article are permitted in all media for any purpose, provided this notice is preserved along with the article's original URL.