Best plane fit from points

[f-frhs]

I added Plane.CreateFittedPlaneFrom(points) method with some unit tests for it.

I would like to discuss the following topic:

• Were testcases I added sufficient?
• implementation
  • which is nice to evaluate?
    a. (z)axis-aligned distance
    b. orthogonal distance

[jkalias]

Thank you for this PR, it seems interesting. Let me know if my comments are not clear enough.

[jkalias]

I would use more descriptive variable names, like throughPoint instead of just c.

[f-frhs]

will do, thanks.

[jkalias]

This is also called the scatter matrix

[f-frhs]

If I understood your comment correctly, it seems to me that the scatter matrix means a different matrix from the matrix A written in the code.

To say that these 2 matrices are different, I checked their sizes(=rows count and columns count).
The scatter matrix is an 3x3 matrix,
and the matrix A is an nx3 matrix,
where $n$ is the count of points.

I agree with you to give the variable a descriptive name.
Do you have any ideas?

[jkalias]

My apologies, you are right. The scatter matrix (M in the paper) would be A^T * A in your notation. Perhaps we can call this relativePointMatrix, since this is exactly what it is.

[jkalias]

What I would do in this case, is to first allocate the necessary matrix space (new DenseMatrix(rows, columns)), and then write a single for-each to populate each row of it. This way we could easily parallelize the loop in case we need it, if for example we have a LOT of points.

[f-frhs]

will do, thanks.
I think that the parallelization is out of scope for me.
Thanks again for your point of view.

[f-frhs]

Does this fix follow the policy of Math.NET Spatial?
I'm afraid the case where you maintainers respect RAII(Resource acquisition is initialization) policy and I would break it.

[jkalias]

I think that the parallelization is out of scope for me.

Sure, I wasn't expecting this to be part of this PR. I don't want to prematurely optimize things in any case.

[jkalias]

Does this fix follow the policy of Math.NET Spatial?
I'm afraid the case where you maintainers respect RAII(Resource acquisition is initialization) policy and I would break it.

I am not sure I understand your question/concern here, can you explain it again please?

[jkalias]

I think you can entirely skip this variable and use var normal = UnitVector3D.OfVector(matV.Column(matV.Cols - 1));. This better captures semantically the intention, that we need the last eigenvector of V.

[f-frhs]

thanks for comment.
I introduced an explaining variable, theIndex, but it didn't work.

I wanted to say "take the corresponding column of matV as the normal vector, to the smallest singular value's index in the diagonal matrix S"

[jkalias]

Fair enough. Perhaps we can call it smallestEigenvalueColumnIndex. You may have noticed that I have a tendency towards long, descriptive names 😃

[jkalias]

var bestFit = new Plane(normal, c);

[f-frhs]

I'm afraid that I didn't define any properties the best Fit plane should have.
Can I use fittedPlane for this?

[f-frhs]

we should evaluate the best fit based on the orthogonal euclidean distance.

Sorry, I couldn't read all your comments at that moment. my fault.
You had already defined "the best fit" for this problem.
Now, we can use var bestFit = new Plane(normal, c);

[jkalias]

It's not really clear to me, what this test is supposed to guarantee.

[f-frhs]

I totally agree with you.
I think I should remove this test.

[jkalias]

We need more tests, which cover what happens when we have less than 3 points (2, 1, or none).

[f-frhs]

Plane.CreateFittedPlaneFrom(points) is not for the testcases with insufficient points.
Plane.CreateFittedPlaneFrom(points) is for over-determined system.
Instead, Plane.Validate(points) will do for this purpose.

It goes a detailed implementation, Plane.Validate(points) will be something like:

public static void Validate(IEnumerable<Point3D> points)
{
    //      [p0x, p0y, p0z]
    //      [p1x, p1y, p1z]
    // mat= [p2x, p2y, p2z]
    //      [p3x, p3y, p3z]
    //      [p4x, p4y, p4z]
    //      [..., ..., ...]

    //this expression will be fixed, according to the discussion
    var mat = CreateMatrix.DenseOfRowVectors(
                points.Select(p => CreateVector.DenseOfArray(new double[] { p.X, p.Y, p.Z })));

    var rankOfMat = mat.Rank();  // the library uses SVD.
    if(rankOfMat <= 0) // DoF=0
    {
        throw new Exception("The same point?");
    }
    if(rankOfMat <= 1) // DoF=1
    {
        throw new Exception("Are they collinear?");
    }

    // do nothing intentionally
    // because there are no problem, we can span 2d-vector-space with given points
}

If few (2 or less) points are given, the method Plane.TryFromPonits() should return false. It is similar to int.TryParse(string)'s behavior, as mentioned in #184.

[jkalias]

I think this is over-engineering the problem. We don't need a Plane.Validate method, we only need to check in this "best fit method" that before we dive in and calculate the SVD of A, its row count is more than or equal to 3. Something like if (A.Rows < 3) { throw ...} (the syntax might not be 100% perfect)

[jkalias]

AssertGeometry.AreEqual(expectedPlane, actualPlane);

[f-frhs]

True. I wish I could have.
I need time to adjust AssertGeometry.AreEqual(PlaneA, PlaneB) to pass this test
where

PlaneA.Normal == -PlaneB.Normal;
PlaneA.RootPoint == PlaneB.RootPoint;

[jkalias]

This is another problem I spotted in the Plane implementation after I took over the maintenance of this repo; the plane is not defined in its Hesse normal form, which I find problematic.

Then the normal direction is uniquely defined and fixed, for any plane. In addition the equality operator should be adjusted, so that we don't check equality of the root points, but that something like other.Conatins(rootPoint), where we check that the root point is contained in the other plane. There is nothing special in the choice of a root point for a plane, any point lying on it is equally valid.

I wanted to open another PR and fix this (breaking) change for the next (v0.8.0) release. Maybe you want to wait with this PR after I am finished (?)

[jkalias]

This seems like a magic number. Perhaps the values from the original example are not as exact as we want them to be. I would adjust the values in our example, so that the default tolerance would pass the test.

[f-frhs]

ok. thank you for your cooperation.

[jkalias]

  * [ ]  which is nice to evaluate?
    a. (z)axis-aligned distance
    b. orthogonal distance

I think your implementation is correct, we should evaluate the best fit based on the orthogonal euclidean distance.

[jkalias]

Fixed accuracy test in #218