Subdividing a 3D surface in MATLAB

Question

I am using MATLAB to create some diagrams to explain my research and am stuck on the following problem. I have created a half-ellipsoid rotated about the x-axis but I would like to demonstrate some internal structure within the ellipsoid which will be rotated around the x- and y-axes. This is roughly shown in the figure below, where I have subdivided the ellipsoid into four parts (with my bad paint skill).

How can I subdivide my ellipsoids into sections using planes of varying rotation about the x- and y-axes? I can create the ellipsoid and intersecting plane surfaces, but from there I don't know how to divide the ellipsoid up and change the face colours.

I have included some basic code below to get things started. I had thought I could just mask parts of the ellipsoid coordinates like I did to cut the ellipsoid in half, but this does not work. I am guessing I need to make some kind of mesh grid,but I don't know how to combine the grid from the intersecting surface and the ellipsoid.

figure
x   = 0;    y   = 0;    z   = 0;
tl  = 10;   tw  = 4;   td  = 2;

% Create ellipsoid
[ex,ey,ez]  = ellipsoid(x, y, z, tl, tw, td,40);
ex          = ex(1:ceil(length(ez)/2),:);   % Remove top half 
ey          = ey(1:ceil(length(ez)/2),:);   % of ellipsoid
ez          = ez(1:ceil(length(ez)/2),:);

% Make some planes
[ySL,zSL] = meshgrid([-10:10],[-2:0.2:2]);
xSL1 = zeros(size(ySL, 1)); % Generate z data
hSL1 = surf(xSL1,ySL,zSL);

hold on
[ySL,zSL] = meshgrid([-10:10],[-3:0.2:1]);
xSL2 = ones(size(ySL, 1)); % Generate z data
hSL2 = surf(xSL2,ySL,zSL);

% rotate(hSL,[1 0 0],5);
rotate([hSL1 hSL2],[0 1 0],-70);

hSurf1  = surf(ex,ey,ez);
set([hSurf1 hSL1 hSL2],'facecolor','blue','facealpha',.2,...
    'edgecolor','none')

% Plot settings
daspect([1 1 0.3]);
hold off
view(-10,6)

Any help is greatly appreciated,

Answer 1

This was an interesting question. So, here is a solution for 1) plotting three dimensional points inside the ellipse, bounded by the two planes, and 2) taking those points to form a smooth looking patch representation.

My basic idea is to build a mesh that covers (as densely as required) the whole X-Y-Z space of interest. What I then do is define the two planes in the parametric form ax+by+cz+d=0. These are the planes between which the data is to be plotted. Using the normal vector of the planes, with information about a point on each of the planes, I can deduce using the dot product whether each point in the X-Y-Z grid is above or below the planes: if the dot product is less than zero, the point is above the plane, and if greater than zero, above the plane.

Using logical indexing I then find which points on the grid satisfy the conditions 1) above plane 1, 2) below plane 2, and 3) inside the ellipsoid and -- optionally -- 4) below Z=0.

Finally, I plot the ellipsoid and the planes, along with the points that fulfil the criteria above.

In this case, I define two planes, but I guess you could do this for any number of them. The code is quick and dirty, but hopefully will give you a way forward!

% Define the ellipse
x   = 0;    y   = 0;    z   = 0;
tl  = 10;   tw  = 4;   td  = 2;

% Create ellipsoid
[ex,ey,ez]  = ellipsoid(x, y, z, tl, tw, td,40);
ex          = ex(1:ceil(length(ez)/2),:);   % Remove top half 
ey          = ey(1:ceil(length(ez)/2),:);   % of ellipsoid
ez          = ez(1:ceil(length(ez)/2),:);


% Define a 3D grid over area of interest
xx = linspace(min(ex(:)),max(ex(:)),50);
yy = linspace(min(ey(:)),max(ey(:)),50);
zz = linspace(min(ez(:)),max(ez(:)),50);
[X, Y, Z] = meshgrid(xx, yy, zz);
V = [X(:) Y(:) Z(:)]; % rearrange

% Define two planes (ax + bx + cz + d = 0)
a = [0; 0];
b = [-1; 1];
c = [-1; 1];
d = [-1; -1];
% Normal vectors of the planes
n = [a b c];
n(1,:) = n(1,:) / norm(n(1,:));
n(2,:) = n(2,:) / norm(n(2,:));

% Find a point (0,0,z) on each plane
zp = [zeros(2) (d- a * 0 - b * 0)./c];

% Define the area of interest.
% Dot product to test if grid points are below or above the planes
% We want: above plane 1, below plane 2.
V1 = sum(bsxfun(@times, bsxfun(@minus, V, -zp(1,:)), n(1,:)),2);
V2 = sum(bsxfun(@times, bsxfun(@minus, V, -zp(2,:)), n(2,:)),2);
between_planes = (V1 < 0) & (V2 < 0);
% ...and the points have to be inside the ellipsoid
in_ellipse = ((X(:) - x)/tl).^2 + ((Y(:)-y)/tw).^2 + ((Z(:)-z)/td).^2 < 1;
% Final AOI
aoi = between_planes & in_ellipse;

% Add this if you want to also have only values with Z < 0
aoi = aoi & (V(:,3) < 0);

figure;
surf(ex, ey, ez, 'facecolor','blue','facealpha',.2,...
    'edgecolor','none')

So, now we are ready to plot some data! First up, let's try just plotting the grid to get a kind of discrete visualization with individual data points showing up:

hold on;
plot3(V(aoi,1), V(aoi, 2), V(aoi, 3), 'r.')

...just to see it works, let's add the planes we defined earlier to the visualization. This is something you probably want to get rid of eventually.

% Draw the planes
[X,Y] = meshgrid(xx,yy);
Z = cell(2,1);
clr = {'r', 'g'};
for k = 1:2
    Z{k} = (a(k) * X + b(k) * Y + d(k))/ (-c(k));
    hs(k) = surf(xx, yy, Z{k},'facecolor',clr{k}, 'facealpha',.2, 'edgecolor','none');
    drawnow;
end
view([-115 10])
legend(hs, 'Plane 1 (above)', 'Plane 2 (below');

Well, finally, if we want to show some smooth looking data, we can use alphaShape to extract a polyhedra from those data points we have. We can visualize the polyhedra as a patch object just by using plot. Finally, we can assign any color we want to each of the vertices in the patch, set a smooth change to the visualized color by setting the 'Facecolor' property, and remove the pesky vertex edges by setting 'EdgeColor' to 'none':

% Moving on to smooth representation...
shp = alphaShape(V(aoi,1), V(aoi, 2), V(aoi, 3),1);
hs = plot(shp);
% Create a colormap to use for the patch (color of each vertex)
cmap = jet(size(hs.Vertices,1));
set(hs, 'FaceVertexCData', cmap);
set(hs, 'FaceColor', 'interp'); % smooth interp coloring
set(hs, 'EdgeColor', 'none'); % Remove ugly edges

So that's it. All that is left is to replace the vertex color values with whatever is appropriate for your data. Also, remember to skip the plot3 part earlier if you don't want those data points to show up in the display along the patch.

Sample end result below, showing the two planes as well.