2 of 150

Aim of Tomographic Reconstruction

It is the task of reconstructing a 2D image (object) from its 1D projections, or a 3D image (object) from its 2D projections.

What is a projection (also called tomographic projection)?

It is defined as the Radon transform of the image in a particular direction (see next slide).

3 of 150

Aim of Tomographic Reconstruction

Involves estimation of inner structural details of an object (eg: internal organs of a person) as opposed to just the surface (as in photographic imaging).

CT scan of a person’s lung

4 of 150

Radon transform

Imagine a line was drawn through the 2D image in a certain direction α, and you integrated the intensity values along that line.

Now you repeat this for lines parallel to the original one but at different offsets.

Each such summation produces a bin of the tomographic projection.

The collection of bins form a 1D array which is called the tomographic projection or the Radon transform of the object in the direction α.

https://en.wikipedia.org/wiki/Radon_transform

5 of 150

Reconstruction from the Radon Transform

The Radon transform computation can be similarly repeated for multiple angles.

The aim of tomographic reconstruction is to reconstruct the 2D image from a collection of its 1D projections in different angles.

6 of 150

Reconstruction of the 3D object is performed from its many 2D Xray projection images, each in a different known direction. This is what happens typically inside a CT machine. The reconstruction quality is better if there are more projection images.

3D object

2D XRay projection image (angle 1)

2D XRay projection image (angle 2)

2D XRay projection image (angle 3)

7 of 150

Where in practice do you see the Radon transform?

The image is a simplification of a set of real biological tissues: example, an organ/tumor surrounded by a background consisting of soft, uniform tissue. The set of tissues is bombarded with an X-Ray beam. The tumor has higher rate of absorption as compared to the surrounding tissue.

The degree of absorption of X-Rays at each point is measured by an X-Ray absorption detector. This detector produces a 1D signal whose amplitude/intensity is directly proportional to the extent of absorption.

Any point in the signal = sum of the absorptivity values across the path of a single ray in the X-Ray beam that spatially maps onto that point.

Xray-beam

detector

8 of 150

Where in practice do you see the Radon transform?

Xray-beam

detector

Sum-total of the two back-projections

Given the 1D signal (called a projection signal), we try to reconstruct the original 2D image by smearing backwards along the direction of projection. This is called as back-projection.

The 1D signal that was measured is duplicated along the columns of the image to be estimated (see the directions marked in yellow).

9 of 150

Introduction

Given projections in K different directions, we can hope to reconstruct the original image by performing back-projection along all these directions, and adding up the results.
Back-projection refers to smearing the 1D projection back across the 2D image, i.e. duplicating the 1D signal across the image in a direction opposite to the direction of projection.
The shape of the object will be approximated better and better as K increases.

10 of 150

Even with many (32) back-projections, there is a blur artifact in the reconstruction. This is called as a “halo effect”.

11 of 150

Reconstructing 3D objects

A 3D object is illuminated with X-Ray beams to produce a projection which is a 2D-image.
Changing the direction of the X-ray beam will produce another image. This set of images when back-projected will yield the 3D volume/object.
However in olden Computed Tomography (CT), each slice of the volume is measured at a time. A slice is a 2D entity obtained by cutting the 3D volume transversely through a plane parallel to the XY plane.
This allows for the employment of a smaller number of detectors at a time, for the same resolution of the measurement.

12 of 150

I₀

X-Ray Image

I0 = intensity of the X-ray beam from the source

I = intensity of the X-ray beam as measured by the detector, given by Beer’s law

13 of 150

Tomography and medical imaging

Tomography is useful in medical imaging, because different tissues (eg: tissues in the bone, lungs, etc.) attenuate the X-rays to different extents.
This allows for good medical diagnosis.

14 of 150

Given a direction ϴ, the source and detector pair move along that direction in fixed steps (i.e. variation in ρ). The distance between source and detector is constant. At each step, the source sends out an X-Ray beam onto the subject, and the projection value is recorded on the detector and stored in a computer. This process is repeated for several values of ϴ. In the end, we record a single 2D slice. Now, the subject is moved in a direction perpendicular to the plane of the source-detector pair, and another 2D slice is recorded. This is called 1^st generation CT.

Source of image: Book by Gonzalez, 3^rd edition

15 of 150

In second-generation CT, the source sends out an X-ray beam in the form of a cone (also called fan). There are more detectors needed, but fewer translations are required to record a single 2D slice (rotations are needed just as before).

In third-generation CT, the source sends out an X-ray beam in the form of a large cone. The number of detectors is long enough to cover the field of view, so no translations are needed to record a single 2D slice (only rotations are needed).

Source of image: Book by Gonzalez, 3^rd edition

16 of 150

In fourth-generation CT, there is a large number of detectors arranged in the form of a circle. The source alone rotates.

Source of image: Book by Gonzalez, 3^rd edition

17 of 150

Newer advancements

Today’s scanners follow cone-beam or helical CT geometry.
A large stack of slices (a sub-volume) is imaged via such cone beams from different directions.
The entire sub-volume is reconstructed at a time.
The object may be moved in a perpendicular fashion to image other sub-volumes.
The reconstruction of different sub-volumes will be stitched/merged together.

18 of 150

Aim of Tomography

To estimate the full 3D structure of an object from its projections.
The projections are directly measured, the 3D structure is estimated.
Applications: medical imaging, industrial applications such as fault detection in machines, observation of plant roots, remote sensing (observation of underground objects or phenomena), seismology, geology.

19 of 150

Representing Projections

The complete set of projections for several different values of the parameters ρ and ϴ gives:

This is called the Radon Transform of f. Its discrete version is:

One single projection vector is obtained with a fixed value of ϴ, but varying ρ.

Dirac delta function

Eqn. of line in “normal form”

Kronecker delta function

20 of 150

Parameterization of a direction/line

The direction of projection is denoted L, and dL is an infinitesimally small element along L. L is parameterized as follows (the “normal representation of the line”):

21 of 150

Example:

Since this is a circular object, the projections are identical in all directions. Hence WLOG, we consider theta = 0, due to which the argument to δ gets simplified

Note: there is no closed form for the Radon transform for most functions f(x,y).

22 of 150

Sinogram (a radon transform plotted as an image in a (ρ,ϴ) grid.

https://en.wikipedia.org/wiki/Radon_transform

23 of 150

Back-projection

Radon transform: obtained by sampling several different angles

Fix the angle ϴ_kand for all x and y, compute the value of ρ. Copy g(ρ, ϴ_k) to hat(f)_ϴk(x,y), which is the image obtained when you back-project along angle ϴ_k.

The back-projection operator is NOT the same as the inverse of the Radon transform! So this does not yield back the true signal f(x,y), but the signal f(x,y) blurred with the kernel (x²+y²)^-0.5. More on this a few slides down, when we do filtered back-projection.

24 of 150

The blur is a painful consequence of (1) discretization of the angle ϴ, and (2) the inherent blurring with the kernel (x²+y²)^-0.5. These images are reconstructed at 0.5 degree changes in ϴ. How do we get rid of this blur? Wait for a few slides!

25 of 150

Fourier transform of the Radon Transform

The Radon transform is given as:

Its 1D Fourier transform w.r.t. ρ (keeping ϴ fixed to some value) is given by:

G(μ, ϴ) is the Fourier transform of the projection of f(x,y) along some direction ϴ.

26 of 150

Fourier transform of the Radon Transform

Its 1D Fourier transform w.r.t. ρ (keeping ϴ fixed) is given by:

The RHS of this equation is a slice of the 2D Fourier transform of f(x,y), i.e. F(u,v), along the angle ϴ in the frequency plane, and passing through the origin

This equation above is called the Projection Slice Theorem or the Fourier Slice Theorem. It states that the Fourier transform of a projection of the 2D object along some direction ϴ (i.e. G(μ, ϴ)) is equal to a slice of the 2D Fourier transform of the object along the same direction ϴ (in the frequency plane), passing through the origin.

27 of 150

The Projection Slice Theorem or the Fourier Slice Theorem states that the following two are equivalent:

Project a 2D object along a certain direction d. Take its 1D Fourier Transform called as F1.
Compute the 2D Fourier transform of the same object. Take a slice of this Fourier transform along a direction parallel to d (but in the frequency plane). Call this slice as F2.

Now F1 = F2.

Source of image: https://en.wikipedia.org/wiki/Projection-slice_theorem

28 of 150

Filtered Back-projection

Consider the 2D inverse Fourier transform of F(u,v), giving us:

Consider u = μ cos(Ѳ), v = μ sin(Ѳ). Then:

Note: we are doing a change of variables from (u,v) to (μ,Ѳ). Hence du dv = μ dμ dѲ.

29 of 150

Filtered Back-projection

By projection slice theorem, this becomes:

Further simplification will give the following (see next slide)

30 of 150

Filtered Back-projection

31 of 150

Filtered Back-projection

This is a 1D Inverse Fourier Transform with an added term |μ| (a ramp filter). But this function is not integrable as |μ| grows unboundedly. Hence the inverse Fourier transform does not exist!

|μ|

32 of 150

Filtered Back-projection versus back-projection

This is a 1D Inverse Fourier Transform with an added term |μ| (a ramp filter). But this function is not integrable as |μ| grows unboundedly. Hence the inverse Fourier transform does not exist!

Now suppose that additional term |μ| were absent. We would then obtain simple back-projection:

By Fourier convolution theorem

33 of 150

Filtered Back-projection algorithm: version 1

The X-Ray or CT machine gives you projections (i.e. Radon transforms) of the data in different directions.
Perform simple back-projection to get a 2D image.
Find the 2D Fourier transform of the back-projected 2D image, multiply it with (u²+v²)^^0.5 and find the 2D inverse Fourier Transform of the product.

34 of 150

Filtered Back-Projection

In practice, the ramp filter is multiplied in the frequency domain by a windowing function (equivalent to convolution with the corresponding kernel in the spatial domain).

It has the advantage of controlling some amount of noise as it band-limits the ramp filter.

The simplest windowing function is a rect filter in the frequency domain.

This formula is called the Ramachandran Lakshminarayanan (Ram-Lak) filter.

http://www.pnas.org/content/68/9/2236

35 of 150

Other windowing functions can also be used!

(eg: Hamming window, which is a truncated cosine). These windowing functions allow lesser ringing artifacts as compared to the rect windowing function.

Ram-Lak filter: ramp multiplied by rect (both in frequency domain)

Freq. domain

36 of 150

Ram-Lak filter with Hamming window

Ram-Lak filter

37 of 150

Filtered Back-projection algorithm: version 2 (faster)

Compute the 1D Fourier transform of each projection of the object (the projections are measured by the X-Ray or CT machine).
Multiply each Fourier transform by the ramp filter |μ| multiplied by a windowing function (such as rect, or Hamming).
Obtain the 1D inverse Fourier transform of the result.
Sum up over all such results (one each per projection angle) from the previous step to give the final image.

38 of 150

Using backprojection – 180 angles

Using filtered backprojection – 180 angles (Ram-Lak filter)

39 of 150

Tomography in MATLAB

The 2D Radon transform can be computed in MATLAB using the function radon.

Filtered back-projection with a choice of various filters – Ram-Lak, Hamming, Hanning, etc., is implemented in MATLAB in a function called iradon.

3D transforms generally work slice by slice for parallel beam projections and are implemented in toolboxes written by various authors: https://web.eecs.umich.edu/~fessler/code/index.html

40 of 150

Tomography and Compressed Sensing

In most tomography applications, the number of angles of projection is limited due to cost, energy and health considerations.

So the problem is of practical interest in the “angle starved” case.

Hence tomography can be considered a compressed sensing problem.

41 of 150

Tomography and Compressed sensing

The filtered backprojection algorithm does not exploit one important property of images – their sparsity or compressibility in standard bases such as DCT.

Instead one can use CS principles to frame the tomographic reconstruction problem as follows:

Vector created by concatenating 1D projections in various angles

Forward model matrix – representing Radon transform computation in those angles

42 of 150

Tomography and Compressed sensing

Instead one can use CS principles to frame the tomographic reconstruction problem as follows:

In practice (say while using ISTA), the R matrix is implemented in function handles that invoke the radon function in MATLAB.

R^T is implemented using iradon.

Vector created by concatenating 1D projections in various angles

Forward model matrix – representing Radon transform computation in those angles

43 of 150

Tomography and Compressed sensing

The R matrix is too large to store in memory.

Likewise the U matrix is never allocated (it is too large to store in memory).

Instead a function handle (idct2 or idwt2) is invoked.

See more about MATLAB function handles here.

Function handles for R, U offer computational speed and do not require too much memory.

44 of 150

Tomography and compressed sensing

The reconstruction results with CS are generally seen to be consistently superior to those by FBP in the angle-starved case.

Next slide shows comparative results: original 240 x 240 cross-sectional images of a walnut (top row), FBP reconstructions (2^nd row), CS reconstructions (3^rd row) and coupled CS reconstructions (4^th row – explained later). All reconstructions with 20 angles.

46 of 150

Original image, size ~ 200 x 200

Reconstructions using FBP (left) and CS (right) using 15 angles (top row), 30 angles (middle row) and 40 angles (bottom row)

47 of 150

Great results, but…

The R matrix as defined earlier on, however, is not known to obey the RIP or exhibit low coherence with standard bases.

Hence the theoretical treatment of CS for tomography has not been fully established.

48 of 150

Augmenting CS for tomography

In multi-slice tomographic reconstruction, one can make use of additional redundancy in the data – for instance, the difference between two consecutive slices is sparse.

This can be used to improve the tomographic reconstruction quality for the same number of measurements – if one chooses different angles for different slices.

49 of 150

Augmenting CS for tomography

The objective function (for two consecutive slices) is as follows:

x₁

x₂

Here x1 and x2 represent two consecutive slices of an organ (each slice is a 2D image), and y1 and y2 represent their tomographic projections expressed as 1D vectors.

50 of 150

Augmenting CS for tomography

See the walnut dataset 3-4 slides earlier for a sample reconstruction to see the advantage of this model.

Note that R₁ and R₂ denote the Radon-based forward models for different angle sets (the number of angles in the two sets may or may not be equal).

51 of 150

Tomography under unknown angles

52 of 150

Problem statement

The previous algorithms for tomographic reconstruction assumed that the angles of Radon projection were accurately known.

In certain applications, this assumption is surprisingly invalid.

This is called tomography under unknown angles.

53 of 150

Problem statement: Applications

Application 1: Patient motion during CT scanning

Application 2: Moving insect tomography

Application 3: Cryo-electron microscopy/tomography

54 of 150

Problem statement: Applications

Application 3: Cryo-electron microscopy

In this, one collects multiple (nearly) identical samples of a structure (such as a virus/macromolecule) which we wish to image. In fact, we are interested in their complete 3D structure.

Each slide contains thousands of macromolecules (i.e. samples) packed in a substrate such as ice.

A tomographic projection is obtained by passing an electron-ray beam through all macromolecules, through some angle (usually vertically downward).

55 of 150

Problem statement: Applications

The electron beam usually destroys the sample, and hence another tomographic projection of a different sample (containing virus particles of the same type) is acquired.

The problem is that each virus or macromolecule copy will be oriented randomly, and all the orientations are unknown!

To make matters worse, the low power of the electron beam produces measurements that are extremely noisy.

In such applications, however several hundred or even thousand projections (all under unknown angles) are acquired.

56 of 150

https://en.wikipedia.org/wiki/Cryogenic_electron_microscopy

Such an image is called a micrograph. It contains tomographic projections of each macromolecule, against a background containing only ice. Each such projection image is called a particle.

57 of 150

https://med.nyu.edu/skirball-lab/stokeslab/phi12.html

58 of 150

https://ki.se/en/research/core-facility-for-electron-tomography-0

59 of 150

Steps in Cryo-em microscopy

Particle picking from noisy micrographs

In some algorithms, similar particles are clustered and averaged to reduce noise

Given the series of particle images, we then seek to solve jointly for the angles of projection and the underlying structure

60 of 150

Cryo-electron microscopy wins the Nobel!

Nobel in Chemistry in 2017

More details here below:

https://www.nobelprize.org/nobel_prizes/chemistry/laureates/2017/advanced-chemistryprize2017.pdf

61 of 150

Two Different Approaches

Moment-based approach

Ordering-based approach

Approach using dimensionality reduction (similar to ordering-based approach).

62 of 150

Moment-based approach

63 of 150

Moment-based approach

We shall restrict ourselves to 2D images and 1D tomographic projections although the theory is extensible to 3D images (and their 2D projections)

The moment of order (p,q) of an image f(x,y) is defined as follows:

64 of 150

Moment-based approach

The moment of order (p,q) where k=p+q of an image f(x,y) is defined as follows:

Note that there can exist multiple pairs of (p,q) which sum up to k, and these are all called order k image moments.

65 of 150

Moment-based approach

The order n moment of a tomographic projection at angle θ is defined as follows:

Substituting the definition of P_θ(s) into M_θ(n):

66 of 150

Moment-based approach: image moments and projection moments

Using the binomial theorem, we have:

We will use this to derive a neat relationship between the tomographic projection moments and the image moments!

See next slide.

67 of 150

Moment-based approach: image moments and projection moments

Image moment of order

(n-l,l)

68 of 150

Moment-based approach: image moments and projection moments

Substituting n = 0, with measurements at one angle.

Substituting n = 1, with measurements at two angles.

69 of 150

Moment-based approach: image moments and projection moments

Substituting n = k, with measurements at k+1 different angles.

70 of 150

Moment-based approach: image moments and projection moments

These equations are called the Helgason-Ludwig consistency conditions (HLCC), and they give relations between image and projection moments.

One can prove that the matrix A is invertible if and only if the projections are acquired at k+1 distinct angles. In fact, unique k+1 angles are necessary and sufficient for estimation of the image moments of order 0 through to order k.

71 of 150

Tomography under unknown angles: rotation ambiguity

In the tomography under unknown angles problem, we would know neither the image moments nor the angles of acquisition.

In such a case, the underlying image can be obtained only up to an unknown rotation.

To understand why, see the next slide.

72 of 150

θ₁

θ₂

θ₃

In the first case you took 3 tomographic projections (denoted p1,p2,p3) of an object respectively at three angles θ₁, θ₂, θ₃

φ+θ₁

φ+ θ₂

φ+ θ₃

In the second case you took tomographic projections (denoted q1,q2,q3) of a version of the same object but rotated by φ respectively at three angles φ +θ₁, φ +θ₂, φ +θ₃

In both cases, the corresponding projections will be identical, ie p1 = q1, p2 = q2, p3 = q3! The parameter φ will always be indeterminate – but this is not a problem in most applications

73 of 150

A surprising result

Given tomographic projections of a 2D image in 8 or more distinct and unknown angles, the image moments of order 1 and 2, as well as the angles can be uniquely recovered – but up to the aforementioned rotation ambiguity.
This result is true for almost any 2D image (i.e. barring a set of very rare “corner case” images).
This result was proved in 2000 by Basu and Bresler at UIUC in a classic paper called “Uniqueness of tomography with unknown view angles”.
In an accompanying paper called “Feasibility of tomography with unknown view angles”, they also proved that these estimates are stable under noise.
The proof of the theorem and the discussion of the corner cases is outside the scope of our course.

74 of 150

A surprising result

In other words, systems of equations of the following form have a unique solution in the angles and the image moments, but modulo the rotation ambiguity:

Image moments

Projection moments

Column vector of image moments of order n

This is the n-th row of a matrix and it represents the linear combination coefficients for moments of order n and at angle θ_i.

75 of 150

Tomography under unknown angles: algorithm

We can now build an algorithm for the aforementioned problem.
Minimize the following objective function in an alternating fashion:

Start with a random initial angle estimate
(&) Compute the image moments by least squares inversion.

Q=#projections

N = max. order of image moments

76 of 150

Tomography under unknown angles: algorithm

($) Next, do an independent brute force search over each angle θ_i.
Repeat the earlier two steps (marked as & and $) until convergence (i.e. Until the estimates of either thetas or the image moments do not change too much)
Perform a multi-start strategy for the best possible results – since this cost function is highly nonconvex.

* For every value of θ_i sampled from 0 to 180, determine the image moments using that value, and hence determine the value of E.

Choose the value of θ_icorresponding to the least value of E.

77 of 150

Tomography under unknown angles: algorithm

Remember: these angles can be estimated only up to a global angular offset φ which is indeterminate.
Following the angle estimates, the underlying image can be reconstructed using FBP.

78 of 150

Image source: Malhotra and Rajwade, “Tomographic reconstruction with unknown view angles exploiting moment-based relationships”

https://www.cse.iitb.ac.in/~ajitvr/eeshan_icip2016.pdf

79 of 150

Summary of moment-based method

Advantage: simple, makes direct use of the consistency conditions and works even if the number of angles is small
Disadvantage: moments are highly sensitive to noise

80 of 150

Ordering-based approach

81 of 150

Approach – the uniform distribution

In this approach, the tomographic projections need to be “sorted” – i.e. arranged in order of increasing angles.
The method relies on the availability of a large number of tomographic projections – all under unknown angles – but sampled independently from a uniform distribution on a unit circle.
Given sufficiently many projections, the angles can be assumed to be equi-spaced.

82 of 150

Approach – the uniform distribution

The problem reduces to projection association – i.e. matching each projection to one of the angles, sampled evenly from the unit circle.
If the angles are spaced sufficiently closely, and the variation of the projections w.r.t. angle, is smooth, then we can employ a nearest neighbour type of heuristic.
This algorithm is summarized on the next slide.

83 of 150

Algorithm

Let us denote the list of available projections as L.
Let q_i be the i-th tomographic projection and s_i be the i-th tomographic projection in the sorted list.
WLOG, define the orientation of q₁to be 0 and set s₁=q₁. Delete s₁ from L.

Find the nearest neighbour of s₁ in L as follows:

Repeat the earlier two steps with s₂.
Continue these steps until L is empty.

84 of 150

Output of the algorithm

At the end of this algorithm, each projection will be associated with an angle at a grid point on the even sampling pattern from 0 to π.
Remember: these angles can be estimated only up to a global angular offset φ which is indeterminate.
Following the angle estimates, the underlying image can be reconstructed using FBP.

85 of 150

Comments on algorithm

Note that essentially, each projection is being assigned an angle ranging from 0 to π in steps of π /Q where Q is the number of projections.

The aforementioned algorithm is an approximate solution for the well-known travelling salesperson problem in computer science.

The role of q_reverse is to flip projections whose angles would have been in [π,2π] because g(s,θ+ π)=g(-s,θ) as per the definition of the Radon transform. If the angles were only from 0 to π, this step wouldn’t be required.

86 of 150

Comments on algorithm

The accuracy of the angle estimate increases as the number of angles increases.
Consider the k-th order statistic θ_(k), i.e. the k-th largest angle in the sorted sequence.
The mean value of this order statistic is πk/Q where Q is the number of angles, assuming uniform distribution of the angles.
The variance of this order statistic can be proved to be O(1/Q²).
This method also relies on the fact that the angles were uniformly distributed, otherwise the angle grid points will need to be chosen differently – as per the order statistics of the other distribution.

87 of 150

Segway: order statistics

Consider a set of N numbers {x₁, …, x_N} drawn from a distribution F(.).
Now suppose we sort in increasing order them to produce the sequence x₍₁₎, x₍₂₎,..., x_(N).
Then x₍₁₎ is called the zero-th order statistic (minimum value), x₍₂₎ is called the first order statistic, and so on.
The order statistics are random variables in their own right, and hence possess their own distribution, expectation, variance, etc.
The expected value of the k-th order statistic for a uniform distribution is k/N given N numbers from the uniform distribution.
For more details, see https://en.wikipedia.org/wiki/Order_statistic

88 of 150

Unknown shifts

Usually, the projection angles are unknown, but so are the projection shifts – due to uncontrolled translation of the object being imaged.
However the unknown shifts can be solved for using a pre-processing step that assumes that the object is surrounded by vacuum.
The tomographic projections can be normalized such that their centroid lies at the origin.
The projection ordering step is performed later.

89 of 150

Using dimensionality reduction

90 of 150

Description of approach

In principle, this method is similar to the ordering-based approach.
However instead of doing pairwise searches, a dimensionality reduction step is applied.
Each tomographic projection q_i (under unknown angles) can be regarded as a high-dimensional vector. There are Q tomographic projections in all.
A dimensionality reduction method is applied to reduce their dimensionality to two (details later), i.e. q_i is mapped to y_i=(φ_1i, φ_2i).

91 of 150

Description of approach

Consequently the tomographic projection q_i will get mapped to the angle:

These angles are then used for creating an ordering.
Given the ordering, the order statistics of the angle distribution is used to assign angles, prior to FBP-based reconstruction.
This is similar to the nearest neighbor based technique.

92 of 150