Simple rotation of a point about another

cormullion commented 3 years ago

I want to write a function that rotates a 3D point by three angles about another point:

function rotateby(pt::Point3D, about::Point3D, angleX, angleY, angleZ)
    ...
    return pt1
end

I think this should be possible (and even made easy) with Rotations.jl, but I could do with a helping hand to get started... (3Blue1Brown's videos about quaternions are great, but they're not in Julia... )

c42f commented 3 years ago

rotates a 3D point by three angles about another point

There's many ways of doing this, depending on the order of rotation axes. Let's assume you want RotXYZ(angleX, angleY, angleZ) for the rotation (reading from right to left — first rotate about Z, then Y, then X). You could equally well choose RotZYX or any of the various other orderings. With RotXYZ, we'd write your function as

using Rotations

function rotateby(point, about, angleX, angleY, angleZ)
    RotXYZ(angleX, angleY, angleZ) * (point - about) + about
end

to get

julia> rotateby(SA[2,0,0], SA[1,0,0], 0, -π/2, 0)
3-element SArray{Tuple{3},Float64,1,3} with indices SOneTo(3):
 1.0
 0.0
 1.0

Another way to do this is to use CoordinateTransformations.recenter(), which produces a transformation functor which can be applied to many points. With this you don't even really need the rotateby utility function:

julia> using CoordinateTransformations

julia> points = rand(SVector{3,Float64}, 4)
4-element Array{SArray{Tuple{3},Float64,1,3},1}:
 [0.9927615884534728, 0.09291848088449939, 0.39064074743025134]
 [0.697136339656516, 0.4664841478024728, 0.4145303821631099]
 [0.9655845551550752, 0.5179123320751253, 0.38367137002584073]
 [0.6205439290284316, 0.7179067253208282, 0.14154719893335987]

julia> xform = recenter(RotXYZ(0, -π/2, 0), SA[1,0,0])
AffineMap([6.123233995736766e-17 -0.0 -1.0; 0.0 1.0 -0.0; 1.0 0.0 6.123233995736766e-17], [0.9999999999999999, 0.0, -1.0])

julia> xform.(points)
4-element Array{SArray{Tuple{3},Float64,1,3},1}:
 [0.6093592525697487, 0.09291848088449939, -0.007238411546527157]
 [0.5854696178368901, 0.4664841478024728, -0.302863660343484]
 [0.6163286299741593, 0.5179123320751253, -0.03441544484492476]
 [0.85845280106664, 0.7179067253208282, -0.3794560709715684]

cormullion commented 3 years ago

Thanks for the help, Chris - much appreciated! Looks like I'd better switch to using StaticArrays for my Point3D stuff...

c42f commented 3 years ago

You're welcome! This should probably all work with any AbstractArray, including plain old Vector — if not, it's likely a bug. Though of course it'll be faster with a static array or other fixed length arrays of some type.

With SArray you get the SA syntax for constructing array literals, which might be nice if that's to your liking.

cormullion commented 3 years ago

I think it's starting to work:

Just have a show problem to solve:

julia> A
3-element Point3D with indices SOneTo(3):
Error showing value of type Point3D:
ERROR: MethodError: no method matching LinearIndices(::Int64)
Closest candidates are:
  LinearIndices(::Tuple{}) at indices.jl:451
  LinearIndices(::R) where {N, R<:Tuple{Vararg{AbstractUnitRange{Int64}, N}}} at indices.jl:448
  LinearIndices(::Tuple{Vararg{AbstractUnitRange{var"#s77"} where var"#s77"<:Integer, N}}) where N at indices.jl:452

Thanks again!

c42f commented 3 years ago

Nice! What is Point3D in this case?

cormullion commented 3 years ago

I'm using:

struct Point3D <: FieldVector{3, Float64}
    x::Float64
    y::Float64
    z::Float64
end

as a replacement for my previous:

struct Point3D
    x::Float64
    y::Float64
    z::Float64
end

and it's just the weird printing errors that I'm trying to isolate...

c42f commented 3 years ago

That is odd. Locally I get

julia> using StaticArrays
[ Info: Precompiling StaticArrays [90137ffa-7385-5640-81b9-e52037218182]

julia> struct Point3D <: FieldVector{3, Float64}
           x::Float64
           y::Float64
           z::Float64
       end

julia> Point3D(1,2,3)
3-element Point3D with indices SOneTo(3):
 1.0
 2.0
 3.0

c42f commented 3 years ago

This very type exists in the StaticArrays tests btw. Defined exactly the same way.

cormullion commented 3 years ago

I know, it's very peculiar - I've obviously done something somewhere that upsets that, because:

julia> A =  Point3D(0, 0, 0); # semicolon!

julia> A.x
0.0

julia> A.y
0.0

julia> A.z
0.0

julia> A
3-element Point3D with indices SOneTo(3):
Error showing value of type Point3D:
ERROR: MethodError: no method matching LinearIndices(::Int64)
Closest candidates are:
  LinearIndices(::Tuple{}) at indices.jl:451
  LinearIndices(::R) where {N, R<:Tuple{Vararg{AbstractUnitRange{Int64}, N}}} at indices.jl:448
  LinearIndices(::Tuple{Vararg{AbstractUnitRange{var"#s77"} where var"#s77"<:Integer, N}}) where N at indices.jl:452
  ...
Stacktrace:
  [1] isassigned
    @ ~/.julia/packages/StaticArrays/LJQEe/src/abstractarray.jl:29 [inlined]
  [2] alignment(io::IOContext{Base.TTY}, X::Point3D, rows::UnitRange{Int64}, cols::UnitRange{Int64}, cols_if_complete::Int64, cols_otherwise::Int64, sep::Int64)
    @ Base ./arrayshow.jl:67
  [3] _print_matrix(io::IOContext{Base.TTY}, X::AbstractVecOrMat{T} where T, pre::String, sep::String, post::String, hdots::String, vdots::String, ddots::String, hmod::Int64, vmod::Int64, rowsA::UnitRange{Int64}, colsA::UnitRange{Int64})
    @ Base ./arrayshow.jl:201
  [4] #invokelatest#2
    @ ./essentials.jl:707 [inlined]
  [5] invokelatest
    @ ./essentials.jl:706 [inlined]
  [6] print_matrix
    @ ./arrayshow.jl:170 [inlined]
  [7] print_array
    @ ./arrayshow.jl:327 [inlined]
  [8] show(io::IOContext{Base.TTY}, #unused#::MIME{Symbol("text/plain")}, X::Point3D)
    @ Base ./arrayshow.jl:368
  [9] (::REPL.var"#38#39"{REPL.REPLDisplay{LineEditREPL}, MIME{Symbol("text/plain")}, Base.RefValue{Any}})(io::Any)
    @ REPL /Users/julia/src/julia/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:220
 [10] with_repl_linfo(f::Any, repl::LineEditREPL)
    @ REPL /Users/julia/src/julia/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:462
 [11] display(d::REPL.REPLDisplay, mime::MIME{Symbol("text/plain")}, x::Any)
    @ REPL /Users/julia/src/julia/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:213
 [12] display(d::REPL.REPLDisplay, x::Any)
    @ REPL /Users/julia/src/julia/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:225
 [13] display(x::Any)
    @ Base.Multimedia ./multimedia.jl:328
 [14] #15
    @ ~/.julia/packages/Media/ItEPc/src/compat.jl:28 [inlined]
 [15] hookless(f::Media.var"#15#16"{Point3D})
    @ Media ~/.julia/packages/Media/ItEPc/src/compat.jl:14
 [16] render
    @ ~/.julia/packages/Media/ItEPc/src/compat.jl:27 [inlined]
 [17] render(x::Point3D)
    @ Media ~/.julia/packages/Media/ItEPc/src/system.jl:160
 [18] display(#unused#::Media.DisplayHook, x::Point3D)
    @ Media ~/.julia/packages/Media/ItEPc/src/compat.jl:9
 [19] display(x::Any)
    @ Base.Multimedia ./multimedia.jl:328
 [20] #invokelatest#2
    @ ./essentials.jl:707 [inlined]
 [21] invokelatest
    @ ./essentials.jl:706 [inlined]
 [22] print_response(errio::IO, response::Any, show_value::Bool, have_color::Bool, specialdisplay::Union{Nothing, AbstractDisplay})
    @ REPL /Users/julia/src/julia/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:247
 [23] (::REPL.var"#40#41"{LineEditREPL, Pair{Any, Bool}, Bool, Bool})(io::Any)
    @ REPL /Users/julia/src/julia/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:231
 [24] with_repl_linfo(f::Any, repl::LineEditREPL)
    @ REPL /Users/julia/src/julia/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:462
 [25] print_response(repl::AbstractREPL, response::Any, show_value::Bool, have_color::Bool)
    @ REPL /Users/julia/src/julia/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:229
 [26] (::REPL.var"#do_respond#61"{Bool, Bool, VSCodeServer.var"#40#41"{LineEditREPL, REPL.LineEdit.Prompt}, LineEditREPL, REPL.LineEdit.Prompt})(s::REPL.LineEdit.MIState, buf::Any, ok::Bool)
    @ REPL /Users/julia/src/julia/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:798
 [27] #invokelatest#2
    @ ./essentials.jl:707 [inlined]
 [28] invokelatest
    @ ./essentials.jl:706 [inlined]
 [29] run_interface(terminal::REPL.Terminals.TextTerminal, m::REPL.LineEdit.ModalInterface, s::REPL.LineEdit.MIState)
    @ REPL.LineEdit /Users/julia/src/julia/usr/share/julia/stdlib/v1.6/REPL/src/LineEdit.jl:2441
 [30] run_frontend(repl::LineEditREPL, backend::REPL.REPLBackendRef)
    @ REPL /Users/julia/src/julia/usr/share/julia/stdlib/v1.6/REPL/src/REPL.jl:1126
 [31] (::REPL.var"#44#49"{LineEditREPL, REPL.REPLBackendRef})()
    @ REPL ./task.jl:406

c42f commented 3 years ago

It looks like maybe you defined size(::Point3D) to be something other than what the AbstractArray interface allows?

c42f commented 3 years ago

You should have

julia> size(Point3D(1,2,3))
(3,)

cormullion commented 3 years ago

Ah, I have:

# for broadcasting
Base.size(::Point3D) = 3

I'll try yours.

c42f commented 3 years ago

Don't overload it — StaticArrays does that for you ;-)

cormullion commented 3 years ago

Haha, you're a genius! Thanks so much... 😅😅😅

c42f commented 3 years ago

Heh :sweat_smile:

You probably shouldn't need any custom methods at all for your Point3D struct. At least, consider deleting any which could overlap with any "array-like" stuff. Otherwise you may run into other problems like this

JuliaGeometry / Rotations.jl

Simple rotation of a point about another #138