close all;
clear all;
clc;
warning off; %#ok<WNOFF>
%_______________________Parameters_________________________________________


Scanner='GE';  %'Phillips' ;  %'Siemens';  % 'GE';  %
addpath D:\MatlabCodes
addr= 'C:\Data\MEG-DTI\EP1157\DTI-Results' ;

addpath addr; 

cd  (addr) 

WholeBrainInput=0;
WholeBrainOutput=1;%1; %0;

if WholeBrainOutput
    passes = true;
else
    passes =false;
end

FiberName='Fx';  %'Genu' %MEG

LoadROIMatfile=1;   %0 MEG
CutFibersOutsideROIs=1;   %1 MEG
SuccessiveLabels=0;
ReverseROI=0;

if WholeBrainInput   
    FiberName=[FiberName '_WholeBrainInput'];
elseif WholeBrainOutput   
    FiberName=[FiberName '_WholeBrainOutput'];
else
%     FiberName=[FiberName '_MultipleROIs'];
end

%--------------------fiber Tracking Parameters-------------
FA_threshold=0.001; %0.01;%0.1;
FT_FA_threshold=0.001; %0.01;%0.1;
EPSILON=0.0000001;
FibersCount=0;
maxNumFibers = 1000000;
minFiberLength = 1;
maxFiberAngle =35*pi/180;
anisoIndexStartThreshold = FT_FA_threshold;
anisoIndexStopThreshold = FT_FA_threshold;
DELTA_FACT = 0.05;

%___________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
load (Inputfile);

cd (pathname);

FA(isnan(FA))=0;
AnisoImage =FA ;
M = size(DWI,1);
N = size(DWI,2);
N_Slice= size(DWI,3);
NVolumes= size(DWI,4);
NrOfDirections=NVolumes-1; % excluding B0
DWI=permute(DWI, [4,1,2,3]);

FA_Image_Mask=find(FA>FA_threshold);
n_vox = length(FA_Image_Mask);
IsPDDsCalculated=zeros(n_vox,1);  %Online

% a=FE(:,:,:,1);
% FE(:,:,:,1)=FE(:,:,:,2);
% FE(:,:,:,2)=a;

% FE(:,:,:,2)=-FE(:,:,:,2);
FE(:,:,:,3)=-FE(:,:,:,3);
% FEFA(:,:,:,1)=-FEFA(:,:,:,1);
% FEFA(:,:,:,3)=-FEFA(:,:,:,3);

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

pds(:,FA <= FA_threshold)=[];

clear DWI;

%___________Load ROI__________________
if LoadROIMatfile
    
    [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);
%     cd (pathname);
        
else
    [filename, pathname] = uigetfile('ROI.img','Load ROI Img File');
    exist_f = isa(filename,'char');
    if exist_f
        FiberFileName = fullfile(pathname,filename);
    else
        return;
    end    
    
    [ROI, FiberImageSize , FiberImagePrecision, FiberVoxelSize] = Read_Analyze(FiberFileName);
    
end

% EP=FiberFileName(end-23:end-17);  %Normals
EP=pathname(end-14:end-9);

FiberName=[FiberName '_' EP];

NubmerOfROILabels=max(ROI(:));   

strFile=[FiberName,'_NOL',num2str(NubmerOfROILabels),'_'];
% strFile=[FiberName,'Hight_Coherence_']; %,num2str(NubmerOfROILabels),'_'];

if WholeBrainInput
    InputROI=ones(size(ROI));
end

dwiSize = MDims';
try
    dwiSize(4)=[];
catch
end
dwiSpacing =  VDims'; %[1 1 1]';
dwiSizeTimesSpacing = dwiSize.*dwiSpacing;
FIBERS = struct([]);

fiberImage = zeros(dwiSize(1),dwiSize(2),dwiSize(3));
r1=1;
r2=1;
NumberOfFibers=zeros(NubmerOfROILabels,NubmerOfROILabels);


if WholeBrainInput
    Seq =1;
    
elseif (ReverseROI && SuccessiveLabels)
    Seq = NubmerOfROILabels:-1:2 ;
elseif (~ReverseROI && SuccessiveLabels)
    Seq = 1:NubmerOfROILabels-1;
elseif (~SuccessiveLabels)
    Seq = 1:NubmerOfROILabels;
end
h_w = waitbar(0,'Finding tract...');
for r1= Seq
%     r1
    waitbar(r1/NubmerOfROILabels); 
    if WholeBrainInput
        
        OutputROI=zeros(size(ROI));
        
        OutputROI(ROI>=1)=1;
        InputROI(OutputROI~=0)=0;
    else
       
        InputROI=zeros(size(ROI));
        InputROI(ROI==r1)=ROI(ROI==r1);
           
        if SuccessiveLabels
            if ReverseROI
                r2=r1-1;
            else
                r2=r1+1;
            end
            OutputROI=zeros(size(ROI));
            OutputROI(ROI==r2)=ROI(ROI==r2);
        else
            OutputROI=ROI;
            OutputROI(InputROI~=0)=0;
        end
    end
    
    
    InputROI_Image_Mask=find(InputROI > 0 );
    
    n_InputROI_Vox = length(InputROI_Image_Mask);
    
    num=0;
    fiber = [];
    indexInputROI=0;
    
%     h_w = waitbar(0,'Finding tract...');
    
    while (indexInputROI < n_InputROI_Vox)
        
        indexInputROI=indexInputROI+1;
        
        
%         waitbar(indexInputROI/n_InputROI_Vox);
        %     [fOD_ODF PDDs] = fOD_ODF_PDDs_inThisVoxel(Signal_Data(InputROI_Image_Mask(indexInputROI)),lmax,Lambda,m,Precision,DIFF_C,e2_e1,eta,ODF_tag);
        [i j k]=ind2sub([M N N_Slice],InputROI_Image_Mask(indexInputROI));
        InitialPoint= [i j k]';
        
        if (  InitialPoint(1) > dwiSize(1) || InitialPoint(2) > dwiSize(2) || InitialPoint(3) > dwiSize(3))
            continue;
        end;
        
        
        
        if ( (isnan(AnisoImage(i,j,k)) ==0 ) && (AnisoImage(i,j,k) > anisoIndexStartThreshold) && (FibersCount < maxNumFibers) )
            
            % found a seed point
            
            
            currentPoint= [i j k]';
            
            if (  currentPoint(1) > dwiSize(1) || currentPoint(2) > dwiSize(2) || currentPoint(3) > dwiSize(3))
                continue;
            end;
            firstPoint = currentPoint;
            
            ind = round(currentPoint);
            if ( (isnan( AnisoImage(ind(1),ind(2),ind(3)) ) ) || (AnisoImage(ind(1),ind(2),ind(3)) < anisoIndexStopThreshold) )
                continue;
            end
            
            index1First= sub2ind([M N N_Slice],ind(1), ind(2), ind(3));
            index2First=find(FA_Image_Mask==index1First);
            
            if ( isempty(index2First))
                break;
            end
            
                if WholeBrainOutput
                    passes = true;
                else
                    passes =false;
                end
                
                pds_vox=pds(:,index2First);
                fiber = [];
                angle = 0;
                innerProd=1;
                currentPoint=firstPoint;
                firstpds_vox=pds_vox;
                oldpds_vox = pds_vox;
                
                %----------------------------------------------------------
                while ( sum(currentPoint < dwiSize)==3 &&  sum(currentPoint >=1)==3  && angle < maxFiberAngle )   %innerProd > 0.5 )  %
                    fiber = [fiber; currentPoint'];
                    
                    tmpPoint= currentPoint;
                    tmpPoint = tmpPoint + pds_vox*DELTA_FACT;
                    while( sum(tmpPoint < round(currentPoint)+1)==3 && sum(tmpPoint > round(currentPoint)-1 )==3 )
                        tmpPoint = tmpPoint + pds_vox*DELTA_FACT;
                    end
                    currentPoint = tmpPoint;
                    
                    ind = round(currentPoint);
                    
                    try
                        index= sub2ind([M N N_Slice],ind(1), ind(2), ind(3));
                    catch %#ok<CTCH>
                        break;
                    end
                    if ( isnan(index))
                        break;
                    end
                    
                    if ( (WholeBrainOutput || WholeBrainInput || OutputROI(index) ~= r1) && ( OutputROI(index) ~= 0) ) %  (OutputROI(index)~=InputROI(index))  )
                        passes = true;
                        r2=OutputROI(index);
                        
                        if CutFibersOutsideROIs
                            fiber = [fiber; ind'];
                            break;
                        end
                    end
                    %                     end
                    
                    if ((isnan( AnisoImage(ind(1),ind(2),ind(3)) )) ||  (AnisoImage(ind(1),ind(2),ind(3)) < anisoIndexStopThreshold) )
                        break;
                    end
                    
                    index1= sub2ind([M N N_Slice],ind(1), ind(2), ind(3));
                    index2=find(FA_Image_Mask==index1);
                    
                    if ( isempty(index2))
                        break;
                    end
                    
                    
                    pds_vox=pds(:,index2);
                    
                    
                    innerProd = pds_vox'* oldpds_vox;
                    
                    if (innerProd < 0)
                        angle = acos(-innerProd - EPSILON);
                        pds_vox = -pds_vox;
                    else
                        angle = acos(innerProd - EPSILON);
                    end
                    oldpds_vox = pds_vox;
                    
                end
                
                if (CutFibersOutsideROIs && ~passes)
                    fiber=[];
                end
                
                %------------------------------The other Direction
                angle = 0;
                currentPoint = firstPoint;
                pds_vox = -firstpds_vox;
                oldpds_vox = pds_vox;
                
                while ( sum(currentPoint < dwiSize)==3 &&  sum(currentPoint >=1)==3 && angle < maxFiberAngle) % && innerProd > 0.5   )  %
                    
                    if (CutFibersOutsideROIs && passes)
                        break;
                    end
                    
                    fiber = [currentPoint'; fiber ];
                    
                    tmpPoint= currentPoint;
                    tmpPoint = tmpPoint + pds_vox*DELTA_FACT;
                    while( sum(tmpPoint < round(currentPoint)+1)==3 && sum(tmpPoint > round(currentPoint)-1 )==3 )
                        tmpPoint = tmpPoint + pds_vox*DELTA_FACT;
                    end
                    currentPoint = tmpPoint;
                    
                    ind = round(currentPoint);
                    
                    try
                        index= sub2ind([M N N_Slice],ind(1), ind(2), ind(3));
                    catch
                        break;
                    end
                    
                    if ( isnan(index))
                        break;
                    end
                    if ( ( WholeBrainOutput || WholeBrainInput || OutputROI(index) ~= r1) && ( OutputROI(index) ~= 0) )%  (OutputROI(index)~=InputROI(index))  )
                        passes = true;
                        r2=OutputROI(index);
                        
                        if CutFibersOutsideROIs
                            fiber = [ind'; fiber ];
                            break;
                        end
                    end
                    
                    if ((isnan( AnisoImage(ind(1),ind(2),ind(3)) )) ||  (AnisoImage(ind(1),ind(2),ind(3)) < anisoIndexStopThreshold) )
                        break;
                    end
                    
                    index1= sub2ind([M N N_Slice],ind(1), ind(2), ind(3));
                    index2=find(FA_Image_Mask==index1);
                    
                    if ( isempty(index2))
                        break;
                    end
                    
                    pds_vox=pds(:,index2);
                    
                    innerProd = pds_vox'* oldpds_vox;
                    if (innerProd < 0)
                        angle = acos(-innerProd - EPSILON);
                        pds_vox = -pds_vox;
                    else
                       angle = acos(innerProd - EPSILON);
                    end
                    
                    oldpds_vox = pds_vox;
                    
                    %------------------------------------------------------
                end
                
                Fiber=[];
                FA_Fiber=[];
                
                if (passes)
                    
                    if (size(fiber,1) >= minFiberLength)
                        FibersCount = FibersCount+1;
                        for mm =1:size(fiber,1)
                            ind = round(fiber(mm,:));
                            try
                                index= sub2ind([M N N_Slice],ind(1), ind(2), ind(3));
%                                 if ind(2) < BranchStart
                                fiberImage(ind(1),ind(2),ind(3))= ceil(256*rand);

                                   
                                Fiber=[Fiber fiber(mm,:)];
                                FA_Fiber=[ FA_Fiber FA(ind(1),ind(2),ind(3))];
                            catch %#ok<CTCH>
                                fprintf('passes loop Index problem %s\n'); %#ok<CTPCT>
                                
                            end
                            
                        end
                        NumberOfFibers(r1,r2)=NumberOfFibers(r1,r2)+1;
                        n=NumberOfFibers(r1,r2);
                        FIBERS(r1,r2,n).Fiber=Fiber;
                        FIBERS(r1,r2,n).FA= FA_Fiber;
                    end
                end            
        end
        
    end
    
%     close(h_w),pause(0.01);
end

close(h_w),pause(0.01);

str=[];
str = [strFile,'Fibers.mat'];
[filename, pathname] = uiputfile(str, 'Save Fiber Structure as...');

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

save(targetfile, 'FIBERS', 'NumberOfFibers');


fiberImageSize=size(fiberImage);
N=size(fiberImage,2);
N_Slice=size(fiberImage,3);


str=[];
str = [strFile,'fibersImage.img'];
[filename, pathname] = uiputfile(str, 'Save fibersImage as...');

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

Analyze = Write_Analyze (imgtargetfile, fiberImage ,'int16', VDims);

% save ([strFile,'seg.mat'], 'fiberImage' );

str = [strFile,'outputFibers.txt'];
[filename, pathname] = uiputfile(str, 'Save Fibers as...');

exist_f = isa(filename,'char');
if exist_f
    targetfile = fullfile(pathname,filename);
else
    return;
end
figure;
imshow(NumberOfFibers,[])

file_1 = fopen(targetfile,'w+');
figure,
num = 0;
linelgh = [];
colors = jet(NubmerOfROILabels*NubmerOfROILabels);

for i = 1:NubmerOfROILabels
    for j = 1:NubmerOfROILabels
%         if NumberOfFibers(i,j) > 80   %40 MEG
        for n=1:NumberOfFibers(i,j)
            
            if isempty(FIBERS(i,j,n).Fiber) == 0
                
                
                i1=[]; j1=[]; k1=[]; Fi=[];
                for k=1:length(FIBERS(i,j,n).Fiber)/3
                    x=FIBERS(i,j,n).Fiber(1+3*(k-1));
                    y=FIBERS(i,j,n).Fiber(2+3*(k-1));
                    z=FIBERS(i,j,n).Fiber(3*k);
                    
                    X=[x; y; z];
                    
                    Fi=[Fi X'];
                    
                    try
                        % FiberImage(ceil(X(1)), ceil(X(2)),ceil(X(3)))=mod(32*i,255);  %floor(255*i/NubmerOfROILabels);  %300
                    catch
                    end
                    
                    i1=[ i1; X(1)];
                    j1=[ j1; X(2)];
                    k1=[ k1; X(3)];
                end
                FIBERS(i,j,n).Fiber=Fi;
                %                 FIBERS(i,j,n).Fiber=repmat(FiberVoxelSize,1,length(FIBERS(i,j,n).Fiber)/3).* FIBERS(i,j,n).Fiber;
                FIBERS(i,j,n).Fiber=repmat(VDims,1,length(FIBERS(i,j,n).Fiber)/3).* FIBERS(i,j,n).Fiber;
                fprintf(file_1,'%4.2f ',length(FIBERS(i,j,n).Fiber));
                fprintf(file_1,'%4.2f ',FIBERS(i,j,n).Fiber);
                linelgh = [linelgh;length(i1)];
                num = num +1;
                lpt(num) = line(i1,j1,k1);
                set(lpt(num),'Color',colors(sub2ind([NubmerOfROILabels NubmerOfROILabels],i, j),:));
                if num == 1
                    hold on
                end
                
            end
        end
        
%         end
    end
    
end
axis equal

fclose(file_1);





