represent a volume given two surfaces

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?

[Walter Roberson]

Do you just need the outside "skin" drawn; or do you need the inside drawn as long blocks with rectangular cross-section; or do you want the inside filled with a bunch of cubes of uniform size ?

[Chad Greene]

Great question Walter, I should have made that more clear. The skin is sufficient.

[Walter Roberson]

Any NaN or inf values? If not then the skin to be drawn would just be based on the rectangle around the X and Y limits, at least if the grids are regular spacing. Are the X Y arrays perhaps not regular spacing ?

[Chad Greene]

Yes, my volume is embedded in a bunch of NaN values on an irregularly-spaced X,Y grid.

Answer 1

[Cedric]

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.

[Chad Greene]

Cedric, thank you so much for this thoughtful response. This solution worked perfectly. You have really helped me out quite a bit. The only edit I had to make was to also include a surface for the top and bottom of my volume. For the top and bottom I used

surf(X,Y,Z1)

and

surf(X,Y,Z2)

Thanks again for your help!

[Cedric]

You're welcome! Note that my example was really for giving you a method for finding the boundary. Once you have it as well as top and bottom surfaces, you could build something a little more powerful than simple patches, based on a 3D mesh for example.

[Pooyan]

dear cedric...is there any code to fill whole between two surface and build a 3d shape...i mean not a hollow shape...i want to do this...