using ImageGeoms: ImageGeom, fovs, MaskCircle
import ImageGeoms # downsample
using ImagePhantoms: shepp_logan, SheppLoganToft, radon, phantom
using MIRTjim: jim, prompt
import Plots
using Sinograms: project_bdd, backproject_bdd
import Sinograms # downsample, _ar, _dso
using Sinograms: SinoFanFlat, rays, plan_fbp, fbp, Window, Hamming, sino_geom_plot!, angles
using Unitful: mm

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

down = 4 # save time
rg = SinoFanFlat( ; nb = 910, d = 1.0239mm, na = 360, dsd = 949mm, dod = 408mm)
rg = Sinograms.downsample(rg, down)

ig = ImageGeom(MaskCircle(); dims=(512,512), deltas = (1mm,1mm))
ig = ImageGeoms.downsample(ig, down)

geo = (DSD = rg.dsd, DS0 = Sinograms._dso(rg), pSize = ig.deltas[1],
    dSize = rg.d, nPix = ig.dims[1], nDet = rg.nb, angle = Sinograms._ar(rg))

pa = jim(axes(ig), ig.mask; prompt=false)
sino_geom_plot!(rg, ig)

prompt()

μ = 0.01 / mm # typical linear attenuation coefficient
ob = shepp_logan(SheppLoganToft(); fovs = fovs(ig), u = (1, 1, μ))
testimage = phantom(axes(ig)..., ob) # Create a phantom image
pt = jim(axes(ig), testimage) # Show the true phantom image

if !@isdefined(sinogramR)
    @time sinogramR = radon(rays(rg), ob)
    # Ensure that sinogram is not truncated
    @assert all(==(0), sum(abs, sinogramR, dims=2)[[1,end]])
end;

if !@isdefined(sinogramB)
    @time sinogramB = project_bdd(reverse(rot180(testimage'), dims=2), geo)
    sinogramB = sinogramB' # todo
end;

p1 = jim(axes(rg), sinogramB; prompt=false,
 title="bdd_2d sinogram", xlabel="r", ylabel="ϕ")
p2 = jim(axes(rg), sinogramR; prompt=false,
 title="Radon test sinogram", xlabel="r", ylabel="ϕ")
p12 = jim(p1,p2)

pd = jim(axes(rg), sinogramR - sinogramB;
 title="Difference sinogram", xlabel="r", ylabel="ϕ")

if !@isdefined(imageB)
    @time imageB = backproject_bdd(sinogramB'/1mm, geo) # todo
    imageB = rotr90(imageB) # todo
end
pb = jim(axes(ig), imageB;
 title="bdd_2d backprojection", xlabel="x", ylabel="y")

plan = plan_fbp(rg, ig; window = Window(Hamming(), 1.0))
fbp_image_b = fbp(plan, sinogramB)
fbp_image_r = fbp(plan, sinogramR);

clim = (0.04, 1.1) .* μ
r1 = jim(axes(ig), fbp_image_b, "FBP image with bdd_2d"; clim)
r2 = jim(axes(ig), fbp_image_r, "FBP image with radon"; clim)
r12 = jim(r1,r2; size=(1000,300))

pe = jim(axes(ig), fbp_image_b - fbp_image_r, "Difference image")