Approximate cubic polynomial with polyline

Cubic polynomial is often used as main building block in some or other algorithms of piecewise function approximation. Often we then want to draw the result of the approximation on screen or printer. But we lack primitives to do it.

To draw the cubic polynomial we have to approximate it with the polyline. The maximum distance from polynomial to polyline should not exceed some "tolerance" value.

Let we have a cubic polynomial defined at [x1, x2] interval.

To approximate it with polyline we should do the following:

  1. Get the deviation polynomial, i.e. the difference between the initial cubic polynomial and the straight line passing through its left and right bound points. This polynomial is either identically equal to zero or has one or two extremum(s) at [x1, x2].
  2. Evaluate the values of deviation polynomial at extremum points. It its absolute values are lower than the tolerance then the initial cubic polynomial can be approximated with a straight line passing through points (x1, y1) and (x2, y2). Otherwise
  3. Split the initial interval  [x1, x2] on two or three subintervals (depending on extremum’s count) and repeat the procedure recursively from (1) for each of subintervals.

The code

/// <summary>
/// Approximating Cubic Polynomial with Polyline.
/// </summary>
public static class CubicPolynomialPolylineApproximation
{
    /// <summary>
    /// Gets the approximation of the polynomial with polyline.
    /// </summary>
    /// <param name="polynomial">The polynomial.</param>
    /// <param name="x1">The abscissas start.</param>
    /// <param name="x2">The abscissas stop.</param>
    /// <param name="tolerance">The tolerance is the maximum distance from the cubic
    /// polynomial to the approximating polyline.</param>
    /// <returns></returns>
    public static Collection<Point> Approximate(Polynomial polynomial
      , double x1, double x2, double tolerance)
    {
        Debug.Assert(x1 <= x2, "x1 <= x2");
        Debug.Assert(polynomial.Order == 3, "polynomial.Order == 3");

        Collection<Point> points = new Collection<Point>();

        // Get difference between given polynomial and the straight line passing its node points.
        Polynomial deviation = DeviationPolynomial(polynomial, x1, x2);
        Debug.Assert(deviation.Order == 3, "diff.Order == 3");
        if (deviation[0] == 0 && deviation[1] == 0 && deviation[2] == 0 && deviation[3] == 0)
        {
            points.Add(new Point(x1, polynomial.GetValue(x1)));
            points.Add(new Point(x2, polynomial.GetValue(x2)));
            return points;
        }

        // Get previous polynomial first derivative
        Polynomial firstDerivative = new Polynomial(new double[] { deviation[1], 2 * deviation[2], 3 * deviation[3] });

        // Difference polynomial extremum’s.
        // Find first derivative roots.
        Complex[] complexRoots = firstDerivative.Solve();
        // Get real roots in [x1, x2].
        List<double> roots = new List<double>();
        foreach (Complex complexRoot in complexRoots)
        {
            if (complexRoot.Imaginary == 0)
            {
                double r = complexRoot.Real;
                if (r > x1 && r < x2)
                    roots.Add(r);
            }
        }
        Debug.Assert(roots.Count > 0, "roots.Count > 0");
        Debug.Assert(roots.Count <= 2, "roots.Count <= 2");

        // Check difference polynomial extremal values.
        bool approximates = true;
        foreach (double x in roots)
        {
            if (Math.Abs(deviation.GetValue(x)) > tolerance)
            {
                approximates = false;
                break;
            }
        }
        if (approximates)
        {// Approximation is good enough.
            points.Add(new Point(x1, polynomial.GetValue(x1)));
            points.Add(new Point(x2, polynomial.GetValue(x2)));
            return points;
        }

        if (roots.Count == 2)
        {
            if (roots[0] == roots[1])
                roots.RemoveAt(1);
            else if (roots[0] > roots[1])
            {// Sort the roots
                // Swap roots
                double x = roots[0];
                roots[0] = roots[1];
                roots[1] = x;
            }
        }
        // Add the end abscissas.
        roots.Add(x2);

        // First subinterval.
        Collection<Point> pts = Approximate(polynomial, x1, roots[0], tolerance);
        // Copy all points.
        foreach (Point pt in pts)
        {
            points.Add(pt);
        }
        // The remnant of subintervals.
        for (int i = 0; i < roots.Count – 1; ++i)
        {
            pts = Approximate(polynomial, roots[i], roots[i + 1], tolerance);
            // Copy all points but the first one.
            for (int j = 1; j < pts.Count; ++j)
            {
                points.Add(pts[j]);
            }
        }
        return points;
    }

    /// <summary>
    /// Gets the difference between given polynomial and the straight line
    /// passing through its node points.
    /// </summary>
    /// <param name="polynomial">The polynomial.</param>
    /// <param name="x1">The abscissas start.</param>
    /// <param name="x2">The abscissas stop.</param>
    /// <returns></returns>
    static Polynomial DeviationPolynomial(Polynomial polynomial, double x1, double x2)
    {
        double y1 = polynomial.GetValue(x1);
        double y2 = polynomial.GetValue(x2);
        double a = (y2 – y1) / (x2 – x1);
        double b = y1 – a * x1;
        if (a != 0)
            return polynomial.Subtract(new Polynomial(new double[] { b, a }));
        else if (b != 0)
            return polynomial.Subtract(new Polynomial(new double[] { b }));
        else
            return polynomial;
    }

}

In the code above I’m using the helper class Polynomial encapsulating operations on polynomials including addition, subtraction, dividing, root finding, etc. It’s ported from "Numerical Recipes in C, 2-nd Edition" book with some additions and bug fixes.

You can download the code sample here:

http://cid-39d56f0c7a08d703.skydrive.live.com/embedrowdetail.aspx/.Public/CubicPolynomialPolylineApproximation

It is Visual Studio 2008 solution targeted to .NET 3.5. It contains WPF Windows Application project designed to demonstrate drawing of some cubic polynomials. You can select one of the curves from Combo Box at the top of the Window, experiment with tolerance and set appropriate XY Scales.

Advertisements

About ovpwp

I am engaged in programming, maintenance and supply of computing technique since 1970. Started with the computers M, BESM, Minsk series and received appropriate education in the least, in which it was then possible in THE USSR. Then he went the usual way - ES series, IBM 360/370, Sun Spark, Ibm Power & PS2, PC. Programming started in the code (machine commands), then were sorts of assemblers, Algol, FORTRAN, PL/1, C, C , VB, C#. It is only the ones that I used the production scale; but there were, of course, others like List, Modula, Pascal, Java, etc. Currently I prefer .NET platform for desktop development. I don't really like web-programming (probably because of the inability for quality design), but I have enough experience in site building in LAMP environment using PHP.
This entry was posted in Uncategorized. Bookmark the permalink.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s