Move tutorials and how-tos to literate folder

Author: nHackel

docs/src/literate/howots/custom.jl

@@ -0,0 +1,19 @@
+# # Block Diagonal Operator
+include("../../util.jl") #hide
+# This operator represents a block-diagonal matrix out of given operators.
+# One can also provide a single-operator and a number of blocks. In that case the given operator is repeated for each block.
+# In the case of stateful operators, one can supply a method for copying the operators.
+blocks = N
+ops = [WeightingOp(fill(i % 2, N)) for i = 1:N]
+dop = DiagOp(ops)
+typeof(dop)
+
+# We can retrieve the operators:
+dop.ops
+
+# And visualize the result:
+fig = Figure()
+plot_image(fig[1,1], image, title = "Image")
+plot_image(fig[1,2], reshape(dop * vec(image), N, N), title = "Block Weighted")
+resize_to_layout!(fig)
+fig

docs/src/literate/howots/gpu.jl

@@ -0,0 +1,19 @@
+# # Fourier Operator
+include("../../util.jl") #hide
+# The Fourier operator and its related operators for the discrete cosine and sine transform are available
+# whenever FFTW.jl is loaded together with LinearOperatorCollection:
+using LinearOperatorCollection, FFTW
+fop = FFTOp(Complex{eltype(image)}, shape = (N, N))
+cop = DCTOp(eltype(image), shape = (N, N))
+sop = DSTOp(eltype(image), shape = (N, N))
+image_frequencies = reshape(fop * vec(image), N, N)
+image_cosine = reshape(cop * vec(image), N, N)
+image_sine = reshape(sop * vec(image), N, N)
+
+fig = Figure()
+plot_image(fig[1,1], image, title = "Image")
+plot_image(fig[1,2], abs.(weighted_frequencies) .+ eps(), title = "Frequency Domain", colorscale = log10)
+plot_image(fig[1,3], image_cosine, title = "Cosine")
+plot_image(fig[1,4], image_sine, title = "Sine")
+resize_to_layout!(fig)
+fig

docs/src/literate/tutorials/diagonal.jl

@@ -0,0 +1,10 @@
+# # Gradient Operator
+include("../../util.jl") #hide
+# This operator computes a direction gradient along one or more dimensions of an array:
+gop = GradientOp(eltype(image); shape = (N, N), dims = 1)
+gradients = reshape(gop * vec(image), :, N)
+fig = Figure()
+plot_image(fig[1,1], image, title = "Image")
+plot_image(fig[1,2], gradients[:, :], title = "Gradient", colormap = :vik)
+resize_to_layout!(fig)
+fig

docs/src/literate/tutorials/fft.jl

@@ -0,0 +1 @@
+include("../../util.jl") #hide

docs/src/literate/tutorials/gradient.jl

@@ -0,0 +1,38 @@
+# # Normal operator
+include("../../util.jl") #hide
+
+# This operator implements a lazy normal operator implementing:
+# ```math
+# \begin{equation}
+#   (\mathbf{A})^*\mathbf{A}
+# \end{equation}
+# ```
+# for some operator $\mathbf{A}$:
+fop = op = FFTOp(ComplexF32, shape = (N, N))
+nop = NormalOp(eltype(fop), parent = fop)
+isapprox(nop * vec(image), vec(image))
+
+# And we can again access our original operator:
+typeof(nop.parent)
+
+# LinearOperatorCollection also provides an opinionated `normalOperator` function which tries to optimize the resulting normal operator.
+# As an example consider the normal operator of a weighted fourier operator:
+weights = collect(range(0, 1, length = N*N))
+wop = WeightingOp(weights)
+pop = ProdOp(wop, fop)
+nop = normalOperator(pop)
+typeof(nop.parent) == typeof(pop)
+# Note that the parent was changed. This is because the normal operator was optimized by initially computing the weights:
+# ```math
+# \begin{equation}
+#   \tilde{\mathbf{W}} = \mathbf{W}^*\mathbf{W}
+# \end{equation}
+# ```
+# and then applying the following each iteration:
+# ```math
+# \begin{equation}
+#   \mathbf{A}^*\tilde{\mathbf{W}}\mathbf{A}
+# \end{equation}
+# ```
+
+# Other operators can define different optimization strategies for the `normalOperator` method.

docs/src/literate/tutorials/nfft.jl

@@ -0,0 +1,46 @@
+# # Getting Started
+# To begin, we first need to load LinearOperatorCollection:
+using LinearOperatorCollection
+
+# If we require an operator which is implemented via a package extensions, we also need to include
+# the package that implements the functionality of the operator:
+using FFTW
+
+# As an introduction, we will construct a two dimensional FFT operator and apply it to an image.
+# To construct an operator we can either call its constructor directly:
+N = 256
+op = FFTOp(ComplexF32, shape = (N, N))
+typeof(op)
+
+# Or we can use the factory method:
+op = createLinearOperator(FFTOp{ComplexF64}, shape = (N, N))
+
+# We will use a Shepp-logan phantom as an example image:
+using ImagePhantoms, ImageGeoms
+image = shepp_logan(N, SheppLoganToft())
+size(image)
+
+# Since our operators are only defined for matrix-vector products, we can't directly apply them to the two-dimensional image.
+# We first have to reshape the image to a vector:
+y = op * vec(image);
+
+# Afterwards we can reshape the result and visualize it with CairoMakie:
+image_freq = reshape(y, N, N)
+using CairoMakie
+function plot_image(figPos, img; title = "", width = 150, height = 150, colorscale = identity)
+  ax = CairoMakie.Axis(figPos[1, 1]; yreversed=true, title, width, height)
+  hidedecorations!(ax)
+  hm = heatmap!(ax, img, colorscale = colorscale)
+  Colorbar(figPos[1, 2], hm)
+end
+fig = Figure()
+plot_image(fig[1,1], image, title = "Image")
+plot_image(fig[1,2], abs.(image_freq), title = "Frequency Domain", colorscale = log10)
+resize_to_layout!(fig)
+fig
+
+# To perform the inverse Fourier transform we can simply use the adjoint of our operator:
+image_inverse = reshape(adjoint(op) * y, N, N)
+plot_image(fig[1,3], real.(image_inverse), title = "Image after Inverse")
+resize_to_layout!(fig)
+fig

docs/src/literate/tutorials/normal.jl

@@ -0,0 +1,32 @@
+# # Product Operator
+include("../../util.jl") #hide
+# This operator describes the product or composition between two operators:
+weights = collect(range(0, 1, length = N*N))
+wop = WeightingOp(weights)
+fop = FFTOp(ComplexF64, shape = (N, N))
+# A feature of LinearOperators.jl is that operator can be cheaply transposed, conjugated and multiplied and only in the case of a matrix-vector product the combined operation is evaluated.
+tmp_op = wop * fop
+tmp_freqs = tmp_op * vec(image)
+
+
+# Similar to the WeightingOp, the main difference with the product operator provided by LinearOperatorCollection is the dedicated type, which allows for code specialisation.
+pop = ProdOp(wop, fop)
+# and the ability to retrieve the components:
+typoef(pop.A)
+# and
+typeof(pop.B)
+
+# Otherwise they compute the same thing:
+pop * vec(image) == tmp_op * vec(image)
+
+# Note that a product operator is not thread-safe.
+
+# We can again visualize our result:
+weighted_frequencies = reshape(pop * vec(image), N, N)
+image_inverse = reshape(adjoint(pop) * y, N, N)
+fig = Figure()
+plot_image(fig[1,1], image, title = "Image")
+plot_image(fig[1,2], abs.(weighted_frequencies) .+ eps(), title = "Frequency Domain", colorscale = log10)
+plot_image(fig[1,3], real.(image_inverse), title = "Image after Inverse")
+resize_to_layout!(fig)
+fig

docs/src/literate/tutorials/overview.jl

@@ -0,0 +1,13 @@
+# # Radon Operator
+include("../../util.jl") #hide
+# The Radon operator is available when loading RadonKA.jl and LinearOperatorCollection:
+using LinearOperatorCollection, RadonKA
+angles = collect(range(0, π, N))
+rop = RadonOp(eltype(image); angles, shape = size(image));
+sinogram = reshape(rop * vec(image), :, N)
+fig = Figure()
+plot_image(fig[1,1], image, title = "Image")
+plot_image(fig[1,2], sinogram, title = "Sinogram")
+plot_image(fig[1,3], reshape(adjoint(rop) * vec(sinogram), N, N), title = "Backprojection")
+resize_to_layout!(fig)
+fig