[C++] Spline/Bezier izracun tock

Naredil sem svojo funkcijo za risanje Spline in Bezier krivulj. Stvar deluje vredu in hitro, vendar vem da ne optimalno. Zanima me, če bi kdo znal najti kakšno boljšo rešitev.

Trenutno je stanje takšno:
- izračunam koeficiente kubičnega polinoma za vsak segment krivulje (od ene do druge točke, za vse pare točk).
- ko dobim funkcijo (koeficiente) za vsak segment je treba nekako izračunat vse točke, ki jih je treba narisat
- pri Spline vzamem 10kratno razdaljo od začetne do končne točke, zato da zagotovo ne bi zmanjkalo pikslov (krivulja lahko vmes kar precej zavije), pri Bezier pa vsoto razdalj med štiri točkami segmenta (krivulja bo vedno znotraj konveksnega "ogrodja")
- delim 1/razdalja in iteriram parameter t kubične funkcije od 0 do 1 s takšnim korakom
- pri vsaki iteraciji izračunam točko in preverim, če je enaka prejšnji ter, če je odvečna (končna linija mora imeti povsod "debelino" 1 pixel, zato odstranim en piksel, če naslednji in prejšnji oba "zavijeta" pod kotom 90 stopinj)

Optimizacija bi bila sigurno možna v zadnjih dveh fazah opisanega algoritma. Rad bi nekako izvedel kakšen je največji korak pri iteraciji t-ja, da bom še vedno dobil vse potrebne piksle in ne bo noben vmes zmanjkal. Poleg tega mogoče obstaja še boljši način za zagotavljanje "čiste" krivulje brez odvečnih pikslov.

moja koda:

void Spline::UpdatePoints()
{
	allPoints.Clear();
	
	if (points.Size() <= 1) return;
	
	// Solve equation system (after Bartels) in order
	// to obtain function derives for every segment
	//------------------------------------------------	
	Matrix<double> leftSides(points.Size(), points.Size());
	Matrix<double> rightSidesX(points.Size(), 1);
	Matrix<double> rightSidesY(points.Size(), 1);
	ArrayList<double> derivesX;
	ArrayList<double> derivesY;
	
	// fill Bartels matrix
	for (int n=0; n<points.Size(); n++)
	{
		for (int m=0; m<points.Size(); m++)
			leftSides.SetData(n,m, 0);
			
		Point pt1, pt2;
		
		if ( n == 0 ){
		
			leftSides.SetData(0,0, 2);
			leftSides.SetData(0,1, 1);
			
			pt1 = points[0];
			pt2 = points[1];
						
		} else
		if ( n == points.Size()-1 ){
		
			leftSides.SetData(n,n-1, 1);
			leftSides.SetData(n,n,   2);
			
			pt1 = points[points.Size()-2];
			pt2 = points[points.Size()-1];
			
		} else {
		
			leftSides.SetData(n,n-1, 1);
			leftSides.SetData(n,n,   4);
			leftSides.SetData(n,n+1, 1);
			
			pt1 = points[n-1];
			pt2 = points[n+1];			
		}
		
		rightSidesX.SetData(n,0, 3*(pt2.x - pt1.x));
		rightSidesY.SetData(n,0, 3*(pt2.y - pt1.y));
	}
	
	// solve
	SolveEquations(&leftSides, &rightSidesX, &derivesX);
	SolveEquations(&leftSides, &rightSidesY, &derivesY);
	
	// Interpolate the other points for each segment
	//------------------------------------------------
	Point ptPrev, ptPrevPrev;
	
	for (int i=0; i<points.Size()-1; i++)
	{
		Point pt1=points[i], pt2=points[i+1];

		//Calculate cubic function coefficients
		double aX = pt1.x;
		double bX = derivesX[i];
		double cX = 3*(pt2.x - pt1.x) - 2*derivesX[i] - derivesX[i+1];
		double dX = 2*(pt1.x - pt2.x) + derivesX[i] + derivesX[i+1];

		double aY = points[i].y;
		double bY = derivesY[i];
		double cY = 3*(pt2.y - pt1.y) - 2*derivesY[i] - derivesY[i+1];
		double dY = 2*(pt1.y - pt2.y) + derivesY[i] + derivesY[i+1];
		
		//Calculate points
		double distance = sqrtf((pt2.x-pt1.x)*(pt2.x-pt1.x) + (pt2.y-pt1.y)*(pt2.y-pt1.y));
		double tStep = 1.0f/(distance*10);

		for (double t=0; t<=1.0f; t+=tStep)
		{
			Point pt;
			pt.x = aX + bX*t + cX*t*t + dX*t*t*t;
			pt.y = aY + bY*t + cY*t*t + dY*t*t*t;
			
			if (allPoints.Size() >= 1)
			{
				ptPrev = allPoints[allPoints.Size()-1];
				
				//We don't want to have multiplied points
				if (pt == ptPrev) continue;
			}
			if (allPoints.Size() >= 2)
			{
				ptPrevPrev = allPoints[allPoints.Size()-2];

				//We must neatly preserve line width of 1 pixel
				if ((pt.x == ptPrev.x) != (ptPrev.x == ptPrevPrev.x) &&
					(pt.y == ptPrev.y) != (ptPrev.y == ptPrevPrev.y))
					allPoints.RemoveAt(allPoints.Size()-1);	
			}

			allPoints.Add(pt);
		}
	}
}

void BezierLine::UpdatePoints()
{
	allPoints.Clear();

	if (points.Size() <= 1) return;
		
	// Interpolate other points for each segment
	//--------------------------------------------
	for (int p=0; p<points.Size()-1; p++)
	{
		Point pt1=points[p], pt2=points[p+1];
		pt1.Offset(points[p].rHandle.x, points[p].rHandle.y);
		pt2.Offset(points[p+1].lHandle.x, points[p+1].lHandle.y);
		Point pts[4] = {points[p], pt1, pt2, points[p+1]};

		//Calculate the cubic function coefficients
		float aX = 3.0f * (pts[1].x - pts[0].x);
		float bX = 3.0f * (pts[2].x - pts[1].x) - aX;
		float cX = pts[3].x - pts[0].x - aX - bX;
		
		float aY = 3.0f * (pts[1].y - pts[0].y);
		float bY = 3.0f * (pts[2].y - pts[1].y) - aY;
		float cY = pts[3].y - pts[0].y - aY - bY;
		
		//Calculate points
		double distance = PtDistance(pts[0], pts[1]) + PtDistance(pts[1], pts[2]) + PtDistance(pts[2], pts[3]);
		double tStep = 1.0f/distance;
		Point ptPrev, ptPrevPrev;
		
		for (double t=0; t<=1.0f; t+=tStep)
		{
			Point pt;
			pt.x = pts[0].x + aX*t + bX*t*t + cX*t*t*t + 0.5f;
			pt.y = pts[0].y + aY*t + bY*t*t + cY*t*t*t + 0.5f;
			
			if (allPoints.Size() >= 1)
			{
				ptPrev = allPoints[allPoints.Size()-1];
				
				//We don't want to have multiplied points
				if (pt == ptPrev) continue;
			}
			if (allPoints.Size() >= 2)
			{
				ptPrevPrev = allPoints[allPoints.Size()-2];

				//We must neatly preserve line width of 1 pixel
				if ((pt.x == ptPrev.x) != (ptPrev.x == ptPrevPrev.x) &&
					(pt.y == ptPrev.y) != (ptPrev.y == ptPrevPrev.y))
					allPoints.RemoveAt(allPoints.Size()-1);	
			}

			allPoints.Add(pt);
		}
	}
}

Samo ena ideja na hitro (prvo kar mi pade na pamet), ki ni nujno da deluje, sploh pa ne da je optimalno.
Od enega do drugega vertexa (se pravi za en segment) ne rišeš stalno z recimo 10x natančnostjo ampak to sproti prilagajaš. Če recimo dobiš ponovljeno točko potem povečaš step (step *= 1.5;). Če dobiš točko ki se ne dotika prejšnje potem ga zmanjšaš (step *= 0.666;) in izračunaš še vmesne točke. Če bi se takole še malo poigral s faktorjem množenja / deljenja stepa bi na lepih odsekih verjetno dobil kar stabilne vrednosti in bi ti to delalo hitro.

Druga stvar, ki bi jo jaz poskusil izboljšati je pa, da bi odpravil čim več preverjanj oz. jih čim bolj poenostavil. Tako da dobiš čimbolj straight-forward algoritem s čim manj ifi.

No zadnje pa v razmislek. Kadar delaš s polinomi so odvodi praktično zastonj. Odvod ti pokaže smer gibanja polinoma (in hitrost). Morda si s tem lahko kaj pomagaš.

rkb2:
tnx, that's the solution!

Gundolf:
Hvala za idejo, ampak mislim da je rešitev ki jo ponuja link rkb2-ja še boljša. Tole z if-i maš pa zelo prav!

Zadnja ideja z odvodi bi lahko bila tudi precej uporabna, še posebej zato ker pri Spline imaš odvode itak že izračunane, ker so del pogojev - sistema enačb (to je to kar se računa z Bartels-ovo matriko), pri Bezier-ju pa bi lahko vzel kar linije od kontrolnih točko do leve in desne točke segmenta (to kar je v CorelDRAW-u predstavljeno kot neke ročice), če se ne motim...

btw: še slikica kako moje krivulje zgledajo v ShmanssyGUI-ju (moj portable GUI):

Heh, pa res nisem pomislil na te probleme pri bisekciji...No vsaj sortiranja na koncu js ne bi imel, ker bi sproti uredil stvari tako da se ustavljajo tam kjer je treba. Če uporabim LinkedList za shranjevanje točko ima vstavljanje komplesnost O(1) in tle ni problema. Če bi še malo premislil bi mogoče našel rešitev tudi za problem sekanja odseka krivulje s samim sabo.

Z vzemanjem večjega korak pri iteraciji in risanja crt za povezovat točke sem pa že sprobal in ne pride lepo. Zelo težko je potem kontrolirat točno kako se bosta črti stikali in včasih pritejo grde napakice, npr. kakšen piksel nekje štrli vn in podobno....vglavnem, je hitra rešitev ampak da vizualno grd rezultat.

Uspelo mi je nardit stvar z metodo bisekcije in to brez potrebe po končnem sortiranju. Za sprotno shranjevanje točk sem uporabil LinkedList in kazalce na njegove "iteratorje" (nosilce elementov), ki se pri vstaljanju ohranjajo (če le dobro paziš kako stvar postaviš, vrstni red kode itd.).

Prvič sem vseskupaj implementiral v celoti z iteracijo, ampak je stvar delala še dvakrat počasneje kot prej, zaradi shranjevanja notranjih parametrov iteracije na sklad itd. Potem sem stvar spremenil da deluje z rekurzivno metodo. Tako sem prišel do rezultata, ki je pomoje še najhitrejši kot se da za to metodo. Nažalost pa tudi rekurzija pobere kar nekaj časa, tako da sem pri časovnih meritvah z bisekcijo dobil le za malo malo počasnješi algoritem od prejšnjega, ko sem računal 10krat več kot je treba (seveda pa govorimo o tisoč osvežitvah krivulje v pol sekunde, tako da vseeno ni nekega velikega problema ;).

Poleg tega da časovno nisem nič pridobil, so se pokazale tudi slabosti metode bisekcije. Ko sem krivuljo zvijal do prav vizualno "bolečih" kotov, so včasih deli krivulje izginili, kar pa se mi s prejšnjim algoritmom nikoli ne dogaja.

Poanta zgodbe: gremo rajši delat na še dodatnem izboljšanju prejšnjega algoritma, kar pomoje lahko privede do kick-ass rezultatov (10000 osvežitve na sekundo bi blo carsko).

Tukaj je še koda za bisekcijo, če vas zanima kako se to izpelje BREZ SORTIRANJA:

double aX, bX, cX, dX;
double aY, bY, cY, dY;

void Spline::CalculatePoint(double t1, double t2, Point *p1, Point *p2, Iterator<Point> *insertPlace)
{	
	static Point p;
	double t;
			
	t = (t1 + t2)/2.0f;
	p.x = aX + bX*t + cX*t*t + dX*t*t*t;	
	p.y = aY + bY*t + cY*t*t + dY*t*t*t;
	
	if (p == *p1 || p == *p2)
		return;
	
	if ((p2->x == p.x) != (p.x == p1->x) &&
		(p2->y == p.y) != (p.y == p1->y))
		return;
	
	newPoints.InsertAt(insertPlace, p);
	
	CalculatePoint(t, t2, &newPoints.ElementAt(insertPlace), p2, insertPlace->next);
	CalculatePoint(t1, t, p1, &newPoints.ElementAt(insertPlace), insertPlace);
}

void Spline::UpdatePoints()
{
	allPoints.Clear();
	newPoints.Clear();
	
	if (points.Size() <= 1) return;
	
	// Solve equation system (after Bartels) in order
	// to obtain function derives for every segment
	//------------------------------------------------	
	Matrix<double> leftSides(points.Size(), points.Size());
	Matrix<double> rightSidesX(points.Size(), 1);
	Matrix<double> rightSidesY(points.Size(), 1);
	ArrayList<double> derivesX;
	ArrayList<double> derivesY;
	
	// fill Bartels matrix
	for (int n=0; n<points.Size(); n++)
	{
		for (int m=0; m<points.Size(); m++)
			leftSides.SetData(n,m, 0);
			
		Point pt1, pt2;
		
		if ( n == 0 ){
		
			leftSides.SetData(0,0, 2);
			leftSides.SetData(0,1, 1);
			
			pt1 = points[0];
			pt2 = points[1];
						
		} else
		if ( n == points.Size()-1 ){
		
			leftSides.SetData(n,n-1, 1);
			leftSides.SetData(n,n,   2);
			
			pt1 = points[points.Size()-2];
			pt2 = points[points.Size()-1];
			
		} else {
		
			leftSides.SetData(n,n-1, 1);
			leftSides.SetData(n,n,   4);
			leftSides.SetData(n,n+1, 1);
			
			pt1 = points[n-1];
			pt2 = points[n+1];			
		}
		
		rightSidesX.SetData(n,0, 3*(pt2.x - pt1.x));
		rightSidesY.SetData(n,0, 3*(pt2.y - pt1.y));
	}
	
	// solve
	SolveEquations(&leftSides, &rightSidesX, &derivesX);
	SolveEquations(&leftSides, &rightSidesY, &derivesY);
	
	// Interpolate the other points for each segment
	//------------------------------------------------
	Point ptPrev, ptPrevPrev;
	LineSegment *whole, *newSeg1, *newSeg2;
	SmartLinkedList<LineSegment*> segments;
	calculations = 0;
	
	for (int i=0; i<points.Size()-1; i++)
	{
		Point pt1=points[i], pt2=points[i+1];

		//Calculate cubic function coefficients
		aX = pt1.x;
		bX = derivesX[i];
		cX = 3*(pt2.x - pt1.x) - 2*derivesX[i] - derivesX[i+1];
		dX = 2*(pt1.x - pt2.x) + derivesX[i] + derivesX[i+1];

		aY = points[i].y;
		bY = derivesY[i];
		cY = 3*(pt2.y - pt1.y) - 2*derivesY[i] - derivesY[i+1];
		dY = 2*(pt1.y - pt2.y) + derivesY[i] + derivesY[i+1];
		
		//Calculate points

		newPoints.Add(pt1);
		newPoints.Add(pt2);
		
		CalculatePoint(0.0f,
				     1.0f,
				     &newPoints.ElementAt(newPoints.First()),
			              &newPoints.ElementAt(newPoints.First()->next),
			              newPoints.First()->next);
	}
}

Vedno lahko izračunaš še eno piko med dvema že enakima, in če je še ta enaka - nehaš. Samo imam občutek da to ne bo blagodejno vplivalo na hitrost

Drugače pa ne, pri bisekciji sem brez idej kako rešiti tole z dt. Zato mi pa ni všeč bisekcija . Lahko pa se poigraš z rekurzijo in ji najprej preč pomečeš čim več parametrov. Potem pa pretvoriš v iteracijo in tisto iteracijo še malo optimiraš. Če boš naredil prav ni hudič, da ne bi dosegel vsaj 1.3x-ne pohitritve.

Aja, hehe, za tole "p.x = aX + bX*t + cX*t*t + dX*t*t*t;" obstaja način izračuna s 4 seštevanji in le 3 množenji. Tule v "t = (t1 + t2)/2.0f;" pa bi lahko deljenje zamenjal z množenjem.

BTW, zakaj pa rabiš tako hitro osvežitev ?

"ko pa gre pod koto 45 stopinj..." - verjetno s tem ne misliš, da bi dejansko računal kot, ampak samo preverjal, če se spremeni le x ali le y, ali če se spremenita oba, je tako? U bistvu bi ti rad, da se korak tStep ne spreminja z fiksno vrednostjo, ampak da bi točno izračunal za kolko ga moraš spremenit, da boš ČIMHITREJ dosegel pravo vrednost tStep. Hmm...že poskušam nardit ;)

BigMan: Zamisli si, recimo, da delaš svoj GUI toolkit in moraš nardit čimhitrejše risalne funkcije. Zdaj bi pa rad imel Spline kot npr. osnovo za nek shape, s katerim boš klippal eno kvadratno teksturo. Recimo da v končni fazi to predstavlja jezerce z neko teksturo "vode" nafilano not in zdaj bi rad da skozi čas to jezerce malo "miga" na robovih, torej je animirano, in vse skupaj izrisuješ v OpenGLju in hočeš met čimvečji FPS. In recimo da imaš na ekranu naenkrat 3 taka jezerca. In recimo da imaš zraven na ekranu še animiran vektorski GUI, ki tudi uporablja Spline za svoje oblike in ravno v istem trenutku hočeš da gumb na nek fancy način popne na sliko (se poveča iz nule na neko velikost npr.)

Ja, sem že prej dojel No, zdaj sem ugotovil da tole z točnim računanje dt-ja ima tudi svoje probleme. Kot prvo, treba je uporablja PointF točke, ker moraš imet float vrednost od prejšnjih točk ko računaš. To pobere dodaten čas pri ostail if stavkih (ko npr. izločaš "odvečne grde točke") ker moraš cast-at nazaj v int. Potem so tu še dodatni izračuni: 2*abs, 1*max, potem pa še dt moraš množit oz. delit, namesto ga prištevat odštevat.

Vse to vzame tolko časa, da že ko sem nardil samo to, da se je dt povečal, če je premajhen (tStep *= (1.0f/diff)), še brez čekiranja če je prevelik, torej z izpuščanjem nekaterih točk, je algoritem že vzel enako časa kot prejšnji.

Potem ko sem dodal še preverjanje in vračanje nazaj, če je dt previlk, je bil že priblžno
100 milisek/10000 osvežitev počasnejši. Poleg tega se algoritem na mestih, kjer krivulja ostro zavije nazaj včasih malo dlje vstavi, ker se dt lovi okrog ene vrednosti, malo več, malo manj, malo več, malo manj, ko so točke ful skupaj natlačene. To se je potem opazlo da ko sem z miško na hitro premikal točke krivulje okrog, gibanje ni zgledalo več tako gladko ampak včasih malo "zobato".

Ne, še prej je treba zafilat skljenjene krivulje Ampak to mislim da niti ne bo tak problem.... Filood fill bi bla rešitev ampak pomoje preve počasna. Boljše bo, če iz točk "sparsam" horizntalna zaporedja notranjih točk in jih zrišem s horizontalnimi črtami. Piece of cake