fillSphericalTriangle

FILLSPHERICALTRIANGLE Fill a triangle on a sphere.

   fillSphericalTriangle(SPHERE, PT1, PT2, PT3);


   See also
   fillSphericalPolygon, drawSphericalTriangle, drawSphere

   ---------
   author : David Legland
   INRA - TPV URPOI - BIA IMASTE
   created the 22/02/2005

0001 function varargout = fillSphericalTriangle(sphere, p1, p2, p3, varargin)
0002 %FILLSPHERICALTRIANGLE Fill a triangle on a sphere.
0003 %
0004 %   fillSphericalTriangle(SPHERE, PT1, PT2, PT3);
0005 %
0006 %
0007 %   See also
0008 %   fillSphericalPolygon, drawSphericalTriangle, drawSphere
0009 %
0010 %   ---------
0011 %   author : David Legland
0012 %   INRA - TPV URPOI - BIA IMASTE
0013 %   created the 22/02/2005
0014 %
0015 
0016 %   HISTORY
0017 %   27/06/2007 manage spheres other than origin
0018 %   30/10/2008 replace intersectPlaneLine by intersectLinePlane
0019 %   2012-02-09 rename as fillSphericalTriangle
0020 
0021 % extract data of the sphere
0022 ori = sphere(:, 1:3);
0023 r   = sphere(4);
0024 
0025 % extract direction vectors for each point
0026 v1  = normalizeVector3d(p1 - ori);
0027 v2  = normalizeVector3d(p2 - ori);
0028 v3  = normalizeVector3d(p3 - ori);
0029 
0030 % create a plane tangent to the sphere containing first point
0031 plane = createPlane(v1, v1);
0032 
0033 % position on the plane of the direction vectors
0034 pp2 = planePosition(intersectLinePlane([ori v2], plane), plane);
0035 pp3 = planePosition(intersectLinePlane([ori v3], plane), plane);
0036 
0037 % create rough parametrization with 2 variables
0038 nTri = 5;
0039 s  = linspace(0, 1, nTri);
0040 t  = linspace(0, 1, nTri);
0041 ns = length(s);
0042 nt = length(t);
0043 s  = repmat(s, [nt, 1]);
0044 t  = repmat(t', [1, ns]);
0045 
0046 % convert to plane coordinates
0047 xp = s * pp2(1) + t .* (1-s) * pp3(1);
0048 yp = s * pp2(2) + t .* (1-s) * pp3(2);
0049 
0050 % convert to 3D coordinates (still on the 3D plane)
0051 x  = plane(1) * ones(size(xp)) + plane(4) * xp + plane(7) * yp - ori(1);
0052 y  = plane(2) * ones(size(xp)) + plane(5) * xp + plane(8) * yp - ori(2);
0053 z  = plane(3) * ones(size(xp)) + plane(6) * xp + plane(9) * yp - ori(3);
0054 
0055 % project on the sphere
0056 norm = hypot(hypot(x, y), z);
0057 xn = x ./ norm * r + ori(1);
0058 yn = y ./ norm * r + ori(2);
0059 zn = z ./ norm * r + ori(3);
0060 
0061 
0062 if nargout == 0
0063     % simply display the patch
0064     surf(xn, yn, zn, 'FaceColor', 'g', 'EdgeColor', 'none', varargin{:});
0065     
0066 elseif nargout == 1
0067     % display the patch and return a handle
0068     h = surf(xn, yn, zn, 'FaceColor', 'g', 'EdgeColor', 'none', varargin{:});
0069     varargout = {h};
0070     
0071 elseif nargout == 3
0072     % If 3 outputs are required, return patch vertex coordinates
0073     varargout = {x, y, z};
0074 end