fix: add notes to documentation about shift #45

kunzaatko · 2023-07-18T14:45:54Z

[Ticket: #44]

kunzaatko · 2023-07-18T14:46:25Z

I left the other formatting changes... They were mostly deletions of end of line spaces.

roflmaostc · 2023-07-18T14:50:15Z

src/deconvolution.jl

@@ -38,6 +38,10 @@ regularizers and mappings.
                     and extended with a mean value (if border regions are used).
 - `debug_f=nothing`: A debug function which must take a single argument, the current reconstruction.

+!!! note
+    If you want to provide your PSF model, ensure that it is not centered. You may need to use `fftshift` for a PSF


maybe put: "ensure that it is centered around the first entry of the array (psf[1]).
And: "You may need to use ifftshift".

Shifting from center size(psf, 1) \div 2 + 1 to 1 requires a `ifftshift

[Ticket: #44]

kunzaatko · 2023-07-18T14:58:42Z

One more thing... The PSF that is created is normalized by its sum... I guess that is all right (although I saw that most literature uses one that has 1 at the center). I found out that if you do not normalize it, the automatic differentiation fails. (I tried it with RL deconvolution, but I think that it would be same elsewhere). Do you know where that originates from?

roflmaostc · 2023-07-18T15:01:19Z

I think it probably should still work. Maybe the initial guess should be divided by a large number to compensate for this.

kunzaatko · 2023-07-18T15:08:54Z

I am getting this error:

In [219]: DeconvOptim.richardson_lucy_iterative(Float32.(lr_imgs[:,:,1]), Float32.(ifftshift(abs.(no_offset_view(psf(tf, 512, Δxy))))))
ERROR: DomainError with -0.015415873:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float32)
    @ Base.Math ./math.jl:33
  [2] sqrt
    @ ./math.jl:677 [inlined]
  [3] sqrt
    @ ~/.julia/packages/ForwardDiff/vXysl/src/dual.jl:240 [inlined]
  [4] #1384
    @ ~/.julia/packages/Zygote/JeHtr/src/lib/broadcast.jl:276 [inlined]
  [5] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
  [6] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
  [7] getindex
    @ ./broadcast.jl:610 [inlined]
  [8] macro expansion
    @ ./broadcast.jl:974 [inlined]
  [9] macro expansion
    @ ./simdloop.jl:77 [inlined]
 [10] copyto!
    @ ./broadcast.jl:973 [inlined]
 [11] copyto!
    @ ./broadcast.jl:926 [inlined]
 [12] copy
    @ ./broadcast.jl:898 [inlined]
 [13] materialize
    @ ./broadcast.jl:873 [inlined]
 [14] broadcast_forward
    @ ~/.julia/packages/Zygote/JeHtr/src/lib/broadcast.jl:282 [inlined]
 [15] _broadcast_generic
    @ ~/.julia/packages/Zygote/JeHtr/src/lib/broadcast.jl:212 [inlined]
 [16] adjoint
    @ ~/.julia/packages/Zygote/JeHtr/src/lib/broadcast.jl:205 [inlined]
 [17] _pullback(__context__::Zygote.Context{false}, 702::typeof(Base.Broadcast.broadcasted), 703::Base.Broadcast.DefaultArrayStyle{2}, f::typeof(sqrt), args::Matrix{Fl
oat32})
    @ Zygote ~/.julia/packages/ZygoteRules/OgCVT/src/adjoint.jl:66
 [18] _apply(::Function, ::Vararg{Any})
    @ Core ./boot.jl:838
 [19] adjoint
    @ ~/.julia/packages/Zygote/JeHtr/src/lib/lib.jl:203 [inlined]
 [20] _pullback
    @ ~/.julia/packages/ZygoteRules/OgCVT/src/adjoint.jl:66 [inlined]
 [21] _pullback
    @ ./broadcast.jl:1311 [inlined]
 [22] _pullback
    @ ~/.julia/packages/DeconvOptim/d984T/src/regularizer.jl:249 [inlined]
 [23] _pullback(ctx::Zygote.Context{false}, f::DeconvOptim.var"#594#595", args::Matrix{Float32})
    @ Zygote ~/.julia/packages/Zygote/JeHtr/src/compiler/interface2.jl:0
 [24] pullback(f::Function, cx::Zygote.Context{false}, args::Matrix{Float32})
    @ Zygote ~/.julia/packages/Zygote/JeHtr/src/compiler/interface.jl:44
 [25] pullback
    @ ~/.julia/packages/Zygote/JeHtr/src/compiler/interface.jl:42 [inlined]
 [26] gradient(f::Function, args::Matrix{Float32})
    @ Zygote ~/.julia/packages/Zygote/JeHtr/src/compiler/interface.jl:96
 [27] #invokelatest#2
    @ ./essentials.jl:816 [inlined]
 [28] invokelatest
    @ ./essentials.jl:813 [inlined]
 [29] (::DeconvOptim.var"#∇reg#149"{DeconvOptim.var"#594#595", Matrix{Float32}})(x::Matrix{Float32})
    @ DeconvOptim ~/.julia/packages/DeconvOptim/d984T/src/lucy_richardson.jl:57
 [30] (::DeconvOptim.var"#iter_with_reg#151"{Float64, DeconvOptim.var"#iter_without_reg#150"{Matrix{Float32}, Matrix{Float32}, Matrix{ComplexF32}, DeconvOptim.var"#co
v#137"{Matrix{ComplexF32}, Matrix{Float32}, AbstractFFTs.ScaledPlan{ComplexF32, FFTW.rFFTWPlan{ComplexF32, 1, false, 2, UnitRange{Int64}}, Float32}, Matrix{ComplexF32
, FFTW.rFFTWPlan{Float32, -1, false, 2, UnitRange{Int64}}}, Matrix{ComplexF32}}, Matrix{Float32}, DeconvOptim.var"#∇reg#149"{DeconvOptim.var"#594#595", Matrix{Float32
}})(rec::Matrix{Float32})
    @ DeconvOptim ~/.julia/packages/DeconvOptim/d984T/src/lucy_richardson.jl:65
 [31] richardson_lucy_iterative(measured::Matrix{Float32}, psf::Matrix{Float32}; regularizer::Function, λ::Float64, iterations::Int64, conv_dims::UnitRange{Int64}, pr
gress::Nothing)
    @ DeconvOptim ~/.julia/packages/DeconvOptim/d984T/src/lucy_richardson.jl:85
 [32] richardson_lucy_iterative(measured::Matrix{Float32}, psf::Matrix{Float32})
    @ DeconvOptim ~/.julia/packages/DeconvOptim/d984T/src/lucy_richardson.jl:33
 [33] top-level scope
    @ REPL[219]:1
 [34] top-level scope
    @ /usr/share/julia/stdlib/v1.9/REPL/src/REPL.jl:1416

My psf at the center is 1

In [220]: psf(tf, 512, Δxy)[0,0]
> 1.0 + 0.0im

If I normalize it, I get an output

In [221]: DeconvOptim.richardson_lucy_iterative(Float32.(lr_imgs[:,:,1]), Float32.(ifftshift(abs.(no_offset_view(psf(tf, 512, Δxy)./sum(psf(tf,512, Δxy)))))))
> 512×512 Matrix{Float32}:
 0.000722969  0.00247124   0.00251804  0.00167182   0.00137161  0.00153345   …  0.000394556  0.000283565  0.000152645  0.000117615  0.00036432   0.00260431
 0.00269969   0.00965069   0.00960698  0.00594005   0.00451054  0.00483012      0.00166718   0.00118713   0.000564906  0.000323449  0.00081945   0.00644365
 0.00323262   0.0119314    0.0120749   0.00726379   0.00505583  0.00487872      0.00262322   0.00195198   0.000861986  0.000330069  0.000406617  0.00166461
 0.00229654   0.00897117   0.00995813  0.00633355   0.00420427  0.00356565      0.00257252   0.00208089   0.000900965  0.000252979  0.00014331   0.000279729
 0.00164897   0.00707362   0.00897682  0.00627268   0.0041557   0.00323875      0.00240401   0.00211108   0.000897628  0.000204787  6.58758f-5   8.38054f-5
 0.00154892   0.00717038   0.00976567  0.00701248   0.00462039  0.00359017   …  0.00257631   0.00237869   0.000971294  0.000193511  4.60927f-5   5.01164f-5
 0.00177654   0.00826314   0.0107331   0.00710855   0.00449392  0.0036923       0.0030169    0.00284355   0.00111592   0.000207322  4.54208f-5   5.16117f-5
 0.00211289   0.00929743   0.0107021   0.00619903   0.00372746  0.0033438       0.00351798   0.00333647   0.00128357   0.000234644  5.41294f-5   7.50505f-5
 ⋮                                                              ⋮            ⋱                                                      ⋮            
 0.000949476  0.00335143   0.00343603  0.00207084   0.0014692   0.00158476      0.00203276   0.0018451    0.00114504   0.000671677  0.000952491  0.00360017
 0.000990412  0.00362271   0.00388912  0.00239712   0.00167872  0.00177523   …  0.00196416   0.00173604   0.00107738   0.000652941  0.000965877  0.00373959
 0.000943852  0.00377063   0.00468583  0.00327644   0.00236232  0.00231858      0.00192046   0.00169956   0.0010891    0.00068204   0.00102001   0.00418354
 0.000652762  0.00278695   0.00407929  0.00339178   0.00262334  0.0023927       0.00160594   0.00147812   0.0010184    0.000649483  0.00087011   0.00327407
 0.000320699  0.00130202   0.00200308  0.00189792   0.0016297   0.00148447      0.00104648   0.0010486    0.000801955  0.000502798  0.000498155  0.00119067
 0.000259589  0.000868155  0.0010141   0.000889455  0.00085538  0.000899868     0.000846093  0.000948417  0.000763956  0.000422859  0.00028609   0.000394601
 0.000887702  0.00282082   0.0018205   0.00106839   0.00114415  0.00166391   …  0.001758     0.00221845   0.00162947   0.000634825  0.000304185  0.000385077
 0.00575295   0.0243073    0.0084146   0.0030446    0.0039014   0.00849665      0.00735271   0.0112588    0.0071264    0.00151111   0.000507424  1.48351f-6

kunzaatko · 2023-07-18T15:10:16Z

I guess that this may be an issue in zygote. I am not going to engage in debugging at the moment :) Should I post an issue?

roflmaostc · 2023-07-18T15:10:25Z

is your image positive?

kunzaatko · 2023-07-18T15:11:45Z

Yes.

In [222]: lr_imgs[:,:,1] |> extrema .|> float
> (0.00067528576f0, 0.021737328f0)

kunzaatko · 2023-07-18T15:12:04Z

Barely, though!

kunzaatko · 2023-07-18T15:13:13Z

Normalized it and the issue is still here

In [224]: DeconvOptim.richardson_lucy_iterative(imadjustintensity(Float32.(lr_imgs[:,:,1])), Float32.(ifftshift(abs.(no_offset_view(psf(tf, 512, Δxy))))))
ERROR: DomainError with -0.5827184:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Stacktrace:
  [1] throw_complex_domainerror(f::Symbol, x::Float32)
    @ Base.Math ./math.jl:33
  [2] sqrt
    @ ./math.jl:677 [inlined]
  [3] sqrt
    @ ~/.julia/packages/ForwardDiff/vXysl/src/dual.jl:240 [inlined]
  [4] #1384
    @ ~/.julia/packages/Zygote/JeHtr/src/lib/broadcast.jl:276 [inlined]
  [5] _broadcast_getindex_evalf
    @ ./broadcast.jl:683 [inlined]
  [6] _broadcast_getindex
    @ ./broadcast.jl:656 [inlined]
  [7] getindex
    @ ./broadcast.jl:610 [inlined]
  [8] macro expansion
    @ ./broadcast.jl:974 [inlined]
  [9] macro expansion
    @ ./simdloop.jl:77 [inlined]
 [10] copyto!
    @ ./broadcast.jl:973 [inlined]
 [11] copyto!
    @ ./broadcast.jl:926 [inlined]
 [12] copy
    @ ./broadcast.jl:898 [inlined]
 [13] materialize
    @ ./broadcast.jl:873 [inlined]
 [14] broadcast_forward
    @ ~/.julia/packages/Zygote/JeHtr/src/lib/broadcast.jl:282 [inlined]
 [15] _broadcast_generic
    @ ~/.julia/packages/Zygote/JeHtr/src/lib/broadcast.jl:212 [inlined]
 [16] adjoint
    @ ~/.julia/packages/Zygote/JeHtr/src/lib/broadcast.jl:205 [inlined]
 [17] _pullback(__context__::Zygote.Context{false}, 702::typeof(Base.Broadcast.broadcasted), 703::Base.Broadcast.DefaultArrayStyle{2}, f::typeof(sqrt), args::Matrix{Fl
oat32})
    @ Zygote ~/.julia/packages/ZygoteRules/OgCVT/src/adjoint.jl:66
 [18] _apply(::Function, ::Vararg{Any})
    @ Core ./boot.jl:838
 [19] adjoint
    @ ~/.julia/packages/Zygote/JeHtr/src/lib/lib.jl:203 [inlined]
 [20] _pullback
    @ ~/.julia/packages/ZygoteRules/OgCVT/src/adjoint.jl:66 [inlined]
 [21] _pullback
    @ ./broadcast.jl:1311 [inlined]
 [22] _pullback
    @ ~/.julia/packages/DeconvOptim/d984T/src/regularizer.jl:249 [inlined]
 [23] _pullback(ctx::Zygote.Context{false}, f::DeconvOptim.var"#615#616", args::Matrix{Float32})
    @ Zygote ~/.julia/packages/Zygote/JeHtr/src/compiler/interface2.jl:0
 [24] pullback(f::Function, cx::Zygote.Context{false}, args::Matrix{Float32})
    @ Zygote ~/.julia/packages/Zygote/JeHtr/src/compiler/interface.jl:44
 [25] pullback
    @ ~/.julia/packages/Zygote/JeHtr/src/compiler/interface.jl:42 [inlined]
 [26] gradient(f::Function, args::Matrix{Float32})
    @ Zygote ~/.julia/packages/Zygote/JeHtr/src/compiler/interface.jl:96
 [27] #invokelatest#2
    @ ./essentials.jl:816 [inlined]
 [28] invokelatest
    @ ./essentials.jl:813 [inlined]
 [29] (::DeconvOptim.var"#∇reg#149"{DeconvOptim.var"#615#616", Matrix{Float32}})(x::Matrix{Float32})
    @ DeconvOptim ~/.julia/packages/DeconvOptim/d984T/src/lucy_richardson.jl:57
 [30] (::DeconvOptim.var"#iter_with_reg#151"{Float64, DeconvOptim.var"#iter_without_reg#150"{Matrix{Float32}, Matrix{Float32}, Matrix{ComplexF32}, DeconvOptim.var"#co
v#137"{Matrix{ComplexF32}, Matrix{Float32}, AbstractFFTs.ScaledPlan{ComplexF32, FFTW.rFFTWPlan{ComplexF32, 1, false, 2, UnitRange{Int64}}, Float32}, Matrix{ComplexF32
, FFTW.rFFTWPlan{Float32, -1, false, 2, UnitRange{Int64}}}, Matrix{ComplexF32}}, Matrix{Float32}, DeconvOptim.var"#∇reg#149"{DeconvOptim.var"#615#616", Matrix{Float32
}})(rec::Matrix{Float32})
    @ DeconvOptim ~/.julia/packages/DeconvOptim/d984T/src/lucy_richardson.jl:65
 [31] richardson_lucy_iterative(measured::Matrix{Float32}, psf::Matrix{Float32}; regularizer::Function, λ::Float64, iterations::Int64, conv_dims::UnitRange{Int64}, pr
gress::Nothing)
    @ DeconvOptim ~/.julia/packages/DeconvOptim/d984T/src/lucy_richardson.jl:85
 [32] richardson_lucy_iterative(measured::Matrix{Float32}, psf::Matrix{Float32})
    @ DeconvOptim ~/.julia/packages/DeconvOptim/d984T/src/lucy_richardson.jl:33
 [33] top-level scope
    @ REPL[224]:1
 [34] top-level scope
    @ /usr/share/julia/stdlib/v1.9/REPL/src/REPL.jl:1416

roflmaostc · 2023-07-18T15:14:25Z

Can you share the data?

Seems like the errors comes from the regularizer?

 [29] (::DeconvOptim.var"#∇reg#149"{DeconvOptim.var"#615#616", Matrix{Float32}})(x::Matrix{Float32})
    @ DeconvOptim ~/.julia/packages/DeconvOptim/d984T/src/lucy_richardson.jl:57

https://github.com/roflmaostc/DeconvOptim.jl/blob/67ffbb8792fcdf96c409c3ee7dab12f5b1e4cbf9/src/lucy_richardson.jl#L57C4-L57C76

kunzaatko · 2023-07-18T15:19:05Z

image.zip

Sorry for the zip. It seems that github only supports compressed image formats. Actually our data is a video with about 600Mb. This is the first frame.

kunzaatko · 2023-07-18T15:20:03Z

I am using a transfer function generated with my package as:

In [230]: tf = IdealOTFwithCurvature(488u"nm", 1.4, 1.0, 0.3)

note: it is a stupid name for a type but this is due to the package being a work in progress

roflmaostc reviewed Jul 18, 2023

View reviewed changes

fix: add notes to documentation about shift

707a6e5

[Ticket: #44]

roflmaostc merged commit 67ffbb8 into roflmaostc:master Jul 18, 2023
4 checks passed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

fix: add notes to documentation about shift #45

fix: add notes to documentation about shift #45

kunzaatko commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

roflmaostc Jul 18, 2023

kunzaatko commented Jul 18, 2023

roflmaostc commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

roflmaostc commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

roflmaostc commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

kunzaatko commented Jul 18, 2023 •

edited

fix: add notes to documentation about shift #45

fix: add notes to documentation about shift #45

Conversation

kunzaatko commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

roflmaostc Jul 18, 2023

Choose a reason for hiding this comment

kunzaatko commented Jul 18, 2023

roflmaostc commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

roflmaostc commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

roflmaostc commented Jul 18, 2023

kunzaatko commented Jul 18, 2023

kunzaatko commented Jul 18, 2023 • edited

kunzaatko commented Jul 18, 2023 •

edited