4

How could I create a 3D binary matrix/image from a surface mesh in Matlab?

For instance, when I create ellipsoid using:

[x, y, z] = ellipsoid(0,0,0,5.9,3.25,3.25,30);

X, Y and X are all 2D matrix with size 31 x 31.

Edited based on suggestion of @Magla:

function Create_Mask_Basedon_Ellapsoid3()
         close all
        SurroundingVol  = [50, 50, 20];
        %DATA
        [MatX,MatY,MatZ] = meshgrid(-24:1:25, -24:1:25, -9:1:10);
        [mask1, x, y, z] = DrawEllipsoid([0, -10, 0], [6, 3, 3], MatX,MatY,MatZ);
        [mask2, x2, y2, z2] = DrawEllipsoid([15, 14, 6], [6, 3, 3], MatX,MatY,MatZ);
        mask =  mask1 + mask2;

        %Surface PLOT
        figure('Color', 'w');
        subplot(1,2,1);

        %help: Ideally I would like to generate surf plot directly from combined mask= mask1 + mask2;
        s = surf(x,y,z); hold on;
        s2 = surf(x2,y2,z2); hold off;     
        title('SURFACE', 'FontSize', 16);
        view(-78,22)

        subplot(1,2,2);
        xslice = median(MatX(:));
        yslice = median(MatY(:));
        zslice = median(MatZ(:));

        %help: Also how do I decide correct "slice" and angles to 3D visualization.
        h = slice(MatX, MatY, MatZ, double(mask), xslice, yslice, zslice)
        title('BINARY MASK - SLICE VOLUME', 'FontSize', 16);
        set(h, 'EdgeColor','none');   
        view(-78, 22)
        %az = 0; el = 90;
        %view(az, el);

    end

    function [mask, Ellipsoid_x, Ellipsoid_y, Ellipsoid_z] = DrawEllipsoid(CenterEllipsoid, SizeEllipsoid, MatX, MatY, MatZ)
        [Ellipsoid_x, Ellipsoid_y, Ellipsoid_z] =  ellipsoid(CenterEllipsoid(1), CenterEllipsoid(2), CenterEllipsoid(3), SizeEllipsoid(1)/2 , SizeEllipsoid(2)/2 , SizeEllipsoid(3)/2 ,30);
        v = [Ellipsoid_x(:), Ellipsoid_y(:), Ellipsoid_z(:)];                       %3D points
        %v = [x(:), y(:), z(:)];                                                    %3D points
        tri = DelaunayTri(v);                                                       %triangulation
        SI = pointLocation(tri,MatX(:),MatY(:),MatZ(:));                            %index of simplex (returns NaN for all points outside the convex hull)
        mask = ~isnan(SI);                                                          %binary
        mask = reshape(mask,size(MatX));                                            %reshape the mask 
    end
Garima Singh
  • 1,410
  • 5
  • 22
  • 46
  • Like you just want plot3(x,y,z) or somethig else ? You can use cat(3,x,y,z) where 3 is for a 3D concat. – Vuwox Oct 17 '13 at 12:41
  • I could not understand your solution. I am asking how can I make a 3D mask from ellipsoid. – Garima Singh Oct 17 '13 at 13:01
  • What you mean by a mask ? Filter ? I'm not sure to undertand. – Vuwox Oct 17 '13 at 13:33
  • By mask, I mean the following: For all coordinates (x,y,z) interior to ellpisoid, mask(x,y,z) = 0; and otherwise mask(x,y,z) = 1. – Garima Singh Oct 17 '13 at 13:49
  • @Alexandre Bizeau You may be interested in the convex hull solution below that works for any kind of volume - a solution from Gnovice that is kind of clever – marsei Oct 17 '13 at 20:30

2 Answers2

4

There you go:

%// Points you want to test. Define as you need. This example uses a grid of 1e6
%// points on a cube of sides [-10,10]:
[x y z] = meshgrid(linspace(-10,10,100)); 
x = x(:);
y = y(:);
z = z(:); %// linearize

%// Ellipsoid data
center = [0 0 0]; %// center
semiaxes = [5 4 3]; %// semiaxes

%// Actual computation:
inner = (x-center(1)).^2/semiaxes(1).^2 ...
      + (y-center(2)).^2/semiaxes(2).^2 ...
      + (z-center(3)).^2/semiaxes(3).^2 <= 1;

For the n-th point of the grid, whose coordinates are x(n), y(n), z(n), inner(n) is 1 if the point lies in the interior of the ellipsoid and 0 otherwise.

For example: draw the interior points:

plot3(x(inner), y(inner), z(inner), '.' , 'markersize', .5)

enter image description here

Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
3

Here is a method for creating a binary mask from an ellipsoid. It creates a corresponding volume and sets to NaN the points outside the ellipsoid (ones inside).

It doesn't take into consideration the formula of the ellipsoid, but uses a convex hull. Actually, it works for any volume that can be correctly described by a 3D convex hull. Here, the convexhulln step is bypassed since the ellipsoid is already a convex hull.

All credits go to Converting convex hull to binary mask

The following plot

enter image description here

is produced by

%DATA
[x, y, z] = ellipsoid(0,0,0,5.9,3.25,3.25,30);

%METHOD
v = [x(:), y(:), z(:)];                       %3D points
[X,Y,Z] = meshgrid(min(v(:)):0.1:max(v(:)));  %volume mesh
tri = DelaunayTri(v);                         %triangulation
SI = pointLocation(tri,X(:),Y(:),Z(:));       %index of simplex (returns NaN for all points outside the convex hull)
mask = ~isnan(SI);                            %binary
mask = reshape(mask,size(X));                 %reshape the mask 

%PLOT
figure('Color', 'w');
subplot(1,2,1);
s = surf(x,y,z);
title('SURFACE', 'FontSize', 16);
view(-78,22)
subplot(1,2,2);
xslice = median(X(:));
yslice = median(Y(:));
zslice = median(Z(:));
h = slice(X, Y, Z, double(mask), xslice, yslice, zslice)
title('BINARY MASK - SLICE VOLUME', 'FontSize', 16);
set(h, 'EdgeColor','none');
view(-78,22)

Several ellipsoids

If you have more than one ellipsoid, one may use this masking method for each of them, and then combine the resulting masks with &.

Choice of slices and angle

"Correct" is a matter of personal choice. You can either

  • create the unrotated mask and rotate it after (Rotate a 3D array in matlab).
  • create a mask on already rotated ellipsoid.
  • create a mask on a slightly rotated ellipsoid (that gives you the choice of a "correct" slicing), and rotate it further to its final position.
Community
  • 1
  • 1
marsei
  • 7,691
  • 3
  • 32
  • 41
  • Thanks. This is exactly what I was looking for. One more question though. How can I have one unit per voxel? i.e. what if I want my ellipsoid to be 6 x 3 x 3 voxel^3 centered at voxel no (i, j, k) within a surrounding volume of 100 x 100 x 20 voxel^3. – Garima Singh Oct 18 '13 at 09:36
  • One way is to define an ellipsoid with `ellipsoid(0,0,0, 6/2 , 3/2 , 3/2 ,30)`, then a mesh with `meshgrid(-49:50, -49:50, -9:10)` (`-49:50` for 100 elements). As each voxel of the mesh is a unit voxel, then your ellipsoid (x_length = 6) will have 6 positive elements in the x. – marsei Oct 18 '13 at 10:12
  • Thanks a lot for your willingness to help. I tried to write my code based on your suggestions. Please see my edited version of original post for the code. However, I still have two questions: 1) If I have two ellipsoids which I converted to 3D data points based on your suggestion. How can use it to display with correct size of surroundings? 2) Also how do I decide correct "slice" and angles to 3D visualization? Both of these questions are marked by %help: in my code. I would appreciate if you can answer those. Both – Garima Singh Oct 19 '13 at 20:43