From ab0797bfb6cdc481b989858cd98589ea204dac51 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Mon, 11 Aug 2025 16:08:04 +0200
Subject: [PATCH 01/11] Prepare documentation structure

---
 .gitignore                     |  5 +++-
 docs/Project.toml              | 11 +++++++
 docs/make.jl                   | 54 ++++++++++++++++++++++++++++++++++
 docs/src/index.md              | 54 ++++++++++++++++++++++++++++++++++
 docs/src/tutorials/overview.jl | 46 +++++++++++++++++++++++++++++
 5 files changed, 169 insertions(+), 1 deletion(-)
 create mode 100644 docs/Project.toml
 create mode 100644 docs/make.jl
 create mode 100644 docs/src/index.md
 create mode 100644 docs/src/tutorials/overview.jl

diff --git a/.gitignore b/.gitignore
index b067edd..e5de18e 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,4 @@
-/Manifest.toml
+Manifest.toml
+docs/build/
+docs/site/
+docs/src/generated/
diff --git a/docs/Project.toml b/docs/Project.toml
new file mode 100644
index 0000000..5fe21cd
--- /dev/null
+++ b/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"
diff --git a/docs/make.jl b/docs/make.jl
new file mode 100644
index 0000000..78ffeaa
--- /dev/null
+++ b/docs/make.jl
@@ -0,0 +1,54 @@
+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
+
+makedocs(
+    format = Documenter.HTML(prettyurls=get(ENV, "CI", "false") == "true",
+        canonical="https://github.com/JuliaImageRecon/LinearOperatorCollection.jl",
+        assets=String[],
+        collapselevel=1,
+    ),
+    modules = [LinearOperatorCollection],
+    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",
+            "Product Operator" => "generated/tutorials/product.md",
+            "Diagonal Operator" => "generated/tutorials/diagonal.md",
+            "Normal Operator" => "generated/tutorials/normal.md",
+            "Gradient Operator" => "generated/tutorials/gradient.md",
+            "FFT Operator" => "generated/tutorials/fft.md",
+            "NFFT Operator" => "generated/tutorials/nfft.md",
+            "Wavelet Operator" => "generated/tutorials/wavelet.md",
+            "Radon Operator" => "generated/tutorials/radon.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/MagneticParticleImaging/MPIReco.jl.git",
+           target = "build")
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 0000000..2850d67
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,54 @@
+# LinearOperatorCollection
+
+*Collection of linear operators for multi-dimensional signal and imaging tasks*
+
+
+## Purpose
+
+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.
diff --git a/docs/src/tutorials/overview.jl b/docs/src/tutorials/overview.jl
new file mode 100644
index 0000000..70ac349
--- /dev/null
+++ b/docs/src/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(ComplexF64, 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

From d4910dfa72b7dda5e599e7698ed9496a019ade01 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Mon, 11 Aug 2025 16:19:14 +0200
Subject: [PATCH 02/11] Add weighting operator docs

---
 docs/src/tutorials/weighting.jl | 22 ++++++++++++++++++++++
 1 file changed, 22 insertions(+)
 create mode 100644 docs/src/tutorials/weighting.jl

diff --git a/docs/src/tutorials/weighting.jl b/docs/src/tutorials/weighting.jl
new file mode 100644
index 0000000..56c3174
--- /dev/null
+++ b/docs/src/tutorials/weighting.jl
@@ -0,0 +1,22 @@
+# # 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)
+
+# Such an operator is also implemented within LinearOperator.jl, however here this operator has a dedicated type on which one can dispatch:
+typeof(op)
+
+# And it is possible to 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
\ No newline at end of file

From 85d9491904916da4395117931cff76ea6537116d Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Mon, 11 Aug 2025 16:37:24 +0200
Subject: [PATCH 03/11] Add prodop documentation

---
 docs/src/tutorials/product.jl | 32 ++++++++++++++++++++++++++++++++
 1 file changed, 32 insertions(+)
 create mode 100644 docs/src/tutorials/product.jl

diff --git a/docs/src/tutorials/product.jl b/docs/src/tutorials/product.jl
new file mode 100644
index 0000000..6bd5f1d
--- /dev/null
+++ b/docs/src/tutorials/product.jl
@@ -0,0 +1,32 @@
+# # Product Operator
+
+# 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
\ No newline at end of file

From 6d843ade78d05e62c279a6b6d3cac470051748e1 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Mon, 11 Aug 2025 17:05:41 +0200
Subject: [PATCH 04/11] Add block diagonal documentation

---
 docs/src/tutorials/diagonal.jl | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)
 create mode 100644 docs/src/tutorials/diagonal.jl

diff --git a/docs/src/tutorials/diagonal.jl b/docs/src/tutorials/diagonal.jl
new file mode 100644
index 0000000..2fe99f8
--- /dev/null
+++ b/docs/src/tutorials/diagonal.jl
@@ -0,0 +1,18 @@
+# # Block Diagonal Operator
+# 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
\ No newline at end of file

From eee43d5e1b032495a0673c9de849a8c401aa81b0 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Mon, 11 Aug 2025 17:20:42 +0200
Subject: [PATCH 05/11] Add normal operator docu

---
 docs/src/tutorials/normal.jl | 38 ++++++++++++++++++++++++++++++++++++
 1 file changed, 38 insertions(+)
 create mode 100644 docs/src/tutorials/normal.jl

diff --git a/docs/src/tutorials/normal.jl b/docs/src/tutorials/normal.jl
new file mode 100644
index 0000000..d4899bc
--- /dev/null
+++ b/docs/src/tutorials/normal.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.
\ No newline at end of file

From e551ea066273975e550a1de551abe5de44d3f36d Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Mon, 11 Aug 2025 17:24:55 +0200
Subject: [PATCH 06/11] Fix deploy url

---
 docs/make.jl | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/docs/make.jl b/docs/make.jl
index 78ffeaa..a7c9c31 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -31,8 +31,8 @@ makedocs(
             "Weighting Operator" => "generated/tutorials/weighting.md",
             "Product Operator" => "generated/tutorials/product.md",
             "Diagonal Operator" => "generated/tutorials/diagonal.md",
-            "Normal Operator" => "generated/tutorials/normal.md",
             "Gradient Operator" => "generated/tutorials/gradient.md",
+            "Normal Operator" => "generated/tutorials/normal.md",
             "FFT Operator" => "generated/tutorials/fft.md",
             "NFFT Operator" => "generated/tutorials/nfft.md",
             "Wavelet Operator" => "generated/tutorials/wavelet.md",
@@ -50,5 +50,4 @@ makedocs(
     warnonly = [:missing_docs]
 )
 
-deploydocs(repo   = "github.com/MagneticParticleImaging/MPIReco.jl.git",
-           target = "build")
+deploydocs(repo   = "github.com/JuliaImageRecon/LinearOperatorCollection.jl")

From 5bfff9fbf69a639c17043fb5de38a22fc355703c Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Tue, 12 Aug 2025 14:03:50 +0200
Subject: [PATCH 07/11] Add gradient operator docs

---
 docs/src/tutorials/gradient.jl | 10 ++++++++++
 1 file changed, 10 insertions(+)
 create mode 100644 docs/src/tutorials/gradient.jl

diff --git a/docs/src/tutorials/gradient.jl b/docs/src/tutorials/gradient.jl
new file mode 100644
index 0000000..caa1e03
--- /dev/null
+++ b/docs/src/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
\ No newline at end of file

From fb5f87242e8a8d8d281a785261b57316a93e5ff4 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Tue, 12 Aug 2025 14:51:22 +0200
Subject: [PATCH 08/11] Add fft op

---
 docs/src/tutorials/fft.jl | 18 ++++++++++++++++++
 1 file changed, 18 insertions(+)
 create mode 100644 docs/src/tutorials/fft.jl

diff --git a/docs/src/tutorials/fft.jl b/docs/src/tutorials/fft.jl
new file mode 100644
index 0000000..3f00def
--- /dev/null
+++ b/docs/src/tutorials/fft.jl
@@ -0,0 +1,18 @@
+# # 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:
+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
\ No newline at end of file

From 7fec43c3841f0c964cfaf4f58f9a5a94971f8943 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Tue, 12 Aug 2025 14:57:38 +0200
Subject: [PATCH 09/11] Fix tmp array type for DCT and DST operator

---
 ext/LinearOperatorFFTWExt/DCTOp.jl | 2 +-
 ext/LinearOperatorFFTWExt/DSTOp.jl | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/ext/LinearOperatorFFTWExt/DCTOp.jl b/ext/LinearOperatorFFTWExt/DCTOp.jl
index 1f7743b..efe1ca7 100644
--- a/ext/LinearOperatorFFTWExt/DCTOp.jl
+++ b/ext/LinearOperatorFFTWExt/DCTOp.jl
@@ -34,7 +34,7 @@ returns a `DCTOpImpl <: AbstractLinearOperator` which performs a DCT on a given
 """
 function LinearOperatorCollection.DCTOp(T::Type; shape::Tuple, dcttype=2)
 
-  tmp=Array{Complex{real(T)}}(undef, shape) 
+  tmp=Array{T}(undef, shape) 
   if dcttype == 2
     plan = plan_dct!(tmp)
     iplan = plan_idct!(tmp)
diff --git a/ext/LinearOperatorFFTWExt/DSTOp.jl b/ext/LinearOperatorFFTWExt/DSTOp.jl
index 4cd4862..32e2ba9 100644
--- a/ext/LinearOperatorFFTWExt/DSTOp.jl
+++ b/ext/LinearOperatorFFTWExt/DSTOp.jl
@@ -32,7 +32,7 @@ returns a `LinearOperator` which performs a DST on a given input array.
 * `shape::Tuple`  - size of the array to transform
 """
 function LinearOperatorCollection.DSTOp(T::Type; shape::Tuple)
-  tmp=Array{Complex{real(T)}}(undef, shape)
+  tmp=Array{T}(undef, shape)
 
   plan = FFTW.plan_r2r!(tmp,FFTW.RODFT10)
   iplan = FFTW.plan_r2r!(tmp,FFTW.RODFT01)

From b6214f53690bca84b4cc9271eb51910eab77bad2 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Tue, 12 Aug 2025 15:28:41 +0200
Subject: [PATCH 10/11] Move tutorials and how-tos to literate folder

---
 docs/src/literate/howots/custom.jl       |  0
 docs/src/literate/howots/gpu.jl          |  0
 docs/src/literate/tutorials/diagonal.jl  | 19 ++++++++++
 docs/src/literate/tutorials/fft.jl       | 19 ++++++++++
 docs/src/literate/tutorials/gradient.jl  | 10 ++++++
 docs/src/literate/tutorials/nfft.jl      |  1 +
 docs/src/literate/tutorials/normal.jl    | 38 ++++++++++++++++++++
 docs/src/literate/tutorials/overview.jl  | 46 ++++++++++++++++++++++++
 docs/src/literate/tutorials/product.jl   | 32 +++++++++++++++++
 docs/src/literate/tutorials/radon.jl     | 13 +++++++
 docs/src/literate/tutorials/sampling.jl  |  2 ++
 docs/src/literate/tutorials/wavelet.jl   | 11 ++++++
 docs/src/literate/tutorials/weighting.jl | 22 ++++++++++++
 13 files changed, 213 insertions(+)
 create mode 100644 docs/src/literate/howots/custom.jl
 create mode 100644 docs/src/literate/howots/gpu.jl
 create mode 100644 docs/src/literate/tutorials/diagonal.jl
 create mode 100644 docs/src/literate/tutorials/fft.jl
 create mode 100644 docs/src/literate/tutorials/gradient.jl
 create mode 100644 docs/src/literate/tutorials/nfft.jl
 create mode 100644 docs/src/literate/tutorials/normal.jl
 create mode 100644 docs/src/literate/tutorials/overview.jl
 create mode 100644 docs/src/literate/tutorials/product.jl
 create mode 100644 docs/src/literate/tutorials/radon.jl
 create mode 100644 docs/src/literate/tutorials/sampling.jl
 create mode 100644 docs/src/literate/tutorials/wavelet.jl
 create mode 100644 docs/src/literate/tutorials/weighting.jl

diff --git a/docs/src/literate/howots/custom.jl b/docs/src/literate/howots/custom.jl
new file mode 100644
index 0000000..e69de29
diff --git a/docs/src/literate/howots/gpu.jl b/docs/src/literate/howots/gpu.jl
new file mode 100644
index 0000000..e69de29
diff --git a/docs/src/literate/tutorials/diagonal.jl b/docs/src/literate/tutorials/diagonal.jl
new file mode 100644
index 0000000..9438aa7
--- /dev/null
+++ b/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:
+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
\ No newline at end of file
diff --git a/docs/src/literate/tutorials/fft.jl b/docs/src/literate/tutorials/fft.jl
new file mode 100644
index 0000000..2da012e
--- /dev/null
+++ b/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.(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
\ No newline at end of file
diff --git a/docs/src/literate/tutorials/gradient.jl b/docs/src/literate/tutorials/gradient.jl
new file mode 100644
index 0000000..caa1e03
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/docs/src/literate/tutorials/nfft.jl b/docs/src/literate/tutorials/nfft.jl
new file mode 100644
index 0000000..b89a53e
--- /dev/null
+++ b/docs/src/literate/tutorials/nfft.jl
@@ -0,0 +1 @@
+include("../../util.jl") #hide
diff --git a/docs/src/literate/tutorials/normal.jl b/docs/src/literate/tutorials/normal.jl
new file mode 100644
index 0000000..d4899bc
--- /dev/null
+++ b/docs/src/literate/tutorials/normal.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.
\ No newline at end of file
diff --git a/docs/src/literate/tutorials/overview.jl b/docs/src/literate/tutorials/overview.jl
new file mode 100644
index 0000000..c9d3538
--- /dev/null
+++ b/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
diff --git a/docs/src/literate/tutorials/product.jl b/docs/src/literate/tutorials/product.jl
new file mode 100644
index 0000000..45983d7
--- /dev/null
+++ b/docs/src/literate/tutorials/product.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
\ No newline at end of file
diff --git a/docs/src/literate/tutorials/radon.jl b/docs/src/literate/tutorials/radon.jl
new file mode 100644
index 0000000..e406318
--- /dev/null
+++ b/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 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
diff --git a/docs/src/literate/tutorials/sampling.jl b/docs/src/literate/tutorials/sampling.jl
new file mode 100644
index 0000000..62b641d
--- /dev/null
+++ b/docs/src/literate/tutorials/sampling.jl
@@ -0,0 +1,2 @@
+# # Sampling
+include("../../util.jl") #hide
diff --git a/docs/src/literate/tutorials/wavelet.jl b/docs/src/literate/tutorials/wavelet.jl
new file mode 100644
index 0000000..59e6943
--- /dev/null
+++ b/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 LinearOperatorCollection, 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
diff --git a/docs/src/literate/tutorials/weighting.jl b/docs/src/literate/tutorials/weighting.jl
new file mode 100644
index 0000000..56c3174
--- /dev/null
+++ b/docs/src/literate/tutorials/weighting.jl
@@ -0,0 +1,22 @@
+# # 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)
+
+# Such an operator is also implemented within LinearOperator.jl, however here this operator has a dedicated type on which one can dispatch:
+typeof(op)
+
+# And it is possible to 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
\ No newline at end of file

From 1bd34e3ad101e26030215805f86f080e06a43819 Mon Sep 17 00:00:00 2001
From: nHackel <nHackel@users.noreply.github.com>
Date: Tue, 12 Aug 2025 18:39:54 +0200
Subject: [PATCH 11/11] Finish initial docs

---
 docs/make.jl                             | 27 ++++++++------
 docs/src/index.md                        |  2 +-
 docs/src/literate/howots/custom.jl       |  0
 docs/src/literate/howots/gpu.jl          |  0
 docs/src/literate/howtos/custom.jl       |  2 ++
 docs/src/literate/howtos/gpu.jl          | 19 ++++++++++
 docs/src/literate/tutorials/diagonal.jl  |  2 +-
 docs/src/literate/tutorials/fft.jl       |  2 +-
 docs/src/literate/tutorials/normal.jl    |  3 +-
 docs/src/literate/tutorials/product.jl   |  7 ++--
 docs/src/literate/tutorials/radon.jl     |  2 +-
 docs/src/literate/tutorials/wavelet.jl   |  2 +-
 docs/src/literate/tutorials/weighting.jl |  6 ++--
 docs/src/references.md                   | 24 +++++++++++++
 docs/src/tutorials/diagonal.jl           | 18 ----------
 docs/src/tutorials/fft.jl                | 18 ----------
 docs/src/tutorials/gradient.jl           | 10 ------
 docs/src/tutorials/normal.jl             | 38 --------------------
 docs/src/tutorials/overview.jl           | 46 ------------------------
 docs/src/tutorials/product.jl            | 32 -----------------
 docs/src/tutorials/weighting.jl          | 22 ------------
 docs/src/util.jl                         | 11 ++++++
 ext/LinearOperatorNFFTExt/NFFTOp.jl      | 23 ++++++------
 ext/LinearOperatorRadonKAExt/RadonOp.jl  | 13 +++++++
 src/ProdOp.jl                            | 10 +-----
 src/SamplingOp.jl                        | 20 +++++------
 26 files changed, 121 insertions(+), 238 deletions(-)
 delete mode 100644 docs/src/literate/howots/custom.jl
 delete mode 100644 docs/src/literate/howots/gpu.jl
 create mode 100644 docs/src/literate/howtos/custom.jl
 create mode 100644 docs/src/literate/howtos/gpu.jl
 create mode 100644 docs/src/references.md
 delete mode 100644 docs/src/tutorials/diagonal.jl
 delete mode 100644 docs/src/tutorials/fft.jl
 delete mode 100644 docs/src/tutorials/gradient.jl
 delete mode 100644 docs/src/tutorials/normal.jl
 delete mode 100644 docs/src/tutorials/overview.jl
 delete mode 100644 docs/src/tutorials/product.jl
 delete mode 100644 docs/src/tutorials/weighting.jl
 create mode 100644 docs/src/util.jl

diff --git a/docs/make.jl b/docs/make.jl
index a7c9c31..689f9e5 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -15,13 +15,19 @@ for (root, dirs, files) in walkdir(INPUT_BASE)
     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 = [LinearOperatorCollection],
+    modules = modules,
     sitename = "LinearOperatorCollection",
     authors = "Tobias Knopp, Niklas Hackelberg and Contributors",
     pages = [
@@ -29,22 +35,23 @@ makedocs(
         "Getting Started" => "generated/tutorials/overview.md",
         "Tutorials" => Any[
             "Weighting Operator" => "generated/tutorials/weighting.md",
-            "Product Operator" => "generated/tutorials/product.md",
+            "FFT Operator" => "generated/tutorials/fft.md",
             "Diagonal Operator" => "generated/tutorials/diagonal.md",
             "Gradient Operator" => "generated/tutorials/gradient.md",
-            "Normal Operator" => "generated/tutorials/normal.md",
-            "FFT Operator" => "generated/tutorials/fft.md",
-            "NFFT Operator" => "generated/tutorials/nfft.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"
+            "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",
+            #"Implement Custom Operators" => "generated/howtos/custom.md",
             "Enable GPU Acceleration" => "generated/howtos/gpu.md",
         ],
-        "Explanations" => Any[
-            "Operator Structure" => "operators.md",
-        ],
+        #"Explanations" => Any[
+        #    "Operator Structure" => "operators.md",
+        #],
         "Reference" => "references.md"
     ],
     warnonly = [:missing_docs]
diff --git a/docs/src/index.md b/docs/src/index.md
index 2850d67..bb9bc72 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -3,7 +3,7 @@
 *Collection of linear operators for multi-dimensional signal and imaging tasks*
 
 
-## Purpose
+## 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.
 
diff --git a/docs/src/literate/howots/custom.jl b/docs/src/literate/howots/custom.jl
deleted file mode 100644
index e69de29..0000000
diff --git a/docs/src/literate/howots/gpu.jl b/docs/src/literate/howots/gpu.jl
deleted file mode 100644
index e69de29..0000000
diff --git a/docs/src/literate/howtos/custom.jl b/docs/src/literate/howtos/custom.jl
new file mode 100644
index 0000000..0931837
--- /dev/null
+++ b/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:
diff --git a/docs/src/literate/howtos/gpu.jl b/docs/src/literate/howtos/gpu.jl
new file mode 100644
index 0000000..d4f8985
--- /dev/null
+++ b/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.
\ No newline at end of file
diff --git a/docs/src/literate/tutorials/diagonal.jl b/docs/src/literate/tutorials/diagonal.jl
index 9438aa7..43c3b05 100644
--- a/docs/src/literate/tutorials/diagonal.jl
+++ b/docs/src/literate/tutorials/diagonal.jl
@@ -9,7 +9,7 @@ dop = DiagOp(ops)
 typeof(dop)
 
 # We can retrieve the operators:
-dop.ops
+typeof(dop.ops)
 
 # And visualize the result:
 fig = Figure()
diff --git a/docs/src/literate/tutorials/fft.jl b/docs/src/literate/tutorials/fft.jl
index 2da012e..40acba6 100644
--- a/docs/src/literate/tutorials/fft.jl
+++ b/docs/src/literate/tutorials/fft.jl
@@ -12,7 +12,7 @@ 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,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)
diff --git a/docs/src/literate/tutorials/normal.jl b/docs/src/literate/tutorials/normal.jl
index d4899bc..38ecfed 100644
--- a/docs/src/literate/tutorials/normal.jl
+++ b/docs/src/literate/tutorials/normal.jl
@@ -8,6 +8,7 @@ include("../../util.jl") #hide
 # \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))
@@ -17,7 +18,7 @@ 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))
+weights = Float32.(collect(range(0, 1, length = N*N)))
 wop = WeightingOp(weights)
 pop = ProdOp(wop, fop)
 nop = normalOperator(pop)
diff --git a/docs/src/literate/tutorials/product.jl b/docs/src/literate/tutorials/product.jl
index 45983d7..c1453ca 100644
--- a/docs/src/literate/tutorials/product.jl
+++ b/docs/src/literate/tutorials/product.jl
@@ -3,7 +3,7 @@ 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))
+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)
@@ -11,8 +11,9 @@ 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:
-typoef(pop.A)
+typeof(pop.A)
 # and
 typeof(pop.B)
 
@@ -23,7 +24,7 @@ pop * vec(image) == tmp_op * vec(image)
 
 # We can again visualize our result:
 weighted_frequencies = reshape(pop * vec(image), N, N)
-image_inverse = reshape(adjoint(pop) * y, 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)
diff --git a/docs/src/literate/tutorials/radon.jl b/docs/src/literate/tutorials/radon.jl
index e406318..0f7062c 100644
--- a/docs/src/literate/tutorials/radon.jl
+++ b/docs/src/literate/tutorials/radon.jl
@@ -1,7 +1,7 @@
 # # Radon Operator
 include("../../util.jl") #hide
 # The Radon operator is available when loading RadonKA.jl and LinearOperatorCollection:
-using LinearOperatorCollection, RadonKA
+using RadonKA
 angles = collect(range(0, π, N))
 rop = RadonOp(eltype(image); angles, shape = size(image));
 sinogram = reshape(rop * vec(image), :, N)
diff --git a/docs/src/literate/tutorials/wavelet.jl b/docs/src/literate/tutorials/wavelet.jl
index 59e6943..9563a38 100644
--- a/docs/src/literate/tutorials/wavelet.jl
+++ b/docs/src/literate/tutorials/wavelet.jl
@@ -1,7 +1,7 @@
 # # Wavelet Operator
 include("../../util.jl") #hide
 # The wavelet operator is available when loading Wavelets.jl together with LinearOperatorCollection:
-using LinearOperatorCollection, Wavelets
+using Wavelets
 wop = WaveletOp(eltype(image), shape = (N, N))
 wavelet_image = reshape(wop * vec(image), N, N)
 fig = Figure()
diff --git a/docs/src/literate/tutorials/weighting.jl b/docs/src/literate/tutorials/weighting.jl
index 56c3174..e7c2d8b 100644
--- a/docs/src/literate/tutorials/weighting.jl
+++ b/docs/src/literate/tutorials/weighting.jl
@@ -4,15 +4,13 @@ include("../../util.jl") #hide
 # 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)
-
-# Such an operator is also implemented within LinearOperator.jl, however here this operator has a dedicated type on which one can dispatch:
 typeof(op)
 
-# And it is possible to retrieve the weights from the operator
+# 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)
+weighted_image = reshape(op * vec(image), N, N);
 
 # To visualize our weighted image, we will again use CairoMakie:
 fig = Figure()
diff --git a/docs/src/references.md b/docs/src/references.md
new file mode 100644
index 0000000..f525825
--- /dev/null
+++ b/docs/src/references.md
@@ -0,0 +1,24 @@
+```@docs
+SamplingOp
+WeightingOp
+DiagOp
+GradientOp
+ProdOp
+NormalOp
+normalOperator
+WaveletOp
+RadonOp
+NFFTOp
+FFTOp
+DCTOp
+DSTOp
+```
+
+
+
+Modules = [LinearOperatorCollection,
+isdefined(Base, :get_extension) ? Base.get_extension(LinearOperatorCollection, :LinearOperatorFFTWExt) : LinearOperatorCollection.LinearOperatorFFTWExt,
+isdefined(Base, :get_extension) ? Base.get_extension(LinearOperatorCollection, :LinearOperatorNFFTWExt) : LinearOperatorCollection.LinearOperatorNFFTWExt,
+isdefined(Base, :get_extension) ? Base.get_extension(LinearOperatorCollection, :LinearOperatoRadonWExt) : LinearOperatorCollection.LinearOperatorRadonExt,
+isdefined(Base, :get_extension) ? Base.get_extension(LinearOperatorCollection, :LinearOperatorWaveletExt) : LinearOperatorCollection.LinearOperatorWaveletExt,
+]
diff --git a/docs/src/tutorials/diagonal.jl b/docs/src/tutorials/diagonal.jl
deleted file mode 100644
index 2fe99f8..0000000
--- a/docs/src/tutorials/diagonal.jl
+++ /dev/null
@@ -1,18 +0,0 @@
-# # Block Diagonal Operator
-# 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
\ No newline at end of file
diff --git a/docs/src/tutorials/fft.jl b/docs/src/tutorials/fft.jl
deleted file mode 100644
index 3f00def..0000000
--- a/docs/src/tutorials/fft.jl
+++ /dev/null
@@ -1,18 +0,0 @@
-# # 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:
-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
\ No newline at end of file
diff --git a/docs/src/tutorials/gradient.jl b/docs/src/tutorials/gradient.jl
deleted file mode 100644
index caa1e03..0000000
--- a/docs/src/tutorials/gradient.jl
+++ /dev/null
@@ -1,10 +0,0 @@
-# # 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
\ No newline at end of file
diff --git a/docs/src/tutorials/normal.jl b/docs/src/tutorials/normal.jl
deleted file mode 100644
index d4899bc..0000000
--- a/docs/src/tutorials/normal.jl
+++ /dev/null
@@ -1,38 +0,0 @@
-# # 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.
\ No newline at end of file
diff --git a/docs/src/tutorials/overview.jl b/docs/src/tutorials/overview.jl
deleted file mode 100644
index 70ac349..0000000
--- a/docs/src/tutorials/overview.jl
+++ /dev/null
@@ -1,46 +0,0 @@
-# # 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(ComplexF64, 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
diff --git a/docs/src/tutorials/product.jl b/docs/src/tutorials/product.jl
deleted file mode 100644
index 6bd5f1d..0000000
--- a/docs/src/tutorials/product.jl
+++ /dev/null
@@ -1,32 +0,0 @@
-# # Product Operator
-
-# 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
\ No newline at end of file
diff --git a/docs/src/tutorials/weighting.jl b/docs/src/tutorials/weighting.jl
deleted file mode 100644
index 56c3174..0000000
--- a/docs/src/tutorials/weighting.jl
+++ /dev/null
@@ -1,22 +0,0 @@
-# # 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)
-
-# Such an operator is also implemented within LinearOperator.jl, however here this operator has a dedicated type on which one can dispatch:
-typeof(op)
-
-# And it is possible to 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
\ No newline at end of file
diff --git a/docs/src/util.jl b/docs/src/util.jl
new file mode 100644
index 0000000..f9fe089
--- /dev/null
+++ b/docs/src/util.jl
@@ -0,0 +1,11 @@
+using CairoMakie, LinearOperatorCollection
+function plot_image(figPos, img; title = "", width = 150, height = 150, colorscale = identity, colormap = :viridis)
+  ax = CairoMakie.Axis(figPos[1, 1]; yreversed=true, title, width, height)
+  hidedecorations!(ax)
+  hm = heatmap!(ax, img, colorscale = colorscale, colormap = colormap)
+  Colorbar(figPos[1, 2], hm)
+end
+N = 256
+using ImagePhantoms, ImageGeoms
+image = shepp_logan(N, SheppLoganToft())
+nothing
\ No newline at end of file
diff --git a/ext/LinearOperatorNFFTExt/NFFTOp.jl b/ext/LinearOperatorNFFTExt/NFFTOp.jl
index 17e185f..9662688 100644
--- a/ext/LinearOperatorNFFTExt/NFFTOp.jl
+++ b/ext/LinearOperatorNFFTExt/NFFTOp.jl
@@ -1,4 +1,15 @@
+"""
+    NFFTOpImpl(shape::Tuple, tr::Trajectory; kargs...)
+    NFFTOpImpl(shape::Tuple, tr::AbstractMatrix; kargs...)
 
+generates a `NFFTOpImpl` which evaluates the MRI Fourier signal encoding operator using the NFFT.
+
+# Arguments:
+* `shape::NTuple{D,Int64}`  - size of image to encode/reconstruct
+* `tr`                      - Either a `Trajectory` object, or a `ND x Nsamples` matrix for an ND-dimenensional (e.g. 2D or 3D) NFFT with `Nsamples` k-space samples
+* (`nodes=nothing`)         - Array containg the trajectory nodes (redundant)
+* (`kargs`)                 - additional keyword arguments
+"""
 function LinearOperatorCollection.NFFTOp(::Type{T};
     shape::Tuple, nodes::AbstractMatrix{U}, toeplitz=false, oversamplingFactor=1.25, 
    kernelSize=3, kargs...) where {U <: Number, T <: Number}
@@ -27,18 +38,6 @@ end
 
 LinearOperators.storage_type(op::NFFTOpImpl) = typeof(op.Mv5)
 
-"""
-    NFFTOpImpl(shape::Tuple, tr::Trajectory; kargs...)
-    NFFTOpImpl(shape::Tuple, tr::AbstractMatrix; kargs...)
-
-generates a `NFFTOpImpl` which evaluates the MRI Fourier signal encoding operator using the NFFT.
-
-# Arguments:
-* `shape::NTuple{D,Int64}`  - size of image to encode/reconstruct
-* `tr`                      - Either a `Trajectory` object, or a `ND x Nsamples` matrix for an ND-dimenensional (e.g. 2D or 3D) NFFT with `Nsamples` k-space samples
-* (`nodes=nothing`)         - Array containg the trajectory nodes (redundant)
-* (`kargs`)                 - additional keyword arguments
-"""
 function NFFTOpImpl(shape::Tuple, tr::AbstractMatrix{T}; toeplitz=false, oversamplingFactor=1.25, kernelSize=3, S = Vector{Complex{T}}, kargs...) where {T}
 
   baseArrayType = Base.typename(S).wrapper # https://github.com/JuliaLang/julia/issues/35543
diff --git a/ext/LinearOperatorRadonKAExt/RadonOp.jl b/ext/LinearOperatorRadonKAExt/RadonOp.jl
index e4d6132..954a201 100644
--- a/ext/LinearOperatorRadonKAExt/RadonOp.jl
+++ b/ext/LinearOperatorRadonKAExt/RadonOp.jl
@@ -1,3 +1,16 @@
+"""
+    RadonOp(::Type{T}; shape::NTuple{N, Int}, angles, geometry = RadonParallelCircle(shape[1], -(shape[1]-1)÷2:(shape[1]-1)÷2), μ = nothing, S = Vector{T}) where {T, N}
+
+Generates a `RadonOp` which evaluates the Radon transform operator and its adjoint (backprojection) for a given geometry and projection angles.
+
+# Arguments:
+* `T`                       - element type for the operator (e.g., `Float64`, `ComplexF32`)
+* `shape::NTuple{N, Int}`   - size of the image
+* `angles`                  - array of projection angles
+* `geometry`                - Radon geometry descriptor (default: parallel beam circle)
+* `μ`                       - optional attenuation map (for attenuated Radon transform)
+* `S`                       - storage type for internal vectors (default: `Vector{T}`)
+"""
 function LinearOperatorCollection.RadonOp(::Type{T}; shape::NTuple{N, Int}, angles,
    geometry = RadonParallelCircle(shape[1], -(shape[1]-1)÷2:(shape[1]-1)÷2), μ = nothing, S = Vector{T}) where {T, N}
   return RadonOpImpl(T; shape, angles, geometry, μ, S)
diff --git a/src/ProdOp.jl b/src/ProdOp.jl
index ccba07e..e33f5c3 100644
--- a/src/ProdOp.jl
+++ b/src/ProdOp.jl
@@ -1,12 +1,4 @@
 export ProdOp
-"""
-    `mutable struct ProdOp{T}`
-
-  struct describing the result of a composition/product of operators.
-  Describing the composition using a dedicated type has the advantage
-  that the latter can be made copyable. This is particularly relevant for
-  multi-threaded code
-"""
 mutable struct ProdOp{T,U,V, vecT <: AbstractVector{T}} <: AbstractLinearOperatorFromCollection{T}
   nrow :: Int
   ncol :: Int
@@ -31,7 +23,7 @@ end
 """
     ProdOp(A,B)
 
-composition/product of two Operators. Differs with * since it can handle normal operator
+composition/product of two Operators. Computes (A * (B * x)). Differs with composition via * since it can handle normal operator and allows for specialisation on `A` and `B`.
 """
 function ProdOp(A, B)
   nrow = size(A, 1)
diff --git a/src/SamplingOp.jl b/src/SamplingOp.jl
index e67997f..e57b0d7 100644
--- a/src/SamplingOp.jl
+++ b/src/SamplingOp.jl
@@ -1,5 +1,15 @@
 export vectorizePattern
 
+"""
+  SamplingOp(pattern::Array{Int}, shape::Tuple)
+
+builds a `LinearOperator` which only returns the vector elements at positions
+indicated by pattern.
+
+# Arguents
+* `pattern::Array{Int}` - indices to sample
+* `shape::Tuple`        - size of the array to sample
+"""
 function LinearOperatorCollection.SamplingOp(::Type{T};
   pattern::P, shape::Tuple=(), S = Vector{T}) where {P} where T <: Number
   if length(shape) == 0
@@ -18,16 +28,6 @@ function vectorizePattern(idx::T, shape::Tuple) where T<:AbstractArray{Int}
   return [ floor(Int,(i-1)/size(idx,1))*shape[1]+idx[i] for i = 1:length(idx) ]
 end
 
-"""
-  SamplingOp(pattern::Array{Int}, shape::Tuple)
-
-builds a `LinearOperator` which only returns the vector elements at positions
-indicated by pattern.
-
-# Arguents
-* `pattern::Array{Int}` - indices to sample
-* `shape::Tuple`        - size of the array to sample
-"""
 function SamplingOpImpl(T::Type{<:Number}, pattern::AbstractArray{Int}, shape::Tuple; S = Vector{T})
   ndims(pattern)>1 ?  idx = vectorizePattern(pattern, shape) : idx = pattern
   return opEye(T,length(idx), S = S)*opRestriction(idx, prod(shape); S = S)