Home > vbmeg > functions > plotfunc > subdirectory > vb_get_cross_section.m

vb_get_cross_section

PURPOSE ^

plot cross section of surface on the slice

SYNOPSIS ^

function [XP,YP,XP1,YP1,XP2,YP2,XP3,YP3,NN] =vb_get_cross_section(V,F,Z0,vdim,dmax,xymode)

DESCRIPTION ^

 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)

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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);

Generated on Mon 22-May-2023 06:53:56 by m2html © 2005