Parker weights divide by 0 bug

[JeffFessler]

    @JeffFessler 

I just tried it out and now get an error:

DomainError with -Inf:

sin(x) is only defined for finite x.

sin_domain_error(::Float64)@trig.jl:28
sin(::Float64)@trig.jl:39
#34@parker.jl:88[inlined]
_broadcast_getindex_evalf@broadcast.jl:670[inlined]
_broadcast_getindex@broadcast.jl:643[inlined]
getindex@broadcast.jl:597[inlined]
macro expansion@broadcast.jl:961[inlined]
macro expansion@simdloop.jl:77[inlined]
copyto!@broadcast.jl:960[inlined]
copyto!@broadcast.jl:913[inlined]
materialize!@broadcast.jl:871[inlined]
materialize!@broadcast.jl:868[inlined]
parker_weight_fan_short(::Int64, ::Int64, ::Float64, ::Float64, ::StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}, ::Vector{Float64}, ::Float64, ::Type{Float32})@parker.jl:98
parker_weight_fan_short@parker.jl:73[inlined]
parker_weight(::Sinograms.SinoFanFlat{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64}, ::Type{Float32})@parker.jl:118
parker_weight@parker.jl:115[inlined]
_fbp_weights(::Sinograms.SinoFanFlat{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64})@plan2.jl:142
#plan_fbp#53@plan2.jl:119[inlined]
reco(::ImageMetadata.ImageMeta{Float32, 2, Matrix{Float32}, Dict{Symbol, Any}}, ::Tuple{Int64, Int64}, ::Tuple{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}}, ::Float64)@[Other: 34](http://localhost:1234/edit?id=315ae8ce-8aa2-11ed-1cdb-d15746ce3525#)
top-level scope@[Local: 13](http://localhost:1234/edit?id=315ae8ce-8aa2-11ed-1cdb-d15746ce3525#)

Apparently there is a division by zero in this line: https://github.com/JuliaImageRecon/Sinograms.jl/blob/main/src/fbp/parker.jl#L97

Not sure how to debug this further.

Originally posted by @tknopp in #47 (comment)

[JeffFessler]

@tknopp sorry about this. Can you summarize what ray geometry you used? I don't need a full MWE probably, just the part where you have a line of code like rg = CtFanArc( ... ).

[tknopp]

Sinograms.SinoFanFlat{Unitful.Quantity{Float64, 𝐋, Unitful.FreeUnits{(cm,), 𝐋, nothing}}, Float64} :
 nb::Int64 1536
 d::Unit{Float64} 0.011067708333333334 cm
 offset::Float32 -6.0
 na::Int64 834
 orbit::Float64 199.15104677834458
 orbit_start::Float64 -5.651046778344568
 source_offset::Unit{Float64} 0.0 cm
 dsd::Unit{Float64} 51.800000000000004 cm
 dod::Unit{Float64} 11.684172500000004 cm

[tknopp]

if I set orbit_start = 0 it works but the reconstruction is rotated. Most of the above parameters were read out of the file. Just the center and orbit_start were hand tuned.

[JeffFessler]

Those were very helpful clues. I have found the problem and it is in this line:

Sinograms.jl/src/fbp/parker.jl

Line 85 in ed697de

bet = ar .- ar[begin] # trick: force 0 start, so this ignores orbit_start!

Surprisingly the first value is slightly (eps) negative. Here is MWE to reproduce it. I will fix it now.

using Unitful: cm
using Sinograms
using Sinograms: _ar

rg = SinoFanFlat(; nb=1536, na = 834,
 d = 0.011067708333333334cm,
 orbit = 199.15104677834458,
 orbit_start = -5.651046778344568,
 dsd =  51.800000000000004cm,
 dod = 11.684172500000004cm,
 offset = -6,
)

ar = _ar(rg) # StepRangeLen
(ar .- ar[begin])[begin] # negative!?
minimum(ar .- minimum(ar)) # negative!?

[JeffFessler]

OK, #50 should fix this and it's so simple that I just merged it without review.
Could you please retest on main branch, and if it is OK then I will tag.
BTW, it makes me happy that the show method for the ray geometry had everything needed to reproduce 😄

[tknopp]

yep, thanks, the issue is fixed!