Another Step Towards Successful Tomographic Imaging in Cancer: Solution of the Problem of Image Reconstruction

Tags: reconstruction, discrete image, image reconstruction, Toeplitz matrix, projection, system of linear equations, parallel rays, breast imaging, cancer, finite number, breast cancer, reconstructed image, random rectangles, parallel projections, pictures of the human body, CT scanners, lung cancer screening, Tomographic imaging, computed tomography, lung cancer, CT images
Content: Another Step Towards Successful Tomographic Imaging in Cancer: Solution of the Problem of Image Reconstruction Artyom M. Grigoryan Department of Electrical and Computer Engineering, The University of Texas at San Antonio, One UTSA Circle, San Antonio, USA ABSTRACT Tomographic imaging is an important tool which is used in cancer to detect tumors and provide information about the extent of disease, help plan treatment, and determine whether the cancer is responding to treatment. Special and complex X-ray equipments are used in tomography to obtain cross-sectional pictures of the human body, which are computed by special and complex algorithms of image reconstruction from the projections. This work describes a new approach for reconstructing images from a finite number of projections. The ray-integrals of the image f(x, y) are transformed uniquely into the ray-sums of the discrete image fn,m on the Cartesian lattice. This transformation allows for calculating the tensor, or paired representation of the discrete image, when the image is considered as the sum of direction images, or splitting-signals carrying the spectral information of the image at frequency-points of different subsets that cover the Cartesian lattice. The proposed approach is presented for parallel projections and the continuous model. The Experimental results of image reconstruction from a finite number of projections are illustrated. Keywords: Tomographic imaging, Fourier slice theorem, reconstruction by projections 1. INTRODUCTION Problem of image reconstruction from a finite number of projections is important in the computed tomography, or computerized tomography (CT) which is used in diagnostic medicine1 and electron microscopy.2,3 The technology of CT scanners is developing and new scanner generations provide high quality pictures. We mention here the third-generation scanners wherein the detector array cover the entire field of view and fan-beam geometry is used instead of parallel-beam geometry.4,5 Computerized tomography which also is called the computerized axial tomography (CAT) is a diagnostic procedure that uses a special X-ray equipment to obtain cross-sectional pictures of the human body. The CT images pictures are called tomograms, images show structures inside the body, which include the internal organs (such as the kidneys, liver or spleen), blood vessels, bones, tumors, and other tissues. Tomograms are used in cancer to detect or confirm the presence of a tumour including a brain tumour, to provide information about the size and location of the tumor and whether it has spread, to guide a biopsy (the removal of cells or tissues for examination under a microscope), and to help plan radiation therapy or surgery. 1.1 CT Images and lung cancer A wide variety of medical conditions can be seen in CT images which are not visible in conventional radiographs. According to the news from the American College of Radiology, CT can be used for lung cancer screening and the recent results of the National Lung Screening Trial have shown that CT lung cancer screening can reduces significantly the number of people who die each year from lung cancer, when CT is performed in the context of careful patient selection and follow-up. Also new research published in the January issue of The Journal of Nuclear Medicine shows that CT and Positron Emission Tomography (PET) imaging offer significant prognostic stratification information at initial staging for patients with locally advanced breast cancer. When compared to conventional imaging, 18F-FDG PET/CT more accurately showed lesions in the chest, abdomen and bones. The study in6 suggests that 18F-FDG PET/CT should be considered as complimentary imaging tools in the pre- and postoperative work-up of patients diagnosed with breast cancer. Further author information: Artyom M. Grigoryan: E-mail: [email protected], Telephone: (210)458-7518.
1.2 CT Images and Brest Cancer As reported in,7 breast cancer is a pathology that causes 40,000 deaths each year in the United States, and is the second leading cause of death in women.8,9 Early detection is necessary to treat patients and improve their chance for survival. Over the last 15 years the number of deaths has been drastically reduced by 30% to 50% because of technological advancements in mammography. Recent studies show the potential benefits of the use of computed tomography technology for breast imaging. The breast CT becomes another alternative to mammography and magnetic resonance imaging for use in breast cancer screening. There are a few complications in mammography, and the major one is the lack of sensitivity with mammography, which makes it impossible to detect some abnormalities, for instance, the breast lesions and breast cancer.10,11 Studies shown that 30% of cancers are not detected by mammography, and 70% to 90% of biopsies conducted as a result of suspicious mammograms are negative, driving up medical costs and causing patients unnecessary stress. Other complications include not being able to visualize the breast in a 3-dimensional (3D) view and how the superimposition of breast tissues with mammography makes it harder to detect abnormalities.12,13 Currently it has been recognized that CT may potentially play a major role in breast imaging, by providing advanced treatment options for patients. Experimental and research studies have centered on using breast CT as another diagnostic tool for breast cancer.14 Studies indicate breast CT can be as valuable as mammograms, if not more so. Using cone-beam technology for breast CT has been shown to benefit image quality. A simulation study conducted by Chen and Ning showed the feasibility of cone-beam dedicated breast CT systems in identifying lesions only a few millimeters in size and calcifications 100 micrometers in size.13 The breast CT can produce better low contrast images which aide in detecting breast masses. A precise measurement of the location and volume visualization of lesions are possible as well as optimal resolution for calcification with high contrast. This is an important aspect of breast CT because physicians are now able to localize the lesions, which will lead to advanced treatment options with CT, such as using breast CT for biopsies. 1.3 Brest cancer with CT, Mammography and MRI The breast CT may become another alternative to MRI for screening women who are high at risk for breast cancer and staging the extent of their disease15-.17 The cost of breast CT will be less expensive than MRI and the potential interpretation time for CT may be shorter than that for MRI. The potential advantages of breast CT systems far outweigh the disadvantages. These advantages include the ability to image the entire volume of the breast, the high spatial resolution, the superior contrast resolution, the elimination of superimposed breast tissues and structures resulting in improved tumor detection, and the fine anatomical detail produced with breast CT. Although dedicated breast CT utilizes ionizing radiation, the radiation dose is equal to or less than the dose used for mammography. Also, unlike mammography, breast CT does not need to use breast compression, which is a major benefit for the comfort of the patient. Although breast CT systems are still in development, the study of breast CT is moving in a positive direction, and it has the potential to gradually change the field of breast imaging. If breast CT becomes commercially available in the future, it will likely play an important role in screening and early diagnosis. Will breast CT replace mammography or MRI? Breast CT may have a vital role in evaluating women with dense breasts, implants, and small cancers in dense breasts. Translational research should help to provide evidence-based guidelines for clinical implementation of breast CT. Lowering patient dose, while increasing the resolution needed to detect angiogenesis at the ductal level, will help to increase survival rates among women. 1.4 Algorithms in CT With new and still developing generations of CT scanners and their effective implementation in tomographic imaging, the complexity of the mathematical methods and approaches for image reconstruction from projections increases. The problem solved by Radon in 1917,18,19 for reconstructing a function from its projections, has been not solved yet for practical applications, because an infinite set of line-integrals (measurements) is not available. Many of existent methods of image reconstruction from a finite number of projections process the projection data on the polar grid of coordinates, and then transformer the data into the Cartesian lattice, without fully covering it, which results in low quality reconstructions. This is the tradition way for calculating the approximation of the reconstructed image, which is shown by the Fourier slice theorem.1,2
In this work we present a new approach for solving the problem of image reconstruction, which differs from the known methods of back-projection, iterative reconstruction, Fourier filtering, Radon filtering, and convolution filtering20-.28 The reconstruction is performed with calculation of the tensor or paired transform29-33 of the image from the projection data. In the tensor representation, a two-dimensional (2-D) discrete image is considered as the sum of direction images, each of which is determined by the corresponding signal. These signals are called the splitting-signals; they also carry the spectral information of the image at frequency-points of different subsets that cover the whole domain of frequencies. The discrete image can be reconstructed by its splitting-signals by calculating the 2-D DFT or directly by performing the inverse 2-D tensor transform. In addition, the tensor representation possesses the following important property. This transform is in the discrete space and all components of the transform are defined as the ray-sums of the discrete image in the Cartesian lattice. However, it is very important to note that these sum-rays can exactly be calculated from the ray-integrals. In other words, the splitting-signals of the discrete image fn,m can be calculated from the projection of the image f(x, y). Each projection is processed by a system of linear equations, or linear convolutions, to calculate the corresponding part of the 2-D tensor or paired representation of the image,33,45,49 and, then, the inverse transform is calculated to obtain the discrete image. The model described for image reconstruction is simple and the reconstruction is exact. The proposed method was implemented on MATLAB and C++, and the experimental results of image reconstruction are illustrated. Our preliminary results also show high quality images, when applying the proposed method with an incomplete set of projections for the angular range scanned down to 10. The rest of the work is organized in the following way. In Section II, we consider the model of image, and, in Section III, we describe briefly the tensor and paired transforms of the discrete image of size N Ч N. Sections IV and V present the method of transferring geometry from the image plane to the lattice for the cases, when N is a prime and a power of two, respectively. The preliminary experimental results of image reconstruction by the proposed method are also given in these sections.
2. MODEL OF THE IMAGE We consider the simple discrete model of image reconstruction when the parallel scanning scheme is used, although the proposed algorithms can be generalized to the fan scanning scheme, too. Figure 1 illustrates the problem of image reconstruction, when assuming that the radiation source and detector represent themselves the points, and the rays spreading between them are straight.
source
· · · ·
·
·
measured data
·
·
original image Figure 1. Image and two sets of parallel rays.
To simplify the calculations, we consider that the 2-D image f(x, y) to be reconstructed is in the unit square [0, 1] Ч [0, 1] with the Cartesian lattice N Ч N on it. The integer N > 1 is the size of the lattice. We denote the lattice by X = XN,N = {(n, m); n, m = 0 : (N - 1)}. The square is divided by small pieces, or image elements (IE) of size x = y = 1/N each, and the knots of the lattice are in the centers of IE. Image elements are numbered by (n, m). In each image element, the value of the image is considered to be constant
fd(x, y)
=
1 (x)2
f(x, y)dxdy, IE
if (x, y) IE.
Thus the image f(x, y) is represented by fd(x, y), and the discrete image to be reconstructed on the lattice from fd(x, y) is fn,m. The values of the discrete image are defined by
fn,m =
f(x, y)dxdy = (x)2fd(x0, y0), (x0, y0) (n, m)-th IE,
(n,m)-th IE
and they are placed in the centers of the corresponding (n, m)-th IE. When N is large and the image f(x, y) is presented by its small elements, the above equations describe a good model of the image on the Cartesian lattice. Such model is considered in the well-known series expansion methods for solving the system of line integrals in the discrete case20-.25
2.1 Line-integrals and ray-sums
The sum of the discrete image along a ray l passing through the knots of the Cartesian lattice is denoted by vl, and we call them the ray-sums,
vl =
fn,m.
(1)
(n,m)l
In the above model, the line-integral along a ray l in the (n, m)-th IE is defined as follows:
wl(n,m) =
(n,m)-th
IE
f (x,
y)dl
=
(l)fd(x0,
y0)
=
(l)
1 (x)2
fn,m
=
(l)N 2fn,m,
(2)
where (x0, y0) is a point in the (n, m)-th IE. We denote by l = ln,m the length of the ray l in the (n, m)-th
IE. The ray is referred to as the line, i.e., the ray with the zero width. The line-integral along the ray l can be
written as
wl = f(x, y)dl =
wl(n,m) = N 2
(ln,m)fn,m ,
(3)
(n,m)l
(n,m)l
where the summation is performing by the image elements through which the ray l passes. Thus we have a large system of linear equations describing the relation between the line-integrals wl and the values fn,m of the discrete image. The solution of this system is unknown; it is difficult to find this solution because, in general, the values of (ln,m) are different even for the parallel rays of the same projection. However, it is important to say, that there exist sets of parallel rays with the same value of length in each IE through which the rays pass, and we will describe these rays in a moment. The rays l of such set may differ from the rays defined the ray-sums in (1) depending on the projection, but in all cases the ray-sums can be uniquely calculated from line-integrals.
2.2 Two types of parallel rays We will describe the rays which pass the unit square and the rays which pass through the knots of the Cartesian lattice. Therefore two co-ordinate systems on the plane are considered. The first system of co-ordinates (x, y) is for the image f(x, y) on the square [0, 1] Ч [0, 1], as shown in Figure 2 in part a, for the image with 13 ellipses with constant intensity each. The second co-ordinate system (n, m), where n and m are integers, is for the lattice XN,N in the square. This system is used for the discrete image fn,m, as shown in part b. Parameters x and n run from left to right, and parameters y and m run from top to down. We now consider the parallel lines on the unit square and lattice, which are parameterized by coordinates of the frequency-points. Given the frequency-point (p, s) XN,N , such that g.c.d.(p, s) = 1, we consider the lines
l(t) = lp,s(t) = {(n, m); pn + sm = t}, t = 0 : (p + s)(N - 1),
on the square lattice XN,N . These lines are referred to as the arithmetical rays. The equations of these lines on the square [0, 1] Ч [0, 1] are
l(t) = lp,s(t) =
(x, y);
px +
sy
=
t N
+
p+s 2N
,
t = 0 : (p + s)(N - 1).
f (x, y) 0
{fn,m}
0.2
50
0.4
100
0.6
150
0.8
200
1
0
0.2 0.4 0.6 0.8
1
(a)
250 50 100 150 200 250 (b)
Figure 2. Two co-ordinate systems for images f (x, y) and fn,m.
These lines are referred to as the geometrical rays, to distinguish the discrete and continuous cases. These two types of rays are denoted by l(t), and the same set of t is considered for the rays, t = 0 : (p + s)(N - 1). The generator (p, s) defines the slop, - tan-1(p/s), of these rays. The set of line-integrals W = Wp,s = {wl(t); l(t) = lp,s(t), t = 0 : (p + s)(N - 1)} is called the (p, s)-projection of the image. We also will consider other sets of (p + s)(N - 1) + 1 shifted lines which are parallel to l(t) and call them the geometric rays, too.
3. IMAGE AND THE SET OF SPLITTING-SIGNALS In this section we describe the images and their two-dimensional discrete Fourier transforms (2-D DFT) by the unique sets of 1-D signals and 1-D DFTs, respectively. The images are considered in the tensor and paired representations which are the 2-D frequency and 1-D time representations of the image29-.35 These representations are unique, invertible, and fast algorithms exist for them. Both tensor and paired transforms of images can be used effectively for calculating the 2-D DFT and other transforms, including 2-D Hartley, Hadamard, and cosine transforms34-.36,48 The tensor transform can be used in image enhancement40-,42 image de-noising,43,49,50 and discrete image reconstruction.44,45 In the tensor representation, the discrete image fn,m is the set of 1-D splitting-signals,
: {fn,m} fTp,s = {fp,s,t; t = 0 : (N - 1)} (p,s)JN,N .
(4)
Here JN,N is a set of frequency-points (p, s), or generators of the splitting-signals which are selected in a way to cover the lattice XN,N of frequency-points (p, s) with a minimum number of the following subsets:
Tp,s = (kp mod N, ks mod N ); k = 0 : (N - 1) . These sets are cyclic groups with generators (p, s). When N is a power of two, N = 2r, r > 1, the set of generators (p, s) contains 3N/2 elements and can be defined as
JN,N = (p, 1); p = 0 : (N - 1) (1, 2s); s = 0 : (N/2 - 1) .
(5)
The components fp,s,t of the splitting-signals fTp,s are the sums of the image fn,m along the parallel lines on the lattice
fp,s,t = p,s,t f =
{fn,m; np + ms = t mod N }, t = 0 : (N - 1).
(6)
(n,m)X
The binary function p,s,t(n, m) takes value 1, when np + ms = t mod N and 0, otherwise. Given (p, s), the components fp,s,t are periodic by t, i.e., fp,s,t+N = fp,s,t, for t = 0 : (N - 1). Figure 3 shows the image in part a, along with the tensor transform in b. The tensor transform is redundant,
{fx,y }
{fp,1,t }
50
50
{f1,2s,t }
20
100
100
40
60
150
150
80
100
120
200
200
50 100 150 200 250
250
250
50 100 150 200 250
50 100 150 200 250
(a)
(b)
Figure 3. (a) Image of size 256 Ч 256 and (b) two parts of the tensor transform of the image.
55
50 f4,1,t 50 100
150 45 200
40
0
50 100 150 200 250
(a)
250 50 100 150 200 250 (b)
Figure 4. (a) Splitting-signal {f4,1,t t = 0 : 255} of the image and (b) its direction image.
and the set of all splitting-signals can be divided by two parts, or matrices. The first part is with the generators {(p, 1); p = 0 : (N -1)} and the second part is for the splitting-signals with the remaining generators {(1, 2s); s = 0 : (N/2 - 1)} of the set J256,256. The splitting-signals are written along the rows in these two matrices. Figure 4 shows the splitting-signal of the image, which is generated by the frequency-point (4, 1) in part a. This splittingsignal is written in the row number 5 in the first part of the tensor transform in Figure 3. The splitting-signal fTp,s carries the spectral information of the 2-D DFT at N frequency-points of the set Tp,s, i.e., the following holhs:
N -1
Fkp mod N,ks mod N =
fp,s,tW kt, k = 0 : (N - 1),
(7)
t=0
where W = WN = exp(-2j/N ). As an example, Figure 5 illustrates the 256 Ч 256 case. The splitting-signal, or image-signal fT2,1 of length 256 is shown in part a, along with the magnitude of the 1-D DFT of the signal in c and the spectrum of the image in d. Two bright parallel lines on the spectrum show the samples at points of the set T2,1 at which the 2-DFT of the image is the 1-DFT of the image-signal. The 2-D DFT at frequency-points of this set has been amplified, in order to see location of the group and directions of the projection along which the components of the tensor are calculated as ray sums. The image after amplifying the 2-DFT at frequency-points of the set T2,1 is shown in part b. The effect of amplifying those spectral components shows the contribution of the direction image. The projections are calculated at angle = 63.4349 and the 1-D DFT is filled the 2-D DFT along three lines at angle = 90 - = 26.5651. To remove the redundancy of the tensor transform, the paired representation of the image is used, when the set of splitting-signals is defined as30,31,33
: {fn,m} fTp,s = fp,s,t; t = 0 : (N/2k - 1), g.c.d.(p, s) = 2k
. (p,s)JN,N
(8)
8
7.5
image-signal f
T
7
2,1
6.5
6
5.5
5
4.5
4 t
3.5
0
50
100 150 200 250
(a) 50
(b) 90°
1-D DFT of f
T
40
2,1

30
amplitude
20
10
0
0
50
100 150 200 250

(c)
(d)
Figure 5. (a) The image-signal corresponding to the set T2,1. (b) The image after amplifying the 2-DFT at frequencypoints of T2,1. (c) Magnitude of the 1-D DFT of the image-signal (zero component shifted to the center and truncated). (d) The 2-D DFT of the image with amplified samples of the set T2,1 .
Here JN,N is a set of 3N - 2 frequency-points (p, s), or generators and it is defined as JN,N 2JN/2,N/2 4JN/4,N/4 . . . {(N/2, N/2)} {(0, 0)}. The components of the paired splitting-signals are defined from the tensor transforms as fp,s,t = fp,s,t - fp,s,t+N/2, when t = 0 : N/2k - 1. The splitting-signals carry information of the 2-D DFT of the image at frequency-points of the subset Tp,s = {(2m + 1)(p, s); m = 0 : N/2k - 1} Tp,s. The 2-D paired transform for the general case, when N is a power of a prime, is defined similarly to the case when such prime is 2.32,49,50
The family of subsets Tp,s, when (p, s) JN,N , covers the entire lattice XN,N , and 3N/2 1-D DFTs of
splitting-signals fTp,s define completely the 2-D DFT of the image. The tensor representation is unique, and
the image can be defined through the 2-D DFT calculated by (7) or directly from the tensor transform. Indeed,
according to the principle of superposition33,46-48 the image fn,m can be composed from splitting-signals as
follows:
fn,m
=
1 r-1 1 2N 2k
f2kp,2ks,t - f2kp,2ks,t+N/2
+
1 N2
f0,0,0
,
(9)
k=0
(p,s)JN/2k ,N/2k
where we denote t = t(2kp, 2ks; n, m) = (n2kp + m2ks)modN. All components f2kp,2ks,t in this equation are defined by (6) but they can also be calculated by using the following recursive formula:
f = f + f , 2kp,2ks,2kt1
2k-1 p,2k-1s,2k-1t1
2k-1p,2k-1 s,2k-1t1+N/2
k = 1, 2, ..., r - 1,
(10)
where p or s equals 1, and t1 = 0 : (N/2k+1 - 1). In other words, these components can be calculated from 3N/2 splitting-signals fTp,s generated by the frequencies (p, s) JN,N . Each image dn,m = fp,s,(np+ms) mod N , n, m = 0 : (N - 1), is the direction image with N values which are located along the parallel lines (np + ms) mod N, that pass through the knots of the lattice. For the case when (p, s) = (3, 1), Figure 4 shows the direction image in part b, that corresponds to the splitting-signal in a.
Thus, the image fn,m is the linear combination of direction images. This composition of the image can be used for reconstructing the image from 3N/2 projections or an incomplete set of projections. It should be mentioned, that the inverse formula in (9) can also be used for other sets of 3N/2 generators (p, s), when the covering of the lattice N Ч N by the groups Tp,s is irreducible. For instance, we can take (2p, 1); p = 0 : (N/2 - 1) (1, s); s = 0 : (N - 1) . In the case, when N is prime, the inverse tensor transform of the N Ч N image fn,m can be composed from (N + 1) tensor splitting-signals as follows:
fn,m
=
1 N
N -1 f1,s,(n+sm) mod N + f0,1,m
- N E[f],
n, m = 0 : (N - 1).
s=0
(11)
The set of generators (p, s) for cyclic groups Tp,s is defined as
JN,N = (1, 0), (1, 1), (1, 2), (1, 3), · · ·, (1, N - 1) (0, 1) .
(12)
This set contains (N + 1) frequency generators. This number shows the number of splitting-signals in the tensor representation of the image.
4. GEOMETRY OF THE PROJECTIONS ON THE LATTICE In this section, we describe the method of transferring the geometry from the image plane to the Cartesian lattice. In other words, we show a way of calculating the ray-sums of the discrete image from the line-integrals of the image f(x, y). The sets of rays for these two geometries, i.e., the sets of arithmetical and geometrical rays, may be the same or different and that will depend on the angle of projections. We first consider a few examples of projections in the case when N is a prime. Example 1 (7 Ч 7). We consider the (1, 2)-projection and the N = 7 case. Figure 6 shows three parallel lines l(2), l(9), and l(16) in the unite square with the Cartesian lattice 7 Ч 7. The generator (p, s) is (1, 2). It determines the slop, arctan(p/s), of the lines. The equations of these three lines on the [0, 1] Ч [0, 1] square are
y and m
l(2) 0
0
0.14
1
l(9)
0.28
2
0.42
3
0.57
4
0.71 5
l(16)
0.85
6
1
0
1
2
3
4
5
6
0
0.14 0.28 0.42 0.57 0.71 0.85 1
x and n
Figure 6. Parallel rays on the square [0, 1] Ч [0, 1] and the 7 Ч 7 lattice.
l(2) = l1,2(2)
=
{(x, y);
x + 2y
=
2 7
+
3 14
=
7 14
},
l(9) = l1,2(9)
=
{(x, y);
x + 2y
=
9 7
+
3 14
=
21 14
},
l(16) = l1,2(16)
=
{(x, y);
x
+
2y
=
16 7
+
3 14
=
35 14
}.
In the discrete lattice 7 Ч 7, these line are described by l(2) = l1,2(2) = {(n, m); n + 2m = 2}, l(9) = l1,2(9) = {(n, m); n + 2m = 9}, and l(16) = l1,2(16) = {(n, m); n + 2m = 16}. The sums of the image along arithmetical rays are denoted by
vp,s(t) = vl(t) =
fn,m.
(n,m)lp,s (t)
The set of 19 geometrical rays that coincide with these arithmetical rays is defined as
l(t) = l1,2(t) =
(x, y); x
+
2y
=
t 7
+
3 14
,
t = 0 : 18,
where (x, y) [0, 1] Ч [0, 1]. The parallel rays n + 2m = t, where t = 0 : 18, are shown in Figure 7 on the left. The number of rays is calculated as (p + s)(N - 1) + 1 = 3· 6 + 1 = 19. To number these rays, we use the set of
y and m
0 0 0.14 1 0.28 2 0.42 3 0.57 4 0.71 5 0.85 6 1 0
l(0) l(1) l(2) l(3) l(4) l(5) l(6) l(7) l(8) l(9) l(10) l(11) l(12) l(13) l(14) l(15) l(16) l(17) l(18)
1
2
3
4
5
6
0· 1· 2· 3· 4· 5· 6·
7· 8·







10·

11· 12·






13·
14·

15· 16·
17· 18·
(for rays at angle -26.56o).
0 0.14 0.28 0.42 0.57 0.71 0.85 1 x and n Figure 7. 19 arithmetical rays for the (1, 2)-projection and the corresponding numbered control points.
control points that are also shown on the right.
19 arithmetical rays of this projection are used to determine the components {f1,2,t; t = 0 : 6} of the tensor transforms. The masks of the basic functions 1,2,t(n, m) have coefficients equal 1 on these parallel line. We consider the masks of the first two basic functions
1 0 0 0 0 0 0
0 1 0 0 0 0 0
0000010


0 0 0 1 0 0 0
0000001


0 0 0 0 1 0 0
[1,2,0]
=

0
1
0
0
0
0
0 ,
[1,2,1]
=

0
0
1
0
0
0
0 .

0
0
0
0
0
0
1


1
0
0
0
0
0
0

0 0 0 0 1 0 0
0 0 0 0 0 1 0
0010000
0001000
These masks and the parallel rays passing trough coefficients 1 are shown in Figure 8. It is not difficult to see, that seven components f1,2,t, t = 0 : 6, of the tensor transform are calculated as
f1,2,0 = v0 + v7 + v14 f1,2,1 = v1 + v8 + v15 f1,2,2 = v2 + v9 + v16 f1,2,3 = v3 + v10 + v17 f1,2,4 = v4 + v11 + v18 f1,2,5 = v5 + v12 f1,2,6 = v6 + v13
(13)
y and m y and m
l(0) 0
0
l(7)
0.14
1
0.28
2
0.42
3
0.57 4
l(14)
0.71
5
0.85
6
1 0123456
l(1) 0
0
0.14
l(8)
1
0.28
2
0.42
3
0.57
4 0.71
l(15)
5
0.85
6
1 0123456
0 0.14 0.28 0.42 0.57 0.71 0.85 1 x and n
0 0.14 0.28 0.42 0.57 0.71 0.85 1 x and n
Figure 8. Masks of functions 1,2,0 (n, m) and 1,2,0 (n, m) with the parallel rays.
where the variables vt = v1,2(t) denote the sums of the discrete image along the arithmetical rays l(t), t = 0 : 18. This system of equations can be written as f1,2,t = vt + vt+7 + vt+14, t = 0 : 6, where vt = 0, if t > 18. To calculate the signal f1,2,t from the projection data, we will derive and, then, simplify the formula describing the linear relation between the integrals w(t) and sums v(t), t = 0 : 18. First we consider the mask of the tensor function 1,2,0(n, m) with coefficients 1 of the rays l(0), l(7), and l(14) shown in Figure 9. It is not difficult to
y and m
l(0) 0 0 0.14 1
l(5) l(6)
~l(7)
l(7)
l(8)
0.28
2
0.42
3
0.57 4
l(14)
0.71
5
0.85
6
1
0
1
2
3
4
5
6
0
0.14 0.28 0.42 0.57 0.71 0.85 1
x and n
Figure 9. A few parallel rays for calculating line-integrals of the (1, 2)-projection.
see, that the line-integral along the ray l = l(6) in image elements can be written as
wl(0,3)
=
l0,3 x
7f0,3,
wl(2,2)
=
l2,2 x
7f2,2,
wl(4,1)
=
l4,1 x
7f4,1,
wl(6,0)
=
l6,0 x
7f6,0.
wl(1,3)
=
l1,3 x
7f1,3,
wl(3,2)
=
l3,2 x
7f3,2,
wl(5,1)
=
l5,1 x
7f5,1,
wl(1,2)
=
l1,2 x
7f1,2,
wl(3,1)
=
l3,1 x
7f3,1,
wl(5,0)
=
l5,0 x
7f5,0,
The lengths of intersection of the rays with image elements equal l0,3 = l2,2 = l4,1 = l6,0 = ( 5/7)/2, and l1,3 = l1,2 = l3,2 = l3,1 = l5,1 = l5,0 = ( 5/7)/4. Therefore, the line-integral wl(6) can be written
as wl(6) =
wl(n,m)
l( n,m)=
= [wl(0,3) + wl(2,2) + wl(4,1) + wl(6,0)] +[wl(1,3) + wl(3,2) + wl(5,1)]+ [wl(1,2) + wl(3,1) + wl(5,0)]
=
7
5 2
[f0,3
+
f2,2
+
f4,1
+
f6,0 ]
+
7 4 5 [f1,3
+
f3,2
+
f5,1]
+
7 4 5 [f1,2
+
f3,1
+
f5,0 ]
=
75 2
v1,2(6) +
1 2
v1,2
(5)
+
1 2
v1,2(7)
.
The similar equation holds for other line-integrals w1,2(t), as well. For instance,

wl(7) =
75 2
v1,2(7)
+
1 2
v1,2(6)
+
1 2
v1,2(8)
.
Different coefficients at the line-sums v1,2(t), v1,2(t - 1), and v1,2(t + 1) complicate the solution of the system of these equations with respect to v1,2(t). To simplify the calculation of the sums v1,2(t), we describe the system of equations, wherein all coefficients at the sums are equal. For that, we first consider the geometrical ray ~l which is located between the rays l(6) and l(7). This ray is denoted by ~l(7) in the figure. One can notice that this ray intersects equally all image elements with numbers (0, 3), (1, 3), (2, 2), (3, 2), (4, 1), (5, 1), and (6, 0). Therefore the integral of the image along the geometrical ray ~l can be calculated from the sums of the discrete image along two arithmetical rays by
w~l =
w~l(n,m) ~l( n,m)=
= wl(0,3) + wl(1,3) + wl(2,2) + wl(3,2) + wl(4,1) + wl(5,1) + wl(6,0)
= [wl(0,3) + wl(2,2) + wl(4,1) + wl(6,0)] +[wl(1,3) + wl(3,2) + wl(5,1)]
=
75 2
[f0,3
+
f2,2
+
f4,1
+
f6,0]
+
7 2 5 [f1,3
+
f3,2
+
f5,1]
=
75 2
v1,2(6) + v1,2(7)
.
We obtain a simple equation for the line-integral along the geometrical ray ~l. This ray is considered as the ray number 7 in the new set of parallel geometrical rays that are defined as
~l(t) = l1,2(t - 0.5) =
(x,
y);
x
+
2y
=
t 7
+
1 7
,
t = 0 : 18,
and illustrated in Figure 10.
y and m
0 0 0.14 1 0.28 2 0.42 3 0.57 4 0.71 5 0.85 6 1
~l(0) ~l(1) ~l(2) ~l(3) ~l(4) ~l(5)
~l(6) ~l(7) ~l(8) ~l(9) ~l(10) ~l(11) ~l(12) ~l(13) ~l(14) ~l(15) ~l(16) ~l(17) ~l(18)
0
1
2
3
4
5
6
0
0.14 0.28 0.42 0.57 0.71 0.85 1
x and n
Figure 10. The set of shifted geometrical rays for calculating 19 line-integrals of the (1, 2)-projection.
The relation between the integrals along the geometrical rays ~l(t) and sums along the arithmetical rays l(t) is described by the following system of linear equations: w1,2(t) = w(~l(t)) = 7 2 5 [v1,2(t - 1) + v1,2(t)], t = 0 : 18,
where v1,2(-1) = 0. These equations can be written in matrix form as
w1,2(0)
1 0 0 0 0 · · · 0 0 v1,2(0)
w1,2 (1)
1 1 0 0 0 ··· 0 0
v1,2 (1)



w
=




w1,2 (2) w1,2 (3) w1,2 (4) ...




=



75 2 Av
=
75 2

0 0 0 ...
1 0 0 ...
1 1 0 ...
0 1 1 ...
0 0 1 ...
··· ··· ··· ...
0 0 0 ...
0

v1,2 (2)

0
0

v1,2 (3) v1,2 (4)


.
...

...

w1,2(17)
0 0 0 0 0 · · · 1 0 v1,2(17)
w1,2(18)
0 0 0 0 0 ··· 1 1
v1,2(18)
The above Toeplitz matrix A of size 19 Ч 19 has the triangle inverse matrix, and all sums along the arithmetical rays are calculated by
v1,2(0)
1 0 0 0 0 · · · 0 0 w1,2(0)
v1,2 (1)

v1,2 (2)

-1 1 0 0 0 · · · 0 0
w1,2(1)

1 -1
1
0
0 ···
0
0

w1,2(2)


v1,2 (3) v1,2 (4) ...


=
2
-1 1

7 5
...
1 -1 ...
-1 1 ...
1 -1 ...
0 ··· 1 ··· ... . ..
0 0 w1,2(3)
0
0
w1,2(4)

.
...
...

...

v1,2(17)
-1 1 -1 1 -1 · · · 1 0 w1,2(17)
v1,2 (18)
1 -1 1 -1 1 · · · -1 1
w1,2 (18)
The
required
sums
v(t)
of
the
discrete
image fn,m
can
be
calculated
directly
by
v
=
2 75
A-1
w,
or
by
the
method
of direct substitution in the following recursive form:
b0 = w(0),
b1 = w(1) - b0,
(14)
bt = w(t) - bt-1, t = 2, 3, . . . , 18,
and, then, v(t) = (2/7 5)bt, where t = 0 : 18. Substituting these values in (13), we obtain all components of the splitting-signal {f1,2,t; t = 0 : 7}.
Example 2 (7 Ч 7). Consider the N = 7 case and (p, s) = (1, 5). The masks of matrices of the tensor functions 1,5,t(n, m) with the first two triplet-numbers of the subset {(1, 5, t); t = 0 : 6} are
1 0 0 0 0 0 0
0010000

0
0
0
0
1
0
0

[1,5,0]
=

0 0
0 1
0 0
0 0
0 0
0 1
1 ,
0

0 0 0 1 0 0 0
0000010
0 1 0 0 0 0 0
0001000

0
0
0
0
0
1
0

[1,5,1]
=

1 0
0 0
0 1
0 0
0 0
0 0
0 .
0

0 0 0 0 1 0 0
0000001
The direction of parallel rays in these masks are defined by the angle arctan(5), as shown in Figure 11 in part a for the first mask. There are six such parallel rays that pass through the 1's in the mask. Instead of parallel rays of the (1, 5)-projection, we consider that (1, -2)-projection which also defines the same set of masks of the tensor functions. Indeed, the equation 1n + 5m = t mod 7 can be written as 1n + (7 - 2)m = t mod 7, or 1n - 2m = t mod 7. The same mask with two parallel rays of the (1, -2)-projection is illustrated in part b.
y and m y and m
0 0 0.14 1 0.28 2 0.42 3 0.57 4 0.71 5 0.85 6 1 0
l(0) l(7) l(14) l(21) l(28) l(35) 0123456 0.14 0.28 0.42 0.57 0.71 0.85 1 x and n
0 0 0.14 1 0.28 2 0.42 3 0.57 4 0.71 5 0.85 6 1 0
l(0) l(-7) 0123456 0.14 0.28 0.42 0.57 0.71 0.85 1 x and n
Figure 11. The mask [1,5,0] with two different sets of parallel rays.
The set of all rays for the (1, -2)-projection is shown in Figure 12 together with the control points numbering these rays. Such change of projection allows for reducing the number of rays from (1 + 5)6 + 1 = 37 to (1 + | - 2|)6 + 1 = 19. The set of 19 arithmetical rays of this projection is described as
l(t) = l1,-2(t) =
(x,
y);
x
-
2y
=
t 7
-
1 14
,
t = 6 : -1 : -12,
(15)
where (x, y) [0, 1] Ч [0, 1]. The masks for the tensor functions 1,5,0(n, m) and 1,5,1(n, m) are shown in Figure 13. It is not difficult to see, that the components of the tensor transform f1,5,t are calculated by
f1,5,0 = v0 + v-7, f1,5,1 = v1 + v-6, f1,5,2 = v2 + v-5 + v-12, f1,5,3 = v3 + v-4 + v-11, f1,5,4 = v4 + v-3 + v-10, f1,5,5 = v5 + v-2 + v-9, f1,5,6 = v6 + v-1 + v-8,
(16)
where the sums of the discrete image along the rays are denoted by v(t) = v1,-2(t), t = 6 : -1 : -12. The above system of linear equations can also be written as follows:
2
2
f1,5,t =
vt-7m =
v1,-2,t-7m,
m=0
m=0
t = 0 : 6,
y and m
0 0 0.14 1 0.28 2 0.42 3 0.57 4 0.71 5 0.85 6 1 0
l(6)
l(5)
l(4)
l(3)
l(2)
l(1)
l(0)
l(-1)
l(-2)
l(-3)
l(-4)
l(-5)
l(-6)
l(-12) l(-11) l(-10) l(-9) l(-8) l(-7)
1
2
3
4
5
6
0· 1· 2· 3· 4· 5· 6·
-2· -1·

-4·
-3·






-6· -5·

-8·
-7·






-10· -9·
-12· -11·
(for rays at angle 26.56o).
0
0.14 0.28 0.42 0.57 0.71 0.85 1
x and n
Figure 12. The set of 19 arithmetical rays for the (1, -2)-projection and control points.
y and m
0 0 0.14 1 0.28 2 0.42 3 0.57 4 0.71 5 0.85 6 1 0 0
1234
0.14
0.28
x and n
l(0) l(-7) 56 0.42
y and m
0 0 0.14 1 0.28 2 0.42 3 0.57 4 0.71 5 0.85 6 1 0 0
1234
0.14
0.28
x and n
l(1) l(-6) 56 0.42
Figure 13. Two masks with parallel rays for the (1, -2)-projection.
where vt = 0, if t < -7. We now consider the shifted set of the rays, l(t) l(t - 0.5), as the set of 19 geometrical
rays
~l(t) = l1,-2(t - 0.5) =
(x, y); x
-
2y
=
t 7
-
1 7
,
t = 6 : -1 : -12.
The set of these G-rays is shown in Figure 14. For this set of G-rays, the system of equations for the line-integrals
is defined as

w(t) = w(~l(t)) =
75 2
v(t) + v(t - 1)
,
t = 6 : -1 : -12,
(17)
where v(-13) = 0. This system can be written in matrix form with the same Toeplitz matrix described for the (1, 2)-projection. Therefore, the sums v(t) of the discrete image fn,m can be determined by v(t) = bt/(7 5/2), where the components bt are calculated recursively by
b-12 = w(-12), b-11 = w(-11) - b-12, b-10 = w(-10) - b-11, bt = w(t) - bt-1, t = -9, -8, -7, -6, · · ·, 6.
(18)
It should be noted that the numbering of the control points can be changed by the transform t 6 - t. The set of 19 geometrical rays of this projection in (15) will then be written as
l(t) = l1,-2(6 - t) =
(x,
y);
x
-
2y
=
-
1 7
+
11 14
,
0
0
0.14
l(6)
1
l(5)
y and m
0.28 2 0.42 3 0.57 4 0.71 5 0.85 6 1
l(-12) l(-11) l(-10) l(-9) l(-8) l(-7)
0
1
2
3
4
5
6
l(4) l(3) l(2) l(1) ~l(0) ~l(-1) ~l(-2) ~l(-3) ~l(-4) ~l(-5) ~l(-6)
0
0.14 0.28 0.42 0.57 0.71 0.85 1
x and n
Figure 14. The set of geometrical rays for the (1, -2)-projection.
and
the
system
of
equation
in
(17)
as
w(t)
=
w(~l(t))
=
7 5/2
v(t)
+ v(t
+ 1)
,
t
=
0
:
18,
where
v(19)
=
0.
The
equation (18) for the coefficients bt = 7 5/2v(t) is considered as
b18 = w(18), b17 = w(17) - b18, b16 = w(16) - b17, bt = w(t) - bt+1, t = 15, 14, · · · , 1, 0.
Therefore, the components of the splitting signal f1,5,t in equation 16 can be calculated as follows:
2
f1,5,t =
v6-t+7m,
m=0
t = 0 : 6,
(vt = 0, t > 18).
4.1 Main Equations for G-rays when N is a prime In general case, when N is prime, the image reconstruction from N + 1 projections is described similar to the N = 7 case. Table 1 provides the general equations of rays and calculating the sums v(t) from the line-integrals w(t), for all projections. Two subsets of generators (p, s) are considered separately. The first subset contains generators (1, s), where s (N - 1)/2. The generators (1, s), where s > (N - 1)/2, compose the second subset. For each generators (1, s) of the second subset, the (1, s)-projection is substituted by the (1, s - N )-projection. In other words, the splitting-signal {f1,s,t; t = 0 : (N - 1)} is calculated from the line-integrals of the (1, s - N )projection, similar to the (1, 5)-projection described above for the N = 7 case. Each such "substitution" saves the number of line-integrals used in the reconstruction. The number of parallel G-rays used in each projection is defined as
n(1, s) = (1 + s)(N - 1) + 1 = N + s(N - 1), s = 0 : (N - 1)/2,
n(1, s) = n(1, N - s) = N 2 - s(N - 1), s = (N + 1)/2 : N - 1,
(19)
n(0, 1) = n(1, 0) = N.
l(t) w(t)
Equations
xp + ys
=
t N
+
1 N
,
where
s
N -1 2
xp + y(s - N ) =
t N
+
s+1-N N
,
where
s>
N -1 2
w(t) = N v(t), where s = 0, or p = 0,
t 0 : (p + s)(N - 1) - 1 N - 1 : -1 : (s - N )(N - 1) 0:N -1
b(t) v(t) fp,s,t
w(t) = K[v(t) + v(t - 1) + · · · + v(t - (s - 1))]
where
s

N -1 2
and
K
=
N s
p2 + s2
w(t) = K[v(t) + v(t - 1) + · · · + v(t - (N - s - 1))]
where
s
>
N
-
12
and
K
=
N N -s
p2 + (N - s)2
b(0) = w(0), bt = w(t) - (bt + · · · + bt-(s-1))
where
s

N -1 2
b(0) = w(0), bt = w(t) - (bt + · · · + bt-(N-s-1))
where
s
>
N -1 2
v(t)
=
bt K
where
s

N -1 2
and
K
=
N s
p2 + s2
v(t)
=
bt K
where
s
>
N -1 2
and
K
=
N N -s
p2 + (N - s)2
m vt+7m,
where
s
N -1 2
m vt-7m,
where
s>
N -1 2
0 : (p + s)(N - 1) - 1 -(N - s)(N - 1) : N - 1 0 : (p + s)(N - 1) - 1 -(N - s)(N - 1) : N - 1 0 : (p + s)(N - 1) - 1 (s - N )(N - 1) : N - 1 0:N -1 0:N -1
Table 1. General equations of rays l(t), linear integrals w(t), sums v(t), and splitting-signals fp,s,t.
4.2 Simulation results on Modeled Images The proposed method of image reconstruction was implemented in MATLAB and C+. For N prime, the blockdiagram of the transform-based method of reconstruction of the image on the Cartesian grid N Ч N from N + 1 projections is shown in Figure 15. The following steps describe the proposed method of image reconstruction. 1. The image N Ч N is composed (for instance, from a few random rectangles or ellipses) on the square [0, 1] Ч [0, 1]. 2. The set JN,N of generators (p, s) is calculated to define N + 1 projections. The considered set is JN,N = {(1, s); s = 0 : (N - 1)} {(0, 1)}. 3. Given frequency-point (p, s), the (p, s)-projection is calculated, i.e., all line-integrals w(t) along the set of geometrical rays ~lp,s(t). 4. The geometry of G-rays of the (p, s)-projection is transformed to A-rays, i.e., the set of sums v(t) of the discrete image fn,m is calculated from the integrals w(t), by solving the system of linear equations with the Toeplitz matrix (M ЧM ), where M = (p+s-1)(N -1)+1, if s (N -1)/2, and M = (p+(N -s)-1)(N -1)+1, if s > (N - 1)/2. The Toeplitz matrix is triangular and the transform of geometry of G-rays to A-rays is fast. 5. The splitting-signal {fp,s,0, fp,s,1, ..., fp,s,N-1} is calculated from v(t). This signal is written into the corresponding row of the matrix N Ч (N + 1) of the 2-D discrete tensor transform. 6. The inverse 2-D tensor transform is calculated, which is the reconstructed discrete image fn,m.
image fd(x, y) Image composition on the unit square 1 ? Line-integrals of the (p, s)-projection 3 ? Transform of G-rays into A-rays, w v 4
2 Set of generators JN,N = {(p, s)}
? Calculation of the splitting-signal {fp,s,t} 5 ? Inverse 2-D discrete tensor transform 6
? reconstructed image fn,m
Figure 15. Block-diagram of the reconstruction of the image f (x, y) = fd(x, y) composed of image elements with constant intensity each on the Cartesian grid N Ч N placed in the square [0, 1] Ч [0, 1].
11 rectangles on the grid 131x131 0
reconstructed image 131x131
20 0.2 40 0.4 60
0.6
80
0.8
100
120
1
0
0.5
1
20 40 60 80 100 120
(a)
(b)
Figure 16. (a) Image with 13 rectangles on the lattice 131 Ч 131 and (b) its reconstruction.
Figure 16 shows the image f(x, y) on the square [0, 1]Ч[0, 1], which is generated by the random rectangles with coordinates on the lattice 131 Ч 131. The reconstruction fn,m of the image by 132 projections is also shown. One can see from the reconstructed image that the reconstruction is exact. The experiments were simulated on the computer with Intel(R) Core(TM) 2Duo E8400 3.00GHz and 3.21GB of RAM. The time for image reconstruction is 32.16 seconds. For the (1, 4)-projection, Figure 17 shows the line-integrals w(t) on the left, the sums v(t) in the middle, and the splitting-signal {f1,4,0, f1,4,1, ..., f1,4,130} on the right. We now consider the 257 Ч 257 example, when 258 projections are required for the image reconstruction. Figure 18 presents the image on the left with the exact reconstruction on the right. Table 2 shows the total time of image processing, when the program was implemented in Ci+. The tested image f(x, y) was composed of 15 random rectangles on the square [0, 1]Ч[0, 1]. The time includes the calculation of all line-integrals w(t), line-sums v(t), and the calculation of the direct and inverse tensor transforms. 5. GEOMETRY FOR THE LATTICE N Ч N , WHEN N IS A POWER OF 2 When reconstructing the image f(x, y) from its projections on the Cartesian lattice of size N Ч N, when N = 2r, r > 1, is a power of two, we can define components of the tensor transform of the image from the line-integrals in a way, similar to the described above case when N is a prime. The effective reconstruction is with the 2-D
6000 5000 4000 3000
w(t), t =0 : 650
1200 1000 800 600
v(t), t =0 : 650
2000
400
f1,4,t , t =0 : 130 26 24 22 20 18
1000
200
16
0
0
14
0
500
0
500
0 50 100
t
t
t
Figure 17. Line-integrals w(t), ray-sums v(t), and splitting-signal for the (1, 4)-projection.
10 rectangles on the grid 257x257 0
reconstructed image 257x257
0.2
50
0.4
100
0.6
150
0.8
200
1
250
0
0.5
1
50 100 150 200 250
(a)
(b)
Figure 18. Image with ten rectangles on the lattice 257 Ч 257 and its reconstruction.
N Ч N projection data processing & reconstruction
67 Ч 67
00 : 00 : 00.08
127 Ч 127
00 : 00 : 00.50
257 Ч 257
00 : 00 : 04.21
521 Ч 521
00 : 00 : 37.24
1031 Ч 1031
00 : 04 : 55.69
2053 Ч 2053
00 : 39 : 40.44
Table 2. Time for the Ci-based program (the 1st version)
paired transform, because it removes the redundancy of the tensor transform. Therefore we consider the image in paired representation and describe the method of calculation of the paired splitting-signals from line-integrals. Given generator (p, s) JN,N , we define the integer k by 2k = g.c.d.(p, s), and k = r - 1 when (p, s) = (0, 0). The components of the 2-D paired transform
fp,s,2kt = p,s,2kt f = fp,s,2kt - fp,s,2kt+N/2, t = 0 : (N/2k+1 - 1), are calculated by the system of orthogonal paired functions which are defined as
1, if np + ms = 2kt mod N p,s,2kt(n, m) = -1, if np + ms = (2kt + N/2) mod N 0, otherwise
where (n, m) XN,N . The set of N 2 triplet-numbers of the paired functions is taken as
r-1
UN,N =
(p, s, 2kt); (p, s) 2kJN/2k,N/2k , t = 0 : (N/2k+1 - 1) {(0, 0, 0)}.
k=0
To describe the method of the transferring geometry from the image plane to the Cartesian lattice, we describe the N = 8 case.
Example 3 (8 Ч 8). We consider the N = 8 case and (p, s) = (4, 1). Four masks of the paired functions 1,4,t(n, m) with the triplet-numbers (1, 4, t), t = 0 : 3, are
1 0 0 0 -1 0 0 0
0 1 0 0 0 -1 0 0
-1 0 0 0 1 0 0 0
0 -1 0 0 0 1 0 0

1
0
0
0 -1
0
0
0


0
1
0
0
0 -1
0
0


[1,4,0]
=


-1 1
0 0
0 0
01 0 -1
0 0
0 0
0 0


,


[1,4,1]
=


0 0
-1 1
0 0
0 0
01 0 -1
0 0
0 0


,


-1
0
0
0
1
0
0
0


0
-1
0
0
0
1
0
0


1
0
0
0 -1
0
0
0


0
1
0
0
0 -1
0
0

-1 0 0 0 1 0 0 0
0 -1 0 0 0 1 0 0
0 0 1 0 0 0 -1 0
0 0 0 1 0 0 0 -1
0 0 -1 0 0 0 1 0
0 0 0 -1 0 0 0 1

0
0
1
0
0
0
-1
0


0
0
0
1
0
0
0
-1


[1,4,2]
=


0 0
0 0
-1 1
0 0
0 0
0 0
1 -1
0 0


,


[1,4,3]
=


0 0
0 0
0 -1 01
0 0
0 0
0 0
1 -1


.


0
0
-1
0
0
0
1
0


0
0
0 -1
0
0
0
1


0
0
1
0
0
0
-1
0


0
0
0
1
0
0
0
-1

0 0 -1 0 0 0 1 0
0 0 0 -1 0 0 0 1
The set of the 36 parallel rays of this projection is shown in Figure 19,
l(t) = l1,4(t) =
(x,
y);
x
+
4y
=
t 8
+
5 16
,
t = 0 : 35,
(x, y) [0, 1] Ч [0, 1].
y and m
0 0 0.25 1 0.25 2 0.375 3 0.5 4 0.625 5 0.75 6 0.875 7 1 0
l(0) l(1) l(2) l(3) l(4) l(5) l(6) l(7) l(8) l(9) l(10) l(11) l(12) l(13) l(14) l(15) l(16) l(17) l(18) l(19) l(20) l(21) l(22) l(23) l(24) l(25) l(26) l(27) l(28) l(29) l(30) l(31) l(32) l(33) l(34) l(35)
0
1
2
3
4
5
6
7
0.125 0.25 0.375
0.5
0.625 0.75 0.875
1
x and n
Figure 19. The sets of parallel rays for the (1, 4)-projection.
(1 + 4)7 + 1 = 36 parallel rays are considered for the (1, 4)-projection. The control points of the set of rays
are defined as shown below:
0· 1· 2· 3· 4· 5· 6· 7·
. . . . 8· 9· 10· 11·

.
.
.
.
12·
13·
14·
15·


.

.
. .
. .
. .
16· 20·
17· 21·
18· 22·
19·

23·

(angle of rays is - tan-1(1/4) = -14.036).

.
.
.
.
24·
25·
26·
27·


.
.
.
.
28·
29·
30·
31·

. . . . 32· 33· 34· 35·
The rays are shown separately on the masks of the paired functions 1,4,t(n, m), t = 0 : 3, in Figure 20.
y and m
0
l(0)
0
l(4) l(8)
0.25 1
l(12)
0.25 2
l(16)
0.375 3
l(20)
0.5 4
l(24)
0.625 5
l(28)
0.75 6
l(32)
0.875
7
1 01234567
0 0.1250.250.3750.50.6250.750.875 1
x and n
y and m
0
l(1)
l(5)
0
l(9)
0.25
1
l(13)
0.25
2
l(17)
0.375
3
l(21)
0.5
4
l(25)
0.625
5
l(29)
0.75
6
l(33)
0.875
7
1 01234567
0 0.1250.250.3750.50.6250.750.875 1
x and n
y and m
0
l(2)
l(6)
0 0.25
l(10)
1 0.25
l(14)
2 0.375
l(18)
3 0.5
l(22)
4 0.625
l(26)
5 0.75
l(30)
6 0.875
l(34)
7
1 01234567
0 0.1250.250.3750.50.6250.750.875 1
x and n
y and m
0
l(3)
l(7)
0
0.25
l(11)
1
0.25
l(15)
2
0.375
l(19)
3
0.5
l(23)
4
0.625
l(27)
5
0.75
l(31)
6
0.875
l(35)
7
1 01234567
0 0.1250.250.3750.50.6250.750.875 1
x and n
Figure 20. Four masks with the set of arithmetical rays for the paired functions with triplet-numbers (1, 4, t), t = 0 : 3.
We denote by vt = v1,4(t), where t = 0 : 35, the sums of the discrete image fn,m along these rays. The components of the paired transform with the triplet-numbers (1, 4, t) can be calculated as follows:
f1,4,0 = v0 - v4 + v8 - v12 + v16 - v20 + v24 - v28 + v32, f1,4,1 = v1 - v5 + v9 - v13 + v17 - v21 + v25 - v29 + v33, f1,4,2 = v2 - v6 + v10 - v14 + v18 - v22 + v26 - v30 + v34, f1,4,3 = v3 - v7 + v11 - v15 + v19 - v23 + v27 - v31 + v35.
(20)
As follows from (10), from the ray-sums of the (1, 4)-projection, we can also calculate the components f2,0,0, f2,0,0, f4,0,0, and f0,0,0. Two masks with rays for the paired functions 2,0,t, t = 0, 2, are shown in Figure 21. The corresponding components of the 2-D paired transform are calculated by
y and m y and m
0 0 0.25 1 0.25 2 0.375 3 0.5 4 0.625 5 0.75 6 0.875
l(0)
l(2)
l(4) l(6)
l(8)
l(10)
l(12)
l(14)
l(16)
l(18)
l(20)
l(22)
l(24)
l(26)
l(28)
l(30)
l(32)
l(34)
7
1 01234567
0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1 x and n
l(1)
l(3)
l(5)
0
l(7)
0
l(9)
0.25
l(11)
1
l(13)
0.25
l(15)
2
l(17)
0.375
l(19)
3
l(21)
0.5
l(23)
4
l(25)
0.625
l(27)
5
l(29)
0.75
l(31)
6
l(33)
0.875
l(35)
7
1 01234567
0 0.125 0.25 0.375 0.5 0.625 0.75 0.875 1 x and n
Figure 21. Two masks with the set of arithmetical rays for the paired functions with triplet-numbers (2, 0, t), t = 0, 1.
f2,0,0 = (v0 + v4 + v8 + v12 + v16 + v20 + v24 + v28 + v32) -(v2 + v6 + v10 + v14 + v18 + v22 + v26 + v30 + v34), f2,0,2 = (v1 + v5 + v9 + v13 + v17 + v21 + v25 + v29 + v33) -(v3 + v7 + v11 + v15 + v19 + v23 + v27 + v31 + v35).
(21)
Similarly, for the remaining two components, we have the following: f4,0,0 = v0 - v1 + v2 - v3 + · · · + v34 - v35 and f0,0,0 = v0 + v1 + v2 + v3 + · · · + v34 + v35. We consider the following set of triplet-numbers:
U (1, 4) = {(1, 4, 0), (1, 4, 1), (1, 4, 2), (1, 4, 3), (2, 0, 0), (2, 0, 2), (4, 0, 0), (0, 0, 0)}
and define by P(1, 4) = {fp1,s1,t; (p1, s1, t) U (1, 4)} the components of the 2-D paired transforms with these triplets. Introducing the new variables
4
4
an =
vn+8m =
v1,4,n+8m,
m=0
m=0
n = 0 : 7,
where vk = 0, if k > 35, the complete system of linear equations for the eight components fp1,s1,t P(1, 4) can
be written as follows:
f1,4,0 = a0 - a4
f1,4,1 = a1 - a5
f1,4,2 = a2 - a6
f1,4,3 = a3 - a7 f2,0,0 = (a0 + a4) - (a2 + a6)
(22)
f2,0,2 = (a1 + a5) - (a3 + a7)
f4,0,0 = (a0 + a4) - (a1 + a5) + (a2 + a6) - (a3 + a7)
f0,0,0 = (a0 + a4) + (a1 + a5) + (a2 + a6) + (a3 + a7).
These calculations correspond to the fast 8-point discrete paired transform (DPT)33,34,50 of the vector a =
(a0, a1, a2, ... , a7) , i.e., P(1, 4) = [8]a.
To express the sums v(t) of the discrete image fn,m through the line-integrals w(t) = w1,4(t), we consider the part of the mask of the paired function with number (1, 4, 0), which is shown in Figure 22. The geometrical ray
l(29) between two points of the discrete image, (0, 7) and (5, 6) passes through six image elements of numbers
(0, 7), (1, 7), (2, 7), (3, 7), (3, 6), and (4, 6). The length of the intersection l0,7 of the ray with the IE of number (0, 7) is twice the length of intersection with each IE of number (3, 7) and (3, 6), but equal to the length of intersection with the IE of numbers (1, 7) and (2, 7). Since l0,7 = 17/4x, we can write the following:
w(29)
=
w1,4(29)
=
l0,7 (x)2
v(28)
+
v(29)
+
v(30)
+
1 2
v(31)
+
1 2
v(27)
.
0.625 5 0.75 6
l(27) l(28) l(29) l(30) l(31) l(32)
y and m
0.875
7
l1 1 l2
0
1
2
3
4
5
6
7
0
0.125 0.25 0.375 0.5 0.625 0.75 0.875
1
x and n
Figure 22. Set of geometrical rays for calculating line-integrals when (p, s) = (1, 4).
Similar equations can be used for calculating other integrals w(t), t = 0 : 35, and then the sums v(t) can be
found. Now we consider the geometrical rays l1 and l2 between rays 28, 29, and 29, 30, respectively. Thelength of the intersection of ray l1 with the image elements of numbers (0, 7), (1, 7), (2, 7), and (3, 6) is the same, 17/4x.
The length of the intersection of ray l2 with each image element of number (0, 7), (1, 7), (2, 7), and (3, 7) is also 17/4x. Therefore, the integrals w1 and w2 of the image along the shifted rays l1 and l2, respectively, can be
written as follows:
w1
=
17/4x (x)2
[v(28)
+
v(29)
+
v(30)
+
v(27)]
=
2 17
[v(27)
+
v(28)
+
v(29)
+
v(30)]
,

w2 = 2 17 [v(28) + v(29) + v(30) + v(31)] .
These equations can be generalized if we denote the integral w1 by w(30), and w2 by w(31). (We can also denote w1 by w(27) and w2 by w(28).) In this case, a simple solution for sums v(t) can be obtained. Indeed, we consider a new set of parallel rays defined by the shifting, t t - 3/2,
~l(t) = l1,4(t - 3/2) =
(x, y);
x
+
4y
=
t 8
+
1 8
,
t = 0 : 35,
which is shown in Figure 23. For the integrals along the set of geometrical rays, we can write the following equations: w(t) = 2 17 [v(t) + v(t - 1) + v(t - 2) + v(t - 3)] , t = 0 : 35,
y and m
0 0 0.25 1 0.25 2 0.375 3 0.5 4 0.625 5 0.75 6 0.875
l(0) l(1) l(2) l(3) l(4) l(5) l(6) l(7) l(9) l(11) l(13) l(15) l(17) l(19) l(21) l(23) l(25) l(27) l(29) l(31) l(33) l(35)
7
1
0
1
2
3
4
5
6
7
0
0.125 0.25 0.375 0.5 0.625 0.75 0.875
1
x and n
Figure 23. The set of geometrical rays for the (1, 4)-projection.
where v(-1) = v(-2) = v(-3) = 0. This system of linear equations can be written in matrix form as
1 0 0 0 0 0 ... 0 0 0 0 0
1 1 0 0 0 0 ... 0 0 0 0 0

1
1
1
0
0
0
...
0
0
0
0
0


1
1
1
1
0
0
...
0
0
0
0
0

w
=
2 17Av
=
2 17


0 0
1 0
1 1
1 1
1 1
0 1
... ...
0 0
0 0
0 0
0 0
0 0
v,
(23)

.
.
.
.
.
. ...
.
.
.
.
.


.
.
.
.
.
. ...
.
.
.
0
0


0
0
0
0
0
0
...
1
1
1
1
0

0 0 0 0 0 0 ... 0 1 1 1 1
where w = (w0, w1, ..., w35) and v = (v0, v1, ..., v35) . The above Toeplitz matrix 36 Ч 36 has the following inverse lower triangle Toeplitz matrix:
1 0 0 0 0 0 0 ... 0 0 0 0
-1 1 0 0 0 0 0 ... 0 0 0 0

0 -1
1
0
0
0
0 ...
0
0
0
0


0
0 -1
1
0
0
0 ...
0
0
0
0


1
0
0 -1
1
0
0 ...
0
0
0
0


-1
1
0
0 -1
1
0 ...
0
0
0
0

A-1
=

0 -1
1
0
0 -1
1 ...
0
0
0
0

.
(24)

0
0 -1
1
0
0 -1 ...
0
0
0
0


...
...
...
...
...
...
... ...
...
...
...
...


1
0
0 -1
1
0
0 ...
1
0
0
0


-1
1
0
0 -1
1
0 ... -1
1
0
0


0 -1
1
0
0 -1
1 ...
0 -1
1
0

0 0 -1 1 0 0 -1 ... 0 0 -1 1
From
the
solution
of
the
equation
b
=
A-1w,
we
obtain
the
sums
v
=
b/(2 17).
The
vector
b
can
be
calculated by using the following recurrent procedure:
b0 = w(0),
b1 = w(1) - b0,
b2 = w(2) - (b1 + b0), b3 = w(3) - (b2 + b1 + b0),
(25)
b4 = w(4) - (b3 + b2 + b1),
bt = w(t) - (bt-1 + bt-2 + bt-3), t = 5, 6, . . . , 35.
Similar calculations hold for other projections, and the ray-sums can be calculated from the line-integrals by solving the system of linear equations described by the triangle Toeplitz matrices (see45 for more detail).
5.1 Algorithm of image reconstruction For the N = 8 case, all 3N/2 = 12 projections that are required for reconstructing the image N Ч N The number of parallel geometrical rays used in the (p, s)-projections is defined by (19). In the case when N is a power of two, the line-integrals are used to calculate the splitting-signals in the paired or tensor representation of the discrete image. Therefore, the block-diagram given in (15) can be used for the reconstruction by the tensor transform. In order to reduce the complexity of calculation which are due to redundancy of the tensor transform, the 2-D paired transform can be calculated, as shown in the above N = 8 example, and then, the reconstruction fn,m is calculated by the inverse paired transform. Each triplet-number (p, s, t) of masks of the paired functions in this table has the form 2k(p1, s1, t1), where 2k=g.c.d.(p, s), and t1 = 0 : N/2k+1 - 1, and its multiplicity (i.e., the number of subsets U covering this triplet-number) equals 2k. Therefore, there are two ways to use this property in image reconstruction. 1. One can use the complete 1-D paired transforms over all vectors a, whose components are defined as
an = vn±Nm = vp,s,n±Nm, n = 0 : N - 1.
m
m
Then, the components fp,s,t of the 2-D paired transform, which have been calculated by 1-D paired transform P(p, s) = [N ]a, are normalized by the factor of 2k = g.c.d.(p, s), when (p, s) = (0, 0), and N, when (p, s) = (0, 0). 2. One can use the incomplete complete 1-D paired transforms over vectors a, to avoid the repetition of calculation of components of the 2-D paired transform of the image. For instance, when N = 8, the complete 1-D paired transform can be used for the first (1, 0)-projection, but only six outputs for the transform when calculating the transform P(1, 2) = [8]a, for the (1, 2)-projection, and half of calculations when calculating the transform P(1, 4) = [8]a, for the (1, 4)-projection, and so on. The use of incomplete 1-D paired transforms reduces the time necessary to calculate the reconstructed image, when compared with the first method.
In the general N = 2r case where r 2, the system of equations for the G-rays is defined in the following way. We first consider the A-rays l(t) = lp,s(t) which are described by the equation
lp,s,t(x, y) =
(x, y);
xp
+
ys
=
t N
+
p+s 2N
.
To reduce the number of line-integrals, we can divide the set of all generators of JN,N by four subsets. The first and second subsets contain all generators (1, s), where s = 0 : N/2 and s = N/2 + 1 : N - 1, respectively. The third and fourth subsets contain all generators (p, 1), where p = 2p1 and p1 = 0 : N/4 and p1 = N/4 +1 : N/2 -1, respectively. For each of these subsets, we can consider the case when the arithmetical rays are shifted to the right, to define the corresponding geometrical rays. Another case is when the geometrical rays are defined from arithmetical rays by shifting to the left. Thus, the geometrical rays are defined from the arithmetical rays by shifting ~l(t) = l(t - t0) (or ~l(t) = l(t + t0). Here the shift t0 is calculated by
t0 =
p
+ 2
s
-
1,
p
+
(N 2
-
s)
-
1,
(N - p) + s - 1,
when 0 < p, s N/2, when p = 1, s > N/2, when p > N/2, s = 1.
(26)
2
5.2 Convolution equations
For all generators (p, s) JN,N , the relations between the sums v(t) = vp,s(t) of the discrete image fn,m and line-integrals w(t) = wp,s(t) of the original unknown image f(x, y) are described in matrix form by the Toeplitz matrices, which are similar to the described N = 8 example. We consider separately the convolution equations of line-integrals for different sets of (p, s)-projections.
Subset I: p = 1 and s N/2.
When s = 0, the convolution equation w = Av is described as

w(t) = N
1
+ s
s2
[v(t)
+
v(t
-
1)
+
v(t
-
2)
+
...
+
v(t
-
s
+
1)],
t = 0 : (1 + s)(N - 1),
(27)
where it is assumed that v(-1) = v(-2) = ... = v(-s + 1) = 0. The inverse transform b = A-1w is calculated by the following recurrent form:
b0 = w(0),
b1 = w(1) - b0,
b2 = w(2) - (b1 + b0),
b3 = w(3) - (b2 + b1 + b0),
(28)
... ...
... ... ... ...
bs-1 = w(s - 1) - (bs-2 + bs-3 + ... + b1 + b0),
bt = w(t) - (bt-1 + bt-2 + ... + bt-s+2 + bt-s+1)), t = s, s + 1, . . . , (1 + s)(N - 1).
The required sums v(t) of the discrete image are calculated by
v(t)
=
bt
N
s 1+
s2
,
t = 0 : (1 + s)(N - 1).
(29)
In the s = 0 case, we have the following simple solution: v(t) = w(t)/N, t = 0 : (N - 1).
Subset II: p = 1 and s > N/2.
The (1, s)-projection is described as the (1, s - N )-projection, i.e., s is considered as s - N. Therefore, the
following convolution equation w = Av is used:

w(t) = N
1
+ sЇ
sЇ2
[v(t)
+
v(t
-
1)
+
v(t
-
2)
+
...
+
v(t
-

+
1)],
t = (N - 1) : -1 : -sЇ(N - 1),
(30)
where sЇ = N - s, and v(-k) = 0, when k > sЇ(N - 1)). We denote by M the negative number -sЇ(N - 1). The inverse transform b = A-1w can be calculated by the following recurrent form:
bM
= w(M ),
bM +1 bM +2
= w(M + 1) - bM , = w(M + 2) - (bM+1 + bM ),
bM +3 ...
= w(M + 3) - (bM+2 + bM+1 + bM ),
...
... ... ... ... ... ...
(31)
bM+sЇ-1 = w(M + sЇ - 1) - (bM+sЇ-2 + bM+sЇ-3 + ... + bM+1 + bM ),
bt
= w(t) - (bt-1 + bt-2 + ... + bt-sЇ+2 + bt-sЇ+1)),
t = M + sЇ, M + sЇ + 1, . . . , (N - 1).
The required sums v(t) of the discrete image are calculated by
v(t)
=
bt
sЇ N 1+
sЇ2
,
t = (N - 1) : -1 : -sЇ(N - 1).
(32)
The remaining two subsets of projections, namely, Subset III, when p N/2 and s = 1, and Subset IV, when p > N/2 and s = 1 are described similarly to Subset I and Subset II, respectively.
5.3 Preliminary results In this section, we illustrate a few results of reconstruction of the image f(x, y) on the lattice of size N Ч N, N = 2r, from 3N/2 projections. As an example, Figure 24 shows the image in part a, along with the graph of all line-integrals w(t) of the (1, 4)-projection in b, and ray-sums v(t) calculated from line-integrals in c. The
image
3000
1276 integrals w(t)
1276 sums v(t)
2500
600
2000 400 1500
1000 200 500
0
0
500
1000
0
0
500
1000
(a)
(b)
(c)
Figure 24. (a) Image, (b) line-integrals, and (c) ray-sums of the (1, 4)-projection.
splitting-signal {f1,4,t; t = 0 : 255} is shown in Figure 25 in part a. This signal is in the tensor representation and it defines the set of paired splitting-signals which are shown in b together. The first signal of length 128 is
tensor splitting-signal f 1,4,t 1700
8 paired splitting-signals 500
1600
1500
0
1400
1300 0
50
100 150 200 250
(a)
-500 0
50
100 150 200 250
(b)
Figure 25. Splitting-signals in (a) tensor and (b) paired representations.
{f1,4,t; t = 0 : 127} and the second signal is of length 64 is {f2,8,2t; t = 0 : 63}, and so on. The first four these splitting-signals are separated in the figure by the vertical dash lines. Figure 26 shows the modeled in MATLAB modified Shepp-Logan phantom image in part a, along with the graph of all line-integrals w(t) of the (1, 2)-projection in b, and ray-sums v(t) calculated from line-integrals in
phantom
15000
766 integrals w(t)
10000
5000
8000 6000 4000 2000
766 sums v(t)
0 0 200 400 600
0 0 200 400 600
(a)
(b)
(c)
Figure 26. (a) Phantom image, (b) line-integrals, and (c) ray-sums of the (1, 2)-projection.
c. The tensor splitting-signal {f1,2,t; t = 0 : 255} is shown in Figure 27 in part a. and the corresponding set of eight paired splitting-signals {f1,2,t; t = 0 : 127}, {f2,4,2t; t = 0 : 63}, {f4,8,4t; t = 0 : 63}, ..., {f64,64,0, f64,64,64},
12000 10000
tensor splitting-signal f 1,2,t
x 104 8 paired splitting-signals 1 0.5
8000
0
6000
-0.5
4000 0
50
100 150 200 250
(a)
-1 0
50
100 150 200 250
(b)
Figure 27. (a) Splitting-signal {f1,2,t} and (b) eight paired splitting-signals.
{f128,128,0}, and {f0,0,0} in b. The first four splitting-signals are separated by the vertical lines. Now we consider the results of the image reconstruction on examples with random rectangles in the square [0, 1] Ч [0, 1], we consider the objects with sharp edges. Figure 28 shows ten rectangles on the square with different intensities on the left. These rectangles are on the grid 128 Ч 128. The image calculated from 192 projections is shown on the right. The reconstruction is exact, and main program required 24 seconds to
10 rectangles and grid 128x128 1
reconstructed image 128x128
20
0.8
1
2
4
40
0.6
6
60
8 10
0.4
3
80
5
7
0.2
9
100
0
0
0.2 0.4 0.6 0.8
1
intensities: 3 5 4 1 2 5 2 5 4 4
120 20 40 60 80 100 120 cpu time/10=23.80sec
Figure 28. Ten random rectangles (left) and the reconstructed image (right).
calculate all projections, i.e., line-integrals w(t), and reconstruct the image. Incomplete 1-D DPT have been used in calculations. The CPU time is calculated as the average time for one rectangle. The process of calculation is organized in such a way, that for each (p, s)-projection, the integrals are calculated for each rectangle separately and, then, the results are added. The inverse 2-D 128 Ч 128-point DPT requires 0.5304 seconds. Figure 29 shows ten rectangles on the square [0, 1] Ч [0, 1] on the left. The coordinates of these rectangles are on the grid 256 Ч 256. The image calculated from 384 projections is shown on the right. The reconstruction is exact, and main program required 4.56 minutes to calculate all projections and reconstruct the image. Incomplete 1-D DPT have been used in calculations. The inverse 2-D 256 Ч 256-point DPT requires 4.4460 seconds. Time characteristics of calculation, when the proposed method was implemented in MATLAB and Ci+ are given in Table 4.5. The time includes the calculation of all line-integrals w(t), line-sums v(t), and the calculation of the 2-D direct and inverse paired transforms. Incomplete 1-D DPTs are used. The data in these tables were obtained from running programs on the personal computer with Intel Dual CPU Processor at 3.20GHz speed.
10 rectangles and grid 256x256 1
reconstructed image 256x256
0.8
50
6
1
4
10
0.6
100
5
2
79
0.4
150
8
0.2
3
200
0
0
0.2 0.4 0.6 0.8
1
intensities: 3 5 3 5 3 5 1 3 5 2
250 50 100 150 200 250 cpu time/10=4.63 min
Figure 29. Ten random rectangles (left) and the reconstructed image (right).
N ЧN 32 Ч 32 64 Ч 64 128 Ч 128 256 Ч 256 512 Ч 512
in MATLAB 00.48s/rec 02.58s/rec 23.80s/rec 04.63m/rec ......
in Ci+ 00.03s 00.11s 00.79s 06.17s 49.76s
Table 3. Time for scanning and reconstruction
This is the first realization of the method on MATLAB and Ci+, which can be improved, in order to achieve a fast reconstruction of images on the large size Cartesian lattices. To achieve the fast reconstruction of the image, we can describe each arithmetical ray, as well as each geometrical ray, as the unique set of IEs, through which the ray passes. This approach allows for calculating the line-integrals of the entire modeled image in one pass. As an example, Figure 30 shows the image of 20 random rectangles with the reconstruction. The total time of reconstruction is 4.98min, which is ten times smaller than
20 rectangles on grid 256x256 1
reconstructed image 256x256
0.8
50
0.6
100
0.4
150
0.2
200
0
0
0.2 0.4 0.6 0.8
1
250 50 100 150 200 250 cpu time=4.98 min
Figure 30. 20 random rectangles (left) and the reconstructed image (right).
the time given in Table 3. Figure 31 shows the result of reconstruction of the image composed by 13 ellipses and circles. The reconstruction is on the Cartesian lattice of size 256 Ч 256. The 256 Ч 256 reconstruction of the modified Shepp-Logan phantom is shown in Figure 32. The reconstruction is exact and all calculations that include the line-integrals, solutions of the convolution equations, and the inverse 2-D paired transform, require about five minutes, when running the code in MATLAB.
13 ellipces on grid 256x256 1
reconstructed image 256x256
0.8
50
0.6
100
0.4
150
0.2
200
0
0
0.2 0.4 0.6 0.8
1
250 50 100 150 200 250 cpu time=5.01 min
Figure 31. Image with 13 ellipses (left) and the reconstructed image (right).
phantom on the grid 256x256
reconstructed image 256x256
0.8
50
0.6
100
0.4
150
0.2
200
0
0.2 0.4 0.6 0.8
1
image f (x,y) d
250 50 100 150 200 250 cpu time 4.98 min
Figure 32. Image (left) and its reconstruction (right).
6. CONCLUSION The well-known model of the image reconstruction has been analyzed, when the image is represented by the set of small image elements with constant intensities each, and the rays are considered with zero width. In the framework of the considered model of the image f(x, y) on the unit square, the solution of the problem of image reconstruction from a finite number of projections is given. The solution is based on the properties of the 2-D tensor and paired transforms, each basis function of which is binary and defined by a set of parallel rays on the lattice. The main task in image reconstruction is the transformation of the geometry from image plane to the Cartesian lattice, on which the reconstructed image is calculated. We solve this problem, by introducing the set of G-rays for each projection and triangle Toeplitz matrix-based system of convolution equations. We describe the cases when image size is N Ч N, and N is a power of two and a prime. We believe that the proposed method can be used in tomographic imaging and the main idea of transferring geometry can be generalized for the parallel rays with not zero widths, as well as, for the fan beam projections. REFERENCES [1] A.C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging, IEEE Press, New York, 1987. [2] G.T. Herman, Image Reconstruction From Projections, New York: Academic Press, 1979. [3] D.E. Olins, A.L. Olins, H.A. Levy, R.C. Durfee, S.M. Margle, E.P. Tinnel, and S.D. Dover, "electron microscope tomography: transcription in three dimensions," Science, vol. 220, pp.498-500, 1983. [4] P. Suetens, Fundamentals of medical imaging, Cambridge University Press, 2002.
[5] Chen B, Ning R, Cone-beam volume CT breast imaging: feasibility study, Med Phys. vol. 29, pp. 755-770,
2002.
[6] Andrei IAGARU, Rinat MASAMED, Sravanthi KEESARA, and Peter S. CONTI, "Breast MRI and 18F
FDG PET/CT in the management of breast cancer," Annals of Nuclear Medicine, vol. 21, No. 1, pp. 3338,
2007.
[7] P. Pettigrew and J.L. Berry,
"The Role of CT in Breast Imaging,"
http://www.eradimaging.com/site/article.cfm?ID=739
[8] Glick S., Breast CT, Annu Rev Biomed Eng., vol. 9, pp. 501-526, 2007.
[9] Lai C, Shaw C, Geiser W, et al. "Comparison of slot scanning digital mammography system with full-field
digital mammography system," Med Phys., vol. 35, pp. 2339-2346, 2008.
[10] Nelson TR, Cervina LI, Boone JM. "Classification of breast computed tomography data," Med Phys.,
vol 35, pp. 1078-1086, 2008.
[11] Gong X, Glick SJ, Liu B, et al., "A comparison simulation study comparing lesion detection accuracy
with digital mammography, breast tomosynthesis, and cone-beam CT breast imaging," Med Phys., vol. 33,
pp. 1041-1052, 2006.
[12] Chen Z., Ning R., "Why should breast tumor detection go three dimensional? Phys Med Biol., vol. 48,
pp. 2217-2228, 2003.
[13] Chen B., Ning R., "Cone-beam volume CT breast imaging: feasibility study, Med Phys., vol. 29, pp. 755-770,
2002.
[14] "Breast CT reaches clinical testing: may improve on mammography, UC Davis health system Web
site. Available at: http://www.ucdmc.ucdavis.edu/newsroom/releases/archives/cancer/2005/breast-ct52005.html. Accessed September 26, 2009.
[15] Shimauchi A, Yamada T, Sato A, et al., "Comparison of MDCT and MRI for evaluating the intraductal
component of breast cancer, AJR Am J Roentgenol, 187, pp. 322-329, 2006.
[16] Inoue M, Sano T, Watai R, et al., " Dynamic multidetector CT of breast tumors: diagnostic features and
comparison with conventional techniques, AJR Am J Roentgenol, 181, pp. 679-686, 2003.
[17] Lindfors KK, Boone JM, Nelson TR, et al., "Dedicated breast CT: initial clinical experience, Radiology,
246, pp. 725-733, 2008.
[18] Radon J. On the determination of function from their integrals along certain manifolds, Ber. Saechs. Akad.
Wiss, Leipzig, Math. Phys. Kl., 69, pp. 262-277, 1917.
[19] S. Helgason, The Radon Transform, Progress in Mathematics 5, BirkhaЁuser Boston, 1980.
[20] R. Gordon, "A tutorial on ART (Algebraic reconstruction techniques)," IEEE Trans. Nuclear Science,
vol. 21, pp. 78-93, 1974.
[21] R. Gordon, R. Bender, and G.T. Herman, "Algebraic reconstruction techniques (ART) for three-dimensional
electron microscopy and X-ray photography," J. Theor. Biol., vol. 29, pp. 471-481, 1970.
[22] G.T. Herman and A. Lent, "Iterative reconstruction algorithms," Comput. Biol. Med., vol. 6, pp. 273-294,
1976.
[23] Y. Censor, "Finite series-expansion reconstruction methods," Proc. of IEEE, 1983, Vol 71, N0 3, pp. 409-419.
[24] R.H.T. Bates and M.J. McDonnell, Image restoration and Reconstruction, Charedon Press-Oxford, 1986.
[25] A.H. Andersen and A.C. Kak, "Simultaneous algebraic reconstruction technique (SART): a superior imple-
mentation of the ART algorithm," Ultrasonic Imaging, vol. 6, pp. 81-94, 1984.
[26] Mathematics and Physics of Emerging Biomedical Imaging, Washington, DC:NRC/IM, National Academic
Press, 1996.
[27] F. Natterer, Mathematical Methods in Image Reconstruction, Philadelphia, PA:SIAM, 2001.
[28] Y. Lifeng and S. Leng, "Image Reconstruction Techniques," 2010, http://www.imagewisely.org/Imaging-
Professionals/Medical-Physicists/Articles/Image- Reconstruction-Techniques.aspx
[29] A. M. Grigoryan, "An algorithm of the two-dimensional Fourier transform," Izvestiya VUZov SSSR,
Radioelectronica, vol. 27, no. 10, pp. 52-57, USSR, Kiev 1984. (translated in http://fasttransforms.com/Art-
USSR-papers/Art- IzvuzovSSSR1984.pdf ) [30] A. M. Grigoryan and M. M. Grigoryan, "Two-dimensional Fourier transform in the tensor presentation and
new orthogonal functions," Avtometria, AS USSR Siberian Section, no. 1, pp. 2127, 1986, Novosibirsk.
[31] A.M. Grigoryan, "New algorithms for calculating discrete Fourier transforms," Journal Vichislitelnoi Matematiki i Matematicheskoi Fiziki, vol. 26, no 9, pp. 1407-1412, AS USSR, Moscow 1986. (translated in http://fasttransforms.com/Art-USSR-papers/Art- JVMATofASUSSR1986.pdf ) [32] A.M. Grigoryan, "An algorithm for computing the two-dimensional Fourier transform of equal orders," Proceedings of the IX All-Union Conference on Coding and Transmission of Information, Part 2, pp. 49-52, USSR, Odessa 1988. (translated in http://fasttransforms.com/Art-USSR-papers/Art-AllUnConf1988.pdf) [33] A.M. Grigoryan, "2-D and 1-D multi-paired transforms: Frequency-time type wavelets," IEEE Trans. on signal processing, vol. 49, no. 2, pp. 344-353, Feb. 2001. [34] A.M. Grigoryan, and S.S. Agaian, "Split manageable efficient algorithm for Fourier and Hadamard transforms," IEEE Trans. on Signal Processing, vol. 48, no. 1, pp. 172-183, Jan. 2000. [35] A.M. Grigoryan and S.S. Agaian, "Shifted Fourier transform based tensor algorithms for 2-D DCT," IEEE Transaction on Signal Processing, vol. 49, no. 9, pp. 2113-2126, Sep. 2001. [36] S.S. Agaian, Hadamard matrices and their applications, lecture notes in Mathematics 1168, Springer ­ Verlag, Berlin, Heitelbery, New York, Tokyo 1985. [37] S.S. Agaian, "Advances and problems of the fast orthogonal transforms for signal-images processing applications (part 1)," Pattern Recognition, Classification, Forecasting. Yearbook, The Russian Academy of Sciences, Nauka, Moscow, Issue 3, pp. 146-215, 1990, (in Russian). [38] S.S. Agaian, "Advances and problems of the fast orthogonal transforms for signal-images processing applications (part 2)," Pattern Recognition, Classification, Forecasting. Yearbook, The Russian Academy of Sciences, Nauka, Moscow, Issue 4, pp. 156-246, 1991, (in Russian). [39] R. Kogan, S.S. Agaian, and K.A. Panetta, "Visualization using rational morphology and zonal magnitude reduction", Proceedings, IS&T/SPIEs Symposium on Electronic Imaging Science & Technology, San Jose. CA February Vol. 3304, IX, pp. 153-163, January 26, 1998. [40] S.S. Agaian, K. Panetta, and A.M. Grigoryan, "Transform-based image enhancement algorithms," IEEE Trans. on Image Processing, vol. 10, no. 3, pp. 367-382, March 2001. [41] S.S. Agaian, K.A. Panetta, and A.M. Grigoryan, "A new measure of image enhancement", Proceedings, IASTED International Conference on Signal Processing and Communications, July 3-6, 2001. [42] A.M. Grigoryan and S.S. Agaian, "Transform-based image enhancement algorithms with performance measure," Advances in Imaging and Electron Physics, Academic Press, vol. 130, pp. 165-242, May 2004. [43] F.T. Arslan and A.M. Grigoryan, "Fast splitting -rooting method of image enhancement: Tensor representation," IEEE Trans. on Image Processing, vol. 15, no. 11, pp. 3375-3384, Nov. 2006. [44] A.M. Grigoryan, "Method of paired transforms for reconstruction of images from projections: Discrete model," IEEE Trans. on Image Processing, vol. 12, no. 9, pp. 985-994, Sep. 2003. [45] A.M. Grigoryan and M.M. Grigoryan, Image Processing: Tensor Transform and Discrete Tomography with MATLAB, CRC Press Taylor and Francis Group, September 25, 2012. [46] A.M. Grigoryan and Nan Du, "Principle of superposition by direction images," IEEE Trans. on Image Processing, vol. 20, no. 9, pp. 2531-2541, September 2011. [47] A.M. Grigoryan and Nan Du, "2-D images in frequency-time representation: Direction images and resolution map," Journal of Electronic Imaging, vol. 19, no. 3, (033012) pp. 1-14, July-September 2010. [48] A.M. Grigoryan, Multidimensional Discrete Unitary Transforms, book-chapter 19 (69 pages) in the Third Edition Transforms and Applications Handbook in The electrical engineering Handbook Series (Editor in Chief, Prof. Alexander Poularikas) CRC Press Taylor and Francis Group, 2010. [49] A. M. Grigoryan and S. S. Agaian, Multidimensional Discrete Unitary Transforms: Representation, Partitioning, and Algorithms, New York: Marcel Dekker, 2003. [50] A. M. Grigoryan and M. M. Grigoryan, Brief Notes in Advanced DSP: Fourier Analysis With MATLAB., CRC Press Taylor and Francis Group, 2009.

File: another-step-towards-successful-tomographic-imaging-in-cancer.pdf
Title: SPIEchapter_Art.dvi
Published: Tue Jan 1 00:00:00 2
Pages: 31
File size: 0.54 Mb


TH Etech, 6 pages, 0.67 Mb
Copyright © 2018 doc.uments.com