Untitled

%---------------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