/Main/src/DynamicDataDisplay/Common/BezierBuilder.cs

  1using System;
  2using System.Collections.Generic;
  3using System.Linq;
  4using System.Text;
  5using System.Windows;
  6
  7namespace Microsoft.Research.DynamicDataDisplay.Charts.NewLine
  8{
  9	// todo check a license
 10	public static class BezierBuilder
 11	{
 12		public static IEnumerable<Point> GetBezierPoints(Point[] points)
 13		{
 14			if (points == null)
 15				throw new ArgumentNullException("points");
 16
 17			Point[] firstControlPoints;
 18			Point[] secondControlPoints;
 19
 20			int n = points.Length - 1;
 21			if (n < 1)
 22				throw new ArgumentException("At least two knot points required", "points");
 23			if (n == 1)
 24			{ // Special case: Bezier curve should be a straight line.
 25				firstControlPoints = new Point[1];
 26				// 3P1 = 2P0 + P3
 27				firstControlPoints[0].X = (2 * points[0].X + points[1].X) / 3;
 28				firstControlPoints[0].Y = (2 * points[0].Y + points[1].Y) / 3;
 29
 30				secondControlPoints = new Point[1];
 31				// P2 = 2P1 – P0
 32				secondControlPoints[0].X = 2 *
 33					firstControlPoints[0].X - points[0].X;
 34				secondControlPoints[0].Y = 2 *
 35					firstControlPoints[0].Y - points[0].Y;
 36
 37				return Join(points, firstControlPoints, secondControlPoints);
 38			}
 39
 40			// Calculate first Bezier control points
 41			// Right hand side vector
 42			double[] rhs = new double[n];
 43
 44			// Set right hand side X values
 45			for (int i = 1; i < n - 1; ++i)
 46				rhs[i] = 4 * points[i].X + 2 * points[i + 1].X;
 47			rhs[0] = points[0].X + 2 * points[1].X;
 48			rhs[n - 1] = (8 * points[n - 1].X + points[n].X) / 2.0;
 49			// Get first control points X-values
 50			double[] x = GetFirstControlPoints(rhs);
 51
 52			// Set right hand side Y values
 53			for (int i = 1; i < n - 1; ++i)
 54				rhs[i] = 4 * points[i].Y + 2 * points[i + 1].Y;
 55			rhs[0] = points[0].Y + 2 * points[1].Y;
 56			rhs[n - 1] = (8 * points[n - 1].Y + points[n].Y) / 2.0;
 57			// Get first control points Y-values
 58			double[] y = GetFirstControlPoints(rhs);
 59
 60			// Fill output arrays.
 61			firstControlPoints = new Point[n];
 62			secondControlPoints = new Point[n];
 63			for (int i = 0; i < n; ++i)
 64			{
 65				// First control point
 66				firstControlPoints[i] = new Point(x[i], y[i]);
 67				// Second control point
 68				if (i < n - 1)
 69					secondControlPoints[i] = new Point(2 * points
 70						[i + 1].X - x[i + 1], 2 *
 71						points[i + 1].Y - y[i + 1]);
 72				else
 73					secondControlPoints[i] = new Point((points
 74						[n].X + x[n - 1]) / 2,
 75						(points[n].Y + y[n - 1]) / 2);
 76			}
 77
 78			return Join(points, firstControlPoints, secondControlPoints);
 79		}
 80
 81		private static IEnumerable<Point> Join(Point[] points, Point[] firstControlPoints, Point[] secondControlPoints)
 82		{
 83			var length = firstControlPoints.Length;
 84			for (int i = 0; i < length; i++)
 85			{
 86				yield return points[i];
 87				yield return firstControlPoints[i];
 88				yield return secondControlPoints[i];
 89			}
 90
 91			yield return points[length];
 92		}
 93
 94		/// <summary>
 95		/// Solves a tridiagonal system for one of coordinates (x or y)
 96		/// of first Bezier control points.
 97		/// </summary>
 98		/// <param name="rhs">Right hand side vector.</param>
 99		/// <returns>Solution vector.</returns>
100		private static double[] GetFirstControlPoints(double[] rhs)
101		{
102			int n = rhs.Length;
103			double[] x = new double[n]; // Solution vector.
104			double[] tmp = new double[n]; // Temp workspace.
105
106			double b = 2.0;
107			x[0] = rhs[0] / b;
108			for (int i = 1; i < n; i++) // Decomposition and forward substitution.
109			{
110				tmp[i] = 1 / b;
111				b = (i < n - 1 ? 4.0 : 3.5) - tmp[i];
112				x[i] = (rhs[i] - x[i - 1]) / b;
113			}
114			for (int i = 1; i < n; i++)
115				x[n - i - 1] -= tmp[n - i] * x[n - i]; // Backsubstitution.
116
117			return x;
118		}
119	}
120}