Calculation of the nearest point on a bézier spline using a subdivision algorithm, http://twitter.com/stemichels/status/485743201668841472

gistfile1.m

void SMSplineSubdivisionIterate_Private(CGPoint p1, CGPoint p2, CGPoint p3, CGPoint p4, CGFloat u, CGFloat du, BOOL *stop, BOOL(^block)(CGPoint p1, CGPoint p2, CGPoint p3, CGPoint p4, CGFloat u, CGFloat du, BOOL *stop)) {
    // Call block for the given spline
    if (!block(p1, p2, p3, p4, u, du, stop)) {
        // don't continue with the subdivisions
        return;
    }
    if (*stop) {
        // stop the whole iteration
        return;
    }
    
    // Calculate all the mid-points of the line segments
    CGPoint p12   = SMLineGetMidPoint(p1, p2);
    CGPoint p23   = SMLineGetMidPoint(p2, p3);
    CGPoint p34   = SMLineGetMidPoint(p3, p4);
    CGPoint p123  = SMLineGetMidPoint(p12, p23);
    CGPoint p234  = SMLineGetMidPoint(p23, p34);
    CGPoint p1234 = SMLineGetMidPoint(p123, p234);
    
    // Continue with the subdivision
    du = du / 2.0;
    SMSplineSubdivisionIterate_Private(p1, p12, p123, p1234, u, du, stop, block);
    if (*stop) {
        return;
    }
    SMSplineSubdivisionIterate_Private(p1234, p234, p34, p4, u + du, du, stop, block);
}

void SMSplineSubdivisionIterate(CGPoint p1, CGPoint p2, CGPoint p3, CGPoint p4, BOOL(^block)(CGPoint p1, CGPoint p2, CGPoint p3, CGPoint p4, CGFloat u, CGFloat du, BOOL *stop)) {
    __block BOOL stop = NO;
    SMSplineSubdivisionIterate_Private(p1, p2, p3, p4, 0.0, 1.0, &stop, block);
}

double SMSplineParameterForShortestDistance(CGPoint p1, CGPoint p2, CGPoint p3, CGPoint p4, CGPoint point) {
    __block CGFloat shortestDistance = CGPointDistanceToPointSquared(p4, point);
    __block double nearstU = 1.0;
    SMSplineSubdivisionIterate(p1, p2, p3, p4, ^BOOL(CGPoint p1, CGPoint p2, CGPoint p3, CGPoint p4, CGFloat u, CGFloat du, BOOL *stop) {
        if (SMSplineIsLinear(p1, p2, p3, p4)) {
            CGFloat currentNearstU = SMLineParameterForShortestDistance(p1, p4, point);
            NSPoint currentNearstPoint = SMLineGetPointAtParameter(p1, p4, currentNearstU);
            CGFloat currentShortestDistance = CGPointDistanceToPoint(point, currentNearstPoint);
            if (currentShortestDistance < shortestDistance) {
                shortestDistance = currentShortestDistance;
                nearstU = currentNearstU * du + u;
            }
            return NO;
        }
        return YES;
    });
    return nearstU;
}

CGPoint SMSplineGetPointAtParameter(CGPoint p1, CGPoint p2, CGPoint p3, CGPoint p4, double u) {
    if (u <= 0.0) {
        return p1;
    }
    if (u >= 1.0) {
        return p4;
    }
    CGFloat b0 = SMSplineTerm1(u);
    CGFloat b1 = SMSplineTerm2(u);
    CGFloat b2 = SMSplineTerm3(u);
    CGFloat b3 = SMSplineTerm4(u);
    
    CGPoint point;
    point.x = p1.x * b0 + p2.x * b1 + p3.x * b2 + p4.x * b3;
    point.y = p1.y * b0 + p2.y * b1 + p3.y * b2 + p4.y * b3;
    return point;
}

BOOL SMSplineIsLinear(CGPoint p1, CGPoint p2, CGPoint p3, CGPoint p4) {
  // see http://www.antigrain.com/research/adaptive_bezier/index.html
    
	CGFloat dx = p4.x - p1.x;
	CGFloat dy = p4.y - p1.y;
	
	CGFloat d2 = fabs(((p2.x - p4.x) * dy - (p2.y - p4.y) * dx));
	CGFloat d3 = fabs(((p3.x - p4.x) * dy - (p3.y - p4.y) * dx));
	
	return (d2 + d3) * (d2 + d3) <= m_distance_tolerance * (dx * dx + dy * dy);
}

CGPoint SMLineGetMidPoint(CGPoint p1, CGPoint p2) {
	return CGPointMake((p1.x + p2.x) * 0.5,
                     (p1.y + p2.y) * 0.5);
}

CGPoint SMLineGetPointAtParameter(CGPoint p1, CGPoint p2, double u) {
    if (u <= 0.0) {
        return p1;
    }
    if (u >= 1.0) {
        return p2;
    }
    if (u == 0.5) {
        return CGPointMake((p1.x + p2.x) * 0.5,
                           (p1.y + p2.y) * 0.5);
    }
    return CGPointMake(p1.x + (p2.x - p1.x) * u,
                       p1.y + (p2.y - p1.y) * u);
}

double SMLineParameterForShortestDistance(CGPoint p1, CGPoint p2, CGPoint point) {
    CGPoint n = CGPointSub(p1, p2);
    return MAX(0.0, MIN(1.0, CGPointDot(n, CGPointSub(p1, point)) / CGPointLengthSquared(n)));
}