1

I'm trying to construct a plane out of three points in 3D. I want to use projective geometry to achieve this.

As far as I know, one can "simply" solve the following to find a plane:

A * x = 0 ,where
A is a 3x4 Matrix - each row being one of the points (x,y,z,1)
x is the plane I want to find

I know that I need to have a constrain. Therefore I want to set x(3) = 1. Can someone please point me to the right method to use?

So far I have the following code:

Eigen::Vector4f p1(0,0,1,1);
Eigen::Vector4f p2(1,0,0,1);
Eigen::Vector4f p3(0,1,0,1);

Eigen::Matrix<float,3,4> A;
A << p1.transpose(), p2.transpose(), p3.transpose();

// Throws compile error
// Eigen::Vector4f Plane = A.jacobiSvd(ComputeThinU | ComputeThinV).solve(Vector4f::Zero()); 

//throws runtime error (row-number do not match)
// Eigen::Vector4f Plane = A.fullPivHouseholderQr().solce(Eigen::Vector4f::Zero());
NewTech
  • 241
  • 1
  • 3
  • 14
  • I would go for a non-matrix method, see https://en.wikipedia.org/wiki/Plane_(geometry)#Method_3 . – Shaana Mar 15 '17 at 10:13
  • The idea is to use projective geometry to have an easy way of testing whether other points belong to this plane by simply computing `Plane.transpose() * otherPoint` – NewTech Mar 15 '17 at 10:17

1 Answers1

3

A 3x4 matrix multiplied by a 4 row vector will give you a 3 row vector. Thus you have to solve for a Vector3f::Zero(). Also, for fixed size matrices you need to compute the full U and V. The final line looks like this:

Vector4f Plane = A.jacobiSvd(ComputeFullU | ComputeFullV).solve(Vector3f::Zero());

Eidt Since this equation system is not fully defined, it might give you the trivial solution of (0,0,0,0). You can solve that by constraining the length of the resulting vector by expanding the matrix to a 4x4, solving for (0,0,0,1) and scaling the result by x(3):

Eigen::Vector4f p1(0,0,1,1);
Eigen::Vector4f p2(1,0,0,1);
Eigen::Vector4f p3(0,1,0,1);
Eigen::Vector4f p4(1,1,1,1);

Eigen::Matrix<float,4,4> A;
A << p1.transpose(), p2.transpose(), p3.transpose(), p4.transpose();

// Throws compile error
Eigen::Vector4f Plane = A.jacobiSvd(Eigen::ComputeFullU | Eigen::ComputeFullV).solve(Vector4f::Unit(3)); 
Plane /= Plane(3);

This will give you the desired solution of (-1, -1, -1, 1).

maddin45
  • 737
  • 7
  • 17
  • Thanks. Now the code compiles but stops at the following line: `eigen_assert(EIGEN_IMPLIES(m_computeThinU || m_computeThinV, MatrixType::ColsAtCompileTime==Dynamic) && "JacobiSVD: thin U and V are only available when your matrix has a dynamic number of columns.")` inside the "JacobiSVD.h" – NewTech Mar 15 '17 at 10:10
  • @NewTech I updated my answer. The problem is that `ComputeThinU` and `ComputeThinV` are not available for fixed-size matrices like `Matrix`. – maddin45 Mar 15 '17 at 10:12
  • Ah. Ok now the code runs perfectly. But only giving the trival solution of (0,0,0,0). I would expect to get (-1, -1, -1, 1). Is there some method to set a constrain? – NewTech Mar 15 '17 at 10:15
  • 1
    @NewTech Constrain the length of your vector to be non-zero, by using a non-zero fourth matrix-row and a fourth entry in the vector you solve for with non-zero value. – maddin45 Mar 15 '17 at 10:28
  • I previously voted up this, but I found input case, when this SVD solution fails. Consider triangle with vertexes (0,0,0), (1,1,0), and (0,0,1) (see https://math.stackexchange.com/questions/305642/how-to-find-surface-normal-of-a-triangle/307155#comment666478_305914). The result using this approach with eigen v 3.3.7 is (-0.408248, -0.408248, -0.816497), which is not correct. Correct result is (-0.707107, 0.707107, 0). I think that better way here to use cross-product method, or this formula: https://math.stackexchange.com/a/307155/852232 – Mykyta Kozlov Jan 23 '22 at 15:09