diff --git a/experimental/Intersection/NearestPoint.cpp b/experimental/Intersection/NearestPoint.cpp
new file mode 100644
index 0000000..f5b2828
--- /dev/null
+++ b/experimental/Intersection/NearestPoint.cpp
@@ -0,0 +1,473 @@
+/*
+Solving the Nearest Point-on-Curve Problem 
+and
+A Bezier Curve-Based Root-Finder
+by Philip J. Schneider
+from "Graphics Gems", Academic Press, 1990
+*/
+
+ /*	point_on_curve.c	*/		
+									
+#include <stdio.h>
+#include <malloc.h>
+#include <math.h>
+#include "GraphicsGems.h"
+
+#define TESTMODE
+
+/*
+ *  Forward declarations
+ */
+Point2  NearestPointOnCurve();
+static	int	FindRoots();
+static	Point2	*ConvertToBezierForm();
+static	double	ComputeXIntercept();
+static	int	ControlPolygonFlatEnough();
+static	int	CrossingCount();
+static	Point2	Bezier();
+static	Vector2	V2ScaleII();
+
+int		MAXDEPTH = 64;	/*  Maximum depth for recursion */
+
+#define	EPSILON	(ldexp(1.0,-MAXDEPTH-1)) /*Flatness control value */
+#define	DEGREE	3			/*  Cubic Bezier curve		*/
+#define	W_DEGREE 5			/*  Degree of eqn to find roots of */
+
+#ifdef TESTMODE
+/*
+ *  main :
+ *	Given a cubic Bezier curve (i.e., its control points), and some
+ *	arbitrary point in the plane, find the point on the curve
+ *	closest to that arbitrary point.
+ */
+main()
+{
+   
+ static Point2 bezCurve[4] = {	/*  A cubic Bezier curve	*/
+	{ 0.0, 0.0 },
+	{ 1.0, 2.0 },
+	{ 3.0, 3.0 },
+	{ 4.0, 2.0 },
+    };
+    static Point2 arbPoint = { 3.5, 2.0 }; /*Some arbitrary point*/
+    Point2	pointOnCurve;		 /*  Nearest point on the curve */
+
+    /*  Find the closest point */
+    pointOnCurve = NearestPointOnCurve(arbPoint, bezCurve);
+    printf("pointOnCurve : (%4.4f, %4.4f)\n", pointOnCurve.x,
+		pointOnCurve.y);
+}
+#endif /* TESTMODE */
+
+
+/*
+ *  NearestPointOnCurve :
+ *  	Compute the parameter value of the point on a Bezier
+ *		curve segment closest to some arbtitrary, user-input point.
+ *		Return the point on the curve at that parameter value.
+ *
+ */
+Point2 NearestPointOnCurve(P, V)
+    Point2 	P;			/* The user-supplied point	  */
+    Point2 	*V;			/* Control points of cubic Bezier */
+{
+    Point2	*w;			/* Ctl pts for 5th-degree eqn	*/
+    double 	t_candidate[W_DEGREE];	/* Possible roots		*/     
+    int 	n_solutions;		/* Number of roots found	*/
+    double	t;			/* Parameter value of closest pt*/
+
+    /*  Convert problem to 5th-degree Bezier form	*/
+    w = ConvertToBezierForm(P, V);
+
+    /* Find all possible roots of 5th-degree equation */
+    n_solutions = FindRoots(w, W_DEGREE, t_candidate, 0);
+    free((char *)w);
+
+    /* Compare distances of P to all candidates, and to t=0, and t=1 */
+    {
+		double 	dist, new_dist;
+		Point2 	p;
+		Vector2  v;
+		int		i;
+
+	
+	/* Check distance to beginning of curve, where t = 0	*/
+		dist = V2SquaredLength(V2Sub(&P, &V[0], &v));
+        	t = 0.0;
+
+	/* Find distances for candidate points	*/
+        for (i = 0; i < n_solutions; i++) {
+	    	p = Bezier(V, DEGREE, t_candidate[i],
+			(Point2 *)NULL, (Point2 *)NULL);
+	    	new_dist = V2SquaredLength(V2Sub(&P, &p, &v));
+	    	if (new_dist < dist) {
+                	dist = new_dist;
+	        		t = t_candidate[i];
+    	    }
+        }
+
+	/* Finally, look at distance to end point, where t = 1.0 */
+		new_dist = V2SquaredLength(V2Sub(&P, &V[DEGREE], &v));
+        	if (new_dist < dist) {
+            	dist = new_dist;
+	    	t = 1.0;
+        }
+    }
+
+    /*  Return the point on the curve at parameter value t */
+    printf("t : %4.12f\n", t);
+    return (Bezier(V, DEGREE, t, (Point2 *)NULL, (Point2 *)NULL));
+}
+
+
+/*
+ *  ConvertToBezierForm :
+ *		Given a point and a Bezier curve, generate a 5th-degree
+ *		Bezier-format equation whose solution finds the point on the
+ *      curve nearest the user-defined point.
+ */
+static Point2 *ConvertToBezierForm(P, V)
+    Point2 	P;			/* The point to find t for	*/
+    Point2 	*V;			/* The control points		*/
+{
+    int 	i, j, k, m, n, ub, lb;	
+    int 	row, column;		/* Table indices		*/
+    Vector2 	c[DEGREE+1];		/* V(i)'s - P			*/
+    Vector2 	d[DEGREE];		/* V(i+1) - V(i)		*/
+    Point2 	*w;			/* Ctl pts of 5th-degree curve  */
+    double 	cdTable[3][4];		/* Dot product of c, d		*/
+    static double z[3][4] = {	/* Precomputed "z" for cubics	*/
+	{1.0, 0.6, 0.3, 0.1},
+	{0.4, 0.6, 0.6, 0.4},
+	{0.1, 0.3, 0.6, 1.0},
+    };
+
+
+    /*Determine the c's -- these are vectors created by subtracting*/
+    /* point P from each of the control points				*/
+    for (i = 0; i <= DEGREE; i++) {
+		V2Sub(&V[i], &P, &c[i]);
+    }
+    /* Determine the d's -- these are vectors created by subtracting*/
+    /* each control point from the next					*/
+    for (i = 0; i <= DEGREE - 1; i++) { 
+		d[i] = V2ScaleII(V2Sub(&V[i+1], &V[i], &d[i]), 3.0);
+    }
+
+    /* Create the c,d table -- this is a table of dot products of the */
+    /* c's and d's							*/
+    for (row = 0; row <= DEGREE - 1; row++) {
+		for (column = 0; column <= DEGREE; column++) {
+	    	cdTable[row][column] = V2Dot(&d[row], &c[column]);
+		}
+    }
+
+    /* Now, apply the z's to the dot products, on the skew diagonal*/
+    /* Also, set up the x-values, making these "points"		*/
+    w = (Point2 *)malloc((unsigned)(W_DEGREE+1) * sizeof(Point2));
+    for (i = 0; i <= W_DEGREE; i++) {
+		w[i].y = 0.0;
+		w[i].x = (double)(i) / W_DEGREE;
+    }
+
+    n = DEGREE;
+    m = DEGREE-1;
+    for (k = 0; k <= n + m; k++) {
+		lb = MAX(0, k - m);
+		ub = MIN(k, n);
+		for (i = lb; i <= ub; i++) {
+	    	j = k - i;
+	    	w[i+j].y += cdTable[j][i] * z[j][i];
+		}
+    }
+
+    return (w);
+}
+
+
+/*
+ *  FindRoots :
+ *	Given a 5th-degree equation in Bernstein-Bezier form, find
+ *	all of the roots in the interval [0, 1].  Return the number
+ *	of roots found.
+ */
+static int FindRoots(w, degree, t, depth)
+    Point2 	*w;			/* The control points		*/
+    int 	degree;		/* The degree of the polynomial	*/
+    double 	*t;			/* RETURN candidate t-values	*/
+    int 	depth;		/* The depth of the recursion	*/
+{  
+    int 	i;
+    Point2 	Left[W_DEGREE+1],	/* New left and right 		*/
+    	  	Right[W_DEGREE+1];	/* control polygons		*/
+    int 	left_count,		/* Solution count from		*/
+		right_count;		/* children			*/
+    double 	left_t[W_DEGREE+1],	/* Solutions from kids		*/
+	   		right_t[W_DEGREE+1];
+
+    switch (CrossingCount(w, degree)) {
+       	case 0 : {	/* No solutions here	*/
+	     return 0;	
+	}
+	case 1 : {	/* Unique solution	*/
+	    /* Stop recursion when the tree is deep enough	*/
+	    /* if deep enough, return 1 solution at midpoint 	*/
+	    if (depth >= MAXDEPTH) {
+			t[0] = (w[0].x + w[W_DEGREE].x) / 2.0;
+			return 1;
+	    }
+	    if (ControlPolygonFlatEnough(w, degree)) {
+			t[0] = ComputeXIntercept(w, degree);
+			return 1;
+	    }
+	    break;
+	}
+}
+
+    /* Otherwise, solve recursively after	*/
+    /* subdividing control polygon		*/
+    Bezier(w, degree, 0.5, Left, Right);
+    left_count  = FindRoots(Left,  degree, left_t, depth+1);
+    right_count = FindRoots(Right, degree, right_t, depth+1);
+
+
+    /* Gather solutions together	*/
+    for (i = 0; i < left_count; i++) {
+        t[i] = left_t[i];
+    }
+    for (i = 0; i < right_count; i++) {
+ 		t[i+left_count] = right_t[i];
+    }
+
+    /* Send back total number of solutions	*/
+    return (left_count+right_count);
+}
+
+
+/*
+ * CrossingCount :
+ *	Count the number of times a Bezier control polygon 
+ *	crosses the 0-axis. This number is >= the number of roots.
+ *
+ */
+static int CrossingCount(V, degree)
+    Point2	*V;			/*  Control pts of Bezier curve	*/
+    int		degree;			/*  Degreee of Bezier curve 	*/
+{
+    int 	i;	
+    int 	n_crossings = 0;	/*  Number of zero-crossings	*/
+    int		sign, old_sign;		/*  Sign of coefficients	*/
+
+    sign = old_sign = SGN(V[0].y);
+    for (i = 1; i <= degree; i++) {
+		sign = SGN(V[i].y);
+		if (sign != old_sign) n_crossings++;
+		old_sign = sign;
+    }
+    return n_crossings;
+}
+
+
+
+/*
+ *  ControlPolygonFlatEnough :
+ *	Check if the control polygon of a Bezier curve is flat enough
+ *	for recursive subdivision to bottom out.
+ *
+ */
+static int ControlPolygonFlatEnough(V, degree)
+    Point2	*V;		/* Control points	*/
+    int 	degree;		/* Degree of polynomial	*/
+{
+    int 	i;			/* Index variable		*/
+    double 	*distance;		/* Distances from pts to line	*/
+    double 	max_distance_above;	/* maximum of these		*/
+    double 	max_distance_below;
+    double 	error;			/* Precision of root		*/
+    double 	intercept_1,
+    	   	intercept_2,
+	   		left_intercept,
+		   	right_intercept;
+    double 	a, b, c;		/* Coefficients of implicit	*/
+    					/* eqn for line from V[0]-V[deg]*/
+
+    /* Find the  perpendicular distance		*/
+    /* from each interior control point to 	*/
+    /* line connecting V[0] and V[degree]	*/
+    distance = (double *)malloc((unsigned)(degree + 1) * 					sizeof(double));
+    {
+	double	abSquared;
+
+	/* Derive the implicit equation for line connecting first *'
+    /*  and last control points */
+	a = V[0].y - V[degree].y;
+	b = V[degree].x - V[0].x;
+	c = V[0].x * V[degree].y - V[degree].x * V[0].y;
+
+	abSquared = (a * a) + (b * b);
+
+        for (i = 1; i < degree; i++) {
+	    /* Compute distance from each of the points to that line	*/
+	    	distance[i] = a * V[i].x + b * V[i].y + c;
+	    	if (distance[i] > 0.0) {
+				distance[i] = (distance[i] * distance[i]) / abSquared;
+	    	}
+	    	if (distance[i] < 0.0) {
+				distance[i] = -((distance[i] * distance[i]) / 						abSquared);
+	    	}
+		}
+    }
+
+
+    /* Find the largest distance	*/
+    max_distance_above = 0.0;
+    max_distance_below = 0.0;
+    for (i = 1; i < degree; i++) {
+		if (distance[i] < 0.0) {
+	    	max_distance_below = MIN(max_distance_below, distance[i]);
+		};
+		if (distance[i] > 0.0) {
+	    	max_distance_above = MAX(max_distance_above, distance[i]);
+		}
+    }
+    free((char *)distance);
+
+    {
+	double	det, dInv;
+	double	a1, b1, c1, a2, b2, c2;
+
+	/*  Implicit equation for zero line */
+	a1 = 0.0;
+	b1 = 1.0;
+	c1 = 0.0;
+
+	/*  Implicit equation for "above" line */
+	a2 = a;
+	b2 = b;
+	c2 = c + max_distance_above;
+
+	det = a1 * b2 - a2 * b1;
+	dInv = 1.0/det;
+	
+	intercept_1 = (b1 * c2 - b2 * c1) * dInv;
+
+	/*  Implicit equation for "below" line */
+	a2 = a;
+	b2 = b;
+	c2 = c + max_distance_below;
+	
+	det = a1 * b2 - a2 * b1;
+	dInv = 1.0/det;
+	
+	intercept_2 = (b1 * c2 - b2 * c1) * dInv;
+    }
+
+    /* Compute intercepts of bounding box	*/
+    left_intercept = MIN(intercept_1, intercept_2);
+    right_intercept = MAX(intercept_1, intercept_2);
+
+    error = 0.5 * (right_intercept-left_intercept);    
+    if (error < EPSILON) {
+		return 1;
+    }
+    else {
+		return 0;
+    }
+}
+
+
+
+/*
+ *  ComputeXIntercept :
+ *	Compute intersection of chord from first control point to last
+ *  	with 0-axis.
+ * 
+ */
+/* NOTE: "T" and "Y" do not have to be computed, and there are many useless
+ * operations in the following (e.g. "0.0 - 0.0").
+ */
+static double ComputeXIntercept(V, degree)
+    Point2 	*V;			/*  Control points	*/
+    int		degree; 		/*  Degree of curve	*/
+{
+    double	XLK, YLK, XNM, YNM, XMK, YMK;
+    double	det, detInv;
+    double	S, T;
+    double	X, Y;
+
+    XLK = 1.0 - 0.0;
+    YLK = 0.0 - 0.0;
+    XNM = V[degree].x - V[0].x;
+    YNM = V[degree].y - V[0].y;
+    XMK = V[0].x - 0.0;
+    YMK = V[0].y - 0.0;
+
+    det = XNM*YLK - YNM*XLK;
+    detInv = 1.0/det;
+
+    S = (XNM*YMK - YNM*XMK) * detInv;
+/*  T = (XLK*YMK - YLK*XMK) * detInv; */
+
+    X = 0.0 + XLK * S;
+/*  Y = 0.0 + YLK * S; */
+
+    return X;
+}
+
+
+/*
+ *  Bezier : 
+ *	Evaluate a Bezier curve at a particular parameter value
+ *      Fill in control points for resulting sub-curves if "Left" and
+ *	"Right" are non-null.
+ * 
+ */
+static Point2 Bezier(V, degree, t, Left, Right)
+    int 	degree;		/* Degree of bezier curve	*/
+    Point2 	*V;			/* Control pts			*/
+    double 	t;			/* Parameter value		*/
+    Point2 	*Left;		/* RETURN left half ctl pts	*/
+    Point2 	*Right;		/* RETURN right half ctl pts	*/
+{
+    int 	i, j;		/* Index variables	*/
+    Point2 	Vtemp[W_DEGREE+1][W_DEGREE+1];
+
+
+    /* Copy control points	*/
+    for (j =0; j <= degree; j++) {
+		Vtemp[0][j] = V[j];
+    }
+
+    /* Triangle computation	*/
+    for (i = 1; i <= degree; i++) {	
+		for (j =0 ; j <= degree - i; j++) {
+	    	Vtemp[i][j].x =
+	      		(1.0 - t) * Vtemp[i-1][j].x + t * Vtemp[i-1][j+1].x;
+	    	Vtemp[i][j].y =
+	      		(1.0 - t) * Vtemp[i-1][j].y + t * Vtemp[i-1][j+1].y;
+		}
+    }
+    
+    if (Left != NULL) {
+		for (j = 0; j <= degree; j++) {
+	    	Left[j]  = Vtemp[j][0];
+		}
+    }
+    if (Right != NULL) {
+		for (j = 0; j <= degree; j++) {
+	    	Right[j] = Vtemp[degree-j][j];
+		}
+    }
+
+    return (Vtemp[degree][0]);
+}
+
+static Vector2 V2ScaleII(v, s)
+    Vector2	*v;
+    double	s;
+{
+    Vector2 result;
+
+    result.x = v->x * s; result.y = v->y * s;
+    return (result);
+}