function SegmentationCCTensor()

addpath D:\MatlabCodes

cd 'D:\EP2026';

%_______________________Parameters__________________
global U Y_Signal Y_ODF;
global M N N_Slice mi ma FA FA_Image_Mask pds pds1 pds2 pds3 pds4 NTSP0 PassHistory M0 N0 Ns0 Speed_F_Thresh Collinearity_Thresh Speed_factor FA_Seg_Thresh
global h1 h2 h3 h4 h5 h6 N_SH DTI_Tag BundleIndex BundleIndexThresh;
global NormalDirection Ten_Sim_Thre;
global notBundleIndex   notBundleIndexThresh
warning off;

FA_threshold=0.1;

%----------------------CC-------------------------
FiberName='CC_';

BundleIndex=1;
BundleIndexThresh= 0.5;

notBundleIndex=2;
notBundleIndexThresh= 0.8; %0.7;


r=120;
Speed_factor=1;
DTI_Tag=1;

Collinearity_Thresh =0.8;

FA_Seg_Thresh=0.1; %0.1;%;
%0.01; %0.05;%0.01;%0.15;%0.3; ok
Speed_F_Thresh =0.1;  %0.1;
%0.03; %0.05;% 0.04;%0.05;%.1;%0.2;

Ten_Sim_Thre=0.1;  %

% % %----------------------ATR-------------------------
% FiberName='ATR_';
%
% BundleIndex=2;
% BundleIndexThresh= 0.8;
%
%
% notBundleIndex=2;
% notBundleIndexThresh= 0.8; %0.7;
%
%
% r=120;
% Speed_factor=1;
% DTI_Tag=1;
%
% Collinearity_Thresh =0.9;
%
% FA_Seg_Thresh=0.05;%;
% %0.01; %0.05;%0.01;%0.15;%0.3; ok
% Speed_F_Thresh =0.1;
% %0.03; %0.05;% 0.04;%0.05;%.1;%0.2;
%
% Ten_Sim_Thre=0.2;   %0.35;%0.4; %
%
%

% % % %----------------------Cg-------------------------
% FiberName='Cg_';
% BundleIndex=2;   %1;
% BundleIndexThresh=0.8; %0.7;  %0.45;% 0.35;%0.6;  ok
%
%  notBundleIndex=1;
% notBundleIndexThresh= 0.3; %0.7;
% %
%
% r=90;
% Speed_factor=1;
% DTI_Tag=1;
%
% Collinearity_Thresh =0.7; %0.95;%0.7;    ok
% FA_Seg_Thresh=0.05;%0.1%0.15;%0.3; ok
% Speed_F_Thresh =0.05;%0.1;%0.05;%.1;%0.2;
%
% Ten_Sim_Thre=0.2;   %0.35;%0.4; %

% % % %----------------------Fornix-------------------------
% FiberName='fx_';
% BundleIndex=2;   %1;
% BundleIndexThresh=0.8;%0.8; %0.7;  %0.45;% 0.35;%0.6;  ok
% r=90;
% Speed_factor=1;
% DTI_Tag=1;
%
% Collinearity_Thresh =0.7; %0.95;%0.7;    ok
% FA_Seg_Thresh=0.2;%0.05;%0.01;%0.05;%0.1%0.15;%0.3; ok
% Speed_F_Thresh =0.2;%0.05;%0.01;%0.05;%0.1;%0.05;%.1;%0.2;
%
% Ten_Sim_Thre=0.2;   %0.35;%0.4; %


% %--------------------CST-----------------------------
% FiberName='CST';
% BundleIndex=3;
% BundleIndexThresh=0.6;
% r=120;
% Speed_factor=1;
% DTI_Tag=1;
% Speed_F_Thresh =0.3;     %0.6:GlauberHydi120NEX1Smoothed;
% Collinearity_Thresh =0.7;     %0.95:GlauberHydi120NEX1Smoothed;
% FA_Seg_Thresh=0.3;
% UseCl=0;
%%-----------------------------------------------------



%___________Load DTI__________________
[filename, pathname] = uigetfile('DTI.mat','Load DTI File');
exist_f = isa(filename,'char');
if exist_f
    Inputfile = fullfile(pathname,filename);
else
    return;
end
cd (pathname);
load (Inputfile);
try
    DT=Full_Rotated_Nei_Tensor;
    FA=Full_Rotated_Nei_FA;
    
    FE=Full_Rotated_Nei_FE;
    
catch
    %     a=FEFA(:,:,:,1);
    %     FEFA(:,:,:,1)=FEFA(:,:,:,2);
    %     FEFA(:,:,:,2)=a;
    %
    %     a=FE(:,:,:,1);
    %     FE(:,:,:,1)=FE(:,:,:,2);
    %     FE(:,:,:,2)=a;
    %
    %     FE(:,:,:,3)=-FE(:,:,:,3);
end



FA(isnan(FA))=0;
AnisoImage =FA ;
[M N N_Slice] = size(FA);

DT=permute(DT, [4,1,2,3]);
TensorAll=double(reshape(DT,9,M*N*N_Slice));
TensorAll(:,FA <= FA_threshold)=[];

FA_Image_Mask=find(FA>FA_threshold);

FEFA(:,:,:,1)=FE(:,:,:,1).*FA;
FEFA(:,:,:,2)=FE(:,:,:,2).*FA;
FEFA(:,:,:,3)=FE(:,:,:,3).*FA;
FEFA=abs(FEFA);

pds=permute(FE, [4,1,2,3]);
pds=double(reshape(pds,3,M*N*N_Slice));

pds(:,FA <= FA_threshold)=[];
RemainingPDDs_All=zeros(12,length(FA_Image_Mask));
RemainingPDDs_All(1:3,:)=pds;
no_RemainingPDDs_All=ones(length(FA_Image_Mask),1);

pds1=pds;
pds2=pds;
pds3=pds;
pds4=pds;

MaxODF_FA=zeros(M*N*N_Slice,3);

for i=1:3
    pds_fa(:,i)=abs(pds(i,:)').*FA(FA_Image_Mask);
end
MaxODF_FA(FA_Image_Mask,:)=pds_fa;
MaxODF_FA=reshape(MaxODF_FA,M,N,N_Slice,3);


%___________Load Input ROI__________________
[filename, pathname] = uigetfile('ROI.mat','Load InputROI File');
exist_f = isa(filename,'char');
if exist_f
    Inputfile = fullfile(pathname,filename);
else
    return;
end

load (Inputfile);

FileName='FiberBundle';

close all;


Seeds=find(ROI~=0);
[ni nj nk]=ind2sub([M N N_Slice],Seeds);
Seeds=[ni nj nk]';
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);


% Make sure we can see the kernel m-files.
% run('../addPathToKernel');

% path(path, 'd:/work/')
%---------------------------------------------------------------------------
% 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


radius = 0.01;
data1 = zeros(size(g.xs{1}));
for i = 1 : g.dim
    data1 = data1 + (g.xs{i} -centers(i,1)).^2;
end
data1 = sqrt(data1) - radius;
data=data1;

for j=2:size(centers,2)
    
    data1 = zeros(size(g.xs{1}));
    for i = 1 : g.dim
        data1 = data1 + (g.xs{i} -centers(i,j)).^2;
    end
    data1 = sqrt(data1) - radius;
    data(data1<=0.01)=data1(data1<=0.01);
end

clear data1

schemeData.TensorAll=TensorAll;

[X,Y,Z] = meshgrid(1:M,1:N,1:N_Slice) ;
[faces,Boundry]=(isosurface(Y,X,Z, data, level));


VolumeBoundry = round(Boundry');

[fq,bin]=hist(VolumeBoundry(1,:),min(VolumeBoundry(1,:)):max(VolumeBoundry(1,:))); [mx,ind]=max(fq); M0=bin(ind);
[fq,bin]=hist(VolumeBoundry(2,:),min(VolumeBoundry(2,:)):max(VolumeBoundry(2,:))); [mx,ind]=max(fq); N0=bin(ind);
[fq,bin]=hist(VolumeBoundry(3,:),min(VolumeBoundry(3,:)):max(VolumeBoundry(3,:))); [mx,ind]=max(fq); Ns0=bin(ind);

neigb=[-1,0,1,-1,0,1,-1,0,1,-1,0,1,-1,1,-1,0,1,-1,0,1,-1,0,1,-1,0,1,0;-1,-1,-1,0,0,0,1,1,1,-1,-1,-1,0,0,1,1,1,-1,-1,-1,0,0,0,1,1,1,0;-1,-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,0];
NeiVolumeBoundry=unique((repmat(VolumeBoundry,1,size(neigb,2))+repmat(neigb,1,size(VolumeBoundry,2)))','rows')';
vz=NeiVolumeBoundry(1,:)<1 | NeiVolumeBoundry(2,:)<1 | NeiVolumeBoundry(3,:)<1 | NeiVolumeBoundry(1,:)>M | NeiVolumeBoundry(2,:)>N | NeiVolumeBoundry(3,:)>N_Slice;
NeiVolumeBoundry(:,vz)=[];

neigb=[0;0;0];
VolumeBoundry=unique((repmat(VolumeBoundry,1,size(neigb,2))+repmat(neigb,1,size(VolumeBoundry,2)))','rows')';
vz=VolumeBoundry(1,:)<1 | VolumeBoundry(2,:)<1 | VolumeBoundry(3,:)<1 | VolumeBoundry(1,:)>M | VolumeBoundry(2,:)>N | VolumeBoundry(3,:)>N_Slice;
VolumeBoundry(:,vz)=[];
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');
h1=figure('Position',[1 1 scrsz(3)/3 scrsz(4)]);

scrsz = get(0,'ScreenSize');
h2=figure('Position', [scrsz(3)/3+5 1 scrsz(3)/3-5 scrsz(4)]);

scrsz = get(0,'ScreenSize');
h3=figure('Position',[2*scrsz(3)/3+5 1 scrsz(3)/3-5 scrsz(4)]);
%
% h4=figure;
% h5=figure;


%---------------------------------------------------------------------------
% Loop until tMax (subject to a little roundoff).
count=1;
startTime = cputime;
NormalDirection=zeros(M,N,N_Slice,3);

DiffData=10000;
itr=1;
while (DiffData > 1)
    %         for itr=1:N_itr
    itr=itr+1
    
    dataOLD=data;
    data = SimilarityMeasure(data, schemeData); %, SHData)
    DiffData=norm(data(:)-dataOLD(:));
    
    [faces,Boundry]=(isosurface(Y,X,Z, data, level));
    
    VolumeBoundry=round(Boundry');
    
    
    
    VolumeBoundryShow =VolumeBoundry;
    
    neigb=[-1,0,1,-1,0,1,-1,0,1,-1,0,1,-1,1,-1,0,1,-1,0,1,-1,0,1,-1,0,1,0;-1,-1,-1,0,0,0,1,1,1,-1,-1,-1,0,0,1,1,1,-1,-1,-1,0,0,0,1,1,1,0;-1,-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,0];
    NeiVolumeBoundry=unique((repmat(VolumeBoundry,1,size(neigb,2))+repmat(neigb,1,size(VolumeBoundry,2)))','rows')';
    vz=NeiVolumeBoundry(1,:)<1 | NeiVolumeBoundry(2,:)<1 | NeiVolumeBoundry(3,:)<1 | NeiVolumeBoundry(1,:)>M | NeiVolumeBoundry(2,:)>N | NeiVolumeBoundry(3,:)>N_Slice;
    NeiVolumeBoundry(:,vz)=[];
    MaxODF_FA=zeros(M*N*N_Slice,3);
    
    for i=1:3
        pds_fa(:,i)=abs(pds(i,:)').*FA(FA_Image_Mask);
    end
    MaxODF_FA(FA_Image_Mask,:)=pds_fa;
    MaxODF_FA=reshape(MaxODF_FA,M,N,N_Slice,3);
    
    %__________For offline PDD estimation---------------------
    
    neigb=[0;0;0];
    VolumeBoundry=unique((repmat(VolumeBoundry,1,size(neigb,2))+repmat(neigb,1,size(VolumeBoundry,2)))','rows')';
    vz=VolumeBoundry(1,:)<1 | VolumeBoundry(2,:)<1 | VolumeBoundry(3,:)<1 | VolumeBoundry(1,:)>M | VolumeBoundry(2,:)>N | VolumeBoundry(3,:)>N_Slice;
    VolumeBoundry(:,vz)=[];
    schemeData.Boundry=VolumeBoundry;
    
    pause(1);
    
end

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


figure(h1);

imshow((im));

%Coronal

offset=fix(N/10);

L2=fix(sqrt((N-offset)/2))+2;

L1=fix((N-offset)/(2*L2))-1;

% L1=fix(sqrt(N/2));
% L2=fix(N/2/L1)+1;


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);
        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

figure(h2);
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


figure(h3);

imshow((im));

ROI=data;
ROI(ROI>0)=r;
ROI(ROI<=0)=0;

FiberImage=ROI;

FiberVoxelSize=VDims;

% [filename, pathname] = uiputfile([ FiberName,'RG_BundleThresh',num2str(  BundleIndexThresh),'Collinearity_Thresh',num2str( Collinearity_Thresh),'FA_Seg_Thresh',num2str(FA_Seg_Thresh),'_Speed_F_Thresh', num2str(Speed_F_Thresh)], 'Save Segmented as...');
[filename, pathname] = uiputfile([ FiberName 'Segmented'], 'Save Segmented as...');

exist_f = isa(filename,'char');
if exist_f
    targetfile = fullfile(pathname,filename);
else
    return;
end

% FiberImage=fiberImage;
Analyze = Write_Analyze (targetfile, FiberImage ,'int16', FiberVoxelSize);


% strFile=[pathname,'/', FiberName,'RG_BundleThresh',num2str(  BundleIndexThresh),'Collinearity_Thresh',num2str( Collinearity_Thresh),'FA_Seg_Thresh',num2str(FA_Seg_Thresh),'_Speed_F_Thresh', num2str(Speed_F_Thresh), '.mat'];
strFile=[pathname,'/', FiberName,'Segmented.mat'];
save(strFile, 'data','FiberImage', 'FiberVoxelSize');
strFigFile=[strFile,'_Axl.Fig'];
hgsave(h1,strFigFile);
strFigFile=[strFile,'_Cor.Fig'];
hgsave(h2,strFigFile);
strFigFile=[strFile,'_Sag.Fig'];
hgsave(h3,strFigFile);

% hgload (strFigFile);

endTime = cputime;
fprintf('Total execution time %g seconds', endTime - startTime);


display('Done!')


end

function out = SimilarityMeasure(data, schemeData)

global M N N_Slice M0 N0 Ns0 FA_Image_Mask NTSP0 pds1 pds2 pds3 pds4 pds FA Collinearity_Thresh Speed_F_Thresh mi ma h3 h4 h5 h6
global N_SH DTI_Tag notBundleIndex BundleIndex BundleIndexThresh notBundleIndexThresh FA_Seg_Thresh Ten_Sim_Thre;
%
DTData=schemeData.TensorAll;
v=2:N_SH;

DataSize=[M N N_Slice];

GPhi=zeros(DataSize(1),DataSize(2),DataSize(3),3);
ImageBoundry=schemeData.Boundry;

data2=zeros(size(data));

data2(find(data<=0))=data(find(data<=0));

for i = 1 : 3
    [ derivL, derivR ]=upwindFirstWENO5a(schemeData.grid, data, i, true); %  upwindFirstWENO5a OK
    GPhi(:,:,:,i) = derivL{3};
end

NormalDirection=zeros(DataSize(1),DataSize(2),DataSize(3),3);
CurrImage=zeros(DataSize(1),DataSize(2),DataSize(3),3);
NextImage=zeros(DataSize(1),DataSize(2),DataSize(3),3);
PrevImage=zeros(DataSize(1),DataSize(2),DataSize(3),3);
SecondPrevImage=zeros(DataSize(1),DataSize(2),DataSize(3),3);
PrevImageP=zeros(DataSize(1),DataSize(2),DataSize(3),3);
SecondPrevImageP=zeros(DataSize(1),DataSize(2),DataSize(3),3);

neig=[-1,0,1,-1,0,1,-1,0,1,-1,0,1,-1,1,-1,0,1,-1,0,1,-1,0,1,-1,0,1,0;-1,-1,-1,0,0,0,1,1,1,-1,-1,-1,0,0,1,1,1,-1,-1,-1,0,0,0,1,1,1,0;-1,-1,-1,-1,-1,-1,-1,-1,-1,0,0,0,0,0,0,0,0,1,1,1,1,1,1,1,1,1,0];
neig(:,end)=[];

NTSP1=NTSP0;
for m=1:size(ImageBoundry,2)
    i=ImageBoundry(1,m);
    j=ImageBoundry(2,m);
    k=ImageBoundry(3,m);
    Normal=reshape(GPhi(i,j,k,:),3,1);
    Normal=-Normal/norm(Normal);
    NormalDirection(i,j,k,:)=Normal;
    ind=sub2ind([M N N_Slice],i,j,k);
    currentPoint=[i j k]' ;
    CurrImage(i,j,k,:)=[1 0 0];
    index2=find(FA_Image_Mask==ind);
    temp_val=0;
    if ~isempty(index2)
        pds_i=pds(:,index2);
        DTData_i=reshape(DTData(:,index2),3,3);
        BoundaryPoints=repmat(currentPoint,1,size(neig,2))+neig;
        for c=1: size(BoundaryPoints,2)
            BoundaryPoint=BoundaryPoints(:,c);
            if sum((BoundaryPoint-currentPoint).*Normal) <=0
                BoundaryPointLabel='Prev';
                PreviousPoint=BoundaryPoint;
            else
                BoundaryPointLabel='Next';
                NextPoint=BoundaryPoint;
            end
            if BoundaryPointLabel == 'Next' & ( (NextPoint(1) <= DataSize(1)) && (NextPoint(2) <= DataSize(2)) && (NextPoint(3) <= DataSize(3)) &&   NextPoint(1) >=1  && NextPoint(2) >=1 &&   NextPoint(3) >=1  ) %#ok<STCMP>
                ind=sub2ind([M N N_Slice],NextPoint(1),NextPoint(2),NextPoint(3));
                index2=find(FA_Image_Mask==ind);
                if ~isempty(index2)
                    NextImage(NextPoint(1),NextPoint(2),NextPoint(3),:)=[0 0 1];
                    DTData_ip1=reshape(DTData(:,index2),3,3);
                    fa_ip1 = FA(index2);
                    DTData_ip1=reshape(DTData(:,index2),3,3);
                    pds1_ip1=pds1(:,index2);
                    tempMax=abs(sum(pds_i.*pds1_ip1));%*(abs(pds1_ip1(1)));
                    DirectionalTerm=tempMax;
                    
                    FA_Sim =FA(i,j,k)*FA(NextPoint(1),NextPoint(2),NextPoint(3));
                    Ten_Sim=trace (DTData_i*DTData_ip1 )/(trace (DTData_i)*trace(DTData_ip1) );
                    
                    
                    SpecificBundleCond =( abs(pds1_ip1(BundleIndex)) > BundleIndexThresh ); % && (abs(pds1_ip1(notBundleIndex)) < notBundleIndexThresh ) ;
                    
                    if  Ten_Sim > Ten_Sim_Thre && sum(isnan(DTData_ip1(:)))==0  && trace (DTData_ip1)~=0 && SpecificBundleCond && (tempMax>=Collinearity_Thresh) && FA(NextPoint(1),NextPoint(2),NextPoint(3)) > FA_Seg_Thresh %            if acos(abs (sum(pds_i.*tempMax ))/(norm(pds_i)*norm(tempMax))) < 10*pi/180
                        
                        
                        temp_val  = temp_val + FA_Sim*Ten_Sim;  %  =1 doesn't matter
                        
                    end
                    
                end
            end
        end
        NTSP1(i,j,k )=temp_val/size(BoundaryPoints,2);
    end
end
if mi>min(NTSP1(:))
    mi=min(NTSP1(:));
end
if ma<max(NTSP1(:))
    ma=max(NTSP1(:));
end
NTSP0=NTSP1;
Speed_F=NTSP1;
Speed_F=( Speed_F- mi )/( ma-mi);
Speed_F(Speed_F < Speed_F_Thresh )=0;
out = Speed_F;

end
