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 Xray equipments are used in tomography to obtain crosssectional 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 rayintegrals of the image f(x, y) are transformed uniquely into the raysums 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 splittingsignals carrying the spectral information of the image at frequencypoints 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 thirdgeneration scanners wherein the detector array cover the entire field of view and fanbeam geometry is used instead of parallelbeam geometry.4,5 Computerized tomography which also is called the computerized axial tomography (CAT) is a diagnostic procedure that uses a special Xray equipment to obtain crosssectional 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 followup. 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, 18FFDG PET/CT more accurately showed lesions in the chest, abdomen and bones. The study in6 suggests that 18FFDG PET/CT should be considered as complimentary imaging tools in the pre and postoperative workup of patients diagnosed with breast cancer. Further author information: Artyom M. Grigoryan: Email:
[email protected], Telephone: (210)4587518.
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 3dimensional (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 conebeam technology for breast CT has been shown to benefit image quality. A
simulation study conducted by Chen and Ning showed the feasibility of conebeam 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 evidencebased 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 lineintegrals (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 backprojection, iterative reconstruction, Fourier filtering, Radon filtering, and convolution filtering20.28 The reconstruction is performed with calculation of the tensor or paired transform2933 of the image from the projection data. In the tensor representation, a twodimensional (2D) discrete image is considered as the sum of direction images, each of which is determined by the corresponding signal. These signals are called the splittingsignals; they also carry the spectral information of the image at frequencypoints of different subsets that cover the whole domain of frequencies. The discrete image can be reconstructed by its splittingsignals by calculating the 2D DFT or directly by performing the inverse 2D 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 raysums of the discrete image in the Cartesian lattice. However, it is very important to note that these sumrays can exactly be calculated from the rayintegrals. In other words, the splittingsignals 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 2D 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 2D 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 wellknown series expansion methods for solving the system of line integrals in the discrete case20.25
2.1 Lineintegrals and raysums
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 raysums,
vl =
fn,m.
(1)
(n,m)l
In the above model, the lineintegral 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 lineintegral 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 lineintegrals 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 raysums in (1) depending on the projection, but in all cases the raysums can be uniquely calculated from lineintegrals.
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 coordinate systems on the plane are considered. The first system of coordinates (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 coordinate 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 frequencypoints. Given the frequencypoint (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 coordinate 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,  tan1(p/s), of these rays. The set of lineintegrals 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 SPLITTINGSIGNALS In this section we describe the images and their twodimensional discrete Fourier transforms (2D DFT) by the unique sets of 1D signals and 1D DFTs, respectively. The images are considered in the tensor and paired representations which are the 2D frequency and 1D 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 2D DFT and other transforms, including 2D Hartley, Hadamard, and cosine transforms34.36,48 The tensor transform can be used in
image enhancement40,42 image denoising,43,49,50 and discrete image reconstruction.44,45 In the tensor representation, the discrete image fn,m is the set of 1D splittingsignals,
: {fn,m} fTp,s = {fp,s,t; t = 0 : (N  1)} (p,s)JN,N .
(4)
Here JN,N is a set of frequencypoints (p, s), or generators of the splittingsignals which are selected in a way to cover the lattice XN,N of frequencypoints (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 splittingsignals 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) Splittingsignal {f4,1,t t = 0 : 255} of the image and (b) its direction image.
and the set of all splittingsignals 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 splittingsignals with the remaining generators {(1, 2s); s = 0 : (N/2  1)} of the set J256,256. The splittingsignals are written along the rows in these two matrices. Figure 4 shows the splittingsignal of the image, which is generated by the frequencypoint (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 splittingsignal fTp,s carries the spectral information of the 2D DFT at N frequencypoints 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 splittingsignal, or imagesignal fT2,1 of length 256 is shown in part a, along with the magnitude of the 1D 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 2DFT of the image is the 1DFT of the imagesignal. The 2D DFT at frequencypoints 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 2DFT at frequencypoints 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 1D DFT is filled the 2D 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 splittingsignals 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
imagesignal 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°
1D DFT of f
T
40
2,1
30
amplitude
20
10
0
0
50
100 150 200 250
(c)
(d)
Figure 5. (a) The imagesignal corresponding to the set T2,1. (b) The image after amplifying the 2DFT at frequencypoints of T2,1. (c) Magnitude of the 1D DFT of the imagesignal (zero component shifted to the center and truncated). (d) The 2D DFT of the image with amplified samples of the set T2,1 .
Here JN,N is a set of 3N  2 frequencypoints (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 splittingsignals 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 splittingsignals carry information of the 2D DFT of the image at frequencypoints of the subset Tp,s = {(2m + 1)(p, s); m = 0 : N/2k  1} Tp,s. The 2D 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 1D DFTs of
splittingsignals fTp,s define completely the 2D DFT of the image. The tensor representation is unique, and
the image can be defined through the 2D DFT calculated by (7) or directly from the tensor transform. Indeed,
according to the principle of superposition33,4648 the image fn,m can be composed from splittingsignals as
follows:
fn,m
=
1 r1 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
2k1 p,2k1s,2k1t1
2k1p,2k1 s,2k1t1+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 splittingsignals 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 splittingsignal 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 splittingsignals 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 splittingsignals 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 raysums of the discrete image from the lineintegrals 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·
9·
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 lineintegrals of the (1, 2)projection.
see, that the lineintegral 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 lineintegral 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 lineintegrals 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 linesums 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 lineintegral 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 lineintegrals 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
A1
w,
or
by
the
method
of direct substitution in the following recursive form:
b0 = w(0),
b1 = w(1)  b0,
(14)
bt = w(t)  bt1, 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 splittingsignal {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 tripletnumbers 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 + v7, f1,5,1 = v1 + v6, f1,5,2 = v2 + v5 + v12, f1,5,3 = v3 + v4 + v11, f1,5,4 = v4 + v3 + v10, f1,5,5 = v5 + v2 + v9, f1,5,6 = v6 + v1 + v8,
(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 =
vt7m =
v1,2,t7m,
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 Grays is shown in Figure 14. For this set of Grays, the system of equations for the lineintegrals
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
b12 = w(12), b11 = w(11)  b12, b10 = w(10)  b11, bt = w(t)  bt1, 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 =
v6t+7m,
m=0
t = 0 : 6,
(vt = 0, t > 18).
4.1 Main Equations for Grays 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 lineintegrals 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 splittingsignal {f1,s,t; t = 0 : (N  1)} is calculated from the lineintegrals 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 lineintegrals used in the reconstruction. The number of parallel Grays 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+1N 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(s1))
where
s
N 1 2
b(0) = w(0), bt = w(t)  (bt + · · · + bt(Ns1))
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 vt7m,
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 splittingsignals 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 transformbased 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 frequencypoint (p, s), the (p, s)projection is calculated, i.e., all lineintegrals w(t) along the set of geometrical rays ~lp,s(t). 4. The geometry of Grays of the (p, s)projection is transformed to Arays, 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+s1)(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 Grays to Arays is fast. 5. The splittingsignal {fp,s,0, fp,s,1, ..., fp,s,N1} is calculated from v(t). This signal is written into the corresponding row of the matrix N Ч (N + 1) of the 2D discrete tensor transform. 6. The inverse 2D tensor transform is calculated, which is the reconstructed discrete image fn,m.
image fd(x, y) Image composition on the unit square 1 ? Lineintegrals of the (p, s)projection 3 ? Transform of Grays into Arays, w v 4
2 Set of generators JN,N = {(p, s)}
? Calculation of the splittingsignal {fp,s,t} 5 ? Inverse 2D discrete tensor transform 6
? reconstructed image fn,m
Figure 15. Blockdiagram 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 lineintegrals w(t) on the left, the sums v(t) in the middle, and the splittingsignal {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 lineintegrals w(t), linesums 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 lineintegrals in a way, similar to the described above case when N is a prime. The effective reconstruction is with the 2D
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. Lineintegrals w(t), raysums v(t), and splittingsignal 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 Cibased 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 splittingsignals from lineintegrals. 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 2D 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 tripletnumbers of the paired functions is taken as
r1
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 tripletnumbers (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  tan1(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 tripletnumbers (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 tripletnumbers (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 raysums 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 2D 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 tripletnumbers (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 tripletnumbers:
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 2D 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 8point 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 lineintegrals 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 lineintegrals 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
A1
=
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
=
A1w,
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)  (bt1 + bt2 + bt3), t = 5, 6, . . . , 35.
Similar calculations hold for other projections, and the raysums can be calculated from the lineintegrals 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 lineintegrals are used to calculate the splittingsignals in the paired or tensor representation of the discrete image. Therefore, the blockdiagram 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 2D 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 tripletnumber (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 tripletnumber) equals 2k. Therefore, there are two ways to use this property in image reconstruction. 1. One can use the complete 1D 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 2D paired transform, which have been calculated by 1D 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 1D paired transforms over vectors a, to avoid the repetition of calculation of components of the 2D paired transform of the image. For instance, when N = 8, the complete 1D 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 1D 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 Grays is defined in the following way. We first consider the Arays 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 lineintegrals, 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 lineintegrals 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 lineintegrals 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 = A1w 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)
... ...
... ... ... ...
bs1 = w(s  1)  (bs2 + bs3 + ... + b1 + b0),
bt = w(t)  (bt1 + bt2 + ... + bts+2 + bts+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

sЇ
+
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 = A1w 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)  (bt1 + bt2 + ... + btsЇ+2 + btsЇ+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 lineintegrals w(t) of the (1, 4)projection in b, and raysums v(t) calculated from lineintegrals 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) lineintegrals, and (c) raysums of the (1, 4)projection.
splittingsignal {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 splittingsignals which are shown in b together. The first signal of length 128 is
tensor splittingsignal f 1,4,t 1700
8 paired splittingsignals 500
1600
1500
0
1400
1300 0
50
100 150 200 250
(a)
500 0
50
100 150 200 250
(b)
Figure 25. Splittingsignals 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 splittingsignals are separated in the figure by the vertical dash lines. Figure 26 shows the modeled in MATLAB modified SheppLogan phantom image in part a, along with the graph of all lineintegrals w(t) of the (1, 2)projection in b, and raysums v(t) calculated from lineintegrals 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) lineintegrals, and (c) raysums of the (1, 2)projection.
c. The tensor splittingsignal {f1,2,t; t = 0 : 255} is shown in Figure 27 in part a. and the corresponding set of eight paired splittingsignals {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 splittingsignal f 1,2,t
x 104 8 paired splittingsignals 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) Splittingsignal {f1,2,t} and (b) eight paired splittingsignals.
{f128,128,0}, and {f0,0,0} in b. The first four splittingsignals 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., lineintegrals w(t), and reconstruct the image. Incomplete 1D 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 2D 128 Ч 128point 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 1D DPT have been used in calculations. The inverse 2D 256 Ч 256point 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 lineintegrals w(t), linesums v(t), and the calculation of the 2D direct and inverse paired transforms. Incomplete 1D 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 lineintegrals 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 SheppLogan phantom is shown in Figure 32. The reconstruction is exact and all calculations that include the lineintegrals, solutions of the convolution equations, and the inverse 2D 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 wellknown 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 2D 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 Grays for each projection and triangle Toeplitz matrixbased 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.498500, 1983. [4] P. Suetens, Fundamentals of
medical imaging, Cambridge University Press, 2002.
[5] Chen B, Ning R, Conebeam volume CT breast imaging: feasibility study, Med Phys. vol. 29, pp. 755770,
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. 501526, 2007.
[9] Lai C, Shaw C, Geiser W, et al. "Comparison of slot scanning digital mammography system with fullfield
digital mammography system," Med Phys., vol. 35, pp. 23392346, 2008.
[10] Nelson TR, Cervina LI, Boone JM. "Classification of breast computed tomography data," Med Phys.,
vol 35, pp. 10781086, 2008.
[11] Gong X, Glick SJ, Liu B, et al., "A comparison simulation study comparing lesion detection accuracy
with digital mammography, breast tomosynthesis, and conebeam CT breast imaging," Med Phys., vol. 33,
pp. 10411052, 2006.
[12] Chen Z., Ning R., "Why should breast tumor detection go three dimensional? Phys Med Biol., vol. 48,
pp. 22172228, 2003.
[13] Chen B., Ning R., "Conebeam volume CT breast imaging: feasibility study, Med Phys., vol. 29, pp. 755770,
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/breastct52005.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. 322329, 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. 679686, 2003.
[17] Lindfors KK, Boone JM, Nelson TR, et al., "Dedicated breast CT: initial clinical experience, Radiology,
246, pp. 725733, 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. 262277, 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. 7893, 1974.
[21] R. Gordon, R. Bender, and G.T. Herman, "Algebraic reconstruction techniques (ART) for threedimensional
electron microscopy and Xray photography," J. Theor. Biol., vol. 29, pp. 471481, 1970.
[22] G.T. Herman and A. Lent, "Iterative reconstruction algorithms," Comput. Biol. Med., vol. 6, pp. 273294,
1976.
[23] Y. Censor, "Finite seriesexpansion reconstruction methods," Proc. of IEEE, 1983, Vol 71, N0 3, pp. 409419.
[24] R.H.T. Bates and M.J. McDonnell, Image restoration and Reconstruction, Charedon PressOxford, 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. 8194, 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/MedicalPhysicists/Articles/Image ReconstructionTechniques.aspx
[29] A. M. Grigoryan, "An algorithm of the twodimensional Fourier transform," Izvestiya VUZov SSSR,
Radioelectronica, vol. 27, no. 10, pp. 5257, USSR, Kiev 1984. (translated in http://fasttransforms.com/Art
USSRpapers/Art IzvuzovSSSR1984.pdf ) [30] A. M. Grigoryan and M. M. Grigoryan, "Twodimensional 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. 14071412, AS USSR, Moscow 1986. (translated in http://fasttransforms.com/ArtUSSRpapers/Art JVMATofASUSSR1986.pdf ) [32] A.M. Grigoryan, "An algorithm for computing the twodimensional Fourier transform of equal orders," Proceedings of the IX AllUnion Conference on Coding and Transmission of Information, Part 2, pp. 4952, USSR, Odessa 1988. (translated in http://fasttransforms.com/ArtUSSRpapers/ArtAllUnConf1988.pdf) [33] A.M. Grigoryan, "2D and 1D multipaired transforms: Frequencytime type wavelets," IEEE Trans. on
signal processing, vol. 49, no. 2, pp. 344353, 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. 172183, Jan. 2000. [35] A.M. Grigoryan and S.S. Agaian, "Shifted Fourier transform based tensor algorithms for 2D DCT," IEEE Transaction on Signal Processing, vol. 49, no. 9, pp. 21132126, 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 signalimages processing applications (part 1)," Pattern Recognition, Classification, Forecasting. Yearbook, The
Russian Academy of Sciences, Nauka, Moscow, Issue 3, pp. 146215, 1990, (in Russian). [38] S.S. Agaian, "Advances and problems of the fast orthogonal transforms for signalimages processing applications (part 2)," Pattern Recognition, Classification, Forecasting. Yearbook, The Russian Academy of Sciences, Nauka, Moscow, Issue 4, pp. 156246, 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. 153163, January 26, 1998. [40] S.S. Agaian, K. Panetta, and A.M. Grigoryan, "Transformbased image enhancement algorithms," IEEE Trans. on Image Processing, vol. 10, no. 3, pp. 367382, 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 36, 2001. [42] A.M. Grigoryan and S.S. Agaian, "Transformbased image enhancement algorithms with
performance measure," Advances in Imaging and Electron Physics, Academic Press, vol. 130, pp. 165242, 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. 33753384, 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. 985994, 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. 25312541, September 2011. [47] A.M. Grigoryan and Nan Du, "2D images in frequencytime representation: Direction images and resolution map," Journal of Electronic Imaging, vol. 19, no. 3, (033012) pp. 114, JulySeptember 2010. [48] A.M. Grigoryan, Multidimensional Discrete Unitary Transforms, bookchapter 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.