# How to find the concave hull for a cloud of points in 3D space?

The existing algorithm for convex hull is not able to capture the feature for a set of 3D points. Moreover, I found few mathematic tools have this function to obtain the concave hull and their responding points.

Given the data of spheres: ...

Any idea?

Update 1:

I found a 2D algorithm, it works fine depending on the threshold value. However, I need the 3D algorithm. • What software are you using, or are you after a mathematical principle? I found using the vector dot product to find the left-leftest the best. Apr 24, 2015 at 0:36
• I use Mathematics in Windows and Octave in Linux. Both coeds don't provide available function for this issue. I am looking for a algorithms based on above coeds.
– KOF
Apr 24, 2015 at 0:42
• en.wikipedia.org/wiki/Convex_hull_algorithms, that's where I started. The simplest (I found) is to use the vector product to find the left of the left then iterate. Apr 24, 2015 at 0:46
• Hi @KOF what algorithm are you using to produce the 2D version? I'm looking for something similar in gis.stackexchange.com/questions/152175/… Jun 25, 2015 at 8:33

A convex hull is unique, whereas there are many possible concave hulls. So you cannot say "the concave hull" but "a concave hull".

There is possibly a minimal volume concave hull, but this is not the case on the example you shown. It is also possible to define various criteria, such as the minimal acceptable concave edge angle, for avoiding deep trenches or pits in the obtained hull.

All hulls on the following picture are valid, depending on the level of "tightness" you are looking for To find a "concave hull" around a set of 3D points, I found that using the marching cube algorithm for volumetric data works best. Here is an example using Python. To run it, you first need to transform your cloud of 3D points into a volumetric dataset. Then, you obtain the points at the border of your cloud and the faces that can be used to make a mesh that can be plotted.

``````from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
import numpy as np
from skimage import measure

# In my case, I have position of voxels in a MRI volume
points = np.array([[-0.025,  0.015, -0.04 ],
[-0.02 ,  0.015, -0.04 ],
[-0.03 ,  0.015, -0.035],
[-0.025,  0.015, -0.035],
[-0.02 ,  0.015, -0.035],
[-0.015,  0.015, -0.035],
[-0.025,  0.01 , -0.03 ],
[-0.025,  0.015, -0.03 ],
[-0.02 ,  0.015, -0.03 ],
[-0.015,  0.015, -0.03 ]])

# And I have a transformation matrix that took the voxels indices
# to these 3D coordinates
trans_mat = np.array([[ 0.005,  0.   ,  0.   , -0.06 ],
[ 0.   ,  0.005,  0.   , -0.075],
[ 0.   ,  0.   ,  0.005, -0.075],
[ 0.   ,  0.   ,  0.   ,  1.   ]])

# The shape of my MRI voxel volume
volume_shape = (25, 32, 26)

# From 3D points to voxel indices
points_4d = np.hstack((points, np.ones((points.shape, 1))))
trans_points = np.linalg.inv(trans_mat)  @  points_4d.T

voxels = np.zeros(volume_shape)
for ix, iy, iz in np.round(trans_points[:3, :]).T.astype(int):
voxels[ix, iy, iz] = 1

vertices, faces = measure.marching_cubes_lewiner(voxels, spacing=(1, 1, 1))[:2]

fig = plt.figure()