diff --git a/Project.toml b/Project.toml index 6831fd5..9f45ae7 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "LocalAnisotropies" uuid = "cb698b3a-2f1d-48b4-a1e9-ddc0b8f648e2" authors = ["Rafael Caixeta and contributors"] -version = "0.6.2" +version = "0.7.0" [deps] Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7" diff --git a/ext/LocalAnisotropiesMakieExt.jl b/ext/LocalAnisotropiesMakieExt.jl index e2ba2fe..cd48b54 100644 --- a/ext/LocalAnisotropiesMakieExt.jl +++ b/ext/LocalAnisotropiesMakieExt.jl @@ -17,13 +17,15 @@ Makie.plottype(::LocalAnisotropy, ::LocalAnisotropies.SpatialData) = function pre_ellipsoid(lpars, geom, f) # pre process coords = LocalAnisotropies.coords_(geom) + mag = lpars.magnitude + mx = maximum D3 = size(coords, 1) >= 3 x = coords[1, :] y = coords[2, :] z = D3 ? coords[3, :] : map(i -> 0, x) s = - D3 ? [f * Makie.Vec3f(i[1], i[2], i[3]) for i in eachcol(lpars.magnitude)] : - [f * Makie.Vec3f(i[1], i[2], i[2]) for i in eachcol(lpars.magnitude)] + D3 ? [f * Makie.Vec3f(i[1]/mx(i), i[2]/mx(i), i[3]/mx(i)) for i in eachcol(mag)] : + [f * Makie.Vec3f(i[1]/mx(i), i[2]/mx(i), i[2]/mx(i)) for i in eachcol(mag)] q = map(q -> Makie.Quaternion(q.q1, q.q2, q.q3, q.q0), lpars.rotation) x, y, z, s, q diff --git a/src/estimators/kriging.jl b/src/estimators/kriging.jl index 47dd5e3..a40f9ba 100644 --- a/src/estimators/kriging.jl +++ b/src/estimators/kriging.jl @@ -40,7 +40,12 @@ krig_estimator(model::KC_SKModel, localpar = nothing) = SimpleKriging(model.γ, neighs_localaniso(model::MWModels, m) = nothing neighs_localaniso(model::KCModels, m) = view(model.hdlocalaniso, m) -function neighs_localaniso(localaniso::LocalAnisotropy, dom, neigh; method::Symbol=:KernelConvolution) +function neighs_localaniso( + localaniso::LocalAnisotropy, + dom, + neigh; + method::Symbol = :KernelConvolution, +) method != :KernelConvolution ? nothing : qmat(nnpars(localaniso, dom, neigh)) end diff --git a/src/estimators/sgs.jl b/src/estimators/sgs.jl index 4472264..dfa619f 100644 --- a/src/estimators/sgs.jl +++ b/src/estimators/sgs.jl @@ -146,7 +146,7 @@ end function local_probmodel(probmodel, localaniso, hdlocalaniso) isnothing(hdlocalaniso) ? MW_SKModel(localaniso, probmodel.γ, probmodel.μ) : - KC_SKModel(localaniso, probmodel.γ, probmodel.μ, hdlocalaniso) + KC_SKModel(localaniso, probmodel.γ, probmodel.μ, hdlocalaniso) end function Meshes._pboxes(::Type{𝔼{N}}, points) where {N} @@ -154,13 +154,13 @@ function Meshes._pboxes(::Type{𝔼{N}}, points) where {N} ℒ = lentype(p) cmin = fill(typemax(ℒ), N) cmax = fill(typemin(ℒ), N) - + for p in points - c = getfield(coords(p), :coords) - for i in 1:N - cmin[i] = min(c[i], cmin[i]) - cmax[i] = max(c[i], cmax[i]) - end + c = getfield(coords(p), :coords) + for i = 1:N + cmin[i] = min(c[i], cmin[i]) + cmax[i] = max(c[i], cmax[i]) + end end Box(withcrs(p, Tuple(cmin)), withcrs(p, Tuple(cmax))) - end \ No newline at end of file +end diff --git a/src/estimators/variograms.jl b/src/estimators/variograms.jl index 29a1dfc..7d20a81 100644 --- a/src/estimators/variograms.jl +++ b/src/estimators/variograms.jl @@ -16,7 +16,8 @@ end qmat(L::LocalAnisotropy, i::Int) = qmat(localpair(L, i)...) qmat(L::LocalGeoData, i::Int) = qmat(localpair(L, i)...) -qmat(L::Union{LocalAnisotropy,Nothing}) = isnothing(L) ? nothing : [qmat(L, i) for i = 1:nvals(L)] +qmat(L::Union{LocalAnisotropy,Nothing}) = + isnothing(L) ? nothing : [qmat(L, i) for i = 1:nvals(L)] function mw_estimator(μ, γ, localpar) # get local Mahalanobis matrix @@ -69,10 +70,3 @@ function kccov(γ::Variogram, xi, xj, Qi::AbstractMatrix, Qj::AbstractMatrix) Cij = (sill(γl) - γl(xi, xj)) (det(Qi)^0.25) * (det(Qj)^0.25) * (det(Qij)^-0.5) * Cij end - -function localball(radii, rot, convention) - P, Λ = rotmat(radii, rot, convention) - Q = P' * Λ * P - Qs = ustrip.(Q ./ collect(radii)' .^ 2) - MetricBall(1.0, Mahalanobis(Symmetric(Qs))) -end diff --git a/src/parametrization/conventions.jl b/src/parametrization/conventions.jl index f3e6b9f..3cd1fe9 100644 --- a/src/parametrization/conventions.jl +++ b/src/parametrization/conventions.jl @@ -98,6 +98,22 @@ function rotmat( # rotation matrix P = angle_to_dcm(angles..., rule.order)[SOneTo(N), SOneTo(N)] - !rule.extrinsic && (P = P') + rule.extrinsic && (P = P') P, Λ end + +function anisodistance( + semiaxes::AbstractVector, + angles::AbstractVector, + convention = :TaitBryanExtr; + rev = false, +) + P, Λ = rotmat(semiaxes, angles, convention; rev) + Q = P' * Λ * P + Mahalanobis(Symmetric(Q)) +end + +function localball(radii, rot, convention) + metric_ = anisodistance(radii, rot, convention) + MetricBall(1.0, metric_) +end diff --git a/src/parametrization/conversions.jl b/src/parametrization/conversions.jl index 2496597..2d66dab 100644 --- a/src/parametrization/conversions.jl +++ b/src/parametrization/conversions.jl @@ -51,18 +51,21 @@ function localanisotropies( if transf rule = convention isa Symbol ? rules[convention] : convention - rule.main == :y && (ranges = ranges[reverse(1:dim, 1, 2)]) + if rule.main == :y + ranges = ranges[reverse(1:dim, 1, 2)] + rule = RotationRule(rule.order, rule.motion, rule.radian, :x, rule.extrinsic) + end end for i = 1:len - xranges = - [Tables.getcolumn(tab, x)[i] for x in ranges] ./ Tables.getcolumn(tab, ranges[1])[i] + vranges = [Tables.getcolumn(tab, x)[i] for x in ranges] + xranges = vranges ./ maximum(vranges) xrot = [Tables.getcolumn(tab, x)[i] for x in rotation] quat = if isquat Quaternion(xrot) else - P, Λ = rotmat(xranges, xrot, rule; rev = false) + P, Λ = rotmat(xranges, xrot, rule) size(P, 1) == 2 && (P = DCM([P[1, 1] P[1, 2] 0; P[2, 1] P[2, 2] 0; 0 0 1])) dcm_to_quat(P) end @@ -196,7 +199,7 @@ function rotmat2angles(dcm::AbstractMatrix, convention::RotConvention) P = N == 2 ? DCM([P[1, 1] P[1, 2] 0; P[2, 1] P[2, 2] 0; 0 0 1]) : dcm rule = convention isa Symbol ? rules[convention] : convention - !rule.extrinsic && (P = P') + rule.extrinsic && (P = P') preangs = dcm_to_angle(DCM(P), rule.order) angles = [preangs.a1, preangs.a2, preangs.a3] diff --git a/test/runtests.jl b/test/runtests.jl index 426f960..a1e21a1 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,18 +1,66 @@ -using LocalAnisotropies +using Distances using GeoStats +using LocalAnisotropies using Test -import LocalAnisotropies: rotmat +import LocalAnisotropies: rotmat, anisodistance @testset "LocalAnisotropies.jl" begin + d₁ = anisodistance([1.0, 1.0], [0.0]) + d₂ = anisodistance([1.0, 2.0], [0.0]) + @test evaluate(d₁, [1.0, 0.0], [0.0, 0.0]) == evaluate(d₁, [0.0, 1.0], [0.0, 0.0]) + @test evaluate(d₂, [1.0, 0.0], [0.0, 0.0]) != evaluate(d₂, [0.0, 1.0], [0.0, 0.0]) + + d₃ = anisodistance([1.0, 0.5, 0.5], [π / 4, 0.0, 0.0]) + @test evaluate(d₃, [1.0, 1.0, 0.0], [0.0, 0.0, 0.0]) ≈ √2 + @test evaluate(d₃, [-1.0, 1.0, 0.0], [0.0, 0.0, 0.0]) ≈ √8 + + # intrinsic conventions + gslib = anisodistance([50.0, 25.0, 5.0], [30.0, -30.0, 30.0], :GSLIB) + tait = anisodistance([25.0, 50.0, 5.0], [-π / 6, -π / 6, π / 6], :TaitBryanIntr) + euler = anisodistance( + [50.0, 25.0, 5.0], + [-deg2rad(78), -deg2rad(41), -deg2rad(50)], + :EulerIntr, + ) + lpf = anisodistance([50.0, 25.0, 5.0], [78.0, 41.0, 50.0], :Leapfrog) + dm = anisodistance([50.0, 25.0, 5.0], [78.0, 41.0, 50.0], :Datamine) + + @test evaluate(gslib, [1.0, 0.0, 0.0], [0.0, 0.0, 0.0]) ≈ 0.1325707358356285 + @test evaluate(gslib, [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]) ≈ 0.039051248379533283 + @test evaluate(gslib, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈ 0.15132745950421558 + + @test evaluate(gslib, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈ + evaluate(tait, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) + @test evaluate(euler, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈ + evaluate(dm, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) + @test evaluate(euler, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈ + evaluate(lpf, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) + @test evaluate(euler, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) - + evaluate(gslib, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) < 10^-3 + + # extrinsic conventions + xtait = anisodistance([50.0, 25.0, 5.0], [π, 0, π / 2], :TaitBryanExtr) + xeuler = anisodistance([50.0, 25.0, 5.0], [-π / 2, -π / 2, -π / 2], :EulerExtr) + + @test evaluate(xtait, [1.0, 0.0, 0.0], [0.0, 0.0, 0.0]) ≈ 0.20 + @test evaluate(xtait, [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]) ≈ 0.04 + @test evaluate(xtait, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈ 0.02 + + @test evaluate(xtait, [1.0, 0.0, 0.0], [0.0, 0.0, 0.0]) ≈ + evaluate(xeuler, [1.0, 0.0, 0.0], [0.0, 0.0, 0.0]) + @test evaluate(xtait, [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]) ≈ + evaluate(xeuler, [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]) + @test evaluate(xtait, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) ≈ + evaluate(xeuler, [0.0, 0.0, 1.0], [0.0, 0.0, 0.0]) + # convert data to LocalAnisotropy dummy = georef( (az = 1:10, r1 = 1:10, r2 = 1:10), PointSet([(i / 2, (i + 1) / 2) for i = 1:2:19]), ) pars = localanisotropies(dummy, [:az], [:r1, :r2], :GSLIB) - @test round(pars.rotation[1][4], digits = 4) ≈ 0.0087 tabpars_quat = to_table(pars) tabpars_angs = to_table(pars, :GSLIB) @@ -30,9 +78,11 @@ import LocalAnisotropies: rotmat # convert between different rotation conventions angs1 = convertangles([30, 30, 30], :GSLIB, :Datamine) angs2 = convertangles(angs1, :Datamine, :GSLIB) - rmat = rotmat([1, 1, 1], [30, 30, 30], :GSLIB) - @test all(rmat .≈ rotmat([1, 1, 1], angs1, :Datamine)) - @test all(rmat .≈ rotmat([1, 1, 1], angs2, :GSLIB)) + rmat, _ = rotmat([1, 1, 1], [30, 30, 30], :GSLIB) + dmat, _ = rotmat([1, 1, 1], angs1, :Datamine) + gmat, _ = rotmat([1, 1, 1], angs2, :GSLIB) + @test all(rmat .≈ dmat) + @test all(rmat .≈ gmat) angs_ = convertangles.(pars.rotation, :GSLIB) @test all([x[1] for x in angs_] .≈ 1:10)