A Support Construction for CT Image Based on K-Means Clustering ()
1. Introduction
For a long time, research regarding image field has become popular. A digitalized image can be analyzed and manipulated in various meanings. The more image quality is improved, the more image explains and provides details. In order to get a better representation of the taken picture or to improve its quality, it is indisputably that taking more pictures from different views or angles is the simple way. This principle made practical use in medical imaging field also, where an accurate internal image is obtained by combining pictures from different views. This well-known method is used not only in nuclear medicine, but also in many fields such as scientific field, engineering field, military field and so on, for example, taking or scanning shape of extremely tiny object. These pictures can be acquired from data called projections and its procedure to combine the projections together to obtain an image called image reconstruction. In this imaging system, it is required to obtain an image when the number of projection profiles is restricted in angle under some situations. One of the most widely used techniques is Computed Tomography. Computed Tomography (CT) is the technique that can generate the internal structure of a target object by using the projections of various angles. Two examples of the limited angle problem are illustrated in Figure 1. They show relationship between object domain and Fourier domain of system. In Figure 1(a), the angular range of the projection is restricted to
by obstacle. Thus, the missing areas of 2D Fourier transform of the projections are also restricted to
shown in Figure 1(b). Although 2D image can be reconstructed from an angular set of 1D projection profiles through many conventional methods, those methods still have weakness in some cases. For example, finding the zero regions in reconstructed image in case there is no any zero value in a set of 1D projection of inverse problem. An example of problem is illustrated in Figure 2 which is shown the comparison between normal case and unusual case. In this paper, we called the zero region as the support of non-object region. The objective of support construction is to find and separate region between object area and non-object area in order to enhance reconstructed image in the next step. An example of appearance of support is demonstrated in Figure 3 to show region between possible object region (white area) and possible non-object region (black area). In order to obtain a good support, firstly a hybrid method is applied for improving the quality of reconstruction from projections when the number of projections is limited or small. The technique consists of 3 traditional methods (ART, FBP and iterative Fourier trans-
Figure 1. (a) Angular range of projections is restricted by shown in object domain, (b) Fourier domain.
Figure 2. Comparing projection profiles.
Figure 3. Support appearance of phantom.
form technique), which transforms repeatedly between object domain and Fourier domain using a prior data in each iteration. The filtered back projection (FBP) [1] is well known method and based on the relationship between the projections and sections on the object and frequency domains, respectively. This technique is developed from simple back projection method. Geometrically, its operation propagates the constructed sonogram back into 2D image plane along the projection paths. Nevertheless, a result image obtained from this simple technique is blurred. Thus, high pass filter will be performed in next step. In other words, filtered back projection is a technique to correct the blurring encountered in simple back projection method. Next procedure applied in this research is one of the most famous iterative methods called ART. The algebraic reconstruction technique (ART) [2] , also called Kaczmarz’s method is often used for solving problem regarding linear system which occurs on the object space, and many developments for the method have been presented. Iterative reconstruction technique has recently become a famous technique and received much attention in this CT research field because it has many advantage compared with conventional techniques. However, each method has advantage and identity individual itself. Thus, many techniques are used and combined their advantages together to get better image quality. All of each method will be more described the deeply details in next section. In the case of all angle projections, the reconstruction problem will be well-solved without ineluctable artifact in the discrete computation.
The projection data is observed from the limited-view angles in the cases of the metallic implants in patients in medical, ocean acoustic tomography, electron microscopy of macromolecules, and so on. Then the problem is to become a kind of ill-posed type. An iterative Fourier method is to use the relationship between the given sections and prior information in Fourier and object domains. Although many kind developments for the algorithms for solving the limited view CT, the stagnation and fail to converge still are unsolved. Goal of our research is to construct support in order to find black region (pixel intensity is 0) in image reconstructed by projection profiles. Consequently, it is necessary to receive a good image as well as possible.
In this paper, we consider the combination of conventional methods for CT and focus on recent developments of the noise filtering method [3] [4] and some traditional techniques (ART, FBP and iterative Fourier method) [2] [5] - [10] , then introduce and discuss the reconstruction method that is based on ART, FBP, iterative Fourier transform and their hybrid with a numerical example.
2. Methodology
2.1. Total Variation
The total variation function is popular in several fields of mathematical image processing. The idea of the total variation has been firstly introduced as a denoising technique by Rudin, Osher and Fatemi [4] in Computer Vision. It is used as a global regularization term. Total variation is not only used for denoising, but also for more general signal restoration problems.
In this paper, Total Variation Denoising (TVD), also known as total variation regularization was used to be an approach for noise reduction which is defined in terms of an optimisation problem. In order to find the output of the TV denoising, output is obtained by minimizing a particular cost function. Although the algorithm can be solved in several different ways, the derivation is based on the min-max property and the majorization-minimization procedure given in [11] . In our experiment, we focus the numerical implementation of total variation denoising for one-dimensional discrete signals which are represented as projection profile of each perspective.
Let us consider a discrete real-valued in projection pro-file of N-point signal
defined on
. The total variation measures how much the signal changes between signal values. There are many ways to define discrete TV by means of finite different signal value, but we used absolute values (l1 norm) because it is one of the simplest ways which is defined as
(1)
Suppose
is a signal consisted by original signal
and additive white Gaussian noise
as the following.
(2)
We want to estimate the target
by using iterative clipping algorithm [3] to minimize the objective function which the general form of its can be defined as
(3)
Smoothness of the signal is controlled by the regularization parameter
.
To find parameter
, the iterative clipping algorithm for TV denoising was applied and it was clearly explained in [3] . In order to enhance image quality after performing TV method, directional Total Variational algorithm in various patterns were performed also. The proposed algorithm consists of 6 patterns of converting 2D image to 1D data shown in Figure 4. In each step of directional total variational algorithm, The TV denoising is processed on 1D data after transforming from 2D image. Then 1D data is converted back to 2D image again. The results are retained the smooth 1D projection and combined the ensemble in final step.
Figure 4. Block diagram patterns of the directional total variational algorithm. (a) top- down horizontal directional pattern based on TV filtering, (b) bottom-up horizontal directional pattern based on TV filtering, (c) top-down vertical directional based on TV filtering, (d) bottom-up vertical directional based on TV filtering, (e) zig-zag directional pattern 1, (f) zig-zag directional pattern 2.
2.2. Iterative FT
An iterative FT presented in this paper is used for combining analytical and algebraic reconstruction techniques. An iterative methods such as Algebraic Reconstruction Technique (ART) and analytical Fourier methods such as Filtered Back projection (FBP) have been considered by several applications. There are many approaches in image reconstruction field using projection data, and almost of those methods were developed from ART algorithm because of its simple principle and good effectiveness. The fundamental equation is presented by the following.
(4)
where each
is a projection along
ray; M is the total number of rays in all projections; and
is a weight of every pixel for all the different rays in the projection;
represents a column vector which contains the values of pixels and N is the total number of pixels. However, in most cases the ray width is often approximately equal to the cell width of the image and a line integral is called a ray sum. The implementation of ART [2] is formed by the Kaczmarz method due to the following.
(5)
where
and
are the current and updated vector respectively;
is the sum of pixels along
ray in
iteration; and
is the sum of projection for the
ray; and
is a coefficient.
Next, we will explain regarding algorithm of FBP solution. The FBP method uses a relationship between the Fourier transform of the projections of a target image and the correspond sections in the Fourier domain. It is based on the Radon transform.
is a polar coordinate of the object domain, and F is the Fourier transform of the target image. The relationship is presented as the following.
(6)
where
is the projection with angle
and site
, then it is represented by
and non-italic
presents imaginary unit.
We started to construct image using ART as an initial input before iterative FT was applied in next step to enhance the quality of image. The concept of iterative revision model was first provided by Gerchberg and Saxton in order to solve the problem of structure in science imaging. A plausible result is provided by the popular methods, ART and FBP, with enough constraints. However, the problem becomes to be an ill-posed type in the case of the limited view. As a good method using the Fourier domain for such the situation, an iterative Fourier transform method was presented in the GP (Gerchberg-Papoulis) algorithm [12] , and the more development was presented in later [5] . The algorithm was presented in Figure 5, and it consists of the following four steps: 1) F is obtained by the Fourier transform of a prior
; 2) replace
with
to satisfy the Fourier domain constraint; 3)
is given by the inverse Fourier transform
and 4) replace
by
that allow it to satisfy the object domain, the prior is updated. The diagram of iterative process is shown in Figure 6.
The sections are overwritten by the given projections in 2), and a prior knowledge of the target is embittered to the
in 4). The following is well used as one of the updated method in 4),
(7)
where D is the region at which
violates the object-domain constraint such the prior knowledge.
Such the updating procedure has been developed in the challenging the phase retrieval problem [8] , and the more refinement are expected using the progress than before. In the limited view of the tomographic problem, each method cannot give a good result in the meaning of the initial condition, parameters, itera-
Figure 5. An illustration of iterative method.
Figure 6. Diagram of iterative procedure.
tions, and falling to the global solution. Therefore, in this presentation, we show a good hybrid of ART, FBP and iterative Fourier transform for the limited view constraint. The algorithm can be explained by following operations:
1) Construct 2D image as initial input which is constructed by traditional method (ART).
2) Construct sinogram from projection project of each angle
.
3) Take Fourier transform of r to obtain 1D-FT.
4) Revise
by replacing 1D-FT from step 2) by 1DFT of original projection value and ramp filtering are applied.
5) Obtain inverse 1D-FT for the filtered projection for each
.
6) Take inverse FT and obtain reconstructed image by back projection method.
7) Revise reconstructed image by applying a prior data and support in image space.
8) Go to step 2).
These eight steps are repeatedly calculated until the process is stopped under some suitable condition of convergence. The criterion which is used to terminate process is the ratio of error between distance of original profiles and profiles of reconstructed image to become sufficiently small. Unfortunately, the number of iteration can be easily fixed but ε or value makes loop terminated cannot be fixed or set only one constant. Different data provide different distance. Thus, appropriate way is to select the result that provides the smallest distance of a set of projection in fixed iteration. In our presentation, the iterative method is connectively used with ART and FBP.
2.3. K-Means Clustering for Support Construction
The definition of object support (outline of the object) in real space is a region of interest where a target object is located. In order to construct support of reconstructed image, K-means clustering for support construction in diffractive imaging approach [7] was used in our work. K-means is an unsupervised clustering algorithm that can classify the input data points into multiple classes. In the beginning, the object domain is divided into two kinds of regions: one is the object support and the other is its complement. When the external shape of the object is clearly determined, the image intensities from areas except for the external dimension are set to be 0 (zero) (real-space constraint conditions). In this experiment, value of K is set equal to 2 because we will focus only 2 regions between object and non-object support. An algorithm of K-means for constructing support can be more clearly explain each step by following:
1) Calculate
from
.
2) Set an initial cluster presented in Figure 7.
3) Calculate center of C1.
a)
= center value of C1.
b)
=
.
where
is number of pixel in region C1.
4) Calculate center of C2.
a)
= center value of C2.
b)
=
.
where
is number of pixel in region C2.
5) For
is number of all pixel.
If
Else
6) Consider region of object by comparing value between C1 and C2.
7) Expand 1 pixel region around object.
8) Go to step 3 and repeat until value of C1 and C2 are not changed.
A simple schematic diagram of K-means clustering method is illustrated in Figure 8. An initial support can be randomly given and then each cluster is updated by K-means iteration which is shown an example of process result in Figure 9. The benefits of using K-means clustering method are effectively to construct support without many setting conditions.
Figure 7. An example of initial cluster.
2.4. Maximum Entropy Thresholding
In image segmentation field, information entropy that occasionally called Shanon’s entropy [12] or maximum entropy is one of the effective techniques to automatically find the optimal threshold. Moreover, it is being increasingly widely used as an optimal technique of image reconstruction in case noisy, incomplete data and elsewhere. In this paper, we applied maximum entropy to construct better support after obtaining from K-means method in previous step. While K-means method is performed to firstly construct initial support, maximum entropy creates more reasonable support.
For the single threshold, an algorithm can be summarized by the following. Suppose that value of a normalized histogram is shown in term of
which
takes integer values from 0 to N and result is converted to 8 bits gray scale image, that is:
. By using
as a threshold value, an entropy of black pixels is defined by
(8)
and an entropy of white pixels is given by
(9)
The optimal threshold maximizing the sum of above two entropies is presented as the following. An example of the result of optimal threshold is shown in Figure 10.
(10)
3. Experimental Result
The performance of the proposed algorithm for an ex-ample of a test signal (projection profile at 0˚) with SNR of 30 dB is shown in Figure 11. The original signal is shown in Figure 12(a). The noisy version of the signal is shown in Figure 12(b) and the output result after performing noise reduction is shown in Figure 12(c). We set the proper regularization parameter lambda. The error results are summarized in Table 1. Then, reconstructed image results were illustrated in Figure 4. In Table 1, we compare error of distance between original projection profiles and projection profiles of reconstructed image. Projection profiles of reconstructed image can be obtained 2 ways between reconstructed image by using FBP under Non-filtering method conditions and proposed method (Total Variation Denoising).
Figure 12. (a) Original projection at 0˚, (b) projection at 0˚ with noise, (c) projection profile at 0˚ after using TV denoising.
Table 1. Angles and corresponding error.
There are typical methods for the reconstruction of an unknown object using the constraints of the object and Fourier domains. For presenting the effectiveness of a hybrid procedure mixing method referred to previous section, a numerical example is settled in the following. The Shepp-Logan phantom (256 × 256) is used for a target object in Figure 12. As the numerical comparison, the result of our pro-posed method is better than traditional methods (ART and FBP) in the meaning of the error calculated by Euclidean distance of each projection profile.
In the next example from Figures 13-16, to test the robustness of the algo-
Figure 13. (a) Original image, (b) Reconstructed image with noise, (c) Reconstructed image after using TV denoising, (d) Reconstructed image after using TV denoising and support construction, (e), (f), (g) and (h), projection profiles of (a), (b), (c) and (d).
Figure 14. (a) Original image, (b) Reconstructed image corrupted by white Guassian noise, (c) Reconstructed image result after performing TV denoising method, (d) Support constructed by K-means clustering, (e) An illustrate of projection profiles at 0˚ of (a), (f) An illustrate of projection profiles at
of (b), (g) An illustrate of projection profiles at
of (c), (h) Final support result constructed by maximum Entropy thresholding and K-means clustering.
Figure 15. (a) Original image, (b) Reconstructed image corrupted by white Guassian noise, (c) Reconstructed image result after performing TV denoising method, (d) Support constructed by K-means clustering, (e) An illustrate of projection profiles at
of (a), (f) An illustrate of projection profiles at
of (b), (g) An illustrate of projection profiles at
of (c), (h) Final support result constructed by maximum Entropy thresholding and K-means clustering.
Figure 16. (a) Original image, (b) Reconstructed image corrupted by white Guassian noise, (c) Reconstructed image result after performing TV denoising method, (d) Support constructed by K-means clustering, (e) An illustrate of projection profiles at
of (a), (f) An illustrate of projection profiles at
of (b), (g) An illustrate of projection profiles at
of (c), (h) Final support result constructed by maximum Entropy thresholding and K-means clustering.
rithm various noise attack are presented to the reconstructed image. Then, we showed the appearance of support constructed by K-means method compare with maximum entropy thresholding and example of projection profiles.
From experimental results, we can deduce that traditional method does not give a plausible image. However, our proposed method provides the better results than traditional method. When the support results are used to combine with the results from iterative FT, these results are a little refined, and such the hybrid gives a stable process for the reconstruction from the limited angles projections. The proposed method will be successful as compare with other traditional methods with noise reduction. However, the results of iterative methods are different in the situation of the stagnation and fail of convergence under the limited constraint of projections. The more refinement of unifying the reconstruction algorithms and developed computation for the ill-posed imaging problems is our future work.
4. Conclusion
To summarize, in this paper we tried to show the results after reconstructing an image using projections from different angles under various noise conditions assumed as Gaussian noise. In order to achieve, we divided into 3 steps for our experiment. Firstly, the topical filter method proposed by [3] was used to reduce noise in projection profiles. Then, an iterative technique was described for tomography image reconstruction in case the number of projections is small and the angular range is limited. Final step is to combine the reconstructed image with support constructed by K-means clustering techniques widely used in diffractive imaging. The effectiveness of our proposed method was confirmed by program simulation. Future perspectives of the proposed work include the application of the implemented algorithm combined with other methods for the determination of optimal projection data.
Acknowledgements
The mathematical analysis and simulation system of our work are supported by MEXT/JSPS Kakenshi 16K00222, JST A-step AS26Z02472H and MP27115663145.