using Sinograms: fbp
using MIRTjim: jim, prompt
using InteractiveUtils: versioninfo

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

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

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

rmax = maximum(r)
fovmask = @. sqrt(abs2(x) + abs2(y)') ≤ rmax
jim(x, y, fovmask, "FOV mask")

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

noisy_fbp_image = fbp(noisy_sinogram; dr)
noisy_fbp_image .*= fovmask # apply FOV mask
jim(x, y, noisy_fbp_image, "Noisy FBP image"; clim=(0,1))