Navigation path calculation algorithm
Virtual colonoscopy (VC) is a diagnostic method enabling the generation of two-dimensional and three-dimensional images of the colon and rectum from the data obtained with relevant imaging modality, usually spiral computed tomography (CT). If CT is used, the method is also called CT colonoscopy, CT colonography, or CT pneumocolon. The main advantages of the VC which support its broader application in medical practice include: limited invasiveness, improved compliance of patients and value for screening for colorectal cancer.
The introduction of virtual endoscopy technique originates from the extended processing options of the data sets obtained with available imaging modalities. The Visible Human Project and related activities (Hong et al., 1996; Hong, 1997) were of key importance for development of the VC.
The quality of first VC images limited the potential of the clinical use of this technique. But with next generations of the imaging equipment and advanced processing algorithms, their applicability in medical practice was established. VC may be performed with all imaging techniques which result in cross-sections of the abdominal cavity. It can be generated both from CT and magnetic resonance imaging (MRI) cross-sections. However, the CT remains the first choice imaging modality for the VC. It is required that spiral acquisition mode with overlapping reconstructions is applied. The quality of the images generated within VC techniques depends strongly on the spatial resolution obtained with imaging modality. Nowadays, multidetector computed tomography (MDCT) is commonly available and this adds to overall quality of assessment of the colon and rectum resulting from VC.
The patient undergoing helical computed tomography with the intent of obtaining VC should undergo complete bowel preparation as for other procedures within abdomen, e.g. endoscopic colonoscopy. The priority is assigned to evacuation of the contents of the colon before CT. For this purpose, many agents are used including ethylene glycol electrolyte solution, magnesium citrate or oral sodium phospate. Netherless, the quality of bowel preparation for VC varies considerably between different centres (Van Uitert et al., 2008).
The trend for the optimisation of the diagnostic procedures and limitation of the burden to the patient resulted also in a strategy focusing on the performing of the optical colonoscopy just after VC, if it is positive for pathological lesions in the colon, in order to avoid repetition of the bowel preparation procedure. A strategy enabling identification of the artefacts resulting from fecal contents in the bowel in the process of generation of VC images was also proposed. This is achieved by labelling it with some type of contrast agent, e.g. barium or meglumine diatrizoate taken orally before the CT (Iannaccone et al., 2004).
The use of MRI imaging for generating of VC is another attempt to avoid bowel preparation and exposure to the radiation related to the CT imaging (Florie et al., 2007). Finally, the techniques enabling “electronic” colon cleansing before generating VC images were developed to allow for less intensive colon preparation procedures (Lakare et al., 2002). VC without laborious bowel cleansing preparation preceding is related to higher acceptance and compliance from patients. This in turn may be a key condition for successful screening strategy in the society.
Growing use of the VC is supported by its lower invasiveness in the comparison to other diagnostic procedures and potential for higher compliance from patients. These features increase the value of the techniques as a screening test for disorders of the colon. The main indications for VC include screening for colonic polyps or cancer and failure or inadequate results of optical colonoscopy due to anatomical conditions or pathological lesions, e.g. obstruction of the colon lumen. Furthermore, the VC enables also for examination of extracolonic structures not accessible during standard colonoscopy. This may be particularly important for these patients in whom pathological lesions were detected inside the colon lumen.
The first studies focusing on the assessment of the sensitivity of the VC reported usually its lower performance in comparison to optical colonoscopy. The sensitivity of the VC depends considerably on the size of the polyps present in the colon and is lower in less advanced lesions. Introduction of MDCT had significant impact on improving the VC efficiency. Most recent studies indicate that VC sensitivity in detecting polyps of size at least 10 mm is comparable with optical colonoscopy and exceeds 90% (Regge et al., 2009; Graser et al., 2009).
The Guidelines issued in 2008 by American Cancer Society, American College of Radiology and US Multi-Society Task Force on Colorectal Cancer included VC within recommended screening tests for colorectal cancer, which should be performed at 5 years intervals in population of at least 50 years or older (Levin et al., 2008). According to the Guidelines, VC should be performed after complete bowel preparation. The detection of a polyp of size >6 mm in VC necessitates the performance of optical colonoscopy, preferably the same day or second complete bowel preparation is needed.
During classical (optical or video-) endoscopy, the endoscope is inserted into the internal space of a tube-like organ. Colonoscopy is performed with flexible endoscope equipped with a camera. A physician performing examination can bend the tip of the endoscope in two orthogonal directions using small wheels which are placed on the head of the instrument in order to navigate and move it through the colon. Colonoscopy may be associated with different level of discomfort in specific patients and depend on the medication used during preparation and the procedure itself. The VC may be an alternative approach to classical endoscopy. In this chapter, technical aspects of the generation of the VC images are explored. Authors describe all steps which are required to create a 3D computer model of the colon from the CT data. A procedure which allows one to compute the VC is presented in figure 1. It includes CT examination of abdomen, electronic colon cleansing, generation of 3D computer model of the colon with appropriate segmentation techniques. Then a navigation path for virtual camera may be calculated in order to simulate the progression of a real endoscope. Finally, views from colon reconstructed from the CT data are calculated based on visualisation methods.
2. Computed tomography data
The imaging modality of computed tomography (CT) was developed in the results of the chain of discoveries starting from the year 1895 when Wilhelm Röntgen invented a new type of electromagnetic radiation, called by him as X-rays. The mathematical principles of CT were first investigated by Johann Radon in 1917 and then extended by Kirillov in 1961. The first CT scanner was presented in 1972 by Allan Cormack and Godfrey Hounsfield, who were awarded the Nobel Prize in medicine in 1979. The detailed description of CT was provided in 1988 (Kak & Slaney, 1988). Nowadays, CT imaging is commonly used throughout medical specialities for diagnostic purposes and support of interventional procedures. Because of easier hardware realization, the fan-beam projection is used in medical CT scanners. In the third-generation CT scanners, the X-ray source and the detector array are rotated around the patient. In helical or spiral cone beam CT scanners, the patient is lying during the procedure on a moving bed and the X-ray source and the detector arrays are rotating around him. The helical scan method enables for quicker scanning and reduction of the radiation dose (necessary for given resolution) to which the patient is exposed during the procedure. Additionally, the number of rows of detectors in helical CT scanners was increasing from several years. Nowadays, the 64 or 128 multi-slice scanners are frequently encountered in clinical applications. Currently, 256-slice cardiac CT scanners enables for scanning of the heart in less than one second. On the other hand, minimizing the radiation risk to the patient, while maintaining satisfactory CT image quality, becomes urgent especially for colon screening with CT (Wang et al., 2008). Dose reduction for CT imaging can be achieved by scanning patient whit low-mAs protocols (less than 100mAs). Unfortunately, these protocols may result in noisy and streak artefacts in the reconstructed images. This effect can be compensated by the optimisation of data acquisition and the application of iterative image reconstruction algorithm (Wang et al., 2008).
The dataset acquired with the CT scanner can by described by number of slices, the number of pixels per slice and the voxel distances. The number of pixels in one slice is also referred to as image matrix. It is usually set of 512 x 512 pixels. The distances between the voxels are differentiated into the slice distance (out-of-plane) and pixel distance (in-plane). In general, the medical data are anisotropic (pixel and slice distances are not equal).
The data resolution influences the noise level. The data having higher resolution are more noisy for the same radiation dose. The computed intensity values represent the densities of the scanned objects that are normalized into Hounsfield units (HUs). This normalization maps the 12-bit data into two bytes (16 bits). Water is mapped to zero and air is mapped to -1000 value.
Medical images are physically stored together with the data essential for their interpretation. This information is highly standardized. The DICOM standard (Digital Imaging and Communications in Medicine) enables the integration of scanners, servers, workstations, printers, and network hardware from multiple manufacturers into a picture archiving and communication system (PACS). It includes a file format definition and a network communications protocol. The DICOM files can be read by many programs and are supported by numerous libraries (VTK, DCMTK, GDCM).
3. Colon cleansing
The first step of the CT data processing is the electronic cleansing (EC) of the colon, especially if other means were not undertaken to remove fecal contents from the colon. The term EC was first introduced by Wax and Liang (Wax et al., 1998). This is a key operation for ensuring correct segmentation being the next step. In figure 2, we can see fluid (contrast) in the colon. Classical threshold operation does not remove voxels lying on the border between air and fluid (see Fig. 2b). The idea of thresholding operation relies on voxels division into 2 groups making use of a predefined threshold value
The threshold value is usually designated on the basis of the histogram of the CT data or general knowledge about values attributed to anatomical structures. A border between contrast and tissues (see fig. 2b), unwanted but usually being obtained after thresholding operation, is a result of limited resolution of detectors and soft property of reconstruction kernel.
3.1. State of the art
In order to remove voxels representing contrast and residual nutrients, many different computer algorithms were proposed for the electronic colon cleansing (ECC). They differ with the image pre-processing steps, local image features (mainly statistical ones, e.g. 23 features in (Lakare et al., 2003) and 35 in (Cai et al., 2011)), the method of reduction of features dimensionality (vector quantization, principal component analysis), applied modelling (of multi-material objects and their edges/gradients), the use of the segmentation techniques (watershed, active contours, level-set) and classification methods (Markov random fields, expectation maximization, support vector machine). Recent publications in this area (Cai et al., 2011), (Serlie at al., 2010) provide the state-of-the-art and historical perspective of the research focused in the EC.
First works on ECC were based on the application of the statistical image features, vector quantization (for dimensionality reduction), image gradient information and Markov random field classification (Chen at al., 2000). Later methods paid special attention to the application of edge modelling during image segmentation, aiming at efficient delineation of tagged regions. The segmentation ray technique (Lakare et al., 2002) is the most important one in this group. In this method the rays were designed to analyse the intensity profile and to detect the intersection between the air and the residual fluid and between the residual fluid and the soft-tissues. The segmentation rays can accurately detect partial volume regions and remove them if necessary. The same authors (Lakare et al., 2003) proposed a method based on vector quantization but it does not assure correct exclusion of the voxels lying near colon wall. In turn, in another method (Zalis et al., 2004) the image gradient is approximated by Sobel mask filtering followed by a morphological dilation. Recently, Wang (Wang et al., 2006) proposed application of statistical features and an expectation-maximization algorithm for distinguishing voxels belonging to multiple materials while (Serlie at al., 2010) built a scale-invariant three-material transition model between air, soft-tissue and tagged material/fluid and used it for classification of each voxel. The most sophisticated and effective algorithm has been proposed recently in (Cai et al. 2011) that consists of several very carefully designed steps making use of many image features (descriptors), two segmentation procedures (watershed transform, level-set method) and very precise SVM classification/cleansing method (sensitivity 97.1%, specificity 85.3%, accuracy 94.6%).
In order to present problems solved with electronic colon cleansing methods and to demonstrate typical existing solutions, the ECC algorithm developed by the authors of this chapter is briefly described below.
3.2. Electronic colon cleansing using non-linear transfer function and morphological operations
The ECC method proposed by us is based on non-linear value transformation combined with morphological voxels processing (Skalski et al., 2007a).
First, if the CT data has HU values without offset, the voxels values are increased by 1024 and unsigned 16-bit fixed-point integer data format is obtained what results in significant reduction of the calculation time. In order to remove voxels representing contrast one has to find them in the CT data. To reach this aim, we compute two binary masks: a fluid mask and a residual mask. The fluid mask is created by thresholding operation described in section 3: voxels having values greater than 1600 are given value 1. In case of the residual mask we are looking for values greater than 1350 and equal or smaller than 1600 since voxels representing stool and fluid remain within this range. Both masks are dilated using regular hexahedron of size 3. Voxels for which the masks are equal 1 are then processed by two transfer functions, shown in fig. 3, representing Gaussian intensity transformation:
with σ = 450 for the fluid mask and 100 for the residual one. This operation is desirable due to necessity to keep smooth changes of intensity on the border between colon and soft tissues.
In the next step a binary data (
The whole operation can be summarized by the following formula:
After subtraction we check intensity profile along normal direction to the body surface (fig. 4) for each voxel equal 1 in the
Exemplary results from application of the proposed algorithm for electronic colon cleansing are shown in figure 5.
Data visualization methods like surface or volume rendering can be used for showing inner structure of the colon. However, inspection of the visualized data requires manual virtual camera movement which make smooth observation difficult. Segmented data may be used for creation of the automatic navigation path for the virtual camera. If we display only segmented structure, we can reduce the time required for visualization. Combined visualization and segmentation algorithms allows for development of 3D models of anatomical structures (fig. 6). Additionally, structures that are external to the colon can be viewed, which improves the assessment of the pathological lesions.
The human abdomen consists mainly of three regions: air, soft tissues and high density materials (bones) and this is reflected by voxels values. Thresholding represents the simplest approach to abdomen segmentation but has many disadvantages, e.g. it does not remove partial volume voxels. For example, voxels near the edge of objects are incorrectly classified when thresholding is used. Therefore, in the VC, segmentation is usually based on region growing (Vilanova et al., 1999; Xie et al., 2003, Sadleir & Whelan, 2005) or active contour methods (Jiang et al., 2005). The idea of the region growing technique is linking thresholding procedure with neighbourhood checking. In first iteration, the algorithm checks membership condition for all voxels of the neighbourhood of voxel being classified. If the voxels pass the membership condition test which can be the same as in classical thersholding procedure, the voxels are added to the object. In next iterations, this process is repeated for all voxels added in the previous iteration until no new voxel can be added. It allows for local operation in contrast to thresholding. Even if in the dataset there are voxels which can pass membership condition, they will not be classified as the object if they have no connection with voxels added before.
Different strategy is applied in the method of active contours, called also “snakes”, proposed by Kass et al. (Kass et al., 1988). The active contours method is a segmentation technique in which the problem of object finding in the analyzed data is formulated as energy minimisation. It is usually calculated in iterative routine where contour evaluation is guided by external constraint forces and influenced by image forces that pull the contour towards lines and edges present in the data. The total energy consists of the internal and the external energies which are responsible, respectively, for contour behaviour and image influence. The active contours method is a parametric technique which is susceptible to parameter tuning and this is a one of its main disadvantages. But even then classification of voxels lying near the colon wall is a source of problems.
The 3D segmentation algorithm based on immersion-based watershed method of Vincent and Soille (Vincent & Soille, 1991) can be also applied. The watershed method exploits topographic and hydrology concepts for the development of region segmentation methods. The image may be seen as a topographic relief, in which the value of a pixel (for 2D images) is interpreted as its altitude in the relief. In case of 2D, the principle of watershed algorithm can be illustrated by an idea of immersing the image from water sources. When the neighbouring catchment basins eventually meet, a dam is created to avoid the water spilling from one basin into the other (Vincent & Soille, 1991). When the water reaches the maximum value, the edges of the union of all dams form the watershed segmentation results (fig. 7). In case of 3D, usage of the algorithm leads to receiving 3D objects separated by the dam. If we use local minima of the image as water sources, oversegmentation problem will appear. In consequence, we receive a huge number of objects in the resultant matrix which do not correspond to data.
One of the solutions is a modified strategy of source selection. We used marker-based Watershed transform, where the immersion processes are started from markers computed from the image.
In order to improve results, the absolute value of gradient of the filtered data (fig. 8) is computed using the 3D Sobel’s mask and then the data are immersed by the watershed algorithm.
Neubauer et al. (Neubauer et al., 2002) proposed manual placing of markers inside the data. On the contrary, we use automatic methods for markers computation: object markers are obtained from voxels after thresholding operation. We know that voxels having intensity below 300 units represent air, and background markers are voxels which have intensity value corresponding to tissues. This approach eliminates the oversegmentation problem.
The 3D watershed segmentation algorithm computes border between the colon and the soft tissue using the gradient map (fig. 8) calculated as:
Thanks to these operations, the colon model, which is traced, reflects details very precisely what we can observe in figure 6.
5. Calculation of navigation path
Fast and accurate navigation path generation is essential for efficient diagnosis using the VC since it allows for simulation of the virtual camera movement inside the segmented structure and the whole colon can be screened by the physician in a short time. The virtual camera can be stopped if a suspicious image is discovered for more careful assessment.
Computation of the colon centerline is not a trivial process. The algorithm should require only minimum operator intervention. Additionally, a centreline approximation of the centre navigation path of the colon must be obtained in reasonable time with acceptable accuracy. Time constraint is a very important factor to evaluate the algorithm especially in a clinical practice.
5.1. State of the art
Centreline calculation methods can be subdivided into three categories since they are mainly based on:
distance transform (e.g Vilanova et al., 1999).
Manual extraction requires manual identification of the centre of colon slices. It does not guarantee that marked points lie in the centres of slices and that they are directly connected. Furthermore, the allocation is difficult because the colon centreline is oriented in different directions.
Methods based on topological thinning and distance transform are automatic usually. The idea of topological thinning is based on peeling off the colon surface points using morphological operations repeatedly until the centreline is obtained. Though the results of this standard algorithm are well-defined they do not always lie in a proper place. Additionally, the algorithm is extremely inefficient computationally. Therefore, other methods were developed that use 3D topological thinning and graph search algorithm (Ge et al. 1999), optimized 3D topological thinning using Look-up Table (Sadleir & Whelan, 2005), distance transform (Zhou & Toga 1999, Van Uitert and Bitter 2007), minimum energy path (Deschamps and Cohen 2001) or Dijkstra’s shortest path algorithm (Bitter et al. 2000).
Some approaches rely on the use of the distance transform. In these methods, the centreline is calculated as a maximum of binary mask representing colon after the distance transformation.
Below an exemplary algorithm of the VC navigation path calculation is presented.
5.2 Exemplary algorithm based on distance transform
In this section we present the colon centreline calculation algorithm based on the distance transform (Skalski et al., 2007b). As an input to the algorithm, the matrix with labelled voxels (label
|compute complement IL of the binary mask L (1 – colon, 0 – others)|
compute distance transform on IL
find maximal value
save location of
set sphere to IL (value=1; centre=(
as a result we receive points saved in the matrix path (3x
sort the path:
find a point which has a minimal value of
for (i=1 to N-2) do:
compute the distance between the last marked point path(
find minimal distance
compute interpolating cubic spline
Firstly, complement of the binary mask resulting from the segmentation process is done. It is a preparation step for the distance transform calculation. The distance transform returns as a result distance to the nearest voxel which belongs to the colon. In order to calculate the distance, the Euclidian metric is usually used:
The final step in the Virtual Colonoscopy is the data visualization. There are two main methods of 3D medical data visualization: indirect (surface rendering) and direct (volume rendering) (Preim & Bartz, 2007). Both techniques can use fully programmable graphical pipeline.
Surface rendering is one of the indirect methods. This technique produces surfaces in the domain of the scalar quantity. Scalar values, contained in 3D medical data, represents tissue properties, like radiodensity in Hounsfield scale or label mask that contains segmentation results. Surface represents a specific scalar value, the so-called isosurface value. In fact, one iso-surface describes only one scalar value. The interior of the object is not described – surface is the boundary of the volume objects. Surface rendering method includes two stages: generation of the 3D surface from 3D data and proper visualisation relying on the image generation by graphics accelerator. There are numerous methods for implementing the surfaces from a discrete set of 3D data (Preim & Bartz, 2007). One of the most useful is the Marching Cubes algorithm (Lorensen & Cline, 1987). This algorithm has many implementations that solve the problem of ambiguities in first cell triangulation method (holes in surface). Possibly the widely used implementation of the Marching Cubes algorithm comes from The Visualization Toolkit (Schroeder et al., 2004). In the algorithm a polygonal mesh of the isosurface is generated from the 3D scalar field. The polygonal mesh is a collection of vertices (points in 3D) connected to triangles. For high resolution data sets the number of generated graphical primitives can be extremely high. To reduce the number of triangles, the mesh can be decimated or smoothed (Schroeder et al., 1992 and 2004).
Surface can be coloured according to the isosurface value or to another scalar field using a texture mapping technique. To increase the perception of the surface shape in the VC visualization, the virtual lighting is used. The standard model in the OpenGL Application Programming Interface is Phong illumination model (Phong, 1975). This is an empirical model of local illumination. It describes the way a surface reflects the light as a combination of the diffuse reflection of rough surfaces with the specular reflection of shiny surfaces. Phong shading includes a model for the reflection of light from surfaces and a compatible method of estimating pixel colours by interpolating surface normals across rasterized polygons. In fact, at each point of the screen full Phong model calculations are performed (per-pixel lighting). Since the interpolation of surface normals is computationally expensive, the Phong shading is slow. In Gouraud shading algorithm the calculating lighting is performed only in vertex. Next, the screen pixel colour on the triangle are bilinearly interpolated from the vertex colour. This method is fast, but the specular highlight will not be rendered correctly if a highlight lies in the middle of a polygon. This limitation can be solved by increasing a number of triangles by mesh tessellation or by increasing of spatial data resolution. The polygonal data can be efficiently processed in modern graphics card. All shading calculations are done in hardware.
Volume rendering is the process of creating a 2D image directly from 3D volumetric data that operates on the actual data sample without creation of intermediate surfaces consisting of triangles (Preim & Bartz, 2007). The purpose of volume rendering is to effectively convey information present within the volumetric data. It is especially important in case of medical data. All direct volume rendering algorithms can be classified into two main groups: object-space and image-space methods. However, many advanced algorithms cannot by easily classified as one or the other, but fuse aspects from both groups into one hybrid algorithm.
Object-space volume rendering techniques use forward mapping scheme where the data is mapped onto the image plane. One of such approach is the Splatting algorithm that projects the data voxels onto image-plane (Westover, 1989). Texture-mapping algorithms are the other widely used object-oriented algorithms. They are supported by computer graphics hardware. In image-order (image-space) algorithms, a backward mapping scheme is used where rays are cast from each pixel in the image plane through the volume data to determine the final pixel colour. The classic direct volume rendering method is the image-space oriented ray casting algorithm. Moreover, some algorithms use domain-base techniques – the spatial volume data is first transformed into an alternative domain, such as frequency or wavelet, and then a projection is generated directly from this domain (Malzbender, 1993).
Modern graphics cards are characterized by immense ability of 3D data processing. They are developed and optimized for processing triangle meshes, which are used for surface rendering. Furthermore, a fully programmable graphics processing unit (GPU) offer new opportunities to use graphics cards for general purpose computing, especially for volume rendering. Ray-casting volume rendering using CPUs is computationally expensive since it requires the interpolation and shading calculations for every sample point along the ray in the data. Interactive volume ray casting was previously restricted to high-end workstations. GPU implementations of ray-casting rendering approaches have received great attention since they enable interactive visualization of volumetric data (Lee et al., 2009).
The most important in virtual colonoscopy visualization is trustworthy surface presentation. In figure 10, examples of applying different rendering methods are shown. The fastest method in interactive visualization is the surface rendering. Unfortunately, the triangle structure of reconstructed surface has an influence on image quality. The colour interpolation (diamond artefact) are visible. Surface in image generated by the direct volume visualization has better quality (fig. 10 b-d). The shape of wrinkles is kept. However, the computational cost is considerably larger. The texture-mapping technique, supported by the GPU acceleration, can be used in interactive visualization. Unfortunately, this method generates images which contain staircase artefacts caused by interpolation and insufficient depth sampling (fig. 10.d).
To extend the field of view in virtual colonoscopy the multi-cameras are used, especially for visual inspection of the colon wall (Serlie I. et al. 2001). The six cameras are located in the same place, but the view directions are different. They are rotated around, to cover 360 degree of view. Each camera has the 90 degree field of view. Images of this cameras can be mapped into the unrolled-box surface. The sample of this technique is shown in the figure 11. Additionally, the light source moves along with the camera and the position of light source is the same as camera position. The light is configured as positional (headlight), and the cone angle corresponds to camera cone angle. To prevent overexposing nearest surfaces, the irregular light intensity along the cone angle was used. Light fading attenuation was used for distance simulation.
Comparison between real colonoscopic image and virtual one is presented in figure 12.
In this chapter virtual colonoscopy has been presented from the point of view of computational sciences. Problems present in the VC software realization have been pointed out and their existing solutions have been cited. For clearity of presentation, to help a reader to understand merits of technical issues associated with the VC, simple examples of computer algorithms have been given, mainly developed by the authors of the chapter. Special attention has been paid to the following technical aspects: electronic colon cleansing, colon lumen segmentation, navigation path calculation and modern 3D visualisation.