Fix ViewState element type, :rotation interpolation, Pause action, and record frame count

.gitignore

@@ -5,3 +5,4 @@
 /docs/Manifest.toml
 /docs/build/
 fly_animation.mp4
+fly_animation_path.mp4

Project.toml

@@ -7,6 +7,7 @@ version = "0.1.3"
 projects = ["test", "docs"]
 
 [deps]
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [weakdeps]
@@ -16,6 +17,7 @@ Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
 FlyThroughPathsMakieExt = "Makie"
 
 [compat]
-Makie = "0.21, 0.22, 0.23, 0.24"
+LinearAlgebra = "1"
+Makie = "0.24"
 StaticArrays = "1"
 julia = "1.10"

docs/src/makie.md

@@ -23,9 +23,10 @@ First, we extract the initial view state from the axis `ax`.
 ```@example simple
 view0 = capture_view(ax)
 ```
-Note that this `ViewState` is a Float32 object, since that's the space
-Makie cameras work in.  If you want this to be Float64, you can 
-simply `convert(ViewState{Float64}, view0)`.
+Note that the element type of this `ViewState` follows the space the Makie
+camera works in, which is Float32 for older versions of Makie and Float64
+for newer ones.  If you want a different element type, you can simply
+`convert(ViewState{Float64}, view0)`.
 
 ### Creating a path
 

ext/FlyThroughPathsMakieExt.jl

@@ -25,27 +25,29 @@ function FlyThroughPaths.set_view!(scene::Scene, view::ViewState)
 end
 FlyThroughPaths.set_view!(axis::Makie.AbstractAxis, view::ViewState) = set_view!(axis.scene, view)
 
-function Makie.record(fig::Makie.FigureLike, file::String, path::Path; framerate = 24, kwargs...)
-    tend = FlyThroughPaths.duration(path)
-    trange = LinRange(0, tend, round(Int, tend / framerate))
-    iterator = path.(trange)
-    return Makie.record(fig, file, iterator; framerate, kwargs...)
-end
-
 # Define the recipe
 import FlyThroughPaths: plotcamerapath, plotcamerapath!
-@recipe(PlotCameraPath, path, time) do scene
-    Attributes(
-        colormap = Makie.inherit(scene, :colormap, :plasma),
-        color = Makie.inherit(scene, :color, :black),
-        linewidth = Makie.inherit(scene, :linewidth, 1.0),
-        linestyle = Makie.inherit(scene, :linestyle, :solid),
-        camera_marker = Makie.inherit(scene, :marker, :none),
-        camera_color = Makie.inherit(scene, :color, :black),
-        camera_markersize = Vec3f(2, 2, 3),
-        density = 30, # points per second
-        cycle = [:color,],
-    )
+"""
+    plotcamerapath(path::Path, [time])
+
+Plot the eye positions along `path` as a line coloured by time, with an arrow showing
+where the camera is looking at `time` (0 by default).
+"""
+@recipe PlotCameraPath (path, time) begin
+    "Colormap for the path, which is coloured by time."
+    colormap = @inherit colormap :plasma
+    color = @inherit color :black
+    linewidth = @inherit linewidth 1.0
+    linestyle = @inherit linestyle :solid
+    camera_color = @inherit color :black
+    """
+    Scales the arrow marking the camera. `automatic` sizes it from the bounding box of
+    the path, which is usually what you want, since a path can span any distance.
+    """
+    camera_markerscale = Makie.automatic
+    "Sampling rate of the path, in points per second of path time."
+    density = 30
+    cycle = [:color]
 end
 
 Makie.convert_arguments(::Type{<: PlotCameraPath}, path::Path, time::Number) = (path, Float64(time))
@@ -57,8 +59,8 @@ function Makie.plot!(plot::PlotCameraPath)
     trange_obs = Observable{LinRange{Float64}}()
     onany(plot, plot.path, plot.density; update = true) do path, density
         tend = FlyThroughPaths.duration(path)
-        trange_obs.val = LinRange(0.0, Float64(tend), round(Int, tend*density))
-        eyepositions_obs.val = Makie.Point3d.(getproperty.(path.(trange_obs.val), :eyeposition))
+        trange_obs.val = LinRange(0.0, Float64(tend), FlyThroughPaths.nframes(path, density))
+        eyepositions_obs.val = Makie.Point3d.(getproperty.(path(trange_obs.val), :eyeposition))
         notify(eyepositions_obs)
         notify(trange_obs)
     end
@@ -90,15 +92,25 @@ function Makie.plot!(plot::PlotCameraPath)
         linewidth = plot.linewidth, 
         linestyle = plot.linestyle,
     )
-    arrows!(
-        plot, 
-        @lift([$eyeposition_obs]), 
+    # The camera arrow has to be sized against the path, not against itself: `automatic`
+    # would scale it by its own bounding box, which is a single unit-length arrow.
+    arrowscale_obs = lift(plot, plot.camera_markerscale, eyepositions_obs) do scale, eyepositions
+        scale isa Makie.Automatic || return Float64(scale)
+        length(eyepositions) < 2 && return 1.0
+        return 0.15 * maximum(Makie.widths(Rect3d(eyepositions)))
+    end
+
+    # `align = :tail` puts the arrow's base at the eye, so it points where the camera looks
+    arrows3d!(
+        plot,
+        @lift([$eyeposition_obs]),
         @lift([$viewdir_obs]);
-        color = plot.camera_color, 
-        arrowsize = plot.camera_markersize, 
-        normalize = true, 
-        shading = Makie.MultiLightShading,
-        align = :headstart,
+        color = plot.camera_color,
+        lengthscale = arrowscale_obs,
+        markerscale = arrowscale_obs,
+        normalize = true,
+        shading = true,
+        align = :tail,
     )
 
 

src/FlyThroughPaths.jl

@@ -1,5 +1,6 @@
 module FlyThroughPaths
 
+using LinearAlgebra
 using StaticArrays
 
 export ViewState, Path, Pause, ConstrainedMove, BezierMove

src/path.jl

@@ -27,16 +27,76 @@ end
 
 duration(path::Path{T}) where T = sum(duration, path.changes; init = zero(T))
 
+"""
+    nframes(path, rate)
+
+Return the number of samples needed to traverse `path` at `rate` samples per second,
+e.g. the number of frames to render at a given framerate.
+
+At least two samples are returned, so that a path shorter than one sampling interval
+still yields a non-degenerate range.
+"""
+nframes(path::Path, rate) = max(2, round(Int, duration(path) * rate))
+
 function (path::Path{T})(t) where T
     view = path.initialview
     tend = zero(T)
     t < tend && return view
     for change in path.changes
         tnext = tend + duration(change)
         if t <= tnext
-            return change(view, t - tend)
+            # `tend` is accumulated separately from `t`, so `t - tend` can land a few ulps
+            # outside `[0, duration(change)]` even though `t` selected this change. Clamp
+            # rather than let `checkt` reject a time we just decided belongs here.
+            return change(view, clamp(t - tend, zero(T), duration(change)))
         end
         tend, view = tnext, filldefaults(target(view, change), view)
     end
     return view
 end
+
+"""
+    (path::Path)(ts::AbstractVector)
+
+Evaluate `path` at every time in `ts`, which must be sorted, and return the resulting
+`Vector{ViewState}`.  The result is identical to `path.(ts)`, elementwise.
+
+Prefer this to broadcasting when sampling a whole path, e.g. once per frame of a video.
+The scalar method walks the path's changes from the beginning on every call, both to find
+the change that owns `t` and to accumulate the `ViewState` that change starts from; this
+method does that walk once and then locates each time by `searchsortedfirst` over the
+segment end times.
+"""
+function (path::Path{T})(ts::AbstractVector) where T
+    issorted(ts) || throw(ArgumentError("`ts` must be sorted; broadcast `path.(ts)` instead"))
+    changes = path.changes
+    nchanges = length(changes)
+    # The start time of each change, the view it starts from, and its end time, all
+    # accumulated exactly as the scalar method accumulates them
+    tstarts = Vector{T}(undef, nchanges)
+    tstops = Vector{T}(undef, nchanges)
+    startviews = Vector{ViewState{T}}(undef, nchanges)
+    tend = zero(T)
+    endview = path.initialview
+    for (i, change) in enumerate(changes)
+        tstarts[i], startviews[i] = tend, endview
+        tstops[i] = tend = tend + duration(change)
+        endview = filldefaults(target(endview, change), endview)
+    end
+    result = Vector{ViewState{T}}(undef, length(ts))
+    i = 1   # the changes are visited in order, since `ts` is sorted
+    for (k, t) in enumerate(ts)
+        if t < zero(T)
+            result[k] = path.initialview
+            continue
+        end
+        i <= nchanges && (i += searchsortedfirst(@view(tstops[i:nchanges]), t) - 1)
+        result[k] = if i > nchanges
+            endview
+        else
+            change = changes[i]
+            change(startviews[i], clamp(t - tstarts[i], zero(T), duration(change)))
+        end
+    end
+    return result
+end

src/pathchange.jl

@@ -38,7 +38,7 @@ end
 
 Pause at the current position for `duration`.
 """
-Pause(duration::T) where T = Pause{T}(duration)
+Pause(duration::T, action=nothing) where T = Pause{T}(duration, action)
 
 Base.convert(::Type{Pause{T}}, p::Pause) where T = Pause{T}(p.duration, p.action)
 Base.convert(::Type{PathChange{T}}, p::Pause) where T = convert(Pause{T}, p)
@@ -145,13 +145,71 @@ Base.@nospecializeinfer function act(@nospecialize(action), t::Real)
     return nothing
 end
 
+"""
+    slerp(vold, vnew, f)
+
+Interpolate between the vectors `vold` and `vnew` at fraction `f`, rotating the direction
+along the great circle joining them while interpolating the length geometrically. The
+length therefore varies monotonically between `norm(vold)` and `norm(vnew)`, and is
+constant when those are equal.
+
+The naive blend `cospi(f/2) * vold + sinpi(f/2) * vnew` is not a rotation: writing
+`d = norm(vold) = norm(vnew)` and `θ` for the angle between the two vectors, its squared
+length is `d^2 * (1 + sinpi(f) * cos(θ))`, which is `d^2` only for `θ = 90°`. At `θ = 0` it
+swells to `d^2*2` halfway through, and at `θ = 180°` it passes through zero, i.e. through
+the point being looked at.
+
+The rotation uses the tangent-vector form `cos(r)*a + sin(r)*dir`, where `a` is the unit
+direction of `vold`, `r = f*θ`, and `dir = normalize((a × b) × a)` is the unit tangent at
+`a` pointing towards the unit direction `b` of `vnew`. This is the formulation used by
+[GeometryOps.jl's `UnitSpherical.slerp`](https://github.com/JuliaGeo/GeometryOps.jl/blob/main/src/utils/UnitSpherical/slerp.jl),
+which adapts it from Google's S2 geometry library. It is preferred to the textbook
+`(sin((1-f)θ)*a + sin(f*θ)*b) / sin(θ)` because that divisor collapses as `θ` approaches
+0 or `π`, whereas the tangent form only has to special-case the two configurations where
+the plane of rotation is genuinely undetermined.
+"""
+function slerp(vold::SVector{3,T}, vnew::SVector{3,T}, f) where T
+    f <= 0 && return vold
+    f >= 1 && return vnew
+    dold, dnew = norm(vold), norm(vnew)
+    # A vector of zero length has no direction to rotate, so interpolate linearly instead
+    (iszero(dold) || iszero(dnew)) && return (1 - f) * vold + f * vnew
+    a, b = vold / dold, vnew / dnew
+    # Geometric interpolation of the length: a constant relative rate of approach reads
+    # more evenly than a linear one when the camera dollies in or out.
+    d = dold * (dnew / dold)^f
+    # `n` is the normal of the plane of rotation and `norm(n)` is `sin(θ)`. `n` carries an
+    # absolute error of a few `eps(T)` however close `a` and `b` are to each other or to
+    # antipodal, so a norm at that level -- and only then -- means the inputs do not
+    # determine a plane at all.
+    n = cross(a, b)
+    s, c = norm(n), dot(a, b)
+    θ = atan(s, c)
+    if s <= 8 * eps(T)
+        # θ ≈ 0: `a` and `b` are the same direction, so only the length changes
+        c > 0 && return d * a
+        # θ ≈ 180°: every plane containing `a` also contains `b`, so the direction of
+        # travel really is arbitrary. Rotate in the plane spanned by `a` and the
+        # coordinate axis it is least aligned with: an arbitrary choice, but a
+        # deterministic and well-conditioned one. (S2, and GeometryOps after it, resolve
+        # this case with exact arithmetic and symbolic perturbation instead, because
+        # their predicates must agree with each other from one call to the next. A camera
+        # fly-through has no such requirement, so that machinery is not reproduced here.)
+        i = argmin(abs.(a))
+        n = cross(a, SVector(ntuple(j -> T(j == i), 3)))
+    end
+    dir = normalize(cross(n, a))
+    r = f * θ
+    return d * normalize(cos(r) * a + sin(r) * dir)
+end
+
 # Compute the view from a PathChange at (relative) time t
 
 function (pause::Pause{T})(view::ViewState{T}, t) where T
     checkt(t, pause)
     action = pause.action
     if action !== nothing
-        tf = t / duration(move)
+        tf = t / duration(pause)
         act(action, tf)
     end
     return view
@@ -173,7 +231,7 @@ function (move::ConstrainedMove{T})(view::ViewState{T}, t) where T
     elseif constraint === :rotation
         vold = eyeposition - lookat
         vnew = eyeposition_new - lookat_new
-        eyeposition = cospi(f/2) * vold + sinpi(f/2) * vnew + lookatf
+        eyeposition = slerp(vold, vnew, f) + lookatf
     end
     upvector = (1 - f) * upvector + f * upvector_new
     fov = (1 - f) * fov + f * fov_new

src/viewstate.jl

@@ -7,7 +7,28 @@ end
 
 ViewState{T}(; eyeposition=nothing, lookat=nothing, upvector=nothing, fov=nothing) where T =
     ViewState{T}(eyeposition, lookat, upvector, fov)
-ViewState(; kwargs...) = ViewState{Float32}(; kwargs...)
+
+_eltype(::Nothing) = Union{}     # contributes nothing to the promotion
+_eltype(x::Number) = typeof(x)
+_eltype(x) = eltype(x)
+
+# `Union{}` (nothing was supplied) is a subtype of everything, hence the first test
+_floattype(::Type{T}) where T = T !== Union{} && T <: AbstractFloat ? T : Float32
+
+"""
+    ViewState(; eyeposition, lookat, upvector, fov)
+
+Construct a `ViewState` whose element type is promoted from the supplied arguments;
+e.g. `Point3d` positions or a `Float64` `fov` give a `ViewState{Float64}`.
+
+Integers express no preference about precision, so they (like omitting an argument
+altogether) keep the `Float32` default, `Float32` being the space Makie cameras work in.
+Use `ViewState{T}(; ...)` to choose the element type explicitly.
+"""
+function ViewState(; eyeposition=nothing, lookat=nothing, upvector=nothing, fov=nothing)
+    T = promote_type(_eltype(eyeposition), _eltype(lookat), _eltype(upvector), _eltype(fov))
+    return ViewState{_floattype(T)}(eyeposition, lookat, upvector, fov)
+end
 
 Base.convert(::Type{ViewState{T}}, v::ViewState) where T =
     ViewState{T}(v.eyeposition, v.lookat, v.upvector, v.fov)

test/glmakie.jl

@@ -26,3 +26,8 @@ tlist = range(0, stop=15, length=31)
 record(fig, "fly_animation.mp4", tlist; framerate=round(Int, length(tlist)/last(tlist))) do t
     set_view!(ax, path(t))
 end
+
+# The path itself can be plotted, in the space it flies through
+f2, a2, p2 = surface(-8..8, -8..8, Makie.peaks())
+FlyThroughPaths.plotcamerapath!(a2, path, 7)
+display(f2)