- Email: [email protected]
A proposed algorithm for adaptive computer tomography
Appl. Radiat. Isot. Vol.48, No. 10-12,pp. 1451-1458, 1997 © 1997ElsevierScienceLtd. All rights reserved Printed in Great Britain PII: S0969-8043(97)00140-1 0969-8043/97 $17.00+ 0.00
Pergamon
A Proposed Algorithm for Adaptive Computer Tomography S. I W A S A K I * , S. O D A N A K A , T. M I T S U I a n d M. K I T A M U R A Department of Quantum Science and Energy Engineering, Tohoku University, Aramaki-aza-Aoba, Aobaku, Sendai 980-77, Japan We propose a quite simple, recursive algorithm of image reconstruction in computer tomography (CT) based only on Bayes' theorem. Object image estimation is carried out with detection or coincidence event-base for emission type CT, and projection event-base for transmission type CT. This event-base method implies, in principle, real-time image reconstruction in CT, and provides a new paradigm called structure adaptive CT, such as parameter adjustable and regions-of-interestCTs. © 1997 ElsevierScience Ltd. All rights reserved
Introduction The technology of computer tomography (CT) has grown rapidly and spread to various applied areas. CT imaging technology is in general based on conventional types of image reconstruction techniques such as FBP (filtered back-projection method, or equivalently the convolution method) and other iterative methods, e.g. the algebraic reconstruction technique (ART), the maximum entropy method, and other methods, and technical development seems to be almost saturated. In this paper, however, we propose a simple recursive algorithm of image reconstruction in CT based only on Bayes' theorem (e.g. see Bernardo and Smith, 1994). In this algorithm, object image estimation is carried out with detection or coincidence event-base for the emission type CT (ECT; S P E C T = s i n g l e photon ECT, PET = positron ECT, etc.), and projection eventbase for the transmission type CT (TCT) using X-ray, y-ray, and neutron beam, etc. This event-base method implies, in principle, real-time image reconstruction in CT, in contrast to all the above conventional methods which are applied to the collected projection data after measurement. In this approach, CT parameters such as beam size and scan mode, etc. can be changed freely and adjusted in an optimum fashion during data collection. This adaptability of the present method provides a new paradigm in CT. For example, it enables us to scan globally an object area with a larger size of mesh and to identify regions of interest (ROI) very fast, e.g. probable defect area in the object in X-ray CT. Then, a more detailed mesh scan *To whom all correspondence should be addressed. Fax: + 81-22-217-7900.
can be made by focusing on the ROI, the defect part of the object. Such ROI setting in the object can not be made in usual reconstruction algorithms. In this adaptive environment, the efficiency of the CT measurement could be increased considerably. Furthermore, the adaptive CT can be combined easily with various advanced techniques related to image recognition (e.g. recognition of the configuration of objects in a container, and identification of defects, etc.); such combined system can be called an intelligent CT. In the present paper, we first give the theoretical base of the Bayesian estimation method, and then show some examples for both ECTs and TCTs. Finally, the simple adaptive CT for the case of the TCT will be shown for a case involving detection of a line type defect in an object.
Image Reconstruction of ECT by the Posterior Accumulation Method Bayes' theorem enables us to make reasonable estimate on the statistical plausibility of a hypothesis as a product of its likelihood and a priori distribution. If we have a reliable likelihood value of the hypotheses and a series of events generated under the hypotheses, recursive estimation of the probability of each hypothesis can be made for each event by cyclic use of Bayes' theorem: the a posteriori probability is used as the new apriori probability for the next event; it usually starts from a prior of an equal probability distribution, known as the non-informative prior. The estimate converges to the true distribution after a sufficient number of events. A typical example of this method was shown in the textbook of data analysis by Smith (1991) in the parameter estimation of a Poisson distribution.
1451
1452
S. lwasaki et al.
However, lwasaki (1996) has shown that such a simple method, i.e. replacement of the prior by the previous posterior failed in more complex problems. He developed an alternative approach called the cumulative method, where the posterior values are accumulated at a rate of constant ~ (the constant is at present an empirical one), re-normalized to unity, and used as the prior values for the next event. This cumulative method has shown remarkable success in the problems of spectrum analysis in radiation measurement and unfolding process, e.g. in the neutron spectrum measurements with a liquid organic scintillator (Iwasaki, 1996). The method can be applied to inverse problems in broad areas of radiation measurement and technology, which are expressed in the form of the Fredholm integral equation of the first kind. In CT, the image reconstruction problem can also be expressed in the same form as done in the A R T method; actually the equation is expressed in matrix form. For example, in ECT, the f: a vector of the concentration distribution of activity in each pixel of an object, which is the target of the problem, is connected to the projection data (detection) vector p with a matrix R, whose element is a contribution of each pixel to a line of the projection (defined by collimator in front of the detector) called response: p = Rf. In the present approach, the concentration distribution f: is regarded as a kind of probability of hypotheses (pixels) in Bayes' theorem as an usual probabilio,. Each element of the response matrix R gives a likelihood as the conditional probability of the contribution of each pixel to the projection beam under a condition of the pixel given. The likelihood can be estimated from the rate of contribution to an event (e.g. coincidence event in PET, and a scan beam line in X-ray CT) prior to the CT measurement.
effects are critically important, but are not for the test of the present method. Further, in image estimation we neglected the decay of the activity of the positron emitters in the object during the PET measurement.
Examph" 1 We assumed a complex two-dimensional distribution of the positron emitters (phantom): a large rectangular wall shaped source, and four spot-like sources with different sizes and intensities in and out of the wall (shown in Fig. 2(a)). The rectangular source is similar to an active area in an organ of human body, and the spot sources simulate abnormally active points, e.g. due to cancers. We simulated the process of positron annihilation in the above distribution of" the sources, just like the calculation of the likelihood, and generated a series of coincidence event data from 1 to total number of events, N. We applied the proposed accumulation method to the above data, event by event, and estimated the source distribution from N = 0 to 100 000. At the beginning of the estimation at N = 0 the uniform probability distribution ( = 1/256) was adopted as the prior as mentioned above. In Fig. 2, some examples of estimated distribution are shown for N = 100, 5000, 50 000 and 100 000. As the number of events increased, the estimate gradually approached the true value. Around 1000 counts the approximate distribution was obtained, and after 50 000 events, the distribution was almost the same as the true one. In an almost converged result (Fig. 2 (f)), the estimate reproduces well the true value in terms of shape and intensity. In the above case the ~ value was 0.0001. Survey of the ~ value dependence of the quality of the estimation is shown in Fig. 3. The figure shows the
Examples of the Present Method in ECTs Let us take an example of ECTs, PET (the proposed method is the same in the other types of ECT). In PET, the event is given a coincidence count of two signals from the detectors around a patient, each of which receives a 511 keV annihilation photon. The two-dimensional PET system is shown in Fig. 1. A total of 64 detectors surround a rectangular object area, which is divided into M × M square meshes, where M = 8 and 16 for the PET cases. The efficiency of each detector against the 511 keV photons is assumed to be 100%. The likelihood of each pixel to the coincidence event between any pair of detectors was evaluated using a Monte-Carlo simulation technique, because the simple geometrical estimation is not so accurate. In the likelihood calculation, 10 000 pairs of annihilation photons were uniformly generated in each pixel. It was assumed that there was no absorption and no scattering of the photons in the media of the object. Practically, these
positron 17 16 a n n i h i l a t i ° ~ t point .
32[ '
~
33
I I I I I . . . . . . . .
III III IIII
k
I I
I
IIIII
!
ectors
II 116
IIII1 IIIli
~1 64
'v ,~
f coincidence line 48
49
Fig. 1. Two-dimensional PET simulation for the test of the algorithm.
An algorithm for adaptive CT
1453
0.1
o~g
phantom
(a)
i::]
,~
,.
0.08
0 0 4 ~
o.o~t
~
(b)
°
0020 20
15
20
x
1
s.s.
o o
15
r~--v
Example 2
5,
In this example, we tested the proposed method in the case of the large dynamic range of the source strengths. In Fig. 4, the phantom of this test is shown; the five point sources distribute randomly in an objective area (M = 8), with their intensities of 1, 0.1, 0.01, 0.001, and 0.0001, respectively. The estimated result is shown in the right side after N = 200 000 events. The original phantom distribution is almost completely reproduced. Note the ordinate of the figure is now in log scale, and the zero value (phantom) or lower values than 10 6 of the probability of each pixel were all set 10 -6 for display purposes.
4.5. 4,
0.04 o o
0.03
N=100
(d)
0.01
2g 15 0"~0. O~08
0.o~
20 o
o
o
o
N=5000
(el
o] lo 0.1 006
$
N=50000
(f)
0.06 0.04
20
15 $
~
standard deviation of the estimated distribution from the true one as a function of N. The quality of the converged estimation is almost the same for the two values, 0.001 and 0.0001. In the former case, convergence is faster, while the quality of the estimation is lower due to the fluctuation. For = 0.01, the quality is unacceptable due to the large fluctuation of the intensity in the wall part in particular. Thus, as the value of ~ increases, the rate of convergence increases as does the size of the fluctuation. Thus, if we can change the a value at certain Ns at turning points as shown in the same figure as 'mix', the estimation can reach a stable solution efficiently: we started first with a = 0.01, changed to ~ = 0.001 at N = 1000 and finally to = 0.0001 at N = 30 000.
~0
o N=100000 Fig. 2. Algorithm test in PET for a large rectangular shape source and four spot-like sources with different intensities. Examples of the estimated distribution are shown at N = 100, 5000, 50 000 and 100 000.
Image Reconstruction of X-ray T C T by a New Recursive Bayesian Method It was assumed that the X-ray T C T problem was treated in the condition of the so-called first generation type one (a set of no overlapping parallel scanning with a narrow beam, and angle rotation) for simplicity. However, the essence of the following discussion is not different in other types of modern scanning modes using the fan or cone beams. Event in T C T is the projection data measurement using a narrow line X-ray beam. F r o m one event, the total linear attenuation coefficient along the beam path on the object area (the sum of the linear attenuation coefficients of all pixels covered by the beam path) can be exactly evaluated using the Radon transform equation. Thus, the total linear attenuation coefficient is not the statistical quantity but the physical one. Further, the beam scan is not a statistical phenomenon like the coincidence event in PET because such beam scan is set by a program of a CT system or an operator. Therefore, the accumulation method adopted in ECT cannot be applied to the T C T case, and a different approach was necessary within the framework of the Bayes' theorem based estimation.
1454
S. Iwasaki et al.
xl0
-3 I
I
I
I
8
-i
7
-.....
- ~ ,,~=o.omj
~
or,=.-o.oool
.'"-~.~..
6
,.i.<3
5
4
.... ;','~=.
e+--mix.
1
m~/i;V"-.~.
.....
....
.;,,_
:...:.. : ..-. :'.'-':+:'.~ . . . . .
.................
l
0
;, _ ' : : 7,
.--.--:-.'..7 ........................ . . . . . . . : . :::.:-:-.:.---. ~-+ .
.._ : - =...
I
I
6
8
I
2
4
",.= .. .~_::.r.-"--.¢~vi:
10
x 10 4
N Fig. 3. Standard deviation of the measured distribution vs the number of data, N, for various values of the variable :(. We again propose a new method for the estimation of the linear attenuation coefficient of every pixel for the T C T after getting data of a beam projection using the following formula:
x (a posteriori probability from Bayes' theorem The above can be expressed as ~,~,K)=
(Estimate of linear attenuation coefficient for a pixel covered by a scan beam line) = (Total linear attenuation coefficient along the beam line from the projection data)
+t
h
+., '°1
A
,o-
It
p,?(w.,,,~
,>)/(~,,_,+w,,.,~l,K
where 6#~,~) is the estimate at the Kth projection given by the Bayes estimation for a part of the attenuation coefficient of the nth pixel from the projection data
+"
A
,o-,i '°'i
.~
,o-i _.j \'~::~/
0
0
phantom
'9
0
Ii It
~'( ?'1 \IL
0
e s t i m a t e at N = 2 0 0 0 0 0
Fig. 4. Test of the present algorithm for the case of a large dynamic range of source strength; their intensities are 1, 0.1, 0.01, 0.001, and 0.0001, respectively. Note that the ordinate axis is in logarithmic scale.
An algorithm for adaptive CT p~N~(ratio of the transmitted X-ray flux to the original one) by the mth beam line; W,mis the area of the pixel partly covered by the mth beam line, which is equal to the likelihood of each pixel to the line; and the estimate of the attenuation coefficient of the nth pixel at the Nth projection line #~N~is given by merging 6pc,m with the original one #~N-~ as ~"~ = (1 - w,~)~,~- ,7 + 6~,~ If the line of a beam covered the full area of the pixel under consideration, i.e. W,,, = 1, the attenuation coefficient of the pixel n would be totally replaced by the new one:
Thus, the linear attenuation coefficient itself cannot be the probability, but its relative value can be used as the probability. The meaning of the above Bayes' estimation formula for 6t~u) is that the sum of the linear attenuation coefficients can be shared rationally between the pixels in proportion to the product of the prior coefficient values and likelihood. Of course, the estimation starts with the uniform prior. It can be imagined that the estimated values approach the true values gradually as new projection data with different angles are given.
Example 3 The above method was applied to a simple two dimensional problem for a phantom having an attenuation coefficient distribution as shown in Fig. 5(a). A square shaped tower is surrounded by a square wall, and the height of the tower is double that of the wall. This example tests the quantification capability of the proposed method by assessing whether it can reproduce the shape of the phantom with the absolute CT values of each pixel, i,e. the heights of the tower and wall, and furthermore, those of the vacant areas between them and outside the wall. The objective area was meshed in 36 × 36, and the size of the beam line for the scan was the same as the pixel size. CT measurement was simulated, and a projection data set was generated under the condition of a special scanning mode, called orthogonal scan, where the order of scanning angle was 0 '~, 90 °, 10°, 100°, 20°, 110°, etc. The statistical fluctuation in the projection data was neglected similarly to the actual X-ray TCT. Figure 5(b) shows the estimated result after scanning at 0 °, parallel to one axis of the graph. Figure 5(c) shows the result after scanning at 90 °, along the other axis. In this figure, the wall is almost reproduced, but the CT values in the wall are strongly affected by the existence of the tower and other sides of the wall. Figure 5(e) shows the result after 36 scans in 10° steps. In this case, the CT values and shape of both the tower and wall become similar to the phantom, but strong fluctuation still exists. The estimated result using the same data of the 36 angle
1455
scans four times is shown in Fig. 5(d), showing almost the same shape and CT values as the phantom. Figure 5(f) is the result obtained using the FBP method, shown for comparison. Although the FBP algorithm used was not tuned optimally to the present example, the quality of the FBP result is not as good as that of the proposed method. It can be seen that the edges of the tower and wall are rather rounded due to the filter used, there are pixels with finite CT values for the vacant areas inside and outside the wall, and there are many line type imaginary structures centering at the four corners, while the present method gives almost zero values for the vacant area. The CT value fluctuation is of the same order of magnitude with the two methods. Thus, the quantification capability of the proposed method has been proved to be high.
ROI Method in TCT As one of the adaptive CT, we present an example of the ROI method. Let an object like a uniform metal plate have a narrow line defect; the size of the defect was assumed to be 1 x 11 pixels in the finest mesh (64 × 64). Our task was to find it as fast as possible. We applied the usual method shown above (with a narrow beam corresponding to 64 × 64 mesh and 64-angle scan) and the ROI method. In the ROI method, we start with the largest mesh size (8 × 8) and eight-angle scan, and set an ROI on the pixels which satisfy the probable defect condition (under a preset value in CT value), and make the mesh size smaller in the binary step: 16 × 16 mesh and 16-angle scan, and scan only around the ROI. We then find new defect candidate pixels and set new ROI, and the mesh size smaller to 32 × 32 and 32-angle scan, and so on, down to 64 × 64 and 64-angle scan. Figure 6 shows the result of the ROI method. In the example, at 8 × 8 mesh, we found two doubtful pixels (line enclosed area in (b)) and set ROI (below 95% of the normal CT value), and changed size to 16 x 16. Scanning only around the ROI pixels made the doubtful pixels near the right upper corner clear, and set ROI on the five pixels (line enclosed area in (c), below 90% of the normal CT value), and change the size to 32 x 32, and so on. Finally, we found the true defect as shown in the figure (f), the darkest 11 pixels in 64 × 64 mesh. It was found that the quality of the final result (f) was almost the same as the ordinary scan method (not shown) in 64 × 64 mesh; both the ROI and ordinary methods succeeded in finding a defect, the location of which coincided with the true one. If we use a smaller mesh size, the number of projections required becomes larger in proportion to the number of mesh times number of angles. However, if we limit the range of the scan to a certain small area (ROI), the total number of projections could be reduced considerably. In this example, the total number of projections needed was 8192 for the
1456
S, Iwasaki et al.
0
E"= e-
()
©
, . 0 e0"~:
o.~
",'-
An algorithm for adaptive CT
(a)
1457
(b)
J
(e)
(c)
(d)
Fig. 6. Test of the ROI method for finding a line type defect in a uniform object in TCT: (a) original phantom having a line type detect; (b) estimate after 8 x 8 mesh scan, and ROI set (enclosed two pixels in 8 x 8 mesh); (c) after 16 x 16 mesh scan around the previous ROI, and new ROI set (enclosed five pixels in 16 x 16 mesh); (d) after 32 x 32 mesh scan around the previous ROI, new ROI set (enclosed 14 pixels in 32 x 32 mesh); (e) final estimate after 64 x 64 mesh scan around the previous ROI.
ordinary method, but only 2206 for the ROI method. Thus, the ROI method gives about a fourfold gain in efficiency in terms of the number of projections required. We are satisfied with the meaningful gain even if the value is not so large.
Conclusions The recursive image reconstruction method based on Bayes' theorem showed successful results in both types of CTs in the simulation studies, as shown above. Features of the method are as follows: • quite simple method, easy to calculate without numerical ill-condition; • event based method, which could realize the real-time CT imaging; • necessary quantities, the likelihood values which can be estimated prior to the data taking by the simulation calculation. Furthermore, in medical imaging using PET and SPECT, it is usual to take C T images of other types ( M R I or X-ray CT) for the same patient. F r o m these, information on the distribution of the mass attenuation coefficients of the patient prior to the ECT application can be obtained. If we adopt the large scale Monte-Carlo calculation to estimate the likelihood, we can take into account the attenuation
and scattering effect naturally and correctly through the evaluation of the likelihood. Thus, the proposed method has the possibility of solving the difficulty of the correction of both effects found in the conventional image reconstruction method (FBP). • A preliminary test of the present method for an ill-conditioned CT, a case of limited angle projection data collection, showed a nice tolerant characteristic to the limitation (not shown in this paper). • Due to the flexibility of the present method, a new concept of CT (structure adaptive CT) was proposed, where the scan beam size and scan mode, etc. could be adjusted to the object structure and quality of the data. As an example of the adaptive CT, the ROI CT tbr defect finding was successful. We are now at a very early phase in the development of the proposed method, and further comprehensive study is needed. First, it is necessary to test the method for more realistic problems in both CTs (ECT and TCT) using actual data from existing CT systems for the various types of phantom. Although the present method is quite simple and is an inherently event-base estimation method, i . e . a true real-time method in principle, the image reconstruction speed using a conventional type workstation with the highest computation power could not reach the true real time mode in the large
1458 scale CT, It would, to pursue algorithm
S. Iwasaki et al. PET, and three dimensional therefore, be interesting in an acceleration method of the which would provide true real
PET, etc. the future estimation time CTs.
References Bernardo, J. M. and Smith, A. F. M. (1994) Bayesian Theory. John Wiley and Sons, Chichester, UK.
Smith D. L. (1991) Probability, Statistics, and Data Uncertainties in Nuclear Science and Technology. OECD/NEA, Nuclear Data Committee Series, Vol. 4, p. 45. American Nuclear Society, Inc., LaGrange Park, 1L. lwasaki, S. (1997) A new approach for unfolding problems based only on the Bayes' Theorem. Proc. Ninth Int. Syrup. on Reactor Dosimetry, 2-6 Sep., 1996, Prague, Czech Republic, pp. 245-252. World Scientific Publishing, Singapore.
Pergamon
A Proposed Algorithm for Adaptive Computer Tomography S. I W A S A K I * , S. O D A N A K A , T. M I T S U I a n d M. K I T A M U R A Department of Quantum Science and Energy Engineering, Tohoku University, Aramaki-aza-Aoba, Aobaku, Sendai 980-77, Japan We propose a quite simple, recursive algorithm of image reconstruction in computer tomography (CT) based only on Bayes' theorem. Object image estimation is carried out with detection or coincidence event-base for emission type CT, and projection event-base for transmission type CT. This event-base method implies, in principle, real-time image reconstruction in CT, and provides a new paradigm called structure adaptive CT, such as parameter adjustable and regions-of-interestCTs. © 1997 ElsevierScience Ltd. All rights reserved
Introduction The technology of computer tomography (CT) has grown rapidly and spread to various applied areas. CT imaging technology is in general based on conventional types of image reconstruction techniques such as FBP (filtered back-projection method, or equivalently the convolution method) and other iterative methods, e.g. the algebraic reconstruction technique (ART), the maximum entropy method, and other methods, and technical development seems to be almost saturated. In this paper, however, we propose a simple recursive algorithm of image reconstruction in CT based only on Bayes' theorem (e.g. see Bernardo and Smith, 1994). In this algorithm, object image estimation is carried out with detection or coincidence event-base for the emission type CT (ECT; S P E C T = s i n g l e photon ECT, PET = positron ECT, etc.), and projection eventbase for the transmission type CT (TCT) using X-ray, y-ray, and neutron beam, etc. This event-base method implies, in principle, real-time image reconstruction in CT, in contrast to all the above conventional methods which are applied to the collected projection data after measurement. In this approach, CT parameters such as beam size and scan mode, etc. can be changed freely and adjusted in an optimum fashion during data collection. This adaptability of the present method provides a new paradigm in CT. For example, it enables us to scan globally an object area with a larger size of mesh and to identify regions of interest (ROI) very fast, e.g. probable defect area in the object in X-ray CT. Then, a more detailed mesh scan *To whom all correspondence should be addressed. Fax: + 81-22-217-7900.
can be made by focusing on the ROI, the defect part of the object. Such ROI setting in the object can not be made in usual reconstruction algorithms. In this adaptive environment, the efficiency of the CT measurement could be increased considerably. Furthermore, the adaptive CT can be combined easily with various advanced techniques related to image recognition (e.g. recognition of the configuration of objects in a container, and identification of defects, etc.); such combined system can be called an intelligent CT. In the present paper, we first give the theoretical base of the Bayesian estimation method, and then show some examples for both ECTs and TCTs. Finally, the simple adaptive CT for the case of the TCT will be shown for a case involving detection of a line type defect in an object.
Image Reconstruction of ECT by the Posterior Accumulation Method Bayes' theorem enables us to make reasonable estimate on the statistical plausibility of a hypothesis as a product of its likelihood and a priori distribution. If we have a reliable likelihood value of the hypotheses and a series of events generated under the hypotheses, recursive estimation of the probability of each hypothesis can be made for each event by cyclic use of Bayes' theorem: the a posteriori probability is used as the new apriori probability for the next event; it usually starts from a prior of an equal probability distribution, known as the non-informative prior. The estimate converges to the true distribution after a sufficient number of events. A typical example of this method was shown in the textbook of data analysis by Smith (1991) in the parameter estimation of a Poisson distribution.
1451
1452
S. lwasaki et al.
However, lwasaki (1996) has shown that such a simple method, i.e. replacement of the prior by the previous posterior failed in more complex problems. He developed an alternative approach called the cumulative method, where the posterior values are accumulated at a rate of constant ~ (the constant is at present an empirical one), re-normalized to unity, and used as the prior values for the next event. This cumulative method has shown remarkable success in the problems of spectrum analysis in radiation measurement and unfolding process, e.g. in the neutron spectrum measurements with a liquid organic scintillator (Iwasaki, 1996). The method can be applied to inverse problems in broad areas of radiation measurement and technology, which are expressed in the form of the Fredholm integral equation of the first kind. In CT, the image reconstruction problem can also be expressed in the same form as done in the A R T method; actually the equation is expressed in matrix form. For example, in ECT, the f: a vector of the concentration distribution of activity in each pixel of an object, which is the target of the problem, is connected to the projection data (detection) vector p with a matrix R, whose element is a contribution of each pixel to a line of the projection (defined by collimator in front of the detector) called response: p = Rf. In the present approach, the concentration distribution f: is regarded as a kind of probability of hypotheses (pixels) in Bayes' theorem as an usual probabilio,. Each element of the response matrix R gives a likelihood as the conditional probability of the contribution of each pixel to the projection beam under a condition of the pixel given. The likelihood can be estimated from the rate of contribution to an event (e.g. coincidence event in PET, and a scan beam line in X-ray CT) prior to the CT measurement.
effects are critically important, but are not for the test of the present method. Further, in image estimation we neglected the decay of the activity of the positron emitters in the object during the PET measurement.
Examph" 1 We assumed a complex two-dimensional distribution of the positron emitters (phantom): a large rectangular wall shaped source, and four spot-like sources with different sizes and intensities in and out of the wall (shown in Fig. 2(a)). The rectangular source is similar to an active area in an organ of human body, and the spot sources simulate abnormally active points, e.g. due to cancers. We simulated the process of positron annihilation in the above distribution of" the sources, just like the calculation of the likelihood, and generated a series of coincidence event data from 1 to total number of events, N. We applied the proposed accumulation method to the above data, event by event, and estimated the source distribution from N = 0 to 100 000. At the beginning of the estimation at N = 0 the uniform probability distribution ( = 1/256) was adopted as the prior as mentioned above. In Fig. 2, some examples of estimated distribution are shown for N = 100, 5000, 50 000 and 100 000. As the number of events increased, the estimate gradually approached the true value. Around 1000 counts the approximate distribution was obtained, and after 50 000 events, the distribution was almost the same as the true one. In an almost converged result (Fig. 2 (f)), the estimate reproduces well the true value in terms of shape and intensity. In the above case the ~ value was 0.0001. Survey of the ~ value dependence of the quality of the estimation is shown in Fig. 3. The figure shows the
Examples of the Present Method in ECTs Let us take an example of ECTs, PET (the proposed method is the same in the other types of ECT). In PET, the event is given a coincidence count of two signals from the detectors around a patient, each of which receives a 511 keV annihilation photon. The two-dimensional PET system is shown in Fig. 1. A total of 64 detectors surround a rectangular object area, which is divided into M × M square meshes, where M = 8 and 16 for the PET cases. The efficiency of each detector against the 511 keV photons is assumed to be 100%. The likelihood of each pixel to the coincidence event between any pair of detectors was evaluated using a Monte-Carlo simulation technique, because the simple geometrical estimation is not so accurate. In the likelihood calculation, 10 000 pairs of annihilation photons were uniformly generated in each pixel. It was assumed that there was no absorption and no scattering of the photons in the media of the object. Practically, these
positron 17 16 a n n i h i l a t i ° ~ t point .
32[ '
~
33
I I I I I . . . . . . . .
III III IIII
k
I I
I
IIIII
!
ectors
II 116
IIII1 IIIli
~1 64
'v ,~
f coincidence line 48
49
Fig. 1. Two-dimensional PET simulation for the test of the algorithm.
An algorithm for adaptive CT
1453
0.1
o~g
phantom
(a)
i::]
,~
,.
0.08
0 0 4 ~
o.o~t
~
(b)
°
0020 20
15
20
x
1
s.s.
o o
15
r~--v
Example 2
5,
In this example, we tested the proposed method in the case of the large dynamic range of the source strengths. In Fig. 4, the phantom of this test is shown; the five point sources distribute randomly in an objective area (M = 8), with their intensities of 1, 0.1, 0.01, 0.001, and 0.0001, respectively. The estimated result is shown in the right side after N = 200 000 events. The original phantom distribution is almost completely reproduced. Note the ordinate of the figure is now in log scale, and the zero value (phantom) or lower values than 10 6 of the probability of each pixel were all set 10 -6 for display purposes.
4.5. 4,
0.04 o o
0.03
N=100
(d)
0.01
2g 15 0"~0. O~08
0.o~
20 o
o
o
o
N=5000
(el
o] lo 0.1 006
$
N=50000
(f)
0.06 0.04
20
15 $
~
standard deviation of the estimated distribution from the true one as a function of N. The quality of the converged estimation is almost the same for the two values, 0.001 and 0.0001. In the former case, convergence is faster, while the quality of the estimation is lower due to the fluctuation. For = 0.01, the quality is unacceptable due to the large fluctuation of the intensity in the wall part in particular. Thus, as the value of ~ increases, the rate of convergence increases as does the size of the fluctuation. Thus, if we can change the a value at certain Ns at turning points as shown in the same figure as 'mix', the estimation can reach a stable solution efficiently: we started first with a = 0.01, changed to ~ = 0.001 at N = 1000 and finally to = 0.0001 at N = 30 000.
~0
o N=100000 Fig. 2. Algorithm test in PET for a large rectangular shape source and four spot-like sources with different intensities. Examples of the estimated distribution are shown at N = 100, 5000, 50 000 and 100 000.
Image Reconstruction of X-ray T C T by a New Recursive Bayesian Method It was assumed that the X-ray T C T problem was treated in the condition of the so-called first generation type one (a set of no overlapping parallel scanning with a narrow beam, and angle rotation) for simplicity. However, the essence of the following discussion is not different in other types of modern scanning modes using the fan or cone beams. Event in T C T is the projection data measurement using a narrow line X-ray beam. F r o m one event, the total linear attenuation coefficient along the beam path on the object area (the sum of the linear attenuation coefficients of all pixels covered by the beam path) can be exactly evaluated using the Radon transform equation. Thus, the total linear attenuation coefficient is not the statistical quantity but the physical one. Further, the beam scan is not a statistical phenomenon like the coincidence event in PET because such beam scan is set by a program of a CT system or an operator. Therefore, the accumulation method adopted in ECT cannot be applied to the T C T case, and a different approach was necessary within the framework of the Bayes' theorem based estimation.
1454
S. Iwasaki et al.
xl0
-3 I
I
I
I
8
-i
7
-.....
- ~ ,,~=o.omj
~
or,=.-o.oool
.'"-~.~..
6
,.i.<3
5
4
.... ;','~=.
e+--mix.
1
m~/i;V"-.~.
.....
....
.;,,_
:...:.. : ..-. :'.'-':+:'.~ . . . . .
.................
l
0
;, _ ' : : 7,
.--.--:-.'..7 ........................ . . . . . . . : . :::.:-:-.:.---. ~-+ .
.._ : - =...
I
I
6
8
I
2
4
",.= .. .~_::.r.-"--.¢~vi:
10
x 10 4
N Fig. 3. Standard deviation of the measured distribution vs the number of data, N, for various values of the variable :(. We again propose a new method for the estimation of the linear attenuation coefficient of every pixel for the T C T after getting data of a beam projection using the following formula:
x (a posteriori probability from Bayes' theorem The above can be expressed as ~,~,K)=
(Estimate of linear attenuation coefficient for a pixel covered by a scan beam line) = (Total linear attenuation coefficient along the beam line from the projection data)
+t
h
+., '°1
A
,o-
It
p,?(w.,,,~
,>)/(~,,_,+w,,.,~l,K
where 6#~,~) is the estimate at the Kth projection given by the Bayes estimation for a part of the attenuation coefficient of the nth pixel from the projection data
+"
A
,o-,i '°'i
.~
,o-i _.j \'~::~/
0
0
phantom
'9
0
Ii It
~'( ?'1 \IL
0
e s t i m a t e at N = 2 0 0 0 0 0
Fig. 4. Test of the present algorithm for the case of a large dynamic range of source strength; their intensities are 1, 0.1, 0.01, 0.001, and 0.0001, respectively. Note that the ordinate axis is in logarithmic scale.
An algorithm for adaptive CT p~N~(ratio of the transmitted X-ray flux to the original one) by the mth beam line; W,mis the area of the pixel partly covered by the mth beam line, which is equal to the likelihood of each pixel to the line; and the estimate of the attenuation coefficient of the nth pixel at the Nth projection line #~N~is given by merging 6pc,m with the original one #~N-~ as ~"~ = (1 - w,~)~,~- ,7 + 6~,~ If the line of a beam covered the full area of the pixel under consideration, i.e. W,,, = 1, the attenuation coefficient of the pixel n would be totally replaced by the new one:
Thus, the linear attenuation coefficient itself cannot be the probability, but its relative value can be used as the probability. The meaning of the above Bayes' estimation formula for 6t~u) is that the sum of the linear attenuation coefficients can be shared rationally between the pixels in proportion to the product of the prior coefficient values and likelihood. Of course, the estimation starts with the uniform prior. It can be imagined that the estimated values approach the true values gradually as new projection data with different angles are given.
Example 3 The above method was applied to a simple two dimensional problem for a phantom having an attenuation coefficient distribution as shown in Fig. 5(a). A square shaped tower is surrounded by a square wall, and the height of the tower is double that of the wall. This example tests the quantification capability of the proposed method by assessing whether it can reproduce the shape of the phantom with the absolute CT values of each pixel, i,e. the heights of the tower and wall, and furthermore, those of the vacant areas between them and outside the wall. The objective area was meshed in 36 × 36, and the size of the beam line for the scan was the same as the pixel size. CT measurement was simulated, and a projection data set was generated under the condition of a special scanning mode, called orthogonal scan, where the order of scanning angle was 0 '~, 90 °, 10°, 100°, 20°, 110°, etc. The statistical fluctuation in the projection data was neglected similarly to the actual X-ray TCT. Figure 5(b) shows the estimated result after scanning at 0 °, parallel to one axis of the graph. Figure 5(c) shows the result after scanning at 90 °, along the other axis. In this figure, the wall is almost reproduced, but the CT values in the wall are strongly affected by the existence of the tower and other sides of the wall. Figure 5(e) shows the result after 36 scans in 10° steps. In this case, the CT values and shape of both the tower and wall become similar to the phantom, but strong fluctuation still exists. The estimated result using the same data of the 36 angle
1455
scans four times is shown in Fig. 5(d), showing almost the same shape and CT values as the phantom. Figure 5(f) is the result obtained using the FBP method, shown for comparison. Although the FBP algorithm used was not tuned optimally to the present example, the quality of the FBP result is not as good as that of the proposed method. It can be seen that the edges of the tower and wall are rather rounded due to the filter used, there are pixels with finite CT values for the vacant areas inside and outside the wall, and there are many line type imaginary structures centering at the four corners, while the present method gives almost zero values for the vacant area. The CT value fluctuation is of the same order of magnitude with the two methods. Thus, the quantification capability of the proposed method has been proved to be high.
ROI Method in TCT As one of the adaptive CT, we present an example of the ROI method. Let an object like a uniform metal plate have a narrow line defect; the size of the defect was assumed to be 1 x 11 pixels in the finest mesh (64 × 64). Our task was to find it as fast as possible. We applied the usual method shown above (with a narrow beam corresponding to 64 × 64 mesh and 64-angle scan) and the ROI method. In the ROI method, we start with the largest mesh size (8 × 8) and eight-angle scan, and set an ROI on the pixels which satisfy the probable defect condition (under a preset value in CT value), and make the mesh size smaller in the binary step: 16 × 16 mesh and 16-angle scan, and scan only around the ROI. We then find new defect candidate pixels and set new ROI, and the mesh size smaller to 32 × 32 and 32-angle scan, and so on, down to 64 × 64 and 64-angle scan. Figure 6 shows the result of the ROI method. In the example, at 8 × 8 mesh, we found two doubtful pixels (line enclosed area in (b)) and set ROI (below 95% of the normal CT value), and changed size to 16 x 16. Scanning only around the ROI pixels made the doubtful pixels near the right upper corner clear, and set ROI on the five pixels (line enclosed area in (c), below 90% of the normal CT value), and change the size to 32 x 32, and so on. Finally, we found the true defect as shown in the figure (f), the darkest 11 pixels in 64 × 64 mesh. It was found that the quality of the final result (f) was almost the same as the ordinary scan method (not shown) in 64 × 64 mesh; both the ROI and ordinary methods succeeded in finding a defect, the location of which coincided with the true one. If we use a smaller mesh size, the number of projections required becomes larger in proportion to the number of mesh times number of angles. However, if we limit the range of the scan to a certain small area (ROI), the total number of projections could be reduced considerably. In this example, the total number of projections needed was 8192 for the
1456
S, Iwasaki et al.
0
E"= e-
()
©
, . 0 e0"~:
o.~
",'-
An algorithm for adaptive CT
(a)
1457
(b)
J
(e)
(c)
(d)
Fig. 6. Test of the ROI method for finding a line type defect in a uniform object in TCT: (a) original phantom having a line type detect; (b) estimate after 8 x 8 mesh scan, and ROI set (enclosed two pixels in 8 x 8 mesh); (c) after 16 x 16 mesh scan around the previous ROI, and new ROI set (enclosed five pixels in 16 x 16 mesh); (d) after 32 x 32 mesh scan around the previous ROI, new ROI set (enclosed 14 pixels in 32 x 32 mesh); (e) final estimate after 64 x 64 mesh scan around the previous ROI.
ordinary method, but only 2206 for the ROI method. Thus, the ROI method gives about a fourfold gain in efficiency in terms of the number of projections required. We are satisfied with the meaningful gain even if the value is not so large.
Conclusions The recursive image reconstruction method based on Bayes' theorem showed successful results in both types of CTs in the simulation studies, as shown above. Features of the method are as follows: • quite simple method, easy to calculate without numerical ill-condition; • event based method, which could realize the real-time CT imaging; • necessary quantities, the likelihood values which can be estimated prior to the data taking by the simulation calculation. Furthermore, in medical imaging using PET and SPECT, it is usual to take C T images of other types ( M R I or X-ray CT) for the same patient. F r o m these, information on the distribution of the mass attenuation coefficients of the patient prior to the ECT application can be obtained. If we adopt the large scale Monte-Carlo calculation to estimate the likelihood, we can take into account the attenuation
and scattering effect naturally and correctly through the evaluation of the likelihood. Thus, the proposed method has the possibility of solving the difficulty of the correction of both effects found in the conventional image reconstruction method (FBP). • A preliminary test of the present method for an ill-conditioned CT, a case of limited angle projection data collection, showed a nice tolerant characteristic to the limitation (not shown in this paper). • Due to the flexibility of the present method, a new concept of CT (structure adaptive CT) was proposed, where the scan beam size and scan mode, etc. could be adjusted to the object structure and quality of the data. As an example of the adaptive CT, the ROI CT tbr defect finding was successful. We are now at a very early phase in the development of the proposed method, and further comprehensive study is needed. First, it is necessary to test the method for more realistic problems in both CTs (ECT and TCT) using actual data from existing CT systems for the various types of phantom. Although the present method is quite simple and is an inherently event-base estimation method, i . e . a true real-time method in principle, the image reconstruction speed using a conventional type workstation with the highest computation power could not reach the true real time mode in the large
1458 scale CT, It would, to pursue algorithm
S. Iwasaki et al. PET, and three dimensional therefore, be interesting in an acceleration method of the which would provide true real
PET, etc. the future estimation time CTs.
References Bernardo, J. M. and Smith, A. F. M. (1994) Bayesian Theory. John Wiley and Sons, Chichester, UK.
Smith D. L. (1991) Probability, Statistics, and Data Uncertainties in Nuclear Science and Technology. OECD/NEA, Nuclear Data Committee Series, Vol. 4, p. 45. American Nuclear Society, Inc., LaGrange Park, 1L. lwasaki, S. (1997) A new approach for unfolding problems based only on the Bayes' Theorem. Proc. Ninth Int. Syrup. on Reactor Dosimetry, 2-6 Sep., 1996, Prague, Czech Republic, pp. 245-252. World Scientific Publishing, Singapore.