Closest point on circle in 3d. What's missing?

Question

I hope I manage to explain this clearly. I am trying to calculate the closest point on a circle in 3D. I found the following solution: http://www.geometrictools.com/Documentation/DistancePoint3Circle3.pdf

My code is below (written in Lua). The main problem is that the projection Q seems to be incorrect, or I don't understand how to calculate it correctly. As you can read in the paper Q should be the projection of the point on the plane of the circle.

For example the normal of the circle is {0,1,0} and its center is located at {3, 3, 3}. My point (p) for which I'm trying to calculate the closest distance to the circle is located at {6, 3, 2}. Then, in my calculation, the projection Q onto the plane of the circle is {6, 0, 2}.

In order to make the algorithm work I seem to have to offset Q with the position of the plane eg the circle center component in the direction of its normal. In this case the y direction so value 3.

I can hack this for normal {0,1,0} because its easy, but once the circle will face any arbitrary position I don't know how to calculate this.

What am I missing and where am I going wrong?

function calculatePointCircleDistance(p, circleCenter, circleNormal, circleRadius)
local C = circleCenter
local R = circleRadius
local Q = projectVectorOntoPlane(p, circleNormal)

-- I need to do a fix like this in order to get the calculations right
-- This for example only works with circleNormal {0,1,0}
-- Adding the y component of the circle position to the projection Q
Q[2] = C[2]

if vec3.equal(Q, C) == 1 then
    print("point exacly aligned with center circle")
    return vec3.mag(vec3.sub(C, p)), C
end

-- the following is calculating X=C+R (Q−C / |Q−C|)
local QminC = vec3.sub(Q, C)
local tmp = vec3.scale(vec3.div(QminC, vec3.mag(QminC)), R)
local X = vec3.add(C, tmp)

-- return distance as |X-p| as well as point X
return vec3.mag(vec3.sub(X, p)), X
end



function projectVectorOntoPlane(v, normal)
-- U = V - (V dot N)N
local vProjected = vec3.sub(v, vec3.scale(normal, vec3.dot(v, normal)))
return vProjected
end

Answer 1

That paper you linked to makes a bit of a meal of this operation, I think.

Your problem is that projectVectorOntoPlane does not actually project the vector onto the plane you want. It projects the vector onto another plane that's parallel to the plane you want, but which passes through the origin. (You're then trying to fix this problem with your Q[2] = C[2], but this just makes things worse.)

A plane can be defined by a normal vector together with some point on the plane, so you could write the projectVectorOntoPlane function like this:

-- Project P onto the plane with normal n containing the point O.
function projectVectorOntoPlane(P, n, O)
    return vec3.sub(P, vec3.scale(n, vec3.dot(vec3.sub(P, O), n)))
end

However, for this problem it's simplest to work all the way through in a coordinate system based on the centre of the circle, so I suggest something like this:

-- Return a point on the circle with center C, unit normal n and radius r
-- that's closest to the point P. (If all points are closest, return any.)
function pointCircleClosest(P, C, n, r)
    -- Translate problem to C-centred coordinates.
    local P = vec3.sub(P, C)

    -- Project P onto the plane containing the circle.
    local Q = vec3.sub(P, vec3.scale(n, vec3.dot(n, P)))

    -- If Q is at the centre, all points on the circle are equally close.
    if vec3.equal(Q, {0,0,0}) then
        Q = perpendicular(n)
    end

    -- Now the nearest point lies on the line through the origin and Q.
    local R = vec3.sub(P, vec3.scale(Q, r / vec3.mag(Q)))

    -- Return to original coordinate system.
    return vec3.add(R, C)
end

-- Return an arbitrary vector that's perpendicular to n.
function perpendicular(n)
    if math.abs(n[1]) < math.abs(n[2]) then
        return vec3.cross(n, {1,0,0})
    else
        return vec3.cross(n, {0,1,0})
    end
end

Oh, and you might find it more convenient to use a nicer vec3 class, perhaps this one, so that you can write P - C instead of the fussy vec3.sub(P, C) and so on.