【问题标题】:Smoothing a hand-drawn curve平滑手绘曲线
【发布时间】:2023-03-09 02:45:01
【问题描述】:

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

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

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

【问题讨论】:

  • 看看这个 SO 帖子(类似的平滑问题,尽管使用 GDI+..不确定您使用的是 OpenGL 还是 GDI+):http://stackoverflow.com/questions/4554698/c-gdi-curve-drawing-issue
  • 我认为这是一个不同的问题。他们似乎在谈论 GDI 绘图例程的一些神器。我的问题仍然适用于所使用的图形平台:如果你在一张纸上手绘曲线,它仍然看起来不太好。我想我需要一些可以智能地选择的东西,比如说,十几个贝塞尔点,并通过这些点拟合一条样条曲线,以平滑近似原始曲线。 (顺便说一句,我正在使用 WPF。)

标签: c# bezier spline smoothing


【解决方案1】:

您可以使用Ramer–Douglas–Peucker algorithm 减少点数,有一个C# 实现here。我使用 WPFs PolyQuadraticBezierSegment 进行了尝试,它显示出少量的改进,具体取决于公差。

经过一番搜索后,来源 (1,2) 似乎表明使用 Philip J Schneider 的 Graphic Gems 的曲线拟合算法效果很好,the C code is availableGeometric Tools 也有一些值得研究的资源。

这是我做的一个粗略的sample,仍然有一些小故障,但它在很多时候都很好用。这是 FitCurves.c 的 quick and dirty 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<Point> FitCurve(Point[] d, double error)
    {
        Vector tHat1, tHat2;    /*  Unit tangent vectors at endpoints */

        tHat1 = ComputeLeftTangent(d, 0);
        tHat2 = ComputeRightTangent(d, d.Length - 1);
        List<Point> result = new List<Point>();
        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<Point> 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 - 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] = 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((a.X*b.X)+(a.Y*b.Y));
    }

}

【讨论】:

  • 非常感谢 - 这个组合看起来最有前途。参考 1 看起来特别有趣——您实际上可以看到它在现实生活中发挥作用。视频说代码可用,但我还没有真正找到它,所以我想我会自己开始实现这些算法......
  • 我添加了一个基于我链接到的 FitCurves.c 的示例。
  • 克里斯,这真是太棒了。谢谢。它确实在 90% 的时间里非常有效。奇怪的是,另外 10% 让事情变得更糟——某个大的随机峰值出现在某处。关于这是否是一个可能的错误,还是算法中更固有的东西有什么想法?
  • 另外,人们也希望允许闭合曲线,但这里的算法似乎移动端点。 (画一条垂直线以清楚地看到这一点)。您认为约束端点位置和切线是封闭曲线的正确方法吗?
  • 这里讨论了 Graphics Gems 代码的数值稳定性,以及如何在边缘情况下改进它:newsgroups.derkeiler.com/Archive/Comp/comp.graphics.algorithms/…>。
【解决方案2】:

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

https://github.com/burningmime/curves

【讨论】:

    【解决方案3】:

    也许这篇基于 WPF+Bezier 的文章是一个好的开始:Draw a Smooth Curve through a Set of 2D Points with Bezier Primitives

    【讨论】:

    • 可能非常有用,但这只是成功的一半。你不会想要使用所有的点,因为有些只是摇摆不定,所以你需要选择一个智能子集。也许使用上面@Kris 提到的 Ramer-Douglas-Peucker 算法,然后在减少的点集上使用它会起作用。我会试试的。
    • 关于 Ramer-Douglas-Peucker 算法的快速说明;它显然不适用于自相交曲线。画一个循环并应用它。
    【解决方案4】:

    嗯,克里斯的工作很有用。

    我意识到他指出的算法由于计算错误以 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<ResultCurve> 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<ResultCurve> 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<ResultCurve> result = new LinkedList<ResultCurve>();
        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<ResultCurve> 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].x-d[first].x);
                if(dist<errorOnLine) {
                    nAcceptableTogether++;
                    if(startPoint==-1) startPoint=i;
                } else {
                    if(startPoint!=-1) {
                        double minLineLength=Double.POSITIVE_INFINITY;
                        if(res!=null) {
                            minLineLength=lineSensitivity * res.ctrl2.distance(d[startPoint]);
                        }
                        double thisFromStart=d[startPoint].distance(d[i]);
                        if(thisFromStart >= 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<errorOnLine) {
                    nAcceptableTogether++;
                    if(startPoint==-1) startPoint=i;
                } else {
                    if(startPoint!=-1) {
                        double thisFromStart=d[startPoint].distance(thisPoint);
                        double minLineLength=Double.POSITIVE_INFINITY;
                        if(res!=null) {
                            minLineLength=lineSensitivity * res.ctrl2.distance(d[startPoint]);
                        }
                        if(thisFromStart >= 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());
            }
        }
    

    这是后者工作的图像,白色是线条,红色是贝塞尔曲线:

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

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

    算法没有做任何减少,由于我的帖子中解释的第一个,我们必须运行一个。这是一个 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<Point> reduce(Point[] points, double tolerance)
        {
            if (points == null || points.length < 3) return Arrays.asList(points);
    
            int firstPoint = 0;
            int lastPoint = points.length - 1;
    
            SortedList<Integer> pointIndexsToKeep;
            try {
                pointIndexsToKeep = new SortedList<Integer>(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<Point> returnPoints = new ArrayList<Point>(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<Integer> 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 实现。对不起,你必须自己做。

    【讨论】:

      【解决方案5】:

      我尚未对其进行测试,但我想到的一种方法是在某个间隔内对值进行采样并创建一条样条线来连接这些点。

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

      可能会出现用户意外画出尖刺等问题,但可能值得一试

      【讨论】:

      • 是的,我认为问题是意外尖峰。 (另外我认为参数化曲线会更通用,即 (x(t), y(t)) 而不是 y(x),以允许闭合曲线等。
      • 只是定期采样,如果您碰巧选择了摆动上的样条点,您可能会使它比原来更糟。看来我们真正想要的可能是将样条点放置在“模糊”高斯外曲率的极值和拐点处。
      【解决方案6】:

      对于 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(b.X, b.Y);
         }
      
         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 = a.X-b.X;
           double y = a.Y-b.Y;
           return Math.Sqrt(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);
        }
      

      【讨论】:

        猜你喜欢
        • 1970-01-01
        • 1970-01-01
        • 2012-06-09
        • 1970-01-01
        • 2017-06-10
        • 2017-01-14
        • 2015-11-18
        • 1970-01-01
        • 2012-01-31
        相关资源
        最近更新 更多