

pds(3,:)=-pds(3,:);

R=[ 1 0 0;
    0 1 0;
    0 0 1 ];
T=[0 ;0 ;0 ];

IntThreshold=1.5;
% % % FiberName='Fx';%Cg';
conf='SuccessiveFiveLables' ; %'EveryOtherSuccessiveFiveLables'  ; %'EveryOtherSuccessiveThreeLables' ; %'SuccessiveFiveLables';  %'AllLables'; %'SuccessiveLables';%


% % % =160;%65;
OnlineTag=0;
DTI_Tag=1;


useOutROI=true;   %false;
LoadROIMatfile=1;  
CutFibersOutsideROIs=1; 

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


NubmerOfROILabels=max(ROI(:));   %for FX     %    /2;  for 2sided Cg
ROI1=ROI;

        
        
strFile=[FiberName,conf,'_NOL',num2str(NubmerOfROILabels),'_'];
switch conf
    case 'AllLables'
        ROIConfigMat=ones(NubmerOfROILabels);
    case 'SuccessiveLables'
        ROIConfigMat=zeros(NubmerOfROILabels);
    case 'SuccessiveFiveLables'; %  'EveryOtherSuccessiveFiveLables'  %'EveryOtherSuccessiveThreeLables'
        ROIConfigMat=zeros(NubmerOfROILabels);
        for i=1 : NubmerOfROILabels
            for j=i-4:1:i+4
                try
                    ROIConfigMat(i,j)=1;
                catch
                end
            end
        end
end
            
for i=1:NubmerOfROILabels
    if i == NubmerOfROILabels-1
        ROIConfigMat(i,i+1)=1;   %-1;
    elseif i < NubmerOfROILabels-1
        ROIConfigMat(i,i+1)=1;
    end
    if i == 2
        ROIConfigMat(i,i-1)=1;  %-1;
    elseif i > 2
        ROIConfigMat(i,i-1)=1;
    end
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);
h_w = waitbar(0,'Finding tracts ...');
for r1= 1:NubmerOfROILabels
    waitbar(r1/NubmerOfROILabels);
    for r2= 1:NubmerOfROILabels
        if ROIConfigMat(r1,r2) ~=0
            
            CutFibersOutsideROIs=ROIConfigMat(r1,r2);
            
            InputROI=zeros(size(ROI1));
            InputROI(ROI1==r1)=ROI1(ROI1==r1);
            
            InputROI_Image_Mask=find(InputROI > 0 );
            n_InputROI_Vox = length(InputROI_Image_Mask);
            
            OutputROI=zeros(size(ROI1));
            OutputROI(ROI1==r2)=ROI1(ROI1==r2);
            
            num=0;
            fiber = [];
            indexInputROI=0;
                          
            while (indexInputROI < n_InputROI_Vox)
                
                indexInputROI=indexInputROI+1;
                 
%                 waitbar(indexInputROI/n_InputROI_Vox);
                [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
                    
                    passes = ~useOutROI;
                    pds_vox=pds(:,index2First);
                    fiber = [];
                    angle = 0;
                    currentPoint=firstPoint;
                    firstpds_vox=pds_vox;
                    oldpds_vox = pds_vox;
%                     fiberFound=0;
                    %----------------------------------------------------------
                    while ( sum(currentPoint < dwiSize)==3 &&  sum(currentPoint >=1)==3  && angle < maxFiberAngle)
                        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) ||~isreal(index))
                            break;
                        end
                        
                        if ( useOutROI && ( OutputROI(index) ~= 0) &&  (OutputROI(index)~=InputROI(index)) && abs(Label(index)-Label(firstPoint(1), firstPoint(2), firstPoint(3) )) < 4 )
                            lastpoint=ind;
                            passes = true;
                            %                             r2=OutputROI(index);
                            fiberFound=1;
                            if CutFibersOutsideROIs==1
                                fiber = [fiber; ind'];
                                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=reshape(pds(:,index2),3,1);
                        
                        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 ~passes    %isempty(fiber) %
                        
                         fiber=[];
                        %------------------------------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)
                            
                            if (CutFibersOutsideROIs==1 && 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) || ~isreal(index))
                                break;
                            end
                            if ( useOutROI && ( OutputROI(index) ~= 0) &&  (OutputROI(index)~=InputROI(index)) && abs(Label(index)-Label(firstPoint(1), firstPoint(2), firstPoint(3) )) < 4 )
                                lastpoint=ind;
                                passes = true;
                                %                             r2=OutputROI(index);
                                
                                if CutFibersOutsideROIs==1
                                    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=reshape(pds(:,index2),3,1);
                            
                            
                            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
                    end
                    
                    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));

                                    
                                    dist_to_Right_Posterior = min(sum((Act_Right_Posterior-repmat(ind' , 1, size(Act_Right_Posterior,2))).^2,1));                     
                                    dist_to_Right_Superior = min(sum((Act_Right_Superior-repmat(ind' , 1, size(Act_Right_Superior,2))).^2,1));
                                    dist_to_Right_Anterior = min(sum((Act_Right_Anterior -repmat(ind' , 1, size(Act_Right_Anterior,2))).^2,1));                                 
                         
                                    dist_to_Left_Posterior = min(sum((Act_Left_Posterior-repmat(ind' , 1, size(Act_Left_Posterior,2))).^2,1));                     
                                    dist_to_Left_Superior = min(sum((Act_Left_Superior-repmat(ind' , 1, size(Act_Left_Superior,2))).^2,1));
                                    dist_to_Left_Anterior = min(sum((Act_Left_Anterior-repmat(ind' , 1, size(Act_Left_Anterior ,2))).^2,1));                                 
                                    
                                    [mmm La]=min([dist_to_Right_Posterior dist_to_Right_Superior dist_to_Right_Anterior dist_to_Left_Posterior dist_to_Left_Superior dist_to_Left_Anterior]);
    
                                    fiberImage(ind(1),ind(2),ind(3))=32*La ;
                                    
                                    Fiber=[Fiber fiber(mm,:)];
                                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;
                        end
                    end
                    
                end
                
            end
            
            
        end
    end
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
cd (pathname);
save(targetfile, 'FIBERS', 'NumberOfFibers');

N=size(fiberImage,2);
N_Slice=size(fiberImage,3);
fiberImageSize=size(fiberImage);
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 = UT_Write_Analyze (imgtargetfile, fiberImage ,'int16', VDims);

data=fiberImage;
save([strFile 'FT_Seg'],'data','Act_Right_Posterior', 'Act_Right_Superior' , 'Act_Right_Anterior','Act_Left_Posterior', 'Act_Left_Superior' , 'Act_Left_Anterior' );
cd (pathname);

str=[];
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
file_1 = fopen(targetfile,'w+');
figure,
num = 0;
linelgh = [];
colors = jet(NubmerOfROILabels*NubmerOfROILabels);

for i = 1:NubmerOfROILabels
    for j = 1:NubmerOfROILabels
        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=R*[x; y; z]+T;
                    Fi=[Fi X'];
                    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(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
axis equal
fclose(file_1);


