plot cross section of surface on the slice
[X,Y] = vb_get_cross_section(V,F,Z0,vdim,dmax,xymode)
[X,Y,X1,Y1,X2,Y2,X3,Y3] = vb_get_cross_section(V,F,Z0,vdim,dmax,xymode)
--- 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]
--- 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)
xymode : 2D plot mode for X-Y (= 0)
= 0 : plot without transpose
= 1 : plot by transposing 2D-image matrix
--- Output
X : List of two X coordinate for line segment : [2 x NP]
Y : List of one Y coordinate for line segment : [2 x NP]
For most cases X & Y are enough for cross section plot
In special cases following lists of lines are given
X1 : List of one X coordinate for points : [1 x NP]
Y1 : List of two Y coordinate for points : [1 x NP]
X2 : List of 4 X coordinate for line segment : [4 x NP]
Y2 : List of 4 Y coordinate for line segment : [4 x NP]
X3 : List of 6 X coordinate for line segment : [6 x NP]
Y3 : List of 6 Y coordinate for line segment : [6 x NP]
--- Calculation of 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)
2007-4-5 M. Sato
Copyright (C) 2011, ATR All Rights Reserved.
License : New BSD License(see VBMEG_LICENSE.txt)0001 function [XP,YP,XP1,YP1,XP2,YP2,XP3,YP3,NN] = ... 0002 vb_get_cross_section(V,F,Z0,vdim,dmax,xymode) 0003 % plot cross section of surface on the slice 0004 % [X,Y] = vb_get_cross_section(V,F,Z0,vdim,dmax,xymode) 0005 % [X,Y,X1,Y1,X2,Y2,X3,Y3] = vb_get_cross_section(V,F,Z0,vdim,dmax,xymode) 0006 % --- Input 0007 % V : voxcel coordinate of the surface [Nvertex x 3] 0008 % F : patch index of the surface [Npatch x 3] 0009 % Z0 : slice index for plot [1 x 1] 0010 % --- Optional Input 0011 % vdim : slice cut direction 0012 % = 'x' : Sagittal cut : Y-Z plane 0013 % = 'y' : Coronal cut : X-Z plane 0014 % = 'z' : Transverse (Axial) cut : X-Y plane 0015 % dmax : search distance from cut slice (= 5) 0016 % xymode : 2D plot mode for X-Y (= 0) 0017 % = 0 : plot without transpose 0018 % = 1 : plot by transposing 2D-image matrix 0019 % --- Output 0020 % X : List of two X coordinate for line segment : [2 x NP] 0021 % Y : List of one Y coordinate for line segment : [2 x NP] 0022 % For most cases X & Y are enough for cross section plot 0023 % In special cases following lists of lines are given 0024 % X1 : List of one X coordinate for points : [1 x NP] 0025 % Y1 : List of two Y coordinate for points : [1 x NP] 0026 % X2 : List of 4 X coordinate for line segment : [4 x NP] 0027 % Y2 : List of 4 Y coordinate for line segment : [4 x NP] 0028 % X3 : List of 6 X coordinate for line segment : [6 x NP] 0029 % Y3 : List of 6 Y coordinate for line segment : [6 x NP] 0030 % 0031 % --- Calculation of intersection point at Z = 0 0032 % Line : V1 = [X1 Y1 Z1] to V2 = [X2 Y2 Z2] 0033 % : V = (V1 - V2) * t + V2 0034 % z = 0 : 0035 % x = (X1*Z2 - X2*Z1)/(Z2 - Z1) 0036 % y = (Y1*Z2 - Y2*Z1)/(Z2 - Z1) 0037 % 0038 % 2007-4-5 M. Sato 0039 % 0040 % Copyright (C) 2011, ATR All Rights Reserved. 0041 % License : New BSD License(see VBMEG_LICENSE.txt) 0042 0043 if ~exist('vdim','var'), vdim = 'z'; end; 0044 if ~exist('dmax','var'), dmax = 5; end; 0045 if ~exist('xymode','var'), xymode = 0; end; 0046 0047 switch vdim 0048 case 'x', 0049 % 'SAG' : Sagittal cut : Y-Z plane 0050 ax = 2; 0051 ay = 3; 0052 az = 1; 0053 case 'y', 0054 % 'COR' : Coronal cut : X-Z plane 0055 ax = 1; 0056 ay = 3; 0057 az = 2; 0058 case 'z', 0059 % 'TRN' : Transverse (Axial) cut : X-Y plane 0060 ax = 1; 0061 ay = 2; 0062 az = 3; 0063 end 0064 0065 switch xymode 0066 case 0, 0067 X = V(:,ax); 0068 Y = V(:,ay); 0069 case 1, 0070 X = V(:,ay); 0071 Y = V(:,ax); 0072 end 0073 0074 % slice axis coordinate 0075 Z = V(:,az) - Z0; 0076 0077 % Flag indicate upper & lower side of slice within dmax 0078 Uix = (Z >= 0 & Z <= dmax); % flag of upper side of slice 0079 Lix = (Z < 0 & Z >= -dmax);% flag of lower side of slice 0080 0081 % Flag for patch index 0082 FUix = Uix(F); 0083 FLix = Lix(F); 0084 0085 % Patch include both upper and lower side of slice within dmax 0086 Flist = find(sum(FUix,2) > 0 & sum(FLix,2) > 0); 0087 Nlist = length(Flist); 0088 0089 XP = zeros(2,Nlist); 0090 YP = zeros(2,Nlist); 0091 XP1 = []; 0092 YP1 = []; 0093 XP2 = []; 0094 YP2 = []; 0095 XP3 = []; 0096 YP3 = []; 0097 0098 NN = []; 0099 m = 0; 0100 m0 = 0; 0101 m1 = 0; 0102 m2 = 0; 0103 m3 = 0; 0104 0105 for j=1:Nlist 0106 n = Flist(j); % patch id 0107 id = F(n,:); % vertex index of patch 0108 up = FUix(n,:); % Upper side flag of each vertex in the n-th patch 0109 dn = FLix(n,:); % Lower side flag of each vertex in the n-th patch 0110 x = []; % x-coordinate to intersect the slice Z = 0 0111 y = []; % y-coordinate to intersect the slice Z = 0 0112 0113 % Check two vertex of edge are in opposite side 0114 % : Check edge intersect the slice Z = 0 0115 % and add intersection points to the list 0116 if (up(1)*dn(2) + up(2)*dn(1)) > 0, 0117 % Two vertex points of edge 0118 X1 = X(id(1)); X2 = X(id(2)); 0119 Y1 = Y(id(1)); Y2 = Y(id(2)); 0120 Z1 = Z(id(1)); Z2 = Z(id(2)); 0121 0122 if Z1 == Z2, 0123 % line is on the Z = 0 slice 0124 x = [x; X1; X2]; 0125 y = [y; Y1; Y2]; 0126 else 0127 % Intersection point at Z = 0 0128 x = [x; (X1*Z2 - X2*Z1)/(Z2 - Z1)]; 0129 y = [y; (Y1*Z2 - Y2*Z1)/(Z2 - Z1)]; 0130 end 0131 end 0132 0133 if (up(2)*dn(3) + up(3)*dn(2)) > 0, 0134 % Two vertex points of edge 0135 X1 = X(id(2)); X2 = X(id(3)); 0136 Y1 = Y(id(2)); Y2 = Y(id(3)); 0137 Z1 = Z(id(2)); Z2 = Z(id(3)); 0138 0139 if Z1 == Z2, 0140 % line is on the Z = 0 slice 0141 x = [x; X1; X2]; 0142 y = [y; Y1; Y2]; 0143 else 0144 % Intersection point at Z = 0 0145 x = [x; (X1*Z2 - X2*Z1)/(Z2 - Z1)]; 0146 y = [y; (Y1*Z2 - Y2*Z1)/(Z2 - Z1)]; 0147 end 0148 end 0149 0150 if (up(3)*dn(1) + up(1)*dn(3)) > 0, 0151 % Two vertex points of edge 0152 X1 = X(id(3)); X2 = X(id(1)); 0153 Y1 = Y(id(3)); Y2 = Y(id(1)); 0154 Z1 = Z(id(3)); Z2 = Z(id(1)); 0155 0156 if Z1 == Z2, 0157 % line is on the Z = 0 slice 0158 x = [x; X1; X2]; 0159 y = [y; Y1; Y2]; 0160 else 0161 % Intersection point at Z = 0 0162 x = [x; (X1*Z2 - X2*Z1)/(Z2 - Z1)]; 0163 y = [y; (Y1*Z2 - Y2*Z1)/(Z2 - Z1)]; 0164 end 0165 end 0166 0167 NP = length(x); 0168 0169 if NP == 2, 0170 m = m + 1; 0171 XP(:,m) = x; 0172 YP(:,m) = y; 0173 elseif NP == 1, 0174 m1 = m1 + 1; 0175 XP1(:,m1) = x; 0176 YP1(:,m1) = y; 0177 elseif NP == 4, 0178 m2 = m2 + 1; 0179 XP2(:,m2) = x; 0180 YP2(:,m2) = y; 0181 elseif NP == 6, 0182 m3 = m3 + 1; 0183 XP3(:,m3) = x; 0184 YP3(:,m3) = y; 0185 else 0186 m0 = m0 + 1; 0187 NN(m0) = NP; 0188 end 0189 end 0190 0191 XP = XP(:,1:m); 0192 YP = YP(:,1:m); 0193 0194 XP1 = XP1(:,1:m1); 0195 YP1 = YP1(:,1:m1); 0196 0197 XP2 = XP2(:,1:m2); 0198 YP2 = YP2(:,1:m2); 0199 0200 XP3 = XP3(:,1:m3); 0201 YP3 = YP3(:,1:m3);