MATLAB: Represent a volume given two surfaces

volume visualization

I have two surfaces given by the grids X, Y, Z1 and X, Y, Z2. Both Z grids share the same X and Y grids. I can plot each Z surface, but I'd like to represent the space between them as a solid mass. Is there a way to represent the space between these two Z surfaces a simple volume akin to a patch or a fill in 2D space?

Best Answer

Saturday EDIT: I took 5 minutes to build an example closer to your case; see the new section after my initial answer.

=== Initial answer.

The issue is to locate the boundary. Once you have it, it's easy: here is an example assuming that the boundary is known (cylinder with octagonal base):

 % - Define (x,y) base and z for the example.
 Xbnd  = [2 1.41 0 -1.41 -2 -1.41 0 1.41 2].' ;
 Ybnd  = [0 1.41 2 1.41 0 -1.41 -2 -1.41 0].' ;
 Z1bnd = rand(size(Xbnd)) ;  Z1bnd(end) = Z1bnd(1) ;
 Z2bnd = 1 + rand(size(Xbnd)) ;  Z2bnd(end) = Z2bnd(1) ;
 % - Define patches nodes (from now on, it is no more specific to the example).
 Xdata = [Xbnd(1:end-1), Xbnd(2:end), Xbnd(2:end), Xbnd(1:end-1)].' ;
 Ydata = [Ybnd(1:end-1), Ybnd(2:end), Ybnd(2:end), Ybnd(1:end-1)].' ;
 Zdata = [Z1bnd(1:end-1), Z1bnd(2:end), Z2bnd(2:end), Z2bnd(1:end-1)].' ;
 % - Plot patches.

 patch(Xdata, Ydata, Zdata, 'm') ;
 rotate3d ;                              % Enable 3D rotations with the mouse 

                                         % in the figure.

As you are dealing with a grid (even if nodes are irregularly spaced, as it doesn't alter the topology), one option for determining the boundary is to work with the image processing toolbox on a BW image generated with ~isnan(Z). Once you get row/col IDs of pixels of the boundary using BWBOUNDARIES [ 1, 2 ], you can use them to select corresponding elements in X, Y, Z1, and Z2.

So you could start trying doing something like (assuming that the nan values are in Z1 and/or Z2):

 BW = ~(isnan(Z1) | isnan(Z2)) ;
 BW_filled = imfill(BW, 'holes') ;
 boundaries = bwboundaries(BW_filled) ;

then maybe filter out unwanted boundaries, or keep the longest one. Finally, you would have something like

 Xbnd = X(sub2ind(size(X), boundary(:,1), boundary(:,2))) ;
 Ybnd = Y(sub2ind(size(Y), boundary(:,1), boundary(:,2))) ;

where boundary equals boundaries{k} for k the index of the longest boundary, or for all k if you iterate over boundaries.

=== Example closer to your case.

 close all ;  clear all ;  clc ; 
 % - Build demo X, Y, Z1, Z2.
 [X, Y] = meshgrid(1:10, 1:10) ;
 X  = X + rand(10) - 0.5 ;
 Y  = Y + rand(10) - 0.5 ;
 Z1 = sin(X) .* cos(Y) ;
 Z2 = X + Y.^2/10 ;
 idNaN = (X-5).^2 + (Y-5).^2 > 3^2 ;    % Set Z1 to NaN around approx. 
 Z1(idNaN) = NaN ;                      % disk shape.
 % - Get boundary.
 BW = ~(isnan(Z1) | isnan(Z2)) ;
 BW_filled = imfill(BW, 'holes') ;
 boundaries = bwboundaries(BW_filled) ;
 boundary = boundaries{1} ;             % Only 1 boundary in this trivial example.
 % - Define x,y,z1,z2 values on the boundary.
 Xbnd  = X(sub2ind(size(X), boundary(:,1), boundary(:,2))) ;
 Ybnd  = Y(sub2ind(size(Y), boundary(:,1), boundary(:,2))) ;
 Z1bnd = Z1(sub2ind(size(Z1), boundary(:,1), boundary(:,2))) ;
 Z2bnd = Z2(sub2ind(size(Z2), boundary(:,1), boundary(:,2))) ;
 % - Define patches nodes.
 Xdata = [Xbnd(1:end-1), Xbnd(2:end), Xbnd(2:end), Xbnd(1:end-1)].' ;
 Ydata = [Ybnd(1:end-1), Ybnd(2:end), Ybnd(2:end), Ybnd(1:end-1)].' ;
 Zdata = [Z1bnd(1:end-1), Z1bnd(2:end), Z2bnd(2:end), Z2bnd(1:end-1)].' ;
 % - Plot patches.
 patch(Xdata, Ydata, Zdata, 'm') ;
 rotate3d ;                             % Enable 3D rotations with the mouse 
                                        % in the figure.

Related Solutions

MATLAB: How to plot 2d & 3d contour with xy coordinates with corresponding z values

If you had the geometry saved, e.g. as .stl file, then you could just pass some constraints to the delunayTriangulation. However, a small problem is that your edges are not straight, so it becomes difficult to define cornerpoints.

Here's another method you could try:

xyz = readmatrix('mesh data.xlsx');
[x,y,z] = deal(xyz(:,1),xyz(:,2),xyz(:,3))
% alphashape to get tight boundary without holes
shp = alphaShape(x,y,'holethreshold',1e5);
% plot edges
edges = boundaryFacets(shp);
plot(x(edges),y(edges),'k','linewidth',2);hold on
%get connectivity
tri = alphaTriangulation(shp);
t = trisurf(tri,x,y,z);
t.EdgeColor = 'none'
view(2)

MATLAB: Filling area between two planes in 3d plot

The isosurface approach that Tim suggested is possible. Here's an example of how you'd do that. It's a bit involved and it uses quite a bit of memory:

%%Start with 2 surfaces
nx = 25;
ny = 25;
nz = 25;
[x,y] = meshgrid(linspace(-2,2,nx),linspace(-2,2,ny));
z1 = peaks(x,y);
z2 = peaks(x,y) + 5;
surf(x,y,z1);
hold on
surf(x,y,z2)
%%Create 3D grid containing distance to closest surface
[y3,x3,z3] = ndgrid(linspace(-2,2,ny),linspace(-2,2,nx),linspace(-10,15,nz));
v = zeros(size(z3));
for r=1:ny
    for c=1:nx
        for s=1:nz
            d1 = z3(r,c,s) - z1(r,c);
            d2 = z2(r,c) - z3(r,c,s);
            if d1 < 0
                v3(r,c,s) = d1;
            elseif d2 < 0
                v3(r,c,s) = d2;
            else
                v3(r,c,s) = min(d1,d2);
            end
        end
    end
end
%%Create isosurface
figure
p = [patch(isosurface(x3,y3,z3,v3,0)), ...
     patch(isocaps(x3,y3,z3,v3,0))];
isonormals(x3,y3,z3,v3,p(1)) 
set(p,'FaceColor','yellow')
set(p,'EdgeColor','none')
set(p,'FaceLighting','gouraud')
view(3)
camlight right

You need both isosurface and isocaps there. The isosurface command creates the parts inside the volume (the same parts surf created) while isocaps creates the parts along the boundary. You can also combine isocaps from that example with your surface objects like this:

figure
surf(x,y,z1);
hold on
surf(x,y,z2)
isocaps(x3,y3,z3,v3,0);

Best Answer

Related Solutions

MATLAB: How to plot 2d & 3d contour with xy coordinates with corresponding z values

MATLAB: Filling area between two planes in 3d plot

Related Question