Is a polygon in a cube - voxelising

[2NywAqkS]

I'm coding a program that converts 3ds files into my own voxel file format how ever part of the process is to test whether any geometry (i.e. polygons) are in a given cubic area. I know the coords of each corner in the cube and I have a list of all the polygons, I just need to know the math to work out whether part of a polygon is in there.

Any help is much appreciated.
Thanks,
Rowan.

[coder777]

Maybe it's time to look at the boost geometry library:

http://www.boost.org/doc/libs/1_50_0/libs/geometry/doc/html/index.html

[2NywAqkS]

I'd much rather work it out mathematically. The way I see it, it's very similar to frustum culling only a cube as opposed to a frustum and there isn't any culling involved. I should probably construct 6 planes and test whether the polygon is in-front of all the planes. However it is this that I don't know how to do.

[Stewbond]

I think I have an answer here (took me a while so sorry for the delay).

#include <iostream> // std::cout, std::endl
#include <cmath>    // sqrt()
typedef double POINT[3];

class PLANE
{
  double A;    // Plane can be defined as
  double B;    // "Ax + By + Cz + D = 0"
  double C;
  double D;
public:
  void init(POINT p1, POINT p2, POINT p3); // Create a plane from 3 points
  bool isAbove(POINT p); // Tests whether a point is above, in front, or to the right of a plane
};

class POLYGON6 // 6 sided polygon (a cube will be our test object)
{
  PLANE plane[6];
public:
  // Defined with 8 points.  Could be less, but I'm just being casual right now.
  POLYGON6(POINT p1, POINT p2, POINT p3, POINT p4, POINT p5, POINT p6, POINT p7, POINT p8);
  bool isInPoly(POINT p);
};

////////////////////////////////////////////////////
///////////// PLANE member functions ///////////////
////////////////////////////////////////////////////

void PLANE::init(POINT p1, POINT p2, POINT p3)
{
  A = p1[1]*(p2[2] - p3[2]) + p2[1]*(p3[2] - p1[2]) + p3[1]*(p1[2] - p2[2]); 
  B = p1[2]*(p2[0] - p3[0]) + p2[2]*(p3[0] - p1[0]) + p3[2]*(p1[0] - p2[0]);
  C = p1[0]*(p2[1] - p3[1]) + p2[0]*(p3[1] - p1[1]) + p3[0]*(p1[1] - p2[1]); 
  D = -1 * (p1[0]*(p2[1]*p3[2] - p3[1]*p2[2]) + p2[0]*(p3[1]*p1[2] - p1[1]*p3[2]) + p3[0]*(p1[1]*p2[2] - p2[1]*p1[2]));
}

bool PLANE::isAbove(POINT p)
{
  // Method:
  // Find the Z-intercept (D) for the parallel plane that the test point lies on.
  // This can be done by using the same A, B and C coefficients
  // If the parallel plane's D is below this one, it's considered below.
  // Take the absolute value of D to avoid problems with upside-down planes.
  double NewD = -1*(A*p[0] + B*p[1] + C*p[2]); 
  return !(abs(NewD) < abs(D));                          
}

////////////////////////////////////////////////////
///////////// POLYGON6 member functions ////////////
////////////////////////////////////////////////////

POLYGON6::POLYGON6(POINT p1, POINT p2, POINT p3, POINT p4, POINT p5, POINT p6, POINT p7, POINT p8)
{
  plane[0].init(p1, p2, p3); // Front 3 
  plane[1].init(p1, p3, p4); 
  plane[2].init(p1, p4, p2);
  plane[3].init(p8, p7, p5); // Back 3
  plane[4].init(p8, p5, p6);
  plane[5].init(p8, p6, p7);
}

bool POLYGON6::isInPoly(POINT p)
{
  int above = 0;
  for (int i = 0; i < 6; ++i)
    if (plane[i].isAbove(p))	above++;
    else                     	above--;

  if (above) return false; // return false if non-zero.  It's not balanced between all planes
  return true; // return true if zero, it's positive of just as many planes as it is negative.
}

////////////////////////////////////////////////////
///////////// main : ie test function //////////////
////////////////////////////////////////////////////

int main()
{
  POINT p1 = {0,0,0};
  POINT p2 = {1,0,0};
  POINT p3 = {0,1,0};
  POINT p4 = {0,0,1};
  POINT p5 = {1,1,0};
  POINT p6 = {1,0,1};
  POINT p7 = {0,1,1};
  POINT p8 = {1,1,1};

  POLYGON6 poly(p1, p2, p3, p4, p5, p6, p7, p8);

  POINT test1 = {0.5, 0.5, 0.5}; // should return true
  POINT test2 = {2.0, 2.0, 2.0}; // should return false

  std::cout << (poly.isInPoly(test1)? "yes" : "no") << std::endl;
  std::cout << (poly.isInPoly(test2)? "yes" : "no") << std::endl;
}

EDIT: the original code didn't work because I forgot to apply the right-hand rule in the POLYGON6 constructor and didn't take the absolute value of D in bool PLANE::isAbove(POINT p).

[Stewbond]

I got it running (I turned that PLANE constructor into a member function). But test1 returns false and test2 returns true (oppsite). There is a problem with my method in bool PLANE::isAbove(POINT p).

[2NywAqkS]

Wow! thanks Stewbond this is really going to help

[Stewbond]

You're very welcome. I was actually working on some linear algebra stuff at work today (not for graphics though) so my mind was in this frameset anyways.

I've updated the code above to make it a little more robust so that the order of points defined in the POLYGON6 constructor doesn't matter (Right-hand-rule isn't crucial anymore). I did that by taking the abs() of D in bool PLANE::isAbove(POINT p).

[mik2718]

Looks like this code will test if a point is inside a cube but will it test if a polygon intersects a cube?

For example, just an edge of the polygon might pass through the cube, in this case the polygon intersects the cube but no polygon vertex is inside the cube.

The polygon can even intersect the cube when no polygon vertices are in the cube and no edges of the polygon intersect the cube.

I have thought about whether there is a simple way of performing this test when the cube is axis aligned but in all the ways I have though of there are some cases where the simple method fails.

I think you will need to use the "method of separating axes"

http://www.cplusplus.com/forum/general/74796/

http://www.codezealot.org/archives/55

so to do this accurately is not simple, it may be the case that a simpler approximation is OK for your program

[helios]

Maybe the rotating calipers algorithm? Basically what you do successively crop the polygon using infinite extensions of each of the sides of the cube. If at the end of the process there's any vertices left at all, it's because the polygon and the cube did indeed intersect.

[Stewbond]

My next thought is to include boundry conditions in the planes by using all 4 points to defined each one. Then determine the intercept point between the vector (ie the edge of your test polygon) and the plane. Then determine if the point is within the boundary conditions.