Create a new version of paste:

%---------------GLAVNI PROGRAM------------------
function bezierjeviZlepki
	%kont=[1 2 4 6;0 4 6 1]
	%tic;
	%narisiBezierCastel(kont,500000,1,1);
	%toc
	%figure(2)
	%tic;
	%narisiBezierCastel(kont,500000,1,1);
	%toc
	%return
	
	zi=linspace(-5,5,70);
	zacetneTocke=[zi;zi.^3]

	%S crno barvo narisemo vse zacetne tocke, katerim se moramo prilagajati
	plot(zacetneTocke(1,:),zacetneTocke(2,:),['.','k'])
	hold on
	
	%Osi grafa dolocimo tako, da bo na grafu  vidna celotna funkcija
	axis([min(zacetneTocke(1,:))-3 max(zacetneTocke(1,:))+3 min(zacetneTocke(2,:))-3 max(zacetneTocke(2,:))+3])

	
	d=length(zacetneTocke);	%stevilo zacetnih tock
	raz=Inf;	%povprecna razdalja krivulje do tock
	stKrivulj=7;	%najvecje stevilo krivulj, ki jih dopustimo da jih algoritem uporabi
	toleranca=0.15;	%maksimalna povprecna razdalja ki jo se dopustimo
	stTockKrivulja=100;	%Stevilo tock ene bezierjeve krivulje
	kontr=startControlPoints(zacetneTocke,[]);	%zacetneTocke
	prejKont=[];	%kontrolne tocke iz prejsnjega koraka
	vseKontr=[];	%vse kontrolne tocke, ki definirajo koncno krivuljo
	prejTocke=[];	%zacetne tocke, ki "pripadajo" dani krivulji
	X=[];	%tocke krivulje
	steviloKrivulj=0;	%koliko krivulj smo uporabili, samo za izpis
	o=1;
	while raz>toleranca
		b=floor(length(zacetneTocke)/o);	%Razdelimo na o krivulj, o tock
		vseKontr=[];
		for g=0:o-1
			if(g*b+b<=length(zacetneTocke) && g!=o-1)	%Tocke razdelimo na o delov, o krivulj
				tocke=zacetneTocke(:,g*b+1:g*b+b);
			else
				tocke=zacetneTocke(:,g*b+1:length(zacetneTocke));
			end
			if( g>0 )
				kontr=startControlPoints(tocke,prejKont);	%dolocimo zacetne tocke dane krivulje
				spr=[0 0 1 0];
				
				[kontr,d]=prilagodiKontrolneTocke(tocke,kontr,spr,0.2);	%in jih prilagodimo
			else
				kontr=startControlPoints(tocke,[]);
				[kontr,d]=prilagodiKontrolneTocke(tocke,kontr,[0 1 1 0],0.2);
			end
			if(length(prejKont)>0)	%ce imamo ze vec kot eno krivuljo, potem lahko vmesne kontrolne tocke prilagajamo(tiste na stičiščih dveh krivulj)
				[lK,dK]=prilagodiPovezovalneKontrolneTocke([prejTocke,tocke],prejKont,kontr,0.2);	%dobimo nove prejsnje in nove zdajsnje kontr. tocke
				vseKontr=[vseKontr(:,1:length(vseKontr)-4),lK,dK];
				X=[X(:,1:length(X)-stTockKrivulja),narisiBezierCastel(lK,stTockKrivulja,0,0),narisiBezierCastel(dK,stTockKrivulja,0,0)];	%tocke krivulje
			else	%v nasprotnem primeru imamo samo eno krivuljo
				vseKontr=[vseKontr,kontr];
				X=[X,narisiBezierCastel(kontr,stTockKrivulja,0,0)];	%tocke krivulje
			end
			prejTocke=tocke;
			prejKont=kontr;
		end
		dist=povpRazdalja(tocke,X);
		if(dist<raz)
			raz=dist
			kontrolneKrivulje=vseKontr;
			steviloKrivulj=o;
		end
		prejTocke=[];
		prejKont=[];
		X=[];
		if(o>=stKrivulj)
			break
		end
		o++;	%Povecamo stevilo krivulj
	end
	raz
	steviloKrivulj
	for i=1:4:length(kontrolneKrivulje)
		newKont=kontrolneKrivulje(:,i:i+3);
		X=narisiBezierCastel(newKont,stTockKrivulja,0,1);
	end

endfunction
%------------------FUNKCIJE----------------------------
%Funkcija vrne zacetne kontrolne tocke krivulje, tako da zracuna povprecje danih tock, ce se ni prejsnje krivulje,
%ce pa prejsnja krivulja ze obstaja, pa moramo paziti na zveznosti in odvedljivost
%Podamo ji 2 argumenta:
%tocke=tocke katerim se prilagajamo
%prejKontr=kontrolne tocke prejsnje krivulje, ki je lahko prazen vektor, kar pomeni, da prejsnje krivulje ni
function kontr=startControlPoints(tocke,prejKontr)
	d=length(tocke); pol=d-floor(d/2);
	avg1=[sum(tocke(1,1:floor(d/2)))/floor(d/2) ; sum(tocke(2,1:floor(d/2)))/floor(d/2)];
	avg2=[sum(tocke(1,floor(d/2)+1:d))/pol ; sum(tocke(2,floor(d/2)+1:d))/pol];
	if(length(prejKontr)==4)	%Imamo prejsnje kontrolne tocke, moramo dolociti odvedljiv zlepek
		%V tem delu funkcija precejkrat napise err, samo zaradi preglednosti, kaj je narobe, te izpise ustvarja funkcija naslednjiKotntrolniTocki
		%Program kljub tem izpisom deluje pravilno
		vmes=naslednjiKontrolniTocki(prejKontr,avg1(1,1),Inf);
		if(length(vmes)==0)	
			vmes=naslednjiKontrolniTocki(prejKontr,Inf,avg1(2,1));
			if(length(vmes)==0)	
				vmes=naslednjiKontrolniTocki(prejKontr,Inf,Inf);	
			end
		end
		kontr=[vmes,tocke(:,length(tocke)-1),tocke(:,length(tocke))];
	else
		kontr=[tocke(:,1),avg1,avg2,tocke(:,d)];
	end
endfunction

%Funkcija, ki prilagodi kontrolne tocke na sticiscu dveh krivulj, kjer moramo zagotavljati zveznost in odvodljivst krivulje
%Dobi 3 parametre:
%tocke=tocke, katerim se na tem odseku prilagajmo
%leveK=kontrolne tocke leve krivulje
%desneK=kontrolne tocke desne krivulje
%korak=korak s katerim popravljamo kontrolne tocke
function [kL,kD]=prilagodiPovezovalneKontrolneTocke(tocke,leveK,desneK,korak)
	%Prilagajamo dve kontrolni tocki
	%po obeh koordinatah
	raz=Inf;
	kL=leveK;
	kd=desneK;
	for k=1:2
		razlika=zeros(2,4);
		razlika(k,3)=korak;
		while(1==1)		%Vsaki tocki poiskusamo pristeti
			kontrD=fliplr(desneK)+razlika;
			kontrL=naslednjiKontrolniTocki(kontrD,leveK(1,3),Inf);	%podamo ji x, da poiscemo novo tocko
			if(length(kontrL)==0)
				kontrL=naslednjiKontrolniTocki(kontrD,Inf,leveK(2,3));	%podamo ji x, da poiscemo novo tocko
				if(length(kontrL)==0)
					kontrL=naslednjiKontrolniTocki(kontrD,Inf,Inf);
				end 
			end
			kontrL=[leveK(:,1:2),fliplr(kontrL)];
			bL=narisiBezierCastel(kontrL,100,0,0);
			bD=narisiBezierCastel(kontrD,100,0,0);
			X=[bL,bD];
			k=povpRazdalja(tocke,X);
			if(k<raz)
				raz=k;
				kL=kontrL;
				kD=kontrD;
			else 
				break;
			end	
		end
		while(1==1)	%oziroma odsteti
			kontrD=fliplr(desneK)-razlika;
			kontrL=naslednjiKontrolniTocki(kontrD,leveK(1,3),Inf);	%podamo ji x, da poiscemo novo tocko
			if(length(kontrL)==0)
				kontrL=naslednjiKontrolniTocki(kontrD,Inf,leveK(2,3));	%podamo ji x, da poiscemo novo tocko
				if(length(kontrL)==0)
					kontrL=naslednjiKontrolniTocki(kontrD,Inf,Inf);
				end 
			end
			kontrL=[leveK(:,1:2),fliplr(kontrL)];
			bL=narisiBezierCastel(kontrL,100,0,0);
			bD=narisiBezierCastel(kontrD,100,0,0);
			X=[bL,bD];
			k=povpRazdalja(tocke,X);
			if(k<raz)
				raz=k;
				kL=kontrL;
				kD=kontrD;
			else 
				break;
			end	
		end
	end
endfunction

%Funkcija, ki spreminja kontrolne tocke tam kjer nam ni treba skrbeti za odvodljivost krivulje
%tocke=zunanje tocke, katerim se prilagajamo
%kontrole=neke dane kotrolne tocke, ki jih bomo popravili
%prilagodi=vektor kjer so nenicelne vrednosti na indeksih tistih tock, ki jih lahko prilagajamo.
%recimo: Če imamo samo 1 krivuljo, lahko prilagajamo (ponavadi) samo 2. in 3. kontrolno tocko
function [newKont,raz]=prilagodiKontrolneTocke(tocke,kontrole,prilagodi,korak)
	pT=find(prilagodi);	%Dobimo ven katere tocke (na indeksih bomo prilagajali)
	raz=povpRazdalja(tocke,narisiBezierCastel(kontrole,200,0,0));
	newKont=kontrole;
	for i=1:10
		for i=1:length(pT)
			for p=1:2	%koordinata x in y
				razlika=zeros(2,4);
				razlika(p,pT(i))=korak;
				while(1==1)		%Vsaki tocki poiskusamo pristeti
					kontrole=newKont+razlika;
					k=povpRazdalja(tocke,narisiBezierCastel(kontrole,200,0,0));
					if(k<raz)
						raz=k;
						newKont=kontrole;
					else 
						break;
					end	
				end
				while(1==1)	%oziroma odsteti
					kontrole=newKont-razlika;
					k=povpRazdalja(tocke,narisiBezierCastel(kontrole,200,0,0));
					if(k<raz)
						raz=k;
						newKont=kontrole;
					else 
						break;
					end	
				end
			end
		end
		korak=korak/2;
	end
endfunction
%FUNKCIJA JE Z RAZVOJEM PROGRAMA "PROPADLA" IN SE V SAMEM PROGRAMU NE UPORABLJA VEC:
%Pridobivanje kontrolnih točk s pomočjo "reverse engineeringa" -> imamo 4 točke neke krivulje
%Če podamo funkciji 4 točke na krivulji, nam rekonstruira kontrolne točke, s tem da sta prva in zadnja začetna oz. končna
%Drugi argument so verjetnosti. Če definiramo 4 točke iz neke krivulje, še vedno lahko skozi te 4 točke potegnemo 
%neskončno mnogo krivulj, verjetnosti pa nam povej položaj na krivulji, seštevek danih verjetnosti mora biti 1
%Zadnji argument je nacin, ce je nacin enak 1 potem podamo 4 tocke na krivulji,
%ce pa je nacin enak 2 potem podamo prvi dve kontrolni tocki in 2 tocki na krivulji (uporabno za nadaljnje krivulje, za zlepke)
function controls = generirajKontrolneTocke(tocke,verjetnost,izris,nacin)
	%imamo tocke t0, t1, t2 in t3 v vektorju tocke
	%POMOC: http://polymathprogrammer.com/2007/06/27/reverse-engineering-narisiBezierCastel-curves/
	
	t0=tocke(:,1); t1=tocke(:,2); t2=tocke(:,3); t3=tocke(:,4);
	k0=t0, k3=t3 %Začnemo in končamo krivuljo v začetni oz. končni točki
	if(izris>0)	%Izrisemo tocke skozi katere mora krivulja
		plot(tocke(1,:),tocke(2,:),['o','r']);	%Narisemo tocke
	end

	u=verjetnost(1); v=verjetnost(2);
	%Kontrolni tocki k0=t0 in k3=t3 imamo ze, rabimo se k1 in k2
	%Uporabimo de Castelauovo formulo
	%t1=(1-u)^3*k0 + 3*(1-u)^2*u*k1 + 3*(1-u)*u^2*k2 + u^3*k3
	%Ce damo k1 in k2 na eno stran, kar iščemo, dobimo sistem 2 neznank k1 in k2:
	%3*(1-u)^2*u*k1 + 3*(1-u)*u^2*k2 = t1 - (1-u)^3*k0 - u^3*k3
	%3*(1-v)^2*v*k1 + 3*(1-v)*v^2*k2 = t2 - (1-v)^3*k0 - v^3*k3
	%Sistem lahko zapisemo v matricni obliki, leva stran:
	if(nacin<2)
		c=t1-(1-u)^3*k0-u^3*k3;
		d=t2-(1-v)^3*k0-v^3*k3;
		L=[3*(1-u)^2*u, 3*(1-u)*u^2; 3*(1-v)^2*v, 3*(1-v)*v^2]'
		D=[c'; d'];
		per=L\D;
		controls=[t0,per',t3];
	%V drugem nacinu imamo prakticno samo eno neznano kontrolno tocko, prvi dve dobimo iz funkcije naslednjiKontrolniTocki, 
	%da je zvezno odvedljiva krivulja, zadnja tocka na krivulji pa je zadnja kontrolna tocka, iscemo torej k2
	elseif(nacin>1)
		k1=t1;
		k2=(t2-((1-v)^3*k0 + 3*(1-v)^2*v*k1 + v^3*k3))/(3*(1-v)*v^2);
		controls=[k0,k1,k2,k3];
	endif
endfunction

%Funkcija izracuna povprecno razdaljo med danimi tockami in bezierjevo krivuljo
%Sprejme 2 argumenta:
%tocke = zacetne tocke, okoli katerih moramo zgraditi bezierjevo krivuljo
%X = tocke, ki predstavljajo bezierjevo krivuljo
%Obe matriki sta sestavljeni iz 2 vrstic, zgornja so vse x koordinate, spodnja vse y (X=[x;y])
function avgDist = povpRazdalja(tocke,X)
	stTock=length(tocke(1,:));	%Stevilo danih tock
	avgDist=0;
	for k=1:stTock
		tocka=tocke(:,k)';
		[dist,i]=min(sqrt((X(1,:)-tocka(1)).^2 + (X(2,:)-tocka(2)).^2));	%Racunanje minimalne razdalje (evklidove)
		avgDist+=dist/stTock;	%Utezena vsota
	end
endfunction

%Funkcija izracuna Bezierjevo krivuljo, ne nujno da 3. reda in jo vrne v obliki N tock
%DEFINIRANA Z BERNSTEINOVIM POLINOMOM
%Podamo ji 4 parametre:
%kTocke = kontrolne tocke: 1. vrstica x koordinate, 2. y koordinate
%N = "resolucija", s koliko tockami bomo predstavili krivuljo
%poligon = Ce je enak 1, potem narisemo poligon, drugace samo krivuljo
function X = narisiBezier(kTocke,N,poligon,krivulj)
	
	n=length(kTocke(1,:));    %Izracunamo stevilo kontrolnih tock
	a=linspace(0,1,N);	 %Izracunamo x koordinate tock na krivulji  
	I=(0:(n-1))'*ones(1,N);    %toliko kot je tock tolikokrat indeks i postavimo od 0-stopnje,
								%ponavadi 0-3, da ne bomo potrebovali for zanke
								
	P=ones(n,1)*a;           %Matrika verjetnosti za binomsko porazdelitev
	B=binopdf(I,n-1,P);  %Izracunamo B kar nad matriko, binopdf vrne binomsko porazdelitev
	X=kTocke*B;                       %Točke na krivulji, računane po 
	X(:,1)=[kTocke(1,1),kTocke(2,1)]';
	X(:,N)=[kTocke(1,n),kTocke(2,n)]';
	
	if(krivulj>0)
		plot(X(1,:), X(2,:));	%Narisemo se krivuljo
		hold on;
	end
	if(poligon==1)
		plot(kTocke(1,:), kTocke(2,:), ['-o','g']); %Narisemo poligon
	end
endfunction


%Funkcija izracuna Bezierjevo krivuljo, ne nujno da 3. reda in jo vrne v obliki N tock
%DEFINIRANA Z DE CASTELJAUOVO METODO  -> Po merjenjih je vec kot 10x hitrejsa od funkcije z Bernsteinovim polinomom
%Podamo ji 4 parametre:
%kTocke = kontrolne tocke: 1. vrstica x koordinate, 2. y koordinate
%N = "resolucija", s koliko tockami bomo predstavili krivuljo
%poligon = Ce je enak 1, potem narisemo poligon, drugace samo krivuljo
function X = narisiBezierCastel(kTocke,N,poligon,krivulj)
	n=length(kTocke(1,:));    %Izracunamo stevilo kontrolnih tock
	u=linspace(0,1,N);	 %Izracunamo x koordinate tock na krivulji 
	%de Casteljauvova formula
	Xx=(1-u) .^3 *kTocke(1,1) + 3*(1-u).^2.*u.*kTocke(1,2) + 3*(1-u).*u.^2.*kTocke(1,3) + u.^3.*kTocke(1,4);
	Xy=(1-u) .^3 *kTocke(2,1) + 3*(1-u).^2.*u.*kTocke(2,2) + 3*(1-u).*u.^2.*kTocke(2,3) + u.^3.*kTocke(2,4);
	X=[Xx;Xy];
	if(krivulj>0)
		plot(X(1,:), X(2,:));	%Narisemo se krivuljo
		hold on;
	end
	if(poligon==1)
		plot(kTocke(1,:), kTocke(2,:), ['-o','g']); %Narisemo poligon
	end
endfunction

%Funkcija izracuna prvi 2 kontrolni tocki naslednje bezierjeve krivulje in jih vrne
%Kot argumente ji podamo prejsnje kontrolne tocke in pa x ali y 2. kontrolne tocke
%Iz prejsnjih kontrolnih tock poiscemo premico, ki gre skozi zadnji dve kontrolni tocki
%iz premice pa glede na podan x ali y izracunamo tocko
%kaj je podano x ali y povemo tako, da je ena vrednost razlicna od Inf (neskoncno)
%V primeru napake, da podamo napacne argumente, sta oba vrnjena vektorja prazna
%Parametra x in y se kontrolirata v bloku PREVERJANJE
function [tocki] = naslednjiKontrolniTocki(kontrole,x,y)
	%koeficient premice: Potrebujemo zadnji dve kontrolni tocki z indeksoma 3 oz. 4
	%k=(y2-y1)/(x2-x1)
	kX3=kontrole(1,3); kY3=kontrole(2,3);
	kX4=kontrole(1,4); kY4=kontrole(2,4);
	
	if(kY4==kY3)	%Zaradi deljenja z 0, ce dobimo vodoravno premico
		k=0;
	elseif(kX4==kX3) %Premica je navpična
		k=Inf;
	elseif(1==1)
		k=(kY4-kY3)/(kX4-kX3);
	endif
	%izracunamo se n iz formule za premico y=k*x+n  -> n=y-k*x
	n=kY4-k*kX4;
	n1=kY3-k*kX3;
	%---------------PREVERJANJE------------- "zaradi nadzora parametrov"
	%seveda veljajo določene omejitve pri izbiri x-a oz. y-a
	%Primer: če je k<0 -> premica pada -> v primeru da je 3. kontrolna točka višje od 4. na prejšnji krivulji:
	%y mora biti manjši od zadnje kontrolne točke in x večji, če ne dobimo nesmiselno situacijo,
	%Ne sme veljati niti enakost, ker sta potem prvi dve kontrolni točki nove krivulje enaki in spet ne dobimo odvedljive krivulje
	dx=abs(kX3-kX4);
	dy=abs(kY3-kY4);
	if (k<0)
		if (kY3<kY4)	%oziroma kX3>kX4
			if(x==Inf&&y==Inf)
				x=kX4-dx;
				"1"
			elseif(x!=Inf && x >= kX4)
				"Neveljavno izbran x, k<0, err 1"
				return
			elseif(y!=Inf && y <= kY4)
				"Neveljavno izbran y, k<0, err 2"
				return
			endif
		
		else
			if(x==Inf&&y==Inf)
				x=kX4+dx;
				"2"
			elseif (x!=Inf && x <= kX4)
				"Neveljavno izbran x, k<0, err 3"
				return
			elseif (y!=Inf && y >= kY4)
				"Neveljavno izbran y, k<0, err 4"
				return
			endif
		endif
	elseif (k > 0 && k!=Inf)
		if (kY3<kY4)	%oziroma kX3<kX4
			if(x==Inf&&y==Inf)
				x=kX4+dx;
				"3"
			elseif (x!=Inf && x <= kX4)
				"Neveljavno izbran x, k>0, err 5"
				return
				
			elseif (y!=Inf && y <= kY4)
				"Neveljavno izbran y, k>0, err 6"
				return
			endif
		
		else
			if(x==Inf&&y==Inf)
				x=kX4-dx;
				"4"
			elseif (x!=Inf && x >= kX4)
				"Neveljavno izbran x, k>0, err 7"
				return
			elseif (y!=Inf && y >= kY4)
				"Neveljavno izbran y, k<0, err 8"
				return
			endif
		endif
	elseif(k==0)
		if(x==Inf&&y==Inf)
			if(kX3>kX4)
				x=kX4-dx;
			else
				x=kX4+dx;
			end
		elseif (y!=Inf)	%Premica je vodoravna, k=0, ne moremo iskati x ker imamo neskončno rešitev, lahko pa iščemo y, trivialna rešitev, y=n
			"Premica je vodoravna, k=0, neznan x, err 9"
			return
		elseif(kX3<kX4 && x<=kX4)
			"vodoravna premica, k=0, slab x, err 10"
			return
		elseif(kX4<kX3 && x>=kX4)
			"vodoravna premica, k=0, slab x, err 11"
			return
		endif	
		
	elseif(k==Inf)	%k==Inf, navpična premica
		if(x==Inf&&y==Inf)
			if(kY4>kY3)
				y=kY4+dy;
			else
				y=kY4-dy;
			end
		elseif(x!=Inf)
			"Navpična premica, x mora biti nedefiniran err 12"
			return
		elseif(kY3<kY4 && y<=kY4)
			"navpična premica, k=0, slab y, err 10"
			return
		elseif(kY4<kY3 && y>=kY4)
			"navpična premica, k=0, slab y, err 11"
			return
		endif	
		
	endif
	%--------------Konec PREVERJANJE-----------
	
	if(k==Inf)	%Če je premica navpična, velja enačba za premico x=...
		x=y;
	elseif(x!=Inf)	
		y=k*x+n;
	elseif (y!=Inf)
		x=(y-n)/k;
	else
		"Napačni parametri!"
	endif
	tocki=[kX4,x;kY4,y];
endfunction