diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000..9016df2 Binary files /dev/null and b/.DS_Store differ diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml new file mode 100644 index 0000000..3a67587 --- /dev/null +++ b/.github/workflows/ci.yml @@ -0,0 +1,64 @@ +name: CI +on: + pull_request: + branches: + - main + - develop + push: + branches: + - main + - develop + tags: '*' +concurrency: + # Skip intermediate builds: always. + # Cancel intermediate builds: only if it is a pull request build. + group: ${{ github.workflow }}-${{ github.ref }} + cancel-in-progress: ${{ startsWith(github.ref, 'refs/pull/') }} +jobs: + test: + name: Julia ${{ matrix.version }} - ${{ matrix.os }} - ${{ matrix.arch }} - ${{ github.event_name }} + runs-on: ${{ matrix.os }} + strategy: + fail-fast: false + matrix: + version: + - '1.9' + os: + - [ubuntu-latest] + arch: + - x64 + steps: + - uses: actions/checkout@v3 + - uses: julia-actions/setup-julia@v1 + with: + version: ${{ matrix.version }} + arch: ${{ matrix.arch }} + - uses: julia-actions/cache@v1 + - name: Add unregistered dependencies + run: julia -e 'using Pkg; Pkg.develop(PackageSpec(url="https://github.com/rainerheintzmann/EvalMultiPoly.jl"))' + + - uses: julia-actions/julia-buildpkg@v1 + - uses: julia-actions/julia-runtest@v1 + - uses: julia-actions/julia-processcoverage@v1 + - uses: codecov/codecov-action@v4 + env: + CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }} + with: + file: lcov.info + + + docs: + name: Documentation + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v2 + - uses: julia-actions/setup-julia@v1 + with: + version: "1.9" + - name: Add unregistered dependencies + run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(url="https://github.com/rainerheintzmann/EvalMultiPoly.jl")); Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()' + + - uses: julia-actions/julia-docdeploy@releases/v1 + env: + GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }} + DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }} \ No newline at end of file diff --git a/.gitignore b/.gitignore index 29126e4..32f88c7 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,19 @@ docs/site/ # committed for packages, but should be committed for applications that require a static # environment. Manifest.toml +examples/Manifest.toml +examples/*.tiff +examples/image_registration_2dpolim.jl + +*.mp4 +*.png +examples/figures/* +examples/*.tif +examples/*.ipynb +*.jpeg +*.gif +*.jpg +*.txt +*.txt +.*tiff +examples/6-4-Rigis.jl diff --git a/Project.toml b/Project.toml index 7c59800..d7cc818 100644 --- a/Project.toml +++ b/Project.toml @@ -4,5 +4,25 @@ authors = ["RainerHeintzmann "] version = "0.1.0" [deps] +EvalMultiPoly = "c78649ec-f9ed-405e-be6d-0472f43586aa" FourierTools = "b18b359b-aebc-45ac-a139-9c0ccbb2871e" Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +Requires = "ae029012-a4dd-5104-9daa-d747884805df" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" + +[compat] +FourierTools = "0.4" +Interpolations = "0.15.1, 0.16" +StaticArrays = "1.9.6, 1.10" +Aqua = "0.8" +Test = "1.9.3" +Zygote = "0.6.70" +julia = "1" + +[extras] +Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595" +Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" + +[targets] +test = ["Aqua", "Test", "Zygote"] diff --git a/README.md b/README.md index 8f912d3..7f863c3 100644 --- a/README.md +++ b/README.md @@ -1,2 +1,6 @@ # DataToFunctions.jl -Represents (measured) data as a function, which supports scaling and shifting. It is intended to be used as a tool for fitting data, where the fitting function is given by measured data. +Represents (measured) data as a function, which supports affine and generally, matrix transformations for arrays. It is intended to be used as a tool for fitting data, where the fitting function is given by measured data. + +It assumes the data as an Interpolation object (function), then it can transform the data to any arbitrary affine or matrix transformations. + +What is important about this package is that it does not assume any pre-defined function as the data for the inverse modeling procedures. \ No newline at end of file diff --git a/docs/Project.toml b/docs/Project.toml new file mode 100644 index 0000000..d2cdc05 --- /dev/null +++ b/docs/Project.toml @@ -0,0 +1,11 @@ +[deps] +DataToFunctions = "64cfdffa-4d02-49ee-ae8b-a805370874f5" +Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +DocumenterTools = "35a29f4d-8980-5a13-9543-d66fff28ecb8" +EvalMultiPoly = "c78649ec-f9ed-405e-be6d-0472f43586aa" +GR = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71" +Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +SyntheticObjects = "e7028c27-0967-45e9-8fdb-dbc10ccb2b0a" diff --git a/docs/make.jl b/docs/make.jl new file mode 100644 index 0000000..f3f5377 --- /dev/null +++ b/docs/make.jl @@ -0,0 +1,23 @@ +import Pkg +Pkg.activate(@__DIR__) + +cd(@__DIR__) # go into `docs` folder + +using Documenter, Literate, DataToFunctions + +# convert tutorial/examples to markdown +Literate.markdown("./src/tutorial.jl", "./src") +# Which markdown files to compile to HTML +# (which is also the sidebar and the table +# of contents for your documentation) + +pages = Any[ + "Introduction" => "index.md", + "Tutorial" => "tutorial.md", + "API" => "api.md", + ] + +# compile to HTML: +makedocs(; sitename="DataToFunctions.jl", pages, modules = [DataToFunctions], warnonly = true) + +deploydocs(repo = "github.com/RainerHeintzmann/DataToFunctions.jl.git", devbranch = "develop") \ No newline at end of file diff --git a/docs/src/api.md b/docs/src/api.md new file mode 100644 index 0000000..2502d57 --- /dev/null +++ b/docs/src/api.md @@ -0,0 +1,20 @@ +# API + +```@docs +get_function +add_dim +red_dim +red_dim_apply +idx_appy +mat_mul +func_transform +func_transform_tup +apply_transform +apply_transform_homogen +apply_transform_affine +get_function_tuple +get_function_svec +get_function_affine +get_function_poly +split_tuple +``` \ No newline at end of file diff --git a/docs/src/index.md b/docs/src/index.md new file mode 100644 index 0000000..fa5f238 --- /dev/null +++ b/docs/src/index.md @@ -0,0 +1,5 @@ +# DataTofunctions.jl + +```@docs +DataToFunctions +``` \ No newline at end of file diff --git a/docs/src/tutorial.jl b/docs/src/tutorial.jl new file mode 100644 index 0000000..2dac15b --- /dev/null +++ b/docs/src/tutorial.jl @@ -0,0 +1,50 @@ +# Tutorial for the DataToFunctions.jl package + +import Pkg +Pkg.add("SyntheticObjects") +Pkg.add("StaticArrays") +Pkg.add("Plots") +Pkg.add("Random") + +# Load the packages +using DataToFunctions +using SyntheticObjects +using StaticArrays +using Random +using Plots + +# We set the Random.seed for reproducibility +Random.seed!(14) + +# Define the data +sz = 128 +data = filaments3D((sz, sz, 1), num_filaments=3)[:, :, 1] + +# define the interpolation function using this package +f_affine = get_function_affine(data) + +# define the transformation parameters +params = [1.0, 2.0, 1.2, 0.8, 0.0, 0.0, 0.0] + +# apply the transformation +data_transformed = f_affine(params) +heatmap(data_transformed, aspect_ratio=1 + , title="Transformed data using a transformation parameters vector" + , titlefontsize=10, size=(500, 500) + , xlabel="X", ylabel="Y") + + +## The next example is to use a transformation matrix +matrix_c = [1.0 0.0 2.0; + 0.0 1.0 3.0; + 0.0 0.0 1.0] + +data_transformed_m = f_affine(SMatrix{3, 3}(matrix_c)) +heatmap(data_transformed_m, aspect_ratio=1 + , title="Transformed data using a transformation matrix" + , titlefontsize=10, size=(500, 500) + , xlabel="X", ylabel="Y") + + +# Now we try to do the transformation using a polynomial function +# we first define the interpolation object using the DataToFunctions package with a polynomial of order 1 diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md new file mode 100644 index 0000000..f1d87d8 --- /dev/null +++ b/docs/src/tutorial.md @@ -0,0 +1,78 @@ +```@meta +EditURL = "tutorial.jl" +``` + +Tutorial for the DataToFunctions.jl package + +````@example tutorial +import Pkg +Pkg.add("SyntheticObjects") +Pkg.add("StaticArrays") +Pkg.add("Plots") +Pkg.add("Random") +```` + +Load the packages + +````@example tutorial +using DataToFunctions +using SyntheticObjects +using StaticArrays +using Random +using Plots +```` + +We set the Random.seed for reproducibility + +````@example tutorial +Random.seed!(14) +```` + +Define the data + +````@example tutorial +sz = 128 +data = filaments3D((sz, sz, 1), num_filaments=3)[:, :, 1] +```` + +define the interpolation function using this package + +````@example tutorial +f_affine = get_function_affine(data) +```` + +define the transformation parameters + +````@example tutorial +params = [1.0, 2.0, 1.2, 0.8, 0.0, 0.0, 0.0] +```` + +apply the transformation + +````@example tutorial +data_transformed = f_affine(params) +heatmap(data_transformed, aspect_ratio=1 + , title="Transformed data using a transformation parameters vector" + , titlefontsize=10, size=(500, 500) + , xlabel="X", ylabel="Y") + + +# The next example is to use a transformation matrix +matrix_c = [1.0 0.0 2.0; + 0.0 1.0 3.0; + 0.0 0.0 1.0] + +data_transformed_m = f_affine(SMatrix{3, 3}(matrix_c)) +heatmap(data_transformed_m, aspect_ratio=1 + , title="Transformed data using a transformation matrix" + , titlefontsize=10, size=(500, 500) + , xlabel="X", ylabel="Y") +```` + +Now we try to do the transformation using a polynomial function +we first define the interpolation object using the DataToFunctions package with a polynomial of order 1 + +--- + +*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).* + diff --git a/examples/PSF_fitting_general.jl b/examples/PSF_fitting_general.jl new file mode 100644 index 0000000..6bc5d8f --- /dev/null +++ b/examples/PSF_fitting_general.jl @@ -0,0 +1,412 @@ +using DataToFunctions +using Optim, StaticArrays, LinearAlgebra +using PointSpreadFunctions +using Zygote +using ForwardDiff, LineSearches, Plots, Printf +using View5D +using Distributions, Rotations +using Plots +using TestImages +using BenchmarkTools +#using InverseModeling +import Random +using Noise, Images, CSV, TiffImages +using ProgressBars + +""" + perform_fit_general(loss_function, fitting_data::AbstractArray) + +Performs a fit to the fitting data using a loss function defined by the user + +# Arguments +`loss_function`: User-defined loss function which is minimized +`fitting_data`: The data which is being fitted + +# Returns +a vector of 7 parameters: 2 for the shift, 2 for the scaling, 2 for shear, and 1 for rotation angle + +# Example +there is an example of this function in the `examples/star_fitting_genaral.jl` + +""" +function perform_fit_general(loss_function, fitting_data::AbstractArray, init_x::AbstractArray{T}) where T + # guess the shift parameters by taking the maximum values of the array and + # centering the positions + ##a, b = Tuple(argmax(fitting_data)) .- size(fitting_data) ./2.0 .- 1.0 + #print("INSIDE!!! hehe") + # assigning the initial parameter estimates + # init_x = vec([0.5, -1.5, 1.0, 1.0, 0.0, 0.0, pi/5]) #ndims(fitting_data)+1, ndims(fitting_data)+1)) + # reshape(Matrix(1.0*I, ndims(fitting_data)+1, ndims(fitting_data)+1), 1, 9)) #[a, b, 1.0, 1.0, 0.001, 0.001, 0.001] + + # setting the lower and upper boundary of the parameter values based on their limits + lower = T[-1*size(fitting_data)[1], -1*size(fitting_data)[2], 0.0, 0.0, -0.01, -0.01, 0.0] + upper = T[size(fitting_data)[1], size(fitting_data)[2], size(fitting_data)[1], size(fitting_data)[2], 0.01, 0.01, pi/2.0] + + # initializing the LBFGS optimizer + inner_optimizer = BFGS()#; m=3, linesearch=LineSearches.BackTracking(order=3)) + + # Computer, Optimize! :D + res = optimize( + loss_function, + + #LBFGS(), + lower, upper, + init_x, + Fminbox(inner_optimizer), + Optim.Options(store_trace = true, extended_trace = true, iterations=5000), + autodiff = :forward + ) + + # return the estimated parameters + return Optim.minimizer(res), res +end + +function perform_fit(loss_function, init_x::AbstractArray{T}) where T + + # Computer, Optimize! :D + res = optimize( + loss_function, + init_x, + #Newton(), + #BFGS(),#; linesearch=LineSearches.BackTracking(order=3)), + LBFGS(),#; linesearch=LineSearches.BackTracking(order=3)), + #lower, upper, + #init_x, + #Fminbox(inner_optimizer), + Optim.Options(store_trace = true, extended_trace = true, iterations=5000), + autodiff = :forward + ) + + # return the estimated parameters + return Optim.minimizer(res), res +end + +""" + apply_transform(;matrix=true, sz=64, dtype=Float32, noise_level=1/20.0) + +This function is designed for applying a transformation to a sample data using either matrix transformations or +parametric transformations based on the provided arguments. + +# Arguments +`matrix`: if true, the fitting is done using a matrix transformation, otherwise, the fitting is done using a parametric transformation +`sz`: the size of the sample data +`dtype`: the data type of the sample data +`noise_level`: the noise level to add to the sample data + +# Returns +a tuple of two arrays: the first array is the fitting data, and the second array is the estimated fitting data + +# Example +apply_transform(matrix=true, sz=64, dtype=Float32, noise_level=1/20.0) + +""" +function apply_transform(;matrix=false, sz=64, dtype=Float32, n_photons=1000, pure_rand=false, from_params=true, plotting=false) + Random.seed!(14) + + # defining the mean and the varixance of the test normal (Gaussian) distribution + μ = [0, 0] + Σ = [sz/10 0.0; + 0.0 sz/10] + + # initializing the multivariate normal distribution + p = MvNormal(μ, Σ) + + # to define the sample array based on a 2D normal distribution + X = -1*sz/2.0:1*sz/2.0 + Y = -1*sz/2.0:1*sz/2.0 + + z = [pdf(p, [x,y]) for y in Y, x in X] + + + # creating a PSF for the widefield microscope + sz_psf = (sz, sz, 100) + sampling = (0.040, 0.040, 0.050) + # simulate a confocal PSF + aberrations = Aberrations([Zernike_HorizontalComa],[0.8]); + pp = PSFParams(0.5,1.4,1.52, method=MethodPropagateIterative, aberrations=aberrations); + + #pp_ex = PSFParams(pp_em; λ=0.488);#, method=MethodPropagateIterative, aplanatic=aplanatic_illumination, aberrations=aberrations); + p_psf_3d = psf(sz_psf, pp, sampling=sampling); + p_psf = p_psf_3d[:, :, 50] + #sample_data = p_psf ./ maximum(p_psf) + + # normalizing the sample data + sample_data = p_psf ./ maximum(p_psf) + # sample_data = dtype.(z[1:sz, 1:sz]./maximum(z)) + # sample_data = dtype.(TestImages.shepp_logan(sz)) + # sample_data = rand(dtype, (sz, sz)) + + # sample_data .+= rand(dtype, (size(sample_data)...)).*noise_level; + p_img = n_photons .* (sample_data);# ./ maximum(sample_data)) + n_img = dtype.(poisson(Float64.(p_img))) + + x_cen, y_cen = (size(n_img) ./ 2.0) + t_to_origin = dtype[1.0 0.0 1*x_cen; 0.0 1.0 y_cen; 0.0 0.0 1.0]; + t_to_center = dtype[1.0 0.0 -1.0*x_cen; 0.0 1.0 -1.0*y_cen; 0.0 0.0 1.0]; + + true_vals = dtype[rand(-4.0:0.001:4.0), rand(-4.0:0.001:4.0), 1.0, 1.0, 0.0, 0.0, 0.0];#rand(0.9:0.001:1.1),rand(0.9:0.001:1.1), 0.0, 0.0, 0.0];#rand(0.001:0.001:pi/2.001)] + + if !pure_rand + + shear_mat = dtype[1.0 true_vals[5] 0.0; true_vals[6] 1.0 0.0; 0.0 0.0 1.0]; + scale_mat = dtype[1.0/true_vals[3] 0.0 0.0; 0.0 1/true_vals[4] 0.0; 0.0 0.0 1.0]; + + shift_mat = dtype[1.0 0.0 true_vals[1]; 0.0 1.0 true_vals[2]; 0.0 0.0 1.0]; + # converting the data to function (DataToFunctions.get_function) + + rot_mat = dtype[cos(true_vals[7]) -1.0*sin(true_vals[7]) 0.0; sin(true_vals[7]) cos(true_vals[7]) 0.0; 0.0 0.0 1.0]; + + matrix_c = (t_to_origin * scale_mat * shear_mat * rot_mat * shift_mat * t_to_center) + else + matrix_c = dtype.(t_to_origin * rand(0.1:0.001:1.0, (3, 3)) * t_to_center) + end + + f_affine_sim_img = get_function_affine(sample_data);#; super_sampling=1);#, extrapolation_bc=0.0); + if matrix + t_img = f_affine_sim_img(SMatrix{3, 3}(matrix_c))#, fitting_data); #.+ dtype.(rand(size(sample_data)...))./5.0; + else + t_img = f_affine_sim_img(true_vals)#, fitting_data); #.+ dtype.(rand(size(sample_data)...))./5.0; + end + fitting_data = dtype.(poisson(Float64.(t_img ./ maximum(t_img) .* n_photons))) #.+= rand(dtype, size(p_img)...).*noise_level; + + heatmap(fitting_data, aspect_ratio=1, size=(600, 600), title="fitting data", titlefont = font(20), legend=:none, axis=([], false)) + annotate!(vec(map(x -> Tuple((reverse(Tuple(x))..., text(@sprintf("%.0f", fitting_data[x]), :center, font(5), :white))), CartesianIndices(sample_data)))) + savefig("figures/fitting/sample_data_1.png") + + # plot(heatmap(sample_data, aspect_ratio=1), heatmap(fitting_data, aspect_ratio=1)) + + contour(sample_data, length=200, fill=false, title="Sample data", titlefont = font(20), legend=:none, aspect_ratio=1, size=(600, 600)) + savefig("figures/fitting/sample_data_1_contour.png") + return sample_data, fitting_data +end + + +function gauss_psf_comp() + x = -10.0:0.01:10.0 + p = Normal(0.0, 1.0) + y = pdf(p, x) + y_psf(x) = (sin(x) /x)^2 + plot(x, y_psf.(x), label="PSF of a circular aperture", title="Comparison of a Gaussian and a PSF", titlefont=20, size=(800, 400)) + plot!(x, y./maximum(y), label="Gaussian with μ=0.0 & σ=1.0") + savefig("figures/fitting/comp_psf_gaussian_1.png") + + plot(x, map(x -> (gradient(y_psf, x)[1]), x), label="Gradient of the PSF") + plot!(x, map(x -> (gradient(x -> pdf(p, x), x)[1]), x), label="Gradient of the Gaussian", title="Comparison of the Gradients", titlefont=20, size=(800, 400)) + savefig("figures/fitting/comp_psf_gaussian_1_gradients.png") +end + +""" + main_fitting(;matrix=true, sz=64, dtype=Float32, iterations=20, noise_level=1/20.0, pure_rand=false, from_params=true) + +This function is designed for performing fitting operations on sample data using either matrix transformations or +parametric transformations based on the provided arguments. + +# Arguments +`matrix`: if true, the fitting is done using a matrix transformation, otherwise, the fitting is done using a parametric transformation +`sz`: the size of the sample data +`dtype`: the data type of the sample data +`iterations`: the number of iterations to perform the fitting +`noise_level`: the noise level to add to the sample data +`pure_rand`: if true, the fitting is done using a random matrix transformation +`from_params`: if true, the fitting is done using the true values as the initial values + +# Returns +a tuple of two arrays: the first array is the fitting data, and the second array is the estimated fitting data + +# Example +x, y = main_fitting(matrix=true, sz=32, dtype=Float32, iterations=10, noise_level=1/20.0, pure_rand=false, from_params=true); + + +""" +function main_fitting(;matrix=true, sz=64, dtype=Float32, iterations=20, use_psf=true, n_photons=1000, pure_rand=false, from_params=true, plotting=false) + Random.seed!(14) + + # defining the mean and the varixance of the test normal (Gaussian) distribution + μ = [0, 0] + Σ = [sz/30 0.0; + 0.0 sz/30] + + # initializing the multivariate normal distribution + p = MvNormal(μ, Σ) + + # to define the sample array based on a 2D normal distribution + X = -1*sz/2.0:1*sz/2.0 + Y = -1*sz/2.0:1*sz/2.0 + + z = [pdf(p, [x,y]) for y in Y, x in X] + + + # creating a PSF for the widefield microscope + sz_psf = (sz, sz, 100) + sampling = (0.040, 0.040, 0.050) + # simulate a confocal PSF + aberrations = Aberrations([Zernike_HorizontalComa],[0.8]); + pp = PSFParams(0.5,1.4,1.52, method=MethodPropagateIterative, aberrations=aberrations); + + #pp_ex = PSFParams(pp_em; λ=0.488);#, method=MethodPropagateIterative, aplanatic=aplanatic_illumination, aberrations=aberrations); + p_psf_3d = psf(sz_psf, pp, sampling=sampling); + p_psf = p_psf_3d[:, :, 50] + #sample_data = p_psf ./ maximum(p_psf) + + # normalizing the sample data + sample_data = p_psf ./ maximum(p_psf) + sample_data_gaussian = dtype.(z[1:sz, 1:sz]./maximum(z)) + # sample_data = dtype.(TestImages.shepp_logan(sz)) + # sample_data = rand(dtype, (sz, sz)) + + # sample_data .+= rand(dtype, (size(sample_data)...)).*noise_level; + p_img = n_photons .* (sample_data);# ./ maximum(sample_data)) + n_img = dtype.(poisson(Float64.(p_img))) + + y = similar(n_img, (size(n_img)..., iterations)) + x = similar(n_img, (size(n_img)..., iterations)) + pos_res = zeros(Float32, iterations, 2) + pos_arr = zeros(Float32, iterations, 2) + + for i in ProgressBar(1:iterations) + # println("iteration: ", i) + + x_cen, y_cen = (size(n_img) ./ 2.0) + t_to_origin = dtype[1.0 0.0 1*x_cen; 0.0 1.0 y_cen; 0.0 0.0 1.0]; + t_to_center = dtype[1.0 0.0 -1.0*x_cen; 0.0 1.0 -1.0*y_cen; 0.0 0.0 1.0]; + + true_vals = dtype[rand(-4.0:0.001:4.0), rand(-4.0:0.001:4.0), 1.0, 1.0, 0.0, 0.0, 0.0];#rand(0.9:0.001:1.1),rand(0.9:0.001:1.1), 0.0, 0.0, 0.0];#rand(0.001:0.001:pi/2.001)] + + if !pure_rand + + shear_mat = dtype[1.0 true_vals[5] 0.0; true_vals[6] 1.0 0.0; 0.0 0.0 1.0]; + scale_mat = dtype[1.0/true_vals[3] 0.0 0.0; 0.0 1/true_vals[4] 0.0; 0.0 0.0 1.0]; + + shift_mat = dtype[1.0 0.0 true_vals[1]; 0.0 1.0 true_vals[2]; 0.0 0.0 1.0]; + # converting the data to function (DataToFunctions.get_function) + + rot_mat = dtype[cos(true_vals[7]) -1.0*sin(true_vals[7]) 0.0; sin(true_vals[7]) cos(true_vals[7]) 0.0; 0.0 0.0 1.0]; + + matrix_c = (t_to_origin * scale_mat * shear_mat * rot_mat * shift_mat * t_to_center) + else + matrix_c = dtype.(t_to_origin * rand(0.1:0.001:1.0, (3, 3)) * t_to_center) + end + + f_affine_sim_img = get_function_affine(sample_data);#; super_sampling=1);#, extrapolation_bc=0.0); + if matrix + t_img = f_affine_sim_img(SMatrix{3, 3}(matrix_c))#, fitting_data); #.+ dtype.(rand(size(sample_data)...))./5.0; + else + t_img = f_affine_sim_img(true_vals)#, fitting_data); #.+ dtype.(rand(size(sample_data)...))./5.0; + end + fitting_data = dtype.(poisson(Float64.(t_img ./ maximum(t_img) .* n_photons))) #.+= rand(dtype, size(p_img)...).*noise_level; + + + if use_psf + f_affine = get_function_affine(p_img);#; super_sampling=1);#, extrapolation_bc=0.0); + else + f_affine = get_function_affine(n_photons .* sample_data_gaussian);#; super_sampling=1);#, extrapolation_bc=0.0); + end + #f_affine = get_function_affine(n_img);#; super_sampling=1);#, extrapolation_bc=0.0); + # defining the loss function based on the gaussian noise + loss_m(p1::AbstractMatrix) = sum(abs2.(f_affine(SMatrix{size(p1)...}(p1)) .- fitting_data)) + # loss_m(p1::AbstractVector) = sum(abs2.(f_affine(p1::AbstractVector) .- fitting_data)) + loss_m(p1::AbstractVector) = sum(abs2.(f_affine([p1[1], p1[2], 0.0, 0.0, 0.0, 0.0, 0.0]) .- fitting_data)) + + + if matrix + if from_params + st_vals = dtype[1.0 0.0 -1.0*(argmax(fitting_data)[1]-size(fitting_data)[1]/2.0); 0.0 1.0 -1.0*(argmax(fitting_data)[1]-size(fitting_data)[2]/2.0); 0.0 0.0 1.0] + else + st_vals = dtype[1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0] + end + elseif from_params + st_vals = dtype[argmax(fitting_data)[1]-size(fitting_data)[1]/2.0, argmax(fitting_data)[2]-size(fitting_data)[2]/2.0, 1.0, 1.0, 0.0, 0.0, 0.0];#pi/8.0] + else + st_vals = dtype[0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0];#pi/8.0] + end + + # perform the main fit to the fitting data by minimizing the loss function + stats = @timed output, res = perform_fit(loss_m, st_vals) + + + if !matrix + pos_arr[i, :] = true_vals[1:2] + pos_res[i, :] = output[1:2] + else + pos_arr[i, :] = matrix_c[1:2, 3] + pos_res[i, :] = output[1:2, 3] + end + + y[:, :, i] = matrix ? f_affine(SMatrix{size(matrix_c)...}(output)) : f_affine(output) + x[:, :, i] = fitting_data + + # plotting the output of the fitting pocedure for further illustration + if plotting + begin + p00 = heatmap(n_img, aspect_ratio=1.0, title="Simulated sample PSF", colormap= :gist_gray); + p01 = heatmap(fitting_data, aspect_ratio=1.0, title="Simulated PSF", colormap= :gist_gray); + p02 = heatmap(matrix ? f_affine(SMatrix{size(matrix_c)...}(output)) : f_affine(output), aspect_ratio=1.0, title="Estimated fit", colormap= :gist_gray); + p03 = heatmap(fitting_data .- (matrix ? f_affine(SMatrix{size(matrix_c)...}(output)) : f_affine(output)), aspect_ratio=1.0, title="Residuals", colormap= :bwr, clim=(-maximum((abs.(fitting_data .- (matrix ? f_affine(SMatrix{size(matrix_c)...}(output)) : f_affine(output))))), maximum((abs.(fitting_data .- (matrix ? f_affine(SMatrix{size(matrix_c)...}(output)) : f_affine(output))))))); + + plot(p00, p01, p03, p02, layout=@layout([A B; C D]), + #framestyle=nothing, + #showaxis=false, + #xticks=false, yticks=false, + size=(1200, 1200), + plot_title=" $(if matrix "Matrix" else "Parametric" end) fitting +True vals: $(map(x -> @sprintf("%.3f",x), (matrix ? matrix_c : true_vals))) +fitted vals: $(map(x -> @sprintf("%.3f",x), output)) +n. of photons: $(@sprintf("%.0f", n_photons)) +time elapsed: $(@sprintf("%.1f", 1000.0*stats.time))ms, loss: $(@sprintf("%.2f", res.trace[1].value)) -> $(@sprintf("%.2f", res.trace[end].value))", + plot_titlevspan=0.14 + ) + savefig("figures/fitting/$(matrix ? "Matrix" : "Parametric")_fitting_$(i).png") + end + end + end + return x, y, pos_arr, pos_res +end + +x, y, pos_arr, pos_res = main_fitting(matrix=false, iterations=1000, sz=65, pure_rand=false, n_photons=10, from_params=true, plotting=false); println(mean(pos_res[:, 1] .- pos_arr[:, 1])); +println(std(pos_res[:, 1] .- pos_arr[:, 1])); + +#, title="Positional errors", markersize=2.0, xlabel="X error (pixels)", ylabel="Y error (pixels)", legend=:none, size=(600, 600), xlim=(-1.0, 1.0), ylim=(-1.0, 1.0), alpha=0.4) + +histogram2d(pos_res[:, 1] .- pos_arr[:, 1], pos_res[:, 2] .- pos_arr[:, 2], title="Positional errors histogram for 100 photons", xlabel="X error (pixels)", ylabel="Y error (pixels)", xlim=(-1.0, 1.0), ylim=(-1.0, 1.0), bins=20, aspect_ratio=1) +# + + +imgg = Gray{N0f16}.(x./maximum(x)); +ff = TiffImages.DenseTaggedImage(imgg); +TiffImages.save("test_4_aberrated.tif", ff); + +res_fiji = CSV.File(open(raw"C:\Users\ho82nat\Desktop\thunderstorm_res_100photons.csv")) +res_fiji_x = res_fiji["x [nm]"] ./ 80.0 .- 65.0 ./ 2.0; +res_fiji_y = res_fiji["y [nm]"] ./ 80.0 .- 65.0 ./ 2.0; + + +scatter(pos_res[:, 1] .- pos_arr[:, 1], pos_res[:, 2] .- pos_arr[:, 2], markershape= :circle, title="Positional errors", markersize=2.0, xlabel="X error (pixels)", ylabel="Y error (pixels)", legend=:none, size=(600, 600), label="DataToFunctions fitting")#, xlim=(-1.0, 1.0), ylim=(-1.0, 1.0), alpha=0.4) +#scatter!(res_fiji_y .- pos_arr[:, 1], res_fiji_x .- pos_arr[:, 2], markershape= :rect, markersize=2.0, alpha=0.2, label="ThunderSTORM fitting") + +println((std(pos_res[:, 1] .- pos_arr[:, 1]), std(pos_res[:, 2] .- pos_arr[:, 2])), (mean(pos_res[:, 1] .- pos_arr[:, 1]), mean(pos_res[:, 2] .- pos_arr[:, 2]))); +#println((std(res_fiji_y .- pos_arr[:, 1]), std(res_fiji_x .- pos_arr[:, 2])), (mean(res_fiji_y .- pos_arr[:, 1]), mean(res_fiji_x .- pos_arr[:, 2]))); + +#anim = @animate for i1 in 1:length(Optim.x_trace(res)) +# +# begin +# p00 = heatmap(sample_data, aspect_ratio=1.0, clim=(0.0, 1.0), title="Sample data", legend = :none); +# p01 = heatmap(fitting_data, aspect_ratio=1.0, clim=(0.0,1.0), title="Fitting data", legend = :none); +# p02 = heatmap(f_general(Optim.x_trace(res)[i1]), aspect_ratio=1.0, clim=(0.0,1.0), title="estimated fit", legend = :none); +# p03 = heatmap(fitting_data .- f_general(Optim.x_trace(res)[i1]), aspect_ratio=1.0, clim=(0.0, 0.3), title="discrepancy", legend = :none); +# +# plot(p00, p01, p02, p03, layout=@layout([A B C D]), +# framestyle=nothing, showaxis=false, +# xticks=false, yticks=false, +# size=(1200, 500), +# plot_title="iteration: $(Int(i1))/$(length(Optim.x_trace(res))), +# estimation: $(Optim.x_trace(res)[i1]) +# true vals : $(true_vals)", +# plot_titlevspan=0.25 +# ) +# end +# +# +#end; +# +#gif(anim, "DataToFunctions.jl/examples/anim_general_generalized.mp4", fps=2) +# diff --git a/examples/PSF_fitting_new.jl b/examples/PSF_fitting_new.jl new file mode 100644 index 0000000..dc9012f --- /dev/null +++ b/examples/PSF_fitting_new.jl @@ -0,0 +1,17 @@ +using PointSpreadFunctions +using Plots + +λ_em = 0.5; NA = 1.4; n = 1.52 +λ_ex = 0.488 # only needed for some PointSpreadFunctions, such as confocal, ISM or TwoPhoton +pp = PSFParams(λ_em, NA, n; pol=pol_x) + +sz = (256, 256, 256) +sampling = (0.020,0.020,0.020) + +aberr_sp = Aberrations([Zernike_VerticalAstigmatism],[1.0]); sz=(256,256,256) +pp_sp = PSFParams(λ_em, NA, n; method=MethodPropagateIterative, aberrations= aberr_sp) +p_sp = psf(sz, pp_sp; sampling=sampling); + +psf_example = sum(p_sp, dims=3)[:,:,1] + +heatmap(p_sp[:, :, 128], aspect_ratio=1) \ No newline at end of file diff --git a/examples/Project.toml b/examples/Project.toml new file mode 100644 index 0000000..5393e21 --- /dev/null +++ b/examples/Project.toml @@ -0,0 +1,44 @@ +[deps] +Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" +BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" +CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298" +DataToFunctions = "64cfdffa-4d02-49ee-ae8b-a805370874f5" +Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f" +Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9" +EvalMultiPoly = "c78649ec-f9ed-405e-be6d-0472f43586aa" +ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" +FourierTools = "b18b359b-aebc-45ac-a139-9c0ccbb2871e" +GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a" +ImageShow = "4e3cecfd-b093-5904-9786-8bbb286a6a31" +Images = "916415d5-f1e6-5110-898d-aaa5f9f070e0" +IndexFunArrays = "613c443e-d742-454e-bfc6-1d7f8dd76566" +Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59" +InverseModeling = "ce844058-9528-415d-a63d-06f3dd08b29f" +LBFGSB = "5be7bae1-8223-5378-bac3-9e7378a2f6e6" +LineSearches = "d3d80556-e9d4-5f37-9878-2ab0fcc64255" +LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a" +NDTools = "98581153-e998-4eef-8d0d-5ec2c052313d" +Noise = "81d43f40-5267-43b7-ae1c-8b967f377efa" +OhMyREPL = "5fb14364-9ced-5910-84b2-373655c76a03" +Optim = "429524aa-4258-5aef-a3af-852621145aeb" +PProf = "e4faabce-9ead-11e9-39d9-4379958e3056" +Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" +Pluto = "c3e4b0f8-55cb-11ea-2926-15256bba5781" +PlutoUI = "7f904dfe-b85e-4ff6-b463-dae2292396a8" +PointSpreadFunctions = "e8810a93-244e-46c5-8da3-35c5dd956001" +Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" +ProgressBars = "49802e3a-d2f1-5c88-81d8-b72133a6f568" +RandomExtensions = "fb686558-2515-59ef-acaa-46db3789a887" +Rotations = "6038ab10-8711-5258-84ad-4b1120ba62dc" +SeparableFunctions = "c8c7ead4-852c-491e-a42d-3d43bc74259e" +StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" +SyntheticObjects = "e7028c27-0967-45e9-8fdb-dbc10ccb2b0a" +TaylorSeries = "6aa5eb33-94cf-58f4-a9d0-e4b2c4fc25ea" +TestImages = "5e47fb64-e119-507b-a336-dd2b206d9990" +TiffImages = "731e570b-9d59-4bfa-96dc-6df516fadf69" +TimerOutputs = "a759f4b9-e2f1-59dc-863e-4aeb61b1ea8f" +View5D = "90d841e0-6953-4e90-9f3a-43681da8e949" +Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" diff --git a/examples/Zygote_adding.jl b/examples/Zygote_adding.jl new file mode 100644 index 0000000..d23ef1d --- /dev/null +++ b/examples/Zygote_adding.jl @@ -0,0 +1,47 @@ +using Zygote +using DataToFunctions + +data = rand(11,10) +f = get_function(data; super_sampling=2); + +loss(p, z) = sum(abs2.(f(p, z) .- data)) + +# Zygote.forwarddiff needs only one input +loss(p2::Vector{Tuple{Float64, Float64}}) = loss(p2[1], p2[2]) + + +for i in 1:20 + pr = [(0.0, -0.1+i/100.0), (1.0, 1.0)] + + # to take derivative of a interpolation process, it is better to use the Zygote.forwarddiff + loss_grad = Zygote.forwarddiff(loss, pr) + + # it can be seen that increasing the scale variable (y direction) leads to increase in the loss function gradient value + println(loss_grad) + +end + +array_scale = Array{Tuple{Float64, Float64}, 2}(undef, 200, 2) +list_loss_grad_scale = Array{Float64, 2}(undef, 200, 200) + +x_1 = 0.9:0.01:1.1 +y_1 = 0.9:0.01:1.1 + + +for i in 1:200 + for j in 1:200 + array_scale[i, :] = [(0.0, 0.0), (0.9 + j/1000.0, 0.9 + i/1000.0)] + + #pr = [(0.0, 0.0), loss_scale[1, :]] + + # to take derivative of a interpolation process, it is better to use the Zygote.forwarddiff + list_loss_grad_scale[i, j] = Zygote.forwarddiff(loss, array_scale[i, :]) + + # it can be seen that increasing the scale variable (y direction) leads to increase in the loss function gradient value + #println(loss_grad) + end + +end + +surface(1:200, 1:200, list_loss_grad_scale) + diff --git a/examples/deform_testimage.jl b/examples/deform_testimage.jl new file mode 100644 index 0000000..3e96f79 --- /dev/null +++ b/examples/deform_testimage.jl @@ -0,0 +1,160 @@ +### A Pluto.jl notebook ### +# v0.19.42 + +using Markdown +using InteractiveUtils + +# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error). +macro bind(def, element) + quote + local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end +end + +# ╔═╡ 28975586-853e-4e19-b9eb-65c41fa61a43 +using Pkg + +# ╔═╡ 0ae2da4f-3f75-47bb-a899-9e89c5c3f17c +Pkg.activate(".") + +# ╔═╡ 4af0c13d-fc42-4fe7-97e6-2248e36b63e2 +Pkg.add("PlutoUI") + +# ╔═╡ a2d75cfb-feab-4130-8439-30c543618d04 +using DataToFunctions, ImageShow, TestImages, PlutoUI, Images + +# ╔═╡ 1c744bec-f085-4812-ab1a-40a32c2ac176 +import PlutoUI: combine + +# ╔═╡ 5ac1123d-5df3-4c9d-aff1-ffe91d931497 +data = Float32.(testimage("resolution_test_512")) + +# ╔═╡ b0c8d15e-1bd3-4e66-b2b9-57885636eb48 +simshow(data) + +# ╔═╡ 2fce0208-732f-4259-a47d-7f78921bfd87 +f = get_interpolated_function(data, AffineMode, super_sampling=1) + +# ╔═╡ dd535378-3f2a-4429-914c-8b7608b99706 +@bind shift_x Slider(-100:0.02:100, default=0) + +# ╔═╡ 3e26d4e8-b5d4-4414-8936-379ec63cb4d2 +@bind shift_y Slider(-100:0.02:100, default=0) + +# ╔═╡ bc3f3659-aa70-4e7f-bef6-4d05229a03c4 +@bind zoom_x Slider(0.2:0.02:4, default=1) + +# ╔═╡ 39f11dc5-351e-4c6d-8fc4-468222e99976 +@bind zoom_y Slider(0.2:0.02:4, default=1) + +# ╔═╡ 64b64d3b-092f-4553-916d-f7db1fdfa428 +simshow(f((shift_x, shift_y, 1/zoom_x, 1/zoom_y, 0.0, 0.0, 0.0))) + +# ╔═╡ c3fe6b25-f4c0-4a80-9d15-9ee30136d43b +typeof(f((shift_x, shift_y, 1/zoom_x, 1/zoom_y, 0.0, 0.0, 0.0))) + +# ╔═╡ ff143f0d-b070-4220-8fa6-0b5a93a56303 +typeof(data) + +# ╔═╡ 6676ba35-3efd-49f5-9819-411ad8f8a95c +md""" +# Polynomial transformations +""" + +# ╔═╡ 13461a95-95ea-4bad-8673-e94e06776254 +md""" +## First order polynomial + +First order polynomial transformation which is as follows: + +``x^{\prime} = c_{1} + c_{2}{x}^{1} + c_{3}{y}^{1}`` + +``y^{\prime} = c_{4} + c_{5}{x}^{1} + c_{6}{y}^{1}`` +""" + +# ╔═╡ 87be45c0-8b2e-4d49-abd1-a274b3c1815e +h = get_interpolated_function(data, PolynomialMode, 1); + +# ╔═╡ 68f771aa-2cde-41cd-990c-9ec7dc2146a4 +function coeffs_input(coeffs::Vector) + + return combine() do Child + + inputs = [ + md""" $(name): $( + Child(name, Slider(-2f0:0.05f0:2f0, default=0, show_value=true)) + )""" + + for name in coeffs + ] + + md""" + #### Transform coefficients + $(inputs) + """ + end +end; + +# ╔═╡ 69331d73-75a3-4727-acda-e79779a2bd03 +@bind c coeffs_input(["c1", "c2", "c3", "c4", "c5", "c6"]) + +# ╔═╡ 74228a9d-6cc2-4aaf-97da-67f32670341e +simshow(h((c.c1, c.c2, c.c3, c.c4, c.c5, c.c6)), cmap=:turbo)#,0f0,0f0,0f0,0f0,0f0,0f0,0f0,0f0,0f0))) + +# ╔═╡ 687a198b-9020-4959-8e27-fd0896d4b1fc +maximum(h((c.c1,c.c2,c.c3,c.c4,c.c5,c.c6))) + +# ╔═╡ 6584aacc-440e-457b-bf52-83f8db40c999 +h((c.c1,c.c2,c.c3,c.c4,c.c5,c.c6)) + +# ╔═╡ 74cf9f27-f559-4fe2-9a30-20cb7cf6fe81 +md""" +## Second order polynomial + +Second order polynomial transformation is as follows: + +``x^{\prime} = c_{1} + c_{2}{x}^{1} + c_{3}{x}^{2} + c_{4}{y}^{1} + c_{5}{x}{y} + c_{6}{y}^{2}`` + +``y^{\prime} = c_{7} + c_{8}{x}^{1} + c_{9}{x}^{2} + c_{10}{y}^{1} + c_{11}{x}{y} + c_{12}{y}^{2}`` +""" + +# ╔═╡ 1f440592-9bbb-429c-8ba3-d018f5b354b0 +h2 = get_interpolated_function(data, PolynomialMode, 2); + +# ╔═╡ 12f59a6d-4600-4379-b536-f3b4701bdbe4 +@bind c2 coeffs_input(["c1", "c2", "c3", "c4", "c5", "c6", "c7", "c8", "c9", "c10", "c11", "c12"]) + +# ╔═╡ fd60df55-2a80-48a9-95ec-eeeb1cfe7491 +simshow(h2((c2.c1, c2.c2, c2.c3, c2.c4, c2.c5, c2.c6, c2.c7, c2.c8, c2.c9, c2.c10, c2.c11, c2.c12)), cmap=:turbo)#,0f0,0f0,0f0,0f0,0f0,0f0,0f0,0f0,0f0))) + +# ╔═╡ Cell order: +# ╠═28975586-853e-4e19-b9eb-65c41fa61a43 +# ╠═0ae2da4f-3f75-47bb-a899-9e89c5c3f17c +# ╟─4af0c13d-fc42-4fe7-97e6-2248e36b63e2 +# ╠═1c744bec-f085-4812-ab1a-40a32c2ac176 +# ╠═a2d75cfb-feab-4130-8439-30c543618d04 +# ╠═5ac1123d-5df3-4c9d-aff1-ffe91d931497 +# ╠═b0c8d15e-1bd3-4e66-b2b9-57885636eb48 +# ╠═2fce0208-732f-4259-a47d-7f78921bfd87 +# ╠═dd535378-3f2a-4429-914c-8b7608b99706 +# ╠═3e26d4e8-b5d4-4414-8936-379ec63cb4d2 +# ╠═bc3f3659-aa70-4e7f-bef6-4d05229a03c4 +# ╠═39f11dc5-351e-4c6d-8fc4-468222e99976 +# ╠═64b64d3b-092f-4553-916d-f7db1fdfa428 +# ╠═c3fe6b25-f4c0-4a80-9d15-9ee30136d43b +# ╠═ff143f0d-b070-4220-8fa6-0b5a93a56303 +# ╟─6676ba35-3efd-49f5-9819-411ad8f8a95c +# ╟─13461a95-95ea-4bad-8673-e94e06776254 +# ╠═87be45c0-8b2e-4d49-abd1-a274b3c1815e +# ╠═69331d73-75a3-4727-acda-e79779a2bd03 +# ╠═68f771aa-2cde-41cd-990c-9ec7dc2146a4 +# ╠═74228a9d-6cc2-4aaf-97da-67f32670341e +# ╠═687a198b-9020-4959-8e27-fd0896d4b1fc +# ╠═6584aacc-440e-457b-bf52-83f8db40c999 +# ╟─74cf9f27-f559-4fe2-9a30-20cb7cf6fe81 +# ╠═1f440592-9bbb-429c-8ba3-d018f5b354b0 +# ╠═12f59a6d-4600-4379-b536-f3b4701bdbe4 +# ╠═fd60df55-2a80-48a9-95ec-eeeb1cfe7491 diff --git a/examples/deform_testimgage backup 1.jl b/examples/deform_testimgage backup 1.jl new file mode 100644 index 0000000..2a37851 --- /dev/null +++ b/examples/deform_testimgage backup 1.jl @@ -0,0 +1,50 @@ +### A Pluto.jl notebook ### +# v0.19.43 + +using Markdown +using InteractiveUtils + +# ╔═╡ 28975586-853e-4e19-b9eb-65c41fa61a43 +using Pkg + +# ╔═╡ a2d75cfb-feab-4130-8439-30c543618d04 +using DataToFunctions, ImageShow, TestImages + +# ╔═╡ 0ae2da4f-3f75-47bb-a899-9e89c5c3f17c +# Pkg.activate(".") + +# ╔═╡ 4af0c13d-fc42-4fe7-97e6-2248e36b63e2 +# Pkg.add("ImageShow") + +# ╔═╡ 5ac1123d-5df3-4c9d-aff1-ffe91d931497 +data = testimage("resolution_test_512", super_sampling=1) + +# ╔═╡ b0c8d15e-1bd3-4e66-b2b9-57885636eb48 + + +# ╔═╡ 2fce0208-732f-4259-a47d-7f78921bfd87 +f = get_function(data) + +# ╔═╡ 1bed34fb-b29a-4042-a493-4835fdb69a9d + + +# ╔═╡ c287fa80-426b-11ef-125e-5fda207e605c +# ╠═╡ disabled = true +#=╠═╡ +using DataToFunctions + ╠═╡ =# + +# ╔═╡ 87be45c0-8b2e-4d49-abd1-a274b3c1815e + + +# ╔═╡ Cell order: +# ╠═28975586-853e-4e19-b9eb-65c41fa61a43 +# ╠═0ae2da4f-3f75-47bb-a899-9e89c5c3f17c +# ╠═4af0c13d-fc42-4fe7-97e6-2248e36b63e2 +# ╠═a2d75cfb-feab-4130-8439-30c543618d04 +# ╠═5ac1123d-5df3-4c9d-aff1-ffe91d931497 +# ╠═b0c8d15e-1bd3-4e66-b2b9-57885636eb48 +# ╠═2fce0208-732f-4259-a47d-7f78921bfd87 +# ╠═1bed34fb-b29a-4042-a493-4835fdb69a9d +# ╠═c287fa80-426b-11ef-125e-5fda207e605c +# ╠═87be45c0-8b2e-4d49-abd1-a274b3c1815e diff --git a/examples/deformed_fitting.jl b/examples/deformed_fitting.jl new file mode 100644 index 0000000..cdd35d2 --- /dev/null +++ b/examples/deformed_fitting.jl @@ -0,0 +1,42 @@ +using Interpolations +using FourierTools + + +function get_function(data::AbstractArray; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) + new_size = super_sampling.*size(data) + upsampled = fftshift(resample(ifftshift(data), new_size)) + @show upsampled + # return upsampled + # itp = LinearInterpolation(axes(upsampled), upsampled, extrapolation_bc=extrapolation_bc); + interpolation = Interpolations.interpolate(upsampled, interp_type) + interpolation = extrapolate(interpolation, extrapolation_bc) + # center of the original data (too keep the axis and number of datapointsi dentical to the original) + center_orig = (size(data) .÷2 .+1) + # create zero-centered original ranges (== axes) + zero_axes = Tuple(ax .- c for (ax, c) in zip(axes(data), center_orig)) + # center of the upsampled data. This is where to access the upsampled data + function zoomed(shift, zoom, theta) + zoom = zoom .* super_sampling + # careful: The center of the original data is not at the expected position! But rather at: + center_upsamp = new_size .÷2 .+1 # ((center_orig .-1) .*super_sampling .+1) # new_size .÷2 .+1 + scaled_axes = ((ax.-myc) .* z .+ cen for (ax, myc, cen, z) in zip(zero_axes, shift, center_upsamp, zoom)) + # @show Tuple(scaled_axes) + return interpolation[scaled_axes...] + # return extrapolate(scale(interpolation, scaled_axes...), extrapolation_bc) + end + zoomed(p) = zoomed([p[1], p[2]], [p[3], p[4]], [p[5]]) + + return zoomed + + # return (pos) -> interp_linear((center .+ pos)...) + # fitp(t) = interp_linear(t...) + # @time res1 = fitp.(tcoords); # 1 sec + # function my_zoom + +end + + +f_d = get_function(sample_data_d; super_sampling=2, extrapolation_bc=0.0) + + +f_d(true_vals); \ No newline at end of file diff --git a/examples/perform_fitting.jl b/examples/perform_fitting.jl new file mode 100644 index 0000000..c12925f --- /dev/null +++ b/examples/perform_fitting.jl @@ -0,0 +1,117 @@ +using DataToFunctions +using Optim, StaticArrays, LinearAlgebra +using Zygote +using ForwardDiff, LineSearches, Plots, Printf + + +Base.show(io::IO, f::Float64) = @printf(io, "%.2f", f) + +true_vals = [1.0, 2.0, 1.0, 1.0] +init_x = [0.5, 1.5, 1.0, 1.0] + +sample_data = rand(11,12) + +f = get_function(sample_data; super_sampling=2, extrapolation_bc=0.0); +#f(p0::Vector{Float64}) = f([p0[1], p0[2], p0[3], p0[4]]) + +fitting_data = f(true_vals) .+ rand(11, 12)./10.0 +#f(p2[1], p2[2]) = f(p2::Vector{Tuple{Float64, Float64}}) +loss(p, z) = sum(abs2.(f(p, z) .- fitting_data)) +loss(p2::Vector{Tuple{Float64, Float64}}) = loss(p2[1], p2[2]) +loss(p3) = loss([p3[1], p3[2]], [p3[3], p3[4]]) + + + + + +#Zygote.forwarddiff(loss, init_x) + +ForwardDiff.gradient(loss, init_x)#true_vals .+ [0.0001, 0.0, 0.0, 0.0]) +#conf = ForwardDiff.GradientConfig(f, init_x, chunk::Chunk = Chunk(init_x)) + + +""" +BFGS(; alphaguess = Optim.LineSearches.InitialStatic(), + linesearch = Optim.LineSearches.HagerZhang(), + initial_invH = nothing, + initial_stepnorm = 0.001, + manifold = Optim.Flat() + ) + +GradientDescent(; alphaguess = 0.01, + linesearch = Optim.LineSearches.HagerZhang(), + P = nothing, + precondprep = (P, x) -> nothing +) +""" +lower = [-1*size(fitting_data)[1], -1*size(fitting_data)[2], 0.0, 0.0] +upper = [size(fitting_data)[1], size(fitting_data)[2], size(fitting_data)[1], size(fitting_data)[2]] +#initial_x = [2.0, 2.0] +# requires using LineSearches +inner_optimizer = LBFGS(; m=1, linesearch=LineSearches.BackTracking(order=2)) +res = optimize( + loss, + lower, upper, + init_x, + Fminbox(inner_optimizer), + Optim.Options(store_trace = true, extended_trace = true, iterations=500), + autodiff = :forward +) + + +p00 = heatmap(sample_data, aspect_ratio=1.0, clim=(0.0, 1.0), title="Sample data", legend = :none); +p01 = heatmap(fitting_data, aspect_ratio=1.0, clim=(0.0,1.0), title="Fitting data", legend = :none); +p02 = heatmap(f(Optim.minimizer(res)), aspect_ratio=1.0, clim=(0.0,1.0), title="estimated fit", legend = :none); +p03 = heatmap(fitting_data .- f(Optim.minimizer(res)), aspect_ratio=1.0, clim=(0.0,1.0), title="discrepancy", legend = :none); + +plot(p00, p01, p02, p03, layout=@layout([A B C D]), + framestyle=nothing, showaxis=false, + xticks=false, yticks=false, + size=(700, 300), + plot_title="True vals: $(true_vals)", plot_titlevspan=0.2 +) + +heatmap(fitting_data .- f(Optim.x_trace(res)[end]), aspect_ratio=1, clim=(0.0, 1.0)) + +""" + +res = optimize( + loss, init_x, + LBFGS(), + Optim.Options(store_trace=true, extended_trace=true, iterations=500), + autodiff = :forward + ) +""" + +trace = Optim.trace(res); +trace + + + + +Optim.minimizer(res) +Optim.f_trace(res) +Optim.x_trace(res) +Optim.converged(res) +Optim.g_norm_trace(res) +Optim.g_calls(res) + +loss(Optim.x_trace(res)[end]) +ForwardDiff.gradient(loss, Optim.x_trace(res)[end]) + + + +anim = @animate for i1 in 1:length(Optim.x_trace(res)) + + heatmap(fitting_data .- f(Optim.x_trace(res)[i1]), + aspect_ratio=1, + clim=(0.0, 1.0), + dpi=300 + ) + title!("iteration: $(Int(i1))/$(length(Optim.x_trace(res))), + estimation: $(Optim.x_trace(res)[i1]) + true vals : $(true_vals)") + +end; + +gif(anim, "anim1.mp4", fps=5) diff --git a/examples/perform_random_optim.jl b/examples/perform_random_optim.jl new file mode 100644 index 0000000..b7a1b69 --- /dev/null +++ b/examples/perform_random_optim.jl @@ -0,0 +1,176 @@ +using DataToFunctions +using Optim, StaticArrays, LinearAlgebra +using Zygote +using ForwardDiff, LineSearches, Plots, Printf +using ProfileView, Profile +using Random, Distributions + +Random.seed!(123) + +# to show all the numbers in 2 decimals format +Base.show(io::IO, f::Float64) = @printf(io, "%.2f", f) + +### +# creating the random matrix in size = (11, 12) +sample_data = rand(11, 12) + +# creating the lower and upper bounds of the fitting variables (shift and scaling) +# shift can not be higher than the size of the array, +# scale can not be lower than zero or higher than size of the data, the latter causes the resulting array to be just one pixel +lower = [-1*size(sample_data)[1], -1*size(sample_data)[2], 0.0, 0.0] +upper = [size(sample_data)[1], size(sample_data)[2], size(sample_data)[1], size(sample_data)[2]] + +# preparing the function to fit +f = get_function(sample_data; super_sampling=2, extrapolation_bc=0.0); + +# assigning a scale (multiplier) to the range of true values, wedo not want that the parameters of +# function to be near the limits and cause strange behavior of the optimization +scale_range = 4.0 + +# true values of the fitting are random for the repeatibility +true_vals = (rand(4) .* (upper .- lower)/scale_range ) .+ (lower/scale_range) #[3.0, 5.5, 1.85, 0.6] + +# create the fitting data and adding random noise with scale of 1/10 +d = Normal() +noise = rand(d, 11, 12) + + +fitting_data = f(true_vals) #.+ noise / 10.0 +heatmap(fitting_data, aspect_ratio=1.0) +savefig("fitting_data_without_noise.png") + +# create the loss function and using the Julia's multiple dispatch +# because the gradients function is required just one input (can be vector) +loss(p, z) = sum(abs2.(f(p, z) .- fitting_data)) +loss(p2::Vector{Tuple{Float64, Float64}}) = loss(p2[1], p2[2]) +loss(p3) = loss([p3[1], p3[2]], [p3[3], p3[4]]) + + + + +# initialization of the LBFGS optimizer +inner_optimizer = LBFGS(; m=1, linesearch=LineSearches.BackTracking(order=2)) + + + +function perform_optim_mthr(loss, n_walkers) + """ + defining a function to survey the parameter space to neglect the local + minima and find the global maxima + + loss: the loss function + n_walkers: number of random initial parameter estimation + """ + # allocating the estimation array: + # consists of four parameters [1:4] and the minimum loss function of them [5] + est_m = zeros(n_walkers, 5) + + #x_tr = Array{Any} + #res = Array{Optim.MultivariateOptimizationResults{}}(undef, n_walkers, 1) + #res = Optim.MultivariateOptimizationResults{} + + # defining the random initial parameter values + walkers = (rand(4, n_walkers) .* (upper .- lower)/scale_range ) .+ (lower/scale_range) + + # main loop to do the optimization for each of the initial parameter values + # it uses the Threads to distribute the for loop to each thread + # note that in the settings.json the Julia is started with 16 threads + Threads.@threads for i in 1:size(walkers)[2] + res = optimize( + loss, + lower, upper, # the limits (simple box constraints) + walkers[:, i], + Fminbox(inner_optimizer), # assigning the limits of fitting (simple constraints) along with the LBFGS + Optim.Options(store_trace = true, extended_trace = true, iterations=500), + autodiff = :forward + ); + #x_tr = Optim.x_trace(res) + + # saving each 4 parameters of the fit and the minimum loss function + est_m[i, 1:4] .= Optim.minimizer(res)#x_tr[end] + est_m[i, 5] = minimum(res)#loss(x_tr[end]) + end + # saving the parameters of the minimum of the loss function + ans_m = est_m[argmin(est_m[:, 5]), :] + return est_m, ans_m +end + +#@profile est_m, ans_m = perform_optim_mthr(loss, 10000) + + +# first time: 4.749537 seconds (31.85 M allocations: 3.468 GiB, 8.50% gc time, 88.55% compilation time) +# 2nd time: 0.645307 seconds (12.68 M allocations: 2.525 GiB, 56.30% gc time) +# changing the fitting_data (noise values) : 0.448227 seconds (9.90 M allocations: 1.984 GiB, 27.55% compilation time: 38% of which was recompilation) +@time est_m, ans_m = perform_optim_mthr(loss, 500); + + +println(string(true_vals) * "\n" * string(ans_m)) + +# [2.23, -1.16, 1.26, 2.59] +# [2.26, -1.15, 1.27, 2.59, 0.41] + + +p00 = heatmap(sample_data, aspect_ratio=1.0, clim=(0.0, 1.0), title="Sample data", legend = :none); +p01 = heatmap(fitting_data, aspect_ratio=1.0, clim=(0.0,1.0), title="Fitting data", legend = :none); +p02 = heatmap(f(ans_m[1:4]), aspect_ratio=1.0, clim=(0.0,1.0), title="estimated fit", legend = :none); +p03 = heatmap(fitting_data .- f(ans_m[1:4]), aspect_ratio=1.0, clim=(0.0,1.0), title="discrepancy", legend = :none); + +plot(p00, p01, p02, p03, layout=@layout([A B C D]), + framestyle=nothing, showaxis=false, + xticks=false, yticks=false, + size=(700, 300), + plot_title="True vals: $(true_vals)", plot_titlevspan=0.2 +) +savefig("Output_mth.png") + + +plot( + heatmap(fitting_data .- f(ans_m[1:4]), aspect_ratio=1.0, clim=(0.0, 0.3), title="discrepancy", legend = :none), + heatmap(noise / 10.0 , aspect_ratio=1.0, clim=(0.0, 0.3), title="noise", legend = :none), + layout=@layout([A B]), + framestyle=nothing, showaxis=false, + xticks=false, yticks=false, + size=(700, 300), + plot_title="Noise comparison", plot_titlevspan=0.2 +) +savefig("Discrepancy.png") + +p12 = plot(fitting_data[1, :], color="black", legend=:none); +t = fitting_data[1, :] +for i in 2:size(fitting_data)[1] + t .+= fitting_data[i, :] + plot!(p12, fitting_data[i, :], color="black", legend=:none) +end +display(p12) + +plot(t) + + +function perform_optim_sthr(loss, n_walkers=1000) + est = Array{Float64, 2}(undef, n_walkers, 5) + + walkers = (rand(4, n_walkers) .* (upper .- lower)/scale_range ) .+ (lower/scale_range) + for i in 1:size(walkers)[2] + res = optimize( + loss, + lower, upper, + walkers[:, i], + Fminbox(inner_optimizer), + Optim.Options(store_trace = true, extended_trace = true, iterations=500), + autodiff = :forward + ); + x_tr = Optim.x_trace(res) + est[i, 1:4] .= x_tr[end] + est[i, 5] = loss(x_tr[end]) + end + ans1 = est[argmin(est[:, 5]), :] + return est, ans1 +end + + +@time est_s, ans_s = perform_optim_sthr(loss, 20000); + +println(string(true_vals) * "\n" * string(ans_s)) +# [1.07, -2.63, 0.78, 1.65] +# [1.02, -2.60, 0.77, 1.66, 0.46] + diff --git a/examples/poly_test_exp.jl b/examples/poly_test_exp.jl new file mode 100644 index 0000000..202a7b3 --- /dev/null +++ b/examples/poly_test_exp.jl @@ -0,0 +1,172 @@ +using Images +using DataToFunctions +using FindShift +using Optim, CUDA +using FourierTools +using View5D, Plots, Statistics + +CUDA.allowscalar(false) + +file1 = raw"D:\Hossein\Programming\Julia\DataToFunctions.jl\examples\test_polim1.jpeg" +file2 = raw"D:\Hossein\Programming\Julia\DataToFunctions.jl\examples\test_polim2.jpeg" + +c1 = Float32.(Gray.(load(file1))) +c2 = Float32.(Gray.(load(file2))) + +#TODO increase the size of the images +img11 = c1[30:1519+30, 30:779+30] +img12 = c1[30:1519+30, 800:779+800] +img21 = c2[30:1519+30, 30:779+30] +img22 = c2[30:1519+30, 800:779+800] + +img11_c = CuArray(img11) +img12_c = CuArray(img12) +img21_c = CuArray(img21) +img22_c = CuArray(img22) + +function resize_img(data, scale) + new_size = size(data) .÷ scale + imresize(data, new_size...) +end + +resample_size = 2 +img1_resampled = img11[1:resample_size:end, 1:resample_size:end] #resize_img(img1, 1) +img2_resampled = img12[1:resample_size:end, 1:resample_size:end] #resize_img(img2, 1) + +#im11 = imfilter(img1_resampled, Kernel.gaussian(5)); +#im12 = imfilter(img2_resampled, Kernel.gaussian(5)); + +f = get_interpolated_function(img12_c, PolynomialMode, 2); +loss_p(p1::AbstractArray) = (sum(abs2.(f(Tuple(p1)) .- img11_c))); +loss_updated(p1::AbstractArray) = mapreduce(abs2, +, f(Tuple(p1)) .- img11_c); +#st_vals = [0.0, 1.00, 0.00, 0.0, 0.00, 1.0] #ones(Float64, 6)./10 +st_vals = Float32[0.0, 1.0, 0, 0, 0, 0, 0.0, 0, 0, 1.0, 0, 0]; +# Float64[9.0, 0, 0, 0, 0, 0, 1, 0, 0, 5.0, 0, 0, 1, 0, 0, 0, 0, 0] +# @vv f(Tuple(st_vals)) +# loss_m(st_vals) + +#a = f(Tuple(st_vals)) +# @time a = f(Tuple(st_vals)); +""" +function do_registeration_step(resample_size, gaussian_kernel_size, img1, img2, st_vals=[0f0, 1f0, 0f0, 0f0, 0f0, 1f0]) + img1_resampled = img1[1:resample_size:end, 1:resample_size:end] #resize_img(img1, 1) + img2_resampled = img2[1:resample_size:end, 1:resample_size:end] #resize_img(img2, 1) + + im1 = imfilter(img1_resampled, Kernel.gaussian(gaussian_kernel_size)); + im2 = imfilter(img2_resampled, Kernel.gaussian(gaussian_kernel_size)); + + im1_c = CuArray(im1) + im2_c = CuArray(im2) + + f = get_interpolated_function(im2_c, PolynomialMode, 1); + loss_updated(p1::AbstractArray) = mapreduce(abs2, +, f(Tuple(p1)) .- im1_c); + + #st_vals = Float32[0.0, 1.0, 0, 0, 0, 0, 0.0, 0, 0, 1.0, 0, 0]; + CUDA.@time res = optimize( + loss_updated, + st_vals, + BFGS(),#; initial_stepnorm = 1f-1),#; linesearch=LineSearches.BackTracking(order=2)), + autodiff = :forward + ) + #return [res.minimizer[1]*resample_size, 1.0, 0, 0, 0, 0, res.minimizer[7]*resample_size, 0, 0, 1.0, 0, 0] + return [res.minimizer[1]*resample_size, 1.0, 0, res.minimizer[4]*resample_size, 0, 1.0] +end + +res_step1 = do_registeration_step(10, 5, img11, img12) +@vt img11 get_interpolated_function(img12, PolynomialMode, 1)(Tuple(res_step1)) img12 +""" + + +""" +function g!(G, x) # (G, x) + G .= gradient(loss_updated, x)[1] +end +od = OnceDifferentiable(loss_updated, g!, st_vals) +""" +aligned_imgs = Array(copy(img11_c)) + + +CUDA.@time res = optimize( + loss_updated, + st_vals, + #Newton(), + BFGS(; initial_stepnorm = 1f-2),#; linesearch=LineSearches.BackTracking(order=2)), + #LBFGS(; linesearch=LineSearches.BackTracking(order=3)), + #lower, upper, + #init_x, + #Fminbox(inner_optimizer), + #Optim.Options(store_trace = true, extended_trace = true, iterations=5000), + autodiff = :forward +) +""" +10.271416 seconds (8.78 M CPU allocations: 589.752 MiB, 1.06% gc time) (2.25 k GPU allocations: 52.251 GiB, 0.30% memmgmt time) +* Status: success + +* Candidate solution +Final objective value: 2.966678e+03 + +* Found with +Algorithm: BFGS + +* Convergence measures +|x - x'| = 1.76e-05 ≰ 0.0e+00 +|x - x'|/|x'| = 3.62e-06 ≰ 0.0e+00 +|f(x) - f(x')| = 0.00e+00 ≤ 0.0e+00 +|f(x) - f(x')|/|f(x')| = 0.00e+00 ≤ 0.0e+00 +|g(x)| = 1.17e+05 ≰ 1.0e-01 + +* Work counters +Seconds run: 5 (vs limit Inf) +Iterations: 73 +f(x) calls: 450 +∇f(x) calls: 450 +""" +aligned_imgs = cat(aligned_imgs, Array(f(Tuple(Optim.minimizer(res)))), dims=3); + + +for img_t in [img21_c, img22_c] + f = get_interpolated_function(img_t, PolynomialMode, 2); + loss_p(p1::AbstractArray) = (sum(abs2.(f(Tuple(p1)) .- img11_c))); + loss_updated(p1::AbstractArray) = mapreduce(abs2, +, f(Tuple(p1)) .- img11_c); + #st_vals = Float64[0.0, 1.00, 0.00, size(img11)[2], 0.00, -1.0] #ones(Float64, 6)./10 + st_vals = Float32[0.0, 1.0, 0, 0, 0, 0, size(img11)[2], 0, 0, -1.0, 0, 0]; + # Float64[9.0, 0, 0, 0, 0, 0, 1, 0, 0, 5.0, 0, 0, 1, 0, 0, 0, 0, 0] + # @vv f(Tuple(st_vals)) + # loss_m(st_vals) + + """ + function g!(G, x) # (G, x) + G .= gradient(loss_p, x)[1] + end + od = OnceDifferentiable(loss_p, g!, st_vals) + """ + CUDA.@time res = optimize( + loss_updated, + st_vals, + #Newton(), + BFGS(; initial_stepnorm = 1f-2),#; linesearch=LineSearches.BackTracking(order=2)), + #LBFGS(),#; linesearch=LineSearches.BackTracking(order=3)), + #lower, upper, + #init_x, + #Fminbox(inner_optimizer), + #Optim.Options(store_trace = true, extended_trace = true, iterations=5000), + autodiff = :forward + ) + + + aligned_imgs = cat(aligned_imgs, Array(f(Tuple(Optim.minimizer(res)))), dims=3); +end + + +aligned_imgs = clamp.(aligned_imgs, 0, 1) +@vv aligned_imgs + +aligned_imgs = permutedims(aligned_imgs, [2, 1, 3]) +save("aligned_all_cuda_polyorder2.gif", Gray.(aligned_imgs)) + + + +raw_imgs = cat(img11, img12, img21, img22, dims=3) +raw_imgs = clamp.(raw_imgs, 0, 1) +raw_imgs = permutedims(raw_imgs, [2, 1, 3]) +save("raw_all_imgs.gif", Gray.(raw_imgs)) diff --git a/examples/poly_test_exp_new.jl b/examples/poly_test_exp_new.jl new file mode 100644 index 0000000..b5f8233 --- /dev/null +++ b/examples/poly_test_exp_new.jl @@ -0,0 +1,173 @@ +using Images +using DataToFunctions +using FindShift +using Optim, CUDA +using FourierTools, Zygote +using View5D, Plots, Statistics, LineSearches + +CUDA.allowscalar(false) + +file1 = raw"D:\Hossein\Programming\Julia\DataToFunctions.jl\examples\markerpen_C1_00001.tif" +file2 = raw"D:\Hossein\Programming\Julia\DataToFunctions.jl\examples\markerpen_C2_00001.tif" + +c1 = Float32.(Gray.(load(file1))) +c2 = Float32.(Gray.(load(file2))) + +wide=true +if wide + img11 = c1[180:1990+180, 160:2050+160] + img12 = c1[160:1990+160, 2180:2050+2180] + img21 = c2[130:1990+130, 140:2050+140] + img22 = c2[130:1990+130, 2160:2050+2160] +end + +@vt img11 img12 img21 img22 + + + +img11_c = CuArray(img11./maximum(img11)) +img12_c = CuArray(img12./maximum(img12)) +img21_c = CuArray(img21./maximum(img21)) +img22_c = CuArray(img22./maximum(img22)) + + +resample_size = 2 +img1_resampled = img11[1:resample_size:end, 1:resample_size:end] #resize_img(img1, 1) +img2_resampled = img12[1:resample_size:end, 1:resample_size:end] #resize_img(img2, 1) + +#im11 = imfilter(img1_resampled, Kernel.gaussian(5)); +#im12 = imfilter(img2_resampled, Kernel.gaussian(5)); + + +""" +resample_size=50 +img1_resampled = img11[1:resample_size:end, 1:resample_size:end] #resize_img(img1, 1) +img2_resampled = img12[1:resample_size:end, 1:resample_size:end] #resize_img(img2, 1) + +im1 = imfilter(img1_resampled, Kernel.gaussian(3)); +im2 = imfilter(img2_resampled, Kernel.gaussian(3)); + +im1_c = CuArray(im1./maximum(im1)) +im2_c = CuArray(im2./maximum(im2)) + +f1 = get_interpolated_function(im2_c, PolynomialMode, 2); +loss_updated1(p1::AbstractArray) = mapreduce(abs2, +, f(Tuple(p1)) .- im1_c); + +st_vals = Float32[0.0, 1.0, 0, 0, 0, 0, 0.0, 0, 0, 1.0, 0, 0]; +CUDA.@time res1 = optimize( + loss_updated1, + st_vals, + BFGS(; initial_stepnorm = 1f-2),#; linesearch=LineSearches.BackTracking(order=2)), + autodiff = :forward +) + +@vt Array(im1_c) Array(f1(Tuple(res1.minimizer))) Array(im2_c) + +function do_registeration_step(resample_size, gaussian_kernel_size, img1, img2, st_vals=Float32[0.0, 1.0, 0, 0, 0, 0, 0.0, 0, 0, 1.0, 0, 0]) + img1_resampled = img1[1:resample_size:end, 1:resample_size:end] #resize_img(img1, 1) + img2_resampled = img2[1:resample_size:end, 1:resample_size:end] #resize_img(img2, 1) + + im1 = imfilter(img1_resampled, Kernel.gaussian(gaussian_kernel_size)); + im2 = imfilter(img2_resampled, Kernel.gaussian(gaussian_kernel_size)); + + im1_c = CuArray(im1) + im2_c = CuArray(im2) + + f = get_interpolated_function(im2_c, PolynomialMode, 2); + loss_updated(p1::AbstractArray) = mapreduce(abs2, +, f(Tuple(p1)) .- im1_c); + + #st_vals = Float32[0.0, 1.0, 0, 0, 0, 0, 0.0, 0, 0, 1.0, 0, 0]; + CUDA.@time res = optimize( + loss_updated, + st_vals, + BFGS(),#; initial_stepnorm = 1f-1),#; linesearch=LineSearches.BackTracking(order=2)), + autodiff = :forward + ) + return [res.minimizer[1]*resample_size, 1.0, 0, 0, 0, 0, res.minimizer[7]*resample_size, 0, 0, 1.0, 0, 0] + #return [res.minimizer[1]*resample_size, 1.0, 0, res.minimizer[4]*resample_size, 0, 1.0] +end + +res_step1 = do_registeration_step(10, 5, img11, img12) +@vt img11 get_interpolated_function(img12, PolynomialMode, 2)(Tuple(res_step1)) img12 +""" + + +order=2 +f = get_interpolated_function(img12_c, PolynomialMode, order); +loss_p(p1::AbstractArray) = (sum(abs2.(f(Tuple(p1)) .- img11_c))); +loss_updated(p1::AbstractArray) = mapreduce(abs2, +, f(Tuple(p1)) .- img11_c); +if order ==1 + st_vals = Float32[0.0, 1.00, 0.00, 0.0, 0.00, 1.0] #ones(Float64, 6)./10 +elseif order == 2 + st_vals = Float32[0.0, 1.0, 0, 0, 0, 0, 0.0, 0, 0, 1.0, 0, 0]; +end + + +aligned_imgs = Array(copy(img11_c)) + + +CUDA.@time res = optimize( + loss_updated, + st_vals, + #Newton(), + BFGS(; initial_stepnorm = 1f-1),#; linesearch=LineSearches.BackTracking(order=2)), + #LBFGS(; linesearch=LineSearches.BackTracking(order=3)), + #lower, upper, + #init_x, + #Fminbox(inner_optimizer), + Optim.Options(store_trace = true, extended_trace = true, iterations=5000, g_tol=1f-2), + autodiff = :forward +) + +open("markerpen_new_data.txt", "w") do f + write(f, "\norder $(order) params = $(res.minimizer)") +end + +aligned_imgs = cat(aligned_imgs, Array(f(Tuple(res.minimizer))), dims=3); +@vv aligned_imgs + + +for img_t in [img21_c, img22_c] + f = get_interpolated_function(img_t, PolynomialMode, 2); + loss_p(p1::AbstractArray) = (sum(abs2.(f(Tuple(p1)) .- img11_c))); + loss_updated(p1::AbstractArray) = mapreduce(abs2, +, f(Tuple(p1)) .- img11_c); + #st_vals = Float64[0.0, 1.00, 0.00, size(img11)[2], 0.00, -1.0] #ones(Float64, 6)./10 + #st_vals = Float32[0.0, 1.0, 0, 0, 0, 0, size(img11)[1], 0, 0, -1.0, 0, 0]; + if order ==1 + st_vals = Float32[0.0, 1.00, 0.00, 0.0, 0.00, 1.0] #ones(Float64, 6)./10 + elseif order == 2 + st_vals = Float32[0.0, 1.0, 0, 0, 0, 0, 0.0, 0, 0, 1.0, 0, 0]; + end + + CUDA.@time res = optimize( + loss_updated, + st_vals, + #Newton(), + BFGS(; initial_stepnorm = 1f-1),#; linesearch=LineSearches.BackTracking(order=2)), + #LBFGS(),#; linesearch=LineSearches.BackTracking(order=3)), + #lower, upper, + #init_x, + #Fminbox(inner_optimizer), + Optim.Options(store_trace = true, extended_trace = true, iterations=5000, g_tol=1f-2), + autodiff = :forward + ) + println(res) + open("markerpen_new_data.txt", "a") do f + write(f, "\norder $(order) params = $(res.minimizer)") + end + aligned_imgs = cat(aligned_imgs, Array(f(Tuple(Optim.minimizer(res)))), dims=3); +end + + +aligned_imgs = clamp.(aligned_imgs, 0, 1) +@vv aligned_imgs + +aligned_imgs = permutedims(aligned_imgs, [2, 1, 3]) +save("aligned_all_cuda_polyorder2_markerpen_new.gif", Gray.(aligned_imgs)) + + + +raw_imgs = cat(img11, img12, img21, img22, dims=3) +raw_imgs = clamp.(raw_imgs, 0, 1) +raw_imgs = permutedims(raw_imgs, [2, 1, 3]) +save("raw_all_imgs_markerpen_new.gif", Gray.(raw_imgs)) diff --git a/examples/polynomial_apply.jl b/examples/polynomial_apply.jl new file mode 100644 index 0000000..501502b --- /dev/null +++ b/examples/polynomial_apply.jl @@ -0,0 +1,93 @@ +using DataToFunctions +using Optim, StaticArrays, LinearAlgebra +using PointSpreadFunctions +using Zygote +using ForwardDiff, LineSearches, Plots, Printf +using View5D +using Distributions, Rotations +using Plots +using TestImages +using BenchmarkTools +#using InverseModeling +import Random +using Noise +using IndexFunArrays + + + +function poly_test(;sz=64, dtype=Float64, n_photons=1000) + + sz_psf = (sz, sz, 100) + sampling = (0.040, 0.040, 0.050) + # simulate a confocal PSF + #aberrations = Aberrations([Zernike_HorizontalComa,Zernike_Tip],[0.8,0.7]); + pp = PSFParams(0.5,1.4,1.52, method=MethodPropagateIterative);#, aberrations=aberrations); + + #pp_ex = PSFParams(pp_em; λ=0.488);#, method=MethodPropagateIterative, aplanatic=aplanatic_illumination, aberrations=aberrations); + p_psf_3d = psf(sz_psf, pp, sampling=sampling); + p_psf = p_psf_3d[:, :, 50] + #sample_data = p_psf ./ maximum(p_psf) + + # normalizing the sample data + sample_data = p_psf ./ maximum(p_psf) * n_photons + + #sample_data = make_grid(); + f_1 = get_interpolated_function(Float64.(sample_data), PolynomialMode, 1); # get_function_affine(sample_data); + true_vals = (1.6, 1.05, 0.1, 1.5, 0.01, 1.02) # dtype[2.0, 1.0, 1.01, 1.0, 0.0, 0.0, 0.0] + # true_vals = dtype[rand(-4.0:0.001:4.0), rand(-4.0:0.001:4.0), rand(0.5:0.001:1.5),rand(0.5:0.001:1.5), 0.0, 0.0, rand(0.001:0.001:pi/2.001)] + + dat2 = f_1(true_vals) + + #sample_data = dtype.(TestImages.shepp_logan(sz)) + f = get_interpolated_function(Float64.(sample_data), PolynomialMode, 1); + + + + loss_p(p1::AbstractArray) = (sum(abs2.(f(Tuple(p1)) .- dat2))) + st_vals = [2.1, 1.00, 0.00, 1.5, 0.00, 1.0] #ones(Float64, 6)./10 + #st_vals = Float64[1.0, 0, 0, 0, 0, 0, 1.0, 0, 0, 1.0, 0, 0, 1.0, 0, 0, 0, 0, 0] + # Float64[9.0, 0, 0, 0, 0, 0, 1, 0, 0, 5.0, 0, 0, 1, 0, 0, 0, 0, 0] + # @vv f(Tuple(st_vals)) + # loss_m(st_vals) + + + function g!(G, x) # (G, x) + G .= gradient(loss_p, x)[1] + end + od = OnceDifferentiable(loss_p, g!, st_vals) + res = optimize( + loss_p, + st_vals, + #Newton(), + BFGS(; initial_stepnorm = 1e-2),#; linesearch=LineSearches.BackTracking(order=2)), + #LBFGS(),#; linesearch=LineSearches.BackTracking(order=3)), + #lower, upper, + #init_x, + #Fminbox(inner_optimizer), + #Optim.Options(store_trace = true, extended_trace = true, iterations=5000, g_tol=1e-3), + autodiff = :forward + ) + + # return the estimated parameters + return true_vals, Optim.minimizer(res), res, f, f_1 +end + +#a, b, c = poly_test() +#@vt f(Tuple(b)) f_affine(a) + + +function make_grid!(arr::AbstractArray) + arr[isinteger.(xx(size(arr))./10) .|| isinteger.(yy(size(arr))./10)] .= 1.0 + return arr + +end +function make_grid(sz::NTuple{N, Int}=(64, 64)) where {N} + arr = zeros(Float64, sz) + make_grid!(arr) + return arr +end +function make_grid(::Type{T}, sz::NTuple{N, Int}=(64, 64)) where {N, T} + arr = zeros(T, sz) + make_grid!(arr) + return arr +end \ No newline at end of file diff --git a/examples/star_fitting.jl b/examples/star_fitting.jl new file mode 100644 index 0000000..0bd63c0 --- /dev/null +++ b/examples/star_fitting.jl @@ -0,0 +1,88 @@ +using DataToFunctions +using Optim, StaticArrays, LinearAlgebra +using Zygote +using ForwardDiff, LineSearches, Plots, Printf +using View5D +using Distributions +using Plots + +using Rotations +using CoordinateTransformations + +Base.show(io::IO, f::Float64) = @printf(io, "%.3f", f) + + +# size of the test array to fit +size_arr = 600.0 + +# defining the mean and the variance of the test normal (Gaussian) distribution +μ = [0, 0] +Σ = [size_arr/2.0 0.0; + 0.0 size_arr/2.0] + +Σ_d = [size_arr/0.5 1.5; + 1.5 size_arr/6.0] + +# initializing the multivariate normal distribution +p = MvNormal(μ, Σ) +p_d = MvNormal(μ, Σ_d) + + + +# this part of the code is to define the sample array based on a 2D normal distribution +X = -1*size_arr/2.0:1*size_arr/2.0 +Y = -1*size_arr/2.0:1*size_arr/2.0 + +z = [pdf(p, [x,y]) for y in Y, x in X] +z_d = [pdf(p_d, [x,y]) for y in Y, x in X] + +#@vv z + +heatmap(z, aspect_ratio=1) +heatmap(z_d, aspect_ratio=1) + +# setting a typical values for the shift (1:2) and scale (3:4) +true_vals = [1.1, -1.5, 0.75, 1.5] + +# normalizing the sample data +sample_data = z./maximum(z) +sample_data_d = z_d./maximum(z_d) + +# converting the data to function (DataToFunctions.get_function) +f = get_function(sample_data; super_sampling=2, extrapolation_bc=0.0); +f_d = get_function(sample_data_d; super_sampling=2, extrapolation_bc=0.0); + +# adding some scaled random noise to the fitting data +fitting_data = f(true_vals) .+ rand(size(z)...)./8.0 + +# @vv fitting_data + +# defining the loss function based on the gaussian noise +loss(p) = sum(abs2.(f(p) .- fitting_data)) + +#loss(p3) = loss([p3[1], p3[2]], [p3[3], p3[4]]) + +heatmap(fitting_data, aspect_ratio=1, title="Fitting data") + +# perform the main fit to the fitting data by minimizing the loss function +output = perform_fit(loss, fitting_data) + +# plotting the output of the fitting pocedure for further illustration +begin + p00 = heatmap(sample_data, aspect_ratio=1.0, clim=(0.0, 1.0), title="Sample data", legend = :none); + p01 = heatmap(fitting_data, aspect_ratio=1.0, clim=(0.0,1.0), title="Fitting data", legend = :none); + p02 = heatmap(f(output), aspect_ratio=1.0, clim=(0.0,1.0), title="estimated fit", legend = :none); + p03 = heatmap(fitting_data .- f(output), aspect_ratio=1.0, clim=(0.0, 0.3), title="discrepancy", legend = :none); + + plot(p00, p01, p02, p03, layout=@layout([A B C D]), + framestyle=nothing, showaxis=false, + xticks=false, yticks=false, + size=(700, 300), + plot_title="True vals: $(true_vals) + est vals: $(output)", + plot_titlevspan=0.25 + ) +end + +# comparing the true values to the best fitting parameters +println(string(true_vals) * "\n" * string(output)) diff --git a/examples/star_rotation.jl b/examples/star_rotation.jl new file mode 100644 index 0000000..5704ec0 --- /dev/null +++ b/examples/star_rotation.jl @@ -0,0 +1,90 @@ +using DataToFunctions +using Optim, StaticArrays, LinearAlgebra, FourierTools +using Zygote +using ForwardDiff, LineSearches, Plots, Printf +using View5D +using Distributions +using Plots + +using Rotations +using CoordinateTransformations + +using Interpolations + +function get_function2(data::AbstractArray; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) + new_size = super_sampling.*size(data) + upsampled = fftshift(resample(ifftshift(data), new_size)) + + # return upsampled + # itp = LinearInterpolation(axes(upsampled), upsampled, extrapolation_bc=extrapolation_bc); + interpolation = Interpolations.interpolate(upsampled, interp_type) + interpolation = extrapolate(interpolation, extrapolation_bc) + # center of the original data (too keep the axis and number of datapointsi dentical to the original) + center_orig = (size(data) .÷2 .+1) + # create zero-centered original ranges (== axes) + zero_axes = Tuple(ax .- c for (ax, c) in zip(axes(data), center_orig)) + # center of the upsampled data. This is where to access the upsampled data + cen_vec = SVector(size(data)./ 2.0) + + function transform_axes(t::SVector, mat::SMatrix, cen_vec::SVector, shift_vec::SVector) + + return (mat * (t - cen_vec)) + cen_vec - shift_vec + + end + + function zoomed(shift_vec, zoom, theta) + zoom = zoom .* super_sampling + # # careful: The center of the original data is not at the expected position! But rather at: + # center_upsamp = new_size .÷2 .+1 # ((center_orig .-1) .*super_sampling .+1) # new_size .÷2 .+1 + # scaled_axes = ((ax.-myc) .* z .+ cen for (ax, myc, cen, z) in zip(zero_axes, shift, center_upsamp, zoom)) + + # #sh_x, sh_y = 0.0, 0.0 + s_x, s_y = zoom + theta = theta[1] + + rot_mat = @SMatrix [cos(theta) -1.0*sin(theta); sin(theta) cos(theta)]; + scale_mat = @SMatrix [s_x 0.0; 0.0 s_y]; + mat = rot_mat * scale_mat + + return interpolation.(transform_axes.(SVector.(Tuple.(CartesianIndices(data))), mat, cen_vec, SVector(shift_vec[1], shift_vec[2]))...) + # return extrapolate(scale(interpolation, scaled_axes...), extrapolation_bc) + end + zoomed(p) = zoomed([p[1], p[2]], [p[3], p[4]], p[5]); + # zoomed([p[1], p[2]], [p[3], p[4]], p[5]) = zoomed[p] + return zoomed + +end + +# defining the mean and the varixance of the test normal (Gaussian) distribution +μ = [0, 0] +Σ = [2 0.0; + 0.0 2] + +Σ_d = [2 1.5; + 1.5 2] + +# initializing the multivariate normal distribution +p = MvNormal(μ, Σ) +p_d = MvNormal(μ, Σ_d) + +# size of the test array to fit +size_arr = 12.0 + +# this part of the code is to define the sample array based on a 2D normal distribution +X = -1*size_arr/2.0:1*size_arr/2.0 +Y = -1*size_arr/2.0:1*size_arr/2.0 + +z = [pdf(p, [x,y]) for y in Y, x in X] +heatmap(z, aspect_ratio=1) + + + + +true_vals = [1.1, -1.5, 0.75, 1.5, pi/2] + + +sample_data = z./maximum(z) +f = get_function2(sample_data; super_sampling=2, extrapolation_bc=0.0); +f(true_vals) + +heatmap(f(true_vals), aspect_ratio=1.0, clim=(0.0,1.0), legend = :none); diff --git a/examples/taylor_test.jl b/examples/taylor_test.jl new file mode 100644 index 0000000..ad771f1 --- /dev/null +++ b/examples/taylor_test.jl @@ -0,0 +1,26 @@ +using DataToFunctions +using TestImages + +function main() + + obj = Float32.(TestImages.shepp_logan(320)); + + f = get_function_poly(obj,1) + c0 = (-10.5, 1.1, 0.1, -20.2, 0.1, 1.1) + @time warped = f(c0); # 1.7 ms + + fa = get_function_affine(obj); #.+ dtype.(rand(size(sample_data)...))./5.0; + ca = [1.5, 1.1, 0.6, 1.2, 0.1, 1.1, 2.0] + @time warpeda = fa(ca); # 1.7 ms + + # non-linear deformation warp + f2 = get_function_poly(obj, 2); #.+ dtype.(rand(size(sample_data)...))./5.0; + c02 = (-150.5, 1.1, 0.1, 0.001, 0.001 ,0.001, 0.001, 0.001, 0.001, + -110.2, 0.1, 1.5, 0.001, 0.0015,-0.001, 0.0014,0.001,-0.001) + @time warped2 = f2(c02); # 3.3 ms + @time warped2 .= f2(c02); + @time f2(c02, warped2); # 3.0 ms + # @vt obj warpeda warped warped2 + +end + diff --git a/src/DataToFunctions.jl b/src/DataToFunctions.jl index fd58afb..af2b99b 100644 --- a/src/DataToFunctions.jl +++ b/src/DataToFunctions.jl @@ -1,55 +1,14 @@ module DataToFunctions -using Interpolations -using FourierTools -export get_function -""" - get_function(data::AbstractArray; super_sampling=2, extrapolation_bc=Flat(), interp_type=Interpolations.BSpline(Linear())) +using Requires -returns a function `dat(shift, zoom)` which generates a shifted and scaled version of the original data. -This is useful for fitting with a function which is itself defined by measured data. +export gpu_or_cpu -# Arguments -`data`: The data to represent by the function `dat` -`super_sampling`: The factor by which the data is internally represented as a supersampled version (Fourier-based upsampling, see `FourierTools.resample`) -`extrapolation_bc`: The extrapolation boundary condition to select for values outside the range. - By default the value 0.0 is used. Other options are `Flat()`, or `Line()`, See the package `Interpolation` for details. -`interp_type`: The type of interpolation to use. See the package `Interpolation` for details. +# to include CUDA +include("requires.jl") -# Example -```jldoctest -``` -""" -function get_function(data::AbstractArray; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) - new_size = super_sampling.*size(data) - upsampled = fftshift(resample(ifftshift(data), new_size)) - @show upsampled - # return upsampled - # itp = LinearInterpolation(axes(upsampled), upsampled, extrapolation_bc=extrapolation_bc); - interpolation = Interpolations.interpolate(upsampled, interp_type) - interpolation = extrapolate(interpolation, extrapolation_bc) - # center of the original data (too keep the axis and number of datapointsi dentical to the original) - center_orig = (size(data) .÷2 .+1) - # create zero-centered original ranges (== axes) - zero_axes = Tuple(ax .- c for (ax, c) in zip(axes(data), center_orig)) - # center of the upsampled data. This is where to access the upsampled data - function zoomed(shift, zoom) - zoom = zoom .* super_sampling - # careful: The center of the original data is not at the expected position! But rather at: - center_upsamp = new_size .÷2 .+1 # ((center_orig .-1) .*super_sampling .+1) # new_size .÷2 .+1 - scaled_axes = ((ax.-myc) .* z .+ cen for (ax, myc, cen, z) in zip(zero_axes, shift, center_upsamp, zoom)) - # @show Tuple(scaled_axes) - return interpolation[scaled_axes...] - # return extrapolate(scale(interpolation, scaled_axes...), extrapolation_bc) - end - - return zoomed - - # return (pos) -> interp_linear((center .+ pos)...) - # fitp(t) = interp_linear(t...) - # @time res1 = fitp.(tcoords); # 1 sec - # function my_zoom - -end - -end # module DataToFunctions +#include("datafunction.jl") +include("transformation_types.jl") +include("transformators.jl") + +end # module DataToFunctions \ No newline at end of file diff --git a/src/datafunction.jl b/src/datafunction.jl new file mode 100644 index 0000000..d0063c1 --- /dev/null +++ b/src/datafunction.jl @@ -0,0 +1,34 @@ +import Base: getindex, setindex!, size, axes, eltype, copy, similar, ndims, iterate, length + + +struct DataFunctionAffine{affineparams{NTuple{N, Number}}, itp, S<:AbstractArray} <: AbstractArray + params::affineparams + interpolation::itp + data::S +end + +function DataFunctionAffine(params::affineparams, interpolation::itp) where {affineparams<:NTuple{N, Number}, itp} + return DataFunctionAffine{affineparams, itp}(params, interpolation) +end + +function getindex(dfa::DataFunctionAffine, I::Vararg{Number, N}) where {N} + return dfa.interpolation(Tuple(I)...) +end + +function copy(s::DataFunctionAffine) + res = similar(s) + res .= s +end + + +function similar(dfa::DataFunctionAffine, ::Type{T}=eltype(dfa)) where {T} + return DataFunctionAffine(dfa.params, similar(dfa.interpolation, T)) +end + +size(dfa::DataFunctionAffine) = size(dfa.data) +axes(dfa::DataFunctionAffine) = axes(dfa.data) +eltype(dfa::DataFunctionAffine) = eltype(dfa.data) +ndims(dfa::DataFunctionAffine) = ndims(dfa.data) +length(dfa::DataFunctionAffine) = length(dfa.data) + + diff --git a/src/requires.jl b/src/requires.jl new file mode 100644 index 0000000..f1eb5dd --- /dev/null +++ b/src/requires.jl @@ -0,0 +1,16 @@ +# from https://github.com/RainerHeintzmann/DeconvOptim.jl/blob/362a741224957155fbc046d3297d43339766721d/src/requires.jl + +isgpu(x) = false +gpu_or_cpu(x) = nothing + +function __init__() + @require CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" begin + @info "DataToFunctions.jl: CUDA.jl is loaded." + + gpu_or_cpu(x) = CUDA.CuArray{Float32} + isgpu(x::CUDA.CuArray) = true + # prevent slow scalar indexing on GPU + CUDA.allowscalar(false); + + end +end \ No newline at end of file diff --git a/src/transformation_types.jl b/src/transformation_types.jl new file mode 100644 index 0000000..119704d --- /dev/null +++ b/src/transformation_types.jl @@ -0,0 +1,4 @@ +abstract type transformation_method end + +struct AffineMode <: transformation_method end +struct PolynomialMode <: transformation_method end \ No newline at end of file diff --git a/src/transformators.jl b/src/transformators.jl new file mode 100644 index 0000000..1aced99 --- /dev/null +++ b/src/transformators.jl @@ -0,0 +1,441 @@ +using Interpolations +using FourierTools +using StaticArrays +using EvalMultiPoly + +export get_interpolated_function, get_function_tuple, get_function_svec, get_function_affine, get_function_poly +export add_dim, red_dim_apply, red_dim, mat_mul, func_transform +export extrapolate, interpolate + +export PolynomialMode, AffineMode + +""" + get_function(data::AbstractArray; super_sampling=2, extrapolation_bc=Flat(), interp_type=Interpolations.BSpline(Linear())) + +returns a function `dat(shift, zoom)` which generates a shifted and scaled version of the original data. +This is useful for fitting with a function which is itself defined by measured data. + +# Arguments +`data`: The data to represent by the function `dat` +`super_sampling`: The factor by which the data is internally represented as a supersampled version (Fourier-based upsampling, see `FourierTools.resample`) +`extrapolation_bc`: The extrapolation boundary condition to select for values outside the range. + By default the value 0.0 is used. Other options are `Flat()`, or `Line()`, See the package `Interpolation` for details. +`interp_type`: The type of interpolation to use. See the package `Interpolation` for details. + + +""" +function get_function_old(data::AbstractArray; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) + new_size = super_sampling.*size(data) + upsampled = fftshift(resample(ifftshift(data), new_size)) + # @show upsampled + # return upsampled + # itp = LinearInterpolation(axes(upsampled), upsampled, extrapolation_bc=extrapolation_bc); + interpolation = Interpolations.interpolate(upsampled, interp_type) + interpolation = extrapolate(interpolation, extrapolation_bc) + # center of the original data (too keep the axis and number of datapointsi dentical to the original) + center_orig = (size(data) .÷2 .+1) + # create zero-centered original ranges (== axes) + zero_axes = Tuple(ax .- c for (ax, c) in zip(axes(data), center_orig)) + # center of the upsampled data. This is where to access the upsampled data + function zoomed(shift, zoom) + zoom = zoom .* super_sampling + # careful: The center of the original data is not at the expected position! But rather at: + center_upsamp = new_size .÷2 .+1 # ((center_orig .-1) .*super_sampling .+1) # new_size .÷2 .+1 + scaled_axes = ((ax.-myc) .* z .+ cen for (ax, myc, cen, z) in zip(zero_axes, shift, center_upsamp, zoom)) + # @show Tuple(scaled_axes) + return interpolation[scaled_axes...] + # return extrapolate(scale(interpolation, scaled_axes...), extrapolation_bc) + end + + return zoomed + + zoomed(p) = zoomed([p[1], p[2]], [p[3], p[4]]) + # return (pos) -> interp_linear((center .+ pos)...) + # fitp(t) = interp_linear(t...) + # @time res1 = fitp.(tcoords); # 1 sec + # function my_zoom + +end + +""" + add_dim(cind) + +adds a dimension to a CartesianIndex + +`cind`: A CartesianIndex +""" +function add_dim(cind) + return SVector.((Tuple(cind))..., 1) +end + +""" + red_dim(svec::SVector{S,T}) + +removes the last dimension of a SVector to convert it from a homogeneous to a Cartesian coordinates + +`svec::SVector{S,T}`: A SVector +""" +@inline function red_dim(svec::SVector{S,T})::SVector{S-1,T} where {S,T} + return @view svec[1:S-1] +end + +""" + red_dim_apply(fct, svec::SVector{S,T}) + +applies a function to a SVector by removing the last dimension to convert it from a homogeneous to a Cartesian coordinates + +`fct`: The function to apply +`svec::SVector{S,T}`: A SVector +""" +@inline function red_dim_apply(fct, svec::SVector{S,T}) where {S,T} + return fct((@view svec[1:S-1])...) +end + +""" + red_dim_apply(fct, tup::NTuple{S,T}) + +applies a function to a Tuple by removing the last dimension to convert it from a homogeneous to a Cartesian coordinates + +`fct`: The function to apply +`tup::NTuple{S,T}`: A Tuple +""" +@inline function red_dim_apply(fct, tup::NTuple{S,T}) where {S,T} + return fct(tup[1:S-1]...) +end + +""" + idx_apply(fct, svec::SVector{S,T}) where {S,T} + +applies a function to a SVector + +`fct`: The function to apply +`svec::SVector{S,T}`: A SVector +""" +@inline function idx_apply(fct, svec::SVector{S,T})::Number where {S,T} + return fct(svec...) +end + +""" + idx_apply(fct, tup::NTuple{S,T}) where {S,T} + +applies a function to a Tuple + +`fct`: The function to apply +`tup::NTuple{S,T}`: A Tuple +""" +@inline function idx_apply(fct, tup::NTuple{S,T})::Number where {S,T} + return fct(tup...) +end + +""" + mat_mul(t::SVector{N, T2}, matrix_c::SMatrix{N,N,T}) + +multiplies a SVector with a SMatrix + +`t::SVector{N, T2}`: The SVector to multiply +`matrix_c::SMatrix{N,N,T}`: The SMatrix to multiply with +""" +@inline function mat_mul(t::SVector{N, T2}, matrix_c::SMatrix{N,N,T})::SVector{N,T} where {N,T, T2} + return matrix_c * t +end + +""" + func_transform(t, coord_transform_func::Function)::SVector + +applies a coordinate transformation function to an array or `CartesianIndex` and returns the transformed array + +`t`: The array or `CartesianIndex` to transform +`coord_transform_func::Function`: The function to apply the transformation +""" +@inline function func_transform(t, coord_transform_func::Function)::SVector + return coord_transform_func(Tuple(t)) +end + +""" + func_transform_tup(t, coord_transform_func::Function) + +applies a coordinate transformation function to a Tuple + +`t`: The array or `CartesianIndex` to transform +`coord_transform_func::Function`: The function to apply the transformation +""" +@inline function func_transform_tup(t, coord_transform_func::Function) + return coord_transform_func(Tuple(t)) +end + + +""" + apply_transform(coord_transf_func::Function, data::AbstractArray{T}, itp) where {T} + +applies a general coordinate transformation function to the indices of an array and returns the transformed array + +`coord_transf_func:Function`: A function that takes a N-+1 dimensional CartesianIndex and returns a new N+1 dimensional SVector or Tuple +`data::AbstractArray{T}`: The data to transform +`itp`: The interpolation object to use +""" +function apply_transform(coord_transf_func::Function, data::AbstractArray{T, N}, itp) where {T, N} #, out::AbstractArray{T}) where {T} + # @info "Applying tuple transformation" + + # return map((it) -> idx_apply(Interpolations.adapt(gpu_or_cpu(nothing), itp), coord_transf_func(Tuple(it))), CartesianIndices(data)) + return idx_apply.(Ref(Interpolations.adapt(gpu_or_cpu(nothing), itp)), coord_transf_func.(Tuple.(CartesianIndices(data)))); + #return idx_apply.(Ref(itp), coord_transf_func.(Tuple.(CartesianIndices(data)))); + # return idx_apply.(Ref(itp), coord_transf_func.(CartesianIndices(data))); +end + +""" + apply_transform_homogen(coord_transf_func::Function, data, itp) +applies a homogeneous coordinate-based coordinate transformation function to the indices of an array and returns the transformed array + +`coord_transf_func::Function`: A function that takes a N-+1 dimensional homogeneous SVector returns a new N+1 dimensional SVector +`data`: The data to transform +`itp`: The interpolation object to use +""" +function apply_transform_homogen(coord_transf_func::Function, data, itp)#, out) + h_coord_transf_func = (c) -> red_dim(coord_transf_func(add_dim(c))) + #@info "Applying homogeneous transformation" + return apply_transform(h_coord_transf_func, data, itp)#, out); + # out .= itp.(red_dim.(coord_transf_func.(add_dim.(CartesianIndices(data))))); + # out .= red_dim_apply.(Ref(itp), coord_transf_func.(add_dim.(CartesianIndices(data)))); +end + +""" + apply_transform_affine(mymat::SMatrix{T}, data, itp) where T + +applies an affine transformation matrix to the indices of an array and returns the transformed array + +`mymat::SMatrix{T}` The affine transformation matrix to apply +`data`: The data to transform +`itp`: The interpolation function (object) to use +""" +function apply_transform_affine(mymat::SMatrix{T}, data, itp) where T #, out) where {T} # The SMatrix spec is important to avoid allocations + # red_dim_apply.(Ref(itp), func_transform.(CartesianIndices(data), Ref(coord_transf_func))); + # return red_dim_apply.(Ref(itp), mat_mul.(add_dim.(CartesianIndices(data)), Ref(mymat))); + # @info "Applying affine transformation" + homogenous_transform = (c) -> mat_mul(c, mymat) + return apply_transform_homogen(homogenous_transform, data, itp)#, out); + #return out +end + +""" + get_function_tuple(data::AbstractArray, fct_tup::Function; super_sampling=2, extrapolation_bc=Flat(), interp_type=Interpolations.BSpline(Linear())) + +returns a function `dat(shift, zoom)` which generates a shifted and scaled version of the original data. +This is useful for fitting with a function which is itself defined by measured data. + +# Arguments +`data`: The data to represent by the function `dat` +`fct_tup`: The function to apply to the data +`super_sampling`: The factor by which the data is internally represented as a supersampled version (Fourier-based upsampling, see `FourierTools.resample`) +`extrapolation_bc`: The extrapolation boundary condition to select for values outside the range. + By default the value 0.0 is used. Other options are `Flat()`, or `Line()`, See the package `Interpolation` for details. +`interp_type`: The type of interpolation to use. See the package `Interpolation` for details. + +# Example + +""" +function get_function_tuple(data::AbstractArray{T, N}, fct_tup::Function; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) where {T, N} + # building the extraplation + interpolation object + itp = extrapolate(interpolate(data, interp_type), extrapolation_bc); + function interpolated(params)#, out = similar(data)) + fct_tup_noparams(ci) = fct_tup(ci, params) + #@show Interpolations.adapt(gpu_or_cpu(1), itp) + return apply_transform(fct_tup_noparams, data, itp); + end + return interpolated +end + +""" + get_function_svec(data::AbstractArray, fct_hom::Function; super_sampling=1, extrapolation_bc=Flat(), interp_type=Interpolations.BSpline(Linear())) + +returns a function `interpolated(params)` which generates a transformed version of the original data parameterized by transform parameters. + +# Arguments +`data`: The data to represent by the function `dat` +`fct_hom`: The function to apply to the data +`super_sampling`: The factor by which the data is internally represented as a supersampled version (Fourier-based upsampling, see `FourierTools.resample`) +`extrapolation_bc`: The extrapolation boundary condition to select for values outside the range. + By default the value 0.0 is used. Other options are `Flat()`, or `Line()`, See the package `Interpolation` for details. +`interp_type`: The type of interpolation to use. See the package `Interpolation` for details. +""" +function get_function_svec(data::AbstractArray{T}, fct_hom::Function; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) where T + # building the extraplation + interpolation object + itp = extrapolate(interpolate(data, interp_type), extrapolation_bc); + function interpolated(params::SVector) #, out = similar(data)) + fct_hom_noparams(c) = fct_hom(c, params) + return apply_transform_homogen(fct_hom_noparams, data, itp)#, out); + # return out; + end + return interpolated +end + +""" + get_function_affine(data::AbstractArray; super_sampling=1, extrapolation_bc=Flat(), interp_type=Interpolations.BSpline(Linear())) + +returns a function `interpolated()` which generates a transformed version of the original data parameterized by transform parameters or by transformation matrix. +This is useful for fitting with a function which is itself defined by measured data. +The returned function supports two ways to be used, with an affine transform matrix `matrix_c` as in input or with a vector `p` of parameters. + + +# Arguments +`data`: The data to represent by the function `dat` +`super_sampling`: The factor by which the data is internally represented as a supersampled version (Fourier-based upsampling, see `FourierTools.resample`) +`extrapolation_bc`: The extrapolation boundary condition to select for values outside the range. + By default the value 0.0 is used. Other options are `Flat()`, or `Line()`, See the package `Interpolation` for details. +`interp_type`: The type of interpolation to use. See the package `Interpolation` for details. + +# Example + +```julia +julia> dat1 = reshape(1:16,(4,4)) +4×4 reshape(::UnitRange{Int64}, 4, 4) with eltype Int64: + 1 5 9 13 + 2 6 10 14 + 3 7 11 15 + 4 8 12 16 + +julia> affine_func = get_function_affine(Float32.(dat1)); + +julia> homogeneous_transform = [1 0 -1; 0 1 1; 0 0 1] # translates by [1,1] +3×3 Matrix{Int64}: + 1 0 -1 + 0 1 1 + 0 0 1 + +julia> affine_func(SMatrix{3,3}(homogeneous_transform)) +4×4 Matrix{Float32}: + 0.0 0.0 0.0 0.0 + 5.0 9.0 13.0 0.0 + 6.0 10.0 14.0 0.0 + 7.0 11.0 15.0 0.0 +``` +""" +function get_function_affine(data::AbstractArray{T}; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) where T + #new_size = super_sampling.*size(data) + #upsampled = fftshift(resample(ifftshift(data), new_size)) + + # building the extraplation + interpolation object + itp = extrapolate(interpolate(data, interp_type), extrapolation_bc); + + function interpolated(matrix_c::SMatrix{T}) where T #, out = similar(data)) where T1 + return apply_transform_affine(matrix_c, data, itp)# , out); + # return out; + end + + function interpolated(p::AbstractVector{T1}) where {T1} #, out = similar(data)) where T1 + x_cen, y_cen = (size(data) .÷ 2.0 .+1) + # x_cen_up, y_cen_up = (size(upsampled) .÷ 2.0 .+ 1.0) + + # creating the matrices of rotation, shear, scale, and shift + rot_mat = @SMatrix T[cos(p[7]) -1.0*sin(p[7]) 0.0; sin(p[7]) cos(p[7]) 0.0; 0.0 0.0 1.0]; + shear_mat = @SMatrix T[1.0 p[5] 0.0; p[6] 1.0 0.0; 0.0 0.0 1.0]; + scale_mat = @SMatrix T[1/p[3] 0.0 0.0; 0.0 1/p[4] 0.0; 0.0 0.0 1.0]; + shift_mat = @SMatrix T[1.0 0.0 -1*p[1]; 0.0 1.0 -1*p[2]; 0.0 0.0 1.0]; + t_to_origin = @SMatrix T[1.0 0.0 1*x_cen; 0.0 1.0 y_cen; 0.0 0.0 1.0]; + t_to_center = @SMatrix T[1.0 0.0 -1.0*x_cen; 0.0 1.0 -1.0*y_cen; 0.0 0.0 1.0]; + # t_orig_upsampled = SMatrix{3, 3}(T[1.0 0.0 -1.0*x_cen_up; 0.0 1.0 -1.0*y_cen_up; 0.0 0.0 1.0]); + + # building the overall transformation matrix + matrix_c = t_to_origin * scale_mat * rot_mat * shear_mat *shift_mat * t_to_center + + return apply_transform_affine(matrix_c, data, itp) #, out); # do not call interolated here for type stability reasons + end + + + function interpolated(p::NTuple{N, T}) where {N, T} #, out = similar(data)) where T1 + x_cen, y_cen = (size(data) .÷ 2.0 .+1) + # x_cen_up, y_cen_up = (size(upsampled) .÷ 2.0 .+ 1.0) + + # creating the matrices of rotation, shear, scale, and shift + rot_mat = @SMatrix T[cos(p[7]) -1.0*sin(p[7]) 0.0; sin(p[7]) cos(p[7]) 0.0; 0.0 0.0 1.0]; + shear_mat = @SMatrix T[1.0 p[5] 0.0; p[6] 1.0 0.0; 0.0 0.0 1.0]; + scale_mat = @SMatrix T[1/p[3] 0.0 0.0; 0.0 1/p[4] 0.0; 0.0 0.0 1.0]; + shift_mat = @SMatrix T[1.0 0.0 -1*p[1]; 0.0 1.0 -1*p[2]; 0.0 0.0 1.0]; + t_to_origin = @SMatrix T[1.0 0.0 1*x_cen; 0.0 1.0 y_cen; 0.0 0.0 1.0]; + t_to_center = @SMatrix T[1.0 0.0 -1.0*x_cen; 0.0 1.0 -1.0*y_cen; 0.0 0.0 1.0]; + # t_orig_upsampled = SMatrix{3, 3}(T[1.0 0.0 -1.0*x_cen_up; 0.0 1.0 -1.0*y_cen_up; 0.0 0.0 1.0]); + + # building the overall transformation matrix + matrix_c = t_to_origin * scale_mat * rot_mat * shear_mat *shift_mat * t_to_center + + return apply_transform_affine(matrix_c, data, itp) #, out); # do not call interolated here for type stability reasons + end + + return interpolated +end + + +""" + get_function_poly(data::AbstractArray, order; super_sampling=1, extrapolation_bc=Flat(), interp_type=Interpolations.BSpline(Linear())) + +returns a function `interpolated(p, [out])` which generates a transformed version of the original data parameterized by transform parameters. +This is useful for fitting with a function which is itself defined by measured data. +The returned function supports two ways to be used, with an affine transform matrix `p` as in input or with a vector `p` of parameters. +The optional argument `out` can be used to store the result of the transformation. + + +# Arguments +`data`: The data to represent by the function `dat` +`order`: The order of the polynomial to use for the transformation +`super_sampling`: The factor by which the data is internally represented as a supersampled version (Fourier-based upsampling, see `FourierTools.resample`) +`extrapolation_bc`: The extrapolation boundary condition to select for values outside the range. + By default the value 0.0 is used. Other options are `Flat()`, or `Line()`, See the package `Interpolation` for details. +`interp_type`: The type of interpolation to use. See the package `Interpolation` for details. +""" +function get_function_poly(data::AbstractArray{T, N}, order; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) where {T, N} + pm = get_multi_poly(Val(ndims(data)), Val(order)) + return get_function_tuple(data, pm; super_sampling= super_sampling, extrapolation_bc=extrapolation_bc, interp_type=interp_type); +end + +""" + get_interpolated_function(data::AbstractArray, ::Type{AffineMode}; super_sampling=2, extrapolation_bc=Flat(), interp_type=Interpolations.BSpline(Linear())) + +returns a function `interpolated(p)` which generates a transformed version of the original data parameterized by transform parameters. +This is useful for fitting with a function which is itself defined by measured data. +The returned function supports two ways to be used, with an affine transform matrix `p` as in input or with a vector or tuple `p` of parameters. + +# Arguments +`data`: The data to represent by the function `dat` +`AffineMode`: The transformation mode to use +`super_sampling`: The factor by which the data is internally represented as a supersampled version (Fourier-based upsampling, see `FourierTools.resample`) +`extrapolation_bc`: The extrapolation boundary condition to select for values outside the range. + By default the value 0.0 is used. Other options are `Flat()`, or `Line()`, See the package `Interpolation` for details. +`interp_type`: The type of interpolation to use. See the package `Interpolation` for details. + +# Returns +A function `interpolated(p)` which generates a transformed version of the original data parameterized by transform parameters +""" +function get_interpolated_function(data::AbstractArray{T, N}, ::Type{AffineMode}; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) where {T, N} + return get_function_affine(data; super_sampling=super_sampling, extrapolation_bc=extrapolation_bc, interp_type=interp_type) +end + +function get_interpolated_function(data::AbstractArray{T, N}; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) where {T, N} + @warn "No transformation mode provided. `AffineMode` is used as default" + return get_interpolated_function(data, AffineMode; super_sampling=super_sampling, extrapolation_bc=extrapolation_bc, interp_type=interp_type) +end + +""" + get_interpolated_function(data::AbstractArray, ::Type{PolynomialMode}, order=nothing; super_sampling=2, extrapolation_bc=Flat(), interp_type=Interpolations.BSpline(Linear())) + +returns a function `interpolated(p)` which generates a transformed version of the original data parameterized by transform parameters. +This is useful for fitting with a function which is itself defined by measured data. +The returned function supports polynomial transformations of the data. + +# Arguments +`data`: The data to represent by the function `dat` +`PolynomialMode`: The transformation mode to use +`order`: The order of the polynomial to use for the transformation +`super_sampling`: The factor by which the data is internally represented as a supersampled version (Fourier-based upsampling, see `FourierTools.resample`) +`extrapolation_bc`: The extrapolation boundary condition to select for values outside the range. + By default the value 0.0 is used. Other options are `Flat()`, or `Line()`, See the package `Interpolation` for details. +`interp_type`: The type of interpolation to use. See the package `Interpolation` for details. + +# Returns +A function `interpolated(p)` which generates a transformed version of the original data parameterized by polynomial transform parameters +""" +function get_interpolated_function(data::AbstractArray{T, N}, ::Type{PolynomialMode}, order=nothing; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) where {T, N} + #TODO from CuArray to CuArray + if isnothing(order) + error("Providing the order of the transformation polynomial is mandatory for the `PolynomialMode`") + end + return get_function_poly(data, order; super_sampling=super_sampling, extrapolation_bc=extrapolation_bc, interp_type=interp_type) +end \ No newline at end of file diff --git a/test/Aqua.jl b/test/Aqua.jl new file mode 100644 index 0000000..e7b0cce --- /dev/null +++ b/test/Aqua.jl @@ -0,0 +1,6 @@ +using Aqua + +Aqua.test_all( + DataToFunctions, + unbound_args=false, + ) \ No newline at end of file diff --git a/test/runtests.jl b/test/runtests.jl index ddd940c..1735829 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,32 +1,36 @@ using Test -# using Zygote +using Zygote using DataToFunctions -@testset "get_function" begin +include("Aqua.jl") + +@testset "get_function_affine" begin data = rand(40,41) - for supersamp = 1:5 - f = get_function(data; super_sampling=supersamp); - @test f((0.0,0.0),(1.0,1.0)) ≈ data - end + f = get_function_affine(data); + @test f([0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0]) ≈ data + end -@testset "gradient" begin +@testset "loss" begin data = rand(11,10) - f = get_function(data; super_sampling=2); - loss(p,z) = sum(abs2.(f(p, z) .- data)) - @test loss((0.0,0.0),(1.0,1.0)) < 1e-20 - @test loss((0.0,0.001),(1.0,1.0)) > 1e-20 - @test loss((0.0,0.0),(1.0001,1.0)) > 1e-20 + f = get_function_affine(data; super_sampling=2); + loss(p) = sum(abs2.(f(p) .- data)) + @test loss([0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0]) < 1e-20 + @test loss([0.001, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0]) > 1e-20 + @test loss([0.0, 0.0, 1.001, 1.0, 0.0, 0.0, 0.0]) > 1e-20 +end +@testset "gradient" begin + data = rand(11,10) + f = get_function_affine(data; super_sampling=2); + loss(p) = sum(abs2.(f(p) .- data)) + st_vals = [0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0] # throws an error... - # Zygote.gradient(loss, (0.0,0.0), (1.0,1.0)) + @test Zygote.gradient(loss, st_vals)[1] ≈ zeros(7) end @testset "keep center" begin data = ones(5,4); data[3,3] = 5.0; - f = get_function(data; super_sampling=5); - @test f((0.0,0.0),(2.0,2.0))[3,3] ≈ 5.0 - - # throws an error... - # Zygote.gradient(loss, (0.0,0.0), (1.0,1.0)) + f = get_function(data); + @test f([0.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0])[3,3] ≈ 5.0 end