Lima J.J.Pedroso, de (ed.). Nuclear Medicine Physics

Подождите немного. Документ загружается.

248 Nuclear Medicine Physics

Ramp

1.0

0.8

0.6

0.4

0.2

0.2 0.4 0.6 0.8 1.0

0.0

Hamming

Hann

ν/ν

w(ν)/ν

FIGURE 6.20

The ramp ﬁlter, truncated in the (spatial) Nyquist frequency vN = 1/2Δxr, and the two most

common apodization windows used, the Hann and the Hamming windows.

methods, which, by assuming that the value of a projection is proportional

to the integral of the activity along the corresponding projection line, force

any data correction procedure to be performed before the reconstruction. The

main disadvantage of the iterative methods compared with the analytical

ones is their slowness, because they are far more intensive from the point of

view of numerical processing (typically they take up to 2 orders of magnitude

times longer than analytical methods to reconstruct the same image). How-

ever, advances in the processing power of modern computers observed in the

lastfew years have madethese methods commonplace inmany circumstances

of clinical practice, and they are nowadays the predominant reconstruction

methods used in commercial cameras. Still, the nonlinear character of itera-

tive methods may, on rare occasions, give rise to erratic behavior capable of

seriously affecting the reconstructed image.

The working principle of iterative methods is simple: Plane by plane,

one generates an estimate of the activity distribution in the patient,

est

(x, y, z)



z=z

, which is iteratively improved by comparing the parallel

projections p

est

, φ, z)



z=z

generated by the estimate with the projections

p(x

, φ, z)



z=z

actually measured. If the two projections are signiﬁcantly dif-

ferent, the estimate of the activity distribution is corrected by an amount that

is a function of the differences found between the projections (Figure 6.21)

[84]. This cyclic procedure is normally initiated with an empty (or sometimes

unitary) activity distribution, and it is repeated until the projections gener-

ated by the estimate and the measured projections are equal (within a preset

tolerance), or a certain condition on the activity distribution is met.

Imaging Methodologies 249

Initial estimate

Current estimate of

the activity distribution

Projections generated

by the estimate

Incluses image

formation models

(attenuation, scattered

radiation, resolition, ...)

est

(x,y,z)

Δf

est

(x,y,z)

Compare

Corrections to the estimate

of the activity distribution

Measure projections

est

(x,ϕ,z)

p(x,ϕ,z)

FIGURE 6.21

Outline of an iterative reconstruction algorithm. The cycle shown is processed plane by plane,

until a stopping condition is met.

Obtaining the projections p

est

, φ, z)



z=z

from the estimate of the activ-

ity distribution f

est

(x, y, z)



z=z

is a fundamental component of the iterative

process. This task is carried out using a system matrix A, which establishes

a relation between the values measured along the different projection lines,

often organized in the shape of a one-dimensional data vector, y, with the

volume elements that constitute the object, usually discretized in voxels

whose values are stored in a vector x:

y = Ax +n. (6.30)

The one-dimensional vector n contains the noise present in the measure-

ment, which can be modeled in many different ways although it is customary

to assume that it obeys a Poisson or Gaussian probability function.

The existing iterative methods mainly differ in how the projections gen-

erated by the activity estimate and the measured projections are compared

and also in how the corrections to the estimate in each iteration are computed

from the differences between those two (this also implies different process-

ing times for different methods). One of the simplest iterative methods is the

Algebraic Reconstruction Technique (ART), which is based on a purely alge-

braic iterative algorithm that improves the activity estimate angle by angle,

comparing each measured projection with the calculated projection using a

simple subtraction [85]. This method, although simple, is slow compared with

250 Nuclear Medicine Physics

most other iterative algorithms for the number of pixels commonly used in

SPECT and PET, and it is, therefore, rarely used.

6.3.1.2.2.1 The MLEM Algorithm The Maximum Likelihood Expectation

Maximization (MLEM) algorithm describes the acquisition of counts in the

parallel projections, y,as



, (6.31)

where x

is the activity in the voxel i, contributing with weight a

for the

projection channel j. The inversion of the matrix A ={a

} in Equation 6.31,

which would lead to the direct solution of x to yield the activity distribu-

tion, is impractical for several reasons: (1) The elements of the system matrix

are not known with enough accuracy, as they depend on the activity x

due to Compton scattering; (2) There might not be enough measured pro-

jections to render Equation 6.31 a fully determined system of equations; and

(3) The statistical nature of radioactive decay and its detection introduces

ambiguities in the system of Equation 6.31. For these reasons, the inversion of

Equation 6.31 has to be done under further assumptions, which leads to a

solution given by an iterative procedure.

The quantitative aspect of that procedure may be deduced by consider-

ing that the probability of observing a given distribution of projected counts,

y = y

obs

, given a hypothetical voxel activity distribution, x

hyp

, can be maxi-

mized if maximum likelihood arguments are applied to those variables. The

probability of success, that is, of detecting the activity of the voxel x

in a given

projection channel y

, is much smaller than the total number of possibilities,

leading to the conclusion that the detection process obeys Poisson statistics.

Further, since a pair of photons interacting in a given pair of crystals cannot

be detected in any other pair (ignoring the existence of Compton scattering

in the crystals), one concludes that, besides obeying Poisson statistics, the

variables are independent as well. It is, therefore, possible to use the Poisson

probability distribution to write the conditional probability, P(y

obs

hyp

), and

maximize it in the sense of having maximum likelihood. The manipulation of

that function [86] shows that it is concave and that the iterative convergence

to the maximum follows the recurrence relation

= x

n−1





n−1

, (6.32)

where n is the iteration number, i = 1, ... and k = 1, ... refer to the voxels

in the image, and j = 1, ... refers to the projection channels. The sum in k in

the denominator of the second factor expresses the projection of the activity

distribution x

n−1

, that is, the measured projection counts if x

n−1

is the true

image. The sum in j in the numerator corresponds to the multiplication of

the measured projections by the transpose of the system matrix A, and it

Imaging Methodologies 251

represents the backprojection of the ratio between the measured data and the

estimated projections in the iteration n. The vector x

is usually a uniform

distribution (x

= 1).

This algorithm can be easily implemented in series, processing the data

directly from list-mode format, as it works only with the projection chan-

nels in the vector y, which contribute to the backprojection, that is, have

nonzerocounts. In terms of its stability, the algorithm converges to the activity

estimate that best approaches the data in terms of a Poisson likelihood. How-

ever, this approximation assuming Poisson noise induces high-frequency

instabilities whose amplitude sharply increase after a certain number of

iterations. This is expressed in high-frequency artifacts in the reconstructed

image, which can be remedied by several techniques. For instance, it is

possible to apply a 3D Gaussian ﬁlter after the image is reconstructed, in

which the intrinsic spatial resolution of the ﬁlter must correspond to the

one that we want to obtain in the image for the chosen SNR. Possibilities

consist of ﬁltering the data before applying the algorithm or even stopping

the reconstruction in a previously chosen iteration. There are several meth-

ods that estimate the most appropriate number of iterations to run [87,88],

although all clinical applications use a number of iterations empirically

determined [89].

6.3.1.2.2.2 Other Methods The Ordered Subset Expectation Maximization

(OSEM) method [90] became very popular, because it signiﬁcantly accelerates

the convergence to a ﬁnal solution relative to the MLEM algorithm, although

it is based on exactly the same principles. The fast speed of the OSEM algo-

rithm is achieved by dividing the data into subsets: in each iteration, the

estimate of the activity distribution obtained from one subset is used as the

initial estimate in the next one. The convergence process of OSEM is swifter

than MLEM; but, on the other hand, it has an increased numerical noise, a

fact that is borne out by the different solutions found by the two methods.

An improvement of the OSEM method often used today is the Attenuation

Weighted OSEM (AWOSEM) [91], which is based on a better modeling of the

image noise to achieve better image quality without signiﬁcant reductions in

the computation speed.

6.3.1.3 3D Reconstruction

In the 3D mode of PET, image reconstruction using analytical methods resorts

to the generalization to three dimensions of the methods developed for the 2D

mode [89,93]. There is, however, an important difference that raises several

problems from the practical point of view: In the 2D mode, one has 2D arrays

corresponding to the projections p

, φ, z)



z=z

of each plane separately;

whereasin the 3D mode all projections have to be simultaneously used, which

forces reconstruction to deal with the projection arrays p

, y

, φ, θ) in the

totality of their four dimensions. In this way, 3D reconstruction exhibits a new

252 Nuclear Medicine Physics

aspect relative to 2D reconstruction: there is information redundancy in the

measured data, so it is possible to reconstruct the imaged object using only

a part of the 3D mode data (e.g., it is possible to reconstruct the object using

only the nonoblique projections, i.e., those whose co-polar angle is 0

◦

, which

corresponds exactly to the 2D mode data).

Additionally, 3D reconstruction implies a drastic step-up in the amount

of data, which now includes oblique projections and makes computational

processing quite a lot more complex and lengthy. The result, however, is

images of better quality than the 2D ones using the same amount of injected

radiotracer, because one has a larger number of detected coincidences in the

3D mode. This fact justiﬁes the increasing interest in the development and

implementation of 3D image reconstruction methods in many PET systems,

in spite of their slowness.

6.3.1.3.1 Analytic Methods

Adapting the analytical methods of the 2D mode to the 3D mode is

straightforward, as the backprojection operation of parallel projections, now

two-dimensional in the coordinates (x

, y

), is still possible for each direc-

tion of space speciﬁed by the angular coordinates (φ, θ). Thus, to reconstruct

the activity distribution, one just needs to backproject all the projections

, y

, φ, θ) along all the directions in space (φ, θ). Just as for the 2D

mode, ﬁltering is still necessary to prevent the oversampling leading to

blurring effects. However, the process of backprojection in the 3D mode

has a fundamental problem that does not exist in the 2D mode: depend-

ing on the acquisition parameters and on the geometrical relation between

the imaged object and the detector, there may be parallel projections that are

incomplete, meaning that the sampling of the object volume is not uniform

(Figure 6.22) [80].

Seg. 0 Seg. +1 Seg. +2

FIGURE 6.22

(a) Axial view of planes deﬁned on a PET system in the 3D mode. The cross-section of the

parallel projections decreases as the inclination θ of the projection lines increases. (b) Drawing of

an incomplete parallel projection due to the limited FOV of the PET system.

Imaging Methodologies 253

The problem of the incomplete projections can be solved by introducing a

previous step of 2D mode reconstruction, which enables an initial estimate

of the activity distribution f

est

(x, y, z) to be obtained. This estimate is then

projected to ﬁll the regions in the (incomplete) oblique sinograms, from which

a set of complete parallel projections p

3Dcompl

, y

, φ, θ) is drawn, and to

which it is possible to apply a 3D FBP reconstruction operation with a chosen

ﬁlter. The whole process is known as 3D-Reprojection (3DRP) [89], and it is

depicted in Figure 6.23.

6.3.1.3.2 Iterative Methods

The philosophy of image reconstruction in PET using iterative methods is the

same as that depicted in Figure 6.21 for SPECT and 2D mode PET. It is a cyclic

process that involves the gradual improvement of an estimate for the activity

distribution by comparing, in each iteration, the projections generated by the

estimate with the measured projections and the computation of a correctionto

the activity distribution estimate from that comparison. The implementation

of this principle is identical for the PET 2D and 3D modes, differing only in

the dimension of the projections used (two-dimensional in the 2D mode, four-

dimensional in the 3D mode) and of the system matrix employed to obtain

the projections corresponding to the estimate of the activity distribution in

each iteration. The different iterative methods available for 3D PET include

ART, MLEM, OSEM, and others.

2D reconstruction

(FBP)

1st estimate

of the object

(noisy)

Reprojection

Estimation of data in lines of

response not measured

Complete truncated

projections

3D Reconstruction

(3D FBP)

Reconstructed

volume

f(x, y, z)

est

(x, y, z)

FIGURE 6.23

Outline of the analytical reconstruction algorithm 3DRP [92]. For the standard analytical recon-

struction algorithm 3D FBP to be correctly applied, the measured projections have to be complete,

though, in general, the oblique projections are truncated and that does not happen. Hence, the

3DRP method aims at completing the oblique projections by projecting the volume obtained from

the reconstruction of the nonoblique projections.

254 Nuclear Medicine Physics

As in SPECT and 2D-mode PET, these methods allow the inclusion of image

formation models which enable the correction of several physical factors that

degrade the image as well as incorporating a priori information about the

image one intends to reconstruct. Some examples of those effects are the cam-

era ﬁnite point spread function (PSF), the existence of scattered radiation, or

the noise embedded in the attenuation correction; whereas a priori informa-

tion may include the setting of mathematical conditions for the formed image,

such as the positivity of the activity in all points of space, or the conﬁnement

of the activity within the frontiers of the object, etc.

Although iterative methods in PET also exhibit nonlinear behavior that is

difﬁcult to foresee and interpret, their advantages over the analytical methods

largely overcome the disadvantages, which is why they are today the most

common approach for image reconstruction in the 2D mode [94–96]. In the 3D

mode, however, it is rare in clinical practice to perform image reconstruction

using pure 3D data, and either the analytical or iterative methods described

in the next section are used. This is due to matters of practical order related to

the large amount of data contained in the parallel projections p

, y

, φ, θ)

and to the fact that they have to be simultaneously dealt with. This makes

the system matrix so big that it is difﬁcult for computers to cope with, both

from the point of view of storage in memory and from the point of view of

processing speed [97]. The iterative reconstruction of 3D data ends up being

mostly employed in research, almost always with the purpose of creating

new computational methodologies capable of making the reconstruction of

3D data feasible in practice.

6.3.1.4 Rebinning Methods

The analytical or iterative reconstruction of PET data acquired in 3D mode

is signiﬁcantly more time consuming than reconstruction in 2D by the same

methods, so it is almost completely excluded from clinical practice unless

dedicated hardware is used. However, the fact that 3D mode data offer a sub-

stantially better sensitivity has prompted the search for methods that, using

all the information contained in the 3D data, would carry out the image recon-

struction in the 2D mode. The solutions found have a simple philosophy that

consists of performing the regrouping, or rebinning, of 3D data into 2D data

(Figure 6.24), adding the counts in oblique planes to the 2D planes. Assuming

that the data areorganized in the sinogram format, this operation corresponds

to the transformation

, φ, z, Δr) → s

, φ, z), (6.33)

where each speciﬁc rebinning method establishes a rule to attribute counts

from oblique directions to the transaxial planes used in the 2D mode.

For the simplest rebinning method, the Single-Slice Rebinning (SSRB) [98],

the counts in an oblique LOR with axial coordinate z are added to the transax-

ial plane with the same coordinate z independently of the LOR inclination Δr.

Imaging Methodologies 255

sinograms

– 1 sinograms

(N direct and N – 1 crossed)

Axial FOV

Rebinning

Δr = N – 1

Δr = 1

Δr = –1

Δr = 0

Δr = N – 2

Δr = –(N – 2)

Δr = –(N – 1)

FIGURE 6.24

Outline of the rebinning process. The 3D sinograms are grouped in a smaller number of 2D

sinograms corresponding only to transaxial planes.

Source on the

detector axis

(a) (b)

Source away from

the detector axis

FIGURE 6.25

Attribution of counts of oblique LORs to transaxial planes in the (a) SSRB and (b) MSRB methods.

The errorlocating the axial coordinateof activity sources that are outside the detector axis in SSRB

is substantially reduced in the MSRB method, as the counts of oblique LORs are distributed over

all the planes crossed by the LORs.

256 Nuclear Medicine Physics

The reduction of the sinograms follows

, φ, z) =



Δr

, φ, z, Δr)

,φ,z) ﬁxed

, (6.34)

implying that the counts of each oblique LOR are added to the middle plane

of the two detectors deﬁning that LOR (Figure 6.25a). The simplicity of this

rebinning rule does, however, have the downside of introducing errors in

the determination of the axial coordinates of activity sources placed away

from the detector axis [80]. The Multi-Slice Rebinning (MSRB) [99] method

is an alternative which addresses that problem and equally distributes the

number of counts in an oblique LOR by all the planes traversed by that LOR

(Figure 6.25b). The MRSB method is fast to process numerically, but it is not

very robust to the presence of noise in the data, as it displaces the activ-

ity of very oblique LORs (which are those with fewer counts) by several

planes.

These two methods are often ignored in favor of a third rebinning method,

known as Fourier Rebinning (FORE) [100]. This method is very stable to

noise and involves the numerical computation of the 2D Fourier transform

, v

, z, Δr) of each oblique sinogram s

, φ, z, Δr) collected in the 3D

mode. The FORE method is based on an approximation, known as frequency–

distance relationship [80], according to which S

, v

, z, Δr) depends only

on the activity distribution f (x, y, z) in a set of discrete points along each LOR.

The transform S

, v

, z, Δr) of each oblique sinogram (z and Δr ﬁxed) is

the data structure suffering the rebinning process, with the value of each coef-

ﬁcient of that transform being added to a sum-sinogram S

, v

, z) in the

z plane coordinate that intersects the LOR in question in the point where

the frequency–distance relationship is satisﬁed. The regrouped 2D sino-

grams s

, φ, z) are ﬁnally obtained by calculating the (two-dimensional)

inverse Fourier transform of S

, v

, z) for each plane individually. Amore

detaileddescription ofthis method, in particular aboutthe frequency-distance

relationship and the region where it is valid, can be found in [80].

Any of these rebinning methods adds numerical noise to the data, leading

to a performance after the 2D reconstruction that is worse than if pure 3D

reconstruction had been performed. However, rebinning makes the image

reconstruction task much faster; and the noise problem can be largely com-

pensated when the associated 2D reconstruction method belongs to the class

of iterative methods, because they are capable of modeling and suppressing

noise; the FORE+OSEM combination, one of the most popular methods of

reconstruction today, is a good example of this type of approach [80].

The combination FORE+AWOSEM has also become very popular, and it

is an improvement over FORE+OSEM, because it allows the reconstruction

algorithm to be applied in conditions nearer to the hypotheses on which they

are based. In fact, the FORE algorithm assumes that the data from which it

starts are corrected for the different physical effects that affect the measure-

ment. However, the carrying out of that correction destroys the Poisson

Imaging Methodologies 257

character of the data, a hypothesis which underpins the OSEM algorithm.

So, to apply the FORE algorithm in its best conditions one would have to

subsequently apply the OSEM algorithm in bad conditions, that is, knowing

that the data output from the FORE procedure is not Poisson distributed,

having a variance that does not equal the average. To avoid this problem, the

solution found with the FORE+AWOSEM procedure consists of ﬁrst apply-

ing the FORE algorithm and afterward correcting the data with a weighting

factor that tries to restore the Poisson character to the data before applying

the OSEM method. After applying OSEM, the correction is reversed, and a

ﬁnal result is eventually produced. As the name indicates, the correction in

the AWOSEM method is based on the attenuation correction coefﬁcients. The

attenuation correction is the one that affects the amplitude of the measured

values the most and due to that, it has a large inﬂuence on removing the

Poisson-like character of the data; it also has the advantage of being easy to

apply, as the attenuation factors are available, being based on the transmis-

sion measurement. Other weighting methods have been proposed, including

other corrections, such as the normalization correction [91], as well as other

criteria directly related to noise [101,102].

Thereare two other speciﬁc rebinningalternatives, named FOREX [103] and

FOREJ [104], which are more precise than the FORE method, which is only

valid up to θ angles of about 20

◦

. These methods are referred to in the liter-

ature as exact methods, in contrast to the approximated methods mentioned

earlier (SSRB, MSRB, and FORE). As the name indicates, they are based on the

frequency-distance relationship as well, but they employ additional correc-

tion terms that allow their validity to be extended to larger θ angles, although

at the cost of a longer computation time and larger implementation complex-

ity. Finally, rebinning methods appropriate for TOF-PET have been recently

put forward [105].

6.3.2 Concepts of Digital Image Processing

6.3.2.1 Image Preprocessing

Image preprocessing concerns image operations aiming at correcting or

improving the image. Point-based processing is a special type of image pre-

processing. Point-based operations are applied pixelwise without taking into

account the values of the pixels in the neighborhood. As a result of these

operations, the histogram of the gray image is changed.

6.3.2.2 Point-Based Preprocessing

In this type of preprocessing, each gray level value is changed into another

value depending on the type of preprocessing operation. Let the set of gray

level values be [0, K − 1]. Given a speciﬁc gray level value u ∈[0, K − 1],itis

changed into another gray level value v ∈[0, K − 1] depending on a function