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