www.gusucode.com > mbcguitools 工具箱 matlab 源码程序 > mbcguitools/mbcInOutTriangles.m
function [insideFaces,outsideFaces,vertices,normals,colors] = mbcInOutTriangles(x,y,z,b,c) %mbcInOutTriangles creates internal and external triangles based on B<0 for xregsurfaceb % [insideFaces,outsideFaces,vertices,normals,colors] = mbcInOutTriangles(x,y,Z,b,c); % x and y are the grid inputs for the output Z, % b is the grid of points to separate the regions based on b<0. % c is the colors specified as a colormap, and or a matrix the same size as % z and b matrices % % [insideFaces,outsideFaces,vertices,normals,colors] = mbcInOutTriangles(x,y,Z,b); % uses colors = z % % B specifies the separating function where B<0 will be inside and B>0 % outside the region. % The vertices, normals and colors are the same for the inside and % outside regions. % % mbcInOutTriangles is used to calculate the patches for xregsurfaceb % % See also xregsurfaceb % Copyright 2014 The MathWorks, Inc. and Ford Global Technologies, Inc. narginchk(5,5) % Find points on the boundary [points,edges] = i_findBoundaryLine(x,y,b); % Triangulate [xx, yy] = ndgrid( x, y ); [DT,TriPoints,inside] = i_triangulation(xx,yy,points,edges,b); % Make Patches zi = interp2( xx.', yy.', z.', TriPoints(:,1), TriPoints(:,2), 'linear' ); vertices = [TriPoints, zi]; [normals, colors] = i_computepatchdata(vertices, x, y, z, c); insideFaces =DT.ConnectivityList(inside,:); outsideFaces =DT.ConnectivityList(~inside,:); end function [points,edges] = i_findBoundaryLine(x,y,b) % points = verticies on the boundary [numpts, 2] % edges = indexes into points [startIndex, endIndex], for every line % segment C = contourc( x, y, b.', [0, 0] ); C = C.'; % Break the contourc matrix C into line segments (find indices) points = zeros( 0, 2 ); % [x, y] numPoints = 0; edges = zeros( 0, 2 ); % [startIdx, endIdx] % Run pointer down the rows of C ptr = 1; while ptr < size( C, 1 ) % The number of vertices in this segment numVertices = C(ptr,2); % Get the list of vertices for this segment indicesForThisSegment = ptr+(1:numVertices); points = [points; C(indicesForThisSegment,:)]; %#ok<AGROW> % Pull out the line segments within this segment newIdx = numPoints+(1:numVertices).'; edges = [edges; [newIdx(1:end-1), newIdx(2:end)]]; %#ok<AGROW> numPoints = numPoints + numVertices; % Move the pointer to the next segment ptr = ptr + numVertices + 1; end end function [DT,TriPoints,inside] = i_triangulation(xx,yy,points,edges,b) % DT = delaunayTriangulation output % points are points on the contour b=0 % inside = boolean vector, true if triangle of DT is inside boundary Tpoints = [points; xx(:),yy(:)]; % scale points to [0,1] based on range of X and Y minPoints = min(Tpoints,[],1); rngPoints = max(Tpoints,[],1)-minPoints; Tpoints = scalePoints(Tpoints,minPoints,rngPoints); warning('OFF', 'MATLAB:delaunayTriangulation:DupPtsConsUpdatedWarnId') DT = delaunayTriangulation( Tpoints, edges ); warning('ON', 'MATLAB:delaunayTriangulation:DupPtsConsUpdatedWarnId') incenterPoints = DT.incenter(); % Need to use unscaled x and y values incenterPoints = unscalePoints(incenterPoints(:,1:2),minPoints,rngPoints); % find interpolated boundary values at center points of trianulation. bT = interp2( xx.', yy.', b.', incenterPoints(:,1), incenterPoints(:,2), 'linear' ); % split patches into inside and outside. % points inside boundary inside = bT < 0; % unscaled triangulation points TriPoints = unscalePoints(DT.Points,minPoints,rngPoints); end function Y = scalePoints(Y,minPoints,rngPoints) %scalePoints scale points from [minPoints , minPoints+rngPoints] to [0,1] Y = bsxfun(@rdivide,bsxfun(@minus,Y,minPoints),rngPoints); end function Y = unscalePoints(Y,minPoints,rngPoints) %unscalePoints unscale points from [0,1] to [minPoints , minPoints+rngPoints] Y = bsxfun(@plus,bsxfun(@times,Y,rngPoints),minPoints); end function [normals, colors] = i_computepatchdata(vertices, x, y, z, c) xs = vertices(:,1); ys = vertices(:,2); % Compute surface normals [dy, dx] = gradient( z, y, x); dx = interp2( y, x, dx, ys, xs, 'linear' ); dy = interp2( y, x, dy, ys, xs, 'linear' ); normals = [dx(:), dy(:), -ones( numel( dx ), 1 )]; nrm = sqrt( sum( normals.^2, 2 ) ); normals(:,1) = normals(:,1)./nrm; normals(:,2) = normals(:,2)./nrm; normals(:,3) = normals(:,3)./nrm; if ismatrix(c) colors = interp2(y, x, c, ys, xs, 'linear'); elseif ndims(c)==3 && size(c,3)==3 % colormap colors = zeros(length(xs), 3); for k = 1:3 colors(:,k) = interp2( y, x, c(:,:,k), ys, xs ); end else error(message('mbc:xregsurfaceb:InvalidArgument1')); end end