source: branches/HeuristicLab.Problems.GrammaticalOptimization/DynamicDataDisplay/Common/BezierBuilder.cs @ 13401

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Windows;

namespace Microsoft.Research.DynamicDataDisplay.Charts.NewLine
{
  // todo check a license
  public static class BezierBuilder
  {
    public static IEnumerable<Point> GetBezierPoints(Point[] points)
    {
      if (points == null)
        throw new ArgumentNullException("points");

      Point[] firstControlPoints;
      Point[] secondControlPoints;

      int n = points.Length - 1;
      if (n < 1)
        throw new ArgumentException("At least two knot points required", "points");
      if (n == 1)
      { // Special case: Bezier curve should be a straight line.
        firstControlPoints = new Point[1];
        // 3P1 = 2P0 + P3
        firstControlPoints[0].X = (2 * points[0].X + points[1].X) / 3;
        firstControlPoints[0].Y = (2 * points[0].Y + points[1].Y) / 3;

        secondControlPoints = new Point[1];
        // P2 = 2P1 – P0
        secondControlPoints[0].X = 2 *
          firstControlPoints[0].X - points[0].X;
        secondControlPoints[0].Y = 2 *
          firstControlPoints[0].Y - points[0].Y;

        return Join(points, firstControlPoints, secondControlPoints);
      }

      // Calculate first Bezier control points
      // Right hand side vector
      double[] rhs = new double[n];

      // Set right hand side X values
      for (int i = 1; i < n - 1; ++i)
        rhs[i] = 4 * points[i].X + 2 * points[i + 1].X;
      rhs[0] = points[0].X + 2 * points[1].X;
      rhs[n - 1] = (8 * points[n - 1].X + points[n].X) / 2.0;
      // Get first control points X-values
      double[] x = GetFirstControlPoints(rhs);

      // Set right hand side Y values
      for (int i = 1; i < n - 1; ++i)
        rhs[i] = 4 * points[i].Y + 2 * points[i + 1].Y;
      rhs[0] = points[0].Y + 2 * points[1].Y;
      rhs[n - 1] = (8 * points[n - 1].Y + points[n].Y) / 2.0;
      // Get first control points Y-values
      double[] y = GetFirstControlPoints(rhs);

      // Fill output arrays.
      firstControlPoints = new Point[n];
      secondControlPoints = new Point[n];
      for (int i = 0; i < n; ++i)
      {
        // First control point
        firstControlPoints[i] = new Point(x[i], y[i]);
        // Second control point
        if (i < n - 1)
          secondControlPoints[i] = new Point(2 * points
            [i + 1].X - x[i + 1], 2 *
            points[i + 1].Y - y[i + 1]);
        else
          secondControlPoints[i] = new Point((points
            [n].X + x[n - 1]) / 2,
            (points[n].Y + y[n - 1]) / 2);
      }

      return Join(points, firstControlPoints, secondControlPoints);
    }

    private static IEnumerable<Point> Join(Point[] points, Point[] firstControlPoints, Point[] secondControlPoints)
    {
      var length = firstControlPoints.Length;
      for (int i = 0; i < length; i++)
      {
        yield return points[i];
        yield return firstControlPoints[i];
        yield return secondControlPoints[i];
      }

      yield return points[length];
    }

    /// <summary>
    /// Solves a tridiagonal system for one of coordinates (x or y)
    /// of first Bezier control points.
    /// </summary>
    /// <param name="rhs">Right hand side vector.</param>
    /// <returns>Solution vector.</returns>
    private static double[] GetFirstControlPoints(double[] rhs)
    {
      int n = rhs.Length;
      double[] x = new double[n]; // Solution vector.
      double[] tmp = new double[n]; // Temp workspace.

      double b = 2.0;
      x[0] = rhs[0] / b;
      for (int i = 1; i < n; i++) // Decomposition and forward substitution.
      {
        tmp[i] = 1 / b;
        b = (i < n - 1 ? 4.0 : 3.5) - tmp[i];
        x[i] = (rhs[i] - x[i - 1]) / b;
      }
      for (int i = 1; i < n; i++)
        x[n - i - 1] -= tmp[n - i] * x[n - i]; // Backsubstitution.

      return x;
    }
  }
}