diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
deleted file mode 100644
index 252cc181..00000000
--- a/.github/workflows/CI.yml
+++ /dev/null
@@ -1,49 +0,0 @@
-name: CI
-on:
-  push:
-    paths-ignore:
-      - '*/*.md'
-      - 'docs/**'
-  pull_request:
-    paths-ignore:
-      - '*/*.md'
-      - 'docs/**'
-
-jobs:
-  test:
-    name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }}
-    runs-on: ${{ matrix.os }}
-    strategy:
-      fail-fast: false
-      matrix:
-        version:
-          - '1.6'
-          - 'nightly'
-        os:
-          - ubuntu-latest
-          - macOS-latest
-          - windows-latest
-        arch:
-          - x64
-    steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
-        with:
-          version: ${{ matrix.version }}
-          arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
-        env:
-          cache-name: cache-artifacts
-        with:
-          path: ~/.julia/artifacts
-          key: ${{ runner.os }}-test-${{ env.cache-name }}-${{ hashFiles('**/Project.toml') }}
-          restore-keys: |
-            ${{ runner.os }}-test-${{ env.cache-name }}-
-            ${{ runner.os }}-test-
-            ${{ runner.os }}-
-      - uses: julia-actions/julia-buildpkg@v1
-      - uses: julia-actions/julia-runtest@v1
-      - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
-        with:
-          file: lcov.info
diff --git a/.github/workflows/CompatHelper.yml b/.github/workflows/CompatHelper.yml
index 8d8c7bb4..98d5bc0a 100644
--- a/.github/workflows/CompatHelper.yml
+++ b/.github/workflows/CompatHelper.yml
@@ -1,4 +1,5 @@
 name: CompatHelper
+
 on:
   schedule:
     - cron: 0 0 * * 0 # weekly
@@ -8,13 +9,25 @@ jobs:
   CompatHelper:
     runs-on: ubuntu-latest
     steps:
-      - name: Pkg.add("CompatHelper")
-        run: julia -e 'using Pkg; Pkg.add("CompatHelper")'
-      - name: CompatHelper.main()
+      - name: "Add the General registry via Git"
+        run: |
+          import Pkg
+          ENV["JULIA_PKG_SERVER"] = ""
+          Pkg.Registry.add("General")
+        shell: julia --color=yes {0}
+      - name: "Install CompatHelper"
+        run: |
+          import Pkg
+          name = "CompatHelper"
+          uuid = "aa819f21-2bde-4658-8897-bab36330d9b7"
+          version = "3"
+          Pkg.add(; name, uuid, version)
+        shell: julia --color=yes {0}
+      - name: "Run CompatHelper"
+        run: |
+          import CompatHelper
+          CompatHelper.main()
+        shell: julia --color=yes {0}
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-          COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }} # need ssh
-        run: julia -e 'using CompatHelper; CompatHelper.main()'
-
-# based on:
-# https://github.com/JuliaRegistries/CompatHelper.jl
+          COMPATHELPER_PRIV: ${{ secrets.DOCUMENTER_KEY }}
diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
index 286c2d52..1123667c 100644
--- a/.github/workflows/Documentation.yml
+++ b/.github/workflows/Documentation.yml
@@ -2,7 +2,11 @@ name: Documentation
 
 on:
   pull_request:
+    paths-ignore:
+      - 'README.md'
   push:
+    paths-ignore:
+      - 'README.md'
     branches:
       - 'main'
       - 'release-'
@@ -15,8 +19,8 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@latest
         with:
-          version: '1.6'
-      - name: Cache artifacts
+          version: '1'
+      - name: CacheArtifacts
         uses: actions/cache@v1
         env:
           cache-name: cache-artifacts
@@ -27,13 +31,13 @@ jobs:
             ${{ runner.os }}-test-${{ env.cache-name }}-
             ${{ runner.os }}-test-
             ${{ runner.os }}-
-      - name: Install dependencies
+      - name: InstallDependencies
         run: |
           julia --project=docs/ -e '
-          ENV["JULIA_PKG_SERVER"] = ""
-          using Pkg
-          Pkg.develop(PackageSpec(path=pwd()))
-          Pkg.instantiate()'
+            ENV["JULIA_PKG_SERVER"] = ""
+            using Pkg
+            Pkg.develop(PackageSpec(path=pwd()))
+            Pkg.instantiate()'
       - name: BuildAndDeploy
         env:
 # https://juliadocs.github.io/Documenter.jl/stable/man/hosting/#Authentication:-GITHUB_TOKEN
diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
new file mode 100644
index 00000000..6f854ddb
--- /dev/null
+++ b/.github/workflows/ci.yml
@@ -0,0 +1,42 @@
+name: CI
+
+on:
+  create:
+    tags:
+  push:
+    branches:
+      - main
+    paths-ignore:
+      - '**.md'
+      - 'docs/**'
+  pull_request:
+    paths-ignore:
+      - '**.md'
+      - 'docs/**'
+  workflow_dispatch:
+
+jobs:
+  test:
+    name: Julia ${{ matrix.version }} - ${{ matrix.os }}
+    runs-on: ${{ matrix.os }}
+    continue-on-error: ${{ matrix.version == 'nightly' }}
+    strategy:
+      fail-fast: false
+      matrix:
+        version: ['1.6', '1', 'nightly']
+        os: [ubuntu-latest, windows-latest, macOS-latest]
+    steps:
+      - uses: actions/checkout@v2
+#     - name: "Set up Julia"
+      - uses: julia-actions/setup-julia@latest
+        with:
+          version: ${{ matrix.version }}
+#     - name: "Unit Test"
+      - uses: julia-actions/julia-buildpkg@latest
+      - uses: julia-actions/julia-runtest@latest
+#     - name: "Cover"
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        if: ${{ matrix.version == '1' && matrix.os == 'ubuntu-latest' }}
+        with:
+          file: lcov.info
diff --git a/.github/workflows/clean.yml b/.github/workflows/clean.yml
new file mode 100644
index 00000000..6d17bcbc
--- /dev/null
+++ b/.github/workflows/clean.yml
@@ -0,0 +1,28 @@
+# https://juliadocs.github.io/Documenter.jl/stable/man/hosting/#gh-pages-Branch
+
+name: DocPreviewCleanup
+
+on:
+  pull_request:
+    types: [closed]
+
+jobs:
+  doc-preview-cleanup:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout gh-pages branch
+        uses: actions/checkout@v2
+        with:
+          ref: gh-pages
+      - name: Delete preview and history + push changes
+        run: |
+            if [ -d "previews/PR$PRNUM" ]; then
+              git config user.name "Documenter.jl"
+              git config user.email "documenter@juliadocs.github.io"
+              git rm -rf "previews/PR$PRNUM"
+              git commit -m "delete preview"
+              git branch gh-pages-new $(echo "delete history" | git commit-tree HEAD^{tree})
+              git push --force origin gh-pages-new:gh-pages
+            fi
+        env:
+            PRNUM: ${{ github.event.number }}
diff --git a/.gitignore b/.gitignore
index 99d5d158..27686f92 100644
--- a/.gitignore
+++ b/.gitignore
@@ -16,6 +16,7 @@ deps/src/
 # Build artifacts for creating documentation generated by the Documenter package
 docs/build/
 docs/site/
+docs/src/generated/
 
 # File generated by Pkg, the package manager, based on a corresponding Project.toml
 # It records a fixed state of all packages used by the project. As such, it should not be
diff --git a/data/trained-cnn-example-6-dl.bson b/data/trained-cnn-example-6-dl.bson
new file mode 100644
index 00000000..add4ebe8
Binary files /dev/null and b/data/trained-cnn-example-6-dl.bson differ
diff --git a/docs/Project.toml b/docs/Project.toml
index a7e7e4df..95d81f11 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,19 +1,18 @@
 [deps]
+BSON = "fbb218c0-5317-5bc6-957e-2ee96dd4b1f0"
 Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Flux = "587475ba-b771-5e3f-ad9e-33799f191a9c"
 ImagePhantoms = "71a99df6-f52c-4da1-bd2a-69d6f37f3252"
+InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearMapsAA = "599c1a8e-b958-11e9-0d14-b1e6b2ecea07"
 Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MIRTjim = "170b2178-6dee-4cb0-8729-b3e8b57834cc"
+NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SPECTrecon = "ab1be465-a7f0-4423-9048-0ee774b70ed9"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 UnitfulRecipes = "42071c24-d89e-48dd-8a24-8a12d9b8861f"
-
-[compat]
-Distributions = "0.25"
-Documenter = "0.27.7"
-ImagePhantoms = "0.0.5"
-LinearMapsAA = "0.7"
-Literate = "2.9.3"
-MIRTjim = "0.13"
+Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
+ZygoteRules = "700de1a5-db45-46bc-99cf-38207098b444"
diff --git a/docs/lit/examples/1-overview.jl b/docs/lit/examples/1-overview.jl
index d0b8fea1..38c14740 100644
--- a/docs/lit/examples/1-overview.jl
+++ b/docs/lit/examples/1-overview.jl
@@ -15,6 +15,7 @@ using LinearAlgebra: mul!
 using LinearMapsAA: LinearMapAA
 using Plots: scatter, plot!, default; default(markerstrokecolor=:auto)
 using Plots # @animate, gif
+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.
@@ -216,3 +217,15 @@ anim = @animate for i in 1:nview
     )
 end
 gif(anim, "views.gif", fps = 8)
+
+
+# ### 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()
diff --git a/docs/lit/examples/3-psf.jl b/docs/lit/examples/3-psf.jl
index 71618240..81e4fc57 100644
--- a/docs/lit/examples/3-psf.jl
+++ b/docs/lit/examples/3-psf.jl
@@ -64,7 +64,7 @@ jim(image, "Original image")
 
 px = 11
 nview = 1 # for simplicity in this illustration
-psf = repeat(psf_gauss( ; ny, px), 1, 1, 1, nview)
+psf = repeat(psf_gauss( ; ny=nx, px), 1, 1, 1, nview)
 jim(psf, "PSF for each of $nx planes")
 
 
@@ -111,7 +111,7 @@ jim(adj, "Adjoint of PSF modeling")
 using LinearMapsAA: LinearMapAA
 
 nx, nz, px = 10, 7, 5 # small size for illustration
-psf3 = psf_gauss( ; ny, px)
+psf3 = psf_gauss( ; ny=nx, px)
 plan = plan_psf( ; nx, nz, px, T)
 idim = (nx,nx,nz)
 odim = (nx,nx,nz)
diff --git a/docs/lit/examples/4-mlem.jl b/docs/lit/examples/4-mlem.jl
index e956f66b..57eb5386 100644
--- a/docs/lit/examples/4-mlem.jl
+++ b/docs/lit/examples/4-mlem.jl
@@ -100,7 +100,7 @@ function mlem(x0, ynoisy, background, A; niter::Int = 20)
     all(>(0), background) || throw("need background > 0")
     x = copy(x0)
     asum = A' * ones(eltype(ynoisy), size(ynoisy))
-	time0 = time()
+    time0 = time()
     for iter = 1:niter
         @show iter, extrema(x), time() - time0
         ybar = A * x .+ background # forward model
@@ -119,7 +119,7 @@ function mlem!(x, ynoisy, background, A; niter::Int = 20)
     ybar = similar(ynoisy)
     yratio = similar(ynoisy)
     back = similar(x)
-	time0 = time()
+    time0 = time()
     for iter = 1:niter
         @show iter, extrema(x), time() - time0
         mul!(ybar, A, x)
diff --git a/docs/lit/examples/6-dl.jl b/docs/lit/examples/6-dl.jl
new file mode 100644
index 00000000..885beb03
--- /dev/null
+++ b/docs/lit/examples/6-dl.jl
@@ -0,0 +1,330 @@
+#---------------------------------------------------------
+# # [SPECTrecon deep learning use](@id 6-dl)
+#---------------------------------------------------------
+
+#=
+This page describes how to end-to-end train unrolled deep learning algorithms
+using the Julia package
+[`SPECTrecon`](https://github.com/JeffFessler/SPECTrecon.jl).
+
+This page was generated from a single Julia file:
+[6-dl.jl](@__REPO_ROOT_URL__/6-dl.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 # [`6-dl.ipynb`](@__NBVIEWER_ROOT_URL__/6-dl.ipynb),
+#md # and opened in [binder](https://mybinder.org/) here:
+#md # [`6-dl.ipynb`](@__BINDER_ROOT_URL__/6-dl.ipynb).
+
+
+# ### Setup
+
+# Packages needed here.
+
+using LinearAlgebra: norm, mul!
+using SPECTrecon
+using MIRTjim: jim, prompt
+using Plots: default; default(markerstrokecolor=:auto)
+using Zygote
+using ZygoteRules: @adjoint
+using Flux: Chain, Conv, SamePad, relu, params, unsqueeze
+using Flux # apparently needed for BSON @load
+using NNlib
+using LinearMapsAA: LinearMapAA
+using Distributions: Poisson
+using BSON: @load, @save
+import BSON # load
+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);
+
+
+#=
+### Overview
+
+Regularized expectation-maximization (reg-EM)
+is a commonly used algorithm for performing SPECT image reconstruction.
+This page considers regularizers of the form ``β/2 * ||x - u||^2``,
+where ``u`` is an auxiliary variable that often refers to the image denoised by a CNN.
+
+### Data generation
+
+Simulated data used in this page are identical to
+[`SPECTrecon ML-EM`](https://jefffessler.github.io/SPECTrecon.jl/stable/examples/4-mlem/).
+We repeat it again here for convenience.
+=#
+
+nx,ny,nz = 64,64,50
+T = Float32
+xtrue = zeros(T, nx,ny,nz)
+xtrue[(1nx÷4):(2nx÷3), 1ny÷5:(3ny÷5), 2nz÷6:(3nz÷6)] .= 1
+xtrue[(2nx÷5):(3nx÷5), 1ny÷5:(2ny÷5), 4nz÷6:(5nz÷6)] .= 2
+
+average(x) = sum(x) / length(x)
+function mid3(x::AbstractArray{T,3}) where {T}
+    (nx,ny,nz) = size(x)
+    xy = x[:,:,ceil(Int, nz÷2)]
+    xz = x[:,ceil(Int,end/2),:]
+    zy = x[ceil(Int, nx÷2),:,:]'
+    return [xy xz; zy fill(average(xy), nz, nz)]
+end
+jim(mid3(xtrue), "Middle slices of xtrue")
+
+
+# ### PSF
+
+# Create a synthetic depth-dependent PSF for a single view
+px = 11
+psf1 = psf_gauss( ; ny, px)
+jim(psf1, "PSF for each of $ny planes")
+
+
+# In general the PSF can vary from view to view
+# due to non-circular detector orbits.
+# For simplicity, here we illustrate the case
+# where the PSF is the same for every view.
+
+nview = 60
+psfs = repeat(psf1, 1, 1, 1, nview)
+size(psfs)
+
+
+# ### SPECT system model using `LinearMapAA`
+
+dy = 8 # transaxial pixel size in mm
+mumap = zeros(T, size(xtrue)) # zero μ-map just for illustration here
+plan = SPECTplan(mumap, psfs, dy; T)
+
+forw! = (y,x) -> project!(y, x, plan)
+back! = (x,y) -> backproject!(x, y, plan)
+idim = (nx,ny,nz)
+odim = (nx,nz,nview)
+A = LinearMapAA(forw!, back!, (prod(odim),prod(idim)); T, odim, idim)
+
+
+# Generate noisy data
+
+if !@isdefined(ynoisy) # generate (scaled) Poisson data
+    ytrue = A * xtrue
+    target_mean = 20 # aim for mean of 20 counts per ray
+    scale = target_mean / average(ytrue)
+    scatter_fraction = 0.1 # 10% uniform scatter for illustration
+    scatter_mean = scatter_fraction * average(ytrue) # uniform for simplicity
+    ynoisy = rand.(Poisson.(scale * (ytrue .+ scatter_mean))) / scale
+end
+jim(ynoisy, "$nview noisy projection views"; ncol=10)
+
+
+# ### ML-EM algorithm
+function mlem!(x, ynoisy, background, A; niter::Int = 20)
+    all(>(0), background) || throw("need background > 0")
+    asum = A' * ones(eltype(ynoisy), size(ynoisy)) # this allocates
+    ybar = similar(ynoisy)
+    yratio = similar(ynoisy)
+    back = similar(x)
+    time0 = time()
+    for iter = 1:niter
+        if isinteractive()
+            @show iter, extrema(x), time() - time0
+        end
+        mul!(ybar, A, x)
+        @. yratio = ynoisy / (ybar + background) # coalesce broadcast!
+        mul!(back, A', yratio) # back = A' * (ynoisy / ybar)
+        @. x *= back / asum # multiplicative update
+    end
+    return x
+end
+
+
+# Apply MLEM algorithm to simulated data
+x0 = ones(T, nx, ny, nz) # initial uniform image
+
+niter = 30
+
+if !@isdefined(xhat1)
+    xhat1 = copy(x0)
+    mlem!(xhat1, ynoisy, scatter_mean, A; niter)
+end;
+
+# Define evaluation metric
+nrmse(x) = round(100 * norm(mid3(x) - mid3(xtrue)) / norm(mid3(xtrue)); digits=1)
+prompt()
+## jim(mid3(xhat1), "MLEM NRMSE=$(nrmse(xhat1))%") # display ML-EM reconstructed image
+
+
+# ### Implement a 3D CNN denoiser
+
+cnn = Chain(
+    Conv((3,3,3), 1 => 4, relu; stride = 1, pad = SamePad(), bias = true),
+    Conv((3,3,3), 4 => 4, relu; stride = 1, pad = SamePad(), bias = true),
+    Conv((3,3,3), 4 => 1; stride = 1, pad = SamePad(), bias = true),
+)
+# Show how many parameters the CNN has
+paramCount = sum([sum(length, params(layer)) for layer in cnn])
+
+
+
+#=
+### Custom backpropagation
+
+Forward and back-projection are linear operations
+so their Jacobians are very simple
+and there is no need to auto-differentiate through the system matrix
+and that would be very computationally expensive.
+Instead, we tell Flux.jl to use the customized Jacobian when doing backpropagation.
+=#
+
+projectb(x) = A * x
+@adjoint projectb(x) = A * x, dy -> (A' * dy, )
+
+backprojectb(y) = A' * y
+@adjoint backprojectb(y) = A' * y, dx -> (A * dx, )
+
+
+# ### Backpropagatable regularized EM algorithm
+# First define a function for unsqueezing the data
+# because Flux CNN model expects a 5-dim tensor
+function unsqueeze45(x)
+    return unsqueeze(unsqueeze(x, 4), 5)
+end
+
+"""
+    bregem(projectb, backprojectb, y, r, Asum, x, cnn, β; niter = 1)
+
+Backpropagatable regularized EM reconstruction with CNN regularization
+-`projectb`: backpropagatable forward projection
+-`backprojectb`: backpropagatable backward projection
+-`y`: projections
+-`r`: scatters
+-`Asum`: A' * 1
+-`x`: current iterate
+-`cnn`: the CNN model
+-`β`: regularization parameter
+-`niter`: number of iteration for inner EM
+"""
+function bregem(
+    projectb::Function,
+    backprojectb::Function,
+    y::AbstractArray,
+    r::AbstractArray,
+    Asum::AbstractArray,
+    x::AbstractArray,
+    cnn::Union{Chain,Function},
+    β::Real;
+    niter::Int = 1,
+)
+
+    u = cnn(unsqueeze45(x))[:,:,:,1,1]
+    Asumu = Asum - β * u
+    Asumu2 = Asumu.^2
+    T = eltype(x)
+    for iter = 1:niter
+        eterm = backprojectb((y ./ (projectb(x) + r)))
+        eterm_beta = 4 * β * (x .* eterm)
+        x = max.(0, T(1/2β) * (-Asumu + sqrt.(Asumu2 + eterm_beta)))
+    end
+    return x
+end
+
+
+# ### Loss function
+# We set β = 1 and train 2 outer iterations.
+
+β = 1
+Asum = A' * ones(T, nx, nz, nview)
+scatters = scatter_mean * ones(T, nx, nz, nview)
+function loss(xrecon, xtrue)
+    xiter1 = bregem(projectb, backprojectb, ynoisy, scatters,
+                    Asum, xrecon, cnn, β; niter = 1)
+    xiter2 = bregem(projectb, backprojectb, ynoisy, scatters,
+                    Asum, xiter1, cnn, β; niter = 1)
+    return sum(abs2, xiter2 - xtrue)
+end
+
+
+# Initial loss
+@show loss(xhat1, xtrue)
+
+# ### Train the CNN
+# Uncomment the following code to train!
+## using Printf
+## nepoch = 200
+## for e = 1:nepoch
+##     @printf("epoch = %d, loss = %.2f\n", e, loss(xhat1, xtrue))
+##     ps = Flux.params(cnn)
+##     gs = gradient(ps) do
+##         loss(xhat1, xtrue) # we start with the 30 iteration EM reconstruction
+##     end
+##     opt = ADAMW(0.002)
+##     Flux.Optimise.update!(opt, ps, gs)
+## end
+
+# Uncomment to save your trained model.
+## file = "../data/trained-cnn-example-6-dl.bson" # adjust path/name as needed
+## @save file cnn
+
+# Load the pre-trained model (uncomment if you save your own model).
+## @load file cnn
+
+#=
+The code below here works fine when run via `include` from the REPL,
+but it fails with the error `UndefVarError: NNlib not defined`
+on the `BSON.load` step when run via Literate/Documenter.
+So for now it is just fenced off with `isinteractive()`.
+=#
+
+if isinteractive()
+    url = "https://github.com/JeffFessler/SPECTrecon.jl/blob/main/data/trained-cnn-example-6-dl.bson?raw=true"
+    tmp = tempname()
+    download(url, tmp)
+    cnn = BSON.load(tmp)[:cnn]
+else
+    cnn = x -> x # fake "do-nothing CNN" for Literate/Documenter version
+end
+
+# Perform recon with pre-trained model.
+xiter1 = bregem(projectb, backprojectb, ynoisy, scatters,
+                Asum, xhat1, cnn, β; niter = 1)
+xiter2 = bregem(projectb, backprojectb, ynoisy, scatters,
+                Asum, xiter1, cnn, β; niter = 1)
+
+clim = (0,2)
+jim(
+    jim(mid3(xtrue), "xtrue"; clim),
+    jim(mid3(xhat1), "EM recon, NRMSE = $(nrmse(xhat1))%"; clim),
+    jim(mid3(xiter1), "Iter 1, NRMSE = $(nrmse(xiter1))%"; clim),
+    jim(mid3(xiter2), "Iter 2, NRMSE = $(nrmse(xiter2))%"; clim),
+)
+
+#=
+For the web-based Documenter/Literate version,
+the three NRMSE values will be the same
+because of the "do-nothing" CNN above.
+But if you download this file and run it locally,
+then you will see that the CNN reduces the NRMSE.
+
+A more thorough investigation
+would compare the CNN approach
+to a suitably optimized regularized approach;
+see [http://doi.org/10.1109/EMBC46164.2021.9630985](http://doi.org/10.1109/EMBC46164.2021.9630985).
+=#
+
+
+# ### 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()
diff --git a/docs/make.jl b/docs/make.jl
index df7dcb4b..2edcd952 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,50 +1,62 @@
-using SPECTrecon
+execute = isempty(ARGS) || ARGS[1] == "run"
+
+org, reps = :JeffFessler, :SPECTrecon
+eval(:(using $reps))
 using Documenter
 using Literate
 
-# based on:
-# https://github.com/jw3126/UnitfulRecipes.jl/blob/master/docs/make.jl
-
 # https://juliadocs.github.io/Documenter.jl/stable/man/syntax/#@example-block
 ENV["GKSwstype"] = "100"
+ENV["GKS_ENCODING"] = "utf-8"
 
-# generate tutorials and how-to guides using Literate
+# generate examples using Literate
 lit = joinpath(@__DIR__, "lit")
 src = joinpath(@__DIR__, "src")
-notebooks = joinpath(src, "notebooks")
+gen = joinpath(@__DIR__, "src/generated")
 
-ENV["GKS_ENCODING"] = "utf-8"
+base = "$org/$reps.jl"
+repo_root_url =
+    "https://github.com/$base/blob/main/docs/lit/examples"
+nbviewer_root_url =
+    "https://nbviewer.org/github/$base/tree/gh-pages/dev/generated/examples"
+binder_root_url =
+    "https://mybinder.org/v2/gh/$base/gh-pages?filepath=dev/generated/examples"
 
-DocMeta.setdocmeta!(SPECTrecon, :DocTestSetup, :(using SPECTrecon); recursive=true)
 
-execute = true # Set to true for executing notebooks and documenter!
-nb = false # Set to true to generate the notebooks
+repo = eval(:($reps))
+DocMeta.setdocmeta!(repo, :DocTestSetup, :(using $reps); recursive=true)
+
+#ENV["JULIA_DEBUG"] = "Literate"
+
 for (root, _, files) in walkdir(lit), file in files
-    splitext(file)[2] == ".jl" || continue
+    splitext(file)[2] == ".jl" || continue # process .jl files only
     ipath = joinpath(root, file)
-    opath = splitdir(replace(ipath, lit=>src))[1]
-    Literate.markdown(ipath, opath, documenter = execute)
-    nb && Literate.notebook(ipath, notebooks, execute = execute)
+    opath = splitdir(replace(ipath, lit => gen))[1]
+    Literate.markdown(ipath, opath, documenter = execute; # run examples
+        repo_root_url, nbviewer_root_url, binder_root_url)
+    Literate.notebook(ipath, opath; execute = false, # no-run notebooks
+        repo_root_url, nbviewer_root_url, binder_root_url)
 end
 
+
 # Documentation structure
 ismd(f) = splitext(f)[2] == ".md"
 pages(folder) =
-    [joinpath(folder, f) for f in readdir(joinpath(src, folder)) if ismd(f)]
+    [joinpath("generated/", folder, f) for f in readdir(joinpath(gen, folder)) if ismd(f)]
 
 isci = get(ENV, "CI", nothing) == "true"
 
 format = Documenter.HTML(;
     prettyurls = isci,
     edit_link = "main",
-    canonical = "https://JeffFessler.github.io/SPECTrecon.jl/stable/",
-#   assets = String[],
+    canonical = "https://$org.github.io/$repo.jl/stable/",
+    assets = ["assets/custom.css"],
 )
 
 makedocs(;
-    modules = [SPECTrecon],
+    modules = [repo],
     authors = "Jeff Fessler & Zongyu Li & contributors",
-    sitename = "SPECTrecon.jl",
+    sitename = "$repo.jl",
     format,
     pages = [
         "Home" => "index.md",
@@ -56,14 +68,12 @@ makedocs(;
 
 if isci
     deploydocs(;
-        repo = "github.com/JeffFessler/SPECTrecon.jl",
+        repo = "github.com/$base",
         devbranch = "main",
         devurl = "dev",
         versions = ["stable" => "v^", "dev" => "dev"],
         forcepush = true,
-        push_preview = true,
-        # see https://JeffFessler.github.io/SPECTrecon.jl/previews/PR##
+#       push_preview = true,
+        # see https://$org.github.io/$repo.jl/previews/PR##
     )
-else
-    @warn "may need to: rm -r src/examples"
 end
diff --git a/docs/src/assets/custom.css b/docs/src/assets/custom.css
new file mode 100644
index 00000000..31e65493
--- /dev/null
+++ b/docs/src/assets/custom.css
@@ -0,0 +1,3 @@
+.docs-sourcelink {
+    opacity:1 !important;
+}
diff --git a/docs/src/index.md b/docs/src/index.md
index ca01de78..dff45ef2 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -37,6 +37,4 @@ The forward projection method works as follows.
 The back-projection method is the exact adjoint
 of the forward projector.
 
-See the
-[Examples](@ref 1-overview)
-tab to the left for usage details.
+See the "Examples" for details.
diff --git a/src/plan-rotate.jl b/src/plan-rotate.jl
index d419d2b1..af8f4f89 100644
--- a/src/plan-rotate.jl
+++ b/src/plan-rotate.jl
@@ -36,18 +36,19 @@ Struct for storing work arrays and factors for 2D square image rotation for one
 - `workmat2 [nx + 2 * padsize, nx + 2 * padsize]` padded work matrix
 - `interp::SparseInterpolator{T, 2, 1}`
 """
-struct PlanRotate{T, R}
+struct PlanRotate{T, R, A2}
     nx::Int
 #   method::Symbol
     padsize::Int
-    workmat1::Matrix{T}
-    workmat2::Matrix{T}
+    workmat1::A2
+    workmat2::A2
     interp::SparseInterpolator{T, 2, 1}
 
     function PlanRotate(
         nx::Int ;
         T::DataType = Float32,
         method::Symbol = :two, # :one is for 1d interpolation, :two is for 2d interpolation
+        arraytype::DataType = Matrix{T},
     )
 
         (method == :one || method == :two) || @warn("method $method") # check interp method
@@ -66,7 +67,7 @@ struct PlanRotate{T, R}
 
         interp = SparseInterpolator(LinearSpline(T), zeros(T, size(workmat1, 1)), size(workmat1, 1))
 
-        new{T, RotateMode{method}}(
+        new{T, RotateMode{method}, arraytype}(
             nx,
             padsize,
             workmat1,
@@ -98,8 +99,9 @@ function plan_rotate(
     T::DataType = Float32,
     method::Symbol = :two,
     nthread::Int = Threads.nthreads(),
+    arraytype::DataType = Matrix{T},
 )
-    return [PlanRotate(nx; T, method) for id = 1:nthread]
+    return [PlanRotate(nx; T, method, arraytype) for id = 1:nthread]
 end
 
 
diff --git a/src/spectplan.jl b/src/spectplan.jl
index d39c1804..7c1ea0d9 100644
--- a/src/spectplan.jl
+++ b/src/spectplan.jl
@@ -29,34 +29,34 @@ Currently code assumes the following:
 * `psf` is symmetric
 * multiprocessing using # of threads specified by `Threads.nthreads()`
 """
-struct SPECTplan{T}
+struct SPECTplan{T,A2,A3,A4,U,R}
     T::DataType # default type for work arrays etc.
     imgsize::NTuple{3, Int}
     px::Int
     pz::Int
-    imgr::Union{Array{T, 3}, Vector{Array{T, 3}}} # 3D rotated image, (nx, ny, nz)
-    add_img::Union{Array{T, 3}, Vector{Array{T, 3}}}
-    mumap::Array{T, 3} # [nx,ny,nz] attenuation map, must be 3D, possibly zeros()
-    mumapr::Union{Array{T, 3}, Vector{Array{T, 3}}} # 3D rotated mumap, (nx, ny, nz)
-    exp_mumapr::Vector{Matrix{T}} # 2D exponential rotated mumap, (nx, ny)
-    psfs::Array{T, 4} # PSFs must be 4D, [px, pz, ny, nview], finally be centered psf
+    imgr::U # 3D rotated image, (nx, ny, nz)
+    add_img::U
+    mumap::A3 # [nx,ny,nz] attenuation map, must be 3D, possibly zeros()
+    mumapr::U # 3D rotated mumap, (nx, ny, nz)
+    exp_mumapr::Vector{A2} # 2D exponential rotated mumap, (nx, ny)
+    psfs::A4 # PSFs must be 4D, [px, pz, ny, nview], finally be centered psf
     nview::Int # number of views
     viewangle::StepRangeLen{T}
     interpmeth::Symbol
     mode::Symbol
     dy::T
     nthread::Int # number of threads
-    planrot::Vector{PlanRotate}
-    planpsf::Vector{PlanPSF}
+    planrot::Vector{PlanRotate{T, R, A2}}
+    planpsf::Vector{PlanPSF} # todo: concrete?
 
     """
         SPECTplan(mumap, psfs, dy; T, viewangle, interpmeth, nthread, mode)
     """
     function SPECTplan(
-        mumap::Array{<:RealU, 3},
-        psfs::Array{<:RealU, 4},
+        mumap::AbstractArray{<:RealU,3},
+        psfs::AbstractArray{<:RealU,4},
         dy::RealU;
-        T::DataType = promote_type(eltype(mumap), Float32),
+        T::DataType = promote_type(eltype(mumap), eltype(psfs), Float32),
         viewangle::StepRangeLen{<:RealU} = (0:size(psfs, 4) - 1) / size(psfs, 4) * T(2π), # set of view angles
         interpmeth::Symbol = :two, # :one is for 1d interpolation, :two is for 2d interpolation
         nthread::Int = Threads.nthreads(),
@@ -105,13 +105,19 @@ struct SPECTplan{T}
             mumapr = Array{T, 3}(undef, nx, ny, nz)
         end
 
+        U = typeof(mumapr)
+        A2 = typeof(mumap[:,:,1])
+        A3 = typeof(mumap)
+        A4 = typeof(psfs)
+
         exp_mumapr = [Matrix{T}(undef, nx, nz) for id = 1:nthread]
 
-        planrot = plan_rotate(nx; T, method = interpmeth)
+        planrot = plan_rotate(nx; T, method = interpmeth, arraytype=A2)
 
         planpsf = plan_psf(; nx, nz, px, pz, nthread, T)
 
-        new{T}(T, # default type for work arrays etc.
+        new{T, A2, A3, A4, U, RotateMode{interpmeth}}(
+            T, # default type for work arrays etc.
             imgsize,
             px,
             pz,
diff --git a/time/rotate3.jl b/time/rotate3.jl
index 8d306424..ec4833d7 100644
--- a/time/rotate3.jl
+++ b/time/rotate3.jl
@@ -29,12 +29,12 @@ if !@isdefined(elapse0)
 end
 
 if !@isdefined(elapse)
-	ntask = 1:20
-	elapse = zeros(length(ntask))
-	for (i,nt) in enumerate(ntask)
-		@show i, nt
-    	elapse[i] = @belapsed imrotate!($out1, $image3, $θ, $plans, $nt)
-	end
+    ntask = 1:20
+    elapse = zeros(length(ntask))
+    for (i,nt) in enumerate(ntask)
+        @show i, nt
+        elapse[i] = @belapsed imrotate!($out1, $image3, $θ, $plans, $nt)
+    end
 end
 
 @show minimum(elapse), elapse0