close all
clear all

global qvalues;

addpath f:\MatlabCodes

cd 'Z:\Epilepsy_ANALYZE\EP1218\Dicoms\EP1218-001\DTI 2.6mm'

%%------------------------Parameters-----------------------------

bvalues=1000;  %0.5;     %sec/mm^2
Delta=0.0334; %sec
delta=0.0322; %sec

qvalues=sqrt(bvalues./(4*pi^2*(Delta-delta/3)))*1e0;   %mm^2/sec

NrOfSamples=1;

NrOfDirections =25; % define number of directions and gradient set (add/change your settings)
NrOfB0=1; 
resampling=1;
smoothing=0;

NVolumes = NrOfDirections + NrOfB0; % +1;  %Philips

%___________________________________________________________

if smoothing
    SmoothingTag='Smo';
else
    SmoothingTag='NSmo';
end

if resampling
    resamplingTag='reSmp';
else
    resamplingTag='NreSmp';
end

current_dir=pwd;

[Input_directory_name] = uigetdir(pwd,'Select directory of Input Dicom Images ...');
exist_d = isa(Input_directory_name,'char');
if ~exist_d
    return;
end

cd (Input_directory_name');

mkdir ('results');

[Output_directory_name] = uigetdir(pwd,'Select directory of Output Resampled Imgs ...');
exist_d = isa(Output_directory_name,'char');
if ~exist_d
    return;
end


cd(Input_directory_name);
dir_struc = dir;
sd = size(dir_struc,1);  

f_name = [Input_directory_name '/' dir_struc(10).name];

d = dicominfo(f_name);

VDims(1)=d.PixelSpacing(1);
VDims(2)=d.PixelSpacing(2);
VDims(3)=d.SpacingBetweenSlices;
MDims = [d.Width d.Height (sd-2)/NVolumes];

h_w = waitbar(0,'Analysing dicom data...');
cn = 0;
for i=1:sd
    f_name = [Input_directory_name '/' dir_struc(i).name];
    try
        d = dicominfo(f_name);
        cn = cn + 1;
        dicom_info{cn} = d;
        waitbar(cn/sd)
    catch
    end
end
close(h_w),pause(0.01);

DWI = zeros([MDims NVolumes],'uint16');%repmat(uint16(0),[MDims NVolumes]);

h_w = waitbar(0,'Reading dicom data...');

z=double(MDims(3));

mi=1000000;
for i=1:length(dicom_info)
    dummy = dicomread(dicom_info{i});
    vv(i) = dicom_info{i}.InstanceNumber;
    if vv(i) < mi
        mi=vv(i);
    end
end
[uu ii]=sort(vv);
for i=1:length(dicom_info)
    dummy = dicomread(dicom_info{i});
    vv = dicom_info{i}.InstanceNumber;    %ii(i); 
    p = mod(vv,z);
    dd = double(floor(vv/z))+1;
    if p==0
        p=MDims(3);
        dd=dd-1;
    end
    
%     disp([p dd]);
       DWI(:,:,p,dd) = dummy';
%     DWI(:,:,MDims(3)+1-p,dd) = dummy';
    waitbar(i/length(dicom_info));
end
close(h_w),pause(0.01);

% N_DWI=zeros(size(DWI));
% 
% for p=1:MDims(3)
%      N_DWI(:,:,p,:)=DWI(:,:,MDims(3)-p+1,:);
% end
% 
% DWI=N_DWI;

%--------For more than 1 B0 images

% % dummy1 = repmat(uint16(0),MDims);
% dummy1=mean(DWI(:,:,:,1:NrOfB0),4);
% DWI(:,:,:,1)=dummy1;
% DWI(:,:,:,2:NrOfB0)=[];
% NVolumes=size(DWI,4);

if MDims(1) > 192
    From_256x256_To_128x128_Tag=1;
else
    From_256x256_To_128x128_Tag=0;
end


if resampling
 
    DWI(DWI==0)=1;

    h = waitbar(0,'Resampling...');
    M1 = double(MDims(1));
    N1 = double(MDims(2));
    N_Slice1= double(MDims(3));  %[256 256 52]
   
     VDims1(3)=VDims(3);
    if From_256x256_To_128x128_Tag    %Philips 1.5 T and GE
        
        VDims2(1)=VDims(1)*2;
        VDims2(2)=VDims(2)*2;
       
        M2=length(1:2:M1); 
        N2=length(1:2:N1);
        
        VolumeThickness=(N_Slice1-1)*VDims1(3);
        N_Slice2=round(VolumeThickness/VDims2(1))+1;
        VDims2(3)=VolumeThickness/(N_Slice2-1);
        [xi,yi,zi] = meshgrid(1:2:M1,1:2:N1,1:VDims2(3)/VDims1(3):N_Slice1);%     1:2*(M1-1)/M1:M1
        
    else
        
        VDims2(1)=VDims(1) ;
        VDims2(2)=VDims(2);
        
        M2=length(1:1:M1);
        N2=length(1:1:N1);
        
        VolumeThickness=(N_Slice1-1)*VDims1(3);
        N_Slice2=round(VolumeThickness/VDims2(1))+1;
        VDims2(3)=VolumeThickness/(N_Slice2-1);
        [xi,yi,zi] = meshgrid(1:1:M1,1:1:N1,1:VDims2(3)/VDims1(3):N_Slice1);
        
    end
    
    im_volume2=zeros(size(xi));
    
    ResampledDWI = repmat(double(0),[M2 N2 N_Slice2 NVolumes]);
    
    pause(0.05)
    
    for i=1:NVolumes
        
        strFile=[Output_directory_name,'\' resamplingTag,'_',SmoothingTag,'_IMG_',num2str(i)];
        waitbar(i/NVolumes);
        im_volume1=single(DWI(:,:,:,i));
        im_volume2= interp3(im_volume1,xi,yi,zi,'cubic');
        
        if smoothing > 0
%             im_volume2 = anisodiff3D(im_volume2, num_iter, delta_t, kappa ,option, VDims);
          im_volume2=smooth3(im_volume2,'gaussian',[3 3 3],smoothing);
        end
        
        ResampledDWI(:,:,:,i) = (im_volume2); %uint16
        
        Analyze = Write_Analyze (strFile, ResampledDWI(:,:,:,i) ,'float', VDims2);
    end
    
    close(h);pause(0.05)
    clear DWI;
    DWI=ResampledDWI;
    % DWI(:,:,:,end)=[];
    clear ResampledDWI;
    MDims = [size(DWI,1) size(DWI,2) size(DWI,3)];
    VDims=VDims2;
else
    
    h = waitbar(0,'Saving...');

    for i=1:NVolumes
        strFile=[Output_directory_name,'/' resamplingTag,'_',SmoothingTag,'_IMG_',num2str(i)];
        waitbar(i/NVolumes);

        Analyze = Write_Analyze (strFile, DWI(:,:,:,i) ,'int16', VDims);
    end
    close(h),pause(0.01);

end
%     cd (current_dir)
clear dicom_info
clear im_volume1 im_volume2 dummy dummy1 uu vv xi yi zi

h = waitbar(0,'Memory allocation...');

mask = sum(DWI(:,:,:,2:end),4);  % Use the smallest redius shell for creating a mask

mask = mat2gray(mask(:));
level = graythresh(mask);
mask = im2bw(mask,level);
mask = reshape(mask,MDims(1),MDims(2),MDims(3));
close(h);pause(0.05)

h = waitbar(0,'Masking ..');

se = strel('disk',6);

for i=1:size(mask,3)
    mask(:,:,i) = bwareaopen(mask(:,:,i),500);
    mask(:,:,i) = imdilate(mask(:,:,i),se);
    mask(:,:,i) = imfill(mask(:,:,i),'holes');
    mask(:,:,i) = imerode(mask(:,:,i),se);
    waitbar(i/size(mask,3));
end
mask = logical(mask);

% mask=ones(MDims);
Vmask=find(mask);
close(h);pause(0.05)

%     figure; imshow(mask(:,:,50),[]);

NumberOfVoxels=length(Vmask);

g=1*GradientTables(NrOfDirections); 
% g(g(:,3)<0 ,:)= -g(g(:,3)<0 ,:);

% tmp=g(:,1);
% g(:,1)=g(:,2);
% g(:,2)=tmp;

%  g(:,[1 2 3]) = g(:,[2 1 3]);
%  g(:,3)=-1*g(:,3);

h = waitbar(0,'Memory allocation...');

if ~exist('b')
    b = bvalues*[g(:,1).^2 2*g(:,1).*g(:,2) 2*g(:,1).*g(:,3) g(:,2).^2 2*g(:,2).*g(:,3) g(:,3).^2];
else
    % Changing to correct Matlab 'x', 'y', 'z' coordinate system
    b(:,[1 2 3 4 5 6]) = b(:,[2 4 6 1 5 3]);
    % factor 2 is for cross terms
    b(:,[2 3 5]) = 2*b(:,[2 3 5]);
    % change sign of z-component (choice of origin for ExploreDTI)
    b(:,[3 5]) = -b(:,[3 5]);
end

DWI(DWI==0)=1;
% DT = repmat(single(0), [MDims(1) MDims(2) MDims(3) 9]);
DT = zeros(MDims(1), MDims(2), MDims(3), 9);

order = [1 2 3 5 6 9];
close(h);pause(0.05)

h = waitbar(0,'Calculating DT...');

for i=1:MDims(1)
    for j=1:MDims(2)
        for k=1:MDims(3)
            
            if mask(i,j,k)
                
                A = log(double(DWI(i,j,k,1))./double(squeeze(DWI(i,j,k,2:end))));
                
                X = b\A;
%                       [X, chi_sq, flag, O_M] = RESTORE_func(double(squeeze(DWI(i,j,k,:))),b);
                
                DT(i,j,k,order) = single(X);
                
            end
            
        end
    end
    waitbar(i/MDims(1));
end

close(h);pause(0.05)


DT(:,:,:,4) = DT(:,:,:,2);
DT(:,:,:,7) = DT(:,:,:,3);
DT(:,:,:,8) = DT(:,:,:,6);

DT(isnan(DT))=0;

[FEFA, FA, FE, SE, eigval,Cl,Cp,Cs,MD] = View3DDTMaps(DT,mask);

FA=FA./sqrt(3);



Lambda1 = eigval(:,:,:,1);
Lambda2 = eigval(:,:,:,2);
Lambda3 = eigval(:,:,:,3);

MD=(Lambda1+Lambda2+Lambda3)/3;
AD=Lambda1;
RD=(Lambda2+Lambda3)/2;

magFactor=100; 
cd (Output_directory_name);

Analyze = Write_Analyze ('Lambda1',magFactor * Lambda1 ,'float', VDims);
Analyze = Write_Analyze ('Lambda2',magFactor * Lambda2 ,'float', VDims);
Analyze = Write_Analyze ('Lambda3',magFactor * Lambda3 ,'float', VDims);
Analyze = Write_Analyze ('AD',magFactor * AD ,'float', VDims);
Analyze = Write_Analyze ('RD',magFactor * RD ,'float', VDims);

Analyze = Write_Analyze ('Cp', Cp ,'float', VDims);
Analyze = Write_Analyze ('Cs', Cs ,'float', VDims);
Analyze = Write_Analyze ('Cl', Cl ,'float', VDims);
Analyze = Write_Analyze ('FA',FA ,'float', VDims);
Analyze = Write_Analyze ('ADC',magFactor * MD ,'float', VDims);
ADC_Map=zeros(MDims(1),MDims(2),MDims(3));

for kk=1:MDims(2)
    
    ADC_Map(:,MDims(2)+1-kk,:) = MD(:,kk,:);
end

Analyze = Write_Analyze ('ADC_Map',magFactor * ADC_Map ,'float', VDims);

figure; imshow(reshape(FEFA(:,:,30,:),MDims(1),MDims(2),3),[]);

fnam=[Output_directory_name,'/DTI_',resamplingTag,'_',SmoothingTag,'.mat'];

save(fnam,'DT','DWI','FE','FA','MD','AD','RD','FEFA','VDims','MDims','eigval','g','SE','Cl','Cp','Cs','Lambda1','Lambda2','Lambda3');

clc

