Geometrie/Algoritmus de Casteljau

Popis

Pro řadu algoritmů s Bézierovou křivkou (výpočet bodu na křivce, zvýšení stupně Bézierovy křivky, rozdělení Bézierovy křivky atd.) bude základem algoritmus de Casteljau. Ukážeme si základní principy tohoto algoritmu.

Vyjádření

Výchozím algoritmem je rekurentní výpočet Bernsteinových polynomů.

$B_{i}^{n}(u)=(1-u)B_{i}^{n-1}+uB_{i-1}^{n-1}$

Vlastní algoritmus de Casteljau vypadá takto:

$P_{i}^{k}=(1-u)P_{i}^{k-1}+uP_{i+1}^{k-1}$ kde k = 1,...,n; i = 0,...,n-k.

Výpočet bodu na Bézierově křivce

Pro výpočet bodu na Bézierově křivce můžeme použít algoritmus de Casteljau.

Když si výpočet graficky znázorníme, zjistíme, že se ve skutečnosti nejedná o nic jiného, než o postupné dělení úseček řídícího polygonu v zadaném poměru. Počet nově vzniklých bodů se v každém kroku zmenšuje o 1 a ve chvíli, kdy zůstane bod jediný, dostaneme hledaný bod křivky. Bod na Bézierově křivce můžeme rovněž vypočítat přímo pomocí vektorové rovnice Bézierovy křivky, kdy použijeme algoritmus pro výpočet Bernsteinových polynomů

Rozdělení Bézierovy křivky

Pro rozdělení Bézierovy křivky v daném bodě na dvě křivky použijeme modifikovaného algoritmu de Casteljau. Pokud chceme rozdělit Bézierovu kubiku, bude dělícím bodem bod $P_{3}^{3}$ a nově vytvořené kubiky budou určeny řídícími polygony:

$P^{1}(u):P_{0}^{0},P_{1}^{1},P_{2}^{2},P_{3}^{3},P_{4}^{4},P_{5}^{5}$

$P^{2}(u):P_{5}^{5},P_{5}^{4},P_{5}^{3},P_{5}^{2},P_{5}^{1},P_{5}^{0}$

Skládání Bézierových křivek

Podmínkou hladkého spojení Bézierových křivek je zajištění požadované spojitosti (parametrické nebo geometrické) ve společném bodě těchto křivek. Mějme dánu Bézierovu křivku $P=P(u),$ u ∈ <0,1> řídícím polygonem $P_{i}$ (i=0,..,n) a Bézierovu křivku $Q=Q(v),$ v ∈ <0,1> řídícím polygonem $Q_{j}(j=0,...,n)$. Předpokládejme, že platí:

$P(1)=Q(0)$ - počáteční bod křivky $Q(v)$ je koncovým bodem křivky $P(u)$.

Pro tečné vektory $P'(1),Q'(0)$platí:

$P'(1)=m(P_{m}-P_{m-1})$ $Q'(1)=n(Q_{1}-Q_{0})$

Pokud chceme pro Bézierovy křivky ve společném bodě zajistit C1 spojitost (totožné tečné vektory), musí platit:

$n(Q_{1}-Q_{0})=m(P_{m}-P_{m-1})$

Pro $P_{m}=Q_{0}$ dostaneme z této rovnice vztah pro výpočet bodu $Q_{1}$:

$Q_{1}=m/n(P_{m}-P_{m-1})+P_{m}$

Pokud nemají obě křivky stejnou parametrizaci, je výpočet vrcholů navazujícího polygonu (který má být C1 nebo C2 spojitý) závislý na poměru parametrizace: Δ $u_{0}$ : Δ $u_{1}$. Stupeň spojitosti určuje počet definovaných vrcholů nově vytvořeného řídícího polygonu. Pro skládání křivek můžeme rovněž použít upraveného algoritmu de Casteljau.

Zvýšení stupně Bézierovy křivky

Častým technickým požadavkem při modelování křivek bývá přidání dalšího bodu do řídícího polygonu tak, aby se výsledný tvar křivky nezměnil (hovoříme o zvýšení stupně aproximační křivky). Pro vyřešení tohoto problému můžeme použít modifikovaný algoritmus de Casteljau.

Pokud si vrcholy původního řídícího polygonu označíme jako $P_{0}^{0},...,P_{n}^{0}$, pak pro vrcholy nového řídícího polygonu platí:

$P_{i}^{1}=u_{i}P_{i-1}^{0}+(1-u_{i})P_{i}^{0}$

kde

$u_{i}=i/(n+1),i=0,...,n+1$

Při bližším pohledu na uvedený postup zjistíme, že neděláme nic jiného, než že původní parametr $u$ jakoby „rozsekáme“ na $n.(n+1)$ částí, kdy původní body jsou umístěny vždy v násobcích $(n+1)$ a nové body umístíme do násobků $n$. Graficky vzato, každou úsečku původního řídícího polygonu rozdělíme na $(n+1)$ stejných částí a pak po každých $n$ dílcích umístíme bod nový.

Algoritmizace

Vrací bod na Bézierově křivce v obecném parametru s

public override Vector GetPoint(double s)
{
	int n = ControlPoly.Count;
	double t = (s - StartParam)/(EndParam-StartParam);
	Vector[] points = new Vector[n];
	for (int i=0;i<n;i++)
	{
		points[i] = ((RichPoint)(ctrlPoly[i])).Locate;
	}
	for (int j=0;j<n-1;j++)
		for (int i=0;i<n-1;i++)
			points[i]=points[i].lerp(t,points[i+1]);
	return points[0];
}

Rozdělí Bézierovu křivku na dvě v bodě určeném prametrem s

// s - parametr, ve kterém se má bod rozdělit
// c1 - první nová křivka
// c2 - druhá nová křivka

public override void Split(double s, out Curve c1, out Curve c2)
{
	int n = ControlPoly.Count;
	if (n>=2)
	{
		double t = (s- StartParam)/(EndParam-StartParam);
		Vector[] pointsc1 = new Vector[n];
		Vector[] pointsc2 = new Vector[n];
		Vector[] points= new Vector[n];
		for (int i=0;i<n;i++)
		{
			points[i] = ((RichPoint)(ctrlPoly[i])).Locate;
		}
		for (int j=0;j<n-1;j++)
		{
			pointsc1[j]= points[0];
			pointsc2[n-1-j]= points[n-1-j];
			for (int i=0;i<n-1;i++)
			{
				points[i]=points[i].lerp(t,points[i+1]);
			}
		}
		pointsc1[n-1]= points[0];
		pointsc2[0]= (Vector)points[0].Clone();
		c1 = new BezierSynteticly(pointsc1);
		c2 = new BezierSynteticly(pointsc2);
	}
	else if (n==1)
	{
		c1 = new BezierSynteticly(ctrlPoly);
		Vector v = (Vector)((RichPoint)ControlPoly[0]).Locate.Clone();
		c2 = new BezierSynteticly(new Vector[]{v});
	}
	else
	{
		c1= new BezierSynteticly(new Vector[]{});
		c2= new BezierSynteticly(new Vector[]{});
	}
}

Vrátí Bezierovu křivku o stupeň vyšší

public void IncreaseDegreeByOne()
{
	double ti=0;
	int n = ControlPoly.Count;
	Vector[] points= new Vector[n];
	for (int i=0;i<n;i++)
	{
		points[i] = ((RichPoint)(ctrlPoly[i])).Locate;
	}
	for (int i=1;i<n;i++)
	{
		ti = ((double)i)/n;
		RichPoint rp = (RichPoint)this.ControlPoly[i];
		rp.Locate = points[i].lerp(ti,points[i-1]);
	}
	this.ControlPoly.Add(new RichPoint(points[n-1]));
}

public override RichPoint CreateRichPoint()
{
	return new RichPoint(new Vector(0,0));
}

Autoři

Tento text vypracovali studenti Univerzity Palackého v Olomouci katedry Matematické informatiky jako zápočtový úkol do předmětu Počítačová geometrie.