Geometrie/Racionální algoritmus de Casteljau

Popis

Jestliže položíme

$P_{i}^{0}=P_{i}$ a $w_{i}^{0}=w_{i}$, pak pro každé $i=0,...,n$ platí:

$P_{i}^{k}=(1-u){\frac {w_{i-1}^{k-1}}{w_{i}^{k}}}P_{i-1}^{k-1}+u{\frac {w_{i}^{k-1}}{w_{i}^{k}}}P_{i}^{k-1}$,kde

$w_{i}^{k}=(1-u)w_{i-1}^{k-1}+uw_{i}^{k-1}$

Jak je vidět, jediným rozdílem oproti klasickému algoritmu de Casteljau je fakt, že do výpočtu zahrnujeme váhové parametry, a pro každý nový bod spočítáme jeho poměrnou váhu vzhledem k bodům předešlým.

Pomocí váhových koeficientů lze měnit tvar Bézierovy křivky, s rostoucím váhovým koeficientem se křivka k danému bodu "přibližuje".

Algoritmizace

Syntetické algoritmy nad racionálními bézierovými křivkami jsou téměř totožné jako algoritmy nad obyčejnými bézierovými křivkami, pouze se počítá v o jednu dimenzi vyšším prostoru. V naší implementaci není takový vektor zahrnut, proto se v následujících metodách počítá s váhovou složkou zvlášť.

Nejprve si popišme pomocné metody a třídy:

Vector

třída pro vektor ve 2D,

RacionalRichPoint

třída pro racionální bod a jeho parametrizaci,

Vector.lerp

metoda pro lineární interpolaci dvou bodů,

ctrlPoly

řídící polygon.

public override Vector GetPoint

Funkce GetPoint slouží k výpočtu polohového vektoru příslušného k bodu na křivce pro daný parametr.

public override Vector GetPoint(double s)
{
	int n = ControlPoly.Count;
	double t = (s - StartParam)/(EndParam-StartParam);
	Vector[] points = new Vector[n];
	double[] weights = new double[n];
	for (int i=0;i<n;i++)
	{
		RacionalRichPoint rp = (RacionalRichPoint)ctrlPoly[i];
		weights[i] = rp.Weight;
		points[i] = rp.Locate*rp.Weight;
	}
	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]);
			weights[i] = weights[i]*(1-t)+weights[i+1]*t;
		}
	return points[0]*(1/weights[0]);
}

public override void Split

Rozdělí Bézierovu křivku na dvě v bodě určeném parametrem s. Používá modifikovaný algoritmus de Casteljau pro výpočet bodu na křivce.

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[] points= new Vector[n];
		double[] weights= new double[n];
		Vector[] pointsc1 = new Vector[n];
		Vector[] pointsc2 = new Vector[n];
		double[] weightsc1 = new double[n];
		double[] weightsc2 = new double[n];
		for (int i=0;i<n;i++)
		{
			RacionalRichPoint rp = (RacionalRichPoint)ctrlPoly[i];
			weights[i] = rp.Weight;
			points[i] = rp.Locate*rp.Weight;
		}
		for (int j=0;j<n-1;j++)
		{
			pointsc1[j]= points[0]*(1/weights[0]);
			weightsc1[j]= weights[0];

			pointsc2[n-1-j]= points[n-1-j]*(1/weights[n-1-j]);
			weightsc2[n-1-j]= weights[n-1-j];
			for (int i=0;i<n-1;i++)
			{
				points[i]=points[i].lerp(t,points[i+1]);
				weights[i] = weights[i]*(1-t)+weights[i+1]*t;
			}
		}
		pointsc1[n-1]= points[0]*(1/weights[0]);
		weightsc1[n-1]= weights[0];
		weightsc2[0]= weights[0];
		pointsc2[0]= points[0]*(1/weights[0]);
		c1 = new RationalBezierSynteticly(pointsc1,weightsc1);
		c2 = new RationalBezierSynteticly(pointsc2,weightsc2);
	}
	else
	if (n==1)
	{
		c1 = new RationalBezierSynteticly(ctrlPoly);
		RacionalRichPoint rp = (RacionalRichPoint)ControlPoly[0];
		c2 = new RationalBezierSynteticly(new Vector[]{(Vector)rp.Locate.Clone()},new double[]{rp.Weight});
	}
	else
	{
		c1= new BezierSynteticly(new ArrayList());
		c2= new BezierSynteticly(new ArrayList());
	}
}

public RationalBezierSynteticly IncreaseDegree

Funkce vrací křivku se zvýšeným stupněm.

public void IncreaseDegree()
{
	for (int i=0;i<DegreeIncrement;i++)
	{
		IncreaseDegreeByOne();
	}
}

public void IncreaseDegreeByOne()
{
	double ti=0;
	int n = ControlPoly.Count;
	Vector[] points= new Vector[n];
	double[] weights= new double[n];
	// naklonovani
	RacionalRichPoint theLast = (RacionalRichPoint)ControlPoly[n-1];
	theLast = new RacionalRichPoint((Vector)theLast.Locate.Clone(),0,theLast.Weight);
	for (int i=0;i<n;i++)
	{
		RacionalRichPoint rp = (RacionalRichPoint)ctrlPoly[i];
		weights[i] = rp.Weight;
		points[i] = rp.Locate*rp.Weight;
	}

	for (int i=1;i<n;i++)
	{
		ti = ((double)i)/n;
		RacionalRichPoint rp = (RacionalRichPoint)this.ControlPoly[i];
		rp.Weight = weights[i]*(1-ti)+weights[i-1]*ti;
		rp.Locate = points[i].lerp(ti,points[i-1])*(1/rp.Weight);

	}
	this.ControlPoly.Add(theLast);
}

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.