CT Scanning & Analysis

Incorporation of chemical and physical data

Computed Tomography is used to analyze the spatial distribution of components of interest in a selected control volume. Initially it was used for diagnostic purposes in the medicinal industry and it was developed as a natural extension of conventional X-ray imaging and analysis. The application of CT has spread into other industries such as endodontics (Metska, 2014) and petroleum (Funk, Jim; Al-Enezi, Sultan; Caliskan, Sinan, 2011).     

Typically a set of attenuated X-rays is measured and inverted, using Fourier Slice Theorem, into a volume of CT numbers. The distribution of CT numbers within the control volume is subsequently analyzed and often converted in to a distribution of density using a linear relation between the density and the CT number. 

The inversion algorithm (Kak & Slaney, 2001) is relatively fast, but doesn’t utilize knowledge of the chemical composition of potential components (Hubbell, J. H.; Seitzer, S. M., 1996) and introduce noise, which blur the spatial distribution and identification of the control volume content.     

Develop a program to

·        read measurements of attenuated X-rays

·        generate distribution of attenuation coefficients within control volume

·        utilize application specific a priori knowledge of potential components

·        visualize input and output

·        validate algorithm against conventional method regarding accuracy and speed

·        speed optimize

·        validate and analyze application of interest

A program InDepAna (In-Depth-Analysis) has been developed in C++ using MFC (Microsoft Foundation Classes) on a 64-bit platform.

reading Attenuated X-rays data

The program reads attenuated X-rays in raw data format (8-bit; 32-bit; 64-bit; little/big endian) and allows the user to assign the instrumentation and specimen configuration.

 For validation purposes the program can generate a data set of attenuated X-rays either from a Shepp & Logan type phantom (Shepp & Logan, 1974) or from an imported data set of attenuation coefficients in 2D or 3D.

The phantom has a default configuration (see Table 1) which can be modified by the user in an Excel sheet subsequently imported by the program.     

Table 1 Default phantom location and attenuation value

Definition
Sub-phantoms	10
Attenuation scaling	1

#	x_0	y_0	z_0	A	B	C	AngleH	AngleV	Attenuation
	m	M	m	m	m	m	deg	deg	dHU
1	0	0	0	0.69	0.92	0.9	0	0	2556.4
2	0	0	0	0.6624	0.874	0.88	0	0	-1202.1
3	-0.22	0	-0.25	0.41	0.16	0.21	108	0	-401.7
4	0.22	0	-0.25	0.31	0.11	0.22	72	0	-1126.0
5	0	0.35	-0.25	0.25	0.21	0.41	0	0	882.9
6	0	0.1	-0.25	0.046	0.046	0.046	0	0	-898.3
7	-0.08	-0.605	-0.25	0.046	0.023	0.02	0	0	-1064.2
8	0.06	-0.605	-0.25	0.046	0.023	0.02	90	0	-1097.9
9	0.06	-0.105	0.625	0.056	0.046	0.1	90	0	-165.9
10	0	0.1	0.625	0.056	0.046	0.1	0	0	-165.9

 

The imported data set of attenuation coefficients can be given in raw data format (8-bit; 32-bit; 64-bit; little/big endian) in units of HU or 1/m. An example of an imported data set and the generated measurements is shown in Figure 1 and Figure 2.

Figure 1 Vertical slice of imported 3D attenuation data set (x186 y312 z210).

Figure 2 Generated attenuated X-rays at detector plate (x433 y433) for a specific angle.

The measurement configuration is set by the user in the instrumentation dialog and the specimen dialog as shown in Figure 3, 4, and 5. The instrumentation and specimen configuration data are key input together with the measurements of the attenuated X-rays.

Figure 3 Instrumentation and specimen volume incl. phantom location shown in 2D and 3D view.

Figure 4 Instrumentation dialog.

Figure 5 Specimen dialog.

generatING Attenuation coefficients volume

From a given set of attenuated X-rays data and the corresponding measurement configuration, the distribution of attenuation coefficients is generated. There are 6 options to choose from for setting the geometry of the measurement configuration.

·        2D line source

·        2D point source – equidistant detector pixels

·        2D point source – equi-angle detector pixels

·        3D line source

·        3D point source – equidistant detector pixels

·        3D point source – equi-angle detector pixels

 

The program generates the distribution of attenuation coefficients using the Fourier Slice Theorem (Kak & Slaney, 2001) and post-processing considering knowledge of chemical composition of potential components.

utilizING a priori chemical knowledge

For most applications there is a limited number of possible components within the control volume. This program has a database with X-ray mass attenuation coefficients for 92 elements at a broad range of photon energies and a default database with density and elemental composition for 64 predefined components commonly occurring in the medicinal and endodontic specimens.

The user can add new components - by assigning the density and the elemental composition - to the database and select which components should be included in the analysis as shown in Figure 6.

The chemical and physical data of the possible components are used indirectly in the post-processing of attenuation coefficients and during analysis of the attenuation distribution as e.g. shown in Figure 1.

Figure 6 Component dialog.

visualization

The distribution of attenuation coefficients can be visualized in 2D, 3D, as iso-surfaces, and in frequency diagrams.

·        2D slice (Figure 1)

·        3D volume (Figure 7 & 8)

·        2D center slice of iso-contour (Figure 9)

·        static iso-surface for cementum (H_37.8C_12.9N_5.4O_25.2P_2.7S_2.1Ca_4.5) (Figure 10)

·        dynamic sequence of iso-surfaces (Figure 11)

·        exported iso-surface for enamel (Figure 12)

·        frequency within selected volume (Figure 13)

·        neighboring frequency (Figure 14)

 

Figure 7 3D view of attenuation data set from top (x186 y312 z210).

 Figure 8 3D view of attenuation data set after deselecting 6 components (x186 y312 z210).

 Figure 9 2D view of iso-contour of same components selection (x186 y312 z210).

 

Figure 10 Iso-surface of dentine [a: 2067.23 HU @277.6 keV].

Figure 11 Four snapshots of dynamic iso-surface sequence [a: 1375.01, 3135.01, 3775.01, 5855.01 HU @277.6 keV].

Figure 12 Iso-surface of cementum visualized with VRML [a: 3155.47 HU @277.6 keV].

Figure 13 Distribution of attenuation values with a frequency greater than 0% and smaller than 16%.

Figure 14 Neighboring frequency with a frequency coloration between 0% (blue) and 100% (red).

The exported iso-surface can be visualized using the WRML, which is a freeware application that can be downloaded from http://freewrl.sourceforge.net/download.html.

Examples of the dynamic sequence of iso-surfaces and an exported WRML file can be accessed through the links below.

·        dynamic sequence of iso-surfaces  http://www.123dtech.com/indepana_dyniso.avi

·        using WRML to visualize exported data http://www.123dtech.com/indepana_exp.wrl 

 

validation

The program can, as mentioned earlier, generate a data set of attenuated X-rays (Figure 17), either from a Shepp & Logan type phantom (Shepp & Logan, 1974) as shown in Figure 15-16, or from an imported data set of attenuation coefficients in 2D or 3D. These data set are subsequently compared to the image reconstruction data (Figure 18-20) in a Line of Perfect Agreement diagram as shown in Figure 21-23.

Figure 15 Slice of 3D phantom attenuation data set (x256 y256 z256).

Figure 16 3D view of phantom attenuation data set (x256 y256 z256).

 

Figure 17 Generated attenuated X-rays at detector plate (x256 y256) for a specific angle.

Figure 18 Slice of 3D reconstructed attenuation data set (x256 y256 z256).

Figure 19 Slice of 3D reconstructed post processed attenuation data set with fill in(x256 y256 z256).

 

Figure 20 3D view of reconstructed and post-processed attenuation data set (x256 y256 z256).

Figure 21 Line of Perfect Agreement diagram for phantom and reconstructed (inversion) attenuation data.

Figure 22 Line of Perfect Agreement diagram for phantom and reconstructed post processed attenuation data without fill in.

Figure 23 Line of Perfect Agreement diagram for phantom and reconstructed post processed attenuation data with fill in.

 

speed optimization

The timing of key modules (Figure 24-25) is reported and will allow for further speed optimization once the prototype has been developed.

Figure 24 Timing of image reconstruction algorithms.

Figure 25 Timing of building measurement data (specimen: x186 y312 z210, detector: x256 y360 z256 - 1 processor) .

The program is currently taking advantage of parallel implementation and fit-for-purpose algorithms when available.

1.     Image reconstruction is implemented with a parallel algorithm utilizing distributed memory (MPICH, 2014).

2.     Inversion algorithm is using Fast Fourier Transformation (Teukolsky, Flannery, Vetterling, & Press, 2007).

3.     Least Square Method reconstruction algorithm is solved with a conjugate gradient iterative solver and implemented in parallel with distributed memory.

4.     The dynamic sequence of iso-surfaces display is using a Marching Cubes algorithm (Young, et al., 2008), implemented in CUDA, and using the GPU for calculation and visualization.  

5.     The calculation of the neighboring frequency diagram is using a Marching Cubes algorithm, implemented in CUDA for accessing the GPU.  

6.     The calculation of fictive measurements of attenuated X-rays given the distribution of attenuation coefficients in the control volume is implemented in CUDA for accessing the GPU.  

7.     The creation of the Least Square Method Harwell-Boeing matrix is implemented in CUDA for accessing the GPU. 

The hardware being used is a Dell M6600, 8 processors, 16 GB RAM, 64-bit, NVIDIA Quadro 4000M running under Windows 7 Ultimate.

 

application validation and analysis

The ultimate validation and analysis is application specific.

endodontics:

1.     Component and composition of interest

a.       Organic

b.      Mostly solids

c.      Age dependent composition

d.      Class (incisors, canine, premolars, molars) dependent composition

2.     Dose of radiation should be kept at a minimum

3.     High resolution (micron ?) to identify cracks

Petroleum:

1.     Component and composition of interest

a.       Inorganic and organic

b.      Mixture of solid and fluid

c.      Phase, temperature, and pressure dependent fluid composition

d.      Mineral dependent composition

4.     Low dose of radiation is not imperative

5.     Efficient turn-around time to allow dynamic reconstruction of fluid flow

6.     Integration with vapor, liquid, (liquid, solid) equilibrium calculations

7.     High resolution (micron) to identify cracks and interfaces between different fluid phases

A sequential use of the inversion technique and a post-processing of the distribution of attenuation coefficients considering chemical composition has been developed and is ready for testing. The inversion technique introduces noise and the Least Square Method is computational quite heavy and requires a high resolution to be accurate. The slowness lies in both setting up the system of equations and solving the equations.

A set of measured raw attenuated X-rays passing through a control volume having organic components is needed to evaluate whether a direct solver reveals a diagnostic friendly pattern in general and in the neighboring frequency diagram in particular. The measured raw data set needs to be supplemented with a description of the measurement configuration. Also it would be preferable if it comes with an inversion image reconstruction data set for comparison purposes.

Extensive testing, debugging and further speed-optimization are also needed.

Finally when and if all of the above is finished satisfactorily, the application specific validation and analysis should be performed.       

I would like to thank Marelli Metska for introducing me to CT-scanning in general, pre- and post-processing software and data handling, and the use of CBCT (Metska, 2014) in endodontics. I also appreciate Jim Funks for introducing me to the use of CT-scanning for core-analysis in petroleum science, and Aon Khamees for showing me the CT scanning lab at EXPEC ARC Saudi Aramco.  

Funk, Jim; Al-Enezi, Sultan; Caliskan, Sinan. (2011, Oct). Core Imaging – Twenty Five Years of Equipment, Techniques, and Applications of X-Ray Computed Tomography (CT) for Core Analysis. Retrieved Oct 23, 2014, from SPE - Los Angeles Section: http://www.laspe.org/pdf/LASPE_Forum(Oct'11).pdf

Hubbell, J. H.; Seitzer, S. M. (1996). Tables of X-Ray Mass Attenuation Coefficients and Mass Energy-Absorption Coefficients from 1 keV to 20 MeV for Elements Z = 1 to 92 and 48 Additional Substances of Dosimetric Interest. Retrieved Apr 10, 2010, from National Institute of Standards and Technology (NIST): http://www.nist.gov/pml/data/xraycoef/

Kak, A. C., & Slaney, M. (2001). Principles of Computerized Tomographic Imaging. Philadelphia: Society of Industrial and Applied Mathematics.

Metska, M. E. (2014). Diagnosis and decision making in endodontics with the use of cone beam computed tomography. Amsterdam: Department of Section Endodontology, Academic Centre for Dentistry (ACTA).

MPICH. (2014). MPICH. Retrieved Oct 24, 2014, from MPICH: http://www.mpich.org/

Shepp, L. A., & Logan, B. F. (1974). The Fourier reconstruction of a head section. IEEE Transactions on Nuclear Science, Vol. NS-21 (3), 21-43.

Teukolsky, S., Flannery, W. T., Vetterling, W. T., & Press, W. T. (2007). Numerical recipes 3rd edition: The art of scientific computing . Cambridge: Cambridge University Press.

Young, P. G., Beresford-West, T. B., Coward, S. R., Notarberardino, B., Walker, B., & Abdul-Aziz, A. (2008, Jun 23). An efficient approach to converting three-dimensional image data into highly accurate computational models. Phil. Trans. R. Soc. A, pp. 3155-3173.