Click here to Skip to main content
15,886,362 members
Articles / Programming Languages / C# 4.0
Tip/Trick

Polynomial class with FindRoots method

Rate me:
Please Sign up or sign in to vote.
5.00/5 (9 votes)
15 Jan 2015CPOL2 min read 21.7K   341   12   12
All the tools to calculate Distance to Bezier curve, find the root of Polynomial, do Complex math

Image 1

Introduction

I needed to find the distance from a Point to a Bezier Segment. Turns out you need lot of Math tools.  A polynomial class which can do simple math, derivate and find its roots. As well as Complex arithmentics.

The source can also be found on GitHub:
https://github.com/superlloyd/Poly

Using the code

There is a Complex.cs class which supports simple arithmetics. ex:

C#
var c1 = 1 + 2 * Complex.I;
var c2 = c1 * c1;
var conjugate = !c1;
Assert.AreEquals(c2 / c1, c1);

 

There is a Polynomial.cs class which supports simple arithmetics, which is defined by listing its coefficient

C#
var X = new Polynomial(0, 1); // poly(x) = X
var x2mx = X * (1 - X) + 1; // poly(x) = 1 + X - X^2
var x3 = new Polynomial(1,0,-1,2); // poly(x) = 1-x^2+x^3

It also support Pow() Derivate() and Integrate()

C#
var Xp1 = 1 + Polynomial.X();
var X2 = Xp1.Pow(2);
Assert.AreEquals(X1, new Polynomial(1,2,1))

// -- other method
var X = new Polynomial(0, 1);
var X2 = 1 + X + (X^3); // unfortunately operator priority is low on '^' hence the parenthesis
Assert.AreEquals(X2.Derivate(), 1+3*(X^2));

// other method
var X = new Polynomial(0, 1);
var X2 = 1 + X;
Assert.AreEquals(X2.Integrate(), X + (X^2)/2);

It also have its most interesting method: FindRoots()

C#
public Complex[] FindRoots(int maxIteration = 10)

Which uses the Durand-Kerner algorithm:
http://en.wikipedia.org/wiki/Durand%E2%80%93Kerner_method

It's a very simple and quick algorithm which find all the root at once. The implementation is only 55 lines long!
It will be used later to calculate distance to a Bezier curves

 

I also added the Bezier curves equations:
http://en.wikipedia.org/wiki/B%C3%A9zier_curve

C#
public static Polynomial LinearBezierCurve(double p0, double p1)
{
 var T = new Polynomial(0, 1);
 return p0 + T * (p1 - p0);
}
public static Polynomial QuadraticBezierCurve(double p0, double p1, double p2)
{
 var T = new Polynomial(0, 1);
 return (1 - T) * LinearBezierCurve(p0, p1) + T * LinearBezierCurve(p1, p2);
}
public static Polynomial CubicBezierCurve(double p0, double p1, double p2, double p3)
{
 var T = new Polynomial(0, 1);
 return (1 - T) * QuadraticBezierCurve(p0, p1, p2) + T * QuadraticBezierCurve(p1, p2, p3);
}

Distance to Bezier Segment

Now that I have all these basic math tools I can finally find the distance from a point to a Bezier segment!

If the curve equations is: BC(t),
The distance from P to the curve is the minimal distance of P to all the points BC(t).
To find the mimums of the distance I need to find the root / 0 of D'(t) (of the derivate of D(t))
Also t is between 0,1 for this reason I should also consider the segment extremities.

Which leave me with this very simple C# implementation of the distance:

C#
public static double DistanceToBezier(this Point p, Point start, Point cp1, Point cp2, Point end)
{
 var BX = Polynomial.CubicBezierCurve(start.X, cp1.X, cp2.X, end.X);
 var BY = Polynomial.CubicBezierCurve(start.Y, cp1.Y, cp2.Y, end.Y);
 var dsquare = (BX - p.X) * (BX - p.X) + (BY - p.Y) * (BY - p.Y);
 var deriv = dsquare.Derivate().Normalize();
 var deriveRoots = deriv.FindRoots();
 return deriveRoots
  .Where(x => Math.Abs(x.Imaginary) < Polynomial.Epsilon && x.Real > 0 && x.Real < 1)
  .Select(x => (double)x.Real)
  .Concat(new double[] { 0, 1 })
  .Select(x => Math.Sqrt(dsquare.Compute(x)))
  .OrderBy(x => x)
  .First();
}

Points of Interest

Bezier curve and Polynomial.FindRoots() are implemented in a very clear and simple fashion. Thanks Wikipedia!

The code code with a Utility libray (including the unit tests) and a WPF application which show live calculated point and calculated distance.

I should thanks Jeremy Alles, as I used his application as the basis of my visual test application.

History

First release fo now...
Code change on GitHub: https://github.com/superlloyd/Poly

License

This article, along with any associated source code and files, is licensed under The Code Project Open License (CPOL)


Written By
Software Developer (Senior) http://www.ansibleww.com.au
Australia Australia
The Australia born French man who went back to Australia later in life...
Finally got over life long (and mostly hopeless usually, yay!) chronic sicknesses.
Worked in Sydney, Brisbane, Darwin, Billinudgel, Darwin and Melbourne.

Comments and Discussions

 
QuestionA question on the problem Pin
Kenneth Haugland15-Jan-15 6:53
mvaKenneth Haugland15-Jan-15 6:53 
AnswerRe: A question on the problem Pin
Super Lloyd15-Jan-15 7:01
Super Lloyd15-Jan-15 7:01 
GeneralRe: A question on the problem Pin
Kenneth Haugland15-Jan-15 7:10
mvaKenneth Haugland15-Jan-15 7:10 
GeneralRe: A question on the problem Pin
Super Lloyd15-Jan-15 17:20
Super Lloyd15-Jan-15 17:20 
GeneralRe: A question on the problem Pin
Kenneth Haugland15-Jan-15 22:45
mvaKenneth Haugland15-Jan-15 22:45 
QuestionDownloads broken Pin
Kenneth Haugland15-Jan-15 6:24
mvaKenneth Haugland15-Jan-15 6:24 
AnswerRe: Downloads broken Pin
Super Lloyd15-Jan-15 6:36
Super Lloyd15-Jan-15 6:36 
GeneralRe: Downloads broken Pin
Kenneth Haugland15-Jan-15 6:43
mvaKenneth Haugland15-Jan-15 6:43 
GeneralRe: Downloads broken Pin
Super Lloyd15-Jan-15 6:50
Super Lloyd15-Jan-15 6:50 
QuestionWhat is wrong with System.Numerics.Complex? Pin
Daniel Pfeffer15-Jan-15 6:23
professionalDaniel Pfeffer15-Jan-15 6:23 
AnswerRe: What is wrong with System.Numerics.Complex? Pin
Super Lloyd15-Jan-15 6:38
Super Lloyd15-Jan-15 6:38 
AnswerRe: What is wrong with System.Numerics.Complex? Pin
Super Lloyd15-Jan-15 6:55
Super Lloyd15-Jan-15 6:55 

General General    News News    Suggestion Suggestion    Question Question    Bug Bug    Answer Answer    Joke Joke    Praise Praise    Rant Rant    Admin Admin   

Use Ctrl+Left/Right to switch messages, Ctrl+Up/Down to switch threads, Ctrl+Shift+Left/Right to switch pages.