Skip to content Skip to navigation

Connexions

You are here: Home » Content » Discretizing Linear Inverse Problems

Navigation

Recently Viewed

This feature requires Javascript to be enabled.
 

Discretizing Linear Inverse Problems

Module by: Justin Romberg. E-mail the author

In many real-world applications, the signal or image we are measuring is a function of a continuous variable (or variables for images). Of course, if we are going to reconstruct the signal/image on a computer, our answer will ultimately be discrete. In this module, we discuss a general way to discretize linear inverse problems using a basis representation. We will start with the particular example of 2D tomography (“reconstruction from projections”), but the framework will be easy to generalize.

The Radon Transform

In the 2D tomographic reconstruction problem, the image f(s,t)f(s,t) we wish to acquire is sampled using line integrals. We can parameterize a line ℓ→ℓ→ using an offset rr and an angle θθ as shown below:

Figure 1
(a)
Figure 1(a) (line-sketch1.png)
(b)
Figure 1(b) (line-sketch2.png)

The line ℓ→ℓ→ is the set of points obeying a linear constraint:

ℓ → = { ( s , t ) : s cos θ + t sin θ = r } ℓ → = { ( s , t ) : s cos θ + t sin θ = r }
(1)

The integral of f(s,t)f(s,t) along ℓ→ℓ→ is given by

R r , θ [ f ] = ∫ f r - t sin θ cos θ , t d t | θ | ≤ π / 4 ∫ f s , r - s cos θ sin θ d s π / 4 < | θ | ≤ π / 2 R r , θ [ f ] = ∫ f r - t sin θ cos θ , t d t | θ | ≤ π / 4 ∫ f s , r - s cos θ sin θ d s π / 4 < | θ | ≤ π / 2
(2)

Of course, these expressions are equal to one another except when θ=0,π/2θ=0,π/2. Note also that the measurements are unique only over a range of ππ, as Rr,θ+π[f]=R-r,θ[f]Rr,θ+π[f]=R-r,θ[f]. It is sometimes convenient to write the line integral as a 2D integral of f(s,t)f(s,t) against a delta ridge:

R r , θ [ f ] = ∫ ∫ f ( s , t ) δ ( s cos θ + t sin θ - r ) d s d t , R r , θ [ f ] = ∫ ∫ f ( s , t ) δ ( s cos θ + t sin θ - r ) d s d t ,
(3)

where δ(·)δ(·) is the Dirac delta function.

The collection of all such line integrals {Rr,θ[x],θ∈[0,π],r∈R}{Rr,θ[x],θ∈[0,π],r∈R} is called the Radon transform of f(s,t)f(s,t). The radon transform is itself a continuous function of two variables. The figure below show an illustrative example: on the left, we see Rr,θRr,θ of a test image as a function of rr for two different fixed values of θθ. On the right is the collection Rr,θRr,θ as a function of both rr and θθ.

Figure 2: Left: Rr,π/4[f]Rr,π/4[f] and Rr,π[f]Rr,π[f] as a function of rr, where f(s,t)f(s,t) is the Shepp-Logan phantom. Right: The Radon transform of the phantom. The rows are indexed by rr and the columns by θθ (in degrees).
(a)
Figure 2(a) (phantom_projections.png)
(b)
Figure 2(b) (phantom_sinogram.png)

Reconstruction from a discrete set of line integrals

Given measurements y[m]y[m], m=1,...,Mm=1,...,M corresponding to line integrals at different different offsets rmrm and angles θmθm (i.e. a finite set of samples of the Radon transform), which have possibly been corrupted by noise, we would like to estimate the underlying image f(s,t)f(s,t). If the measurements are dense in (r,θ)(r,θ) space, the natural approach to this problem is to use filtered backprojection. Our focus here will be setting this problem up as finite linear inverse problem

y = A x + noise , y ∈ R M , x ∈ R N y = A x + noise , y ∈ R M , x ∈ R N
(4)

so that it can be attacked with the general set of tools for solving such problems (e.g. least-squares).

We start by choosing a finite-dimensional space VV in which to perform the reconstruction that comes equipped with a set of NN basis vectors {ψγ(s,t)}{ψγ(s,t)}. We will use the general index γ∈Γγ∈Γ where ΓΓ is a set of size NN as, depending on the basis, it may be convenient to index the basis in different ways (i.e. by integers, pairs of integers over the same range, pairs of integers over different ranges, etc.).

For example, if f(s,t)f(s,t) is non-zero only for (s,t)∈[0,1]2(s,t)∈[0,1]2, we might take our reconstruction space VV to be the set of all “pixellated” images — images that are piecewise-constant on squares of side length 1/n1/n for some integer nn. A natural basis for this space is the set of indicator functions on these squares:

ψ j , k ( s , t ) = 1 s ∈ [ j / n , ( j + 1 ) / n ] , t ∈ [ k / n , ( k + 1 ) / n ] 0 otherwise ψ j , k ( s , t ) = 1 s ∈ [ j / n , ( j + 1 ) / n ] , t ∈ [ k / n , ( k + 1 ) / n ] 0 otherwise
(5)

Using our general index notation, we can write any f(s,t)∈Vf(s,t)∈V as

f ( s , t ) = ∑ γ ∈ Γ x [ γ ] ψ γ ( s , t ) , f ( s , t ) = ∑ γ ∈ Γ x [ γ ] ψ γ ( s , t ) ,
(6)

where Γ={(j,k):j,k=0,1,...,n-1}Γ={(j,k):j,k=0,1,...,n-1} with size N=n2N=n2, and the x[γ]∈Rx[γ]∈R are the basis expansion coefficients, which are, in this case, the pixel values. (Another natural basis for VV would be the two-dimensional Haar basis we encountered earlier). The point is that knowing the discrete set of coefficients x[γ]x[γ] for all γ∈Γγ∈Γ is the same as knowing the continuous-space function f(s,t)f(s,t).

We can also write the measurements of an f(s,t)∈Vf(s,t)∈V in terms of the basis functions:

y [ m ] = R r m , θ m ∑ γ ∈ Γ x [ γ ] ψ γ ( s , t ) = ∑ γ ∈ Γ x [ γ ] R r m , θ m ψ γ ( s , t ) (since R r , θ [ · ] is linear) = ∑ γ ∈ Γ A [ m , γ ] x [ γ ] where A [ m , γ ] = R r m , θ m ψ γ ( s , t ) , y [ m ] = R r m , θ m ∑ γ ∈ Γ x [ γ ] ψ γ ( s , t ) = ∑ γ ∈ Γ x [ γ ] R r m , θ m ψ γ ( s , t ) (since R r , θ [ · ] is linear) = ∑ γ ∈ Γ A [ m , γ ] x [ γ ] where A [ m , γ ] = R r m , θ m ψ γ ( s , t ) ,
(7)

which can be written in more compact form as

y = A x . y = A x .
(8)

The entries of the M×NM×N matrix AA contain the results of each of the MM measurements functionals Rrm,θm[·]Rrm,θm[·] applied to each of the NN basis functions ψγ(s,t)ψγ(s,t), the NN-vector xx contains the expansion coefficients for f(s,t)f(s,t) in the basis {ψγ}{ψγ}, and yy contains the MM measurements. This is illustrated in the figure below. As we can see, not too many of the ℓ→mℓ→m pass through a given pixel, meaning that the matrix AA will be very sparsely populated.

Figure 3: Left: A sketch of one of the basis functions ψγ(s,t)ψγ(s,t) from the discussion above. Right: The entries of AA in the column indexed by γγ will be the result of measuring the basis function ψγ(s,t)ψγ(s,t): A[m,γ]=Rrm,θm[ψγ]A[m,γ]=Rrm,θm[ψγ].
(a)
Figure 3(a) (pixel-sketch1.png)
(b)
Figure 3(b) (pixel-sketch2.png)

Of course, the true underlying image will in general not lie in the chosen finite-dimensional subspace VV. This means that even when there is no measurement noise, there will still be some inherent error in our calculations. But solving Equation 8 will in some sense find the member of VV that best explains the measurements that have been observed. If the true image can be closely approximated by a member of VV, then we will not lose much through this discretization. A major consideration in choosing the space VV is how well we can use it to approximate images we expect to encounter.

General linear operators

The technique above can be very naturally generalized to different kinds of measurement operators that map signals of a continuous variable(s) into RMRM. In general, the measurements y[m]y[m] will consist of inner products of f(s,t)f(s,t) against different “test functions” φm(s,t)φm(s,t) for m=1,...,Mm=1,...,M:

y [ m ] = 〈 f , φ m 〉 = ∫ ∫ f ( s , t ) φ m ( s , t ) d s d t . y [ m ] = 〈 f , φ m 〉 = ∫ ∫ f ( s , t ) φ m ( s , t ) d s d t .
(9)

In the tomography example above, we took

φ m ( s , t ) = δ ( s cos θ m + t sin θ m - r m ) φ m ( s , t ) = δ ( s cos θ m + t sin θ m - r m )
(10)

(see Equation 3 above). Just as before, for f(s,t)∈Vf(s,t)∈V we can write

y [ m ] = 〈 ∑ γ ∈ Γ x [ γ ] ψ γ , φ m 〉 = ∑ γ ∈ Γ x [ γ ] 〈 ψ γ , φ m 〉 , y [ m ] = 〈 ∑ γ ∈ Γ x [ γ ] ψ γ , φ m 〉 = ∑ γ ∈ Γ x [ γ ] 〈 ψ γ , φ m 〉 ,
(11)

and so

y = A x , y = A x ,
(12)

where the M×NM×N matrix AA has entries

A [ m , γ ] = 〈 ψ γ , φ m 〉 = ∫ ∫ ψ γ ( s , t ) φ m ( s , t ) d s d t . A [ m , γ ] = 〈 ψ γ , φ m 〉 = ∫ ∫ ψ γ ( s , t ) φ m ( s , t ) d s d t .
(13)

Notice that since the entries of AA do not depend on ff, they can be pre-computed.

Content actions

Download module as:

PDF | EPUB (?)

What is an EPUB file?

EPUB is an electronic book format that can be read on a variety of mobile devices.

Downloading to a reading device

For detailed instructions on how to download this content's EPUB to your specific device, click the "(?)" link.

| More downloads ...

Add module to:

My Favorites (?)

'My Favorites' is a special kind of lens which you can use to bookmark modules and collections. 'My Favorites' can only be seen by you, and collections saved in 'My Favorites' can remember the last module you were on. You need an account to use 'My Favorites'.

| A lens I own (?)

Definition of a lens

Lenses

A lens is a custom view of the content in the repository. You can think of it as a fancy kind of list that will let you see content through the eyes of organizations and people you trust.

What is in a lens?

Lens makers point to materials (modules and collections), creating a guide that includes their own comments and descriptive tags about the content.

Who can create a lens?

Any individual member, a community, or a respected organization.

What are tags? tag icon

Tags are descriptors added by lens makers to help label content, attaching a vocabulary that is meaningful in the context of the lens.

| External bookmarks