To show the solution I'll work with a torus, which is very similar to a donut. The first thing to do will be to create a Torus as reference. The mathematical equation for a torus is given by:
Creating the torus can be done via the following equation:
def create_torus(xlim, ylim, zlim, step, R, r, epsilon):
"""
:param xlim: limits of the space in the x axis, tuple
:param ylim: limits of the space in the y axis, tuple
:param zlim: limits of the space in the z axis, tuple
:param step: precision of the 3d space
:param R: distance from the center of the tube to the center of the torus, int
:param r: radius of the tube, int
:param epsilon: since this is discrete, stating an epsilon for numerical calculations
:return: 3D space with voxels indicating the shape of a torus
Torus equation in cartesian coordinates:
"""
# ==================================================================================================================
# Local variables
# ==================================================================================================================
x_vec = np.arange(xlim[0], xlim[1], step)[:, None]
y_vec = np.arange(ylim[0], ylim[1], step)[None, :]
z_vec = np.arange(zlim[0], zlim[1], step)[None, None, :]
space = np.zeros([len(x_vec), len(y_vec[0, :]), len(z_vec[0, 0, :])])
# ==================================================================================================================
# Computing torus function
# ==================================================================================================================
torus_func = ((np.sqrt(x_vec**2 + y_vec**2) - R)**2)[:, :, None] + z_vec**2
# ==================================================================================================================
# Thresholding and allocating voxel location
# ==================================================================================================================
threshold = np.abs(torus_func - r**2) < epsilon
space[threshold] = 1
return space
an Example of such a torus is given below (R = 4, r = 1, epsilon = 0.1):
Now That we have the input created, the solution is given as follows:
Solution Explanation
A common engineering method will be to break a problem into a set of smaller problems which is easier to solve. The 3D problem can be reduced into 3D via slicing the 3D space into plains. Assuming that the donut is symmetric w.r.t the z axis(can be expanded to align the donut as well but this is not the topic), one could slice the 3D space into XY plains, and check for a set of 2d "islands" with specific characteristics. An XY slice of a torus will have either 1, 2 or 4 "islands". ignoring the case of 1 island, if in any XY plain, one would fine ONLY 2 or 4 "islands" then the 3D "island" is a Donut.
Let me note that any method of solving the 3D case can be accepted, since in 2D slice of a 3D torus in the XY plain, there will be 1, 2 or 4 circles meaning the the "floodfill" method of solving the 2D case also works here.
Anyway, the solution I made for the function is as follows:
def isDonut(space):
"""
:param space: Voxel space where '0' is nothing and '1' is somethings
:return: function returns boolean indicating if the voxel space contains a donut (torus)
Assumptions:
1. only 1 "island" in the voxel space
2. if there is a donut, it is aligned with the z axis (can be expanded to align the torus but for simplification...)
The method of doing this is to reduce the problem into another, easier problem. Function slices the 3D space into
several XY planes, checking if there are exactly 2 closed contours or 4 closed contours. If there are, the 3d "island" is a torus
"""
z_length = space.shape[2]
# ==================================================================================================================
# Operating on slices
# ==================================================================================================================
for ii in range(z_length):
# ----------------------------------------------------------------------------------------------------------
# Extracting plain
# ----------------------------------------------------------------------------------------------------------
plain = space[:, :, ii]
if np.max(plain) == 0:
continue
# ----------------------------------------------------------------------------------------------------------
# Getting contours
# ----------------------------------------------------------------------------------------------------------
plain_uint8 = (plain * 255).astype(np.uint8)
contours, hierarchies = cv2.findContours(plain_uint8, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) == 2 or len(contours) == 4:
return True
return False
Note that this solution does not solve for the case of a cylinder with a "dent" inside, since a 2D slice might meet the conditions mentioned above, but for most of the cases it will be sufficient. To extend the solution above to the more general case one would see that each 2D slice should contain at least 1 2d Donut, meaning that there will be '0' inside the contour and '0' outside the contour.
More general solution
A more general approach will be to find the center of mass Y coordinate for each of the contours in the XY plain, and check the number of non-overlapping contours in the XZ plain. If there is a Donut, there will be 2 non-overlapping contours in this plain. A true full solution will be to check all XZ plains but this relaxation of checking the center of mass will usually be enough. The implantation if this is a follows:
def isDonut(space):
"""
:param space: Voxel space where '0' is nothing and '1' is somethings
:return: function returns boolean indicatin if the voxel space contains a donut (torus)
Assumptions:
1. only 1 "island" in the voxel space
2. if there is a donut, it is aligned with the z axis (can be expanded to align the torus but for simplification...)
The method of doing this is to reduce the problem into another, easier problem. Function slices the 3D space into
several XY planes, checking if there are exactly 2 closed contours or 4 closed contours. If there are, the 3d "island" is a torus
"""
z_length = space.shape[2]
# ==================================================================================================================
# Operating on slices
# ==================================================================================================================
for ii in range(z_length):
# ----------------------------------------------------------------------------------------------------------
# Extracting plain
# ----------------------------------------------------------------------------------------------------------
xy_plain = space[:, :, ii]
if np.max(xy_plain) == 0:
continue
# ----------------------------------------------------------------------------------------------------------
# Getting contours
# ----------------------------------------------------------------------------------------------------------
xy_plain_uint8 = (xy_plain * 255).astype(np.uint8)
contours, hierarchies = cv2.findContours(xy_plain_uint8, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours) == 2 or len(contours) == 4:
for c in contours:
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Getting the center of mass for the Y coordinate
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
M = cv2.moments(c)
cY = int(M["m01"] // M["m00"])
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Getting the XZ slice of the space
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
xz_plain = space[:, cY, :]
xz_plain_uint8 = (xz_plain * 255).astype(np.uint8)
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Checking if there are two non-overlapping contours in the XZ plain according to hierarchies
# [Next, Previous, First_Child, Parent]
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
contours_temp, hierarchies_temp = cv2.findContours(xz_plain_uint8, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)
if len(contours_temp) > 1:
for jj in range(hierarchies_temp.shape[1]): # number of contours
if all(hierarchies_temp[0, jj][2:] == -1):
return True
return False
