%Anton Ullrich
%multiphase MCF via Voronoi
function MCFVoronoi
    it=200;
    fig = uifigure('Name', 'Voronoi');
    fig.Position(3)=1500;
    fig.Position(4)=768;
    global X;
    global Y;
    global h;
    h=100;%200 weiter als 100
    n=1000;
    [X,Y] = meshgrid(-1:.01:1);
    hAxes=axes('Parent',fig);
    set(hAxes,'Position',[.05 .05 .9 .9]);
    %set(hAxes,'nextplot','replacechildren'); 
    vid = VideoWriter('voronoi.avi');
    open(vid);
    Z=zeros(size(X));
    points=2*rand(n+4,2)-1;
    points(1:4, :)=[-10 -10; -10 10; 10 -10; 10 10];
    [v,c] = voronoin(points);
    m=length(c)-4;
    A = zeros([m+1,size(X)]);
    A(m+1,:,:)=ones(size(X));
    for i = 1:m
        v1 = v(c{i+4},1) ; 
        v2 = v(c{i+4},2) ;
        A(i,:,:)=inpolygon(X,Y,v1,v2);
        A(m+1,A(i,:,:)==1)=0;
        %Z(A(i,:,:)==1)=i;
    end
    for i = 1:m
        %A(i,:,:)=zeros(size(X));
        %A(i,Z==i)=1;
        XT=X(A(i,:,:)==1);
        YT=Y(A(i,:,:)==1);
        k = boundary(XT,YT,.9);
        hold(hAxes,'on');
        plot(hAxes,XT(k),YT(k),'black');
        hold(hAxes,'off');
    end
    Adiff = zeros([m+1,size(X)]);
    Adiff(m+1,:,:)=A(m+1,:,:);
    Adiff(m+1,1:10,:)=2;
    Adiff(m+1,end-10:end,:)=2;
    Adiff(m+1,:,1:10)=2;
    Adiff(m+1,:,end-10:end)=2;
    for i=1:it
        cla(hAxes);
        for j = 1:m
            if(~isempty(A(j,11:end-11,11:end-11)))
                Adiff(j,:,:)=reshape(convolve(A(j,:,:)),[1, size(X)]);%charis(Z, j)
            end
        end
        [M,index]=max(Adiff);
        %index
        for j = 1:m
            A(j,:,:)=zeros(size(X));
            A(j,index==j)=1;
            XT=X(A(j,:,:)==1);
            %YT=Y(index==j+1);
            YT=Y(A(j,:,:)==1);
            k = boundary(XT,YT,.9);
            hold(hAxes,'on');
            plot(hAxes,XT(k),YT(k),'black');
            hold(hAxes,'off');
        end
        frame = getframe(hAxes);
        writeVideo(vid,frame);
    end
    close(vid);
    5
end

function Z=convolve(Zold)
    global h;
    global X;
    global Y;
    Zold=reshape(Zold, size(X));
    Z = abs(ifft2(fftshift(exp((-(X).^2-(Y).^2).*h/2).*fftshift(fft2(Zold)))));%./(2*pi)
    %real or abs nicht notwendig aber um ggf. translationsfehler zu
    %beheben, da nichtnegativ
end

function Z=charis(Zold,i)
    Z=zeros(size(Zold));
    oneArray=ones(size(Zold));
    Z(Zold==i)=oneArray(Zold==i);
end