Untitled

%---------------GLAVNI PROGRAM------------------
zacetneTocke=[1 1.5 2 3 4 5 4; 1 1.5 2 2.5 3 2 0.5] %Zgornja vrstica so x koordinate spodnja pa y
%Narisemo zacetne tocke
%plot(zacetneTocke(1,:),zacetneTocke(2,:),['o','b'])
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])

kT=[1 0 2 3; 1 4 4 0]
X=bezier(kT,200,0)
%aD=avgDistance(zacetneTocke,X)
%hold off
%l=bezierGetControlPoints([1 2 5 4; 1 4 5 1],[0.001 1-0.001],1)
%hold on
%X=bezier(l,200,1)
%hold on
[pl,gh]=nextControlPoint(kT,4,Inf);
%plot(gh(1,:),gh(2,:),['-o','r'])
newKont=[gh,[4,6; -1 0.5]]
X=bezier(newKont,200,0)

%---------------------------------
%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 = avgDistance(tocke,X)
        stTock=length(tocke(1,:));      %Stevilo danih tock
        avgDist=0;
        for k=1:stTock
                tocka=tocke(:,k)';
%----------------!!!PREVERI CE JE PRAVILNO, mislim da nekaj ne stima!!!-------------
                [dist,i]=min(sqrt((X(1,:)-tocka(1)).^2 + (X(2,:)-tocka(2)).^2));        %Racunanje minimalne razdalje (evklidove)
                avgDist+=dist/stTock;   %Utezena vsota
        end

        %plot([tocka(1),X(1,i)],[tocka(2),X(2,i)],['-o','r']); %
endfunction

%Funkcija izracuna Bezierjevo krivuljo, ne nujno da 3. reda in jo vrne v obliki N tock
%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 = bezier(kTocke,N,poligon)
        if(poligon==1)
                plot(kTocke(1,:), kTocke(2,:), ['-o','g']); %Narisemo poligon
                hold on;
        end
        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)]';

        plot(X(1,:), X(2,:));   %Narisemo se krivuljo

endfunction

%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
function controls = bezierGetControlPoints(tocke,verjetnost,izris)
        %imamo tocke t0, t1, t2 in t3 v vektorju tocke
        %POMOC: http://polymathprogrammer.com/2007/06/27/reverse-engineering-bezier-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
        razd=tocke(1,1)+tocke(length(tocke))
        a=[0 1/5 4/5 1]
        u=verjetnost(1); v=verjetnost(2);
        %Kontrolni tocki k0=t0 in k3=t3 imamo ze, rabimo se k1 in k2
        %Uporabimo de Casteljejevo 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:
        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];

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)
function [kont,tocki] = nextControlPoint(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);

        k=(kY4-kY3)/(kX4-kX3);
        %izracunamo se n iz formule za premico y=k*x+n  -> n=y-k*x
        n=kY4-k*kX4;
        n1=kY3-k*kX3;
        if(n==n1)
                "n-ja sta enaka, vse je kul"
        endif
        if(x!=Inf)
                y=k*x+n;
        elseif(y!=Inf)
                x=(y-n)/k;
        else
                "NAPAKA"
        endif
        tocki=[kX4,x;kY4,y];
        kont=[kX3,kX4;kY3,kY4];
endfunction