Geometrie/C2 kubická interpolace

Z Wikiknih

Popis[editovat | editovat zdroj]

Splajny byly studovány už začátkem 20. století. Označení spline použil poprvé v roce 1946 Shoenberg. Původně byla tímto názvem označována speciální pružná dřevěná nebo kovová pravítka, kterých používali konstruktéři, když při vytváření návrhů trupů lodí potřebovali zadanými body proložit hladkou křivku.

C2 kubická interpolace patří do skupiny interpolací křivek po obloucích, tj. každý úsek mezi dvěma opěrnými body se interpoluje zvlášť. Pro C2 interpolační křivku musí být zabezpečena C2 spojitost, tj. až do řádu dva.

Spline křivka řádu k
pro výpočet interpolačních křivek se nejčastěji používají spline křivky. Spline křivka je křivka složená z polynomu stupně k, přičemž v opěrných bodech je zajištěna spojitost až do (k – 1) derivací. V praxi se nejčastěji používají kubické spline křivky.

Vyjádření[editovat | editovat zdroj]

Kubická spline křivka
Mějme dáno (n + 1) opěrných bodů P0,...,Pn a (n + 1) reálných parametrů u0 <...< un. Interpolační kubická spline křivka P(u) je složena z kubických polynomů Pi(u), které musí splňovat tyto podmínky:

Volbou čísel u0,...,un se určuje parametrizace spline křivky. Pro parametrizaci obvykle volíme ui=i (Uniformní parametrizace).

Pro výpočet kubické spline křivky lze použít kubické polynomy ve tvaru:

Tyto polynomy dosadíme do podmínek pro spline křivku a vyřešíme soustavu 4 × n rovnic.

Spline křivka je určena 4 × n koeficienty ai , bi , ci a di. Ze zadání (n + 1) opěrných bodů a z podmínek pro kubický spline však dostáváme pouze (4×n - 2) podmínek. Aby byla soustava jednoznačně řešitelná, musíme ještě určit další dvě podmínky. Nejčastěji se volí jedna z následujících možností:

  1. Vektory prvních derivací v krajních bodech 0 a n.
  2. Vektory druhých derivací v krajních bodech P´´0 a P´´n. Pokud použijete P´´0=0 a P´´n=0 hovoříme o přirozené spline křivce.
  3. Křivka je uzavřená.

Při výpočtu kubické spline křivky se obvykle postupuje tak, že nejdříve vypočteme tečné vektory 0,...,P´n. Jednotlivé kubiky Pi(u) pak můžeme vyjádřit pomocí Hermitovy interpolace.

Jako ukázku implementace jsme zvolili druhou možnost (přirozený kubický splajn). Z podmínek P´´(u 0)=0 a P´´(u n)=0 dostaneme soustavu (n+1) rovnic pro výpočet tečných vektorů 0,...,P´n. Zvolíme-li uniformní parametrizaci, pak platí:

Vlastnosti[editovat | editovat zdroj]

U C2 kubických splajnů je zajištěna spojitá derivace až do stupně 2. Používají se polynomy lichého stupně z důvodu vzniku pouze pásové matice s kterou se podstatně lépe pracuje, ale hlavně kvůli složitosti výpočtu.

Algoritmizace[editovat | editovat zdroj]

Uveďme si nejdříve nějaké poznámky k proměnným a použitým metodám.

Vector
třída popisující vektor,
Richpoint
třída popisující bod,
ctrlPoly
body řídícího polygonu.

A nyní již k vlastnímu algoritmu.

Metoda Vector GetPoint(double t)[editovat | editovat zdroj]

K zadanému parametru vrátí polohový vektor bodu na křivce.

public override Vector GetPoint(double t)
{
	// preventivně setřídím body
	ctrlPoly.Sort();
	// definice krajních bodů a tečných vektorů
	RichPoint b1, b2;
	Vector db1, db2;
	b1 = b2 = null;
	db1 = db2 = null;
	// předpočítání derivací
	Vector[] deriv = derivace();
	// zkoumání polohy parametru pro určení krajních bodů
	for (int i=0; i<ctrlPoly.Count; i++)
		if (t <= ((RichPoint)ctrlPoly[i]).T)
		{
			// je-li parametrem 1. bod, beru tento a následný bod za krajní
			// jinak beru body o 1 zpětně abych přesně pokryl poslední krajní bod
			if(i==0)
			{
				b1 = (RichPoint)ctrlPoly[i];
				b2 = (RichPoint)ctrlPoly[i+1];
				db1 = deriv[i];
				db2 = deriv[i+1];
				break;
			}
			else
			{
				b1 =(RichPoint)ctrlPoly[i-1];
					b2 =(RichPoint)ctrlPoly[i];
				db1 = deriv[i-1];
				db2 = deriv[i];
				break;
			}
		}
	Vector f = hermitovaKubika (b1,b2,db1,db2,t);
	return f;
}

Metoda Vector[] derivace()[editovat | editovat zdroj]

Hledání derivací pro přirozený spline. Kód je spíše přípravou 1. matice pro řešení rovnic, pokračování tedy v metodě "reseniRovnice".

 [2 1       ] [D[0]]   [3(x[1] - x[0])  ]
 |1 4 1     | |D[1]|   |3(x[2] - x[0])  |
 |  1 4 1   | | .  | = |      .         |
 |    ..... | | .  |   |      .         |
 |     1 4 1| | .  |   |3(x[n] - x[n-2])|
 [       1 2] [D[n]]   [3(x[n] - x[n-1])]

Výsledkem je pole vektorů derivaci D[i] (X,Y,0).

private Vector[] derivace()
{
	// n = pocet bodu kontrolniho polygonu
	int n = ctrlPoly.Count;
	// definice prvku pro 3-diag. matici a matice s deriv bodu X a Y
	double[] ma = new double[n];
	double[] mb = new double[n];
	double[] mx = new double[n];
	double[] my = new double[n];
	double[] dx = new double[n];
	double[] dy = new double[n];
	// definice vypoctenych derivaci (vektor(dX,dY,0))
	Vector[] deriv = new Vector[n];
	for (int i=1; i<n-1; i++)
	{
		// vytvoreni prvku 3-diagonalni matice A
		ma[i] = 1;
		mb[i] = 4;
		// vytvoreni prvku sloupcove derivaci bodu pro X
		mx[i] = 3 * ( ((RichPoint)ctrlPoly[i+1]).Locate.X - ((RichPoint)ctrlPoly[i-1]).Locate.X );
		// vytvoreni prvku sloupcove derivaci bodu pro Y
		my[i] = 3 * ( ((RichPoint)ctrlPoly[i+1]).Locate.Y - ((RichPoint)ctrlPoly[i-1]).Locate.Y );
	}
	// doplneni kraju diag. matice A
	ma[0] = ma[n-1] = 1;
	mb[0] = mb[n-1] = 2;
	// doplneni kraju diag. matice s derivacemi bodu pro X
	mx[0] = 3 * ( ((RichPoint)ctrlPoly[1]).Locate.X - ((RichPoint)ctrlPoly[0]).Locate.X );
	mx[n-1] = 3 * ( ((RichPoint)ctrlPoly[n-1]).Locate.X - ((RichPoint)ctrlPoly[n-2]).Locate.X );
	// doplneni kraju diag. matice s derivacemi bodu pro Y
	my[0] = 3 * ( ((RichPoint)ctrlPoly[1]).Locate.Y - ((RichPoint)ctrlPoly[0]).Locate.Y );
	my[n-1] = 3 * ( ((RichPoint)ctrlPoly[n-1]).Locate.Y - ((RichPoint)ctrlPoly[n-2]).Locate.Y );
	// vypocet derivaci z rovnic 3-diag. matice
	dx = reseniRovnic (ma, mb, ma, mx);
	dy = reseniRovnic (ma, mb, ma, my);
	for (int i=0; i<n; i++)
	{
		deriv[i] = new Vector(dx[i],dy[i],0);
	}
	return deriv;
}

Metoda double[] reseniRovnic (...)[editovat | editovat zdroj]

Řešení rovnic z daných matic - viz metoda "derivace". Parametry:

a
1 tj. pravé okolí diagonály 3-diag. matice,
b
2 resp 4 tj. diagonála 3-diag. matice,
c
1 tj. pravé okolí diagonaly 3-diag. matice,
d
1-sloupcová matice z pravé strany rovnice.

Metoda vrací pole derivací pro daný směr X nebo Y.

private double[] reseniRovnic (double[] a, double[] b, double[] c, double[] d)
{
	// n = pocet bodu kontrolniho polygonu
	int n = ctrlPoly.Count;
	double[] deriv = new double[n];
	double[] e = new double[n];
	double[] f = new double[n];
	// převod matice na horní trojúhelníkovou
	e[0] = b[0];
	f[0] = d[0];
	for (int i=1; i<n; i++)
	{
		e[i] = b[i] - (c[i-1] * a[i] / e[i-1]);
		f[i] = d[i] - (f[i-1] * a[i] / e[i-1]);
	}
	// zpětná substituce a tedy výpočet rovnic
	deriv[n-1] = f[n-1] / e[n-1];
	for (int j=n-2; j>=0; j--)
	{
		deriv[j] = (f[j] - c[j] * deriv[j+1]) / e[j];
	}
	return deriv;
}

Metoda Vector hermitovaKubika(...)[editovat | editovat zdroj]

Spočítá a do vnitřní proměnné (pole richpointů) uloží derivace v bodech křivky. Derivace jsou počítány pomocí FMILL a v krajních bodech pomocí BESSEL. Body musejí mít UNIFORMNÍ PARAMETIZACI.

Vector hermitovaKubika(RichPoint P0, RichPoint P1, Vector P0C, Vector P1C, double t)
{
	/// POZOR
	/// při násobení zlomkem (3/2) je nutno UVÉST TYP (3f/2), jinak celočíselné dělení
	/// POZOR
	/// definice koeficientů
	Vector a,b,c,d;
	double h,u0;
	// kontrola parametrů bodů
	h = P1.T - P0.T;
	if (h==0)
		throw new Exception("Two or more points cannot have the same parameter");
	u0=P0.T;
	// koeficienty viz skripta - vzorec [J.DRDLA - Počítačová grafika, Olomouc 1991]
	a = P0.Locate;
	b = P0C;
	c = (3f/Math.Pow(h,2)) * (P1.Locate + ((-1)*P0.Locate)) +
		((-1f/h) * ((2*P0C) + P1C));
	d = (2f/Math.Pow(h,3)) * (P0.Locate + ((-1)*P1.Locate)) +
		(1f/Math.Pow(h,2)) * (P0C + P1C);
	// dosazení koeficientů do rovnice
	Vector vysledek;
	vysledek = a+(t-u0)*b + Math.Pow((t-u0),2)*c + Math.Pow((t-u0),3)*d;
	return vysledek;
}

Metoda Vector GetPoint(double t)[editovat | editovat zdroj]

Přetěžovaná metoda třídy Curve. V zadaném parametru vrátí polohový vektor bodu na křivce.

public override Vector GetPoint(double t)
{
	double w = ((RichPoint) ctrlPoly[0]).T;
	double w1 = ((RichPoint) ctrlPoly[0]).T;
	for (int q=0; q<ctrlPoly.Count;q++)
	{
	if (((RichPoint) ctrlPoly[q]).T>w) w=((RichPoint) ctrlPoly[q]).T;
	if (((RichPoint) ctrlPoly[q]).T<w1) w1=((RichPoint) ctrlPoly[q]).T;
	}
	Vector ret = new Vector();
	Vector v;
	double inter = w-w1;
	if (inter==0) throw new Exception("Sme vyrizeny. A vo tom ani tralala.");
	for (int i = 0;i<ctrlPoly.Count;i++)
	{
	v = ((RichPoint)ctrlPoly[i]).Locate;
	double nt = (t-w1)/inter;
	ret = ret + BernsteinPolynom(ctrlPoly.Count-1,i,nt)*v;
	}
	return ret;
}

Metoda RichPoint CreateRichPoint()[editovat | editovat zdroj]

Přetěžovaná metoda třídy Curve. Vytvoření řídícího bodu pro křivku. Oproti rodiči má navíc rozšíření o váhový koeficient.

public override RichPoint CreateRichPoint()
{
	double w = ((RichPoint) ctrlPoly[0]).T;
	for (int q=0; q<ctrlPoly.Count;q++)
	{
	if (((RichPoint) ctrlPoly[q]).T>w) w=((RichPoint) ctrlPoly[q]).T;
	}
	RichPoint tw = new RichPoint(new Vector(), w+0.1);
	return tw;
}

Autoři[editovat | editovat zdroj]

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.

Podívejte se také na[editovat | editovat zdroj]