The use of transport and diffusion equations in the three-dimensional reconstruction of computerized tomographic images

O uso de equações de transporte e difusão na reconstrução de imagens de tomografia computadorizada em três dimensões

ARTIGO ORIGINAL

The use of transport and diffusion equations in the three-dimensional reconstruction of computerized tomographic images

O uso de equações de transporte e difusão na reconstrução de imagens de tomografia computadorizada em três dimensões

Sandrerley Ramos PiresI,^*; Edna Lúcia Flôres^III; Celia Aparecida Zorzo Barcelos^II; Dulcinéia Gonçalves Ferreira Pires^III; Gilberto Arantes Carrijo^III; Antônio Cláudio Paschoarelli Veiga^III

^IEscola de Engenharia Elétrica e de Computação – EEEC, Universidade Federal de Goiás – UFG, Av. Universitária, n. 1488, Setor Leste Universitário, CEP 74605-010, Goiânia, GO, Brasil

^IIFaculdade de Matemática, Universidade Federal de Uberlândia – UFU, Av. João Naves de Ávila, 2121, Campus Santa Mônica, CEP 38048-100, Uberlândia, MG, Brasil

^IIIFaculdade de Engenharia Elétrica – FEELT, Universidade Federal de Uberlândia – UFU, Av. João Naves de Ávila, 2160, Campus Santa Mônica, CEP 38400-902, Uberlândia, MG, Brasil

ABSTRACT

The visualization of a computerized tomographic (TC) exam in 3D increases the quality of the medical diagnosis and, consequently, the success probability in the treatment. To obtain a high quality image it is necessary to obtain slices which are close to one another. Motivated towards the goal of reaching an improved balance between quantity of slices and visualization quality, this research work presents a digital inpainting technique of 3D interpolation for CT slices used in the visualization of human body structures. The inpainting is carried out via non-linear partial differential equations (PDE). The PDE's have been used, in the image-processing context to fill in the damaged regions in a digital 2D image. Inspired by this idea, this article proposes an interpolation method for the filling in of the empty regions between the CT slices. To do it, considering the high similarity between two consecutive real slice, the first step of the proposed method is to create the virtual slices. The virtual slices contain all similarity between the intercaleted slices and, when there aren't similarities between real slices, the virtual slices will contain indefinite portions. In the second step of the proposed method, the created virtual slices will be used together with the real slices images, in the reconstruction of the structure in three dimensions, mapped onto the exam. The proposed method is capable of reconstructing the curvatures of the patient's internal structures without using slices that are close to one another. The experiments carried out show the proposed method's efficiency.

Keywords: 3D image interpolation, Computerized tomography, Differential equations.

RESUMO

A visualização do resultado de uma tomografia computadorizada (CT) em 3D aumenta a qualidade do diagnóstico médico e, consequentemente, a probabilidade de sucesso no tratamento. Para obter uma maior qualidade na imagem 3D, é necessário obter fatias próximas umas das outras. Motivados pela meta de obter um apropriado balanceamento entre a quantidade de fatias e a qualidade da visualização, este trabalho apresenta uma técnica de retoque digital por interpolação 3D para fatias de CT usadas na visualização de estruturas do corpo humano. O retoque digital é feito através de equações diferenciais parciais não lineares (EDP). As EDPs têm sido utilizadas no contexto de processamento de imagens para preencher regiões defeituosas de imagens em 2D. Inspirado por esta ideia, este artigo propõe um método de interpolação para preencher as regiões vazias entre as fatias de TC. Para isto, considera-se a grande semelhança entre duas fatias reais consecutivas, onde o primeiro passo do método proposto é a criação de fatias virtuais. As fatias virtuais contêm todas as regiões similares entre as duas fatias intercaladas, sendo que onde não há semelhança entre as fatias reais, as fatias virtuais conterão regiões indefinidas. No segundo passo do método proposto, as fatias virtuais criadas serão utilizadas junto com as fatias reais para a reconstrução da estrutura em três dimensões, mapeadas no exame. O método é capaz de reconstruir as curvaturas das estruturas internas do paciente sem a necessidade de utilizar fatias muito próximas uma das outras. Os experimentos realizados mostram a eficiência do método proposto.

Palavras-chave: Interpolação de imagens em 3D, Tomografia computadorizada, Equações diferenciais.

Introduction

The visualization of a computerized tomography (CT) exam in 3D increases the quality of the medical diagnosis and thus the efficiency of the treatment. This type of visualization relieves the medical professional from the task of mentally reconstructing the human body's internal structures, thus aiding him/her in arriving at a more accurate diagnostic analysis (Hernandez et al., 2006; Piccand et al., 2008).

Traditionally, the CT results show images that are transversal cuts (slices) from the analyzed structure. In an exam these images are distanced one from the other by a few millimeters, in general between 1 to 10 mm. Nowadays, using a CT scanner with a spiral approach, it is possible to provide the 3D visualization. In this case, slices are only subsets of the existing data; however, normally the examination results are still slices. Thus, the 3D visualization is achieved through the placing of the slices side by side and the filling in of the space between them, through the use of an interpolation process.

Even though, many software packages (E-Film Solutions, 2011; 3D Slicer, 2010; Nolf, 2003) using a variety of techniques for image visualization are available, it is generally appreciated that to present a high quality image, it is necessary to obtain slices which are close to one another. To illustrate this, we present in Figure 1 the reconstruction of a human jaw, through the use of an interpolation process using different quantities of slices to visualize the same structure.

Though there are tomographs able to produce slices very close to each other (1 mm or less), they are not common in daily use. Another possibility to reduce the processing time is to reduce the amount of slices (Genu, 2007). Thus, it is still common to find computerized tomography exams where the distance between two consecutive slices is more than five times the distance between two adjacent pixels in the 2D image. Such low resolutions demand the use of better interpolation methods, rather than a simple linear interpolation, as proposed for Tbaina et al. (1997). According to Wang et al. (2001) the linear interpolation produces a weak approximation of the real artifacts. The produced result is similar to a "jagged steps", as shown in Figure 1b. Despite Figure 1a also presents the same effect, it is possible to observe a better quality in the image, because the "jagged steps" effects is reduced, when compared to the second case.

Motivated by the goal of reaching an improved balance between slice's quantity and visualization quality, this work presents digital inpainting as a form of 3D interpolation for CT slices used in the visualization of structures from the human body. This inpainting is carried out via non-linear partial differential equations, aiming to produce much smoother and better defined forms for the structure in question, without the need of a large quantity of slices. Initially, digital inpainting was used for the recovery of photographic images, as well as old films, as described by Bertalmio et al. (2000) which uses non-linear PDEs for the recovery of 2D images. This research uses a similar approach, but for 3D interpolation.

Correlated works

This subsection presents some previous research work on 3D image interpolation, with added comments showing how the cited research contributed to the method proposed in this article.

Goshtasby et al. (1992) proposed a linear interpolation method between two slices, where a similarity measure defined for a pixel in a certain slice will be the corresponding pixel in the other slice. For this to be completed adequately, a comparison is made between the reference pixels in one slice to all the neighboring of corresponding pixels in another. This connection between the pixels defines the path for the linear interpolation. The same idea is used in our method. However, a tolerance threshold is used in the comparison process between the real pixels, it is possible that there exist pixels, which hold no correspondence to another, thus creating regions with no information.

Wang et al. (2001) proposed an approach where the free spaces between the real slices were divided into small pixel cubes. From the analysis of the alterations in the gray level on the cube's corners, the direction of the interpolation process is defined. The obtained results are similar to those obtained by linear interpolation. The authors also explore this idea. The pixel neighborhoods are used to define the connections between the slices, thus increase the space for the identification of similarities.

Shuying et al. (2005) proposed a non linear interpolation method, where the pixel intensity of a virtual slice is determined by a weighed average derived from the pixels of the adjacent real slices. However, the definition of the weights is reached through the use of a non-linear function such as a cosine function. That work is suported by the hypothesis that the change of the internal organ's form is smooth and non linear. This idea is also explored in this research work.

Another approach for the generation of 3D images with spaced apart slices was presented by Köstler et al. (2006). The interpolation process is carried out on the sinogram image, before the reconstruction of the original image. The best result obtained by Köstler was through the use of anisotropic diffusion, which results in a much smoother encounter of the different flows represented on the sinogram. The anisotropic diffusion is used in both stages of the interpolation process proposed in this research work, with the goal of reaching a greater preservation of the image contours.

The methods of Bors et al. (2002), Goshtasby et al. (1992) and Wang et al. (2001) consider the hypothesis that there exist no abrupt changes along the contours of human body tissue. In this way, they consider the possibility of the continuity of a determined tissue from a real slice onto the adjacente slice. This is the fundamental hypothesis to the method proposed in this article.

There also exist other approaches for carrying out 3D reconstruction. Among them we can mention the use of mesh as in Cateriano and Nonato et al. (2005) and of Taillandier and Deriche (2002) or even the use of pixel cubes as in Harald and Scherzer (2003).

Materials and Methods

With the aim of building 3D images with a small quantity of CT slices, which are denoted here as real slices, the authors propose to create virtual slices between two real slices. The creation of these virtual slices is carried out in two phases. The first phase consists of creating an initial representation of the virtual slices by analyzing the similarities. This can be done, because the images of the internal human body structures represented on CT slices present great similarity to one another (Goshtasby et al. 1992), as shown in Figure 2.

With the initial representation complete, the second phase consists of updating the existing information in this representation by transporting the 3D information to the regions, that do not receive information during the first phase. The non-similar regions which still have no information at the end of the first phase, will define the inpainting domain D in the second phase, according to the proposed method.

Initial representation of the 3D image

The joining of real and virtual slices composes the initial representation of the 3D image. Between the real slices it is possible to add one or more copies of the virtual slices. The initial quantity of inserted virtual slices between two adjacent real slices is obtained through the distance and the resolution of these images. All the virtual slices between two adjacent real slices will contain the same information after the first phase.

To create a virtual slice between two adjacent real slices, F1 and F2, the similarity between them needs to be analyzed. The initial representation of virtual slices is composed by similar portions between F1 and F2. Non similar slices portions, denoted as A* in Equation 1, will be created in second phase of the process.

This similarity is determined in the following way: let be F1 and F2 two adjacent real slices and Fv a virtual slice, which will be created between them. A particular pixel from Fv, denoted as Fv (i, j), is defined by:

where:

and

• T_g - tolerance threshold;

• A* - attributed value given simply for the definition of the pixels which are going to define the inpainting domain D.

To determine the value of the threshold T_g, the analysis of the real slice's gray level standard deviation σF_n involved in the process was adopted. The standard deviation indicates the average gray level difference of the tissues contained in the real slice F_n. Besides, a coefficient K (0 < K < 1), which defines a level of tolerance, is used. The value of T_g is given by:

where:

• σF_n - standard deviation between the levels of gray from the slice F_n;

• K - constant value, 0 < K < 1, defined in the experiments.

The objective behind the initial representation is to decrease the size of the inpainting domain, which in essence can be the whole virtual slice. By decreasing the size of the domain D via PDEs, the processing time will also be decreased. In this manner, the more similar two adjacent real slices, are, the smaller the inpainting domain will be. In order to increase the similarity of two real slices, an anisotropic diffusion is applied, with edge preservation of the 3D image, which is composed through the joining of the real slices before the similarity analysis.

Some pioneering work for anisotropic diffusion was proposed by Malik and Perona (1988) for the removal of noise from images. Following this work, various other researches appeared in the scientific literature with the aim of carrying out diffusion, maintaining edge preservation (Alvarez et al., 1992; Barcelos and Batista, 2007). In this work we use the mean curvature flow (Ecker, 2004), described by:

where I₀ is the original image and ∇I the gradient of I.

The proposed model uses a diffusion term and a forcing term. The balance between these terms is carried out in a selective manner on the image, where the edge and the uniform region are treated in different manners. This approach was presented by Barcelos et al. (2003) and described by:

where g = g (|∇ I|) and I₀ is the original image.

The function g(s) > 0 is a nonincreasing function and satisfies the conditions: g(0) = 1 and g(s) → 0, when s → ∞. The function g most commonly used is:

where K_g is a constant.

The Equation 7 balances, in a selective manner, the image's edge and the uniform region, that will receive the diffusion action in a differentiated manner. The function g stresses the action in regions away from the edge and reduces the action in regions near the edge.

As described earlier, this work aims at increasing the similarity level between two adjacent real slices by using the diffusion method on a 3D image composed through the joining of the CT slices. This procedure makes the convergence of the intensity of tissues' gray levels on the slices.

Figure 3 illustrates the virtual slices obtained with and without the diffusion step. Figures 3c and 3d show that the diffusion step aids in producing an initial representation with smaller regions containing no information and much better definition.

Interpolation through the use of inpainting in 3D

The inpainting process aims the transportation of information from regions containing information of the 3D image into regions which contain no information, that is, into D inpainting domain. Bertalmio et al. (2000) proposed an inpainting method split into two stages. The first consists of the transportation of information from the external region to the interior of the region containing no information; the second consists of the diffusion of the information already transported. Barcelos and Batista (2007) proposed a modification of the Bertalmio inpainting process, where the diffusion step is extended across the whole image, however, in the portion of the image outside of the inpainting region, the edge preservation criterion is used. In this manner, the process carries out the inpainting and noise elimination in only one process.

This work uses the approach proposed by Barcelos and Batista (2007), but in this case, the inpainting domain consists of the regions which contain no information of the 3D image obtained from the initial representation, that is, the inpainting process occurs in domains containing no information. The transportation step is carried out by:

where:

• I: Ω ⊂ R3 → [a,b] ⊂ R, D ⊂ R3,

• (I) = ΔI = div (∇I),

• ∇ ┴ I - ortogonal vector to ∇I;

• I₀ - original image.

The diffusion process is carried out differently inside and outside the domain D. For the points outside of the domain D the diffusion is carried out using Equation 7, for the points inside the domain D, the diffusion is carried out by Equation 5, the mean curvature flow.

The application of this inpainting method permits that the region D receives the information transported from the region near to the D domain and not only from the border of the D domain, thus resulting in more precise curves and the need of a small quantity of slices for interpolation.

Implementation details

In this section the implementation details are presented for interpolation in 3D. Each element from I_ijk is a value which corresponds to an image's intensity gray level I(x_i , y_j , z_k ). Let I(x_i , y_j , z_k , t_n ) = Iⁿ_ijk, where t_n = nΔt. The derivative of I in relation to time t (scale parameter), that is, I_t calculated in the space (x_i , y_j , z_k , t_n ) is approximated by the Euler method:

By using the Neumann boundary conditions, calculated for , n = 1, 2, ..., N, through the discretization of Equation 8 for the transportation and the Equations 5 and 7 for the diffusion, one obtains:

For the transport of information

where P_t(u) is obtained by the discretization of Equation 8.

• K_t - constant of attenuation;

• D - inpainting domain,

• D^C - complement of D.

For diffusion

where is obtained by the discretization using central difference techniques:

and g is the function described in Equation 7.

The numerical implementation of the transport and diffusion terms were carried out by using central difference techniques (Haralick and Shapiro, 1992) in the 3D image. The transport and diffusion terms are given respectively by:

Results

Two experiments were carried out using CT exams containing 52 and 50 slices, respectively. The slices used in the tests were images of 512 x 512 pixels with tomograph spatial resolution, which was configured so that each pixel represents 0.34 mm of the image's mapped structure. In this way, the quantity of virtual slices can be obtained by:

where:

• d(F_a,F_b) - distance between the slices F_a and F_b;

• Trunc(x) - returns the largest integer less than or equal to x.

For each of these exams, sub-sets were made with 50, 25 and 12.5% of the original size, always maintaining a constant space between the slices of each of the sub-groups. The objective behind the difference in information quantity used in each sub-group is to measure the method's efficiency when the quantity of slices is decreased for representing the same structure. For all the tests two transference functions were used to accentuate the bone structure as well as the patient's face.

For the first experiment, Figure 4 presents the images which resulted from the proposed method. In these CT exam, the slices are separated by 1 mm. To obtain the initial virtual slices the parameters used were: K = 0.15, N_d = 50 and K_g = 0.85. The parameters for the transportation and diffusion equations were: M = 25, M_t = 3, M_d = 10, K_t = 0.002, K_d = 0.05 and K_g = 0.6.

Figures 4a and 4d show the results obtained using only the even numbered slices from the original data set. In this case, the data set was reduced by 50% and a new distance between two adjacent real slices was set at 2 mm. With this in mind, and taking into consideration the equation 14, it became necessary to generate 5 virtual slices between two adjacent real slices.

Figures 4b and 4e show the results after the removal of three real slices from each group of four real slices of the original set. In this case, the data set was reduced to 25% of its original size and a new distance between each adjacent real slice is 4 mm, which meant that 11 virtual slices were necessary between two adjacent real slices.

Figures 4c and 4f show the results after the removal of seven real slices from each group of eight real slices of the original set. In this case, the data set was reduced to 12.5% of its original size and a new distance between two adjacent real slices is 8 mm, which meant that now 24 virtual slices were needed between two real adjacent slices.

Discussion

To verify the efficiency of the proposed method, the authors compared its performance to the Goshtasby et al. (1992) method and to the linear interpolation method. For the first experiment, Figures 5 and 6 present the images which resulted from the linear interpolation in the first column, Goshtasby et al. (1992) in the second column and the proposed method in the third column.

Through a visual analysis, it is possible to observe that the three methods produce a satisfactory visualization. However, the proposed method presents much more precise and smoother contours, in bone structure as well as the visualization of the patient's face.

Through the use of this data set, the 3D visualization, using the 13 real slices and the virtual slices interpolated by the proposed method, shows itself to be satisfactory. The same does not happen when interpolation is used with the other two methods (see Figures 5d, 5e, 6d and 6e).

One can observe that even though very few, the proposed method still permits structure recognition, whereas neither of the other two allowed any type of structure recognition. The visual analysis allows one to verify that the greater the distance between the real slices, the more evident the superiority of the proposed method becomes. The time spent by each algorithm was: 09 sec (Linear); 30 sec (Goshtasby et al.; 1992) and; 285 sec (proposed method).

Other experiments were carried out, where the same standard of behavior was obtained for the tested methods. So, in the second experiment only the visual results for both the proposed method and linear interpolation method are presented.

Numerical comparison

Besides the visual results, statistical correlation measurement as well as the PSNR were carried out in the obtained results. As the real slices removed from the group are substituted by virtual slices by the interpolation process, the numerical comparisons were made by using one removed real slice and comparing it to its corresponding virtual slice. Figures 7 and 8 show the graphs with the measurements realized in the first experiment. The graphs show that the transportation and diffusion process achieved better correlation and PSNR results in all the three cases studied. Besides, these measurements showed that the proposed method's behavior is maintained even after altering the distances between the real slices.

Conclusion

The use of non-linear differential equations for 3D interpolation produces well-defined and smooth contours, allowing for a perfect reconstruction of structures, through the use of virtual slices inserted between real slices.

The diffusion process made it possible to model satisfactorily the non-existence of abrupt changes in human body tissues, which validated the hypothesis used in this work.

The approximation made possible a satisfactory reconstruction even after using only 5% of real slices and 95% of interpolated slices as in the cases shown in Figures 5 and 6, which allowed for a significant reduction of the amount of real slices from the CT scans.

Recebido: 30/12/2011

Aceito: 06/04/2012

• Alvarez L, Lions P, Morel J. Image selective smoothing and edge detection by nonlinear diffusion. SIAM Journal on Numerical Analysis. 1992; 29(3):845-66. http://dx.doi.org/10.1137/0729052
• Barcelos CMZ, Batista M. Image restoration using digital inpainting and noise. Image and Vision Computing. 2007; 25(1):61-9. http://dx.doi.org/10.1016/j.imavis.2005.12.008
• Barcelos CAZ, Boaventura M, Silva Junior EC. A well-balanced flow equation for noise removal and edge detection. IEEE Transactions On Image Processing. 2003; 12(7):751-63. PMid:18237950. http://dx.doi.org/10.1109/TIP.2003.814242
• Bertalmio M, Sapiro G, Caselles V, Ballester C. Image inpainting. In: SIGGRAPH: Proceedings of the 27th Annual Conference on Computer Graphics; 2000 July 23-28; New Orleans. ACM; 2000. p. 417-24.
• Bors GA, Kechagias L, Pitas L. Binary morphological shape-based interpolation applied to 3D tooth reconstruction. Transactions on Medical Imagining. 2002; 21(2):100-8. PMid:11929098. http://dx.doi.org/10.1109/42.993129
• E-Film Solutions. Canada Merge Healthcare; 2011. [cited 2011 Oct 04]. Available from: www.merge.com/na/index.aspx
• Ecker K. Regularity theory for mean curvature flow. Basel: Birkhäuser Book; 2004. http://dx.doi.org/10.1007/978-0-8176-8210-1
• Genu AM. Tomografia computadorizada multidetectores na avaliação do tromboembolismo pulmonar: uso de reformatações em projeção de intensidade máxima [tese]. São Paulo: Universidade de São Paulo; 2007.
• Goshtasby A, Turner D, Ackerman L. Matching of tomographic slices for Interpolation. IEEE Transaction on Medical Imaging. 1992; 11(4):507-16. PMid:18222892. http://dx.doi.org/10.1109/42.192686
• Harald G, Scherzer O. Using the complex ginzburg-landau equation for digital inpainting in 2D and 3D. Scale Space Methods in Computer Vision. 2003; 2695(1):225-36.
• Haralick RM, Shapiro L. Computer and robot vision. Boston: Addison-Wesley; 1992. 630 p.
• Hernandez J, Aguilar J, Lara A. Volumetric reconstruction from dicom format in magnetic resonance imaging and 3D visualization. In: MEP 2006: Proceedings of the Multiconference on Electronics and Photonics; 2006; Guanajuato, México. Guanajuato: University of Guanajuato; 2006. p. 163-7.
• Köstler H, Prummer M, Rude U, Hornegger J. Adaptive variational sinogram interpolation of sparsely sampled CT data. In: ICPR06: Proceedings of the 18th International Conference on Pattern Recognition, 2006 Aug 20-24; Hong Kong. Hong Kong; 2006. p. 778-81.
• Malik J, Perona P. Scale-space and edge detection using anisotropic diffusion. IEEE Transactions on Pattern Analisys and Machine Intelligence. 1988; 12(7):12-49.
• Nolf E. XMedCon - An open-source medical image conversion toolkit. European Journal of Nuclear Medicine. 2003; 30(2):46.
• Nonato LG, Castelo A, Campos JP, Bíscaro HH, Minghim R. Topological tetrahedron characterization with application in volume reconstruction. International Journal of Shape Modeling. 2005; 11(2):189-216. http://dx.doi.org/10.1142/S0218654305000773
• Shuying Z, Jinzhu Y, Liying C, Guojun S, Xinhe X. An algorithm of nonlinear interpolation segmentation based on distance mask applied to human brain 3D recognition. In: Proceedings of the 2th International IEEE EMBS Conference on Neural Engineering; 2005 Mar 16-19; Arlington, VA. Arlington: IEEE; 2005. p. 5-8.
• 3D Slicer. A Multi-plataform, free and open source software package for visualization and medical image computing. [Internet]. 2010. [cited 2011 Oct 04]. Avaiable from: www.slicer.org
• Taillandier F, Deriche R. Reconstruction of 3D linear primitives from multiples views for urban areas modelisation. In: PCV02: Proceedings of the Symposium on Photogrammetric Computer Vision; 2002 Sep 9-13; Graz, Austria. Graz; 2002.
• Tbaina A, Prado A, Bueno J. 3D reconstruction of tomographic images applied to largely spaced slices. In: CBMS'97: Proceedings of the 10th IEEE International Conference on Computer Based Medical Systems; 1997; Maribor, Slovenia. IEEE; 1997. p. 353-67.
• Wang Y, Zhang Z, Guo B. 3D image interpolation based on directional. In: MMBIA2001: IEEE Workshop on Mathematical Methods in Biomedical Image Analysis; 2001; Kauai, Hawaii. IEEE; 2001. p. 195-202.

Publication Dates

• Publication in this collection
  07 Dec 2012
• Date of issue
  Sept 2012

History

• Received
  30 Dec 2011
• Accepted
  06 Apr 2012