平滑手绘曲线

我有一个允许用户绘制曲线的程序。 但是这些曲线看起来不太好 – 它们看起来摇摇晃晃,手绘。

所以我想要一种能够自动平滑它们的算法。 我知道平滑过程中存在固有的模糊性,因此每次都不会是完美的,但是这些算法确实存在于几个绘图包中并且它们工作得很好。

是否有类似这样的代码示例? C#会很完美,但我可以翻译其他语言。

你可以使用Ramer-Douglas-Peucker算法减少点数, 这里有一个C#实现。 我尝试使用WPFs PolyQuadraticBezierSegment,并根据容差显示出少量改进。

经过一些搜索后,源( 1,2 )似乎表明使用Philip J Schneider的Graphic Gems的曲线拟合算法效果很好, C代码可用 。 几何工具也有一些值得研究的资源。

这是我制作的一个粗略的样本 ,仍然有一些故障,但它在很多时候都很好用。 这是FitCurves.c的快速而脏的 C#端口。 其中一个问题是,如果不减少原始点,则计算出的误差为0,并且它会提前终止,样本会事先使用点减少算法。

/* 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) { Vector 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, Vector tHat1, Vector 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 */ Vector 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, Vector tHat1, Vector tHat2) { int i; Vector[,] A = new Vector[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; Vector 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++) { Vector 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 = ((Vector)d[first + i] - ( ((Vector)d[first] * B0(uPrime[i])) + ( ((Vector)d[first] * B1(uPrime[i])) + ( ((Vector)d[last] * B2(uPrime[i])) + ((Vector)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].X = (Q[i+1].X - Q[i].X) * 3.0; Q1[i].Y = (Q[i+1].Y - Q[i].Y) * 3.0; } /* Generate control vertices for Q'' */ for (i = 0; i <= 1; i++) { Q2[i].X = (Q1[i+1].X - Q1[i].X) * 2.0; Q2[i].Y = (Q1[i+1].Y - Q1[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 - PX) * (Q1_u.X) + (Q_u.Y - PY) * (Q1_u.Y); denominator = (Q1_u.X) * (Q1_u.X) + (Q1_u.Y) * (Q1_u.Y) + (Q_u.X - PX) * (Q2_u.X) + (Q_u.Y - PY) * (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] = 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 Vector ComputeLeftTangent(Point[] d,int end) { Vector tHat1; tHat1 = d[end+1]- d[end]; tHat1.Normalize(); return tHat1; } static Vector ComputeRightTangent(Point[] d,int end) { Vector tHat2; tHat2 = d[end-1] - d[end]; tHat2.Normalize(); return tHat2; } static Vector ComputeCenterTangent(Point[] d,int center) { Vector V1, V2, tHatCenter = new Vector(); 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 */ Vector 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(Vector a,Vector b) { return((aX*bX)+(aY*bY)); } } 

Kris的答案是原始C#的一个非常好的端口,但性能并不理想,并且有些地方浮点不稳定会导致一些问题并返回NaN值(在原始代码中也是如此)。 我创建了一个包含我自己的端口以及Ramer-Douglas-Peuker的库,不仅应该使用WPF点,还应该使用新的支持SIMD的矢量类型和Unity 3D:

https://github.com/burningmime/curves

也许这篇基于WPF + Bezier的文章是一个好的开始: 通过一组带有Bezier基元的2D点绘制平滑曲线

好吧,Kris的工作非常有用。

我意识到他指出的算法因为错误计算的错误以0结尾而早先终止算法的问题是由于重复一个点并且计算的切线是无限的。

我已经根据Kris的代码完成了对Java的翻译,它运行正常我相信:

编辑:

我仍在努力并尝试在算法上获得更好的行为。 我意识到,在非常尖锐的角度上,贝塞尔曲线的表现并不好。 所以我试着将Bezier曲线与Lines结合起来,这就是结果:

 import java.awt.Point; import java.awt.Shape; import java.awt.geom.CubicCurve2D; import java.awt.geom.Line2D; import java.awt.geom.Point2D; import java.util.LinkedList; import java.util.List; import java.util.concurrent.atomic.AtomicInteger; import javax.vecmath.Point2d; import javax.vecmath.Tuple2d; import javax.vecmath.Vector2d; /* An Algorithm for Automatically Fitting Digitized Curves by Philip J. Schneider from "Graphics Gems", Academic Press, 1990 */ public class FitCurves { /* Fit the Bezier curves */ private final static int MAXPOINTS = 10000; private final static double epsilon = 1.0e-6; /** * Rubén: * This is the sensitivity. When it is 1, it will create a line if it is at least as long as the * distance from the previous control point. * When it is greater, it will create less lines, and when it is lower, more lines. * This is based on the previous control point since I believe it is a good indicator of the curvature * where it is coming from, and we don't want long and second derived constant curves to be modeled with * many lines. */ private static final double lineSensitivity=0.75; public interface ResultCurve { public Point2D getStart(); public Point2D getEnd(); public Shape getShape(); } public static class BezierCurve implements ResultCurve { public Point start; public Point end; public Point ctrl1; public Point ctrl2; public BezierCurve(Point2D start, Point2D ctrl1, Point2D ctrl2, Point2D end) { this.start=new Point((int)Math.round(start.getX()), (int)Math.round(start.getY())); this.end=new Point((int)Math.round(end.getX()), (int)Math.round(end.getY())); this.ctrl1=new Point((int)Math.round(ctrl1.getX()), (int)Math.round(ctrl1.getY())); this.ctrl2=new Point((int)Math.round(ctrl2.getX()), (int)Math.round(ctrl2.getY())); if(this.ctrl1.x<=1 || this.ctrl1.y<=1) { throw new IllegalStateException("ctrl1 invalid"); } if(this.ctrl2.x<=1 || this.ctrl2.y<=1) { throw new IllegalStateException("ctrl2 invalid"); } } public Shape getShape() { return new CubicCurve2D.Float(start.x, start.y, ctrl1.x, ctrl1.y, ctrl2.x, ctrl2.y, end.x, end.y); } public Point getStart() { return start; } public Point getEnd() { return end; } } public static class CurveSegment implements ResultCurve { Point2D start; Point2D end; public CurveSegment(Point2D startP, Point2D endP) { this.start=startP; this.end=endP; } public Shape getShape() { return new Line2D.Float(start, end); } public Point2D getStart() { return start; } public Point2D getEnd() { return end; } } public static List FitCurve(double[][] points, double error) { Point[] allPoints=new Point[points.length]; for(int i=0; i < points.length; i++) { allPoints[i]=new Point((int) Math.round(points[i][0]), (int) Math.round(points[i][1])); } return FitCurve(allPoints, error); } public static List FitCurve(Point[] d, double error) { Vector2d tHat1, tHat2; /* Unit tangent vectors at endpoints */ int first=0; int last=d.length-1; tHat1 = ComputeLeftTangent(d, first); tHat2 = ComputeRightTangent(d, last); List result = new LinkedList(); FitCubic(d, first, last, tHat1, tHat2, error, result); return result; } private static void FitCubic(Point[] d, int first, int last, Vector2d tHat1, Vector2d tHat2, double error, List result) { PointE[] bezCurve; /*Control points of fitted Bezier curve*/ double[] u; /* Parameter values for point */ double[] uPrime; /* Improved parameter values */ double maxError; /* Maximum fitting error */ int nPts; /* Number of points in subset */ double iterationError; /*Error below which you try iterating */ int maxIterations = 4; /* Max times to try iterating */ Vector2d tHatCenter; /* Unit tangent vector at splitPoint */ int i; double errorOnLine=error; iterationError = error * error; nPts = last - first + 1; AtomicInteger outputSplitPoint=new AtomicInteger(); /** * Rubén: Here we try to fit the form with a line, and we mark the split point if we find any line with a minimum length */ /* * the minimum distance for a length (so we don't create a very small line, when it could be slightly modeled with the previous Bezier, * will be proportional to the distance of the previous control point of the last Bezier. */ BezierCurve res=null; for(i=result.size()-1; i >0; i--) { ResultCurve thisCurve=result.get(i); if(thisCurve instanceof BezierCurve) { res=(BezierCurve)thisCurve; break; } } Line seg=new Line(d[first], d[last]); int nAcceptableTogether=0; int startPoint=-1; int splitPointTmp=-1; if(Double.isInfinite(seg.getGradient())) { for (i = first; i <= last; i++) { double dist=Math.abs(d[i].xd[first].x); if(dist= minLineLength) { splitPointTmp=i; startPoint=-1; break; } } nAcceptableTogether=0; startPoint=-1; } } } else { //looking for the max squared error for (i = first; i <= last; i++) { Point thisPoint=d[i]; Point2D calculatedP=seg.getByX(thisPoint.getX()); double dist=thisPoint.distance(calculatedP); if(dist= minLineLength) { splitPointTmp=i; startPoint=-1; break; } } nAcceptableTogether=0; startPoint=-1; } } } if(startPoint!=-1) { double minLineLength=Double.POSITIVE_INFINITY; if(res!=null) { minLineLength=lineSensitivity * res.ctrl2.distance(d[startPoint]); } if(d[startPoint].distance(d[last]) >= minLineLength) { splitPointTmp=startPoint; startPoint=-1; } else { nAcceptableTogether=0; } } outputSplitPoint.set(splitPointTmp); if(nAcceptableTogether==(last-first+1)) { //This is a line! System.out.println("line, length: " + d[first].distance(d[last])); result.add(new CurveSegment(d[first], d[last])); return; } /*********************** END of the Line approach, lets try the normal algorithm *******************************************/ if(splitPointTmp < 0) { if(nPts == 2) { double dist = d[first].distance(d[last]) / 3.0; //sqrt((last.x-first.x)^2 + (last.y-first.y)^2) / 3.0 bezCurve = new PointE[4]; bezCurve[0] = new PointE(d[first]); bezCurve[3] = new PointE(d[last]); bezCurve[1]=new PointE(tHat1).scaleAdd(dist, bezCurve[0]); //V2Add(&bezCurve[0], V2Scale(&tHat1, dist), &bezCurve[1]); bezCurve[2]=new PointE(tHat2).scaleAdd(dist, bezCurve[3]); //V2Add(&bezCurve[3], V2Scale(&tHat2, dist), &bezCurve[2]); result.add(new BezierCurve(bezCurve[0],bezCurve[1],bezCurve[2],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, outputSplitPoint); if(maxError < error) { result.add(new BezierCurve(bezCurve[0],bezCurve[1],bezCurve[2],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, outputSplitPoint); if(maxError < error) { result.add(new BezierCurve(bezCurve[0],bezCurve[1],bezCurve[2],bezCurve[3])); return; } u = uPrime; } } } /* Fitting failed -- split at max error point and fit recursively */ tHatCenter = ComputeCenterTangent(d, outputSplitPoint.get()); FitCubic(d, first, outputSplitPoint.get(), tHat1, tHatCenter, error,result); tHatCenter.negate(); FitCubic(d, outputSplitPoint.get(), last, tHatCenter, tHat2, error,result); } //Checked!! static PointE[] GenerateBezier(Point2D[] d, int first, int last, double[] uPrime, Vector2d tHat1, Vector2d tHat2) { int i; Vector2d[][] A = new Vector2d[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; PointE[] bezCurve = new PointE[4]; /* RETURN bezier curve ctl pts */ nPts = last - first + 1; /* Compute the A's */ for (i = 0; i < nPts; i++) { Vector2d v1=new Vector2d(tHat1); Vector2d v2=new Vector2d(tHat2); v1.scale(B1(uPrime[i])); v2.scale(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] += A[i][0].dot(A[i][0]); //C[0][0] += V2Dot(&A[i][0], &A[i][0]); C[0][1] += A[i][0].dot(A[i][1]); //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][0] = C[0][1] C[1][1] += A[i][1].dot(A[i][1]); //C[1][1] += V2Dot(&A[i][1], &A[i][1]); Tuple2d scaleLastB2=new Vector2d(PointE.getPoint2d(d[last])); scaleLastB2.scale(B2(uPrime[i])); // V2ScaleIII(d[last], B2(uPrime[i])) Tuple2d scaleLastB3=new Vector2d(PointE.getPoint2d(d[last])); scaleLastB3.scale(B3(uPrime[i])); // V2ScaleIII(d[last], B3(uPrime[i])) Tuple2d dLastB2B3Sum=new Vector2d(scaleLastB2); dLastB2B3Sum.add(scaleLastB3); //V2AddII(V2ScaleIII(d[last], B2(uPrime[i])), V2ScaleIII(d[last], B3(uPrime[i])) Tuple2d scaleFirstB1=new Vector2d(PointE.getPoint2d(d[first])); scaleFirstB1.scale(B1(uPrime[i])); //V2ScaleIII(d[first], B1(uPrime[i])) Tuple2d sumScaledFirstB1andB2B3=new Vector2d(scaleFirstB1); sumScaledFirstB1andB2B3.add(dLastB2B3Sum); //V2AddII(V2ScaleIII(d[first], B1(uPrime[i])), V2AddII(V2ScaleIII(d[last], B2(uPrime[i])), V2ScaleIII(d[last], B3(uPrime[i]))) Tuple2d scaleFirstB0=new Vector2d(PointE.getPoint2d(d[first])); scaleFirstB0.scale(B0(uPrime[i])); //V2ScaleIII(d[first], B0(uPrime[i]) Tuple2d sumB0Rest=new Vector2d(scaleFirstB0); sumB0Rest.add(sumScaledFirstB1andB2B3); //V2AddII(V2ScaleIII(d[first], B0(uPrime[i])), V2AddII( V2ScaleIII(d[first], B1(uPrime[i])), V2AddII(V2ScaleIII(d[last], B2(uPrime[i])), V2ScaleIII(d[last], B3(uPrime[i])))))); Vector2d tmp=new Vector2d(PointE.getPoint2d(d[first + i])); tmp.sub(sumB0Rest); X[0] += A[i][0].dot(tmp); X[1] += A[i][1].dot(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].distance(d[last]); //(d[first] - d[last]).Length; double epsilonRel = epsilon * segLength; if (alpha_l < epsilonRel || alpha_r < epsilonRel) { /* fall back on standard (probably inaccurate) formula, and subdivide further if needed. */ double dist = segLength / 3.0; bezCurve[0] = new PointE(d[first]); bezCurve[3] = new PointE(d[last]); Vector2d b1Tmp=new Vector2d(tHat1); b1Tmp.scaleAdd(dist, bezCurve[0].getPoint2d()); bezCurve[1] = new PointE(b1Tmp); //(tHat1 * dist) + bezCurve[0]; Vector2d b2Tmp=new Vector2d(tHat2); b2Tmp.scaleAdd(dist, bezCurve[3].getPoint2d()); bezCurve[2] = new PointE(b2Tmp); //(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] = new PointE(d[first]); bezCurve[3] = new PointE(d[last]); Vector2d alphaLTmp=new Vector2d(tHat1); alphaLTmp.scaleAdd(alpha_l, bezCurve[0].getPoint2d()); bezCurve[1] = new PointE(alphaLTmp); //(tHat1 * alpha_l) + bezCurve[0] Vector2d alphaRTmp=new Vector2d(tHat2); alphaRTmp.scaleAdd(alpha_r, bezCurve[3].getPoint2d()); bezCurve[2] = new PointE(alphaRTmp); //(tHat2 * alpha_r) + bezCurve[3]; return (bezCurve); } /* * Reparameterize: * Given set of points and their parameterization, try to find * a better parameterization. * */ static double[] Reparameterize(Point2D[] d,int first,int last,double[] u, Point2D[] 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(Point2D[] Q, Point2D P, double u) { double numerator, denominator; Point2D[] Q1 = new Point2D[3]; //Q' Point2D[] Q2 = new Point2D[2]; //Q'' Point2D 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++) { double qXTmp=(Q[i+1].getX() - Q[i].getX()) * 3.0; //Q1[i].x = (Q[i+1].x - Q[i].x) * 3.0; double qYTmp=(Q[i+1].getY() - Q[i].getY()) * 3.0; //Q1[i].y = (Q[i+1].y - Q[i].y) * 3.0; Q1[i]=new Point2D.Double(qXTmp, qYTmp); } /* Generate control vertices for Q'' */ for (i = 0; i <= 1; i++) { double qXTmp=(Q1[i+1].getX() - Q1[i].getX()) * 2.0; //Q2[i].x = (Q1[i+1].x - Q1[i].x) * 2.0; double qYTmp=(Q1[i+1].getY() - Q1[i].getY()) * 2.0; //Q2[i].y = (Q1[i+1].y - Q1[i].y) * 2.0; Q2[i]=new Point2D.Double(qXTmp, qYTmp); } /* 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.getX() - P.getX()) * (Q1_u.getX()) + (Q_u.getY() - P.getY()) * (Q1_u.getY()); denominator = (Q1_u.getX()) * (Q1_u.getX()) + (Q1_u.getY()) * (Q1_u.getY()) + (Q_u.getX() - P.getX()) * (Q2_u.getX()) + (Q_u.getY() - P.getY()) * (Q2_u.getY()); 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 Point2D BezierII(int degree, Point2D[] V, double t) { int i, j; Point2D Q; /* Point on curve at parameter t */ Point2D[] Vtemp; /* Local copy of control points */ /* Copy array */ Vtemp = new Point2D[degree+1]; for (i = 0; i <= degree; i++) { Vtemp[i] = new Point2D.Double(V[i].getX(), V[i].getY()); } /* Triangle computation */ for (i = 1; i <= degree; i++) { for (j = 0; j <= degree-i; j++) { double tmpX, tmpY; tmpX = (1.0 - t) * Vtemp[j].getX() + t * Vtemp[j+1].getX(); tmpY = (1.0 - t) * Vtemp[j].getY() + t * Vtemp[j+1].getY(); Vtemp[j].setLocation(tmpX, tmpY); } } 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 Vector2d ComputeLeftTangent(Point[] d, int end) { Vector2d tHat1=new Vector2d(PointE.getPoint2d(d[end+1])); tHat1.sub(PointE.getPoint2d(d[end])); tHat1.normalize(); return tHat1; } static Vector2d ComputeRightTangent(Point[] d, int end) { //tHat2 = V2SubII(d[end-1], d[end]); tHat2 = *V2Normalize(&tHat2); Vector2d tHat2=new Vector2d(PointE.getPoint2d(d[end-1])); tHat2.sub(PointE.getPoint2d(d[end])); tHat2.normalize(); return tHat2; } static Vector2d ComputeCenterTangent(Point[] d ,int center) { //V1 = V2SubII(d[center-1], d[center]); Vector2d V1=new Vector2d(PointE.getPoint2d(d[center-1])); V1.sub(new PointE(d[center]).getPoint2d()); //V2 = V2SubII(d[center], d[center+1]); Vector2d V2=new Vector2d(PointE.getPoint2d(d[center])); V2.sub(PointE.getPoint2d(d[center+1])); //tHatCenter.x = (V1.x + V2.x)/2.0; //tHatCenter.y = (V1.y + V2.y)/2.0; //tHatCenter = *V2Normalize(&tHatCenter); Vector2d tHatCenter=new Vector2d((V1.x + V2.x)/2.0, (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].distance(d[i]); } 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(Point2D[] d, int first, int last, Point2D[] bezCurve, double[] u, AtomicInteger splitPoint) { int i; double maxDist; /* Maximum error */ double dist; /* Current error */ Point2D P; /* Point on curve */ Vector2d v; /* Vector from point to curve */ int tmpSplitPoint=(last - first + 1)/2; maxDist = 0.0; for (i = first + 1; i < last; i++) { P = BezierII(3, bezCurve, u[i-first]); v = new Vector2d(P.getX() - d[i].getX(), P.getY() - d[i].getY()); //P - d[i]; dist = v.lengthSquared(); if (dist >= maxDist) { maxDist = dist; tmpSplitPoint=i; } } splitPoint.set(tmpSplitPoint); return maxDist; } /** * This is kind of a bridge between javax.vecmath and java.util.Point2D * @author Ruben * @since 1.24 */ public static class PointE extends Point2D.Double { private static final long serialVersionUID = -1482403817370130793L; public PointE(Tuple2d tup) { super(tup.x, tup.y); } public PointE(Point2D p) { super(p.getX(), p.getY()); } public PointE(double x, double y) { super(x, y); } public PointE scale(double dist) { return new PointE(getX()*dist, getY()*dist); } public PointE scaleAdd(double dist, Point2D sum) { return new PointE(getX()*dist + sum.getX(), getY()*dist + sum.getY()); } public PointE substract(Point2D p) { return new PointE(getX() - p.getX(), getY() - p.getY()); } public Point2d getPoint2d() { return getPoint2d(this); } public static Point2d getPoint2d(Point2D p) { return new Point2d(p.getX(), p.getY()); } } 

这里是后者工作的图像,白色是线条,红色是Bezier:

在此处输入图像描述

使用这种方法,我们使用更少的控制点,更准确。 可以通过lineSensitivity属性调整线条创建的灵敏度。 如果您根本不想使用线条,只需将其设置为无限。

我相信这可以改进。 随意贡献:)

算法没有做任何减少,因为在我的post中首先解释我们必须运行一个。 这是一个DouglasPeuckerReduction实现,对于我来说,在某些情况下,比其他FitCurves更有效地工作(存储点数更少,呈现速度更快)

 import java.awt.Point; import java.util.ArrayList; import java.util.Arrays; import java.util.LinkedList; import java.util.List; public class DouglasPeuckerReduction { public static List reduce(Point[] points, double tolerance) { if (points == null || points.length < 3) return Arrays.asList(points); int firstPoint = 0; int lastPoint = points.length - 1; SortedList pointIndexsToKeep; try { pointIndexsToKeep = new SortedList(LinkedList.class); } catch (Throwable t) { t.printStackTrace(System.out); ErrorReport.process(t); return null; } //Add the first and last index to the keepers pointIndexsToKeep.add(firstPoint); pointIndexsToKeep.add(lastPoint); //The first and the last point cannot be the same while (points[firstPoint].equals(points[lastPoint])) { lastPoint--; } reduce(points, firstPoint, lastPoint, tolerance, pointIndexsToKeep); List returnPoints = new ArrayList(pointIndexsToKeep.size()); for (int pIndex : pointIndexsToKeep) { returnPoints.add(points[pIndex]); } return returnPoints; } private static void reduce(Point[] points, int firstPoint, int lastPoint, double tolerance, List pointIndexsToKeep) { double maxDistance = 0; int indexFarthest = 0; Line tmpLine=new Line(points[firstPoint], points[lastPoint]); for (int index = firstPoint; index < lastPoint; index++) { double distance = tmpLine.getDistanceFrom(points[index]); if (distance > maxDistance) { maxDistance = distance; indexFarthest = index; } } if (maxDistance > tolerance && indexFarthest != 0) { //Add the largest point that exceeds the tolerance pointIndexsToKeep.add(indexFarthest); reduce(points, firstPoint, indexFarthest, tolerance, pointIndexsToKeep); reduce(points, indexFarthest, lastPoint, tolerance, pointIndexsToKeep); } } } 

我在这里使用我自己的SortedList和Line的实现。 抱歉,你必须自己做。

我还没有对它进行过测试,但我想到的一种方法是在某个时间间隔内采样值并创建一个连接点的样条线。

例如,假设曲线的x值从0开始并在10结束。因此,您在x = 1,2,3,4,5,6,7,8,9,10处对y值进行采样并创建样条曲线从点(0,y(0)),(1,y(1)),…(10,y(10))

它可能会出现诸如用户意外尖峰等问题,但可能值得一试

对于Kris回答的Silverlight用户,Point是蹒跚而且Vector不存在。 这是一个支持代码的最小Vector类:

 public class Vector { public double X { get; set; } public double Y { get; set; } public Vector(double x=0, double y=0) { X = x; Y = y; } public static implicit operator Vector(Point b) { return new Vector(bX, bY); } public static Point operator *(Vector left, double right) { return new Point(left.X * right, left.Y * right); } public static Vector operator -(Vector left, Point right) { return new Vector(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; } } } 

还必须解决使用Length和+, – 运算符的问题。 我选择只为FitCurves类添加函数,并重写编译器抱怨的用法。

  public static double Length(Point a, Point b) { double x = aX-bX; double y = aY-bY; return Math.Sqrt(x*x+y*y); } public static Point Add(Point a, Point b) { return new Point(aX + bX, aY + bY); } public static Point Subtract(Point a, Point b) { return new Point(aX - bX, aY - bY); }