Automated Feature Extraction with Machine Learning and Image Processing

PD Stefan Bosse

University of Siegen - Dept. Maschinenbau
University of Bremen - Dept. Mathematics and Computer Science

1 / 45

PD Stefan Bosse - AFEML - Module Y:

X-ray Imaging: Computer Tomography

Principles of Computer Tomography

From Projections to reconstruction of object slices. Algorithms and beyond...

Quality, Noise, artifacts, and other issues with CT reconstruction

2 / 45

PD Stefan Bosse - AFEML - Module Y: Further Readings

Further Readings

  1. J. Baruchel, Jean-Yves, B. Eric, M. P. Merle, and G. Peix, X-Ray Tomography in Material Science. Hermes, 2000.
  2. M. Soleimani and T. Pengpen, “Introduction: a brief overview of iterative algorithms in X-ray computed tomography,” Phil.Trans.R.Soc.A, vol. 373, 2015.
  3. R. Cierniak, X-Ray Computed Tomography in Biomedical Engineering. Springer, 2011.
  4. BDG, BDG – Richtlinie P 203: Porositätsanalyse und -beurteilung mittels industrieller Röntgen- Computertomographie (CT). BDG, Hansaallee 203, 40549 Düsseldorf, 2019.
3 / 45

PD Stefan Bosse - AFEML - Module Y: Computer Tomography (CT)

Computer Tomography (CT)

Find an image I(x,y) from a set of rotated line projections p(s,φ)

More general:

Find an image set (volume) V(x,y,z) from a set of rotated image projections P(x,y,φ)

Definitions:

  • Projections: Input data (intensity images, line profiles, typically a radial projection set)
  • Volume: Output data as a set of images forming a 3-dim (or 2-dim) cartesian space
  • Voxel: A discrete 3-dim (or 2-dim) element of the output volume/image set (Voxel size typically same for all coordinates)
4 / 45

PD Stefan Bosse - AFEML - Module Y: Computer Tomography (CT): Examples

Computer Tomography (CT): Examples

Baruchel,2000 (Left) 3D rendered view of a tomographic image of a composite material with 400 yon glass balls inside an organic matrix; Voexl size 42 μm (Right) 3D rendered view of a tomographic image of an aluminium foam (density 0.06); Voxel size: 150 μm

5 / 45

PD Stefan Bosse - AFEML - Module Y: CT Beam Geometries

CT Beam Geometries

Cierniak, 2011

Shapes of X-ray beams used in CT scanner projection systems: (a) a parallel beam of radiation, (b) a fan beam of radiation, (c) a beam in the form of a cone

6 / 45

PD Stefan Bosse - AFEML - Module Y: CT Geometries

CT Geometries

Cierniak, 2011 (a) Darkening of a photographic film by X-rays ⇒ Inverse Attenuation / Intensity (b) Obtaining one-dimensional projections using a parallel beam of X-rays (c) Projections carried out at an angle α

7 / 45

PD Stefan Bosse - AFEML - Module Y: CT Reconstructions: Basics

CT Reconstructions: Basics

CT bases on projections. A photo (or image) is a projection of a 3-dim object onto a 2-dim plane!

Zvolský, 2014 Example: Two trees in a park, make 2 pictures from east and south, try to create a map of the park.

8 / 45

PD Stefan Bosse - AFEML - Module Y: CT Reconstructions: Basics

CT Reconstructions: Basics

Zvolský, 2014 Other configuration: If you see two separate trees on both views, can you uniquely reconstruct the map of trees? Here you cannot reconstruct the position of both trees.

If we take another picture at 45◦, we are able to solve the ambiguity.

9 / 45

PD Stefan Bosse - AFEML - Module Y: CT Reconstructions: Basics

CT Reconstructions: Basics

We now consider line projections only for the sake of simplicity. There are projections p(s,φ) at angle φ with s as the coordinate on detector, which is a line integral of a photo.

Zvolský, 2014 (Left) Projection p(s) the same for any φ (Right) Projection p(s,φ) depends on orientation
10 / 45

PD Stefan Bosse - AFEML - Module Y: CT Reconstructions: Basics

CT Reconstructions: Basics

Projections are angle dependent.

A simple example should demonstrate this: A source point on the y axis is viewed under different angles φ.

  • The location s of the spike on the 1D detector is given by:

s=rsinφ

where r is the distance of the point from the origin (measuring system) and φ the viewing angle.

  • The projection p(s,φ) in the s-φ-coordinate system is a sine function.
11 / 45

PD Stefan Bosse - AFEML - Module Y: CT Reconstructions: Basics

CT Reconstructions: Basics

Projections and calculation of the projected point position under angle φ

12 / 45

PD Stefan Bosse - AFEML - Module Y: CT Reconstructions: Basics

CT Reconstructions: Basics

A sinogram is a representation of the projections on the s-φ plane.

A sinogram combines all line profiles in a diagram providing a s-φ coordinate system (semi-polar). The point example creates a sine wave diagram.

13 / 45

PD Stefan Bosse - AFEML - Module Y: Analytic Image Reconstruction

Analytic Image Reconstruction

  1. Backprojection
  2. Filtered Backprojection

What is wrong with analytic image reconstruction from projections? Why is filtering required, and why is a filter never a good idea?

These methods are based on an analytical expression of the inversion of the Radon transform.

14 / 45

PD Stefan Bosse - AFEML - Module Y: Projection

Projection

  • Create an example projection set

We create four projections from an object consisting of 4 different densities (μ values) by calculating the sum of the contributions along a line of response (LOR)

15 / 45

PD Stefan Bosse - AFEML - Module Y: Backprojection

Backprojection

  • Placing a value of p(s,φ) back into the position of the appropriate LOR
  • But the knowledge of where the values came from was lost in the projection step
  • The best we can do is to place a constant value into all elements along the line

Summing up all p(φ) along a line of response (LOR)

16 / 45

PD Stefan Bosse - AFEML - Module Y: Backprojection

Backprojection

Get back the density distribution via I(x,y) from the projections...

17 / 45

PD Stefan Bosse - AFEML - Module Y: Backprojection

Backprojection

Sum all up all projection values under respective angles.

18 / 45

PD Stefan Bosse - AFEML - Module Y: Backprojection

Backprojection

Subtract the total projection sum Σ=a+b+c+d from all backprojected entries.

19 / 45

PD Stefan Bosse - AFEML - Module Y: Backprojection

Backprojection

We are done. In theory ...

Divide by number of projections −1 = 3

20 / 45

PD Stefan Bosse - AFEML - Module Y: Backprojection

Backprojection

With few views. But with a high number of views and continuous distributions the following is happening...

Many angles → Tall and broadened spike at the location of the point source

21 / 45

PD Stefan Bosse - AFEML - Module Y: Backprojection

Backprojection

Baruchel, 2000 (a) Original object (b) Some projections (c) Backprojection in a point of the object (d) sinogram or set of projections over 180° (e) Backprojection of the sinogram

22 / 45

PD Stefan Bosse - AFEML - Module Y: Backprojection

Backprojection

The backprojected image compared with original object is blurred. As a result of the backprojection process, each pixel contains information about what the object really contains at the pixel location, but this information is added to a blurred version of the rest of the object.

23 / 45

PD Stefan Bosse - AFEML - Module Y: Backprojection

Backprojection

The backprojected image compared with original object is blurred. As a result of the backprojection process, each pixel contains information about what the object really contains at the pixel location, but this information is added to a blurred version of the rest of the object.

An exact mathematical correction of the backprojection smoothing effect can be done by an appropriate pre-filtering of the projections, as in the Filtered Backprojection (FBP) algorithm.

This can be demonstrated based on Fourier considerations.

24 / 45

PD Stefan Bosse - AFEML - Module Y: Radon Transformation

Radon Transformation

  • The Rdaon transformation is the base function for back projection, and given for a distribution f(x, y) by:

p(s,ϕ)=∫∞−∞f(x,y)δ(xcosϕ+ysinϕ−s)dxdy

  • Due to the δ function the integrand is zero except on the Line L(s,φ)

  • The backprojected image is given by the integration over 189° (no more information in the other half, really?):

b(x,y)=∫π0p(s,ϕ)∣s=xcosϕ+ysinϕdϕ

  • The image b is blurred by: b(x,y)=f(x,y)/√(x2+y2)
25 / 45

PD Stefan Bosse - AFEML - Module Y: Radon Transformation

Radon Transformation

There is a close relationship between Radon and the Fourier transformations!

F{p(s)}=P(ω)=12π∫∞−∞p(s)e−iωsds

with:

  • Summation of lines causes duplication in the center
  • Oversampling in the center of the Fourier space
26 / 45

PD Stefan Bosse - AFEML - Module Y: Radon Transformation

Radon Transformation

  • Via the central slice theorem we get relationship between projection and image slice space:

F1{p(s,ϕ´)}=F2{f(x,y)}∣ϕ=ϕ´

F1: Take a 2D function f(x,y), project it onto a line, and do a FT of that projection ⇔ F2: Do a 2D FT of f(x,y) first, and then take a slice through origin parallel to the projection line.

27 / 45

PD Stefan Bosse - AFEML - Module Y: Filtered backprojection

Filtered backprojection

  • We saw the relationship between projection space and spatial frequency space via the Fourier transformation:

f(x,y)=F−12{F(vx,vy)}

with vx=ω cos(φ), vy=ω sin(φ), and dvxdvy=ωdωdφ

  • For the reconstruction we get finally:

f(x,y)=∫π0dϕ[∫∞−∞dω|ω|P(ω)e2πiωs)]s=xcosϕ+ysinϕ

28 / 45

PD Stefan Bosse - AFEML - Module Y: Filtered backprojection

Filtered backprojection

That measn for the FBP "algorithm":

  • P: FT of projection p(s,φ)
  • Multiply by frequency filter |ω|
  • Inverse-transform this product
  • This filtered projection is backprojected
  • Then sum over all filtered projections

A complete set of 1D projections allows the reconstruction of the original 2D distribution without loss of information (theoretically=

29 / 45

PD Stefan Bosse - AFEML - Module Y: Filtered backprojection

Filtered backprojection

Filters

Baruchel, 2000

  1. Ramp filter
  2. Sine wave filter
  3. Ramp + Hanning function
  4. ...

? What is the right filter? Is there a right filter?

! All filter functions have advantages and disadvantages!

Transfer functions of different filters |ω| in the frequency space

30 / 45

PD Stefan Bosse - AFEML - Module Y: Filtered backprojection

Filtered backprojection

Filters

Ramp filter:

  • Blurring (low ω) is minimized
  • Contrast (high ω) is accentuated
  • But noise is amplified (high ω), too!

Hanning/Sine filters:

  • Blurring (low ω) is minimized (improved)
  • Contrast (mid ω) is accentuated
  • Noise (high ω) is reduced
  • Better compromise between blurring and noise!
31 / 45

PD Stefan Bosse - AFEML - Module Y: Filtered backprojection

Filtered backprojection

  • Analytical Reconstruction works exactly for lim N → ∞! It bases on Poisson distribution, i.e., Signal-to-noise ratio ≈ √N

Zvolský, 2014 The reconstruction quality increases with the number of projections

32 / 45

PD Stefan Bosse - AFEML - Module Y: Filtered backprojection

Filtered backprojection

Zvolský, 2014 BP (left) vs FBP (right)

33 / 45

PD Stefan Bosse - AFEML - Module Y: Reconstruction and Number of Projections

Reconstruction and Number of Projections

http://www.impactscan.org/slides/impactday/basicct The quality of the reconstruction and artifacts (noise) depends on the number of used projections

34 / 45

PD Stefan Bosse - AFEML - Module Y: Reconstruction Errors

Reconstruction Errors

Reality is more complex:

  • Data is discrete: FBP only precise if all angles are available
  • Data is noisy: Measurements follow a probability distribution
  • Detectors are unprecise: mis-positioning of photons
  • Detector geometry may not provide complete data
  • Not all photons travel along straight lines: scatter, absorption

Other reconstruction methods:

  • Itertative Image Reconstruction
  • Statistical methods (using forward- and backward projection learning)
  • ML-based (data-driven) reconstruction models
35 / 45

PD Stefan Bosse - AFEML - Module Y: Reconstruction Errors

Baruchel, 2000

Much attention must also be paid to the noise of the camera or, more precisely, to its dynamic range. When the inspection's issue is the determination of the accurate size of some internal feature, or the local characterization of materials (density measurement for instance), then an increased attention must be paid to the reconstruction artifacts. They create artificial patterns inside the reconstructed slice (streak artifacts), or they locally modify the pixels values (cupping effect), and hence the quantitative result.

36 / 45

PD Stefan Bosse - AFEML - Module Y: Reconstruction Errors

Reconstruction Errors

Artifacts and errors in slice reconstruction due to:

Baruchel, X-Ray Tomography in Material Science. 2000, pp 23

  • Beam hardening
  • Detector saturation
  • Aliasing, Line and ring artifacts
  • Scattered photons
  • Ill corrected detector
  • Spatial distortion of the detector
  • Centering error
37 / 45

PD Stefan Bosse - AFEML - Module Y: Aliasing

Aliasing

  • High (spatial) frequencies are encountered in the signal corresponding to every projection.

  • They are due to the steep edges which are eventually present in the object. As the detector samples the signal (all along the projection) with a non-zero step, high frequencies corrupt the data, within the Fourier domain. Streaks are generated.

  • Ill corrected detector: The signal delivered by every sensitive cell of the detector must be linearly spread between the offset level (corresponding to the absence of photons) and the gain level (corresponding to the non-attenuated flux).

38 / 45

PD Stefan Bosse - AFEML - Module Y: Aliasing

Aliasing

  • A bad correction of one cell will generate, in the reconstructed image a "ring artifact", i.e. the image of a ring, centered on the pixel corresponding to the location of the rotation axis, or "line artifacts" by a dead pixel.

  • Saturation noise: Different pixels have different saturation levels ⇒ spatial patterns!

39 / 45

PD Stefan Bosse - AFEML - Module Y: Aliasing

Aliasing

BDG – Richtlinie P 203 Line artifact due to a defect detector pixel

40 / 45

PD Stefan Bosse - AFEML - Module Y: Beam Hardening

Beam Hardening

X-ray emission is commonly polychromatic with a continuous energy spectrum. But μ and attenuation depends on energy! Beam hardening is the correction of this phenomena using physical energy or numerical filters.

BDG – Richtlinie P 203 (Left) Without beam hardening (boundary intensification) (Right) With beam hardening

41 / 45

PD Stefan Bosse - AFEML - Module Y: Rotation Axis Error

Rotation Axis Error

CT reconstruction with sine-wave filtered back projection from 800 projections (image width=400 pixels) of an Aluminum die casted plate with pores (40 mm width, 5 mm thickness), no X-ray energy filter (a) Centered rotation axis, error: ± 1 pixel accuracy (b) Rotation axis error 2 pixels! (c) Rotation axis error: 28 pixels

42 / 45

PD Stefan Bosse - AFEML - Module Y: Normal CT versa μCT

Normal CT versa μCT

  • "Normal" CT: FSD=0.5 mm, M=2, Detecor pixel size 200 μm, without X-ray energy filtration using classical sine-wave filtered backprojection algorithm;
  • μCT: FSD=50 μm, using Zeiss Xradia and proprietary reconstruction software and X-ray energy filtration

Specimen: Aluminum die casted plate with pores, 800-1000 projections (a) μCT reconstruction, 120 kV (b) Normal CT, 60 kV

43 / 45

PD Stefan Bosse - AFEML - Module Y: Do not trust CT

Do not trust CT

BDG – Richtlinie P 203 Comparison of a μCT slice (voxel size 20μm) and nearly the micrograph slice at nearly the same depth. There are differences and artifacts. Find them!

44 / 45

PD Stefan Bosse - AFEML - Module Y: Conclusions

Conclusions

  1. CT commonly records a set of radial 1-dim or 2-dim projections of a 3-dim object
  2. Backprojection is used to transform the radial 2-dim image projection set to a set of slice images forming a 3-dim volume
  3. We have: X-ray intensity images; We want: μ distributions!
  4. Backprojection introduces errors; blurring is compensated by frequency filtering.
  5. Do not trust reconstructed CT volume slices if accuracy is important! Too many artifacts ...
45 / 45