What is the algorithm to determine whether a cluster of 3D voxels is a donut shape?

Question

I will start off with some background for this problem. It is an extension of several easier problems. Strangely the final extension leads to it being almost impossible to solve.

1. Given a matrix of integers with 1s, and 0s. 1s representing land and 0s representing ocean count the amount of islands in the matrix which is clusters of 1s separated from each other.

2. Extend the above problem to donuts. Given the same thing as above only count donuts. That means clusters of 1s with one or more holes of "water" or 0s in it. There can be donuts within donuts.

3. Extend the above problem to 3D. When you extend the problem above to 3D it becomes hard enough that I will move the question away from "counting donuts" and more towards "is donut?" So in other words, instead of counting clusters of 1s, now you are told that in this 3D grid space there is one and only one cluster of voxels. That cluster is either a donut or it is not a donut. Which means it has a hole (or several holes) going through it or it does not. Write an algorithm to identify this.

Each question is a more challenging extension of the other. With (1.) being the simplest and (3.) being the hardest.

The first 2 questions are quite straight forward. 1. is a classic interview question. (2.) is an extension of that; simply color all separate bodies of water via flood fill with a separate "color" (aka number) and all "islands" touching 2 or more colors is a donut (the hole touches one body of water, and the outside touches a different body of water).

However the 3rd question is challenging. I cannot come up with a way. My coworkers at 2 jobs... nobody could find a way. So I post it here. The isDonut algorithm question.

from typing import List
#3D_space is gaunteed to have one and only one cluster of pixels in it. 
def isDonut(3d_space: List[List[List[int]]]) -> bool:
    #implement this code

There are many solutions that seem correct but actually fail under specific circumstances. Be careful if you decide to answer.

Edit: For clarity I will define what a voxel donut is:

The above is a voxel donut. A cluster of voxels with a hole going through it. I can only define it in high level english terms. You know it when you see it.

A formal definition of what this is in terms of voxels is the solution to this problem and is described in terms of a programming algorithm. I therefore am unable to describe it, as it's basically my question. Essentially you can think of this question as isomorphic to this one: "what is the formal definition of a voxel donut"?

Edit 2: I'm getting some vague answers and people using advanced math that are hard to understand. Let me put it this way. If you have a straightforward answer you should be able to finish off that python function signature above. You do that the answer is correct, whether or not anyone fully understands your reasoning.

Answer 1

Consider the boundary of your shape: that is, all faces that lie between a filled voxel and an empty one.

Let V be the number of vertices on the boundary, E the number of edges on the boundary, and F the number of faces on the boundary. These are easy to compute; just be careful not to count edges & vertices that belong to multiple boundary faces more than once.

A shape is a donut if and only if (1) the boundary faces are connected, and (2) V-E+F=0.

For more information on this strange and magical second condition, see Euler characteristic.

Answer 2

As you already know #1,#2 then lets focus on #3 (detect if 3D voxel cluster has one ore more holes). After some thinking I revised the original algo a bit:

1. mark border voxels

   so any voxel equal to 1 set to 2 if its neigbors any voxel with 0. After this 0 is empty space, 1 is interior, 2 is surface.

2. use growth fill to create SDR map of your object

   so mark all voxels which are set to 1 to 3 if they neighboring voxel set to 2. Then mark with 4 those which neighbors 3 and so on until no voxel set to 1 is left. This will create something like SDR map (distance to surface).

3. find and count number of local maximums

   for objects without holes there should be just one local max however with holes there would be more of them. In edge case few local max voxels could group to small voxel so count those as one.

Here small C++/OpenGL/VCL example:

//---------------------------------------------------------------------------
#include <vcl.h>
#pragma hdrstop
#include "Unit1.h"
#include "gl_simple.h"
//---------------------------------------------------------------------------
#pragma package(smart_init)
#pragma resource "*.dfm"
TForm1 *Form1;
//---------------------------------------------------------------------------
const int n=40;     // voxel map resolution
int map[n][n][n];   // voxel map color
int max_pos[n][3];  // position of local max
int max_cnt=0;      // number of local max
int max_dis=0;      // number of distinct local max
int pal[32]=        // 0xAABBGGRR
    {
    0x00808080,
    0x00707070,
    0x00606060,
    0x00505050,
    0x00404040,
    0x00303030,
    0x00202020,
    0x00101010,

    0x00800000,
    0x00700000,
    0x00600000,
    0x00500000,
    0x00400000,
    0x00300000,
    0x00200000,
    0x00100000,

    0x00008000,
    0x00007000,
    0x00006000,
    0x00005000,
    0x00004000,
    0x00003000,
    0x00002000,
    0x00001000,

    0x00000080,
    0x00000070,
    0x00000060,
    0x00000050,
    0x00000040,
    0x00000030,
    0x00000020,
    0x00000010,
    };
//---------------------------------------------------------------------------
void TForm1::draw()
    {
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
    glEnable(GL_CULL_FACE);
    glEnable(GL_DEPTH_TEST);
    glEnable(GL_LIGHTING);
    glEnable(GL_LIGHT0);
    glEnable(GL_COLOR_MATERIAL);

    // center the view around map[][][]
    float a;
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
    glTranslatef(0.0,0.0,-80.0);
    a=-0.4*float(n); glTranslatef(a,a,a);
    a=16.0/float(n); glScalef(a,a,a);
//  glRotatef( 15.0,1.0,0.0,0.0);

    // render map[][][] as cubes (very slow old api version for simplicity)
    int x,y,z,i,j;
    for (x=0;x<n;x++)
     for (y=0;y<n;y++)
      for (z=0;z<n;z++)
       if (map[x][y][z])
        {
        glPushMatrix();
        glTranslatef(x+x,y+y,z+z);
        glColor4ubv((BYTE*)&(pal[map[x][y][z]&31]));
        glBegin(GL_QUADS);
        for (i=0;i<3*24;i+=3)
            {
            glNormal3fv(vao_nor+i);
            glVertex3fv(vao_pos+i);
            }
        glEnd();
        glPopMatrix();
        }

    // local max
    glDisable(GL_DEPTH_TEST);
    glColor4f(0.9,0.2,0.1,1.0);
    for (j=0;j<max_cnt;j++)
        {
        x=max_pos[j][0];
        y=max_pos[j][1];
        z=max_pos[j][2];
        glPushMatrix();
        glTranslatef(x+x,y+y,z+z);
        glBegin(GL_QUADS);
        for (i=0;i<3*24;i+=3)
            {
            glNormal3fv(vao_nor+i);
            glVertex3fv(vao_pos+i);
            }
        glEnd();
        glPopMatrix();
        }

    glFlush();
    SwapBuffers(hdc);
    }
//---------------------------------------------------------------------------
__fastcall TForm1::TForm1(TComponent* Owner):TForm(Owner)
    {
    gl_init(Handle);

    // init map[][][]
    int x,y,z,xx,yy,zz,c0,c1,c2,e;
    int x0=n/2,y0=n/2,z0=n/2,rr0=(n/2)-3; rr0*=rr0; // ball
    int x1=n/3,y1=n/2,rr1=(n/5); rr1*=rr1;          // cylinder hole
    for (x=0;x<n;x++)
     for (y=0;y<n;y++)
      for (z=0;z<n;z++)
        {
        // clear map
        map[x][y][z]=0;
        // ball
        xx=x-x0; xx*=xx;
        yy=y-y0; yy*=yy;
        zz=z-z0; zz*=zz;
        if (xx+yy+zz<=rr0) map[x][y][z]=1;
        // hole
        xx=x-x1; xx*=xx;
        yy=y-y1; yy*=yy;
        if (xx+yy<=rr1) map[x][y][z]=0;
        }
    // palette
//  for (x=0;(x<n)&&(x<32);x++) map[x][n-1][n-1]=x;

    // SDR growth fill
    c0=0;   // what to neighbor
    c1=1;   // what to fill
    c2=2;   // recolor to
    for (e=1,c0=0,c1=1,c2=2;e;c0=c2,c2++)
     for (e=0,x=1;x<n-1;x++)
      for (y=1;y<n-1;y++)
       for (z=1;z<n-1;z++)
        if (map[x][y][z]==c1)
         if ((map[x-1][y][z]==c0)
           ||(map[x+1][y][z]==c0)
           ||(map[x][y-1][z]==c0)
           ||(map[x][y+1][z]==c0)
           ||(map[x][y][z-1]==c0)
           ||(map[x][y][z+1]==c0)){ map[x][y][z]=c2; e=1; }

    // find local max
    max_cnt=0;
    max_dis=0;
    for (x=1;x<n-1;x++)
     for (y=1;y<n-1;y++)
      for (z=1;z<n-1;z++)
        {
        // is local max?
        c0=map[x][y][z];
        if (map[x-1][y][z]>=c0) continue;
        if (map[x+1][y][z]>=c0) continue;
        if (map[x][y-1][z]>=c0) continue;
        if (map[x][y+1][z]>=c0) continue;
        if (map[x][y][z-1]>=c0) continue;
        if (map[x][y][z+1]>=c0) continue;
        // is connected to another local max?
        for (e=0;e<max_cnt;e++)
         if (abs(max_pos[e][0]-x)+abs(max_pos[e][1]-y)+abs(max_pos[e][2]-z)==1)
          { e=-1; break; }
        if (e>=0) max_dis++;
        // add position to list
        max_pos[max_cnt][0]=x;
        max_pos[max_cnt][1]=y;
        max_pos[max_cnt][2]=z;
        max_cnt++;
        }

    Caption=AnsiString().sprintf("local max: %i / %i",max_dis,max_cnt);     
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormDestroy(TObject *Sender)
    {
    gl_exit();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormResize(TObject *Sender)
    {
    gl_resize(ClientWidth,ClientHeight);
    draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::FormPaint(TObject *Sender)
    {
    draw();
    }
//---------------------------------------------------------------------------
void __fastcall TForm1::Timer1Timer(TObject *Sender)
    {
    draw();
    }
//---------------------------------------------------------------------------

Just ignore the VCL and OpenGL stuff (they are not important) and focus on the stuff marked with // SDR growth fill and // find local max comments...

Here preview for ball without and with hole:

The local max are rendered without depth testing in orange color and their count (distinct / all) are printed in window Caption...

Answer 3

Edit: This definition does not work. See Rawling's comment below.

Here is an attempt to define a donut, first in a continuous world, through a few observations:

A convex set cannot be a donut. Let S be potential donut, and let conv(S) be the convex hull of S. We define the hole to be H := S \ conv(S). Then S is a donut if H has exactly two disjoint contact surfaces with R^3 \ conv(S). (See below for definitions of "conv()" and "".)

Now, in a discrete voxel world. We can do pretty much the same, except that there are some ambiguities. However, since "donut" is rather informal, they can be resolved according to your personal preferences.

We first need to compute conv(S). There are multiple valid answers here. For example, voxels that partially intersect the continuous conv(S) could be considered part or not part of the discrete convex hull. The construction of H is straightforward, and so are the contact surfaces. The second ambiguity concerns the two disjoint surfaces, specifically what constitutes contiguous voxel faces. A restrictive definition would count 12 neighbors for each voxel face (must have a cube edge in common). But this can be extended to many more if adjacent cube vertices are considered enough.

Note that here I considered that if H is shaped like a Y, then S is not a donut. But this could be up for discussion too.

Disclaimer: not a topologist, my vocabulary may be off. Links to definitions:

Convex hull conv(S): https://en.wikipedia.org/wiki/Convex_hull

S \ conv(S): "set complement" / "Boolean subtraction": https://en.wikipedia.org/wiki/Complement_(set_theory)#Relative_complement / https://en.wikipedia.org/wiki/Constructive_solid_geometry

Edit: Here is an illustration. Yellow: donut. Blue: convex hull. Green: hole. Red: surfaces. Generated with https://evanw.github.io/csg.js/