Kaczmarz extended algorithm for tomographic image reconstruction from limited-data

https://doi.org/10.1016/j.matcom.2004.01.021Get rights and content

Abstract

The algebraic reconstruction technique (ART), based on the well known algorithm proposed by Kaczmarz in 1937, is one of the most important class of solution methods for image reconstruction problems. But unfortunately, almost all the methods from the ART class give satisfactory results only in the case of consistent problems. In the inconsistent case (and unfortunately this is what happens in real applications, because of measurement errors) they give only more or less “controllable” versions of the exact solutions. This is exactly the case that we analyze in the present paper. We start with a theoretical analysis of the classical Kaczmarz’s projection method in the case of an inconsistent linear least-squares problem and we prove that the approximations so obtained are at a certain distance from the set of exact least-squares solutions. This distance is controlled by the component of the right hand side of the problem lying in the orthogonal complement of the range of problem’s matrix, i.e. exactly the component that makes the problem inconsistent. For overcoming this difficulty we consider an extended version of Kaczmarz’s algorithm, previously analyzed by one of the authors. In the numerical experiments described in the last part of the paper we compare the above mentioned extension with two well known (ART) type algorithms for image reconstruction in two electromagnetic geotomography problems. The results indicate that the extended Kaczmarz algorithm gives much better results than the other two.

Section snippets

Kaczmarz’s algorithm for inconsistent problems

We consider the linear least-squares problem: find xRn such that ||Ax−b||=min!,where A is an m×n real matrix and b∈Rm a given vector (by ||·|| and 〈·,·〉 we shall denote the Euclidean norm and scalar product on some space Rq). Let LSS(A,b) be the set of all its least-squares solutions and xLS its minimal norm one. Also Bt,R(B),N(B) will denote the transpose, range and null space of a matrix B, respectively, and PS the orthogonal projection onto a closed convex set S⊂Rq. By aiRn, we shall

Kaczmarz’s algorithm with relaxation parameters

The results from , in the previous section tell us that, when applying classical Kaczmarz algorithm , to an inconsistent problem (1) for which the “perturbation from consistency” bA is big (see also Remark 5), the obtained approximation x(x0) can be enough far from the solution set. This can explain the troubles appearing in image reconstruction problems for which (because of measurement errors) the right hand side b changes and (1) becomes inconsistent (see also the comments in [1], [3],

The extended Kaczmarz algorithm

In [14] (see also [13], [15]) one of the authors proposed a different extension of Kaczmarz’s algorithm , , for dealing with the inconsistent case for (1). For this, we have to consider, beside the orthogonal projections fi(b;x) from (7) onto the hyperplanes generated by the rows ai of A, similar ones but with respect to the columns, denoted here by αj≠0,j=1,…,n. They are defined by ϕj(y)=y−〈y,αj||αj||2αj,j=1,…,n,Φ(y)=(ϕ1,…,ϕn)(y).Then the Extended Kaczmarz algorithm is defined as follows.

Algorithm KE

Let

Right hand side perturbation

As we have already mentioned at the beginning of Section 2, in practical applications perturbations can appear in the right hand side of (1) because of measurement errors in the process of scanning the given object. Let us suppose that the (initial) problem (1) is consistent, i.e. (11) holds. The perturbed one will be ||Ax*b̃||=min!,where b̃=b+δb,δb=δbA+δbA,δbA∈R(A),δbA∈N(At).According to the considerations and results from Section 1, if x0 is the initial approximation, the limit of the

Conclusions

The results in Section 4.4 prove the theoretical considerations discussed in Section 4.1. Fig. 15, Fig. 17, Fig. 19 show that for KERP applied to noisy data with δbA=0, z*(x0=0) in (71) converges to xLS, and hence the distance in (72) approaches to zero. Increasing a number of iterations in KERP for δbA=0, we observe a gradual improvement in the quality of the reconstructed images (see Fig. 3, Fig. 9), whereas this does not hold for KR and CEG (see Fig. 5, Fig. 7, Fig. 11, Fig. 13). This is in

Acknowledgements

The paper was supported by the NATO Grant PST.CLG.977924 and the DAAD grant that one of the authors had as a visiting professor at the Friedrich-Alexander-University Erlangen-Nuremberg, Germany, in the period October 2002 to August 2003.

Rafal Zdunek was born on 26 June 1972, Brzesko, Poland. Presently he is a Polish. He is a Lecturer in the Institute of Telecommunications and Acoustics (ITA), Wroclaw University of Technology (WUT), Wybrzeze Wyspianskiego 27, 50-370 Wroclaw, Poland; Tel.: +48-71-3203215; fax: +48-71-3223664; E-mail address: [email protected]. He received his MSc degree in 1997 in broadcast telecommunication and PhD degree in 2002 in telecommunication (specialization in numerical methods). From 2001–2002,

References (17)

There are more references available in the full text version of this article.

Cited by (74)

View all citing articles on Scopus

Rafal Zdunek was born on 26 June 1972, Brzesko, Poland. Presently he is a Polish. He is a Lecturer in the Institute of Telecommunications and Acoustics (ITA), Wroclaw University of Technology (WUT), Wybrzeze Wyspianskiego 27, 50-370 Wroclaw, Poland; Tel.: +48-71-3203215; fax: +48-71-3223664; E-mail address: [email protected]. He received his MSc degree in 1997 in broadcast telecommunication and PhD degree in 2002 in telecommunication (specialization in numerical methods). From 2001–2002, he worked as a research assistant at the ITA, (WUT). From 2002 to till date, he is a lecturer at the ITA, (WUT). He has published 18 scientific papers in proceedings of national and international conferences. His current research interests include electromagnetic geotomography (limited-data tomography), algorithms for tomographic image reconstruction, inverse problems–regularization techniques. Since 2003, he has been a member of SEP (Society of Polish Electrical Engineers), Poland. He can speak two foreign languages: English and French.

Constantin Popa was born on 10 October 1956, Bucharest, Romania. Presently he is a Romanian. He is a professor of applied mathematics and head of the Department of Computer Science and Numerical Methods, Faculty of Mathematics and Computer Science, “OVIDIUS” University, Blvd. Mamaia 124, 900527 Constanta, Romania; Tel.: +40-241-618070; fax: +40-241-618372; E-mail address: [email protected]. From September 1996 to September 1997, he was a visiting scientist (postdoc grant) at Department of Applied Mathematics and Computer Science, Weizmann Institute of Science, Rehovot, Israel. From 1 October 2002 to 30 August 2003 he was a visiting professor (DAAD Gastdozent) at Lehrstuhl für Informatik 10 (Systemsimulation), “Friedrich-Alexander” University Erlangen-Nürnberg, Germany. He has published 3 books and more than 40 papers in national and international refereed journals. His current research interests include preconditioning techniques for finite element and finite differences discretizations of boundary value problems, iterative methods for least-squares formulations of linear systems of equalities and inequalities (projection algorithms), inverse problems—regularization techniques and methods for approximating the minimal norm solution of first kind integral equations. He was a member of GAMM Society, Germany (from 1992). From February 2001 to July 2002, he was a secretary for applied mathematics of the Romanian Section of GAMM. From December 2001, he was a reviewer at Zentralblatt fuer Mathematik, Germany. He can speak three foreign languages: English, French, Italian.

1

On leave from“Ovidius” University of Constanţa, Romania.

View full text