19 replies to this topic

[Vexator]

I have an arbitrary number of points which are all located on a plane in 3D space. I need to find the four corner points of the rectangle which tightly encloses all points. how can i do that?

thank you!

[pantaloons]

I have an arbitrary number of points which are all located on a plane in 3D space. I need to find the four corner points of the rectangle which tightly encloses all points. how can i do that?

thank you!

"Tightly encloses" is ill-defined. If you mean the rectangle with minimal area, this is called the minimal enclosing rectangle problem and can be solved in linear time. See http://cgm.cs.mcgill.ca/~orm/maer.html

For other definitions of "tightly encloses" it may be significantly easier or significantly harder.

[Vexator]

thanks for the link! I'm having trouble computing the extreme points, though. If it was in 2D I'd know how to do it but on an arbitrary plane in 3D, I'm not sure

[Kryzon]

You need to view those points in an axis-aligned space.

Multiply each point's transform by the inverse of the plane's transform.

[Vexator]

Multiply each point's transform by the inverse of the plane's transform.

How exaclty would I do that? I store my planes as normal plus distance to origin.

[Álvaro]

You should use a change of coordinates. Take your normal vector, find two unit vectors that are perpendicular to it and to each other and use that as your base. In these new coordinates, your plane is z = constant and you can deal with the x and y coordinates as if you were in 2D.

[Kryzon]

How exaclty would I do that? I store my planes as normal plus distance to origin.

Your 'normal' describes a rotation, and the 'distance to origin' describes a position. With these two elements you can build a transformation matrix that represents your plane's space.
If you want to view your points in an axis-aligned manner, multiply their XYZ coordinates by the inverse of this matrix. This can be viewed as aligning the plane with your world's origin - you are removing any transformations from this plane, and bringing the points with it. Any distance these points had towards the plane's origin, now they will have the same towards the world's origin.

After you've multiplied them all with the plane's inverse matrix and stored their new position results, find out the enclosing rectangle with the algorithm of your choice - run it on these stored positions you have. This should give you the 4 points for the enclosing rectangle. Multiply each of these 4 points' XYZ coordinates by the plane's matrix to bring them back into the plane's space, which would be akin to having the rectangle get aligned with that arbitrary plane, whichever rotation it might have. This is your transformed rectangle.

For me the most difficult part here would be to extract a rotation from the plane's normal vector - I don't have any resources at hand where this is shown easily. I use an API with utilitary functions for this.

[quasar3d]

I have an arbitrary number of points which are all located on a plane in 3D space. I need to find the four corner points of the rectangle which tightly encloses all points. how can i do that?

thank you!

"Tightly encloses" is ill-defined. If you mean the rectangle with minimal area, this is called the minimal enclosing rectangle problem and can be solved in linear time. See http://cgm.cs.mcgill.ca/~orm/maer.html

For other definitions of "tightly encloses" it may be significantly easier or significantly harder.

Are you sure it can be solved in linear time? The link you posted only seems to solve it in linear time for convex polygons. I can't really think of quick way to prove that it's not possible in linear time for an arbitary set of points, but I have a suspicion that its not possible.

[solenoidz]

Is this kind of approach would work (linear time) ?
Pseudo code follows.
Basically, you try to expand the rectangle, starting with zero dimentions.

struct rectangle
{
vector2d upperLeft ;
vector2d lowerRight ;
} ;
rectangle boundRec;
boundRec.upperLeft  = boundRec.lowerRight  = points[0] ;

for (int i = 0 ;  i < points.size() ; i++)
  {
    if(points[i].X > boundRec.lowerRight.X)  boundRec.lowerRight.X = points[i].X ;
   if(points[i].Y > boundRec.lowerRight.Y)  boundRec.lowerRight.Y = points[i].Y ;
  if(points[i].X > boundRec.upperLeft  .X)  boundRec.upperLeft  .X = points[i].X ;
  if(points[i].Y > boundRec.upperLeft  .Y)  boundRec.upperLeft  .Y = points[i].Y ;
  }

[quasar3d]

Is this kind of approach would work (linear time) ?
Pseudo code follows.
Basically, you try to expand the rectangle, starting with zero dimentions.

struct rectangle
{
vector2d upperLeft ;
vector2d lowerRight ;
} ;
rectangle boundRec;
boundRec.upperLeft  = boundRec.lowerRight  = points[0] ;

for (int i = 0 ;  i < points.size() ; i++)
  {
    if(points[i].X > boundRec.lowerRight.X)  boundRec.lowerRight.X = points[i].X ;
   if(points[i].Y > boundRec.lowerRight.Y)  boundRec.lowerRight.Y = points[i].Y ;
  if(points[i].X > boundRec.upperLeft  .X)  boundRec.upperLeft  .X = points[i].X ;
  if(points[i].Y > boundRec.upperLeft  .Y)  boundRec.upperLeft  .Y = points[i].Y ;
  }

That always gives axis aligned rectangle. The minimal one will often not be axis aligned.

[solenoidz]

Yes, Kryzon gave an explanation how to switch betwen these two spaces.

[quasar3d]

Yes, Kryzon gave an explanation how to switch betwen these two spaces.

Kryzon gave an explanation how to switch a space which has the plane with the points aligned on the x,y plane. This doesn't turn the bounding rect into an axis aligned one though.

[solenoidz]

Yes, Kryzon gave an explanation how to switch betwen these two spaces.

Kryzon gave an explanation how to switch a space which has the plane with the points aligned on the x,y plane. This doesn't turn the bounding rect into an axis aligned one though.

Now I can see what you mean - thanks. Reputation ++

[pantaloons]

I have an arbitrary number of points which are all located on a plane in 3D space. I need to find the four corner points of the rectangle which tightly encloses all points. how can i do that?

thank you!

"Tightly encloses" is ill-defined. If you mean the rectangle with minimal area, this is called the minimal enclosing rectangle problem and can be solved in linear time. See http://cgm.cs.mcgill.ca/~orm/maer.html

For other definitions of "tightly encloses" it may be significantly easier or significantly harder.

Are you sure it can be solved in linear time? The link you posted only seems to solve it in linear time for convex polygons. I can't really think of quick way to prove that it's not possible in linear time for an arbitary set of points, but I have a suspicion that its not possible.

You are right, its loglinear for nonconvex polygons. Any algorithm which enumerates the convex hull trivially cannot be linear time.

[Vexator]

thanks for the input everybody!
here's what i've come up with. results are not as expected, though.

original portal vertices:

[0]	{ x=-1.5885940 y=-1.2763360 z=4.7325292 }
[1]    { x=-1.5885930 y=1.8201370 z=4.7325320 }
[2]    { x=-3.3042150 y=1.8201380 z=3.7420130 }
[3]    {x=-3.3042150 y=-0.56992221 z=3.7420115 }
[4]    { x=-2.6216710 y=-1.2763354 z=4.1360798 }

vertices transformed by plane matrix inverse:

[0]	{  x=6.1884351 y=-1.2763298 z=-7.5414600 }
[1]	{  x=6.1884346 y=1.8201432 z=-7.5414615 }
[2]	{  x=7.1789441 y=1.8201436 z=-5.8258338 }
[3]	{  x=7.1789441 y=-0.56991661 z=-5.8258338 }
[4]	{  x=6.7848792 y=-1.2763295 z=-6.5083799 }

computing a matrix from the portal's plane and multiplying its vertices with the inverse matrix:

const vec3 v1 = portal->getVertex(0);
const vec3 v2 = portal->getVertex(1);
const vec3 v3 = portal->getVertex(2);

const vec4 portalPlane = computePlane( v1, v2, v3 );

vec3 normal = vec3( portalPlane );

vec3 up( 0, 1, 0 );

vec3 f( glm::normalize(normal) );
vec3 s( glm::normalize(glm::cross(f, up)) );

vec3 u( glm::normalize(glm::cross(s, f)) );

mat4 planeMatrix(
	s[0], 	u[0], 	-f[0], 	0.0,
	s[1], 	u[1], 	-f[1], 	0.0,
	s[2], 	u[2], 	-f[2], 	0.0,
	0.0f, 	0.0f, 	0.0f,      1.0f);

vec3 v( normal*(-portalPlane.w) );

for (unsigned i = 0; i < 3; ++i)
{
	double tmp = v[i];

	if (tmp == 0)
		continue;

	planeMatrix[3][0] += tmp*planeMatrix[i][0];
	planeMatrix[3][1] += tmp*planeMatrix[i][1];
	planeMatrix[3][2] += tmp*planeMatrix[i][2];
	planeMatrix[3][3] += tmp*planeMatrix[i][3];
}

mat4 planeMatrixInverse = glm::inverse( planeMatrix );

vector<vec3> newVertices( portal->getVertexCount() );

for( uint i = 0; i < portal->getVertexCount(); i++ )
{
	const vec3 &position = portal->getVertex(i);

	newVertices[i] = vec3( planeMatrixInverse * vec4(position, 1.0f) );
}

what do you think?

[solenoidz]

Can't you fill the translation part of the matrix directly. is you roation matrix orthogonal, can you transform and render a mesh with it to see if it looks right.
Also, are you sure vec3 normal = vec3 (portalPlane) returns normalized normal, because you use it later on normal *(-portalPlane.w)

[quasar3d]

I'm not that familiar with glm, but I just took a quick look at the source, and it seems that the mat4 constructor takes its parameters in column major order, so the matrix you create is in fact the transpose of the matrix you want.

[Vexator]

thank you! i adopted your corrections:

const vec3 v1 = portal->getVertex(0);
const vec3 v2 = portal->getVertex(1);
const vec3 v3 = portal->getVertex(2);

const vec4 portalPlane = computePlane( v1, v2, v3 );

// FIX: normalize normal
vec3 normal = glm::normalize( vec3(portalPlane) );

vec3 up( 0, 1, 0 );

vec3 f( glm::normalize(normal) );
vec3 s( glm::normalize(glm::cross(f, up)) );

vec3 u( glm::normalize(glm::cross(s, f)) );

mat4 planeMatrix(
	s[0], 	u[0], 	-f[0], 	0.0,
	s[1], 	u[1], 	-f[1], 	0.0,
	s[2], 	u[2], 	-f[2], 	0.0,
	0.0f, 	0.0f, 	0.0f,      1.0f);

vec3 v( normal*(-portalPlane.w) );

for (unsigned i = 0; i < 3; ++i)
{
	double tmp = v[i];

	if (tmp == 0)
		continue;

	planeMatrix[3][0] += tmp*planeMatrix[i][0];
	planeMatrix[3][1] += tmp*planeMatrix[i][1];
	planeMatrix[3][2] += tmp*planeMatrix[i][2];
	planeMatrix[3][3] += tmp*planeMatrix[i][3];
}

// FIX: transpose matrix
mat4 planeMatrixInverse = glm::inverse( glm::transpose(planeMatrix) );

vector<vec3> newVertices( portal->getVertexCount() );

for( uint i = 0; i < portal->getVertexCount(); i++ )
{
	const vec3 &position = portal->getVertex(i);

	newVertices[i] = vec3( planeMatrixInverse * vec4(position, 1.0f) );
}

portal = shared_ptr<Portal>( new Portal(newVertices) );

now the resulting vertices are aligned on a 2D plane (see z-coordinate). however, they are still not at the origin:

[0]	{x=-0.99051094 y=-1.2763327 z=-4.8927865 }
[1]	{x=-0.99051315 y=1.8201402 z=-4.8927865 }
[2]	{x=0.99051857 y=1.8201412 z=-4.8927865 }
[3]	{x=0.99051929 y=-0.56991905 z=-4.8927870 }
[4]	{x=0.20238473 y=-1.2763321 z=-4.8927870 }

[quasar3d]

The matrix you create is the matrix that maps the xz plane with z = 0 to your original plane (the inverse of the matrix you eventually transform your vertices with). The rotation part of this matrix yields a plane through the origin, parallel to portalPlane. Then to let it map to portalPlane, you have to shift this plane by -portalPlane.w into the direction of its normal. This means that you can simply set the last column of planeMatrix to normal * -portalPlane.w. There's no need to transform this by the rotation part of your matrix, as you are doing now.

Also, right now you do the transpose after you add the translation part, even though you do use the last column instead of the last row for it. To fix this, either add the rotation to the last row, or do the transpose before you add the translation.

[Vexator]

thank you! I adopted your suggestions and it seems to work just fine. I was just surprised that I had to change the bounding vertices' order; otherwise, the clip planes' normals were facing the wrong way.

here's the final code:

const vec3 v1 = portal->getVertex(0);
const vec3 v2 = portal->getVertex(1);
const vec3 v3 = portal->getVertex(2);

// compute plane reflecting portal's orientation
const vec4 portalPlane = computePlane( v1, v2, v3 );

// compute plane matrix

vec3 normal = vec3( portalPlane );
vec3 position( normal * (-portalPlane.w) );

vec3 up( 0, 1, 0 );

vec3 xAxis( glm::normalize(glm::cross(normal, up)) );
vec3 yAxis( glm::normalize(glm::cross(xAxis, normal)) );

mat4 planeMatrix;

planeMatrix[0] = vec4(  xAxis, 0.0f );
planeMatrix[1] = vec4(  yAxis, 0.0f );
planeMatrix[2] = vec4( -normal, 0.0f );
planeMatrix[3] = vec4(  position, 1.0f );

mat4 planeMatrixInverse = glm::inverse( planeMatrix );

// un-project portal vertices

uint vertexCount = portal->getVertexCount();

vector<vec3> projectedVertices( vertexCount );

for( uint i = 0; i < vertexCount; i++ )
{
	const vec3 &position = portal->getVertex(i);

	projectedVertices[i] = vec3( planeMatrixInverse * vec4(position, 1.0f) );
}

// compute bounding vertices
// from projected portal vertices
					
vec3 minimum = projectedVertices.at(0);
vec3 maximum = projectedVertices.at(0);

for( uint i = 0; i < projectedVertices.size(); i++ )
{
	const vec3 &vertex = projectedVertices.at(i);

	if( vertex.x < minimum.x )	minimum.x = vertex.x;
	if( vertex.y < minimum.y )	minimum.y = vertex.y;
	if( vertex.x > maximum.x )	maximum.x = vertex.x;
	if( vertex.y > maximum.y )	maximum.y = vertex.y;
}

// re-project bounding vertices

vector<vec3> boundingVertices( 4 );

boundingVertices[3] = vec3( planeMatrix * vec4(maximum.x, minimum.y, 0, 1) );
boundingVertices[2] = vec3( planeMatrix * vec4(maximum.x, maximum.y, 0, 1) );
boundingVertices[1] = vec3( planeMatrix * vec4(minimum.x, maximum.y, 0, 1) );
boundingVertices[0] = vec3( planeMatrix * vec4(minimum.x, minimum.y, 0, 1) );

// create new portal from bounding vertices
portal = shared_ptr<Portal>( new Portal(boundingVertices) );