Find the 3-D grid coordinates points index that are interior of a 3D surface boundary

The aim of the code is to use the given data array , to create a isosurface at value 0.25. The data is scattered 3D array with values for each point in xyz co-ordinate. So we create a 3D grid mesh array , select the grid points inside the data surface boundary , the interpolate the values [data(4,:)] into the 3D grid array and the get the isosurface of the values.

I am stuck in the step to find or search the index of 3-D grid coordinates array for all points which lie inside a volume defined by data points.

Is there any function to find the index of the xx,yy,zz which are inside the surface boundary by data(1,:), data(2,:),data(3,:). That is to find the index of 3-D grid coordinates using a 3D array.

The main code steps are

◾Create boundary of data points [ data(1,:)=x,data(2,:)=y,data(3,:)=z, data(4,:)=value]

◾mesh grid created xx yy zz

◾Missing step - to select index of 3-D grid co-ordinate which is interior of boundary of data points...

◾DATA - 3D array created for holding interpolated values

◾create iso surface using meshgrid DATA

% return index for boundary of all points sequence which is inside volume
j = boundary(data(1,:)',data(2,:)',data(3,:)',0.97);
hold on; trisurf(j,data(1,:)',data(2,:)',data(3,:)',0.97,'FaceColor','red','FaceAlpha',0.1)

% PREPARATION
% ===========
% GENERATE A GRID for x y z co-ordinate
[xx yy zz] = meshgrid(linspace(xx_min,xx_max,20),linspace(yy_min,yy_max,20),linspace(0,zz_max,20));  % replace, 0 1, 10 with range of your values
% GENERATE RANDOM DATA for interplot 3D array value v - grid
DATA = zeros(20,20,20);
%Need to interplot the grid points [xx,yy,zz] inside the volume???
DATA= griddata(data(1,:)',data(2,:)',data(3,:)',data(4,:)',xx,yy,zz); 
% EXTRACT THE  SURFACE
SURF = isosurface(xx,yy,zz,DATA,0.25);  % all x,y,z,v are 3D array 20X20X20
% PLOT THE MASK SURFACE AND DATA
hold on; axis square; axis([0 1 0 1 0 1]); view(3); camlight
p = patch(SURF); isonormals(xx,yy,zz,MASK,p) ;p.FaceColor = 'red'; p.EdgeColor = 'none';

I have used the tsearch using tri boundary points , but this donot work for 3D grid matrix. for example:

% to located interior points in 3D grid 
xyz3D=[xx(:) yy(:) zz(:)]; %  Make 3D grid to a 3D array for using tsearchn
plot3(xyz3D(:,1),xyz3D(:,2),xyz3D(:,3),'*');view(20,30); % plot 3D grid as array
tri = delaunayn([data(1,j)' data(2,j)' data(3,j)']); % Generate delaunay triangulization
[tn p] = tsearchn([data(1,j)' data(2,j)' data(3,j)'],tri,xyz3D) ; % how to use 3D arry for double
tnn=(tn>1);
plot3(xyz3D(tnn,1),xyz3D(tnn,2),xyz3D(tnn,3),'*');view(20,30); % plot inside grid points

temp= griddata(data(1,:)',data(2,:)',data(3,:)',data(15,:)',xyz3D,xyz3D,xyz3D); %interpolate data points
DATA2=reshape(temp,[20,20,20]); % Interpolated data as 3D grid array

I have two things to ask:

• The tnn index [delaunayn] here is not exactly holding all interior points - some are outside the boundary trisuf.
• How to get the volume of the region with a value (say value 0.25) inside boundary surface.

The image show

• the surface boundary with grid points
• the surface with interior points [tnn] from delaunayn

Answer 1

[jonas]

Seems your isosurface is not a convex hull, which is why you end up with some points on the outside. You can try a combination of these two FEX functions/TB's

To find a tight non-convex mesh:

https://se.mathworks.com/matlabcentral/fileexchange/63730-surface-reconstruction-from-scattered-points-cloud

To find the points inside this mesh:

https://se.mathworks.com/matlabcentral/fileexchange/48041-in_polyhedron

[sam l]

Hi thanks a lot for helping out

I tried the two funcions: The triangulation is better now in top regions (to avoid the points getting outside actual surface data ). some issue near bottom as not hollow (no triangles)

However the MyRobustCrust return a NX3 array , while we needed NX4 traingulation to input the Polyhedron . Is there something missing this...

% Robust Crust 
[t,tnorm]=MyRobustCrust([data(1,j)' data(2,j)' data(3,j)']);%t are points id contained in triangles, nx3 array .
trisurf(t,data(1,j)',data(2,j)',data(3,j)','facecolor','c','edgecolor','b')%plot della superficie trattata

% Polyhedron function
TR = triangulation(double(t),[data(1,j)' data(2,j)' data(3,j)']);  %TR:4969X4 triangulation;tet:4969X4 double; X :1456X3 double
[S.faces, S.vertices] = freeBoundary(TR);
in1 = in_polyhedron(S, xyz3D); % 1X1 struct : faces 2174X3 , vertices 1079X3, points (2000X3)

[jonas]

The crust function returns nx3 pts, which does not seem to work with tsearchn. The in_polyhedron works with nx3 face list so that's why I suggested to use it instead.

This is untested code, because I dont have your data:

[t,tnorm]=MyRobustCrust([data(1,j)' data(2,j)' data(3,j)']);%t are points id contained in triangles, nx3 array .
trisurf(t,data(1,j)',data(2,j)',data(3,j)','facecolor','c','edgecolor','b')%plot della superficie trattata
tri.faces = t; 
tri.vertices = [data(1,j)' data(2,j)' data(3,j)';
in1 = in_polyhedron(tri, xyz3D);

This is a working example, using data that you can download from another Q (link)

load seg_info.mat
tri.vertices = seg_info(:,1:3);
tri.faces = MyRobustCrust(tri.vertices);
bounds = [min(tri.vertices);max(tri.vertices)];
[X,Y,Z] = meshgrid(bounds(1,1):0.2:bounds(2,1),...
    bounds(1,2):0.2:bounds(2,2),...
    bounds(1,3):0.2:bounds(2,3));
XYZ = [X(:),Y(:),Z(:)];
in = in_polyhedron(tri,[X(:),Y(:),Z(:)]);
ht = trisurf(tri.faces,tri.vertices(:,1),tri.vertices(:,2),tri.vertices(:,3),'facealpha',0);hold on
hp = scatter3(XYZ(in,1),XYZ(in,2),XYZ(in,3),[],'r');

[sam l]

Thanks it works now..

[jonas]

My pleasure!