From 4e117fb99beb6c3440a102f5fb8575b71113a167 Mon Sep 17 00:00:00 2001 From: hzarei4 Date: Tue, 23 Jul 2024 13:26:18 +0200 Subject: [PATCH] modified to use EvalMultiPoly.jl --- Project.toml | 1 + examples/Project.toml | 2 + examples/deform_testimage.jl | 96 ++++++++++--------- examples/polynomial_apply.jl | 25 +++-- src/DataToFunctions.jl | 3 +- src/polynomials.jl | 173 ----------------------------------- src/transformation_types.jl | 4 + src/transformators.jl | 83 ++++++++++++++++- src/utils.jl | 18 ---- 9 files changed, 151 insertions(+), 254 deletions(-) delete mode 100644 src/polynomials.jl create mode 100644 src/transformation_types.jl delete mode 100644 src/utils.jl diff --git a/Project.toml b/Project.toml index 9f7abd6..eb883bc 100644 --- a/Project.toml +++ b/Project.toml @@ -4,6 +4,7 @@ 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" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" diff --git a/examples/Project.toml b/examples/Project.toml index a25d76a..5a2e1dd 100644 --- a/examples/Project.toml +++ b/examples/Project.toml @@ -4,6 +4,8 @@ CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" 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" FindShift = "643ec891-bf64-479f-8088-26ff5ce1b396" ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210" FourierTools = "b18b359b-aebc-45ac-a139-9c0ccbb2871e" diff --git a/examples/deform_testimage.jl b/examples/deform_testimage.jl index d2ac90f..274c0ae 100644 --- a/examples/deform_testimage.jl +++ b/examples/deform_testimage.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.19.43 +# v0.19.42 using Markdown using InteractiveUtils @@ -20,11 +20,14 @@ 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 -# ╔═╡ 4af0c13d-fc42-4fe7-97e6-2248e36b63e2 -# Pkg.add("PlutoUI") +# ╔═╡ 1c744bec-f085-4812-ab1a-40a32c2ac176 +import PlutoUI: combine # ╔═╡ 5ac1123d-5df3-4c9d-aff1-ffe91d931497 data = Float32.(testimage("resolution_test_512")) @@ -33,7 +36,7 @@ data = Float32.(testimage("resolution_test_512")) simshow(data) # ╔═╡ 2fce0208-732f-4259-a47d-7f78921bfd87 -f = get_function(data, super_sampling=1) +f = get_interpolated_function(data, AffineMode, super_sampling=1) # ╔═╡ dd535378-3f2a-4429-914c-8b7608b99706 @bind shift_x Slider(-100:0.02:100, default=0) @@ -48,57 +51,68 @@ f = get_function(data, super_sampling=1) @bind zoom_y Slider(0.2:0.02:4, default=1) # ╔═╡ 64b64d3b-092f-4553-916d-f7db1fdfa428 -f((shift_x, shift_y), (1/zoom_x, 1/zoom_y)) +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))) +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) -# ╔═╡ 1bed34fb-b29a-4042-a493-4835fdb69a9d -g =DataToFunctions.get_function_affine(data, super_sampling=1) +# ╔═╡ 6676ba35-3efd-49f5-9819-411ad8f8a95c +md""" +# Polynomial transformations +""" -# ╔═╡ c287fa80-426b-11ef-125e-5fda207e605c -# ╠═╡ disabled = true -#=╠═╡ -g() - ╠═╡ =# - -# ╔═╡ 87be45c0-8b2e-4d49-abd1-a274b3c1815e -h = get_function_poly(data, 1) +# ╔═╡ 13461a95-95ea-4bad-8673-e94e06776254 +md""" +First order polynomial transformation which is as follows: -# ╔═╡ 473d635e-d06d-4cdb-991f-22c82a06b491 -@bind c1 Slider(-20f0:0.05f0:20f0, default=0) +``x^{\prime} = c_{1} + c_{2}{x}^{1} + c_{3}{y}^{1}`` -# ╔═╡ 7fb3ef17-d0a0-4919-b4d5-8e66b4a0fe60 -@bind c2 Slider(0.2f0:0.05f0:2f0, default=1) +``y^{\prime} = c_{4} + c_{5}{x}^{1} + c_{6}{y}^{1}`` +""" -# ╔═╡ 200bab4d-b444-47c6-b4d8-d5d3c450e5f9 -@bind c3 Slider(-2f0:0.05f0:2f0, default=0) - -# ╔═╡ f17402a6-ba63-44e9-8e22-da027b07ffc3 -@bind c4 Slider(-20f0:0.05f0:20f0, default=0) - -# ╔═╡ d913278e-1658-4bee-9e12-ad778e530c1b -@bind c5 Slider(-2f0:0.05f0:2f0, default=0) - -# ╔═╡ cee301da-2b3e-432f-9551-0fe83bf6c8ec -@bind c6 Slider(0.2f0:0.05f0:2f0, default=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""" + #### Transfrorm 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((c1, c2, c3, c4, c5, c6)), cmap=:turbo)#,0f0,0f0,0f0,0f0,0f0,0f0,0f0,0f0,0f0))) +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((c1,c2,c3,c4,c5,c6))) +maximum(h((c.c1,c.c2,c.c3,c.c4,c.c5,c.c6))) # ╔═╡ 6584aacc-440e-457b-bf52-83f8db40c999 -h((c1,c2,c3,c4,c5,c6)) +h((c.c1,c.c2,c.c3,c.c4,c.c5,c.c6)) # ╔═╡ 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 @@ -110,15 +124,11 @@ h((c1,c2,c3,c4,c5,c6)) # ╠═64b64d3b-092f-4553-916d-f7db1fdfa428 # ╠═c3fe6b25-f4c0-4a80-9d15-9ee30136d43b # ╠═ff143f0d-b070-4220-8fa6-0b5a93a56303 -# ╠═1bed34fb-b29a-4042-a493-4835fdb69a9d -# ╠═c287fa80-426b-11ef-125e-5fda207e605c +# ╟─6676ba35-3efd-49f5-9819-411ad8f8a95c +# ╟─13461a95-95ea-4bad-8673-e94e06776254 # ╠═87be45c0-8b2e-4d49-abd1-a274b3c1815e -# ╠═473d635e-d06d-4cdb-991f-22c82a06b491 -# ╠═7fb3ef17-d0a0-4919-b4d5-8e66b4a0fe60 -# ╠═200bab4d-b444-47c6-b4d8-d5d3c450e5f9 -# ╠═f17402a6-ba63-44e9-8e22-da027b07ffc3 -# ╠═d913278e-1658-4bee-9e12-ad778e530c1b -# ╠═cee301da-2b3e-432f-9551-0fe83bf6c8ec -# ╠═74228a9d-6cc2-4aaf-97da-67f32670341e +# ╟─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 diff --git a/examples/polynomial_apply.jl b/examples/polynomial_apply.jl index 39b565d..501502b 100644 --- a/examples/polynomial_apply.jl +++ b/examples/polynomial_apply.jl @@ -29,25 +29,22 @@ function poly_test(;sz=64, dtype=Float64, n_photons=1000) #sample_data = p_psf ./ maximum(p_psf) # normalizing the sample data - sample_data = p_psf ./ maximum(p_psf) + sample_data = p_psf ./ maximum(p_psf) * n_photons #sample_data = make_grid(); - f_1 = get_function_poly(Float64.(sample_data), 1); # get_function_affine(sample_data); - true_vals = dtype[2.1, 1.05, 0.02, 1.5, 0.05, 1.02] # dtype[2.0, 1.0, 1.01, 1.0, 0.0, 0.0, 0.0] + 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(Tuple(true_vals)) - p_img = n_photons .* (dat2 ./ maximum(dat2)) - n_img = dtype.(poisson(Float64.(p_img))) + dat2 = f_1(true_vals) - s_data = Float64.(sample_data .* n_photons) #sample_data = dtype.(TestImages.shepp_logan(sz)) - f = get_function_poly(Float64.(s_data), 1); + f = get_interpolated_function(Float64.(sample_data), PolynomialMode, 1); - loss_m(p1::AbstractVector) = sum(abs2.(f(Tuple(p1)) .- n_img)) - st_vals = dtype[2.1, 1.01, 0.01, 1.5, 0.01, 0.9] #ones(Float64, 6)./10 + 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)) @@ -55,11 +52,11 @@ function poly_test(;sz=64, dtype=Float64, n_photons=1000) function g!(G, x) # (G, x) - G .= gradient(loss_m, x)[1] + G .= gradient(loss_p, x)[1] end - od = OnceDifferentiable(loss_m, g!, st_vals) + od = OnceDifferentiable(loss_p, g!, st_vals) res = optimize( - od, + loss_p, st_vals, #Newton(), BFGS(; initial_stepnorm = 1e-2),#; linesearch=LineSearches.BackTracking(order=2)), @@ -67,7 +64,7 @@ function poly_test(;sz=64, dtype=Float64, n_photons=1000) #lower, upper, #init_x, #Fminbox(inner_optimizer), - Optim.Options(store_trace = true, extended_trace = true, iterations=5000, g_tol=1e-3), + #Optim.Options(store_trace = true, extended_trace = true, iterations=5000, g_tol=1e-3), autodiff = :forward ) diff --git a/src/DataToFunctions.jl b/src/DataToFunctions.jl index 93857cd..77019a6 100644 --- a/src/DataToFunctions.jl +++ b/src/DataToFunctions.jl @@ -1,7 +1,6 @@ module DataToFunctions -include("utils.jl") -include("polynomials.jl") +include("transformation_types.jl") include("transformators.jl") end # module DataToFunctions \ No newline at end of file diff --git a/src/polynomials.jl b/src/polynomials.jl deleted file mode 100644 index 969dc3f..0000000 --- a/src/polynomials.jl +++ /dev/null @@ -1,173 +0,0 @@ -export get_multi_poly, get_num_poly_vars, get_num_multipoly_vars -export polynomial - -""" - polynomial(::Val{0}, ::T1, c::T2, ::Val{cstart}=Val(1), ::Val{numvars}=Val(1)) where {NV, TS, T1 <: NTuple{NV, Integer}, T2 <: NTuple{TS, Float32}, cstart, numvars} - -Create a polynomial of order 0 with numvars variables. - -returned is a function that takes a tuple of variables and a tuple of coefficients and returns the value of the polynomial - and the number of coefficients required (here 1). -""" -function polynomial(::Val{0}, ::T1, c::T2, ::Val{cstart}=Val(1), ::Val{numvars}=Val(1)) where {NV, TS, T1 <: NTuple{NV, Integer}, T2 <: NTuple{TS, Float32}, cstart, numvars} - # println("c: $(c) $(length(c))"); - return c[cstart] -end #, (t,c) -> ntuple(n->c[1], Val(numvars)) - -""" - polynomial(::Val{N}, t::T1, c::T2, ::Val{cstart}=Val(1), ::Val{numvars}=Val(length(t)))::Float32 where {N, NV, TS, T1 <: NTuple{NV, Integer}, T2 <: NTuple{TS, Float32}, cstart, numvars} - -Represents a polynomial of order N with numvars variables (also implicitely defined via the length of the NTuple `t`). -Note that `numvars` is needed for the internal workings of the polynomial generator, but notmally not by the user. - -E.g. to represent a polynomial of order 1 with 2 variables, the coefficients are ordered as follows: -c = (c0, c1, c2) where the polynomial is: c0 + c1*x + c2*y -or for a polynomial of order 2 with 2 variables: -c = (c0, c1, c2, c3, c4, c5) where the polynomial is c0 + c1*x + c2*x^2 + c3*y + c4*x*y + c5*y^2 -Note that the coefficients are ordered not by the multiples in which they appear in the polynomial, but by the order of the variables. - -Example: -```jldoctest ->julia polynomial(Val(2), (2, 20), (1f0, 1f0, 1f0, 1f0, 1f0, 1f0)) -467.0f0 ->julia cs = Tuple(Float32.(collect(1:27))) - (1.0f0, 2.0f0, 3.0f0, 4.0f0, 5.0f0, 6.0f0, 7.0f0, 8.0f0, 9.0f0, 10.0f0, 11.0f0, 12.0f0, 13.0f0, 14.0f0, 15.0f0, 16.0f0, 17.0f0, 18.0f0, 19.0f0, 20.0f0, 21.0f0, 22.0f0, 23.0f0, 24.0f0, 25.0f0, 26.0f0, 27.0f0) ->julia res = zeros(Float32, 200,200) ->julia @time res .= polynomial.(Ref(Val(2)), Tuple.(CartesianIndices((200,200))), Ref(cs)); - 0.054017 seconds (147.45 k allocations: 10.168 MiB, 99.75% compilation time) ->julia @time res .= polynomial.(Ref(Val(2)), Tuple.(CartesianIndices((200,200))), Ref(cs)); - 0.000089 seconds (3 allocations: 184 bytes) -``` -""" -@generated function polynomial(::Val{N}, t::T1, c::T2, ::Val{cstart}=Val(1), ::Val{numvars}=Val(length(t)))::Float32 where {N, NV, TS, T1 <: NTuple{NV, Integer}, T2 <: NTuple{TS, Float32}, cstart, numvars} - quote - c_start = cstart - res = c[c_start] - c_start += 1 - # iterate through the polynomial variables: (but this leads to dynamic memory allocation!) - # for n = 1:numvars # eachindex(t) # 1:length(t) # eachindex(t) # 1:length(t) - # # c_end = c_start + get_num_poly_vars(Val(n), Val(N-1)) - 1 - # res += t[n] * polynomial(Val(N-1), t, c, Val(c_start), Val(n)) - # c_start += get_num_poly_vars(Val(n), Val(N-1)) # c_end + 1 - # end - # # does not work: - # @macroexpand Base.Cartesian.@nexprs 4 n -> begin - # if (numvars >= n) - # res += t[n] * polynomial(Val(N-1), t, c, Val(c_start), Val(n)) - # c_start += get_num_poly_vars(Val(n), Val(N-1)) # c_end + 1 - # end - # end - - # @generated function myvarsum(t, c, res, c_start, ::Val{numvars}) where {N} - # quote - # c_start = cstart - # res = c[c_start] - # c_start += 1 - # Base.Cartesian.@nexprs $N i A n -> begin - # res += t[n] * polynomial(Val(N-1), t, c, Val(c_start), Val(n)) - # c_start += get_num_poly_vars(Val(n), Val(N-1)) # c_end + 1 - # end - # res, c_start - # end - # end - - Base.Cartesian.@nexprs $numvars n -> begin - res += t[n] * polynomial(Val(N-1), t, c, Val(c_start), Val(n)) - c_start += get_num_poly_vars(Val(n), Val(N-1)) # c_end + 1 - end - # if numvars > 5 - # error("Only up to 5 dimensions are currently supported for polynomials") - # end - # this simply unrolls the loop by hand (up to 4D input variables): - # if numvars >= 1 - # res += t[1] * polynomial(Val(N-1), t, c, Val(c_start), Val(1)) - # c_start += get_num_poly_vars(Val(1), Val(N-1)) # c_end + 1 - # end - # if numvars >= 2 - # res += t[2] * polynomial(Val(N-1), t, c, Val(c_start), Val(2)) - # c_start += get_num_poly_vars(Val(2), Val(N-1)) # c_end + 1 - # end - # if numvars >= 3 - # res += t[3] * polynomial(Val(N-1), t, c, Val(c_start), Val(3)) - # c_start += get_num_poly_vars(Val(3), Val(N-1)) # c_end + 1 - # end - # if numvars >= 4 - # res += t[4] * polynomial(Val(N-1), t, c, Val(c_start), Val(4)) - # c_start += get_num_poly_vars(Val(4), Val(N-1)) # c_end + 1 - # end - # if numvars >= 5 - # error("Only up to 4 dimensions are currently supported for polynomials") - # end - return res - end -end - - -# function unroll_loop(res::Float32, ::Var{0}, ::Var{numvars}) where{n} -# return 0f0; -# end - -# function unroll_loop(res::Float32, t, ::Var{c_start},::Var({n}, ::Var{numvars}) where{c_start, n, numvars} -# res += t[n] * polynomial(Val(N-1), t, c, Val(c_start), Val(n)) -# c_start += get_num_poly_vars(Val(n), Val(N-1)) # c_end + 1 -# return unroll_loop(res::Float32, ::Var({n}, ::Var{numvars}) + -# end - - -function get_multi_poly(::Val{numvars}, ::Val{N}) where {numvars, N} - # cs_per_comp = ((numvars+1)^N) - @info "Creating polynomials with $(numvars) variables of order , $(N). Required constants: $(numvars*((numvars+1)^N))" - p = (t,c) -> polynomial(Val(N), Tuple.(t), c) - # return p - function mpol(t,c)#::NTuple{numvars, T} where T - return ntuple(n->p(t, split_tuple(c, Val(numvars))[n]), Val(numvars)) - - # println("N: $(N), c: $(c) $(length(c))"); - # return Tuple(p(t, c[1+(n-1)*((numvars+1)^N):n*((numvars+1)^N)]) for n=1:numvars) - # return ntuple(n->p(t, c[1+(n-1)*cs_per_comp:n*cs_per_comp]), Val(numvars)) - # return ntuple(n->p(t, c[1+(n-1)*((numvars+1)^N):n*((numvars+1)^N)]), Val(numvars)) - # return (p(t, c[1+(1-1)*((numvars+1)^N):1*((numvars+1)^N)]), p(t, c[1+(2-1)*((numvars+1)^N):2*((numvars+1)^N)])) - end - return mpol # (t, c)->ntuple(n->p(t, c[1+(n-1)*((numvars+1)^N):n*((numvars+1)^N)]), Val(numvars)) - # return (t, c) -> Tuple(p(t,c[1+(n-1)*((numvars+1)^N):n*((numvars+1)^N)]) for n=1:numvars) -end - -function get_num_poly_vars(::Val{numvars}, ::Val{N}) where {numvars, N} - return binomial(numvars+N, N) -end - -function get_num_multipoly_vars(::Val{numvars}, ::Val{N}) where {numvars, N} - return numvars*get_num_poly_vars(Val(numvars), Val(N)) -end - -function test_poly_allocations() - # get_num_poly_vars(Val(2), Val(3)) # 10 indices - # p = get_polynomial(Val(2), Val(3)) - # @time p.(Tuple.(CartesianIndices((200,200))),Ref((1.1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27))); - - polynomial(Val(2), (2, 20), (1f0, 1f0, 1f0, 1f0, 1f0, 1f0)) == 467 - - cids = Tuple.(CartesianIndices((200,200))) - # cfds = map((t)->Tuple(Float32.([t...])), cids) - cs = Tuple(Float32.(collect(1:27))) - res = zeros(Float32, 200,200) - get_num_poly_vars(Val(2), Val(2)) # 6 indices required - @time res .= polynomial.(Ref(Val(2)), cids, Ref(cs)); # 2 orders, two variables - # 0.000122 seconds (3 allocations: 168 bytes) - - get_num_poly_vars(Val(3), Val(2)) # 10 indices required - @time res .= polynomial.(Ref(Val(3)), cids, Ref(cs)); # 3 orders, two variables - # 0.003710 seconds (240.00 k allocations: 13.428 MiB) - - get_num_poly_vars(Val(4), Val(2)) # 15 indices required - @time res .= polynomial.(Ref(Val(4)), cids, Ref(cs)); # 2 orders, two variables - # 0.008149 seconds (480.00 k allocations: 26.856 MiB) - - get_num_poly_vars(Val(5), Val(2)) # 21 indices required - @time res .= polynomial.(Ref(Val(5)), cids, Ref(cs)); # 2 orders, two variables - #0.014299 seconds (1.08 M allocations: 60.425 MiB, 23.18% gc time) - - @time polynomial.(Ref(Val(0)), cids, Ref(cs)); - # 0.000076 seconds (5 allocations: 156.461 KiB) - -end 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 index 820d91f..9a72e7b 100644 --- a/src/transformators.jl +++ b/src/transformators.jl @@ -1,11 +1,14 @@ using Interpolations using FourierTools using StaticArrays +using EvalMultiPoly -export get_function, get_function_tuple, get_function_svec, get_function_affine, get_function_poly +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())) @@ -21,7 +24,7 @@ This is useful for fitting with a function which is itself defined by measured d """ -function get_function(data::AbstractArray; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) +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 @@ -177,7 +180,7 @@ function apply_transform(coord_transf_func::Function, data::AbstractArray{T}, it out[it] = idx_apply(itp, coord_transf_func(it)) end """ - return map((it) -> idx_apply(itp, coord_transf_func(it)), (CartesianIndices(data))) + return map((it) -> idx_apply(itp, coord_transf_func(it)), Tuple.(CartesianIndices(data))) # out .= idx_apply.(Ref(itp), coord_transf_func.(CartesianIndices(data))); # return idx_apply.(Ref(itp), coord_transf_func.(CartesianIndices(data))); end @@ -311,7 +314,26 @@ function get_function_affine(data::AbstractArray{T}; super_sampling=2, extrapola 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 - #return out; + 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 @@ -338,4 +360,57 @@ The optional argument `out` can be used to store the result of the transformatio function get_function_poly(data::AbstractArray{T}, order; super_sampling=2, extrapolation_bc=zero(eltype(data)), interp_type=Interpolations.BSpline(Linear())) where T 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} + 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/src/utils.jl b/src/utils.jl deleted file mode 100644 index b04e514..0000000 --- a/src/utils.jl +++ /dev/null @@ -1,18 +0,0 @@ -export split_tuple - -""" - split_tuple(t::NTuple{S,T},::Val{numvars}) where {S,T,numvars} - -Split a tuple into `numvars` parts packed into a tuple of tuples. The tuple `t` is assumed to have a length that is a multiple of `numvars`. - -Example: -```juliadoc -julia> t = (1,2,3,4,5,6) -julia> split_tuple(t, Val{2}) -((1,2), (3,4), (5,6)) -``` -""" -function split_tuple(t::NTuple{S,T},::Val{numvars}) where {S,T,numvars} - return ntuple(n->t[1+(n-1)*(S÷numvars):n*(S÷numvars)], Val(numvars)) -end -