diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index 9c45dbe..d9f67a2 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.6","generation_timestamp":"2024-11-25T21:32:09","documenter_version":"1.2.1"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.7","generation_timestamp":"2025-01-15T17:25:32","documenter_version":"1.2.1"}} \ No newline at end of file diff --git a/dev/benchmark/index.html b/dev/benchmark/index.html index 71bb5e9..beda658 100644 --- a/dev/benchmark/index.html +++ b/dev/benchmark/index.html @@ -172,4 +172,4 @@ %time astra.algorithm.run(alg_id, 1) rec = astra.data3d.get(rec_id) -np.copy(rec)

+np.copy(rec)

diff --git a/dev/functions/index.html b/dev/functions/index.html index f7d37c8..3984cf2 100644 --- a/dev/functions/index.html +++ b/dev/functions/index.html @@ -42,7 +42,7 @@ 3.7320508075688767 0.0 3.7320508075688767 - 2.0source

Backproject (adjoint Radon)

RadonKA.backproject — Function
backproject(sinogram, θs; <kwargs>)

Conceptually the adjoint operation of radon. Intuitively, the backproject smears rays back into the space. See also radon.

For filtered backprojection see backproject_filtered.

Example

julia> arr = zeros((5,2)); arr[2,:] .= 1; arr[4, :] .= 1
+ 2.0
source

Backproject (adjoint Radon)

RadonKA.backproject — Function
backproject(sinogram, θs; <kwargs>)

Conceptually the adjoint operation of radon. Intuitively, the backproject smears rays back into the space. See also radon.

For filtered backprojection see backproject_filtered.

Example

julia> arr = zeros((5,2)); arr[2,:] .= 1; arr[4, :] .= 1
 2-element view(::Matrix{Float64}, 4, :) with eltype Float64:
  1.0
  1.0
@@ -69,5 +69,5 @@
  0.0  0.0  0.0       1.08455   0.0        0.0      0.0        1.08455   0.0       0.0
  0.0  0.0  0.591307  0.493247  0.0        0.0      0.0        0.493247  0.591307  0.0
  0.0  0.0  0.700352  0.0       0.0        0.0      0.0        0.0       0.700352  0.0
- 0.0  0.0  0.0       0.0       0.0        0.0      0.0        0.0       0.0       0.0
source
RadonKA.backproject_filtered — Function
backproject_filtered(sinogram, θs; 
-                        geometry, μ, filter)

Reconstruct the image from the sinogram using the filtered backprojection algorithm.

  • filter=nothing: The filter to be applied in Fourier space. If nothing, a ramp filter is used. filter should be a 1D array of the same length as the sinogram.

See radon for the explanation of the keyword arguments

source

Specifying Geometries

RadonKA.RadonGeometry — Type
abstract type RadonGeometry end

Abstract supertype for all geometries which are supported by radon and backproject.

List of geometries:

See radon and backproject how to apply.

source
RadonKA.RadonParallelCircle — Type
RadonParallelCircle(N, in_height, weights)

  • N is the size of the first and second dimension of the input array for radon.

  • in_height is a vector or a range indicating the positions of the detector with respect to the midpoint which is located at N ÷ 2 + 1. The rays travel along straight parallel paths through the array. Maximum and minimum values are (N+1) ÷ 2 - 1 and -(N+1) ÷ 2 + 1 respectively.

For example, for an array of size N=10 the default definition is: RadonParallelCircle(10, -4:4) So the resulting sinogram has the shape: (9, length(angles), size(array, 3)).

  • weights can weight each of the rays with different strength. Per default weights = 0 .* in_height .+ 1

See radon and backproject how to apply.

source
RadonKA.RadonFlexibleCircle — Type
RadonFlexibleCircle(N, in_height, out_height, weights)
  • N is the size of the first and second dimension of the input for radon.
  • in_height is a vector or range indicating the vertical positions of the rays entering the circle with respect to the midpoint which is located at N ÷ 2 + 1. Maximum and minimum values are (N+1) ÷ 2 - 1 and -(N+1) ÷ 2 + 1 respectively.
  • out_height is a vector or range indicating the vertical positions of the rays exiting the circle with respect to the midpoint which is located at N ÷ 2 + 1. Maximum and minimum values are (N+1) ÷ 2 - 1 and -(N+1) ÷ 2 + 1 respectively.

One definition could be: RadonFlexibleCircle(10, -4:4, zeros((9,))) It would describe rays which enter the circle at positions -4:4 but all of them would focus at the position 0 when leaving the circle. This is an extreme form of fan beam tomography.

  • weights can weight each of the rays with different strength. Per default weights = 0 .* in_height .+ 1

See radon and backproject how to apply.

source
+ 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0source
RadonKA.backproject_filtered — Function
backproject_filtered(sinogram, θs; 
+                        geometry, μ, filter)

Reconstruct the image from the sinogram using the filtered backprojection algorithm.

  • filter=nothing: The filter to be applied in Fourier space. If nothing, a ramp filter is used. filter should be a 1D array of the same length as the sinogram.

See radon for the explanation of the keyword arguments

source

Specifying Geometries

RadonKA.RadonGeometry — Type
abstract type RadonGeometry end

Abstract supertype for all geometries which are supported by radon and backproject.

List of geometries:

See radon and backproject how to apply.

source
RadonKA.RadonParallelCircle — Type
RadonParallelCircle(N, in_height, weights)

  • N is the size of the first and second dimension of the input array for radon.

  • in_height is a vector or a range indicating the positions of the detector with respect to the midpoint which is located at N ÷ 2 + 1. The rays travel along straight parallel paths through the array. Maximum and minimum values are (N+1) ÷ 2 - 1 and -(N+1) ÷ 2 + 1 respectively.

For example, for an array of size N=10 the default definition is: RadonParallelCircle(10, -4:4) So the resulting sinogram has the shape: (9, length(angles), size(array, 3)).

  • weights can weight each of the rays with different strength. Per default weights = 0 .* in_height .+ 1

See radon and backproject how to apply.

source
RadonKA.RadonFlexibleCircle — Type
RadonFlexibleCircle(N, in_height, out_height, weights)
  • N is the size of the first and second dimension of the input for radon.
  • in_height is a vector or range indicating the vertical positions of the rays entering the circle with respect to the midpoint which is located at N ÷ 2 + 1. Maximum and minimum values are (N+1) ÷ 2 - 1 and -(N+1) ÷ 2 + 1 respectively.
  • out_height is a vector or range indicating the vertical positions of the rays exiting the circle with respect to the midpoint which is located at N ÷ 2 + 1. Maximum and minimum values are (N+1) ÷ 2 - 1 and -(N+1) ÷ 2 + 1 respectively.

One definition could be: RadonFlexibleCircle(10, -4:4, zeros((9,))) It would describe rays which enter the circle at positions -4:4 but all of them would focus at the position 0 when leaving the circle. This is an extreme form of fan beam tomography.

  • weights can weight each of the rays with different strength. Per default weights = 0 .* in_height .+ 1

See radon and backproject how to apply.

source
diff --git a/dev/geometries/index.html b/dev/geometries/index.html index 6954b7f..547a8c0 100644 --- a/dev/geometries/index.html +++ b/dev/geometries/index.html @@ -36,4 +36,4 @@ simshow(μ)

projected_1 = backproject(sinogram, angles, μ=μ);
 projected_2 = backproject(sinogram, angles .+ π / 4, μ=μ);
 
-[simshow(projected_1) simshow(projected_2)]

+[simshow(projected_1) simshow(projected_2)]

diff --git a/dev/index.html b/dev/index.html index 3acefe1..413c69c 100644 --- a/dev/index.html +++ b/dev/index.html @@ -20,4 +20,4 @@ eprint={2402.06283}, archivePrefix={arXiv}, primaryClass={physics.optics} -} +} diff --git a/dev/mathematical/index.html b/dev/mathematical/index.html index e706422..228bb26 100644 --- a/dev/mathematical/index.html +++ b/dev/mathematical/index.html @@ -1,2 +1,2 @@ -Mathematical Background · RadonKA.jl
GitHub

Mathematial background

Introduction

The Radon transform [1] is an integral transform which is the foundation behind computed tomography (CT). It projects the values of a function along straight lines onto a detector. These projections are calculated for a set of different observation angles around the function. Such a dataset allows to reconstruct the original function distribution.

Definition

The general definition of the parallel attenuated Radon transform (ART) $\mathcal{R}_{\mu}$ is

\[g(s, \theta) = (\mathcal{R}_{\mu} f)(s, \theta) = \int_{-\infty}^{\infty} f(s\,\pmb{\omega} + t\,\pmb{\omega}^{\perp}) \exp\left(- \int_t^{\infty} \mu\left(s\,\pmb{\omega} + \tau\,\pmb{\omega}^{\perp}\right)\,\,\mathrm{d}\tau\right)\,\, \mathrm{d}t\]

where $\mu$ is the position dependent absorption coefficient. The vectors are defined as $\omega=(\cos \theta, \sin \theta)$ and $\omega^{\perp}=(-\sin \theta, \cos \theta)$.

In case of zero absorption ($\mu=0$) the ART reduces to the ordinary Radon transform.

Adjoint

If $\mu$ is known, the adjoint of the ART with respect to $f$ is given by

\[(\mathcal{R}_{\mu}^* g)(x, y) = \int_{0}^{2\pi} g( \pmb{\omega} \cdot [x,y]^{\intercal}, \theta) \exp\left(- \int_{[x,y]^{\intercal} \cdot \pmb{\omega}^{\perp}}^{\infty} \mu\left([x,y]^{\intercal} \cdot \pmb{\omega} + \tau\,\pmb{\omega}^{\perp}\right)\,\,\mathrm{d}\tau\right)\,\, \mathrm{d}\theta\]

where $\mu$ is the position dependent absorption coefficient. Qualitatively, the Radon transform calculates the shadow of a sample if illuminated with plane waves.

Discretization

On the other hand, the adjoint Radon transform smears the shadows (projection patterns) back into space and deposits values at each point. The pair of $\mathcal{R}_{\mu}$ and $\mathcal{R}_{\mu}^*$ can be used in gradient descent based optimizers.

In the figure below we show the geometry to the presented equations. Further, we show how the grid is discretized. RadonKA.jl traces a ray through the pixel grid. The function value it accumulates (Radon transform) or deposits (adjoint Radon transform) depends on the intersection length with each pixel.

The geometry of the Radon transform.

Applications

The ordinary Radon transform is the mathematical concept behind CT [2]. CT allows to recover a 3D absorption map of objects such as human bodies. One application of the ART is single photon emission computerized tomography (SPECT) [3]. In SPECT a tracer is injected in the specimen (e.g. blood of human bodies). The trace is radioactive and emits gamma rays which are detected. Since the gamma rays are created inside the body, its intensity is weakened from the starting point of the emission until hitting the detector. A more recent usage of the ART is within the context of Tomographic Volumetric Additive Manufacturing (TVAM) [4], [5].

In TVAM a photosensitive resin is illuminated with projection patterns from different angular positions. The light is absorbed while propagating through that medium and reacting with the photo-initiator. Once a certain voxel in space received a sufficient energy dose, it is going to polymerize. The mathematical process behind this printer is given by the adjoint of the ART.

References

[1]
J. Radon. Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten. Akad. Wiss. 69, 262–277 (1917).
[2]
T. M. Buzug. Computed Tomography. In: Springer Handbook of Medical Technology, edited by R. Kramme, K.-P. Hoffmann and R. S. Pozos (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011); pp. 311–342.
[3]
[4]
B. E. Kelly, I. Bhattacharya, H. Heidari, M. Shusteff, C. M. Spadaccini and H. K. Taylor. Volumetric additive manufacturing via tomographic reconstruction. Science 363, 1075–1079 (2019).
[5]
D. Loterie, P. Delrot and C. Moser. High-resolution tomographic volumetric additive manufacturing. Nature Communications 11, 852 (2020).
+Mathematical Background · RadonKA.jl
GitHub

Mathematial background

Introduction

The Radon transform [1] is an integral transform which is the foundation behind computed tomography (CT). It projects the values of a function along straight lines onto a detector. These projections are calculated for a set of different observation angles around the function. Such a dataset allows to reconstruct the original function distribution.

Definition

The general definition of the parallel attenuated Radon transform (ART) $\mathcal{R}_{\mu}$ is

\[g(s, \theta) = (\mathcal{R}_{\mu} f)(s, \theta) = \int_{-\infty}^{\infty} f(s\,\pmb{\omega} + t\,\pmb{\omega}^{\perp}) \exp\left(- \int_t^{\infty} \mu\left(s\,\pmb{\omega} + \tau\,\pmb{\omega}^{\perp}\right)\,\,\mathrm{d}\tau\right)\,\, \mathrm{d}t\]

where $\mu$ is the position dependent absorption coefficient. The vectors are defined as $\omega=(\cos \theta, \sin \theta)$ and $\omega^{\perp}=(-\sin \theta, \cos \theta)$.

In case of zero absorption ($\mu=0$) the ART reduces to the ordinary Radon transform.

Adjoint

If $\mu$ is known, the adjoint of the ART with respect to $f$ is given by

\[(\mathcal{R}_{\mu}^* g)(x, y) = \int_{0}^{2\pi} g( \pmb{\omega} \cdot [x,y]^{\intercal}, \theta) \exp\left(- \int_{[x,y]^{\intercal} \cdot \pmb{\omega}^{\perp}}^{\infty} \mu\left([x,y]^{\intercal} \cdot \pmb{\omega} + \tau\,\pmb{\omega}^{\perp}\right)\,\,\mathrm{d}\tau\right)\,\, \mathrm{d}\theta\]

where $\mu$ is the position dependent absorption coefficient. Qualitatively, the Radon transform calculates the shadow of a sample if illuminated with plane waves.

Discretization

On the other hand, the adjoint Radon transform smears the shadows (projection patterns) back into space and deposits values at each point. The pair of $\mathcal{R}_{\mu}$ and $\mathcal{R}_{\mu}^*$ can be used in gradient descent based optimizers.

In the figure below we show the geometry to the presented equations. Further, we show how the grid is discretized. RadonKA.jl traces a ray through the pixel grid. The function value it accumulates (Radon transform) or deposits (adjoint Radon transform) depends on the intersection length with each pixel.

The geometry of the Radon transform.

Applications

The ordinary Radon transform is the mathematical concept behind CT [2]. CT allows to recover a 3D absorption map of objects such as human bodies. One application of the ART is single photon emission computerized tomography (SPECT) [3]. In SPECT a tracer is injected in the specimen (e.g. blood of human bodies). The trace is radioactive and emits gamma rays which are detected. Since the gamma rays are created inside the body, its intensity is weakened from the starting point of the emission until hitting the detector. A more recent usage of the ART is within the context of Tomographic Volumetric Additive Manufacturing (TVAM) [4], [5].

In TVAM a photosensitive resin is illuminated with projection patterns from different angular positions. The light is absorbed while propagating through that medium and reacting with the photo-initiator. Once a certain voxel in space received a sufficient energy dose, it is going to polymerize. The mathematical process behind this printer is given by the adjoint of the ART.

References

[1]
J. Radon. Über die Bestimmung von Funktionen durch ihre Integralwerte längs gewisser Mannigfaltigkeiten. Akad. Wiss. 69, 262–277 (1917).
[2]
T. M. Buzug. Computed Tomography. In: Springer Handbook of Medical Technology, edited by R. Kramme, K.-P. Hoffmann and R. S. Pozos (Springer Berlin Heidelberg, Berlin, Heidelberg, 2011); pp. 311–342.
[3]
[4]
B. E. Kelly, I. Bhattacharya, H. Heidari, M. Shusteff, C. M. Spadaccini and H. K. Taylor. Volumetric additive manufacturing via tomographic reconstruction. Science 363, 1075–1079 (2019).
[5]
D. Loterie, P. Delrot and C. Moser. High-resolution tomographic volumetric additive manufacturing. Nature Communications 11, 852 (2020).
diff --git a/dev/tutorial/index.html b/dev/tutorial/index.html index 3aa39df..c905f6a 100644 --- a/dev/tutorial/index.html +++ b/dev/tutorial/index.html @@ -54,4 +54,4 @@ @time radon(volume, angles); # 2.527153 seconds (8.95 k CPU allocations: 363.703 KiB) (2 GPU allocations: 498.027 MiB, 0.06% memmgmt time) -CUDA.@time CUDA.@sync radon(volume_c, angles_c); +CUDA.@time CUDA.@sync radon(volume_c, angles_c);