clear all;
close all;
img = zeros(128,128);

[fileima,pname] = uigetfile('\*.ima','Please input Image filename');
dashpoint=findstr(fileima,'-');
filehead=fileima(1:dashpoint(2));
firstnumber=str2num(fileima(dashpoint(2)+1:findstr(fileima,'.')-1));
tempstr=input('Please input Slice number?');
nslice=tempstr;
dwith=5; % DWI threshold

qwrite=input('Write DTI images(y/n)?','s');


if qwrite=='y' | qwrite=='Y'
    [outfilehead,pfile] = uiputfile('\*.*','Please input Image filename');    
    %qinter=input('Write images with interpolation(y/n)?','s');
    qinter='y';
    
end
disp('Loading DTI images....');
for cslice=1:nslice
    disp(['Loading DTI images....Slice' num2str(cslice)]);
    %load B0
    fid=fopen(strcat(pname,filehead,num2str(firstnumber+(cslice-1)),'.ima'),'r','b');
    %fid2=fopen(strcat(pname,'342-9-',num2str(firstnumber+(cslice-1)+7*16),'.ima'),'r','b');
    fseek(fid,6144,-1);
    %fseek(fid2,6144,-1);
    dtiall(cslice,1,:,:)=fread(fid,[128 128],'uint16')';
    % dtiall(cslice,1,:,:)=(fread(fid,[128 128],'uint16')'+fread(fid2,[128 128],'uint16')')/2;
    fclose(fid);
    for gd = 1:6;
        %load 6 directions
        fid=fopen(strcat(pname,filehead,num2str(firstnumber+nslice+(cslice-1)*6+gd-1),'.ima'),'r','b');
        % fid2=fopen(strcat(pname,'342-9-',num2str(firstnumber+nslice+(cslice-1)*6+gd-1+7*16),'.ima'),'r','b');
        fseek(fid,6144,-1);
        % fseek(fid2,6144,-1);
        dtiall(cslice,gd+1,:,:)=fread(fid,[128 128],'uint16')';
        %dtiall(cslice,gd+1,:,:)=(fread(fid,[128 128],'uint16')'+fread(fid2,[128 128],'uint16')')/2;
        fclose(fid);
    end
end

disp('Processing DTI images....');
%Thresholding DWIs
dtiall(find(dtiall<dwith))=1;
%B Matrix
M = [1, 0, 1, 0, 0, 2;...
        1, 0, 1, 0, 0,-2;...
        0, 1, 1, 0, 2, 0;...
        0, 1, 1, 0,-2, 0;...
        1, 1, 0, 2, 0, 0;...
        1, 1, 0,-2, 0, 0];
Diffusion_tensor = zeros(128,128);
Diffusion_tensor1 = zeros(128,128);
Diffusion_tensor2 = zeros(128,128);
Diffusion_tensor3 = zeros(128,128);
ori_ani=zeros(128,128);
ori_ani2 = zeros(128,128);
ori=zeros(128,128);
ori2=zeros(128,128);
fa=zeros(128,128);
ra=zeros(128,128);
vr=zeros(128,128);
for cslice=1:nslice
    disp(['Processing DTI images....Slice' num2str(cslice)]);
    ADCimg1=log(reshape(dtiall(cslice,2,:,:),128,128)./reshape(dtiall(cslice,1,:,:),128,128))/-1000;
    ADCimg2=log(reshape(dtiall(cslice,3,:,:),128,128)./reshape(dtiall(cslice,1,:,:),128,128))/-1000;
    ADCimg3=log(reshape(dtiall(cslice,4,:,:),128,128)./reshape(dtiall(cslice,1,:,:),128,128))/-1000;
    ADCimg4=log(reshape(dtiall(cslice,5,:,:),128,128)./reshape(dtiall(cslice,1,:,:),128,128))/-1000;
    ADCimg5=log(reshape(dtiall(cslice,6,:,:),128,128)./reshape(dtiall(cslice,1,:,:),128,128))/-1000;
    ADCimg6=log(reshape(dtiall(cslice,7,:,:),128,128)./reshape(dtiall(cslice,1,:,:),128,128))/-1000;
    
    
    for i = 1:128;
        for j = 1:128;
            if dtiall(cslice,1,i,j) > 100
                ADC = [ADCimg1(i,j);...
                        ADCimg2(i,j);...
                        ADCimg3(i,j);...
                        ADCimg4(i,j);...
                        ADCimg5(i,j);...
                        ADCimg6(i,j)];
                
                rADC = inv(M)*ADC;
                
                diffusion = [rADC(1,1), rADC(4,1), rADC(6,1);...
                        rADC(4,1), rADC(2,1), rADC(5,1);...
                        rADC(6,1), rADC(5,1), rADC(3,1)];
                
                [Evector,Evalue] = eig(diffusion,'nobalance');
                
                meantrace=trace(abs(Evalue))/3;
                Diffusion_tensor1(i,j) = (abs(Evalue(1,1))-meantrace).^2;
                Diffusion_tensor2(i,j) = (abs(Evalue(2,2))-meantrace).^2;
                Diffusion_tensor3(i,j) = (abs(Evalue(3,3))-meantrace).^2;
                Diffusion_tensor(i,j) = Evalue(1,1).^2+Evalue(2,2).^2+Evalue(3,3).^2;
                
                fa(i,j)=(3^0.5)/(2^0.5)*(((Diffusion_tensor1(i,j)+Diffusion_tensor2(i,j)+Diffusion_tensor3(i,j))./Diffusion_tensor(i,j))^0.5);
                faimage(i,j,1:3)=[fa(i,j) fa(i,j) fa(i,j)];
                tracetemp=(rADC(1,1)+rADC(2,1)+rADC(3,1))*100;
                traceimage(i,j,1:3)=[tracetemp tracetemp tracetemp];
                
                ra(i,j)=(((Diffusion_tensor1(i,j)+Diffusion_tensor2(i,j)+Diffusion_tensor3(i,j))./(sqrt(3)*meantrace)));
                
                vr(i,j)=Evalue(1,1)*Evalue(2,2)*Evalue(3,3)/meantrace^3;
                
                V2(:,i,j) = Evector(1:3,2);
                % [temp ori(i,j)]=max(abs(Evector(:,3)),[],1);
                ei=abs(Evector(1:3,3)*Evalue(3,3));
                ori2(i,j,1:3)=(abs(ei)./sqrt(sum(ei.^2)));
                ori(i,j,1:3)=abs(Evector(1:3,3))./sqrt(sum(abs(Evector(1:3,3)).^2));
                
                ori_ani(i,j,1:3)=fa(i,j)*ori(i,j,1:3);
                ori_ani2(i,j,1:3)=fa(i,j)*ori2(i,j,1:3);
                
            else          
                Diffusion_tensor1(i,j) = 0;
                Diffusion_tensor2(i,j) = 0;
                Diffusion_tensor3(i,j) = 0;
                Diffusion_tensor(i,j) = 1;
                ori_ani(i,j,1:3)=[0 0 0];
                ori_ani2(i,j,1:3)=[0 0 0];
                ori(i,j,1:3)=[0 0 0];
                faimage(i,j,1:3)=[0 0 0];
                traceimage(i,j,1:3)=[0 0 0];
            end
        end
    end
    % figure(1000+cslice);imagesc(faimage),colormap(gray);axis off;axis equal;title('FA');
    figure(4000);
    ori_ani2=ori_ani2*1.5;
    ori_ani2(find(ori_ani2>1))=1;
    traceimage(find(traceimage>1))=1;
    traceimage(find(traceimage<0))=0;
    %image(ori_ani2);axis off;axis equal;
    if cslice==1
        tempmontage=[traceimage faimage ori_ani2];
        imgmontage=tempmontage;
    elseif cslice==2 | cslice==3
        tempmontage=[tempmontage traceimage faimage ori_ani2];
        imgmontage=tempmontage;
    elseif mod(cslice,3)==1 
        rowtempmontage=[traceimage faimage ori_ani2];
       imgmontage=[tempmontage;rowtempmontage zeros(128,128*6,3)];
    elseif mod(cslice,3)==2
        rowtempmontage=[rowtempmontage traceimage faimage ori_ani2];
     imgmontage=[tempmontage;rowtempmontage zeros(128,128*3,3)];
    else
        rowtempmontage=[rowtempmontage traceimage faimage ori_ani2];
        imgmontage=[tempmontage;rowtempmontage];
        tempmontage=imgmontage;
    end
    imshow(imgmontage);
    
    if qwrite=='y' | qwrite=='Y'
        % imwrite([traceimage faimage ori_ani2],strcat(pfile,outfilehead,sprintf('%03d',cslice),'.jpg'),'jpg');
        if qinter=='y' | qwrite=='Y'
            imwrite([imresize(traceimage,[256 256],'bilinear') imresize(faimage,[256 256],'bilinear') imresize(ori_ani2,[256 256],'bilinear')],strcat(pfile,outfilehead,sprintf('%03d',cslice),'.jpg'),'jpg','Quality',95);
        else
            imwrite([imresize(traceimage,[256 256]) imresize(faimage,[256 256]) imresize(ori_ani2,[256 256])],strcat(pfile,outfilehead,sprintf('%03d',cslice),'.jpg'),'jpg','Quality',95);
        end
    end
end