I am optimizing cross correlation of two arrays taken over there transfer function intersection. My function looks like this:
otf_lr_1 = otf(tf, fft_lr_1, Δxy) # does not change with Δϕ
fxs = fftfreq(size(fft_lr_2, 1), 1 / Δxy) .* ones(size(fft_lr_2, 2))'
fys = ones(size(fft_lr_2, 1)) .* fftfreq(size(fft_lr_2, 2), 1 / Δxy)'
function xcorr(Δϕ)
otf_lr_2 = otf.(tf,
fxs .+ (Δϕ[1] / (Δxy * (size(fft_lr_2, 1)))),
fys .+ (Δϕ[2] / (Δxy * (size(fft_lr_2, 2))))
) # non-shifted
fft_lr_2_sΔϕ = fft(
ifft(fft_lr_2) .* exp.(-2im * π * (Δϕ[1] * fftshift(fxs * Δxy) + Δϕ[2] * fftshift(fys * Δxy)))
) # shifted by Δϕ (fourier shift theorem)
Δϕr = hypot.(fxs .* Δxy, fys .* Δxy) ./ splat(hypot)(Δϕ ./ size(fft_lr_2_sΔϕ))
intersection = dropdims(all(isone, cat(
otf_lr_1 .> alg.MTF_min,
otf_lr_2 .> alg.MTF_min,
Δϕr .> alg.Δϕr_min,
Δϕr .< 1 - alg.Δϕr_min;
dims=3); dims=3); dims=3)
fft_lr_1_i = zeros(eltype(fft_lr_1), size(fft_lr_1)) # at intersection
fft_lr_1_i[intersection] = fft_lr_1[intersection] ./ otf_lr_1[intersection] # transfer normalized
fft_lr_2_i = zeros(eltype(fft_lr_2_sΔϕ), size(fft_lr_2_sΔϕ)) # at intersection
fft_lr_2_i[intersection] = fft_lr_2_sΔϕ[intersection] ./ otf_lr_1[intersection] # same OTF... Imitating the image formation where `lr_img_1` is at the center
lr_1_i = ifft(fft_lr_1_i) # spatial from intersection
lr_2_i = ifft(fft_lr_2_i)
return -abs(sum(lr_2_i .* conj(lr_1_i)) / sum(abs.(lr_1_i) .^ 2))
end
When using Optim.jl and after waiting for about 10 minutes, the optimization did not finish running.
res = Optim.optimize(xcorr,
[(alg.Δϕ_0 .- box_size)...],
[(alg.Δϕ_0 .+ box_size)...],
[float.(alg.Δϕ_0)...], Fminbox(LBFGS(; m=4))
)
if Optim.converged(res)
return Optim.minimizer(res)
else
@error "Optimization did not converge"
end
Then after changing to scipy, I got the result (correct results) in about 10 seconds with
using PyCall
optimize = pyimport("scipy.optimize")
res = optimize.minimize(
xcorr,
alg.Δϕ_0;
bounds=((alg.Δϕ_0[1] .+ (-box_size, +box_size)),
(alg.Δϕ_0[2] .+ (-box_size, +box_size)))
)
if res["success"]
return res["x"] |> Tuple
else
@error "Optimization did not converge: $(res["message"])"
end
It seems to me that the parameters of the optimization is the same in both cases.
Here are the packages that I am using and the serialized data that is input
using FFTW
using Unitful
using TransferFunctions
using OffsetArrays
using OffsetArrays: centered
using TransferFunctions
using TransferFunctions: otf_support
using Serialization
using Optim
fft_lr_1 = deserialize("file_one.ser")
fft_lr_2 = deserialize("file_two.ser")
struct Dummy
Δϕ_0
Δϕr_min
MTF_min
end
alg = Dummy((-176,-24), 0.15, 0.05)
tf = IdealOTFwithCurvature(488u"nm", 1.4, 1.0, 0.3)
Δxy = 61u"nm"
box_size = Inf
otf_lr_1 = otf(tf, fft_lr_1, Δxy) # does not change with Δϕ
fxs = fftfreq(size(fft_lr_2, 1), 1 / Δxy) .* ones(size(fft_lr_2, 2))'
fys = ones(size(fft_lr_2, 1)) .* fftfreq(size(fft_lr_2, 2), 1 / Δxy)'
The TransferFunctions package is of my creation and is here on GH.
I hope I provided everything that is necessary. Thank you very much for helping me untangle this!