using Plots: plot, gui # these 2 must precede Sinograms for Requires to work!
using Unitful: cm
using Sinograms: CtFanArc, CtFanFlat # CtPar
using Sinograms: rays, plan_fbp, Window, Hamming, fdk, ct_geom_plot3
using ImageGeoms: ImageGeom, MaskCircle, fovs
using ImagePhantoms: ellipsoid_parameters, ellipsoid
using ImagePhantoms: radon, phantom
using MIRTjim: jim, prompt

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

ig = ImageGeom(MaskCircle(); dims=(64,62,30), deltas = (1,1,2) .* 0.4cm)

μ = 0.1 / cm # typical linear attenuation coefficient
params = ellipsoid_parameters( ; fovs = fovs(ig), u = (1, 1, μ))
ob = ellipsoid(params)

clim = (0.95, 1.05) .* μ

oversample = 3
true_image = phantom(axes(ig)..., ob, oversample)
jim(axes(ig), true_image, "True 3D Shepp-Logan phantom image"; clim)

p = (ns = 130, ds = 0.3cm, nt = 80, dt = 0.4cm, na = 50, dsd = 200cm, dod = 40cm)
cg = CtFanArc( ; p...)

ct_geom_plot3(cg, ig)

prompt()

proj_arc = radon(rays(cg), ob)
jim(cg.s, cg.t, proj_arc ;
    title="Shepp-Logan projections (arc)", xlabel="s", ylabel="t")

plan = plan_fbp(cg, ig; window = Window(Hamming(), 1.0))
fdk_arc = fdk(plan, proj_arc)
jim(axes(ig), fdk_arc, "FDK image (arc)"; clim)

err_arc = fdk_arc - true_image
jim(axes(ig), err_arc, "Error image (arc)"; clim = (-1,1) .* (0.05μ))

cg = CtFanFlat( ; p...)
proj_flat = radon(rays(cg), ob)
jim(cg.s, cg.t, proj_flat ;
    title="Shepp-Logan projections (flat)", xlabel="s", ylabel="t")

plan = plan_fbp(cg, ig; window = Window(Hamming(), 1.0))
fdk_flat = fdk(plan, proj_flat)
jim(axes(ig), fdk_flat, "FDK image (flat)"; clim)

err_flat = fdk_flat - true_image
jim(axes(ig), err_flat, "Error image (flat)"; clim = (-1,1) .* (0.05μ))