diff --git a/demo/fitCurve.html b/demo/fitCurve.html
new file mode 100644
index 0000000..6c341ff
--- /dev/null
+++ b/demo/fitCurve.html
@@ -0,0 +1,54 @@
+
+
+
+
+ 三次bezier曲线拟合
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/demo/js/fitCurve.js b/demo/js/fitCurve.js
new file mode 100644
index 0000000..77c793a
--- /dev/null
+++ b/demo/js/fitCurve.js
@@ -0,0 +1,161 @@
+/**
+ * @file bezier曲线拟合
+ * @author mengke01(kekee000@gmail.com)
+ */
+
+define(
+ function(require) {
+
+ var fitCurve = require('graphics/image/fitCurve');
+
+
+
+ /**
+ * 计算bezier曲线B0参数
+ *
+ * @param {number} u t值
+ * @return {number}
+ */
+ function B0(u) {
+ var tmp = 1.0 - u;
+ return (tmp * tmp * tmp);
+ }
+
+ /**
+ * 计算bezier曲线B1参数
+ *
+ * @param {number} u t值
+ * @return {number}
+ */
+ function B1(u) {
+ var tmp = 1.0 - u;
+ return (3 * u * (tmp * tmp));
+ }
+
+ /**
+ * 计算bezier曲线B2参数
+ *
+ * @param {number} u t值
+ * @return {number}
+ */
+ function B2(u) {
+ var tmp = 1.0 - u;
+ return (3 * u * u * tmp);
+ }
+
+ /**
+ * 计算bezier曲线B3参数
+ *
+ * @param {number} u t值
+ * @return {number}
+ */
+ function B3(u) {
+ return (u * u * u);
+ }
+
+ var entry = {
+
+ /**
+ * 初始化
+ */
+ init: function () {
+
+ var canvas = $('#canvas').get(0);
+ var ctx = canvas.getContext('2d');
+ var width = canvas.offsetWidth;
+ var height = canvas.offsetHeight;
+
+ var points =[{"x":167,"y":449},{"x":233,"y":331},{"x":431,"y":332},{"x":481,"y":417}];
+
+ $(points).each(function(index, item) {
+ $('[data-index="'+index+'"]').css({
+ left: points[index].x,
+ top: points[index].y
+ })
+ });
+
+ var point;
+
+ $('.point').on('mousedown', function(e) {
+ point = $(this);
+ });
+
+ $(document.body).on('mousemove', onMove);
+ $(document.body).on('mouseup', function(e) {
+ onMove.call(this, e);
+ point = null;
+ });
+
+ function onMove(e) {
+ var px = e.pageX;
+ var py = e.pageY;
+ if(point) {
+ point.css({
+ left: px,
+ top: py
+ });
+ var p = points[+point.attr('data-index')];
+ p.x = px;
+ p.y = py;
+ draw();
+ }
+ }
+
+ function draw() {
+ ctx.clearRect(0, 0, width, height);
+
+ ctx.beginPath();
+ ctx.strokeStyle='gray';
+ ctx.moveTo(points[0].x, points[0].y);
+ ctx.bezierCurveTo(points[1].x,points[1].y,points[2].x,points[2].y,points[3].x,points[3].y);
+ ctx.stroke();
+
+ var curvePoints = [];
+ var count = 40;
+ var tSegment = 1 / count;
+ for (var i = 0; i < count; i++) {
+ var t = tSegment * i;
+ curvePoints.push({
+ x: points[0].x * B0(t) + points[1].x * B1(t) + points[2].x * B2(t) + points[3].x * B3(t),
+ y: points[0].y * B0(t) + points[1].y * B1(t) + points[2].y * B2(t) + points[3].y * B3(t)
+ });
+ }
+
+ ctx.fillStyle='green';
+ for (var i = 0; i < count; i++) {
+ ctx.fillRect(curvePoints[i].x, curvePoints[i].y, 2, 2);
+ }
+
+ var result = fitCurve(curvePoints, 100);
+ result.unshift(points[0]);
+ console.log(result);
+
+ ctx.strokeStyle='blue';
+ for (var i = 3; i < result.length - 1; i+=3) {
+ var p = {
+ x: (result[i].x + result[i + 1].x) / 2,
+ y: (result[i].y + result[i + 1].y) / 2
+ };
+ result.push(p);
+ i++;
+ }
+
+ ctx.moveTo(result[0].x, result[0].y);
+ for (var i = 0; i < result.length - 1; i+=4) {
+ ctx.bezierCurveTo(result[i + 1].x, result[i + 1].y, result[i + 2].x, result[i + 2].y, result[i + 3].x, result[i + 3].y);
+ ctx.moveTo(result[i + 3].x, result[i + 3].y);
+ }
+
+ ctx.stroke();
+
+ }
+
+ draw();
+ }
+ };
+
+ entry.init();
+
+ return entry;
+ }
+);
diff --git a/src/graphics/image/fitCurve.js b/src/graphics/image/fitCurve.js
new file mode 100644
index 0000000..51729cb
--- /dev/null
+++ b/src/graphics/image/fitCurve.js
@@ -0,0 +1,534 @@
+/**
+ * @file 三次bezier曲线点拟合
+ *
+ * An Algorithm for Automatically Fitting Digitized Curves
+ * by Philip J. Schneider
+ * from "Graphics Gems", Academic Press, 1990
+ *
+ * @author mengke01(kekee000@gmail.com)
+ */
+
+
+define(
+ function (require) {
+
+ var MAX_ITERATIONS = 4; // 最大的迭代次数
+
+ /**
+ * 将点转换成余弦值
+ *
+ * @param {Object} p 点
+ * @return {Object} 余弦值
+ */
+ function normalize(p) {
+ var factor = 1 / Math.sqrt(p.x * p.x + p.y * p.y);
+ p.x = p.x * factor;
+ p.y = p.y * factor;
+ return p;
+ }
+
+ /**
+ * 求两点距离
+ *
+ * @param {Object} p0 p0
+ * @param {Object} p1 p1
+ * @return {number} 距离值
+ */
+ function dist(p0, p1) {
+ return Math.sqrt(Math.pow(p0.x - p1.x, 2) + Math.pow(p0.y - p1.y, 2));
+ }
+
+ /**
+ * 计算左侧正切点
+ *
+ * @param {Array} points 点集合
+ * @param {number} end 点索引
+ * @return {Object} 点
+ */
+ function computeLeftTangent(points, end) {
+ var x = points[end + 1].x - points[end].x;
+ var y = points[end + 1].y - points[end].y;
+ return normalize({
+ x: x,
+ y: y
+ });
+ }
+
+ /**
+ * 计算右侧正切点
+ *
+ * @param {Array} points 点集合
+ * @param {number} end 点索引
+ * @return {Object} 点
+ */
+ function computeRightTangent(points, end) {
+ var x = points[end - 1].x - points[end].x;
+ var y = points[end - 1].y - points[end].y;
+ return normalize({
+ x: x,
+ y: y
+ });
+ }
+
+ /**
+ * 计算中间正切点
+ *
+ * @param {Array} points 点集合
+ * @param {number} end 点索引
+ * @return {Object} 点
+ */
+ function computeCenterTangent(points, center) {
+ var x = (points[center - 1].x - points[center + 1].x) / 2;
+ var y = (points[center - 1].y - points[center + 1].y) / 2;
+ return normalize({
+ x: x,
+ y: y
+ });
+ }
+
+ /**
+ * 计算bezier曲线B0参数
+ *
+ * @param {number} u t值
+ * @return {number}
+ */
+ function B0(u) {
+ var tmp = 1.0 - u;
+ return (tmp * tmp * tmp);
+ }
+
+ /**
+ * 计算bezier曲线B1参数
+ *
+ * @param {number} u t值
+ * @return {number}
+ */
+ function B1(u) {
+ var tmp = 1.0 - u;
+ return (3 * u * (tmp * tmp));
+ }
+
+ /**
+ * 计算bezier曲线B2参数
+ *
+ * @param {number} u t值
+ * @return {number}
+ */
+ function B2(u) {
+ var tmp = 1.0 - u;
+ return (3 * u * u * tmp);
+ }
+
+ /**
+ * 计算bezier曲线B3参数
+ *
+ * @param {number} u t值
+ * @return {number}
+ */
+ function B3(u) {
+ return (u * u * u);
+ }
+
+ /**
+ * 计算向量内积
+ *
+ * @param {Object} a a点
+ * @param {Object} b b点
+ */
+ function v2Dot(a, b) {
+ return (a.x * b.x) + (a.y * b.y);
+ }
+
+ /*
+ * ChordLengthParameterize :
+ * Assign parameter values to digitized points
+ * using relative distances between points.
+ * 使用线段法计算各个点上的t值
+ *
+ * @param {Array} points 点集合
+ * @param {number} first 点索引
+ * @param {number} last 点索引
+ * @return {Array} t值集合
+ */
+ function chordLengthParameterize(points, first, last)
+ {
+ var i;
+ var u = [0];
+ for (i = first + 1; i <= last; i++) {
+ u[i - first] = u[i - first - 1] + dist(points[i - 1], points[i]);
+ }
+ var length = u[last - first];
+ for (i = first + 1; i <= last; i++) {
+ u[i - first] = u[i - first] / length;
+ }
+
+ return u;
+ }
+
+ /*
+ * Bezier :
+ * Evaluate a Bezier curve at a particular parameter value
+ *
+ */
+ /**
+ * 计算bezier点在特定位置起始点位置，使用连接点法，参考bezier曲线定义
+ *
+ * @param {number} degree 深度
+ * @param {Array} bezierPoints bezier点集合
+ * @param {number} t t值
+ */
+ function bezierII(degree, bezierPoints, t) {
+ var i;
+ var j;
+ var temp = []; // Local copy of control points
+
+ for (i = 0; i <= degree; i++)
+ {
+ temp[i] = {
+ x: bezierPoints[i].x,
+ y: bezierPoints[i].y
+ };
+ }
+
+ /* Triangle computation */
+ for (i = 1; i <= degree; i++)
+ {
+ for (j = 0; j <= degree - i; j++)
+ {
+ temp[j].x = (1.0 - t) * temp[j].x + t * temp[j + 1].x;
+ temp[j].y = (1.0 - t) * temp[j].y + t * temp[j + 1].y;
+ }
+ }
+
+ return temp[0];
+ }
+
+
+ /**
+ * 计算最大错误点的位置和错误值
+ *
+ * @param {Array} points 点集合
+ * @param {number} first 起始点
+ * @param {number} last 结束点
+ * @param {Array} bezCurve bezier点参数
+ * @param {Array} u t参数集合
+ *
+ * @return {Object} 最大错误点的位置和错误值
+ */
+ function computeMaxError(points, first, last, bezCurve, u, splitPoint) {
+ var maxDist = 0; // Maximum error
+
+ var splitPoint = Math.floor((last - first + 1) / 2);
+ for (var i = first + 1; i < last; i++) {
+ var p = bezierII(3, bezCurve, u[i - first]);
+ var x = p.x - points[i].x;
+ var y = p.y - points[i].y;
+ var dist = x * x + y * y;
+ if (dist >= maxDist) {
+ maxDist = dist;
+ splitPoint = i;
+ }
+ }
+
+ return {
+ dist: maxDist,
+ index: splitPoint
+ };
+ }
+
+ /**
+ * 牛顿迭代法求优化bezier优化点
+ *
+ * @param {Array} Q bezier参数点
+ * @param {Object} P 当前点
+ * @param {number} u 当前点t值
+ * @return {u} 优化后的t值
+ */
+ function newtonRaphsonRootFind(Q, P, u) {
+ var numerator;
+ var denominator;
+ var Q1 = []; // Q'
+ var Q2 = []; // Q''
+ var Q_u; // u evaluated at Q, Q', & Q''
+ var Q1_u;
+ var Q2_u;
+ var i;
+
+ Q_u = bezierII(3, Q, u);
+
+ /* Generate control vertices for Q' */
+ for (i = 0; i <= 2; i++) {
+ Q1[i] = {
+ x: (Q[i + 1].x - Q[i].x) * 3.0,
+ y: (Q[i + 1].y - Q[i].y) * 3.0
+ };
+ }
+
+ /* Generate control vertices for Q'' */
+ for (i = 0; i <= 1; i++) {
+ Q2[i] = {
+ x: (Q[i + 1].x - Q[i].x) * 2.0,
+ y: (Q[i + 1].y - Q[i].y) * 2.0
+ };
+ }
+
+ /* Compute Q'(u) and Q''(u) */
+ Q1_u = bezierII(2, Q1, u);
+ Q2_u = bezierII(1, Q2, u);
+
+ /* Compute f(u)/f'(u) */
+ numerator = (Q_u.x - P.x) * (Q1_u.x) + (Q_u.y - P.y) * (Q1_u.y);
+ denominator = (Q1_u.x) * (Q1_u.x)
+ + (Q1_u.y) * (Q1_u.y)
+ + (Q_u.x - P.x) * (Q2_u.x)
+ + (Q_u.y - P.y) * (Q2_u.y);
+
+ if (denominator === 0) {
+ return u;
+ }
+
+ /* u = u - f(u)/f'(u) */
+ return u - (numerator / denominator);
+ }
+
+ /**
+ * 对当前的bezier参数和t值集合进行优化获取新的t值集合
+ *
+ * @param {Array} points 点集合
+ * @param {number} first 起始索引
+ * @param {number} last 结束索引
+ * @param {Array} u t值集合
+ * @param {Array} bezCurve bezier曲线点集
+ * @return {Array} 优化后的t值集合
+ */
+ function reparameterize(points, first, last, u, bezCurve) {
+ var uPrime = [];
+ for (var i = first; i <= last; i++) {
+ uPrime[i - first] = newtonRaphsonRootFind(bezCurve, points[i], u[i - first]);
+ }
+ return uPrime;
+ }
+
+
+ /**
+ * 根据当前参数集合生成bezier曲线点
+ *
+ * @param {Array} points 点集合
+ * @param {number} first 起始点
+ * @param {number} last 结束点
+ * @param {Array} uPrime t参数集合
+ * @param {Object} tHat1 左侧切线点
+ * @param {Object} tHat2 右侧切线点
+ * @return {Array} bezier曲线点
+ */
+ function generateBezier(points, first, last, uPrime, tHat1, tHat2) {
+
+ var A = [];
+ var C = [[0, 0], [0, 0]]; // Matrix C
+ var X = [0, 0]; // Matrix X
+ var tmp; // Utility variable
+ var u;
+ var nPts = last - first + 1; // number of pts in sub-curve
+
+ /* Compute the A's */
+ for (i = 0; i < nPts; i++) {
+ var v1 = {
+ x: tHat1.x,
+ y: tHat1.y
+ };
+ var v2 = {
+ x: tHat2.x,
+ y: tHat2.y
+ };
+ u = uPrime[i];
+ v1.x *= B1(u);
+ v1.y *= B1(u);
+ v2.x *= B2(u);
+ v2.y *= B2(u);
+ A[i] = [v1, v2];
+ }
+
+ var firstPoint = points[first]; // 起始点
+ var lastPoint = points[last]; // 结束点
+ for (var i = 0; i < nPts; i++) {
+ C[0][0] += v2Dot(A[i][0], A[i][0]);
+ C[0][1] += v2Dot(A[i][0], A[i][1]);
+ C[1][0] = C[0][1];
+ C[1][1] += v2Dot(A[i][1], A[i][1]);
+
+ var pi = points[first + i]; // 当前点
+ u = uPrime[i];
+ tmp = {
+ x: pi.x - (firstPoint.x * B0(u) + firstPoint.x * B1(u) + lastPoint.x * B2(u) + lastPoint.x * B3(u)),
+ y: pi.y - (firstPoint.y * B0(u) + firstPoint.y * B1(u) + lastPoint.y * B2(u) + lastPoint.y * B3(u))
+ };
+
+ X[0] += v2Dot(A[i][0], tmp);
+ X[1] += v2Dot(A[i][1], tmp);
+ }
+
+ /* Compute the determinants of C and X */
+ var det_C0_C1 = C[0][0] * C[1][1] - C[1][0] * C[0][1];
+ var det_C0_X = C[0][0] * X[1] - C[1][0] * X[0];
+ var det_X_C1 = X[0] * C[1][1] - X[1] * C[0][1];
+
+ /* Finally, derive alpha values */
+ var alpha_l = (det_C0_C1 == 0) ? 0.0 : det_X_C1 / det_C0_C1;
+ var alpha_r = (det_C0_C1 == 0) ? 0.0 : det_C0_X / det_C0_C1;
+
+ /* If alpha negative, use the Wu/Barsky heuristic (see text) */
+ /* (if alpha is 0, you get coincident control points that lead to
+ * divide by zero in any subsequent NewtonRaphsonRootFind() call. */
+ var segLength = Math.sqrt(Math.pow(firstPoint.x - lastPoint.x, 2) + Math.pow(firstPoint.y - lastPoint.y, 2));
+ var epsilon = 1.0e-6 * segLength;
+ var bezCurve = [];
+ bezCurve[0] = firstPoint;
+ bezCurve[3] = lastPoint;
+
+ if (alpha_l < epsilon || alpha_r < epsilon) {
+ /* fall back on standard (probably inaccurate) formula, and subdivide further if needed. */
+ var dist = segLength / 3.0;
+ bezCurve[1] = {
+ x: tHat1.x * dist + firstPoint.x,
+ y: tHat1.y * dist + firstPoint.y
+ };
+ bezCurve[2] = {
+ x: tHat2.x * dist + lastPoint.x,
+ y: tHat2.y * dist + lastPoint.y
+ };
+
+ return bezCurve;
+ }
+
+ /* First and last control points of the Bezier curve are */
+ /* positioned exactly at the first and last data points */
+ /* Control points 1 and 2 are positioned an alpha distance out */
+ /* on the tangent vectors, left and right, respectively */
+
+ bezCurve[1] = {
+ x: tHat1.x * alpha_l + firstPoint.x,
+ y: tHat1.y * alpha_l + firstPoint.y
+ };
+
+ bezCurve[2] = {
+ x: tHat1.x * alpha_r + firstPoint.x,
+ y: tHat1.y * alpha_r + firstPoint.y
+ };
+
+ return bezCurve;
+ }
+
+
+
+ /**
+ * 三次bezier曲线拟合
+ *
+ * @param {Array} points 点集合
+ * @param {number} first 起始点
+ * @param {number} last 结束点
+ * @param {Object} tHat1 起始切线点点
+ * @param {Object} tHat2 右侧切线点
+ * @param {number} error 最大错误
+ * @param {Array} result 结果点集合
+ */
+ function fitCubic(points, first, last, tHat1, tHat2, error, result) {
+ var bezCurve = []; // Control points of fitted Bezier curve
+ var nPts; // Number of points in subset
+ var i;
+
+ nPts = last - first + 1;
+
+ var firstPoint = points[first]; // 起始点
+ var lastPoint = points[last]; // 结束点
+ var segLength = Math.sqrt(Math.pow(firstPoint.x - lastPoint.x, 2) + Math.pow(firstPoint.y - lastPoint.y, 2));
+
+ // Use heuristic if region only has two points in it
+ if (nPts == 2) {
+ var dist = segLength / 3.0;
+
+ bezCurve[0] = firstPoint;
+ bezCurve[3] = lastPoint;
+ bezCurve[1] = {
+ x: tHat1.x * dist + firstPoint.x,
+ y: tHat1.y * dist + firstPoint.y
+ };
+ bezCurve[2] = {
+ x: tHat2.x * dist + lastPoint.x,
+ y: tHat2.y * dist + lastPoint.y
+ };
+
+ result.push(bezCurve[1]);
+ result.push(bezCurve[2]);
+ result.push(bezCurve[3]);
+ return;
+ }
+
+ /* Parameterize points, and attempt to fit curve */
+ var uPrime; // Improved parameter values
+ var maxError; // Maximum fitting error
+ var iterationError = error * error; // Error below which you try iterating
+ var u = chordLengthParameterize(points, first, last); // Parameter values for point
+
+ bezCurve = generateBezier(points, first, last, u, tHat1, tHat2);
+ /* Find max deviation of points to fitted curve */
+ var maxError = computeMaxError(points, first, last, bezCurve, u);
+
+ if (maxError.dist < error) {
+ result.push(bezCurve[1]);
+ result.push(bezCurve[2]);
+ result.push(bezCurve[3]);
+ return;
+ }
+
+
+ /* If error not too large, try some reparameterization */
+ /* and iteration */
+ if (maxError.dist < iterationError) {
+ for (i = 0; i < MAX_ITERATIONS; i++) {
+ uPrime = reparameterize(points, first, last, u, bezCurve);
+ bezCurve = generateBezier(points, first, last, uPrime, tHat1, tHat2);
+ maxError = computeMaxError(points, first, last, bezCurve, uPrime);
+
+ if (maxError.dist < error) {
+ result.push(bezCurve[1]);
+ result.push(bezCurve[2]);
+ result.push(bezCurve[3]);
+ return;
+ }
+
+ u = uPrime;
+ }
+ }
+
+ /* Fitting failed -- split at max error point and fit recursively */
+ var splitPoint = maxError.index;
+ var tHatCenter = computeCenterTangent(points, splitPoint); // Unit tangent vector at splitPoint
+ fitCubic(points, first, splitPoint, tHat1, tHatCenter, error, result);
+ tHatCenter.x = -tHatCenter.x;
+ tHatCenter.y = -tHatCenter.y;
+ fitCubic(points, splitPoint, last, tHatCenter, tHat2, error, result);
+ }
+
+
+
+ /**
+ * 三次bezier曲线点拟合点集
+ *
+ * @param {Array} points 点集合
+ * @param {number} error 最大错误距离
+ *
+ * @return {Array} 结果点集
+ */
+ function fitCurve(points, error) {
+ var last = points.length - 1;
+ var tHat1 = computeLeftTangent(points, 0);
+ var tHat2 = computeRightTangent(points, last);
+ var result = [];
+ fitCubic(points, 0, last, tHat1, tHat2, error, result);
+ return result;
+ }
+
+
+ return fitCurve;
+ }
+);
diff --git a/test/fitcurve.cs b/test/fitcurve.cs
new file mode 100644
index 0000000..53715e3
--- /dev/null
+++ b/test/fitcurve.cs
@@ -0,0 +1,499 @@
+using System;
+using System.Collections.Generic;
+using System.Drawing;
+using System.Linq;
+using System.Text;
+using System.Threading.Tasks;
+
+namespace Bezier
+{
+
+
+ public class Point
+ {
+ public double X { get; set; }
+ public double Y { get; set; }
+
+ public double Length
+ {
+ get
+ {
+ return Math.Sqrt(X * X + Y * Y);
+ }
+ }
+ public Point(double x = 0, double y = 0)
+ {
+ X = x;
+ Y = y;
+ }
+
+ public Point(Point p)
+ {
+ X = p.X;
+ Y = p.Y;
+ }
+
+ public static Point operator *(Point left, double right)
+ {
+ return new Point(left.X * right, left.Y * right);
+ }
+
+ public static Point operator /(Point left, double right)
+ {
+ return new Point(left.X / right, left.Y / right);
+ }
+
+ public static Point operator +(Point left, Point right)
+ {
+ return new Point(left.X + right.X, left.Y + right.Y);
+ }
+
+ public static Point operator -(Point left, Point right)
+ {
+ return new Point(left.X - right.X, left.Y - right.Y);
+ }
+
+ internal void Negate()
+ {
+ X = -X;
+ Y = -Y;
+ }
+
+ internal void Normalize()
+ {
+ double factor = 1.0 / Math.Sqrt(LengthSquared);
+ X *= factor;
+ Y *= factor;
+ }
+
+ public double LengthSquared { get { return X * X + Y * Y; } }
+
+
+ public static Point Add(Point a, Point b)
+ {
+ return new Point(a.X + b.X, a.Y + b.Y);
+ }
+
+ public static Point Subtract(Point a, Point b)
+ {
+ return new Point(a.X - b.X, a.Y - b.Y);
+ }
+ }
+
+
+ /*
+ An Algorithm for Automatically Fitting Digitized Curves
+ by Philip J. Schneider
+ from "Graphics Gems", Academic Press, 1990
+ */
+ public static class FitCurves
+ {
+ /* Fit the Bezier curves */
+
+ private const int MAXPOINTS = 10000;
+ public static List FitCurve(Point[] d, double error)
+ {
+ Point tHat1, tHat2; /* Unit tangent vectors at endpoints */
+
+ tHat1 = ComputeLeftTangent(d, 0);
+ tHat2 = ComputeRightTangent(d, d.Length - 1);
+ List result = new List();
+ FitCubic(d, 0, d.Length - 1, tHat1, tHat2, error, result);
+ return result;
+ }
+
+
+ private static void FitCubic(Point[] d, int first, int last, Point tHat1, Point tHat2, double error, List result)
+ {
+ Point[] bezCurve; /*Control points of fitted Bezier curve*/
+ double[] u; /* Parameter values for point */
+ double[] uPrime; /* Improved parameter values */
+ double maxError; /* Maximum fitting error */
+ int splitPoint; /* Point to split point set at */
+ int nPts; /* Number of points in subset */
+ double iterationError; /*Error below which you try iterating */
+ int maxIterations = 4; /* Max times to try iterating */
+ Point tHatCenter; /* Unit tangent vector at splitPoint */
+ int i;
+
+ iterationError = error * error;
+ nPts = last - first + 1;
+
+ /* Use heuristic if region only has two points in it */
+ if (nPts == 2)
+ {
+ double dist = (d[first] - d[last]).Length / 3.0;
+
+ bezCurve = new Point[4];
+ bezCurve[0] = d[first];
+ bezCurve[3] = d[last];
+ bezCurve[1] = (tHat1 * dist) + bezCurve[0];
+ bezCurve[2] = (tHat2 * dist) + bezCurve[3];
+
+ result.Add(bezCurve[1]);
+ result.Add(bezCurve[2]);
+ result.Add(bezCurve[3]);
+ return;
+ }
+
+ /* Parameterize points, and attempt to fit curve */
+ u = ChordLengthParameterize(d, first, last);
+ bezCurve = GenerateBezier(d, first, last, u, tHat1, tHat2);
+
+ /* Find max deviation of points to fitted curve */
+ maxError = ComputeMaxError(d, first, last, bezCurve, u, out splitPoint);
+ if (maxError < error)
+ {
+ result.Add(bezCurve[1]);
+ result.Add(bezCurve[2]);
+ result.Add(bezCurve[3]);
+ return;
+ }
+
+
+ /* If error not too large, try some reparameterization */
+ /* and iteration */
+ if (maxError < iterationError)
+ {
+ for (i = 0; i < maxIterations; i++)
+ {
+ uPrime = Reparameterize(d, first, last, u, bezCurve);
+ bezCurve = GenerateBezier(d, first, last, uPrime, tHat1, tHat2);
+ maxError = ComputeMaxError(d, first, last,
+ bezCurve, uPrime, out splitPoint);
+ if (maxError < error)
+ {
+ result.Add(bezCurve[1]);
+ result.Add(bezCurve[2]);
+ result.Add(bezCurve[3]);
+ return;
+ }
+ u = uPrime;
+ }
+ }
+
+ /* Fitting failed -- split at max error point and fit recursively */
+ tHatCenter = ComputeCenterTangent(d, splitPoint);
+ FitCubic(d, first, splitPoint, tHat1, tHatCenter, error, result);
+ tHatCenter.Negate();
+ FitCubic(d, splitPoint, last, tHatCenter, tHat2, error, result);
+ }
+
+ static Point[] GenerateBezier(Point[] d, int first, int last, double[] uPrime, Point tHat1, Point tHat2)
+ {
+ int i;
+ Point[,] A = new Point[MAXPOINTS, 2];/* Precomputed rhs for eqn */
+
+ int nPts; /* Number of pts in sub-curve */
+ double[,] C = new double[2, 2]; /* Matrix C */
+ double[] X = new double[2]; /* Matrix X */
+ double det_C0_C1, /* Determinants of matrices */
+ det_C0_X,
+ det_X_C1;
+ double alpha_l, /* Alpha values, left and right */
+ alpha_r;
+ Point tmp; /* Utility variable */
+ Point[] bezCurve = new Point[4]; /* RETURN bezier curve ctl pts */
+ nPts = last - first + 1;
+
+ /* Compute the A's */
+ for (i = 0; i < nPts; i++)
+ {
+ Point v1, v2;
+ v1 = tHat1;
+ v2 = tHat2;
+ v1 *= B1(uPrime[i]);
+ v2 *= B2(uPrime[i]);
+ A[i, 0] = v1;
+ A[i, 1] = v2;
+ }
+
+ /* Create the C and X matrices */
+ C[0, 0] = 0.0;
+ C[0, 1] = 0.0;
+ C[1, 0] = 0.0;
+ C[1, 1] = 0.0;
+ X[0] = 0.0;
+ X[1] = 0.0;
+
+ for (i = 0; i < nPts; i++)
+ {
+ C[0, 0] += V2Dot(A[i, 0], A[i, 0]);
+ C[0, 1] += V2Dot(A[i, 0], A[i, 1]);
+ /* C[1][0] += V2Dot(&A[i][0], &A[i][9]);*/
+ C[1, 0] = C[0, 1];
+ C[1, 1] += V2Dot(A[i, 1], A[i, 1]);
+
+ tmp = ((Point)d[first + i] -
+ (
+ ((Point)d[first] * B0(uPrime[i])) +
+ (
+ ((Point)d[first] * B1(uPrime[i])) +
+ (
+ ((Point)d[last] * B2(uPrime[i])) +
+ ((Point)d[last] * B3(uPrime[i]))))));
+
+
+ X[0] += V2Dot(A[i, 0], tmp);
+ X[1] += V2Dot(A[i, 1], tmp);
+ }
+
+ /* Compute the determinants of C and X */
+ det_C0_C1 = C[0, 0] * C[1, 1] - C[1, 0] * C[0, 1];
+ det_C0_X = C[0, 0] * X[1] - C[1, 0] * X[0];
+ det_X_C1 = X[0] * C[1, 1] - X[1] * C[0, 1];
+
+ /* Finally, derive alpha values */
+ alpha_l = (det_C0_C1 == 0) ? 0.0 : det_X_C1 / det_C0_C1;
+ alpha_r = (det_C0_C1 == 0) ? 0.0 : det_C0_X / det_C0_C1;
+
+ /* If alpha negative, use the Wu/Barsky heuristic (see text) */
+ /* (if alpha is 0, you get coincident control points that lead to
+ * divide by zero in any subsequent NewtonRaphsonRootFind() call. */
+ double segLength = (d[first] - d[last]).Length;
+ double epsilon = 1.0e-6 * segLength;
+ if (alpha_l < epsilon || alpha_r < epsilon)
+ {
+ /* fall back on standard (probably inaccurate) formula, and subdivide further if needed. */
+ double dist = segLength / 3.0;
+ bezCurve[0] = d[first];
+ bezCurve[3] = d[last];
+ bezCurve[1] = (tHat1 * dist) + bezCurve[0];
+ bezCurve[2] = (tHat2 * dist) + bezCurve[3];
+ return (bezCurve);
+ }
+
+ /* First and last control points of the Bezier curve are */
+ /* positioned exactly at the first and last data points */
+ /* Control points 1 and 2 are positioned an alpha distance out */
+ /* on the tangent vectors, left and right, respectively */
+ bezCurve[0] = d[first];
+ bezCurve[3] = d[last];
+ bezCurve[1] = (tHat1 * alpha_l) + bezCurve[0];
+ bezCurve[2] = (tHat2 * alpha_r) + bezCurve[3];
+ return (bezCurve);
+ }
+
+ /*
+ * Reparameterize:
+ * Given set of points and their parameterization, try to find
+ * a better parameterization.
+ *
+ */
+ static double[] Reparameterize(Point[] d, int first, int last, double[] u, Point[] bezCurve)
+ {
+ int nPts = last - first + 1;
+ int i;
+ double[] uPrime = new double[nPts]; /* New parameter values */
+
+ for (i = first; i <= last; i++)
+ {
+ uPrime[i - first] = NewtonRaphsonRootFind(bezCurve, d[i], u[i - first]);
+ }
+ return uPrime;
+ }
+
+
+
+ /*
+ * NewtonRaphsonRootFind :
+ * Use Newton-Raphson iteration to find better root.
+ */
+ static double NewtonRaphsonRootFind(Point[] Q, Point P, double u)
+ {
+ double numerator, denominator;
+ Point[] Q1 = new Point[3], Q2 = new Point[2]; /* Q' and Q'' */
+ Point Q_u, Q1_u, Q2_u; /*u evaluated at Q, Q', & Q'' */
+ double uPrime; /* Improved u */
+ int i;
+
+ /* Compute Q(u) */
+ Q_u = BezierII(3, Q, u);
+
+ /* Generate control vertices for Q' */
+ for (i = 0; i <= 2; i++)
+ {
+ Q1[i] = new Point((Q[i + 1].X - Q[i].X) * 3.0,(Q[i + 1].Y - Q[i].Y) * 3.0 );
+ }
+
+ /* Generate control vertices for Q'' */
+ for (i = 0; i <= 1; i++)
+ {
+ Q2[i] = new Point((Q[i + 1].X - Q[i].X) * 2.0, (Q[i + 1].Y - Q[i].Y) * 2.0);
+ }
+
+ /* Compute Q'(u) and Q''(u) */
+ Q1_u = BezierII(2, Q1, u);
+ Q2_u = BezierII(1, Q2, u);
+
+ /* Compute f(u)/f'(u) */
+ numerator = (Q_u.X - P.X) * (Q1_u.X) + (Q_u.Y - P.Y) * (Q1_u.Y);
+ denominator = (Q1_u.X) * (Q1_u.X) + (Q1_u.Y) * (Q1_u.Y) +
+ (Q_u.X - P.X) * (Q2_u.X) + (Q_u.Y - P.Y) * (Q2_u.Y);
+ if (denominator == 0.0f) return u;
+
+ /* u = u - f(u)/f'(u) */
+ uPrime = u - (numerator / denominator);
+ return (uPrime);
+ }
+
+
+
+ /*
+ * Bezier :
+ * Evaluate a Bezier curve at a particular parameter value
+ *
+ */
+ static Point BezierII(int degree, Point[] V, double t)
+ {
+ int i, j;
+ Point Q; /* Point on curve at parameter t */
+ Point[] Vtemp; /* Local copy of control points */
+
+ /* Copy array */
+ Vtemp = new Point[degree + 1];
+ for (i = 0; i <= degree; i++)
+ {
+ Vtemp[i] = new Point(V[i]);
+ }
+
+ /* Triangle computation */
+ for (i = 1; i <= degree; i++)
+ {
+ for (j = 0; j <= degree - i; j++)
+ {
+ Vtemp[j].X = (1.0 - t) * Vtemp[j].X + t * Vtemp[j + 1].X;
+ Vtemp[j].Y = (1.0 - t) * Vtemp[j].Y + t * Vtemp[j + 1].Y;
+ }
+ }
+
+ Q = Vtemp[0];
+ return Q;
+ }
+
+
+ /*
+ * B0, B1, B2, B3 :
+ * Bezier multipliers
+ */
+ static double B0(double u)
+ {
+ double tmp = 1.0 - u;
+ return (tmp * tmp * tmp);
+ }
+
+
+ static double B1(double u)
+ {
+ double tmp = 1.0 - u;
+ return (3 * u * (tmp * tmp));
+ }
+
+ static double B2(double u)
+ {
+ double tmp = 1.0 - u;
+ return (3 * u * u * tmp);
+ }
+
+ static double B3(double u)
+ {
+ return (u * u * u);
+ }
+
+ /*
+ * ComputeLeftTangent, ComputeRightTangent, ComputeCenterTangent :
+ *Approximate unit tangents at endpoints and "center" of digitized curve
+ */
+ static Point ComputeLeftTangent(Point[] d, int end)
+ {
+ Point tHat1;
+ tHat1 = d[end + 1] - d[end];
+ tHat1.Normalize();
+ return tHat1;
+ }
+
+ static Point ComputeRightTangent(Point[] d, int end)
+ {
+ Point tHat2;
+ tHat2 = d[end - 1] - d[end];
+ tHat2.Normalize();
+ return tHat2;
+ }
+
+ static Point ComputeCenterTangent(Point[] d, int center)
+ {
+ Point V1, V2, tHatCenter = new Point();
+
+ V1 = d[center - 1] - d[center];
+ V2 = d[center] - d[center + 1];
+ tHatCenter.X = (V1.X + V2.X) / 2.0;
+ tHatCenter.Y = (V1.Y + V2.Y) / 2.0;
+ tHatCenter.Normalize();
+ return tHatCenter;
+ }
+
+
+ /*
+ * ChordLengthParameterize :
+ * Assign parameter values to digitized points
+ * using relative distances between points.
+ */
+ static double[] ChordLengthParameterize(Point[] d, int first, int last)
+ {
+ int i;
+ double[] u = new double[last - first + 1]; /* Parameterization */
+
+ u[0] = 0.0;
+ for (i = first + 1; i <= last; i++)
+ {
+ u[i - first] = u[i - first - 1] + (d[i - 1] - d[i]).Length;
+ }
+
+ for (i = first + 1; i <= last; i++)
+ {
+ u[i - first] = u[i - first] / u[last - first];
+ }
+
+ return u;
+ }
+
+
+
+
+ /*
+ * ComputeMaxError :
+ * Find the maximum squared distance of digitized points
+ * to fitted curve.
+ */
+ static double ComputeMaxError(Point[] d, int first, int last, Point[] bezCurve, double[] u, out int splitPoint)
+ {
+ int i;
+ double maxDist; /* Maximum error */
+ double dist; /* Current error */
+ Point P; /* Point on curve */
+ Point v; /* Vector from point to curve */
+
+ splitPoint = (last - first + 1) / 2;
+ maxDist = 0.0;
+ for (i = first + 1; i < last; i++)
+ {
+ P = BezierII(3, bezCurve, u[i - first]);
+ v = P - d[i];
+ dist = v.LengthSquared;
+ if (dist >= maxDist)
+ {
+ maxDist = dist;
+ splitPoint = i;
+ }
+ }
+ return maxDist;
+ }
+
+ private static double V2Dot(Point a, Point b)
+ {
+ return ((a.X * b.X) + (a.Y * b.Y));
+ }
+
+ }
+
+}