curve_pick Subroutine

    subroutine curve_pick(a, frame, cam, m, hits, n)
        class(CurveArtist), intent(in), target :: a !! the artist
        class(ReferenceFrame), intent(in), target :: frame  !! the frame drawing it (for the hit)
        type(PickCamera), intent(in) :: cam         !! the pick camera
        real(wp), intent(in) :: m(4,4)              !! the frame's transform to the root frame
        type(PickHit), allocatable, intent(inout) :: hits(:)  !! hits are added to hits(1:n)
        integer, intent(inout) :: n                 !! number of hits
        real(wp), dimension(:,:), pointer :: pos
        real(wp) :: extra, lo(3), hi(3), c(3), p0(3), p1(3), px, s, depth, best_px, best_depth, best_s, t0, t1, frac
        integer :: j, first, count, i1, i2, i, best_i, best_j, nd, d1, d2
        logical :: xyz
        type(PickHit) :: h
        if (.not. associated(a%traj) .or. .not. a%chain%is_built()) return
        nd = a%n_drawn
        if (nd < 1) return
        d1 = a%first_drawn_                                     ! (the points drawn: d1 .. d2)
        d2 = d1 + nd - 1
        extra = 0.5_wp * a%width_ * cam%pixel_ratio
        xyz = is_position_xyz(a%sources_) .and. .not. transformed(a)   ! (transformed: art_point)
        if (xyz) call a%traj%positions(pos)
        best_px = huge(1.0_wp)
        best_depth = huge(1.0_wp)
        best_i = 0
        best_j = 0                                              ! (the hit is between points best_i and best_j)
        best_s = 0.0_wp
        if (nd == 1) then                                       ! (a single point)
            if (point_pick(cam, root_point(d1), px, depth, extra)) then
                best_px = px
                best_depth = depth
                best_i = d1
            end if
        end if
        do j = 1, a%chain%num_chunks()                          ! the segments drawn, chunk by chunk
            call a%chain%chunk_bounds(j, first, count, lo, hi)
            i1 = max(first, d1)
            i2 = min(first + count - 1, d2)
            if (count < 2 .or. i2 <= i1) cycle
            c = matmul(m(1:3, 1:3), 0.5_wp * (lo + hi)) + m(1:3, 4)
            if (.not. sphere_near(cam, c, 0.5_wp * norm2(hi - lo), extra)) cycle       ! (a chunk far off: skipped)
            p0 = root_point(i1)
            do i = i1, i2 - 1
                p1 = root_point(i + 1)
                if (segment_pick(cam, p0, p1, px, s, depth, extra)) then
                    if (better(px, depth, best_px, best_depth)) then
                        best_px = px
                        best_depth = depth
                        best_i = i
                        best_j = i + 1
                        best_s = s
                    end if
                end if
                p0 = p1
            end do
        end do
        ! trace mode: the tail from the latest point passed to the point at the time drawn
        if (a%trace .and. a%tail_shown) then
            t0 = a%traj%time(a%tail_from)
            t1 = a%traj%time(a%tail_to)
            frac = 0.0_wp
            if (t1 /= t0) frac = max(0.0_wp, min(1.0_wp, (a%tail_time - t0) / (t1 - t0)))
            p0 = root_point(a%tail_from)
            p1 = p0 + frac * (root_point(a%tail_to) - p0)
            if (segment_pick(cam, p0, p1, px, s, depth, extra)) then
                if (better(px, depth, best_px, best_depth)) then
                    best_px = px
                    best_depth = depth
                    best_i = a%tail_from
                    best_j = a%tail_to
                    best_s = s * frac
                end if
            end if
        end if
        if (best_i == 0) return
        h%kind = PICK_CURVE
        h%frame => frame
        h%artist => a
        h%trajectory => a%traj
        h%index = best_i
        if (best_j > 0) h%index = min(best_i, best_j)           ! (the point the hit is after, in point order)
        h%local = local_point(best_i)
        h%time = a%traj%time(best_i)
        if (best_j > 0) then
            h%local = h%local + best_s * (local_point(best_j) - h%local)
            h%time = h%time + best_s * (a%traj%time(best_j) - h%time)
        end if
        h%position = matmul(m(1:3, 1:3), h%local) + m(1:3, 4)
        h%distance = best_depth
        h%pixels = best_px
        call add_pick_hit(hits, n, h)

    contains

        !> Point i in the frame's coordinates (straight from the positions if they are the sources).
        function local_point(i) result(v)
            integer, intent(in) :: i                        !! the trajectory point
            real(wp) :: v(3)
            if (xyz) then
                v = pos(:, i)
            else
                v = art_point(a, i, a%sources_)
            end if
        end function local_point

        !> Point i in the root frame.
        function root_point(i) result(v)
            integer, intent(in) :: i                        !! the trajectory point
            real(wp) :: v(3)
            v = matmul(m(1:3, 1:3), local_point(i)) + m(1:3, 4)
        end function root_point

    end subroutine curve_pick