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

29 views (last 30 days)
sam l on 20 Aug 2020
Commented: jonas on 22 Aug 2020
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

jonas on 21 Aug 2020
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:
To find the points inside this mesh:

Show 1 older comment
jonas on 22 Aug 2020
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);
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 on 22 Aug 2020
Thanks it works now..
jonas on 22 Aug 2020
My pleasure!