CREATEPLANE Create a plane in parametrized form. PLANE = createPlane(P1, P2, P3) creates a plane containing the 3 points PLANE = createPlane(PTS) The 3 points are packed into a single 3-by-3 array. PLANE = createPlane(P0, N); Creates a plane from a point P0 and a normal N to the plane. The parameter N is given either as a 3D vector (1-by-3 row vector), or as [THETA PHI], where THETA is the colatitute (angle with the vertical axis) and PHI is angle with Ox axis, counted counter-clockwise (both given in radians). PLANE = createPlane(P0, Dip, DipDir); Creates a plane from a point and from a dip and dip direction angles of the plane. Parameters Dip and DipDir angles are given as numbers. Dip : maximum inclination to the horizontal. DipDir : direction of the horizontal trace of the line of dip, measured clockwise from north. The created plane data has the following format: PLANE = [X0 Y0 Z0 DX1 DY1 DZ1 DX2 DY2 DZ2], with - (X0, Y0, Z0) is a point belonging to the plane - (DX1, DY1, DZ1) is a first direction vector - (DX2, DY2, DZ2) is a second direction vector The 2 direction vectors are normalized and orthogonal. See also planes3d, medianPlane
0001 function plane = createPlane(varargin) 0002 %CREATEPLANE Create a plane in parametrized form. 0003 % 0004 % PLANE = createPlane(P1, P2, P3) 0005 % creates a plane containing the 3 points 0006 % 0007 % PLANE = createPlane(PTS) 0008 % The 3 points are packed into a single 3-by-3 array. 0009 % 0010 % PLANE = createPlane(P0, N); 0011 % Creates a plane from a point P0 and a normal N to the plane. The 0012 % parameter N is given either as a 3D vector (1-by-3 row vector), or as 0013 % [THETA PHI], where THETA is the colatitute (angle with the vertical 0014 % axis) and PHI is angle with Ox axis, counted counter-clockwise (both 0015 % given in radians). 0016 % 0017 % PLANE = createPlane(P0, Dip, DipDir); 0018 % Creates a plane from a point and from a dip and dip direction angles 0019 % of the plane. Parameters Dip and DipDir angles are given as numbers. 0020 % Dip : maximum inclination to the horizontal. 0021 % DipDir : direction of the horizontal trace of the line of dip, 0022 % measured clockwise from north. 0023 % 0024 % The created plane data has the following format: 0025 % PLANE = [X0 Y0 Z0 DX1 DY1 DZ1 DX2 DY2 DZ2], with 0026 % - (X0, Y0, Z0) is a point belonging to the plane 0027 % - (DX1, DY1, DZ1) is a first direction vector 0028 % - (DX2, DY2, DZ2) is a second direction vector 0029 % The 2 direction vectors are normalized and orthogonal. 0030 % 0031 % See also 0032 % planes3d, medianPlane 0033 0034 % ------ 0035 % Author: David Legland 0036 % E-mail: david.legland@inrae.fr 0037 % Created: 2005-02-18 0038 % Copyright 2005-2024 INRA - TPV URPOI - BIA IMASTE 0039 0040 if isscalar(varargin) 0041 % If a single input is provided, it can be: 0042 % * an array of three points belonging to the plane 0043 % * a cell array -> one plane per array element is created 0044 0045 var = varargin{1}; 0046 if iscell(var) 0047 plane = zeros([length(var) 9]); 0048 for i = 1:length(var) 0049 plane(i,:) = createPlane(var{i}); 0050 end 0051 0052 elseif size(var, 1) >= 3 0053 % 3 points in a single array 0054 p1 = var(1,:); 0055 p2 = var(2,:); 0056 p3 = var(3,:); 0057 0058 % create direction vectors 0059 v1 = p2 - p1; 0060 v2 = p3 - p1; 0061 0062 % create plane 0063 plane = normalizePlane([p1 v1 v2]); 0064 0065 else 0066 error ('MatGeom:createPlane', 'In case of single a.'); 0067 end 0068 0069 elseif length(varargin) == 2 0070 % Two arguments -> correspond to plane origin and plane normal 0071 0072 % plane origin 0073 p0 = varargin{1}; 0074 0075 % second parameter is either a 3D vector or a 3D angle (2 params) 0076 var = varargin{2}; 0077 if size(var, 2) == 2 0078 % normal is given in spherical coordinates 0079 n = sph2cart2([var ones(size(var, 1))]); 0080 elseif size(var, 2) == 3 0081 % normal is given by a 3D vector 0082 n = normalizeVector3d(var); 0083 else 0084 error ('MatGeom:createPlane', 'Can not parse plane normal.'); 0085 end 0086 0087 % ensure same dimension for parameters 0088 if size(p0, 1) == 1 0089 p0 = repmat(p0, [size(n, 1) 1]); 0090 end 0091 if size(n, 1) == 1 0092 n = repmat(n, [size(p0, 1) 1]); 0093 end 0094 0095 % find a vector not colinear to the normal, in the direction of [1 0 0] 0096 % first try with vector [0 0 1] 0097 v0 = repmat([0 0 1], [size(p0, 1) 1]); 0098 % if vectors are close to colinear, use vector [0 -1 0] 0099 inds = vectorNorm3d(cross(n, v0, 2)) < 1e-12; 0100 v0(inds, :) = repmat([0 -1 0], [sum(inds) 1]); 0101 0102 % create direction vectors 0103 v1 = normalizeVector3d(cross(n, v0, 2)); 0104 v2 = -normalizeVector3d(cross(v1, n, 2)); 0105 0106 % concatenate results in the array representing the plane 0107 plane = [p0 v1 v2]; 0108 0109 elseif length(varargin) == 3 0110 % Three input arguments: 0111 % * three points (as 1-by-3 or 1-by-N numeric arrays) 0112 % * center, Dip and DipDir (?) 0113 0114 var1 = varargin{1}; 0115 var2 = varargin{2}; 0116 var3 = varargin{3}; 0117 if size(var1, 2) == 3 && size(var2, 2) == 3 && size(var3, 2) == 3 0118 % input arguments are three points 0119 p1 = var1; 0120 p2 = var2; 0121 p3 = var3; 0122 0123 % create direction vectors 0124 v1 = p2 - p1; 0125 v2 = p3 - p1; 0126 0127 plane = normalizePlane([p1 v1 v2]); 0128 0129 elseif size(var1, 2) == 3 && size(var2, 2) == 1 && size(var3, 2) == 1 0130 p0 = var1; 0131 n = [sin(var2)*sin(var3) sin(var2)*cos(var3) cos(var2)]; 0132 0133 % find a vector not colinear to the normal 0134 v0 = repmat([1 0 0], [size(p0, 1) 1]); 0135 inds = vectorNorm3d(cross(n, v0, 2))<1e-14; 0136 v0(inds, :) = repmat([0 1 0], [sum(inds) 1]); 0137 0138 % create direction vectors 0139 v1 = normalizeVector3d(cross(n, v0, 2)); 0140 v2 = -normalizeVector3d(cross(v1, n, 2)); 0141 0142 % concatenate result in the array representing the plane 0143 plane = [p0 v1 v2]; 0144 0145 else 0146 error('MatGeom:createPlane', 'Wrong argument in "createPlane".'); 0147 end 0148 else 0149 error ('MatGeom:createPlane', 'Wrong number of input parameters.'); 0150 end 0151