tmporarily added some details for gaussian fitting.

Author: RainerHeintzmann

examples/Project.toml

@@ -2,6 +2,7 @@
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 ComponentArrays = "b0b7db55-cfe3-40fc-9ded-d10e2dbeff66"
 InverseModeling = "ce844058-9528-415d-a63d-06f3dd08b29f"
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 Noise = "81d43f40-5267-43b7-ae1c-8b967f377efa"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 SeparableFunctions = "c8c7ead4-852c-491e-a42d-3d43bc74259e"

examples/gauss_fit.jl

@@ -6,66 +6,110 @@ using Noise
 using CUDA
 
 # simulate a gaussian blob with Poisson noise and fit it with a Gaussian function
-sz = (7,7) # (1600, 1600)
-many_fits = true
+sz = (9, 9) # (1600, 1600)
+many_off = true
+many_int = true
+many_bg = true
+many_sig = true
+
 use_cuda = false
-N = 1_000
-hyperplanes = many_fits ? rand(Float32, (1, N)) : 0
-hp_zeros = many_fits ? zeros(Float32, (1, N)) : 0
+DType = Float32
+N = 10_000
+hp_off = many_off ? 2 .*rand(DType, (1, N)) : 0
+hp_sig = many_sig ? zeros(DType, (1, N)) : 0
+hp_int = many_int ? 1 .+ rand(DType, (1, N)) : 1
 
-off = [3.2f0, 3.5f0] .+ hyperplanes
-sigma = [1.4f0, 1.1f0] .+ hp_zeros
-hyperplanes = many_fits ? rand(Float32, (1, N)) : 0
-intensity = [50f0] .* (1 .+ hyperplanes)
-vec_true = ComponentVector(;bg=10.0f0 .+ hp_zeros, intensity=intensity, off = off, args = sigma)
-vec_true = Float64.(vec_true)
+off = [5.2, 4.5] .+ hp_off
+sigma = [1.4, 1.1] .+ hp_sig
+intensity = [50] .* hp_int
+bg = many_bg ? 10.2 .+ zeros(DType, (1, N)) : 10.2
+vec_true = DType.(ComponentVector(;bg=bg, intensity=intensity, off = off, args = sigma))
 
+# create the perfect spots:
 pdat = gaussian_vec(sz, vec_true)
-dat = Float32.(poisson(Float64.(pdat)))
+dat = DType.(poisson(Float64.(pdat)))
+
+startvals = DType.(ComponentVector(gauss_start(dat, 0.2, length(sz))));
+sum(abs2.(collect(startvals.off) .- vec_true.off))
 
 pdat = (use_cuda) ? CuArray(pdat) : pdat
 dat = (use_cuda) ? CuArray(dat) : dat
-
-qdat = copy(pdat)
-qdat .= qdat[:,:,1]
-# now prepare the fitting:
-# myfg! = get_fg!(pdat, gaussian_raw, length(sz); loss=loss_anscombe_pos, bg=7f0);
-myfg! = get_fg!(qdat, gaussian_raw, length(sz); loss=loss_gaussian);
-shyperplanes = many_fits ? zeros(Float32, (1, size(dat)[end])) : 0
-soff = [4.0f0, 4.0f0] .+ shyperplanes
-bg = [0.5f0] .+ shyperplanes
-intensity = [45f0] .+ shyperplanes
-sigma = [3.0f0, 2.0f0] .+ shyperplanes
 if (use_cuda)
-    bg = CuArray([bg])
-    intensity = CuArray(intensity)
-    soff = CuArray(soff)
-    sigma = CuArray(sigma)
+    startvals = ComponentVector(;bg=CuArray(startvals.bg), intensity=CuArray(startvals.intensity), off = CuArray(startvals.off), args = CuArray(startvals.args))
 end
-startvals = ComponentVector(;bg=bg, intensity=intensity, off = soff, args = sigma)
-startvals = Float64.(startvals)
-opt = Optim.Options(iterations = 50); #
-odo = OnceDifferentiable(Optim.NLSolversBase.only_fg!(myfg!), startvals);
+
+# @vt pdat gaussian_vec(sz, startvals)
+# qdat = copy(pdat)
+# qdat .= qdat[:,:,1]
+# now prepare the fitting:
+# myfg! = get_fg!(dat, gaussian_raw, length(sz); loss=loss_anscombe_pos, bg=0.1f0);
+# myfg! = get_fg!(dat, gaussian_raw, length(sz); loss=loss_gaussian);
+myfg! = get_fg!(dat, gaussian_raw, length(sz); loss=loss_poisson_pos);
+# hp_off2 = many_off ? zeros(DType, (1, size(dat)[end])) : 0
+# soff = [4.0, 4.0] .+ hp_off2
+# bg = [0.5] .+ shyperplanes
+# intensity = [45.0] .+ hp_int
+# sigma = [3.0, 2.0] .+ hp_sig
+# if (use_cuda)
+#     bg = CuArray([bg])
+#     intensity = CuArray(intensity)
+#     soff = CuArray(soff)
+#     sigma = CuArray(sigma)
+# end
+# startvals = DType.(ComponentVector(;bg=bg, intensity=intensity, off = soff, args = sigma))
+opt = Optim.Options(iterations = 1500); #
 
 if (false)
     G = copy(startvals)
     myfg!(1, G, startvals)
 
-    myfg2! = get_fg!(pdat[:,:,1], gaussian_raw, length(sz); loss=loss_anscombe_pos, bg=7f0);
+    myfg2! = get_fg!(pdat[:,:,1], gaussian_raw, length(sz); loss=loss_anscombe_pos, bg=0.1f0);
     sv = ComponentVector{Float32}(bg=startvals.bg[1], intensity=startvals.intensity[1], off = startvals.off[:,1], args = startvals.args[:,1])
     G2 = copy(sv)
     myfg2!(1, G2, sv)
     G2
 end
 
 # and perform the fit
-@time reso = Optim.optimize(odo, startvals, Optim.LBFGS(), opt);
+svb = copy(startvals)
+# svb.args = svb.args .* 1.2f0
+odo = OnceDifferentiable(Optim.NLSolversBase.only_fg!(myfg!), svb);
+@time reso = Optim.optimize(odo, svb, Optim.LBFGS(), opt);
+# 14 sec, CUDA: 2.5 sec
 # 2 sec, 5k fits/s (44.25 k allocations: 1.546 GiB, 7.35% gc time)
 # with intensity variations: 26.833106 seconds (532.47 k allocations: 20.251 GiB, 5.99% gc time)
 # in Cuda: 
 reso.f_calls # 61   # 1766 für 10_000 fits, 155, 2.2 sec for 10_000 fits with all entries being vectors
 reso.minimum # 
-@vt pdat gaussian_vec(sz, startvals) gaussian_vec(sz, reso.minimizer)
+@vt pdat dat (gaussian_vec(sz, startvals).-dat) (gaussian_vec(sz, reso.minimizer).-dat)
+
+success = sum(abs.(collect(startvals.off) .- vec_true.off), dims=1) .< 0.5
+success = success .&& sum(abs.(collect(reso.minimizer.off) .- vec_true.off), dims=1) .< 0.5
+sum(.!success)
+# findfirst(.!success)
+@vt collect(dat)[:,:,.!success[:]] collect(pdat)[:,:,.!success[:]] collect(gaussian_vec(sz, startvals))[:,:,.!success[:]] collect(gaussian_vec(sz, reso.minimizer))[:,:,.!success[:]]
+
+ff = findfirst(.!success)[2]
+sv = ComponentVector(off=startvals.off[:,ff], bg = startvals.bg[:,ff], intensity=startvals.intensity[:,ff], args = startvals.args[:,ff])
+sdat = dat[:,:,ff]
+afg! = get_fg!(sdat, gaussian_raw, length(sz); loss=loss_gaussian); # loss_poisson_pos
+odo = OnceDifferentiable(Optim.NLSolversBase.only_fg!(afg!), sv);
+@time reso = Optim.optimize(odo, sv, Optim.LBFGS(), opt);
+@vt sdat gaussian_vec(sz, sv) gaussian_vec(sz, reso.minimizer)
+
+trueoff = vec_true.off[:,collect(success[:])]
+startoff = collect(startvals.off[:,success[:]])
+minoff = collect(reso.minimizer.off[:,success[:]])
+sum(abs2.(startoff .- trueoff), dims=2) # 0 failed, 68.29712, 42.612137
+sum(abs2.(minoff .- trueoff), dims=2) # 21 failed, 57.92457, 33.990776 # Anscombe: 5 failed, 58.806953, 34.21795
+@vt abs.(startoff .- trueoff) abs.(minoff .- trueoff)
+
+sum(abs2.(collect(startvals.intensity) .- vec_true.intensity))
+sum(abs2.(collect(reso.minimizer.intensity) .- vec_true.intensity))
+sum(abs2.(collect(startvals.bg) .- vec_true.bg))
+sum(abs2.(collect(reso.minimizer.bg) .- vec_true.bg))
+sum(abs2.(collect(startvals.args) .- vec_true.args), dims=2)
+sum(abs2.(collect(reso.minimizer.args) .- vec_true.args), dims=2)
 
 odo = OnceDifferentiable(Optim.NLSolversBase.only_fg!(myfg!), startvals);
 if isa(dat, CuArray)
@@ -97,3 +141,32 @@ res1
 # 5 ms (39192 allocations: 4.35 MiB)
 
 # @btime Optim.optimize($loss, $off_start, $sigma_start, LBFGS(); autodiff = :forward); # 1.000 ms (10001 allocations: 1.53 MiB)
+
+# here is the LsqFit example from https://github.com/JuliaNLSolvers/LsqFit.jl
+using LsqFit
+
+x = collect(range(0, stop=200, length=201))
+y = collect(range(0, stop=200, length=201))
+
+xy = hcat(x, y)
+
+function twoD_Gaussian(xy, p)
+    amplitude, xo, yo, sigma_x, sigma_y, theta, offset = p
+    a = (cos(theta)^2)/(2*sigma_x^2) + (sin(theta)^2)/(2*sigma_y^2)
+    b = -(sin(2*theta))/(4*sigma_x^2) + (sin(2*theta))/(4*sigma_y^2)
+    c = (sin(theta)^2)/(2*sigma_x^2) + (cos(theta)^2)/(2*sigma_y^2)
+
+    # creating linear meshgrid from xy
+    x = xy[:, 1]
+    y = xy[:, 2]
+    g = offset .+ amplitude .* exp.( - (a.*((x .- xo).^2) + 2 .* b .* (x .- xo) .* (y .- yo) + c * ((y .- yo).^2)))
+    return g[:]
+end
+
+p0 = Float64.([3, 100, 100, 20, 40, 0, 10])
+data = twoD_Gaussian(xy, p0)
+
+# Noisy data
+data_noisy = data + 0.2 * randn(size(data))
+
+fit = LsqFit.curve_fit(twoD_Gaussian, xy, data_noisy, p0)

src/SeparableFunctions.jl

@@ -47,13 +47,15 @@ export kwargs_to_args
 
 export loss_anscombe, loss_anscombe_pos, loss_gaussian, loss_poisson, loss_poisson_pos 
 export get_vec_dim
+export gauss_start
 
 DefaultResElType = Float32
 DefaultArrType = Array{DefaultResElType}
 DefaultComplexArrType = Array{complex(DefaultResElType)}
 
 include("utilities.jl")
 include("losses.jl")
+include("gauss_params.jl")
 include("specific.jl")
 include("general.jl")
 include("radial.jl")

src/gauss_params.jl

@@ -0,0 +1,89 @@
+"""
+    get_com(dat::AbstractArray{T, N}, mask, prod_dims=N) where {T, N}
+
+returns the center of mass of a data array `dat` with a mask `mask`. 
+"""
+function get_com!(dst::AbstractMatrix{T}, dat::AbstractArray{T, N}, mask, prod_dims=N) where {T, N}
+    ax = axes(dat)[1:prod_dims]
+    ax = ntuple((d) -> reorient(ax[d], Val(d), Val(prod_dims)), Val(prod_dims))
+    sum_dat_mask = sum(dat.*mask, dims=1:prod_dims)
+    for d in 1:prod_dims
+        dst[d,:] .= (sum(ax[d].*dat.*mask, dims=1:prod_dims) ./ sum_dat_mask)[:]
+    end
+end
+function get_com(dat::AbstractArray{T, N}, mask, prod_dims=N) where {T, N}
+    dst = similar(dat, (prod_dims, size(dat)[prod_dims+1:end]...))
+    get_com!(dst, dat, mask, prod_dims)
+    return dst
+end
+
+"""
+    get_std(dat::AbstractArray{T, N}, t_ctr, mask) where {T, N}
+
+returns the center of mass of a data array `dat` with a mask `mask`. 
+"""
+function get_std!(dst::AbstractMatrix{T}, dat::AbstractArray{T, N}, mask, t_ctr, prod_dims=N) where {T, N}
+    sum_dat_mask = sum(dat.*mask, dims=1:prod_dims)
+    for d in 1:prod_dims
+        ax = reorient(axes(dat)[d], Val(d), Val(prod_dims))
+        ctrs = reshape(t_ctr[d,:], (ones(Integer, prod_dims)..., size(t_ctr)[2:end]...))
+        dst[d,:] .= (sqrt.(sum(dat.*mask.*(ax .- ctrs).^2, dims=1:prod_dims) ./ sum_dat_mask))[:]
+    end
+end
+function get_std(dat::AbstractArray{T, N}, mask, t_ctr, prod_dims=N) where {T, N}
+    dst = similar(dat, (prod_dims, size(dat)[prod_dims+1:end]...))
+    get_std!(dst, dat, mask, t_ctr, prod_dims)
+    return dst
+end
+
+function get_bg_sum(dat::AbstractArray{T, N}, prod_dims=N) where {T, N}
+    sz = size(dat)[1:prod_dims]
+    isz = max.(1, sz.-2)
+    inner = select_region_view(dat, isz)
+    sdat = sum(dat, dims=1:prod_dims)
+    sinner = sum(inner, dims=1:prod_dims)
+    ninner = prod(isz)
+    nedge = prod(sz) - ninner
+    bg = (sdat .- sinner)/nedge
+    return bg, (sinner .- ninner.*bg) # minimum(dat)
+end
+
+function get_intensity(dat::AbstractArray{T, N}, prod_dims=N) where {T, N}
+    sz = size(dat)[1:prod_dims]
+    # select a 3x3 pixel region and calculate the average intensity in it
+    psz = min.(sz, 3)
+    peak = select_region_view(dat, psz)
+    npeak = prod(psz)
+    return sum(peak, dims=1:prod_dims) / npeak # maximum(meas) - offset
+end
+
+function gauss_start(meas::AbstractArray{T, N}, rel_thresh=0.2, prod_dims=N; has_covariance=false) where {T, N}
+    sz = size(meas)[1:prod_dims]
+    bg, sinner = get_bg_sum(meas, prod_dims)
+    # just a fletch factor 1.25 for typical spotsize
+    i0 = get_intensity(meas, prod_dims) * 1.25
+    meas = meas .- bg .- i0 .* rel_thresh
+    mymask = meas .> 0
+
+    t_ctr = get_com(meas, mymask, prod_dims)
+    σ =  get_std(meas, mymask, t_ctr, prod_dims) .* 1.22 ./ (1-rel_thresh)
+    μ =  t_ctr #  .- (sz.÷2 .+1)
+    # @show μ
+    # @show σ
+    # sumpix = sum(meas .* mymask)
+    # tosum(apos, dat) = apos .* dat
+    # pos = idx(sz)
+    # μ = tuple_sum(tosum.(pos, meas.*mymask)) ./ sumpix
+    # mysqr(apos, dat) = abs2.(apos) .* dat 
+    # pos = idx(size(meas), offset=size(pos).÷2 .+1 .+ μ)
+    # σ = 1.0 .* max.(1.0, sqrt.(tuple_sum(mysqr.(pos, meas.*mymask )) ./ sumpix))
+    # if has_covariance
+    #     σ = [σ..., zeros(((length(μ))*(length(μ)-1))÷2)...]
+    #     @show σ
+    # end
+    maxint = reshape(0.178 .* sinner[:] ./ prod(σ, dims=1)[:], (1, prod(size(sinner))))
+    bg = reshape(bg[:], size(maxint))
+    start_params = (bg=bg, intensity = maxint, off=μ, args=σ)
+    # @show start_params
+    return start_params # Fixed(), Positive
+end