
% Integration parameters.
tMax = 1.0;                  % End time.
plotSteps = 6;               % How many intermediate plots to produce?
t0 = 0;                      % Start time.
singleStep = 0;              % Plot at each timestep (overrides tPlot).
% Period at which intermediate plots should be produced.
tPlot = (tMax - t0) / (plotSteps - 1);
% How close (relative) do we need to get to tMax to be considered finished?
small = 10 * eps;
%---------------------------------------------------------------------------
% What level set should we view?
level = 0;
% Pause after each plot?
pauseAfterPlot = 0;
% Delete previous plot before showing next?
deleteLastPlot = 0;
% Plot in separate subplots (set deleteLastPlot = 0 in this case)?
useSubplots = 0;
%---------------------------------------------------------------------------
% Use periodic boundary conditions?
periodic = 0;

mi=Inf;
ma=0;

g.dim = 3;
g.min = -1;
g.N=[M N N_Slice]';

NTSP0=zeros(M,N,N_Slice);
PassHistory=zeros(M,N,N_Slice);

if(periodic)
    g.max = (1 - g.dx);
    g.bdry = @addGhostPeriodic;
else
    g.max = +1;
    g.bdry = @addGhostExtrapolate;
end
g = processGrid(g);

%---------------------------------------------------------------------------
% What kind of display?
if(nargin < 3)
    switch(g.dim)
        case 1
            displayType = 'plot';
        case 2
            displayType = 'contour';
        case 3
            displayType = 'surface';
        otherwise
            error('Default display type undefined for dimension %d', g.dim);
    end
end
% Seeds(:,3)=[];%=Seeds(:,1:2);
centers=[];
centers(1,:)=2/(M-1)*Seeds(1,:)-(M+1)/(M-1);
centers(2,:)=2/(N-1)*Seeds(2,:)-(N+1)/(N-1);
centers(3,:)=2/(N_Slice-1)*Seeds(3,:)-(N_Slice+1)/(N_Slice-1);

% data = zeros(size(g.xs{1}));
[X,Y,Z] = meshgrid(1:M,1:N,1:N_Slice) ;

schemeData.TensorAll=TensorAll;

radius = 0.01;
% neigb=[ -1 -1 -1 0 0 1 1 1;0 0 0 0 0 0 0 0; -1 0 1 -1 1 -1 0 1];
neig=[ -1 -1 -1 0 0 1 1 1; -1 0 1 -1 1 -1 0 1; 0 0 0 0 0 0 0 0];
neigb=[0;0;0];
% neigb=[0;0;0];
VolumeBoundry=[];
data=zeros(size(g.xs{1}));
LS=[];
for j=1:size(centers,2)
    
    LS(j).data = zeros(size(g.xs{1}));
    for i = 1 : g.dim
        LS(j).data = LS(j).data + (g.xs{i} -centers(i,j)).^2;
    end
    LS(j).data = sqrt(LS(j).data) - radius;
    data=data+LS(j).data;
    Boundry=[];
    [faces,Boundry]=(isosurface(Y,X,Z, LS(j).data, level));
    %  data(data1<=0.01)=data1(data1<=0.01);
    % imshow(reshape(data(:,65,:),128,68),[]);
    
    VB = round(Boundry');
    %     VB(2,:)=round(mean(VB(2,:)));   %Seeds(2,1)
    
    VB(3,:)=round(mean(VB(3,:)));
    
    %     VB=unique((repmat(VB,1,size(neigb,2))+repmat(neigb,1,size(VB,2)))','rows')';
    
    VB=repmat(VB,1,size(neigb,2))+repmat(neigb,1,size(VB,2));
    VB=unique(VB','rows')';
    vz=VB(1,:)<1 | VB(2,:)<1 | VB(3,:)<1 | VB(1,:)>M | VB(2,:)>N | VB(3,:)>N_Slice;
    VB(:,vz)=[];
    LS(j).VolumeBoundry=VB;
    LS(j).NTSP0=zeros(M,N,N_Slice);
    if isempty(VolumeBoundry)
        VolumeBoundry=LS(j).VolumeBoundry;
    else
        VolumeBoundry=[VolumeBoundry   LS(j).VolumeBoundry];
    end
    
end

% schemeData.Boundry=VolumeBoundry;
schemeData.grid = g;

%---------------------------------------------------------------------------
% Set up motion in the normal direction (derivative choice is set below).
schemeFunc = @termNormal;
schemeData.grid = g;
% schemeData.Data3d_SH=SH_from_MAT; %CrossingPattern_SH8;


scrsz = get(0,'ScreenSize');

%---------------------------------------------------------------------------
% Loop until tMax (subject to a little roundoff).

startTime = cputime;
NormalDirection=zeros(M,N,N_Slice,3);

DiffData=10000;
itr=1;

offs=50;
A=offs;
B=offs;
C=scrsz(3)*0.6/2;
D=scrsz(4)*0.8;

while (DiffData > 0.1)
    dataOLD=data;
    
    itr=itr+1
    data=zeros(size(g.xs{1}));
    VolumeBoundry=[];
    for j=1:size(centers,2)
        %         LS(j).dataOLD(:)=LS(j).data(:);
        
        try
            %             tmp= SimilarityMeasureCoronal(LS(j).data, schemeData ,  LS(j).VolumeBoundry , LS(j).NTSP0);
            tmp= SimilarityMeasure_7_11_11(LS(j).data, schemeData ,  LS(j).VolumeBoundry , LS(j).NTSP0,neig);
            LS(j).NTSP0 = tmp;
            %             data=data+tmp;
            [faces,Boundry]=isosurface(Y,X,Z, tmp,level );
            VB=round(Boundry');
            %             VB(2,:)=round(mean(VB(2,:)));   %Seeds(2,1)
            VB(3,:)=round(mean(VB(3,:)));   %Seeds(2,1)
            
            VB=unique((repmat(VB,1,size(neigb,2))+repmat(neigb,1,size(VB,2)))','rows')';
            vz=VB(1,:)<1 | VB(2,:)<1 | VB(3,:)<1 | VB(1,:)>M | VB(2,:)>N | VB(3,:)>N_Slice;
            VB(:,vz)=[];
            
            LS(j).VolumeBoundry=VB;
            LS(j).data=tmp;
            
            LS(j).diff= norm(LS(j).data(:)-dataOLD(:));
            
            if isempty(VolumeBoundry)
                VolumeBoundry=LS(j).VolumeBoundry;
            else
                VolumeBoundry=[ VolumeBoundry LS(j).VolumeBoundry];
            end
        catch
            
        end
    end
      
    
    DiffData=norm(data(:)-dataOLD(:));
       
end


% Axial

offset=fix(N_Slice/8);
L2=fix(sqrt((N_Slice-offset)/1));%+2;
L1=fix((N_Slice-offset)/L2)-1;

im=zeros(L1*M,L2*N,3);
impp=zeros(M,N,3);

for i1=1:L1
    for i2=1:L2
        
        vz=VolumeBoundry(3,:)==(i1-1)*L2+i2+offset; imp=reshape(FEFA(:,:,(i1-1)*L2+i2+offset,:),M,N,3);
        a=sub2ind([M,N,3],VolumeBoundry(1,vz),VolumeBoundry(2,vz),repmat(1,1,length(VolumeBoundry(3,vz)) ));
        imp(ind2sub([M,N,3], a))=2*max(imp(:));
        b=sub2ind([M,N,3],VolumeBoundry(1,vz),VolumeBoundry(2,vz),repmat(2,1,length(VolumeBoundry(3,vz)) ));
        imp(ind2sub([M,N,3], b))=2*max(imp(:));
        c=sub2ind([M,N,3],VolumeBoundry(1,vz),VolumeBoundry(2,vz),repmat(3,1,length(VolumeBoundry(3,vz)) ));
        imp(ind2sub([M,N,3], c))=2*max(imp(:));
        for i3=1:3
            impp(:,:,i3)=(imp(:,:,i3)');
        end
        
        im((i1-1)*M+1:i1*M,(i2-1)*N+1:i2*N,:)=impp;
        
    end
end


h1=figure('Position',[A B C D]);

imshow((im));

%Coronal

L2=11;
L1=5;
offset=36;

im=zeros(L2*N_Slice,L1*M,3);
impp=zeros(N_Slice,M,3);


for i1=1:L1
    for i2=1:L2
        
        %             vz=VolumeBoundry(2,:)==(2*i1-1)*L2+2*i2+offset; imp=reshape(FEFA(:,(2*i1-1)*L2+2*i2+offset,:,:),M,N_Slice,3);
        imp=reshape(FEFA(:,(1*i1-1)*L2+1*i2+offset,:,:),M,N_Slice,3);
        vz=VolumeBoundry(2,:)==(1*i1-1)*L2+1*i2+offset;
        a=sub2ind([M,N_Slice,3],VolumeBoundry(1,vz),VolumeBoundry(3,vz),repmat(1,1,length(VolumeBoundry(2,vz)) ));
        imp(ind2sub([M,N_Slice,3], a))=2*max(imp(:));
        b=sub2ind([M,N_Slice,3],VolumeBoundry(1,vz),VolumeBoundry(3,vz),repmat(2,1,length(VolumeBoundry(2,vz)) ));
        imp(ind2sub([M,N_Slice,3], b))=2*max(imp(:));
        c=sub2ind([M,N_Slice,3],VolumeBoundry(1,vz),VolumeBoundry(3,vz),repmat(3,1,length(VolumeBoundry(2,vz)) ));
        imp(ind2sub([M,N_Slice,3], c))=2*max(imp(:));
        for i3=1:3
            impp(:,:,i3)=flipud(imp(:,:,i3)');
        end
        im((i2-1)*N_Slice+1:i2*N_Slice,(i1-1)*M+1:i1*M,:)=impp;
    end
end


h2=figure('Position', [A+C+offs B C D ]);
imshow((im));

%Sagittal

offset=fix(M/3)+5;

L2=fix(sqrt((M-offset)/1))+0;
L1=fix((M-1*offset)/(2*L2));

im=[];
im=zeros(L2*N_Slice,L1*N,3);
impp=zeros(N_Slice,N,3);


for i1=1:L1
    for i2=1:L2
        
        vz=VolumeBoundry(1,:)==(1*i1-1)*L2+1*i2+offset; imp=reshape(FEFA((1*i1-1)*L2+1*i2+offset,:,:,:),N,N_Slice,3);
        a=sub2ind([N,N_Slice,3],VolumeBoundry(2,vz),VolumeBoundry(3,vz),repmat(1,1,length(VolumeBoundry(1,vz)) ));
        imp(ind2sub([N,N_Slice,3], a))=2*max(imp(:));
        b=sub2ind([N,N_Slice,3],VolumeBoundry(2,vz),VolumeBoundry(3,vz),repmat(2,1,length(VolumeBoundry(1,vz)) ));
        imp(ind2sub([N,N_Slice,3], b))=2*max(imp(:));
        c=sub2ind([N,N_Slice,3],VolumeBoundry(2,vz),VolumeBoundry(3,vz),repmat(3,1,length(VolumeBoundry(1,vz)) ));
        imp(ind2sub([N,N_Slice,3], c))=2*max(imp(:));
        for i3=1:3
            impp(:,:,i3)=flipud(imp(:,:,i3)');
        end
        
        im((i2-1)*N_Slice+1:i2*N_Slice,(i1-1)*N+1:i1*N,:)=impp;
    end
end

h3=figure('Position',[A+2*C+2*offs B C D  ]);
imshow((im));

pause(1);

% Seeds=[Right_Posterior Right_Anterior Left_Posterior Left_Anterior];


for j=1:size(centers,2)
    tmp=LS(j).data;
    
    %     [j length(find(tmp>0))]
    
    if length(find(tmp>0)) <1000
        
        if j<= size( Right_Posterior,2)    %--------------------------------------------------Right  Posterior
            ROI(tmp>0)=j;
            Label(tmp>0)=1;
            
        elseif j<= size( Right_Posterior,2)+size(Right_Anterior,2)    %------------------------Right  Anterior
            ROI(tmp>0)=size( Right_Superior,2) +j;
            Label(tmp>0)=3;
            
        elseif j<= size( Right_Posterior,2)+size(Right_Anterior,2)+size(Left_Posterior,2)     %Left  Posterior
            ROI(tmp>0)=size(Right_Superior,2)+size(Right_Inferior,2)+j;
            Label(tmp>0)=5;
            
        else  %---------------------------------------------------------------------------     % Left Anterior
            ROI(tmp>0)=size(Right_Superior,2)+size(Right_Inferior,2)+size(Left_Superior,2) +j;
            Label(tmp>0)=7;
            %             Analyze = UT_Write_Analyze (ROItargetfile, ROI ,'int16', VDims);
            
        end
    end
end

