Geometrie/C2 kubická interpolace

Popis

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í

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

$P_{i}(u_{i})=P_{i-1}(u_{i})$

$P'_{i}(u_{i})=P'_{i-i}(u_{i})$

$P''_{i}(u_{i})=P''_{i-1}(u_{i}),i=1,\dots ,n-1$

Volbou čísel u₀,...,u_n se určuje parametrizace spline křivky. Pro parametrizaci obvykle volíme u_i=i (Uniformní parametrizace).

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

$P_{i}(u)=a_{i}+b_{i}(u-u_{i})+c_{i}(u-u_{i})^{2}+d_{i}(u-u_{i})^{3},i=0,...,n-1$

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 a_i , b_i , c_i a d_i. 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í:

Vektory prvních derivací v krajních bodech P´₀ a P´_n.
Vektory druhých derivací v krajních bodech P´´₀ a P´´_n. Pokud použijete P´´₀=0 a P´´_n=0 hovoříme o přirozené spline křivce.
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 P´₀,...,P´_n. Jednotlivé kubiky P_i(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 a P´´(u _n)=0 dostaneme soustavu (n+1) rovnic pro výpočet tečných vektorů P´₀,...,P´_n. Zvolíme-li uniformní parametrizaci, pak platí:

${\begin{bmatrix}2&1&0&&&\ldots &0\\1&4&1&0&&\ldots &0\\0&1&4&1&0&\ldots &0\\\vdots &&&&&&\vdots \\0&\ldots &&0&1&4&1\\0&\ldots &&&0&1&2\end{bmatrix}}\cdot {\begin{bmatrix}P'_{0}\\P'_{1}\\\vdots \\P'_{n}\end{bmatrix}}={\begin{bmatrix}3\cdot (P_{1}-P_{0})\\3\cdot (P_{2}-P_{0})\\\vdots \\3\cdot (P_{n}-P_{n-2})\\3\cdot (P_{n}-P_{n-1})\end{bmatrix}}$

Vlastnosti

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

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)

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()

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 (...)

Ř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(...)

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)

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()

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

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