Add initial documenter docs

.gitignore

@@ -1 +1,4 @@
-/Manifest.toml
+Manifest.toml
+docs/build/
+docs/site/
+docs/src/generated/

docs/Project.toml

@@ -0,0 +1,11 @@
+[deps]
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+ImageGeoms = "9ee76f2b-840d-4475-b6d6-e485c9297852"
+ImagePhantoms = "71a99df6-f52c-4da1-bd2a-69d6f37f3252"
+LinearOperatorCollection = "a4a2c56f-fead-462a-a3ab-85921a5f2575"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+NFFT = "efe261a4-0d2b-5849-be55-fc731d526b0d"
+RadonKA = "86de8297-835b-47df-b249-c04e8db91db5"
+Wavelets = "29a6e085-ba6d-5f35-a997-948ac2efa89a"

docs/make.jl

@@ -0,0 +1,60 @@
+using Documenter, Literate # Documentation
+using RadonKA, Wavelets, NFFT, FFTW # Extensions
+using CairoMakie, ImageGeoms, ImagePhantoms # Documentation Example Packages
+
+# Generate examples
+OUTPUT_BASE = joinpath(@__DIR__(), "src", "generated")
+INPUT_BASE = joinpath(@__DIR__(), "src", "literate")
+for (root, dirs, files) in walkdir(INPUT_BASE)
+    for dir in dirs
+        OUTPUT = joinpath(OUTPUT_BASE, dir)
+        INPUT = joinpath(INPUT_BASE, dir)
+        for file in filter(f -> endswith(f, ".jl"), readdir(INPUT))
+            Literate.markdown(joinpath(INPUT, file), OUTPUT)
+        end
+    end
+end
+
+modules = [LinearOperatorCollection,
+            isdefined(Base, :get_extension) ? Base.get_extension(LinearOperatorCollection, :LinearOperatorFFTWExt) : LinearOperatorCollection.LinearOperatorFFTWExt,
+            isdefined(Base, :get_extension) ? Base.get_extension(LinearOperatorCollection, :LinearOperatorNFFTExt) : LinearOperatorCollection.LinearOperatorNFFTExt,
+            isdefined(Base, :get_extension) ? Base.get_extension(LinearOperatorCollection, :LinearOperatorRadonKAExt) : LinearOperatorCollection.LinearOperatorRadonKAExt,
+            isdefined(Base, :get_extension) ? Base.get_extension(LinearOperatorCollection, :LinearOperatorWaveletExt) : LinearOperatorCollection.LinearOperatorWaveletExt]
+
+makedocs(
+    format = Documenter.HTML(prettyurls=get(ENV, "CI", "false") == "true",
+        canonical="https://github.com/JuliaImageRecon/LinearOperatorCollection.jl",
+        assets=String[],
+        collapselevel=1,
+    ),
+    modules = modules,
+    sitename = "LinearOperatorCollection",
+    authors = "Tobias Knopp, Niklas Hackelberg and Contributors",
+    pages = [
+        "Home" => "index.md",
+        "Getting Started" => "generated/tutorials/overview.md",
+        "Tutorials" => Any[
+            "Weighting Operator" => "generated/tutorials/weighting.md",
+            "FFT Operator" => "generated/tutorials/fft.md",
+            "Diagonal Operator" => "generated/tutorials/diagonal.md",
+            "Gradient Operator" => "generated/tutorials/gradient.md",
+            #"Sampling Operator" => "generated/tutorials/sampling.md",
+            #"NFFT Operator" => "generated/tutorials/nfft.md",
+            "Wavelet Operator" => "generated/tutorials/wavelet.md",
+            "Radon Operator" => "generated/tutorials/radon.md",
+            "Product Operator" => "generated/tutorials/product.md",
+            "Normal Operator" => "generated/tutorials/normal.md",
+        ],
+        "How to" => Any[
+            #"Implement Custom Operators" => "generated/howtos/custom.md",
+            "Enable GPU Acceleration" => "generated/howtos/gpu.md",
+        ],
+        #"Explanations" => Any[
+        #    "Operator Structure" => "operators.md",
+        #],
+        "Reference" => "references.md"
+    ],
+    warnonly = [:missing_docs]
+)
+
+deploydocs(repo   = "github.com/JuliaImageRecon/LinearOperatorCollection.jl")

docs/src/index.md

@@ -0,0 +1,54 @@
+# LinearOperatorCollection
+
+*Collection of linear operators for multi-dimensional signal and imaging tasks*
+
+
+## Introduction
+
+This package contains a collection of linear operators that are particularly useful for multi-dimensional signal and image processing tasks. Linear operators or linear maps behave like matrices in a matrix-vector product, but aren't necessarily matrices themselves. They can utilize more effective algorithms and can defer their computation until they are multiplied with a vector.
+
+All operators provided by this package extend types and methods [LinearOperators.jl](https://github.com/JuliaSmoothOptimizers/LinearOperators.jl). For example this package
+provides operators for the FFT (Fast Fourier Transform) and its non-equidistant variant (NFFT), the DCT (Discrete Cosine Transform), and the Wavelet transform. This package, however, does not implement
+these transformation itself but uses established libraries for them.
+
+LinearOperatorCollection's main purpose is provide a wrapper around low-level libraries like FFTW.jl and NFFT.jl, which allows using the transformations as linear operators, i.e., implementing `Op * x`, `adjoint(Op) * x` and the `mul!` based in-place variants of the former.
+
+## Installation
+
+Within Julia, use the package manager to install this package:
+```julia
+using Pkg
+Pkg.add("LinearOperatorCollection")
+```
+This will install `LinearOperatorCollection` and a subset of the available operators.
+To keep the load time of this package low, many operators are implemented using package extensions.
+For instance, in order to get the `FFTOp`, one needs to install not only `LinearOperatorCollection` but also `FFTW` and load both in a Julia sessiong:
+```julia
+Pkg.add("FFTW")
+using LinearOperatorCollection, FFTW
+```
+Small operators are implemented in LinearOperatorCollection directly.
+
+
+## License / Terms of Usage
+
+The source code of this project is licensed under the MIT license. This implies that
+you are free to use, share, and adapt it. However, please give appropriate credit
+by citing the project.
+
+## Contact
+
+If you have problems using the software, find mistakes, or have general questions please use
+the [issue tracker](hthttps://github.com/JuliaImageRecon/LinearOperatorCollection.jl/issues) to contact us.
+
+## Related Packages
+
+There exist many related packages which also implement efficient and/or lazy operators:
+
+* [LinearOperators.jl](https://github.com/JuliaSmoothOptimizers/LinearOperators.jl)
+* [LinearMaps.jl](https://github.com/JuliaLinearAlgebra/LinearMaps.jl)
+* [LazyArrays.jl](https://github.com/JuliaArrays/LazyArrays.jl)
+* [BlockArrays.jl](https://github.com/JuliaArrays/BlockArrays.jl)
+* [Kronecker.jl](https://github.com/MichielStock/Kronecker.jl)
+
+Generally, it should be possible to combine operators and arrays from various packages.

docs/src/literate/howtos/custom.jl

@@ -0,0 +1,2 @@
+# # Implement Custom Operators
+# There are two different ways one can implement a custom LinearOperator. The first one is to directly implement an operator as a LinearOperator from LinearOperators.jl:

docs/src/literate/howtos/gpu.jl

@@ -0,0 +1,19 @@
+# # GPU Acceleration
+include("../../util.jl") #hide
+# GPU kernels generally require all their arguments to exist on the GPU. This is not ncessarily the case for matrix-free operators as provides LinearOperators or LinearOperatorCollection.
+# In the case that a matrix free operator is solely a function call and contains no internal array state, the operator is GPU compatible as long as the method has a GPU compatible implementation.
+
+# If the operator has internal fields required for its computation, such as temporary arrays for intermediate values or indices, then it needs to move those to the GPU.
+# Furthermore if the operator needs to create a new array in its execution, e.g. it is used in a non-inplace matrix-vector multiplication or it is combined with other operators, then the operator needs to specify
+# a storage type. LinearOperatorCollection has several GPU compatible operators, where the storage type is given by setting a `S` parameter:
+# ```julia
+# using CUDA # or AMDGPU, Metal, ...
+# image_gpu = cu(image)
+# ```
+using LinearOperatorCollection.LinearOperators
+image_gpu = image #hide
+storage = Complex.(similar(image_gpu, 0))
+fop = FFTOp(eltype(image_gpu), shape = (N, N), S = typeof(storage))
+LinearOperators.storage_type(fop) == typeof(storage)
+
+# GPU operators can be used just like the other operators. Note however, that a GPU operator does not necessarily work with a CPU vector.

docs/src/literate/tutorials/diagonal.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:
+typeof(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/tutorials/fft.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.(image_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/gradient.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/nfft.jl

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

docs/src/literate/tutorials/normal.jl

@@ -0,0 +1,39 @@
+# # 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}$:
+using FFTW
+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 = Float32.(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/overview.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/product.jl

@@ -0,0 +1,33 @@
+# # 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)
+typeof(pop)
+# and the ability to retrieve the components:
+typeof(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) * vec(weighted_frequencies), 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/radon.jl

@@ -0,0 +1,13 @@
+# # Radon Operator
+include("../../util.jl") #hide
+# The Radon operator is available when loading RadonKA.jl and LinearOperatorCollection:
+using 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

docs/src/literate/tutorials/sampling.jl

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

docs/src/literate/tutorials/wavelet.jl

@@ -0,0 +1,11 @@
+# # Wavelet Operator
+include("../../util.jl") #hide
+# The wavelet operator is available when loading Wavelets.jl together with LinearOperatorCollection:
+using Wavelets
+wop = WaveletOp(eltype(image), shape = (N, N))
+wavelet_image = reshape(wop * vec(image), N, N)
+fig = Figure()
+plot_image(fig[1,1], image, title = "Image")
+plot_image(fig[1,2], abs.(wavelet_image) .+ eps(), title = "Wavelet Image", colorscale = sqrt)
+resize_to_layout!(fig)
+fig

docs/src/literate/tutorials/weighting.jl

@@ -0,0 +1,20 @@
+# # Weighting Operator
+include("../../util.jl") #hide
+# The weighting operator implements a diagonal matrix which multiplies a vector index-wise with given weights.
+# Such an operator is also implemented within LinearOperator.jl, however here this operator has a dedicated type on which one can dispatch:
+weights = collect(range(0, 1, length = N*N))
+op = WeightingOp(weights)
+typeof(op)
+
+# One can retrieve the weights from the operator
+op.weights == weights
+
+# Note that we didn't need to specify the element type of the operator here. In this case the eltype was derived from the provided weights.
+weighted_image = reshape(op * vec(image), N, N);
+
+# To visualize our weighted image, we will again use CairoMakie:
+fig = Figure()
+plot_image(fig[1,1], image, title = "Image")
+plot_image(fig[1,2], weighted_image, title = "Weighted Image")
+resize_to_layout!(fig)
+fig