Geometrie/C1 kubická interpolace

Popis

C1 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ášť. Každý úsek C1 interpolační křivky bude kubický polynom spočítaný pomocí kubické Hermitovy interpolace. Pro C1 interpolační křivku musí být zabezpečena C1 spojitost. Tzn. v každém opěrném bodě mají sousední křivky stejný tečný vektor. K tomu potřebujeme tečné vektory v opěrných bodech vypočítat.

Vyjádření

Pro výpočet tečných vektorů opěrných bodů jsme použili metodu BESSEL. Vytvoříme interpolační parabolu $G(u)=au^{2}+bu+c$ určenou posloupností tří po sobě jdoucích opěrných bodů ( $P_{i-1},P_{i},P_{i+1}$ ). Tečný vektor prostředního bodu $P_{i}$ počítáme vždy jako derivaci interpolační paraboly tomto bodě. Zvolíme parametrizaci(obvykle $G(0)=P_{i-1},G(1)=P_{i},G(2)=P_{i+1}$ ) opěrných bodů paraboly, za levé strany těchto tří rovnic dosadíme konkrétní tvary (tj. po dosazení za parametr u) pravých stran paraboly a vypočítáme koeficienty a, b, c ( $a={1 \over 2}(P_{i-1}-2P_{i}+P_{i+1}),b={1 \over 2}(-3P_{i-1}+4P_{i}-P_{i+1}),c=P_{i-1}$ ). Tyto koeficienty dosadíme do derivace paraboly $G'(u)=2au+b$ . Tečný vektor je tedy roven $P'_{i}=G'(u)=2au+b$ , kde $u=1$ . V případě prvního (resp. posledního) opěrného bodu vytvoříme parabolu z prvních (resp. posledních) tří opěrných bodů a parametr $u$ volíme 0 (resp. 2).

Vlastnosti

Pro výpočet tečných vektorů by měla být splněna podmínka uniformní parametrizace.

Algoritmizace

Metoda HermitovaKubika pro zadaný oblouk křivky (P0 a P1 a jejich tečné vektory P0C a P1C) a parametr vrátí polohový vektor bodu na křivce. Polohový vektor je spočítán pomocí Hermitových kubik. Oblouk je tvořen dvěma po sobě jdoucími body a parametr musí spadat mezi ně.

Vector: třída popisující vektor,
Richpoint: třída popisující bod,
ctrlPoly: body řídícího polygonu,
P0: počáteční bod oblouku,
P1: koncový bod oblouku,
P0C: derivace v počátečním bodu oblouku,
P1C: derivace v koncovém bodu oblouku,
t: parametr (musí náležet do daného oblouku).

Nejprve si popišme pomocnou funkci VratDerivace, jež vrací pole derivací v opěrných bodech.

RichPoint[] VratDerivace(ArrayList Body) {
 RichPoint[] Derivace=new RichPoint[Body.Count];
 int n=Body.Count;
 int u, j;
 Vector a,b;
 for (int i=0; i<=n-1; i++) {
    // volba j - pro bod na 1. oblouku nutvo volit parabolu prvních
    // 3 bodů a pro poslední oblouk nutno volit parabolu poslední tří bodů
    if(i==0) {
        u=0; // zvolená parametrizace pro první bod paraboly
        j=i+1; // pro první bod počítám derivaci v prvním bodě paraboly
    } else if(i==n-1) {
        u=2; // zvolená parametrizace pro poslední bod paraboly
        j=i-1; // pro poslední bod počítám derivaci v posledním bodě dané paraboly
    } else {
        u=1; // zvolená parametrizace pro vrchol paraboly
        j=i; // jinak
    }
    a = 1f/2*( ((RichPoint)Body[j-1]).Locate -
        2*((RichPoint)Body[j]).Locate + ((RichPoint)Body[j+1]).Locate );
    b = 1f/2*( -3*((RichPoint)Body[j-1]).Locate +
        4*((RichPoint)Body[j]).Locate + (-1)*((RichPoint)Body[j+1]).Locate );
    Derivace[i]=new RichPoint(2*a*u+b);
    }
 return Derivace;
}

A nyní již k vlastnímu algoritmu.

Vector HermitovaKubika(RichPoint P0, RichPoint P1, RichPoint P0C,
  RichPoint P1C, double t) {
 Vector a,b,c,d; //viz skripta - vzorec [J.DRDLA - Počítačová grafika, Olomouc 1991]
 double h,u0;
 h=P1.T-P0.T;
 if(h==0)
  throw new Exception("Two or more points cannot have the same parameter");
 u0 = P0.T;
 a = P0.Locate;
 b = P0C.Locate;
 c = (3f/Math.Pow(h,2))*(P1.Locate+((-1)*P0.Locate)) 
     +(-1f/h)*((2*P0C.Locate)+P1C.Locate));
 d = (2f/Math.Pow(h,3))*(P0.Locate+((-1)*P1.Locate)) +
     (1f/Math.Pow(h,2))*(P0C.Locate+P1C.Locate);
 Vector vVysledek;
 vVysledek=a+(t-u0)*b + Math.Pow((t-u0),2)*c + Math.Pow((t-u0),3)*d;
 return vVysledek;
}

Metoda Hermit pro daný parametr zvolí oblouk křivky a vrátí polohový vektor bodu pro zadaný parametr.

Vector Hermit(double t) {
 this.prpDerivace=VratDerivace(this.ctrlPoly);
 // slouží pro výběr oblouku křivky
 RichPoint PredchoziBod, DalsiBod, DerivacePredchoziho 
 DerivaceDalsiho;
 PredchoziBod=DalsiBod=DerivacePredchoziho=DerivaceDalsiho=null;
 // vyberu oblouk křivky pro daný parametr
 for (int i=0; i<ctrlPoly.Count; i++)
 if (t <= ((RichPoint)ctrlPoly[i]).T) {
   ...
 }
  ...
 return
 HermitovaKubika(PredchoziBod, DalsiBod, 
   DerivacePredchoziho, DerivaceDalsiho, t);
}

Pomocná metoda.

public override Vector GetPoint(double t) {
 return Hermit(t);
}

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.