plot3d

 PLOT3D - produces an image of a 3D object
 
 To use: [b,d]=plot3d(a,alfa,beta);
 
 This function produces an image of a 3D object
 defined by matrix a(l,m,n) in terms of voxels.
 The image is a view after rotating the object 
 by angles alfa and beta (in degrees).
 b is the image and d is its distance to the viewer matrix
 The first figure depicts the object using only its gray level values
 The second image depicts the object using some lighting effect
 rotate3d may be used for reorientation but the obtained image is planar
 Kindly be patient the prog. is very slow!
 July-9-1998 Dr. H Azhari

 Example:
 x=zeros(40,40,40);
 x(10:30,10:30,10:30)=ones(21,21,21);
 x(15:25,15:25,:)=zeros(11,11,40);
 plot3d(x,30,20);
 colormap gray

0001 function [b,d]=plot3d(a,alfa,beta);
0002 % PLOT3D - produces an image of a 3D object
0003 %
0004 % To use: [b,d]=plot3d(a,alfa,beta);
0005 %
0006 % This function produces an image of a 3D object
0007 % defined by matrix a(l,m,n) in terms of voxels.
0008 % The image is a view after rotating the object
0009 % by angles alfa and beta (in degrees).
0010 % b is the image and d is its distance to the viewer matrix
0011 % The first figure depicts the object using only its gray level values
0012 % The second image depicts the object using some lighting effect
0013 % rotate3d may be used for reorientation but the obtained image is planar
0014 % Kindly be patient the prog. is very slow!
0015 % July-9-1998 Dr. H Azhari
0016 %
0017 % Example:
0018 % x=zeros(40,40,40);
0019 % x(10:30,10:30,10:30)=ones(21,21,21);
0020 % x(15:25,15:25,:)=zeros(11,11,40);
0021 % plot3d(x,30,20);
0022 % colormap gray
0023 
0024 
0025 time0=clock;
0026 %-----------------
0027 % Intial values
0028 [l,m,n]=size(a);
0029 M=max(size(a));
0030 M1=round(sqrt(2)*M)+2;
0031 b=zeros(M1,M1);
0032 d=-M1*ones(M1,M1);
0033 alfa=alfa*pi/180;
0034 beta=beta*pi/180;
0035 T=zeros(3,3); % Rotation matrix
0036 T(1,1)=cos(alfa);
0037 T(1,2)=sin(alfa);
0038 T(1,3)=0;
0039 T(2,1)=-cos(beta)*sin(alfa);
0040 T(2,2)=cos(beta)*cos(alfa);
0041 T(2,3)=sin(beta);
0042 T(3,1)=sin(beta)*sin(alfa);
0043 T(3,2)=-sin(beta)*cos(alfa);
0044 T(3,3)=cos(beta);
0045 xyz=[1,1,1]';
0046 %-------------------------
0047 for i=1:l,
0048     for j=1:m,
0049     for k=1:n,
0050     xyz=round(T*[i-l/2,j-m/2,k-n/2]');
0051     x1=round(M1/2+xyz(1)+1);
0052     y1=round(M1/2+xyz(2)+1);
0053     z1=round(M1/2+xyz(3)+1);
0054         if a(i,j,k)>0,
0055             if d(x1,y1)<z1,
0056             b(x1,y1)=a(i,j,k);
0057             d(x1,y1)=z1;
0058             end
0059         end
0060 
0061 
0062     end
0063     end
0064 end
0065 
0066 imagesc(b)
0067 figure
0068 [fx,fy]=gradient(d);
0069 h=(b.*(d.^1.5+mean(mean(d/2))*(cos((atan(fy)+atan(fx))/2)).^2));
0070 pcolor(flipud(h));shading 'interp' ;
0071 %pcolor(flipud((b.*d).^1.5));
0072 %shading 'interp'
0073 
0074 %------------------------------
0075 Elapsed_Time=etime(clock,time0)