DTIguicode

0001 tic
0002 
0003 
0004 
0005 getDicomsgui           %uses my script to load an entire Dicom directory.
0006 
0007 disp('Calculating tensor...')
0008 
0009 noise=evalin('base','noise');                                   %noise threshold
0010 doFAmap=evalin('base','doFAmap');                               %write DICOMs?
0011 dotADC=evalin('base','dotADC');
0012 dolambda1=evalin('base','dolambda1');
0013 dolambda2=evalin('base','dolambda2');
0014 dolambda3=evalin('base','dolambda3');
0015 docm=evalin('base','docm');
0016 counter=evalin('base','counter'); %counter comes from getDicomsgui
0017 extension=evalin('base','extension');
0018 
0019 
0020 slices=counter/7;
0021 
0022 
0023 FAmaps=[];
0024 tADCmaps=[];
0025 lambda1_maps=[];
0026 lambda2_maps=[];
0027 lambda3_maps=[];
0028 cmaps=[];
0029 xyzmaps=[];
0030 
0031 
0032 
0033 
0034 
0035 
0036 
0037 warning off MATLAB:divideByZero                                 %supress matlab warning
0038 
0039 
0040 for numslice=1:slices
0041 
0042 bo=tensorvol(:,:,numslice,1);%bo=bo';                            %transposes the images to correct for getDicoms
0043 b1=tensorvol(:,:,numslice,2);%b1=b1';
0044 b2=tensorvol(:,:,numslice,3);%b2=b2';
0045 b3=tensorvol(:,:,numslice,4);%b3=b3';
0046 b4=tensorvol(:,:,numslice,5);%b4=b4';
0047 b5=tensorvol(:,:,numslice,6);%b5=b5';
0048 b6=tensorvol(:,:,numslice,7);%b6=b6';
0049 
0050 
0051 ADC101=log(im2double(b1)./im2double(bo))./(-1000);              %Calculate ADC maps for each direction
0052 ADC_101=log(im2double(b2)./im2double(bo))./(-1000);
0053 ADC011=log(im2double(b3)./im2double(bo))./(-1000);
0054 ADC01_1=log(im2double(b4)./im2double(bo))./(-1000);
0055 ADC110=log(im2double(b5)./im2double(bo))./(-1000);
0056 ADC_110=log(im2double(b6)./im2double(bo))./(-1000);
0057 
0058 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
0059 
0060 xx=zeros(256);                                                  %Create empty tensor elements
0061 yy=zeros(256);
0062 zz=zeros(256);
0063 xy=zeros(256);
0064 xz=zeros(256);
0065 yz=zeros(256);
0066 lambda1=zeros(256);                                             %create empty images
0067 lambda2=zeros(256);
0068 lambda3=zeros(256);
0069 r=zeros(256);
0070 g=zeros(256);
0071 b=zeros(256);
0072 FAmap=zeros(256);
0073 tADC=zeros(256);
0074 
0075 for i=1:256                                                     %Fill in tensor elements pixel by pixel.
0076     for j=1:256
0077         if (bo(i,j)>noise)                                      %noise threshold.  Change this to zero if you want the whole image to be calculated (takes longer)
0078             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/
0079             ADCe=inv(transform)*ADCm;
0080             xx(i,j)=ADCe(1,1);                                  %get tensor elements
0081             yy(i,j)=ADCe(2,1);
0082             zz(i,j)=ADCe(3,1);
0083             xy(i,j)=ADCe(4,1);
0084             xz(i,j)=ADCe(5,1);
0085             yz(i,j)=ADCe(6,1);
0086             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
0087             [V,D]=eig(ten);
0088             D=eig(ten);                                         %get eigenvalues
0089             D=abs(D);                                           
0090             E=D;                                                %get the eigenvalues before re-ordering
0091             D=sort(D);                                          %sort the eigenvalues (upwards)
0092             e1=D(3,1);e1=e1*2;
0093             e2=D(2,1);e2=e2*2;
0094             e3=D(1,1);e3=e3*2;
0095             lambda1(i,j)=e1;
0096             lambda2(i,j)=e2;
0097             lambda3(i,j)=e3;
0098             %V=abs(V);       %get the x,y,z components of e3
0099             largest=find(E==D(3,1));                            %find the value of the largest eigenvalue in the original matrix
0100             sz=length(largest);
0101             r(i,j)=V(1,largest(sz));                                %in order to get its corresponding eigenvector components in x, y and z
0102             g(i,j)=V(2,largest(sz));
0103             b(i,j)=V(3,largest(sz));
0104             trace=(e1+e2+e3)/3;
0105             %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
0106             FA=(sqrt(3*((e1-trace)^2+(e2-trace)^2+(e3-trace)^2)))/(sqrt(2*(e1^2+e2^2+e3^2)));  %Le Bihan 2001
0107             tADC(i,j)=trace;
0108             FAmap(i,j)=FA;
0109         end;
0110     end;
0111 end;
0112 
0113 
0114 
0115 
0116 rgb=cat(3,r,g,b);
0117 xyz=cat(3,r,(g*-1),b);
0118 FAmap3=cat(3,FAmap,FAmap,FAmap);                                %needed to combine with the colormap, need 3D.
0119 cm=abs(rgb).*FAmap3;                                            %FA weighting for colormap.  Has absolute values for vector components
0120 
0121 %eval (['rgb' num2str(numslice) '=rgb;']);
0122 %nameofmap=['rgb' num2str(numslice)];
0123 %assignin('base',nameofmap,rgb);
0124 
0125 
0126 FAmaps=cat(3,FAmaps,FAmap);
0127 tADCmaps=cat(3,tADCmaps,tADC);
0128 lambda1_maps=cat(3,lambda1_maps,lambda1);
0129 lambda2_maps=cat(3,lambda2_maps,lambda2);
0130 lambda3_maps=cat(3,lambda3_maps,lambda3);
0131 cmaps=cat(4,cmaps,cm);
0132 xyzmaps=cat(4,xyzmaps,xyz);
0133 
0134 
0135 
0136 
0137 
0138 
0139 %eval (['xyz' num2str(numslice) '=rgb;']);                       %xyz is the same as cm (colormap) but without abs.  But, apparently, doesn't make a difference.
0140 %nameofmap=['xyz' num2str(numslice)];
0141 %assignin('base',nameofmap,rgb);
0142 
0143 %image (cm);axis image
0144 %eval (['FA' num2str(numslice) '=FAmap;']);
0145 %nameofmap=['FA' num2str(numslice)];
0146 %assignin('base',nameofmap,FAmap);
0147 if (doFAmap==1)
0148     number=num2str(numslice);                                   %FOUR LINES FOR FILE WRITING, REPEAT FOR THE OTHER IMAGES
0149     st1='FA';st2='.dcm';
0150     nameforfile=[st1 number st2];
0151     dicomwrite(FAmap,nameforfile);
0152     disp('writing FA map to DICOM file')
0153 end
0154 %eval (['tADC' num2str(numslice) '=tADC;']);
0155 %nameofmap=['tADC' num2str(numslice)];
0156 %assignin('base',nameofmap,tADC);
0157 if (dotADC==1)
0158     number=num2str(numslice);                                   %FOUR LINES FOR FILE WRITING, REPEAT FOR THE OTHER IMAGES
0159     st1='tADC';st2='.dcm';
0160     nameforfile=[st1 number st2];
0161     dicomwrite(tADC,nameforfile);
0162     disp('writing tADC map to DICOM file')
0163 end
0164 %eval (['lambda1' num2str(numslice) '=lambda1;']);
0165 %nameofmap=['lambda1_' num2str(numslice)];
0166 %assignin('base',nameofmap,lambda1);
0167 if (dolambda1==1)
0168     number=num2str(numslice);                                   %FOUR LINES FOR FILE WRITING, REPEAT FOR THE OTHER IMAGES
0169     st1='lambda1_';st2='.dcm';
0170     nameforfile=[st1 number st2];
0171     dicomwrite(lambda1,nameforfile);
0172     disp('writing lambda1 map to DICOM file')
0173 end
0174 %eval (['lambda2' num2str(numslice) '=lambda2;']);
0175 %nameofmap=['lambda2_' num2str(numslice)];
0176 %assignin('base',nameofmap,lambda2);
0177 if (dolambda2==1)
0178     number=num2str(numslice);                                   %FOUR LINES FOR FILE WRITING, REPEAT FOR THE OTHER IMAGES
0179     st1='lambda2_';st2='.dcm';
0180     nameforfile=[st1 number st2];
0181     dicomwrite(lambda2,nameforfile); 
0182     disp('writing lambda2 map to DICOM file')
0183 end
0184 %eval (['lambda3' num2str(numslice) '=lambda3;']);
0185 %nameofmap=['lambda3_' num2str(numslice)];
0186 %assignin('base',nameofmap,lambda3);
0187 if (dolambda3==1)
0188     number=num2str(numslice);                                   %FOUR LINES FOR FILE WRITING, REPEAT FOR THE OTHER IMAGES
0189     st1='lambda3_';st2='.dcm';
0190     nameforfile=[st1 number st2];
0191     dicomwrite(lambda3,nameforfile);
0192     disp('writing lambda3 map to DICOM file')
0193 end
0194 %eval (['cm' num2str(numslice) '=cm;']);
0195 %nameofmap=['cm' num2str(numslice)];
0196 %assignin('base',nameofmap,cm);
0197 if (docm==1)
0198     number=num2str(numslice);                                   %FOUR LINES FOR FILE WRITING, REPEAT FOR THE OTHER IMAGES
0199     st1='cm';st2='.dcm';
0200     nameforfile=[st1 number st2];
0201     dicomwrite(cm,nameforfile); 
0202     disp('writing colormap to DICOM file')
0203 end
0204 
0205 mes=[num2str(numslice) ' of ' num2str(slices) ' finished'];
0206 disp(mes)
0207 end
0208  
0209 %assignin('base','ten',ten)
0210 clear FA FAmap tADC lambda1 lambda2 lambda3 cm FAmap3 rgb bo b1 b2 b3 b4 b5 b6 ADC101 ADC_101 ADC011;
0211 clear ADC01_1 ADC110 ADC_110 r g b xx yy zz xy xz yz trace ten tensorvol transform ADCe ADCm D V;
0212 assignin('base','tensorvol',0)
0213 
0214 
0215 assignin('base','FAmaps',FAmaps);
0216 assignin('base','tADCmaps',tADCmaps);
0217 assignin('base','lambda1_maps',lambda1_maps);
0218 assignin('base','lambda2_maps',lambda2_maps);
0219 assignin('base','lambda3_maps',lambda3_maps);
0220 assignin('base','cmaps',cmaps);
0221 assignin('base','xyzmaps',xyzmaps);
0222 
0223 toc