How to determine a point is inside a cube?

Question

Given a cube with 8 vertex in 3D space. How could I determine the myPoint is inside the cube?

cube[0] = (x0, y0, z0);
cube[1] = (x1, y1, z1);
cube[2] = (x2, y2, z2);
cube[3] = (x3, y3, z3);
cube[4] = (x4, y4, z4);
cube[5] = (x5, y5, z5);
cube[6] = (x6, y6, z6);
cube[7] = (x7, y7, z7);

myPoint = (x, y, z);

I found a solution for this on Stack overflow(below). This sample solution returns the indexes of all points outside the cube. I tried to remake it to make the function return the indexes of the points that are inside the cube but failed. I've been sitting on this for 2 days. Is anyone able to help?

import numpy as np

def inside_test(points , cube3d):
    """
    cube3d  =  numpy array of the shape (8,3) with coordinates in the clockwise order. first the bottom plane is considered then the top one.
    points = array of points with shape (N, 3).

    Returns the indices of the points array which are outside the cube3d
    """
    b1,b2,b3,b4,t1,t2,t3,t4 = cube3d

    dir1 = (t1-b1)
    size1 = np.linalg.norm(dir1)
    dir1 = dir1 / size1

    dir2 = (b2-b1)
    size2 = np.linalg.norm(dir2)
    dir2 = dir2 / size2

    dir3 = (b4-b1)
    size3 = np.linalg.norm(dir3)
    dir3 = dir3 / size3

    cube3d_center = (b1 + t3)/2.0

    dir_vec = points - cube3d_center

    res1 = np.where( (np.absolute(np.dot(dir_vec, dir1)) * 2) > size1 )[0]
    res2 = np.where( (np.absolute(np.dot(dir_vec, dir2)) * 2) > size2 )[0]
    res3 = np.where( (np.absolute(np.dot(dir_vec, dir3)) * 2) > size3 )[0]

    return list( set().union(res1, res2, res3) )

When you say it doesn’t work, could you give some examples? Does it crash or does it produce logical errors? Also what have you changed about the original working snippet in your attempt to invert the logic? — lordnoob, Sep 28 '22 at 22:04
+1 to @Woodford but if you wanted to do this ‘properly’ based on my vague understanding of what this code is doing (projecting the vectors between the cube centre and point coordinates along the cube axes and checking the result is greater than the cube side length) it should also work to flip the > sign to < if you want to check if they’re inside. I am not sure about how np sets work and don’t know why the res1,2,3 lines end with [0] though — lordnoob, Sep 28 '22 at 22:11
@lordnob Swapping the > sign to < doesn't help. My problem with this is logical in nature. When changing the sign, the problem is edge values. For example: p(1, 1, 1,) will be recognized as a point outside the cube which is wrong, in that case it helps to add the equal sign <= but then for example p(1.5, 1.5, 1) which is outside the cube will be return as the one that belongs to the cube because of the Z coordinate. res1 = np.where( (np.absolute(np.dot(dir_vec, dir1)) * 2) <= size1)[0] the "where" method will return true in this case. — elfzorik, Sep 28 '22 at 22:23
@elfzorik You need to change the equal sign to <=, but also need to change `union` to `intersection` — AndrzejO, Sep 28 '22 at 22:29

AndrzejO · Accepted Answer · 2022-09-28T22:36:29.913

According to De Morgan's laws:

The complement of the union of two sets is the same as the intersection of their complements

The function returns the outside points, we want the inside – the complement set of points. The function returns the union of three sets. So, in order to get the complement, according to De Morgan's law we need to compute the intersection of the complements of these three sets.

So you need to reverse the unequality signs, and change union to intersection. This should work:

    ...
    res1 = np.where( (np.absolute(np.dot(dir_vec, dir1)) * 2) <= size1 )[0]
    res2 = np.where( (np.absolute(np.dot(dir_vec, dir2)) * 2) <= size2 )[0]
    res3 = np.where( (np.absolute(np.dot(dir_vec, dir3)) * 2) <= size3 )[0]

    return list( set().intersection(res1, res2, res3) )

score 0 · Answer 2 · answered Sep 28 '22 at 23:28

I don't mean to demean the question, but this can be seen as a learning opportunity with regards to understanding how to use the right function for the right job. Like other commenters have noted, the function that should be called in the last line is intersection, not union. To quickly describe the difference and its relevance:

Union: The union of two sets A and B is defined as the set of all elements that belong to A, B, or both.

Intersection: The intersection of two sets A and B is defined as the set of all elements that belong to both A and B at the same time.

Assuming the code to generate the final sets being compared is correct, the problem lies in the last line. Python's intersection function (link here for examples) needs to be called from a set onto another set. What this means, going back to sets A and B, is that an intersection call on these two sets would look like:

A.intersection(B)

There are several things that are important to note here. First, intersection can accept multiple sets to be compared to, so in the case of the original question, the line would look like:

res1.intersection(res2, res3)

The second thing to note is that numpy's where function returns as a tuple, and so the res1, res2, and res3 tuples need to be converted to a set. Thus, the last line would be replaced with the following:

res1 = set(res1)
res2 = set(res2)
res3 = set(res3)
return list(res1.intersection(res2, res3))

score 0 · Answer 3 · answered Sep 29 '22 at 06:41

Included for the general case (based on my answer here). For any n-dimensional convex shape with vertices in any order:

import scipy.spatial

def inside_test2(points, vertices):    
    deln = scipy.spatial.Delaunay(vertices) 
    out_idx = np.flatnonzero(deln.find_simplex(points) + 1)
    return out_idx

Testing (using non-ordered points for cube3d because . . . seriously make real test data as part of your question! Nobody wants to waste time making your unit tests for you.)

import itertools as it

cube3d = np.array(list(it.product([-1,1], repeat=3)))  
# not specified order for original question
points = np.random.rand(10,3) * 1.2

inside_test(points, cube3d))
Out[]: []  #doesn't work like this, because order-dependent

inside_test2(points, cube3d)
Out[]: [1 3 7 8 9]

How to determine a point is inside a cube?

3 Answers3