%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%TENSORCALC2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Written by Enzyme. Feb, 2004%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Matlab implementation of the instructions given by Dr. Maj Hedehus
in: http://www-radiology.stanford.edu/majh/
This program calculates the tensor for an MRI slice and creates maps for
Fractional Anisotropy (FAmap), traceADC (tADC), principal eigenvector
(lambda1) and its orthogonals (lambda2 and lambda3), and an FA-weighted
colormap for the principal eigenvector (cm).
It also creates the tensor elements (xx,yy,zz,xy,yz,xz).
HOW TO USE:
1. Load the image slices to variables bo (note it is bo, letter 'o'),
b1, b2, b3, b4, b5 and b6. Use the function dicomread to do this i.e:
>>bo=dicomread('image.dcm');
2. Run this script
>>DTI
3. The maps and images are created with the names as described above. To
view a particular map or image type
>>imagesc(FAmap);colormap(gray);axis image;
Please note this program requires the image processing toolbox with dicomread
function. If you do not have it, you can load raw data images with the 'load'
function, availabe in any matlab.
VERY IMPORTANT!!!! This program assumes you only have one b-value and it is
set to b=1000.
Also: feel free to make any modifications to this code.
ENJOY!0001 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0002 %%%%%%%%%%%%%%%%%%%%%%TENSORCALC2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0004 %Written by Enzyme. Feb, 2004%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0006 %Matlab implementation of the instructions given by Dr. Maj Hedehus 0007 %in: http://www-radiology.stanford.edu/majh/ 0008 %This program calculates the tensor for an MRI slice and creates maps for 0009 %Fractional Anisotropy (FAmap), traceADC (tADC), principal eigenvector 0010 %(lambda1) and its orthogonals (lambda2 and lambda3), and an FA-weighted 0011 %colormap for the principal eigenvector (cm). 0012 %It also creates the tensor elements (xx,yy,zz,xy,yz,xz). 0013 % 0014 % 0015 %HOW TO USE: 0016 % 0017 %1. Load the image slices to variables bo (note it is bo, letter 'o'), 0018 % b1, b2, b3, b4, b5 and b6. Use the function dicomread to do this i.e: 0019 % >>bo=dicomread('image.dcm'); 0020 %2. Run this script 0021 % >>DTI 0022 %3. The maps and images are created with the names as described above. To 0023 % view a particular map or image type 0024 % >>imagesc(FAmap);colormap(gray);axis image; 0025 % 0026 % 0027 % 0028 %Please note this program requires the image processing toolbox with dicomread 0029 %function. If you do not have it, you can load raw data images with the 'load' 0030 %function, availabe in any matlab. 0031 % 0032 % 0033 %VERY IMPORTANT!!!! This program assumes you only have one b-value and it is 0034 % set to b=1000. 0035 %Also: feel free to make any modifications to this code. 0036 % 0037 %ENJOY! 0038 0039 0040 0041 0042 0043 0044 ADC101=log(im2double(b1)./im2double(bo))./(-1000); %Calculate ADC maps for each direction 0045 ADC_101=log(im2double(b2)./im2double(bo))./(-1000); 0046 ADC011=log(im2double(b3)./im2double(bo))./(-1000); 0047 ADC01_1=log(im2double(b4)./im2double(bo))./(-1000); 0048 ADC110=log(im2double(b5)./im2double(bo))./(-1000); 0049 ADC_110=log(im2double(b6)./im2double(bo))./(-1000); 0050 0051 transform=[1 1 0 2 0 0;1 0 1 0 2 0;0 1 1 0 0 2;1 1 0 -2 0 0;1 0 1 0 -2 0;0 1 1 0 0 -2] %the transformation matrix to obtain tensor elements 0052 0053 xx=zeros(256); %Create empty tensor elements 0054 yy=zeros(256); 0055 zz=zeros(256); 0056 xy=zeros(256); 0057 xz=zeros(256); 0058 yz=zeros(256); 0059 lambda1=zeros(256); %create empty images 0060 lambda2=zeros(256); 0061 lambda3=zeros(256); 0062 0063 r=zeros(256); 0064 g=zeros(256); 0065 b=zeros(256); 0066 FAmap=zeros(256); 0067 tADC=zeros(256); 0068 0069 for i=1:256 %Fill in tensor elements pixel by pixel. 0070 for j=1:256 0071 if (bo(i,j)>60) %noise threshold. Change this to zero if you want the whole image to be calculated (takes longer) 0072 ADCm=[ADC110(i,j);ADC101(i,j);ADC011(i,j);ADC_110(i,j);ADC_101(i,j);ADC01_1(i,j)]; %see Hedehus, http://www-radiology.stanford.edu/majh/ 0073 ADCe=inv(transform)*ADCm; 0074 xx(i,j)=ADCe(1,1); 0075 yy(i,j)=ADCe(2,1); 0076 zz(i,j)=ADCe(3,1); 0077 xy(i,j)=ADCe(4,1); 0078 xz(i,j)=ADCe(5,1); 0079 yz(i,j)=ADCe(6,1); 0080 ten=[xx(i,j) xy(i,j) xz(i,j);xy(i,j) yy(i,j) yz(i,j);xz(i,j) yz(i,j) zz(i,j)]; %the tensor itself 0081 [V,D]=eig(ten); 0082 D=eig(ten); %get eigenvalues 0083 D=abs(D); 0084 D=sort(D); %sort the eigenvalues (upwards) 0085 e1=D(3,1);e1=e1*2; 0086 e2=D(2,1);e2=e2*2; 0087 e3=D(1,1);e3=e3*2; 0088 lambda1(i,j)=e1; 0089 lambda2(i,j)=e2; 0090 lambda3(i,j)=e3; 0091 V=abs(V); %get the x,y,z components of e3 0092 r(i,j)=V(1,3); 0093 g(i,j)=V(2,3); 0094 b(i,j)=V(3,3); 0095 trace=(e1+e2+e3)/3; 0096 %FA=(sqrt(3/2)*sqrt((1/3)*((e1-e2)^2+(e2-e3)^2+(e3-e1)^2)))/(sqrt(e1^2+e2^2+e3^3)); %Another formula for FA 0097 FA=(sqrt(3*((e1-trace)^2+(e2-trace)^2+(e3-trace)^2)))/(sqrt(2*(e1^2+e2^2+e3^2))); %Le Bihan 2001 0098 tADC(i,j)=trace; 0099 FAmap(i,j)=FA; 0100 end; 0101 end; 0102 end; 0103 0104 0105 0106 0107 rgb=cat(3,r,g,b); 0108 FAmap3=cat(3,FAmap,FAmap,FAmap); %needed to combine with the colormap, need 3D. 0109 cm=rgb.*FAmap3; %FA weighting for colormap 0110 figure, image (cm);axis image 0111 0112 0113 0114 0115 0116 0117 0118