(4)
Using these notations, equation (3) with Tp = 1 is replaced by the following discrete relationship
(4)
which, by introducing the matrix [A], defined by
(5)
and called in the following the imaging matrix, can be written in the synthetic form
(5)
where wj is a random vector generated
by a white Gaussian process.
The imaging matrix [A] is rectangular, with N
rows and N + 2K columns. Since this matrix is ill-conditioned,
this solution can be corrupted by an amplified propagation of the data
noise, so that regularization methods must be used for controlling this
noise propagation.
Taking into account that restored images must be positive and not corrupted by noise amplification, we have implemented a particular version of so-called projected Landweber method (Eicke 1992), proposing the following iterative method:
(6)
where:
denotes the
transposed of the matrix [A];
(7)
is a positive parameter,
known as relaxation parameter, which, for our problem, must be smaller
than 0.125 (in our code 0.1125).
between the j-th column of the measured data and the j-th column of the
data computed by means of the k-th iterate as the root means square (r.m.s.)
value of the vector 
:
(8)
Even if the value of 
is used for all columns, the number of iterations in general is changing
from column to column: the number of iterations is small if the column
is characterized by a low value of the signal-to-noise ratio S/N and is
larger if S/N is higher. If one does not expect the ratio S/N to change
dramatically from column to column it may be more convenient to use a second
stopping rule which provides the same number of iterationssame number of
iterations for all columns. To this purpose we define the average relative
discrepancy as follows
(9)
is
smaller than some estimated value 
of the relative r.m.s. error affecting the image.