#---------------------------------------------------------
# # [Tomography](@id 01-tomography)
#---------------------------------------------------------

#=
This page gives an overview of the Julia package
[`Sinograms.jl`](https://github.com/JeffFessler/Sinograms.jl).

This page was generated from a single Julia file:
[01-tomography.jl](@__REPO_ROOT_URL__/01-tomography.jl).
=#

#md # In any such Julia documentation,
#md # you can access the source code
#md # using the "Edit on GitHub" link in the top right.

#md # The corresponding notebook can be viewed in
#md # [nbviewer](http://nbviewer.jupyter.org/) here:
#md # [`01-tomography.ipynb`](@__NBVIEWER_ROOT_URL__/01-tomography.ipynb),
#md # and opened in [binder](https://mybinder.org/) here:
#md # [`01-tomography.ipynb`](@__BINDER_ROOT_URL__/01-tomography.ipynb).


# ### Setup

# Packages needed here.

using Sinograms #: todo
using MIRTjim: jim, prompt
using InteractiveUtils: versioninfo


# The following line is helpful when running this example.jl file as a script;
# this way it will prompt user to hit a key after each figure is displayed.

isinteractive() ? jim(:prompt, true) : prompt(:draw);


#=
### Tomography

[Tomography](https://en.wikipedia.org/wiki/Tomography)
is the process of imaging cross sections of an object
without actually slicing the object.
There are many forms of tomography;
the description here
focuses on
[X-ray computed tomography (CT scans)](https://en.wikipedia.org/wiki/CT_scan).
(See
[SPECTrecon.jl](https://github.com/JeffFessler/SPECTrecon.jl)
for a Julia package related to
[SPECT](https://en.wikipedia.org/wiki/Single-photon_emission_computed_tomography)
imaging.)

In an
[X-ray CT imaging system](http://en.wikipedia.org/wiki/File:Ct-internals.jpg),
X-rays emitted from an X-ray source
are transmitted through an object
(e.g., a patient in medical CT)
towards a detector array.
The source and detector
[rotate rapidly]("https://www.youtube.com/watch?v=2CWpZKuy-NE)
around the object.
The signal intensities recorded by the detector
are related
to the sizes and densities
of the object materials
between the source and each detector element.

### Radon transform

The mathematical foundation
for 2D X-ray CT imaging
is the
[Radon transform](https://en.wikipedia.org/wiki/Radon_transform),
an
[Integral transform](https://en.wikipedia.org/wiki/Integral_transform)
that maps a 2D function
``f(x,y)``
into the collection of line integrals
through that function.
Here we describe each line
by its angle ``\phi``
measured counter-clockwise from the ``y`` axis,
and by the radial distance ``r`` of the line from the origin.
The collection of integrals
is called a sinogram.

Mathematically,
the Radon transform ``p(r,\phi)`` of ``f(x,y)`` is defined by
```math
p(r,\phi) = \int_{-\infty}^{\infty}
f(r \cos \phi - l \sin \phi,
r \sin \phi + l \cos \phi) \, \mathrm{d} l.
```

The Radon transform of the object that is 1 inside the unit disk
and 0 elsewhere is given by
```math
p_1(r,\phi) = 2 \sqrt{1 - r^2} \mathbb{1}_{\abs{r} < 1},
```
where ``\mathbb{1}`` denotes the
[indicator function](https://en.wikipedia.org/wiki/Indicator_function),
i.e.,
it is 1 for ``\abs{r} < 1`` and 0 elsewhere.

By the Radon transform's translation property,
the Radon transform
of a disk of radius ``r_0``
centered at ``(x_0,y_0)``
is given by
```math
p(r,\phi) = r_0 p_1
\left( \frac{r - (x_0 \cos \phi + y_0 \sin \phi)}{r_0}, ϕ \right)
```

### Sinogram example

Here is an display of that Radon transform
for a disk of radius ``3``
centered at coordinate ``(5,1)``.
Note that maximum value is approximately ``6``,
the length of the longest chord
through a disk of radius ``3``.
=#

nr = 128
dr = 20 / nr
r = ((-(nr-1)/2):((nr-1)/2)) * dr # radial sample locations
na = 130
ϕ = (0:(na-1))/na * π # angular samples
proj1 = r -> abs(r) < 1 ? 2 * sqrt(1 - r^2) : 0.
proj2 = (r, ϕ, x, y, r0) -> r0 * proj1(r/r0 - (x * cos(ϕ) + y * sin(ϕ))/r0)
sino = proj2.(r, ϕ', 5, 1, 3)
jimsino = (sino, title) -> jim(
    r, ϕ, sino; title, aspect_ratio=:none,
    xlabel = "r", ylabel = "ϕ", ylims=(0,π), yticks=([0, π], ["0", "π"]),
	yflip=false, xticks = [-10, 0, 2, 5, 8, 10],
	clim = (0, 6),
)
jimsino(sino, "Sinogram for one disk")


#=
As the preceding figure illustrates,
the Radon transform of a unit disk
has a somewhat sinusoidal shape.
Indeed every point in the ``(x,y)`` plane
traces out a distinct sinusoid
in the sinogram.

The mapping from the object ``f(x,y)``
to data like the sinogram ``p(r,\phi)``
is called the "forward problem."


### Image reconstruction

[Image reconstruction](http://en.wikipedia.org/wiki/Image_reconstruction)
is the process of solving the
[inverse problem](https://en.wikipedia.org/wiki/Inverse_problem)
of recovering the object ``f(x,y)``
from a measured sinogram,
i.e.,
(usually noisy) samples of ``p(r,\phi)``.


### FBP

A simple image reconstruction method
is called the
"filtered back-projection" (FBP) approach.
It works by filtering each row of the sinogram
with a filter,
called the
[ramp filter](https://en.wikipedia.org/wiki/Radon_transform#Radon_inversion_formula),
whose frequency response
is roughly ``\abs{\nu}``,
where ``\nu`` is the spatial frequency variable
(units cycles / m),
followed by a
[back-projection](https://en.wikipedia.org/wiki/Radon_transform#Dual_transform)
step that is the
[adjoint](https://en.wikipedia.org/wiki/Hermitian_adjoint)
of the Radon transform.


### FBP example

Here is an illustration
of using the `fbp` method
in this package
to perform image reconstruction
from the preceding sinogram.

=#

image = fbp(sino; dr)
(nx,ny) = size(image)
dx = dr # default
x = ((-(nx-1)/2):((nx-1)/2)) * dr
y = x
jim(x, y, image, "FBP image",
    xtick=[-10, 0, 2, 5, 8, 10],
    ytick=[-10, 0, -2, 1, 4, 10],
)


#=

The FBP reconstructed image
looks pretty similar
to a disk of radius ``3``
centered at ``(5,1)``
as expected.
However,
there are quite a few ripples;
these are
[aliasing](https://en.wikipedia.org/wiki/Aliasing)
artifacts
due to the finite sampling
in ``r`` and ``\phi``.


### Noise effects

Simulating the effects of measurement noise in sinogram
leads to even worse FBP results.

=#

noisy_sinogram = sino + 0.1 * randn(size(sino))
jimsino(noisy_sinogram, "Noisy sinogram")

#
noisy_fbp_image = fbp(noisy_sinogram; dr)
rmax = maximum(r)
fovmask = @. sqrt(abs2(x) + abs2(y)') ≤ rmax
jim(x, y, fovmask)
noisy_fbp_image .*= fovmask
jim(x, y, noisy_fbp_image, "Noisy FBP image"; clim=(0,1))
throw(0)


#=

Many computational imaging methods use system models
that are too large to store explicitly
as dense matrices,
but nevertheless
are represented mathematically
by a linear mapping `A`.

Often that linear map is thought of as a matrix,
but in imaging problems
it often is more convenient
to think of it as a more general linear operator.

The `LinearMapsAA` package
can represent both "matrix" versions
and "operator" versions
of linear mappings.
This page illustrates both versions
in the context of single-image
[super-resolution](https://en.wikipedia.org/wiki/Super-resolution_imaging)
imaging,
where the operator `A` maps a `M × N` image
into a coarser sampled image of size `M÷2 × N÷2`.

Here the operator `A` is akin to down-sampling,
except, rather than simple decimation,
each coarse-resolution pixel
is the average of a 2 × 2 block of pixels in the fine-resolution image.
=#

# ### System operator (linear mapping) for down-sampling

# Here is the "forward" function needed to model 2× down-sampling:

down1 = (x) -> (x[1:2:end,:] + x[2:2:end,:])/2 # 1D down-sampling by 2×
down2 = (x) -> down1(down1(x)')'; # 2D down-sampling by factor of 2×

# The `down2` function is a (bounded) linear operator
# and here is its adjoint:
down2_adj(y::AbstractMatrix{<:Number}) = kron(y, fill(0.25, (2,2)));


#=
Mathematically, and adjoint is a generalization
of the (Hermitian) transpose of a matrix.
For a (bounded) linear mapping `A` between
inner product space X
with inner product <.,.>_X
and inner product space Y
with inner product <.,.>_Y,
the adjoint of `A`, denoted `A'`,
is the unique bound linear mapping
that satisfies
<A x, y>_Y = <x, A' y>_X
for all x ∈ X and y ∈ Y.
One can verify that the `down2_adj` function
satisfies that equality
for the usual inner product
on the space of `M × N` images.
=#


# ### LinearMap as an operator for super-resolution

#=
We now pick a specific image size
and define the linear mapping
using the two functions above:
=#

nx, ny = 200, 256
A = LinearMapAA(down2, down2_adj, ((nx÷2)*(ny÷2), nx*ny);
    idim = (nx,ny), odim = (nx,ny) .÷ 2)

#=
The `idim` argument specifies
that the input image is of size `nx × ny`
and
the `odim` argument specifies
that the output image is of size `nx÷2 × ny÷2`.
This means that when we invoke
`y = A * x`
the input `x` must be a 2D array of size `nx × ny`
(not a 1D vector!)
and the output `y` will have size `nx÷2 × ny÷2`.
This behavior is a generalization
of what one might expect
from a conventional matrix-vector expression,
but is quite appropriate and convenient
for imaging problems.

Here is an illustration.
We start with a 2D test image.
=#

image = shepp_logan(ny, SheppLoganToft())[(ny-nx)÷2 .+ (1:nx),:]
jim(image, "SheppLogan")


# Apply the operator `A` to this image to down-sample it:

down = A * image
jim(down, title="Down-sampled image")


# Apply the adjoint of `A` to that result to "up-sample" it:
up = A' * down
jim(up, title="Adjoint: A' * y")


# That up-sampled image does not have the same range of values
# as the original image because `A'` is an adjoint, not an inverse!


# ### AbstractMatrix version

#=
Some users may prefer that the operator `A` behave more like a matrix.
We can implement approach from the same ingredients
by using `vec` and `reshape` judiciously.
The code is less elegant,
but similarly efficient
because `vec` and `reshape` are non-allocating operations.
=#

B = LinearMapAA(
        x -> vec(down2(reshape(x,nx,ny))),
        y -> vec(down2_adj(reshape(y,Int(nx/2),Int(ny/2)))),
        ((nx÷2)*(ny÷2), nx*ny),
    )

#=
To apply this version to our `image`
we must first vectorize it
because the expected input is a vector in this case.
And then we have to reshape the vector output
as a 2D array to look at it.
(This is why the operator version with `idim` and `odim` is preferable.)
=#

y = B * vec(image) # This would fail here without the `vec`!
jim(reshape(y, nx÷2, ny÷2)) # Annoying reshape needed here!


#=
Even though we write `y = A * x` above,
one must remember that internally `A` is not stored as a dense matrix.
It is simply a special variable type
that stores the forward function `down2` and the adjoint function `down2_adj`,
along with a few other values like `nx,ny`,
and applies those functions when needed.
Nevertheless,
we can examine elements of `A` and `B`
just like one would with any matrix,
at least for small enough examples to fit in memory.
=#


# Examine `A` and `A'`

nx, ny = 8,6
idim = (nx,ny)
odim = (nx,ny) .÷ 2
A = LinearMapAA(down2, down2_adj, ((nx÷2)*(ny÷2), nx*ny); idim, odim)

# Here is `A` shown as a Matrix:
jim(Matrix(A)', "A")

# Here is `A'` shown as a Matrix:
jim(Matrix(A')', "A'")


#=
When defining the adjoint function of a linear mapping,
it is very important to verify
that it is correct (truly the adjoint).

For a small problem we simply use the following test:
=#
@assert Matrix(A)' == Matrix(A')

# For some applications we must check approximate equality like this:
@assert Matrix(A)' ≈ Matrix(A')


# Here is a statistical test that is suitable for large operators.
# Often one would repeat this test several times:
T = eltype(A)
x = randn(T, idim)
y = randn(T, odim)

@assert sum((A*x) .* y) ≈ sum(x .* (A'*y)) # <Ax,y> = <x,A'y>


# ### Reproducibility

# This page was generated with the following version of Julia:

io = IOBuffer()
versioninfo(io)
split(String(take!(io)), '\n')


# And with the following package versions

import Pkg; Pkg.status()
#using ImagePhantoms: shepp_logan, SheppLoganToft
