Examining an Image Reconstruction Method in Infrared Emission Tomography

TL;DR

This paper introduces RISE, a new infrared tomography reconstruction method that adapts photon emission tomography techniques using Monte Carlo simulations, validated through simulated and experimental infrared data.

Contribution

The paper presents RISE, a novel reconstruction method applying statistical physics and Monte Carlo techniques to infrared tomography, expanding established emission tomography methods to infrared radiation.

Findings

RISE performs comparably to ART and MLEM in simulations.

Experimental infrared projections successfully reconstructed with RISE.

Validated the applicability of photon emission tomography principles to infrared imaging.

Abstract

We present and evaluate the application of the "Reconstructed Image from Simulations Ensemble" (RISE), a novel tomographic image reconstruction method, in infrared tomography. We demonstrate that established methods of photon emission tomography, widely used with penetrating ionizing radiation, are applicable to infrared radiation. RISE, the method of choice, employs statistical physics concepts and utilizes Monte Carlo techniques to construct the imaged object from its infrared planar projections. The validity of the InfraRed Emission Tomographic (IRET) method is demonstrated, and the efficacy of RISE is evaluated with A) simulated data and B) experimental sets of infrared projections obtained from a thermal phantom with an infrared camera. For the simulation studies presented, the reconstructed images obtained with RISE and the well - known Algebraic Reconstruction Technique (ART) and Maximum Likelihood Expectation Maximization (MLEM) method were evaluated using well-established metrics.

Tables

Table 1. Table 1: Descriptors of the Shepp-Logan variant - ”Phantom C” - shown in Figure 4 .

”Object”	Temperature $T_0$	Position $(x,y)$	Size $(u,v)$	Orientation $\theta$
	(${}_{}{}^{\circ }C$)	(pixel units)	(pixel units)	(degrees)
a	2	(42.2, 34.0)	(7.9, 2.8)	72
b	2	(22.0, 34.0)	(10.5, 4.1)	108
c	3	(32.1, 48.0)	(6.4, 5.4)	0
d	4	(32.1, 30.0)	(1.2, 1.2)	-
e	3	(32.1, 9.8)	(3.0, 1.7)	0

Table 2. Table 2: CC, NMSE, SSIM, PSNR and CNR scores of the reconstructed images presented in Figure 5 .

		Phantom A			Phantom B
	RISE	ART	MLEM	RISE	ART	MLEM
CC	0.99	0.55	0.79	0.95	0.52	0.70
NMSE	0.01	0.81	0.49	0.05	1.19	0.75
SSIM	0.77	0.30	0.44	0.55	0.22	0.39
PSNR	33.25	16.53	12.78	32.33	16.41	11.24
CNR	6.52	2.24	3.17	13.55	1.97	3.11

Table 3. Table 3: CC, NMSE, SSIM, PSNR and CNR scores of the RISE, ART, and MLEM images (shown in Figure 6 ) reconstructed in the second simulation study (”Resolving of Temperature Differences”).

$T_0$ ${(}^{o}C)$		RISE	ART	MLEM
3.0	CC	0.96	0.55	0.80
	NMSE	0.03	0.89	0.47
	SSIM	0.65	0.22	0.45
	PSNR	28.60	16.93	12.90
	CNR	6.32	2.33	3.30
2.0	CC	0.95	0.48	0.80
	NMSE	0.05	1.18	0.48
	SSIM	0.62	0.21	0.43
	PSNR	26.78	17.95	13.24
	CNR	6.28	2.07	3.41
1.0	CC	0.86	0.46	0.77
	NMSE	0.16	1.35	0.53
	SSIM	0.52	0.22	0.42
	PSNR	23.80	18.98	14.01
	CNR	5.04	1.71	2.50

Table 4. Table 4: CC, NMSE, SSIM, PSNR, and CNR scores of the reconstructed images presented in Figure 7 .

Phantom		RISE	ART	MLEM
C	CC	0.99	0.67	0.92
	NMSE	0.02	0.50	0.13
	SSIM	0.94	0.42	0.60
	PSNR	29.68	15.31	21.24
	CNR	5.39	1.66	3.28
D	CC	0.98	0.56	0.83
	NMSE	0.01	0.46	0.18
	SSIM	0.82	0.36	0.49
	PSNR	27.51	12.64	16.76
	CNR	2.00	0.18	0.08
E	CC	0.99	0.73	0.96
	NMSE	0.01	0.27	0.20
	SSIM	0.77	0.44	0.66
	PSNR	29.46	13.98	15.31
	CNR	3.68	0.33	0.43

Equations

$$I(T(x,y),\lambda ,x,y)=ϵ(\lambda )\cdot \frac{c_1}{\pi {\lambda }^5}\cdot (exp(\frac{c_2}{\lambda T(x,y)})-1{)}^{-1}$$

$$t=x\cdot \cos (\theta )+y\cdot sin(\theta )$$

$$s=-x\cdot \sin (\theta )+y\cdot cos(\theta )$$

$$R(\lambda,t,\theta)=\int_{s_{1}}^{s_{2}}I(T(s,t),\lambda,s,t)~{}e^{\bigg{(}{-\int_{s}^{s_{2}}{k(\lambda,s^{\prime},t)ds^{\prime}}\bigg{)}}}~{}ds$$

$$R_0(t,\theta )=\Delta \tau \Delta S{\int }_{{\lambda }_1}^{{\lambda }_2}R(\lambda ,t,\theta )d\lambda$$

$$R_i=\sum _{j=1}^{N^2}{P}_{ij}{F}_j$$

$$P_{ij}={\tilde{P}}_{ij}exp(-k_u{r}_{ij})$$

$$T(x,y)=\sum _{i=0}^{M_z}\sum _{j=-i}^i{C}_i^j{Z}_i^j(x,y)+\sum _{i=1}^{N_s}{A}_i(exp(\frac{r_i-r_{0i}}{s_i{r}_{0i}})+1{)}^{-1}$$

$$r_0(\theta )=\frac{u\cdot v}{\sqrt{(v\cdot cos(\varphi -\theta ){)}^2+(u\cdot sin(\varphi -\theta ){)}^2}}$$

$${\chi }^2=\sum _i\frac{(R_i-Y_i{)}^2}{{ϵ}_i^2}$$

$$\text{BIC}={\chi }_{min}^2(N_s)+(7\cdot N_s+1)\cdot log(M)$$

$$CC=\frac{{\sum }_{i=1}^{N^2}(F_i^a-\overline{F^a})(F_i^b-\overline{F^b})}{\sqrt{{\sum }_{i=1}^{N^2}(F_i^a-{\bar{F}}^a{)}^2{\sum }_{i=1}^{N^2}(F_i^b-{\bar{F}}^b{)}^2}}$$

$$NMSE=\frac{{\sum }_{i=1}^{N^2}(F_i^a-F_i^b{)}^2}{{\sum }_{i=1}^{N^2}{F_i^b}^2}$$

$$PSNR=10log_{10}\bigg{(}\frac{N^{2}\cdot max(F^{b})^{2}}{\sum_{i=1}^{N^{2}}(F^{a}_{i}-F^{b}_{i})^{2}}\bigg{)}$$

$$CNR=\frac{\bar{T}-\bar{B}}{{\sigma }_B}$$

$$SSIM=l\cdot c\cdot s$$

$$l=\frac{2{\bar{F}}^a{\bar{F}}^b+C_1}{{\overline{F^a}}^2+{\overline{F^b}}^2+C_1}$$

$$c=\frac{2{\sigma }_a{\sigma }_b+C_2}{{\sigma }_a^2+{\sigma }_b^2+C_2}$$

$$s=\frac{{\sigma }_{ab}+C_3}{{\sigma }_a{\sigma }_b+C_3}$$

$$T(r)=T_{0}\cdot\bigg{(}\theta(R_{0}-r)+\theta(r-R_{0})\cdot\frac{R_{0}}{r}\bigg{)}$$

$$\theta (x)=\{\begin{array}{ll}1,& \text{for }x\ge 0\\ 0,& \text{for }x<0\end{array}$$

$$T(r)=T_{0}\cdot exp\big{(}-\frac{1}{2}\frac{r^{2}}{R^{2}_{0}}\big{)}$$

Full text

Infrared Imaging, Tomography, Image Reconstruction, AMIAS, RISE

I Introduction

Thermal Imaging provides a diagnostic imaging modality that can be used to visualize the temperature distribution in semi-opaque media (such as tissue and skin layer) from sets of infrared radiation measurements Ring and Ammer (2012); Baddour (2006). Utilizing non-ionizing radiation this imaging modality allows non-invasive measurements of temperature fluctuations in the human body caused by an abnormal generation of heat in the examined volume of interest. These fluctuations which are related to potential abnormalities of the cells’ metabolic activity can be correlated to other symptoms and medical imaging observations to indicate various pathogenies Lahiri et al. (2012).

Currently, the technique of thermal imaging is being explored and is considered potentially valuable for a number of medical applications including the monitoring of chemotherapy Keyserlingk et al. (2000), the diagnosis of inflammatory arthritis Ring and Collins (1970), osteoarthritis Glehr et al. (2011); Romano et al. (2011); Varju et al. (2004); Mayr (1995) and the assessment of peripheral circulation Vardasca (2012). Much effort has been extended for the detection of malignant diseases Herman and Cetingul (2011); Cruz et al. (2009); Maillard and Hessler (1969) and in particular for breast cancer Kirubha et al. (2015); Kennedy et al. (2009); Ng and Kee (2008); Lapayowker et al. (1973). Melanomas can be identified and localized within the skin layers of the patient by examining the different patterns of temperature distribution Ring and Ammer (2012). Various studies report that the heat generated due to the presence of cancerous lesions within the body can be up to ten times higher than that generated by healthy tissue Cetingul and Herman (2010); Song et al. (1984). This deviation is caused by the abnormally higher metabolic activity of the malignant cells Kirubha et al. (2015).

Thermal imaging, despite the promise, has not been yet accepted as a technique which can be used to provide accurate and reliable diagnostic information. The main limiting factor derives from the fact that the infrared radiation emitted from the inner organs is diffused, absorbed and scattered by the intermediate tissues before exiting the skin surface. As a result, the obtained infrared images from a source below a certain depth in the body are blurred and of low-resolution Snyder et al. (2000). The results mentioned above highlighting the promise of thermal imaging as a medical diagnostic tool refer to single planar imaging (thermography). The obvious next step for further development that of going from planar imaging to tomographic imaging has hardly been explored presumably because of the same reasons. Few successful attempts utilizing the traditional image reconstruction algorithms of emission tomography have been reported for 3D reconstruction of the flame temperature distribution Hossain et al. (2013); Goyal et al. (2014). In these studies, the tomographic images were produced in ”idealized” conditions using static infrared projections acquired in the absence of an intermediate absorber or scatterer. To date, the tomographic image reconstruction of the temperature distribution in an absorbing medium from the set of its infrared image projections remains a challenging problem with untapped potential.

II Objectives

The overarching objective of the work presented here is to demonstrate that infrared emission tomographic imaging is feasible and promising when implemented with the “Reconstructed Image from Simulations Ensemble (RISE)” technique. In particular this manuscript endeavors to:

1. I

Extend the well established and widely used methodology of emission tomography with ionizing radiation (Positron Emission Tomography (PET) or Single Photon Emission Computed Tomography (SPECT)) to Infrared Radiation Emission Tomography (IRET) utilizing infrared radiation. Enable the extension of IR imaging from the single thermal – planar imaging ( thermography) to a tomographic modality visualizing the 3D temperature distribution in an absorbing and scattering medium. 2. II

For the first time reconstruct tomographic images from hotspots embedded in a semitransparent medium. 3. III

Evaluate the success of implementing the proposed methodology employing the well established tomographic reconstruction techniques ART, MLEM, and RISE by utilizing software and hardware phantoms. Quantify the quality of the constructed tomographic images using well-accepted metrics and demonstrate that RISE is best suited for IRET.

The argumentation and the results supporting the achievement of the objectives as mentioned above are presented as follows: the fundamental aspects of the IRET problem are reviewed, and the mathematical foundations of IRET as an inverse tomographic problem are presented in Section III (”IRET Methodology”). The MLEM, the ART, and the RISE tomographic reconstruction methods are introduced in Section IV (”Image Reconstruction”). In Section V (”Simulation Studies”), software phantoms are utilized to benchmark the performance of the various reconstruction methods in IRET along with the appropriate metrics quantifying the quality of the reconstructed images. In Section VI (”Thermal Phantom”), an elaborate hardware phantom is employed to showcase IRET in achieving the stated objective. A summary of the presented results and an assessment of achieving the stated objectives along with suggestions for future work is presented in Section VII (”Conclusions”).

III InfraRed Emission Tomography (IRET) Methodology

The IRET problem concerns the reconstruction of tomographic images depicting the temperature distribution of the imaged structure derived from its infrared planar projections obtained at different angles.

IRET can be formulated as an inverse proble by assuming that the temperature distribution in the tomographic plane is represented by the 2D function $T(x,y)$ and that the intensity of the radiative power per unit wavelength ($\lambda$) emitted from a point $(x,y)$ can be calculated through the Planck’s formula Lienhard IV and Lienhard V (2008); Modest (2013):

[TABLE]

where $\epsilon(\lambda)$ is the wavelength depended emissivity of the imaged material and $c_{1}=3.7417\cdot 10^{-16}Wm^{2}$, $c_{2}=1.4387\cdot 10^{-2}mK$ are the first and second radiation constants respectively.

Utilizing a more convenient coordinate system $(t,s)$ :

[TABLE]

where $\theta$ is the angle of the camera rotation, the intensity of the radiative power per unit wavelength which is transmitted along the line path $L$ and being detected at $(t,\theta)$ by the detection unit is given by:

[TABLE]

where $k(\lambda,s,t)$ is the wavelength dependent attenuation coefficient of the imaged object. The geometry of the above problem is best explained with the help of Figure 1.

For an infrared detector (camera) operating in the spectral range ($\lambda_{1}$,  $\lambda_{2}$), the total radiative energy being detected is:

[TABLE]

where $\Delta\tau$ is the exposure time and $\Delta S$ is the area of the camera’s sensor.

For the infrared window ($\lambda_{1}$ = 7.5 $\mu m$,  $\lambda_{2}$ = 13.0 $\mu m$) and the temperature range (300 K - 340 K) examined in this work, the relation between the detected radiative energy $R_{0}$ and the temperature distribution $T(x,y)$ of the imaged object can be approximated by a linear function (Figure 2). We need to highlight that the linearization of the emitted radiative energy $R_{0}(t,\theta)$ as a function of temperature allows the formulation of IRET as a conventional emission tomographic problem. The IRET problem can then be formulated as a system of equations described by a projection matrix. In such a formulation of the infrared tomographic problem, a bin measurement $R_{i}$ corresponding to the total energy transferred by a finite width beam can be estimated from the sum:

[TABLE]

where, $F_{j}$ is the vector representation of the tomographic image of the temperature distribution, $N^{2}$ is the number of image elements, and $P_{ij}$ is the projection matrix weighting the contribution in transmitted photons of the $j^{th}$ image element to the $i^{th}$ bin (pixelated) measurement.

In emission tomography Lalush and Wernick (2004); Natterer and Wübbeling (2001), the Projection Matrix $P_{ij}$ not only provides the geometric weight with which the $F_{j}$ element (pixel/voxel) contributes to $R_{i}$, it also incorporates all the attenuation and scattering effects that the emitted photons undergo along their path to the detector. Unlike the x-ray or $\gamma$-ray regions, in the infrared region, the dominant light-tissue interaction is elastic scattering, and the ability to accurately incorporate its effects is important. The incorporation of attenuation and scattering into $P_{ij}$ may entail detailed and cumbersome cascade calculations. In this generalized approach, the association of $R(\lambda,t,\theta)$ to $I(T,\lambda,s,t)$ via the Projection Matrix can no longer be interpreted in terms of simple geometrical optics as in Figure 1, but rather as the response of the detector element to the source of radiation regardless of how the photon reached it (directly or via scattering). Methods for modeling light propagation in tissue  Ishimaru (1978); Bohren and Huffman (2008); Nissilä et al. (2005) are well developed, and the methodology above allows their implementation for specific cases. In real samples (e.g. in infrared tomographic medical imaging), knowledge about the morphology (structure) and the optical properties (attenuation and scattering cross sections) will be taken into account for the construction of the Projection Matrix.

In the examples presented in this methodological paper, we examine cases where uniform attenuation suffices. In such cases,

[TABLE]

where $\tilde{P}_{ij}$ is the ”bare” Projection Matrix which takes into account only the geometrical weight of each pixel/voxel, $k_{u}$ is a uniform attenuation coefficient characterizing the imaged medium, and, $r_{ij}$ is the distance of the $j^{th}$ image element from the boundaries of the imaged volume, measured in the direction towards the $i^{th}$ pixelated detector of the thermal camera system.

The simplified uniform attenuation map employed in this work provides a good approximation of the attenuation effects characterizing a homogeneous medium; such is the case of the hardware phantom examined in this study. Different and more complex attenuation and scattering models can be employed if knowledge on the morphology of the medium is available, as it is expected to be the case in applying IRET to medical imaging.

IV Image Reconstruction

IV.1 Iterative Reconstruction Techniques

Two widely used reconstruction techniques are employed in this study to solve the image reconstruction problem defined in Equation 6 and provide reference images for comparison are the Algebraic Reconstruction Technique (ART) and the Maximum Likelihood Expectation Maximization (MLEM) method. ART Gordon et al. (1970); Gilbert (1972) and MLEM Shepp and Vardi (1982); Lange et al. (1984) are well established and widely used methods in emission tomography, although MLEM invariably yields superior results and it is widely regarded as representing the ”State of the Art” method. Both methods are capable of incorporating the projection model of Equation 6 and can provide reconstruction results using sets containing a limited number of planar projections Wieczorek (2010). In the present work, the Newton-Raphson variant of ART Angeli and Stiliaris (2009) and an accelerated version of MLEM Hwang and Zeng (2006) were used for the reconstruction of the tomographic images.

IV.2 The Reconstructed Image from Simulations Ensemble (RISE)

RISE provides an alternative method to ART and MLEM for tomographic image reconstruction. Based on the AMIAS framework Alexandrou et al. (2015); L. Markou et al. (2018); Papanicolas and Stiliaris (2012), RISE produces images by parametrizing the physical characteristics of the imaged ”target” and simulating a large ensemble of image configurations, each one corresponding to a solution of the tomographic problem. In the case of IRET, the term ”target” is used to represent a source of thermal radiation. In AMIAS and RISE, each of the numerous randomly generated solutions describing the characteristics of the target is assigned a weight quantifying its ”goodness” of representing the projection data. The image is reconstructed by statistically weighting the entire ensemble of the simulated solutions. The formulation of RISE and the mathematical concepts of AMIAS are presented elsewhere  Alexandrou et al. (2015); L. Markou et al. (2018); Papanicolas and Stiliaris (2012); Papanicolas, C. N. and Koutsantonis, L. and Stiliaris, E. (2018). In this work, we summarize the essential features of the method and the specific choices required to implement it in the IRET problem; they are described in the next paragraphs.

IV.2.1 Modeling the Imaged Object

A model parametrized by a set of parameters is chosen to represent the physical characteristics of the target. Prior knowledge can be incorporated at this step to endow the scheme with Bayesian capabilities. The model employed in RISE to describe the temperature distribution is that of a summation of elementary sources (shown in Figure 3) having an ellipsoidal shape and embedded in a slowly varying background:

[TABLE]

where:

• •

$Z^{j}_{i}(x,y)$ is a set of Zernike polynomials Lencrerot et al. (2009); Born and Wolf (1959) representing the background tempereture distribution and $C^{j}_{i}$ are their associated amplitudes,

• •

$M_{z}$ defines the total number of such polynomial required to approximate the background distribution,

• •

$N_{s}$ is the total number of elementary ellipsoidal targets,

• •

$A_{i}$ is the amplitude of temperature at the center of the $i^{th}$ target,

• •

$r_{i}$ is the euclidean distance of a point $(x,y)$ lying on the tomographic plane from the center $(x_{i},y_{i})$ of the $i^{th}$ target,

• •

$s_{i}$ is a coefficient defining the ”sharpness” of the temperature distribution in the surrounding by the source medium, and,

• •

$r_{0}(\theta)$ is a geometrical factor given as a function of the semi-major $u$ and semi-minor $v$ axes of an ellipse respectively:

[TABLE]

where the angle $\phi$ is the above parametric equation defines the orientation of the target in the tomographic plane.

IV.2.2 Composing the Ensemble of Solutions

A Monte-Carlo procedure is used to sample the parameters of interest. Each set of sampled parameter values constituting a solution to the tomographic problem is used in Equation 8 to construct a tomographic image configuration. A set of projections is simulated from the tomographic image by employing the forward projection model presented in Equation 6. For each solution, the model estimated projections are compared to the measured projections through the $\chi^{2}$ criterion:

[TABLE]

where $\epsilon_{i}$ is the associated error of the $i^{th}$ projection measurement. The $\chi^{2}$ value is assigned to the solution to quantify its goodness of representing the data. The procedure is repeated to construct a large ensemble of possible solutions (typically 105 to 106), which as prescribed by the RISE method, is statistically weighted to determine the parameters of interest.

One of the distinct features of RISE and the underlying AMIAS, which sets it apart from the standard reconstruction methods, is that it can incorporate any forward model simulating the propagation of the radiative energy from the source to the detector. This forward model can be a function or a functional describing a very complex propagation process such as the entire cascade of non-uniform absorption and re-scattering.

IV.2.3 Determining the parameters of interest

The RISE reconstruction process allows through different techniques Alexandrou et al. (2015); L. Markou et al. (2018); Papanicolas, C. N. and Koutsantonis, L. and Stiliaris, E. (2018) the determination of the total number of terms ($N_{s}$) which are needed in the model (Equation 8) to represent the thermal sources. The number $N^{opt}_{s}$ best describing the data is selected as the one minimizing the Bayesian Information Criterion (BIC)Neath and Cavanaugh (2012) given by the formula:

[TABLE]

where the $\chi^{2}_{min}(N_{s})$ is the minimum $\chi^{2}$ value in the ensemble of solutions constructed for the model of $N_{s}$ terms and $M$ is the number of ray projections used for the reconstruction of the tomographic image.

Having defined the number $N^{opt}_{s}$, the solutions in the corresponding ensemble are weighted with the probability value $exp(-\frac{1}{2}\chi^{2})$ to derive the Probability Distribution Functions (PDFs) of the parameters $(A_{i},x_{i},y_{i},u_{i},v_{i},\phi_{i},s_{i})$. Mean values and associated uncertainties are derived from the PDFs and used in Equation 8 to reconstruct the tomographic image of the temperature distribution.

V Simulation Studies

The validation of IRET and in particular the success of RISE in IRET has been accomplished with simulated data from numerical phantoms and real data from an experimental study with a hardware phantom.

In three case studies employing five software phantoms, the ability of ART, MLEM, and RISE to reconstruct the image of the ”true” temperature distribution has been examined using simulated sets of noisy infrared projections.

V.1 Image Quality Metrics

The objective assessment of a reconstruction method requires the calculation of different metrics quantifying the quality of the produced images. For the simulation studies presented, the metrics employed for this task are the Correlation Coefficient (CC), the Normalized Mean Square Error (NMSE), the Peak Signal-to-Noise Ratio (PSNR), the Contrast-to-Noise Ratio (CNR) and the Structural SIMilarity (SSIM) index.

CC, which provides a spatial similarity measure between two images Hill et al. (2001), is given by:

[TABLE]

where $N^{2}$ is the total number of image pixels, $F^{a}_{i}$ and $F^{b}_{i}$ are the pixel values of the reconstructed and ”true” image respectively, and $\bar{F^{a}}$, $\bar{F^{b}}$ are their corresponding average pixel values. The CC as a metric expressing the spatial similarity between two images it was used to evaluate the capability of ART, MLEM and RISE to adequately resolve the geometrical characteristics of the imaged ”targets”.

The NMSE and PSNR metrics were used to calculate the overall reconstruction error and to provide, respectively, a global measure of image contrast. The two metrics are defined as:

[TABLE]

[TABLE]

CNR is meant to objectively evaluate the detectability of a ”target” in a noisy background. It provides a local measure of contrast; it was calculated on a Region of Interest (ROI) extracted from the reconstructed image through the formula:

[TABLE]

where $\bar{T}$ is the average reconstructed temperature value in the target ROI, $\bar{B}$ and $\sigma_{B}$ is the average and the standard deviation, respectively, of image elements corresponding to the background ROI.

The SSIM index is calculated on the reconstructed images with respect to the ”true” phantom images to quantify their visual similarity. SSIM is defined as a multiplicative combination of three indexes, namely the luminance index ($l$), the contrast index ($c$) and the similarity index ($s$) Wang et al. (2004):

[TABLE]

where the indexes $l$, $c$ and $s$ are defined in terms of the average pixel values $\bar{F^{a}}$ and $\bar{F^{b}}$, the standard deviations $\sigma_{a}$ and $\sigma_{b}$, and the covariance $\sigma_{ab}$ of the reconstructed ($a$) and ”true” ($b$) images:

[TABLE]

[TABLE]

[TABLE]

The coefficients $C_{1}$, $C_{2}$ and $C_{3}$ are constants introduced to prevent numerical instabilities when the denominators in the above equations are close to zero Wang et al. (2004). Unlike the CC index evaluating the spatial characteristics of an image, the SSIM index accounts for the additional comparisons of pixel intensities to quantify the similarity in luminance and contrast between the reconstructed and ”true” image.

V.2 Numerical Phantoms and Projection Data Generation

Five numerical phantoms, shown in Figure 4, were used to generate the sets of infrared projections. The phantoms simulating the presence of thermal sources (”hotspots”) were defined on rectangular grids of $64\times 64$ pixels.

In the first phantom (”A”), the two depicted sources were generated using the radial temperature distribution:

[TABLE]

where $T_{0}$ is a constant coefficient representing the temperature within the area of each hotspot, $R_{0}$ is the radius of the hotspot, and:

[TABLE]

Phantom ”B” has the same geometry as phantom ”A”; however, the radial temperature distribution in the medium surrounding each source exhibits a Gaussian profile, given by the equation:

[TABLE]

where, in this radial profile, the constant coefficient $T_{0}$ represents the temperature at the center of the hotspot and $R_{0}$ is a constant defining the diffusion of temperature in the surrounding medium. For both phantoms (”A” and ”B”), the coefficient $R_{0}$ was set to $2$ pixel units.

Phantom ”C” is a variant of the Shepp-Logan mathematical phantom; it consists five ”hotspots” of different size and orientation having a Gaussian radial profile (see Equation 21). The parameters values defining the five hotspots are given in Table 1.

Phantom ”D” consists the same five hotspots of phantom ”C” embedded in a disk having a non-uniform and non-symmetric temperature distribution. The temperature distribution of the disk introducing physical effects of background into the projection data was defined as a second order polynomial of the $(x,y)$ coordinates. The maximum temperature of the background was set to 2.2*∘*C. Both ”C” and ”D” phantoms present challenging cases having five hotspots of different size, temperature, and orientation. Given that the projection data were generated by simulating attenuation conditions, the reconstruction of all hotspots and especially of the one in the middle, lying between the two larger ones, is particularly difficult.

Phantom ”E” was also generated to present a non-uniform - non-symmetric temperature distribution. It presents two small hotspots, each one having a Gaussian radial profile (Equation 21) and a radius of 2 pixel units. The temperature $T_{0}$ of the two hotspots was set to 2*∘C, whereas the maximum temperature of the background was equal to 2.2∘*C.

Sets of 24 infrared projections were generated from each phantom using the forward projection model of Equation 6. For each simulation case, the set of 24 projections was generated in the full $360^{\circ}$ angular range with a constant angular step of $15^{\circ}$. The detected infrared radiation $R_{i}$ resulting from the true temperature distribution $T_{j}$ was calculated using the linearization of the problem as provided in Equation 6 for the spectral range $[7.5\mu m,13\mu m]$ assuming that $\epsilon(\lambda)=1$. The uniform attenuation coefficient of the medium $k_{u}$ was set to the value of $0.1/w_{p}$ where $w_{p}$ is the size of a pixel. Each projection was generated as a vector of $91$ rays $(NR=91)$ and further randomized with Gaussian distributed noise $(n\sim\mathcal{N}(0,0.1^{2}))$.

The simulated sets of projections were used in three distinct studies examining ART, MLEM and especially RISE in IRET:

1. A

Model Capacity: The model formulated in RISE (Equation 8) is examined in its capability to provide sufficient reconstruction results when the ”true” temperature distribution has a general and not necessarily the same radial temperature profile. Both of the two phantoms (”A”, ”B”) defined by setting the temperature coefficient $T_{0}$ to $4^{o}C$ were used in this study. 2. B

Resolving of Temperature Differences: The efficacy of ART, MLEM, and RISE to reconstruct images revealing small differences in the temperature distribution is examined by the use of the phantom ”A”. In three simulation cases, the temperature coefficient $T_{0}$ was varied from $1^{o}C$ to $3^{o}C$ with a step size of $1^{o}C$ to produce different realizations of the temperature distribution. 3. C

Non-uniform Background Distribution: This study was conducted to examine the capability of the three reconstruction methods to identify hotspots in the strong background. The effects of the non-uniform - non-symmetric background distribution on the reconstruction quality were examined using the numerical phantoms ”C” and ”D”. The phantom ”E”, also presenting a non-uniform - non-symmetric background distribution was used to assess the detectability of small-sized targets in the produced reconstructions.

V.3 Reconstruction Results

V.3.1 Model Capacity

Figure 5 shows the reconstructed images of the two simulated phantoms as they were produced in the first simulation study with RISE, ART and MLEM. ART and MLEM reconstructions were obtained by performing two and three cycles of iterations respectively.

From a visual inspection of the reconstructed images, it can be seen that RISE images exhibit higher contrast and less noise as compared to ART and MLEM reconstructions. Line-profiles extracted from the reconstructed images across the diagonal line connecting the centers of the two hotspots are depicted in the same figure. As seen, the Fermi-like model employed in RISE yielded images that well-resolve the two simulated temperature distributions.

The image quality metrics (CC, NMSE, PSNR, CNR, and SSIM) comparing the reconstructed images are presented in Table 2. The results confirm and quantify the visual observations. For both phantoms, while all the methods lead to acceptable image quality metrics, RISE led to superior CC and SSIM values (the highest among the three methods) indicating its ability to resolve the image of the temperature distribution adequately. Moreover, the RISE images show an improved contrast and hotspots detectability as indicated by the calculated PSNR and CNR scores.

V.3.2 Resolving of Temperature Differences

RISE, ART and MLEM were used to reconstruct the images of the phantom ”A” for three simulation cases varying the temperature coefficient $T_{0}$ in Equation 20 ($T_{0}=1^{o}C$, $T_{0}=2^{o}C$, $T_{0}=3^{o}C$). RISE reconstructions were produced using the Fermi-like model, the appropriateness of which is validated in the previous simulation study. For the case of ART and MLEM, the reconstructions were obtained by performing two and three grand iterations respectively.

The reconstructed images from the three methods are shown in Figure 6. In all three simulation cases, streak artifacts reducing the detectability of the target are presented in ART reconstructions. RISE and MLEM images exhibit less amount of noise and provide reconstructions better revealing the simulated distribution of temperature. By visually examining these images it can be seen that the boundaries of the two hotspots are well shown in RISE and MLEM images, whereas, it is hard to detect them in the image produced with ART.

The CC, NMSE, PSNR, CNR and SSIM scores evaluating the quality of the reconstructed images by the three methods are shown in Table 3. The calculated scores indicate superior image quality for the images reconstructed with RISE. As expected, the reconstruction quality of the three methods is degraded as the temperature difference $T_{0}$ in the simulated phantom is decreased. For the case exhibiting the smallest difference in temperature ($T_{0}=1^{o}C$), the case with the most interest in this simulation study, the highest CC and SSIM values obtained for RISE show an improvement in image similarity as compared to the corresponding CC and SSIM scores of MLEM and ART images respectively. For the same simulation case, the RISE reconstruction is shown to be superior concerning the hotspots detectability (CNR). RISE led to a CNR which is $192\%$ and $100\%$ higher than those calculated for the images reconstructed with ART and MLEM respectively. Likewise, the NMSE value obtained for RISE is more than two times less than those obtained for ART and MLEM.

V.3.3 Non-uniform Background Distribution

Figure 7 compares the reconstructed images of the phantoms ”C”, ”D”, and ”E” obtained with RISE, ART and MLEM in the third simulation study. The calculated metrics (CC, NMSE, PSNR, CNR, SSIM) comparing the reconstructed images of the three methods are shown in Table 4.

In the absence of background (phantom ”C”), as indicated by the CC and SSIM scores, the RISE image shows the highest spatial and structural similarity with the true image. In this ”ideal” case of zero background, RISE presents improved image contrast and hotspots detectability as quantified via the PSNR and CNR respectively. From Figure 7, it can be seen that the reconstruction quality of the three methods is affected by the introduction of the non-uniform background in phantom ”D”. As the image quality metrics also indicate it in Table 4, both the structural similarity (SSIM) and the hotspots detectability (CNR) are reduced for all three methods. However, as compared to ART and MLEM, the scores obtained for RISE show a superior performance which leads to an acceptable hotspots detectability (CNR) and a sufficiently high structural similarity (SSIM). Visually, the four of the five hotspots of ”Phantom D” can be identified in the RISE image, whereas it is hard to be separated from the background in the images produced with ART and MLEM.

Differences in image characteristics between the three methods are also apparent in the reconstructions of phantom ”E”. The combination of the non-uniform background with the small size of the hotspots simulated in this case makes the reconstruction problem difficult. As seen in Figure 7, the two hotspots cannot be observed in ART and MLEM reconstructions. The RISE image clearly shows the boundaries of the two hotspots which can be easily separated from the background. From Table 4, it can be seen that RISE yielded the highest spatial (CC) and structural (SSIM) similarity scores while it also led to the highest contrast (PSNR) and hotspots detectability (CNR).

Overall, the results from the three simulation studies validate the reconstruction efficacy of RISE and indicate good performance in cases of non-uniform - non-symmetric background distributions. The model which is employed in the method (Equation 8) to represent the imaged targets proves to be adequate to describe hotspots exhibiting different radial intensity profiles. This feature of the method proves to be beneficial in reducing the amount of noise and improving the hotspots detectability of the reconstructed image.

VI Experiments with a Thermal Phantom

The evaluation of RISE in IRET was further explored using sets of thermal images captured from a thermal phantom. In two experimental studies reported in Rapsomanikis et al. (2012), sets of 24 infrared images were acquired from the phantom using an infrared camera. In the present study, RISE and the two conventional methods ART and MLEM are examined in their capacity to reconstruct the tomographic images of the phantom temperature distribution. The derived tomographs by the three methods were assembled in stacks to visualize the 3D distribution of temperature in the volume of interest.

VI.1 The Thermal Phantom

The thermal phantom was constructed out of three pairs of $R_{A}=100\Omega$ and $R_{B}=47\Omega$ resistors, alternately configured to form the capital letter ’M’, with the two middle resistors oriented out of the plane (Figure 8). The resistors were connected in series and supplied with direct current. A thermocouple monitoring the temperature of the phantom was adjusted at the upper-left corner of the configuration as it is shown in Figure 8.

Two experimental studies Rapsomanikis et al. (2012) were carried out with the implemented hardware phantom:

1. A

Zero Absorption Case. The thermal radiation was captured as it was emitted directly from the resistors with no intervening medium. The temperature of the phantom measured by a thermocouple was comparable to the core temperature of humans ($38^{o}C$), and the environmental temperature measured at a distance of 5 cm from the phantom was $25^{o}C$. The set of infrared projections obtained in this study was used to provide reconstructed images for visual comparisons. 2. B

Absorption Case. The same configuration of six resistors was placed in a conical vessel filled with agarose gel of about $1$ gr/ml concentration. The inner diameter of the vessel was $40$ mm (small diameter) at its top surface and $65$ mm at its bottom surface (big diameter). The environmental temperature was $25^{o}C$, and so was the temperature of the upper-left resistor, kept at $38^{o}C$. In this study, we examine the ability of RISE to reconstruct tomographic images in absorption conditions degrading the quality of the planar thermal images.

VI.2 Data Acquisition

In both experiments, the phantom was placed on a rotating table at a distance of $15$ cm from the front of the camera. A set of $24$ infrared images of the phantom, captured using a thermal camera (Thermovision 570, AGEMA Infrared Systems), were obtained in the full range of $360^{\circ}$ with an angular step of $15^{\circ}$. The thermal camera has a $24^{o}\times 18^{o}$ Field of View (FOV) and is characterized by $0.1^{o}C$ thermal sensitivity; it operates in the spectral range of 7.5 to 13 $\mu m$. The 24 planar images (thermographs) were further sliced to provide a set of 15 sinograms, each one corresponding to a tomographic level at a specific vertical offset. The extracted sinograms were used as inputs in RISE, ART, and MLEM to reconstruct the sets of $64\times 64$ tomographic images of the phantom.

VI.3 Reconstruction Results

The reconstructed images obtained in the two experimental studies, with and without absorption conditions, are visually compared in Figures 9 and 10. The capability of each method to reproduce the structure of the imaged ’M’ shaped object and to visualize the difference in temperature distribution from the alternatively positioned resistors were the two primary criteria for the evaluation of the images.

All the three methods yielded reliable reconstructions in the absorption free case (Figure 9) validating the IRET methodology as described in Section III. From a visual inspection of the tomographic images reconstructed by the three methods (shown at the bottom panel of Figure 9), it can be seen that the alternation in temperature between the $100\Omega$ and $47\Omega$ resistors is visualized in the RISE image. This alternation in temperature cannot be observed in the images reconstructed with ART or MLEM.

The RISE, ART and MLEM 3D volumetric images of the phantom (shown at the top panel of Figure 9) obtained by assembling the 2D tomographic images in stacks indicate that all three methods successfully resolved the geometry of the phantom. In the case of RISE, the structure of the phantom was visualized by surfacing the volumetric data at the determined model radius parameter $r_{0}$ (see Equation 9). In the case of ART and MLEM, isothermal surfaces were extracted at half of the maximum temperature to resolve the 3D shape of the thermal phantom.

The tomographic images of the phantom produced in the second experimental study (”Absorption Case”) with the three methods are shown in Figure 10. For all three reconstruction methods, images were obtained by using a uniform attenuation coefficient $\kappa_{u}=0.23mm^{-1}$ (Equation 6). From a visual inspection of the images produced in this case study, it can be observed that the results obtained with RISE are superior. ART images exhibit noisy artifacts and provide limited information about the structure of the phantom. Compared to MLEM, the RISE images exhibit higher contrast and provide better approximations of the thermal sources. In the images produced by the three methods, the difference in temperature from the alternatively positioned resistors is not observed.

The 3D images of the phantom composed from the sets of the reconstructed images by following the same procedure used in the previous case (”Zero Absorption Case”) are shown at the top panel of Figure 10. RISE yields more reliable imaging of the ’M’ shaped phantom and to a lesser degree MLEM, while ART fails. The amplitude of temperature at the center of the phantom, indicating the presence of the middle pair of resistors, is more visible in the RISE result.

An additional 3D image of the phantom visualizing the RISE reconstruction is shown in Figure 11. This volumetric image depicts the distribution of temperature on the adjacent surface lying at a distance $r_{0}$ from the center of sources. The 3D reconstruction reveals a positive gradient in temperature which is observed across the vertical axis ($z$-axis). This gradient is caused by the difference in the produced thermal energy between the upper and bottom part of the phantom, comprising four and two heaters respectively. It can also be observed that compared to the RISE image obtained in the ”Zero Absorption Case,” the image obtained in this study with the same method presents an over-sized representation of the thermal sources. This difference can be understood as resulting from the heated medium (agarose gel) surrounding the resistors. A more precise image of the heaters, if such a result is desired, would require the implementation of a more sophisticated model that takes into account heat transport effects. Such a forward model can be implemented in the general framework of RISE to accommodate the diffusion of infrared radiation within the absorbing material.

The tomographic images of the hardware phantom produced by RISE compared to those of MLEM and ART, demonstrate the suitability of RISE for Infrared Tomography. These images should also be compared with the planar thermal images (thermographs) of the phantom (shown in panel D of Figure 8), which do not allow any conclusions to be drawn for the thermal distribution within the gel.

VII Conclusions

In the work presented in this paper, the Infrared Emission Tomography (IRET) technique is revisited, and it is shown that the use of the methodology employed in ionizing radiation tomography (e.g., PET or SPECT modalities) is both justified and implementable. The IRET methodology implemented with the widely used MLEM and ART techniques yields satisfactory results, especially in cases where attenuation and scattering effects are minimal. The RISE technique yields superior tomographic images even in cases where medium modifications are present.

Images from all three reconstruction methods were evaluated using well-established metrics of the field. All yielded acceptable tomographic results according to these metrics although RISE images exhibit a significantly higher structural similarity (SSIM) to the simulated distributions and improved hotspot detectability (CNR). The three methods were also used to reconstruct tomographic images from experimental data obtained from a thermal phantom imaged in nil and high absorption conditions. In the case of absorption, which represents a semi-realistic case simulating attenuation conditions in the human body, RISE yielded an image that reliably resolves the geometry of the phantom and visualizes the gradient of the temperature distribution within the image medium.

The robust methodology presented in this work specifies how further refinements can be accomplished with the inclusion of detailed modeling of non-uniform attenuation and scattering effects. Following the reported successful implementation of the RISE method in IRET further expanded experimentation with phantoms and possibly small animals is warranted to ascertain the potential of this modality in medical imaging.

Acknowledgements.

This work was supported by the Graduate School of The Cyprus Institute and the Cy-Tera Project ”NEA IPODOMI/STRATI/0308/31”, which is co-funded by the European Regional Development Fund and the Republic of Cyprus through the Research Promotion Foundation.

Conflicts of interest

We have no conflict of interest to declare.