// (c) Tuomas J. Lukka #include namespace Vob { namespace Geom { /** Project a point through another point to a circle. * @param pt the point to project * @param p the projection center * @param ctr The center of the circle * @param rad The radius of the circle * @return The projected point, or if not possible, the original pt. */ Vec project2circle(Vec pt, Vec p, Vec ctr, float rad) { Vec ao = pt - ctr; Vec ap = pt - p; // Coefficients of the 2nd degree equation float a = dot(ap, ap); float b = 2*dot(ap, ao); float c = dot(ao, ao) - rad * rad; // determinant of the equation float det = b*b - 4*a*c; float ans = (det > 0 ? (-b + sqrt(det)) / (2*a) : 0); return pt + ans * ap; } } }