Alisar una curva dibujada a mano

Tengo un progtwig que permite a los usuarios dibujar curvas. Pero estas curvas no se ven bien, se ven tambaleantes y dibujadas a mano.

Entonces quiero un algoritmo que los suavice automáticamente. Sé que hay ambigüedades inherentes en el proceso de suavizado, por lo que no siempre será perfecto, pero tales algoritmos parecen existir en varios paquetes de dibujo y funcionan bastante bien.

¿Hay ejemplos de código para algo como esto? C # sería perfecto, pero puedo traducir desde otros idiomas.

Puede reducir la cantidad de puntos usando el algoritmo Ramer-Douglas-Peucker. Aquí hay una implementación de C #. Lo probé utilizando WPFs PolyQuadraticBezierSegment y mostró una pequeña cantidad de mejora en función de la tolerancia.

Después de buscar fonts ( 1 , 2 ) que indican que el algoritmo de ajuste de curvas de Gemas de gráficos de Philip J Schneider funciona bien, el código C está disponible . Geometric Tools también tiene algunos recursos que valdría la pena investigar.

Esta es una muestra aproximada que hice, todavía hay algunas fallas técnicas pero funciona bien la mayor parte del tiempo. Aquí está el puerto C # rápido y sucio de FitCurves.c. Uno de los problemas es que si no reduce los puntos originales, el error calculado es 0 y termina antes, la muestra usa el algoritmo de reducción de puntos de antemano.

/* 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)); } } 

La respuesta de Kris es un muy buen puerto del original a C #, pero el rendimiento no es ideal y hay algunos lugares donde la inestabilidad del punto flotante puede causar algunos problemas y devolver valores de NaN (esto también es cierto en el código original). Creé una biblioteca que contiene mi propio puerto así como Ramer-Douglas-Peuker, y debería trabajar no solo con los puntos de WPF, sino también con los nuevos tipos de vectores habilitados para SIMD y Unity 3D:

https://github.com/burningmime/curves

Tal vez este artículo basado en WPF + Bezier sea un buen comienzo: dibuje una curva suave a través de un conjunto de puntos 2D con primitivos Bezier

Bueno, el trabajo de Kris fue muy útil.

Me di cuenta de que el problema que señaló sobre el algoritmo que terminaba antes debido a un error mal calculado que terminaba en 0, se debe al hecho de que un punto se repite y la tangente calculada es infinita.

He hecho una traducción a Java, basada en el código de Kris, funciona bien, creo:

EDITAR:

Todavía estoy trabajando e intentando obtener un mejor comportamiento en el algoritmo. Me di cuenta de que en angularjs muy puntiagudos, las curvas de Bezier simplemente no se comportan bien. Así que traté de combinar las curvas de Bezier con Lines y este es el resultado:

 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()); } } 

Aquí una imagen de este último trabajando, blanco son líneas, y rojo son Bezier:

enter image description here

Usando este enfoque usamos menos puntos de control y más precisos. La sensibilidad para la creación de líneas se puede ajustar mediante el atributo lineSensitivity. Si no quiere que se usen líneas en absoluto, simplemente configúrelo en infinito.

Estoy seguro de que esto se puede mejorar. Siéntase libre de contribuir 🙂

El algoritmo no está haciendo ninguna reducción, y debido a la primera explicación en mi publicación, tenemos que ejecutar uno. Aquí hay una implementación de DouglasPeuckerReduction, que para mí funciona en algunos casos de manera aún más eficiente (menos puntos para almacenar y más rápido de renderizar) que un FitCurves adicional

 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); } } } 

Estoy usando aquí mi propia implementación de una SortedList y de una línea. Tendrás que hacerlo tú mismo, lo siento.

No lo he probado, pero un enfoque que me viene a la mente sería muestrear valores en algún intervalo y crear una spline para conectar los puntos.

Por ejemplo, supongamos que el valor x de tu curva comienza en 0 y termina en 10. Entonces muestras los valores y en x = 1,2,3,4,5,6,7,8,9,10 y creas una spline desde los puntos (0, y (0)), (1, y (1)), … (10, y (10))

Probablemente tenga problemas tales como picos accidentales dibujados por el usuario, pero puede valer la pena intentarlo

Para los usuarios de Silverlight de la respuesta de Kris, Point está obstaculizado y Vector no existe. Esta es una clase Vector mínima que admite el código:

 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; } } } 

También tuvo que abordar el uso de Longitud y +, – operadores. Elegí simplemente agregar funciones a la clase FitCurves y reescribir sus usos donde se quejaba el comstackdor.

  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); } 
    Intereting Posts