plot cross section of surface on the slice
vb_plot_cross_section(V,F,Z0)
vb_plot_cross_section(V,F,Z0,vdim,dmax,Msize,Mtype,mode)
--- Input
V : voxcel coordinate of the surface [Nvertex x 3]
F : patch index of the surface [Npatch x 3]
Z0 : slice index for plot [1 x 1]
--- Output
h : plotted line handles
--- Optional Input
vdim : slice cut direction
= 'x' : Sagittal cut : Y-Z plane
= 'y' : Coronal cut : X-Z plane
= 'z' : Transverse (Axial) cut : X-Y plane
dmax : search distance from cut slice (= 5)
Msize : Size of plot line (= 1)
Mtype : Color of plot line (= 'y')
mode : 2D plot mode for X-Y (= 0)
= 0 : plot without transpose
= 1 : plot by transposing 2D-image matrix
--- Related function
vb_plot_vertex(V, vdim, Z0, dmax, Msize, Mtype, mode)
--- Intersection point at Z = 0
Line : V1 = [X1 Y1 Z1] to V2 = [X2 Y2 Z2]
: V = (V1 - V2) * t + V2
z = 0 :
x = (X1*Z2 - X2*Z1)/(Z2 - Z1)
y = (Y1*Z2 - Y2*Z1)/(Z2 - Z1)
2006-10-20 M. Sato
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function h = vb_plot_cross_section(V,F,Z0,vdim,dmax,Msize,Mtype,mode) 0002 % plot cross section of surface on the slice 0003 % vb_plot_cross_section(V,F,Z0) 0004 % vb_plot_cross_section(V,F,Z0,vdim,dmax,Msize,Mtype,mode) 0005 % --- Input 0006 % V : voxcel coordinate of the surface [Nvertex x 3] 0007 % F : patch index of the surface [Npatch x 3] 0008 % Z0 : slice index for plot [1 x 1] 0009 % --- Output 0010 % h : plotted line handles 0011 % --- Optional Input 0012 % vdim : slice cut direction 0013 % = 'x' : Sagittal cut : Y-Z plane 0014 % = 'y' : Coronal cut : X-Z plane 0015 % = 'z' : Transverse (Axial) cut : X-Y plane 0016 % dmax : search distance from cut slice (= 5) 0017 % Msize : Size of plot line (= 1) 0018 % Mtype : Color of plot line (= 'y') 0019 % mode : 2D plot mode for X-Y (= 0) 0020 % = 0 : plot without transpose 0021 % = 1 : plot by transposing 2D-image matrix 0022 % --- Related function 0023 % vb_plot_vertex(V, vdim, Z0, dmax, Msize, Mtype, mode) 0024 % 0025 % --- Intersection point at Z = 0 0026 % Line : V1 = [X1 Y1 Z1] to V2 = [X2 Y2 Z2] 0027 % : V = (V1 - V2) * t + V2 0028 % z = 0 : 0029 % x = (X1*Z2 - X2*Z1)/(Z2 - Z1) 0030 % y = (Y1*Z2 - Y2*Z1)/(Z2 - Z1) 0031 % 0032 % 2006-10-20 M. Sato 0033 % 0034 % Copyright (C) 2011, ATR All Rights Reserved. 0035 % License : New BSD License(see VBMEG_LICENSE.txt) 0036 0037 if ~exist('vdim','var'), vdim = 'z'; end; 0038 if ~exist('dmax','var'), dmax = 5; end; 0039 if ~exist('Mtype','var'), Mtype = 'y-'; end; 0040 if ~exist('Msize','var'), Msize = 1; end; 0041 if ~exist('mode','var'), mode = 0; end; 0042 0043 h = []; 0044 0045 switch vdim 0046 case 'x', 0047 % 'SAG' : Sagittal cut : Y-Z plane 0048 ax = 2; 0049 ay = 3; 0050 az = 1; 0051 case 'y', 0052 % 'COR' : Coronal cut : X-Z plane 0053 ax = 1; 0054 ay = 3; 0055 az = 2; 0056 case 'z', 0057 % 'TRN' : Transverse (Axial) cut : X-Y plane 0058 ax = 1; 0059 ay = 2; 0060 az = 3; 0061 end 0062 0063 switch mode 0064 case 0, 0065 X = V(:,ax); 0066 Y = V(:,ay); 0067 case 1, 0068 X = V(:,ay); 0069 Y = V(:,ax); 0070 end 0071 0072 % slice axis coordinate 0073 Z = V(:,az) - Z0; 0074 0075 % Flag indicate upper & lower side of slice within dmax 0076 Uix = (Z >= 0 & Z <= dmax); % flag of upper side of slice 0077 Lix = (Z < 0 & Z >= -dmax);% flag of lower side of slice 0078 0079 % Flag for patch index 0080 FUix = Uix(F); 0081 FLix = Lix(F); 0082 0083 % Patch include both upper and lower side of slice within dmax 0084 Flist = find(sum(FUix,2) > 0 & sum(FLix,2) > 0); 0085 Nlist = length(Flist); 0086 0087 for j=1:Nlist 0088 n = Flist(j); % patch id 0089 id = F(n,:); % vertex index of patch 0090 up = FUix(n,:); % Upper side flag of each vertex in the n-th patch 0091 dn = FLix(n,:); % Lower side flag of each vertex in the n-th patch 0092 x = []; % x-coordinate to intersect the slice Z = 0 0093 y = []; % y-coordinate to intersect the slice Z = 0 0094 0095 % Check two vertex of edge are in opposite side 0096 % : Check edge intersect the slice Z = 0 0097 % and add intersection points to the list 0098 if (up(1)*dn(2) + up(2)*dn(1)) > 0, 0099 % Two vertex points of edge 0100 X1 = X(id(1)); X2 = X(id(2)); 0101 Y1 = Y(id(1)); Y2 = Y(id(2)); 0102 Z1 = Z(id(1)); Z2 = Z(id(2)); 0103 0104 if Z1 == Z2, 0105 % line is on the Z = 0 slice 0106 x = [x; X1; X2]; 0107 y = [y; Y1; Y2]; 0108 else 0109 % Intersection point at Z = 0 0110 x = [x; (X1*Z2 - X2*Z1)/(Z2 - Z1)]; 0111 y = [y; (Y1*Z2 - Y2*Z1)/(Z2 - Z1)]; 0112 end 0113 end 0114 0115 if (up(2)*dn(3) + up(3)*dn(2)) > 0, 0116 % Two vertex points of edge 0117 X1 = X(id(2)); X2 = X(id(3)); 0118 Y1 = Y(id(2)); Y2 = Y(id(3)); 0119 Z1 = Z(id(2)); Z2 = Z(id(3)); 0120 0121 if Z1 == Z2, 0122 % line is on the Z = 0 slice 0123 x = [x; X1; X2]; 0124 y = [y; Y1; Y2]; 0125 else 0126 % Intersection point at Z = 0 0127 x = [x; (X1*Z2 - X2*Z1)/(Z2 - Z1)]; 0128 y = [y; (Y1*Z2 - Y2*Z1)/(Z2 - Z1)]; 0129 end 0130 end 0131 0132 if (up(3)*dn(1) + up(1)*dn(3)) > 0, 0133 % Two vertex points of edge 0134 X1 = X(id(3)); X2 = X(id(1)); 0135 Y1 = Y(id(3)); Y2 = Y(id(1)); 0136 Z1 = Z(id(3)); Z2 = Z(id(1)); 0137 0138 if Z1 == Z2, 0139 % line is on the Z = 0 slice 0140 x = [x; X1; X2]; 0141 y = [y; Y1; Y2]; 0142 else 0143 % Intersection point at Z = 0 0144 x = [x; (X1*Z2 - X2*Z1)/(Z2 - Z1)]; 0145 y = [y; (Y1*Z2 - Y2*Z1)/(Z2 - Z1)]; 0146 end 0147 end 0148 0149 % Plot intersection line of triangle with slice 0150 h = [h ; plot(x,y,Mtype,'MarkerSize',Msize)]; 0151 hold on 0152 end