3
\$\begingroup\$

I have a 3D-laser scanner which gives me a series of points of a real-world object. The points are related by a common origin.

I know I can derive the equation of a plane from 3 points and that the "maximum material condition", "flatness" or "thickness" of the blob will be defined by a 4th point which has the greatest height-magnitude relative to the plane of the object.

How do I choose which 3 points from a shape (which happens to look like a CD, as in it's disc-shaped) define the plane of the object?

To account for any ambiguity in the question, I already know the diameter of this disc so even though that is technically the widest part of my blob, I actually want to know how "flat" it is. Basically if this "disc" were resting on a perfectly flat surface, what would be it's highest point?

EDITS: Question has been significantly edited, because I originally asked the "wrong" question. I'm stuck at the part where I'm try to find 3 points that define the plane this disc will rest on, before I find the 4th point that indicates how "thick" it is...

I can guarantee the disc will not rest on it's thin edge as it should almost be level or maybe tilted up to 15° when scanned. Essentially a disc may be scanned at a slight angle/tilt and that should not influence it's maximal thickness.

Note: My points come in as cylindrical coordinates initially, but I already wrote functions to convert back and forth to Cartesian3D and Spherical systems. I know which axis will be X, Y and Z but I have intentionally omitted that information as it should be arbitrary to any code.

Conceptual answers are very welcome, I can google to figure out implementation details.

using Common.FluentValidation;
using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using Common.Extensions;

namespace Common.Mathematics
{
    // Coordinate system conversions http://en.wikipedia.org/wiki/Spherical_coordinate_system#Cartesian_coordinates
    // Cylindrical Coordinates http://en.wikipedia.org/wiki/Cylindrical_coordinate_system
    // http://gamedev.stackexchange.com/questions/81713/how-do-i-translate-a-spherical-coordinate-to-a-cartesian-one
    // http://gamedev.stackexchange.com/questions/44738/spherical-to-cartesian-coordinates

    // World Frame
    /**********************************************************
     *                                                        *
     *                      +Height                           *
     *                   (Polar-Axis)                         *
     *                 (Elevation = 90°)                     *
     *                         |                              *
     *                         |  / +Depth (Rotation = 90°)   *
     *                         | /                            *
     *                         |/                             *
     *  -Length ---------------*--------------- +Length       *
     * (Rotation = 180°)      /|            (Rotation = 0°)   *
     *                       / |            (Elevation = 0°)  *
     *               -Depth /  |                              *
     *     (Rotation = 270°)   |                              *
     *                      -Height                           *
     *                 (Elevation = -90°)                     *
     *                                                        *
     **********************************************************/

    public static class Coordinates
    {
        // Config
        public static double Cartesian_Max_Resolution = 0.0000001;
        public static double Degree_Max_Resolution = 0.0000001;
    }

    public class Cylindrical
    {
        // Public Fields
        /// <summary>
        /// Horizontal vector distance from polar axis to point.
        /// </summary>
        public readonly double Radius;

        /// <summary>
        /// Height along polar axis (Graphically -Down, +Up).
        /// </summary>
        public readonly double Height;

        /// <summary>
        /// Angle (theta) of rotation (+ is counter clockwise) about the polar axis. Rotation of 0° is colinear with the positive Length-axis; 90° is colinear with the positive Depth-axis.
        /// </summary>
        public readonly double Rotation;

        // Constructor
        public Cylindrical(double p_radius, double p_height, double p_rotation)
        {
            // Init
            Radius = p_radius;
            Height = p_height;
            Rotation = p_rotation;

            // Correct Inverted Radius
            if (Radius < 0)
            {
                Radius = -Radius;
                Height = -Height;
                Rotation += 180;
            }

            // Limit Rotation from 0° to 360°
            Rotation = Rotation % 360;
            Rotation = (Rotation < 0) ? Rotation + 360 : Rotation;
        }
        public Cylindrical(Cartesian3D p_coordinate)
        {
            // Init
            double L = p_coordinate.Length;
            double D = p_coordinate.Depth;
            Height = p_coordinate.Height;

            // Horizontal Vector
            Radius = Math.Sqrt(L * L + D * D);

            // Atan2 returns -PI to +PI Radians
            // Adjust angle to be between 0° and 360°
            double DomainValue = Math.Atan2(D, L) * Constants.DEGREES_PER_RADIAN;
            Rotation = ((DomainValue >= 0) ? DomainValue : DomainValue + 360);
        }
        public Cylindrical(Spherical p_coordinate)
        {
            // Init
            double M = p_coordinate.Magnitude;
            double E = p_coordinate.Elevation * Constants.RADIANS_PER_DEGREE;
            Rotation = p_coordinate.Rotation % 360;

            // Calc
            Radius = M * Math.Cos(E);
            Height = M * Math.Sin(E);

            // Limit Rotation from 0° to 360°
            Rotation = Rotation % 360;
            Rotation = (Rotation < 0) ? Rotation + 360 : Rotation;
        }

        // Overrides
        public override bool Equals(object obj)
        {
            if (this == (obj as Cylindrical))
                return true;

            return false;
        }
        public static bool operator ==(Cylindrical A, Cylindrical B)
        {
            // Both null
            if (Object.Equals(A, null) && Object.Equals(B, null))
                return true;

            // Either Null
            if (Object.Equals(A, null) || Object.Equals(B, null))
                return false;

            // Compare fields/properties
            if (Math.Abs(A.Radius - B.Radius) > Coordinates.Cartesian_Max_Resolution)
                return false;

            if (Math.Abs(A.Height - B.Height) > Coordinates.Cartesian_Max_Resolution)
                return false;

            if (Math.Abs(A.Rotation - B.Rotation) > Coordinates.Degree_Max_Resolution)
                return false;

            return true;
        }
        public static bool operator !=(Cylindrical A, Cylindrical B)
        {
            return !(A == B);
        }
        public override int GetHashCode()
        {
            // http://primes.utm.edu/lists/small/1000.txt
            int SomeLowValuePrimeA = 503;
            int SomeLowValuePrimeB = 947;

            // http://stackoverflow.com/a/720282/1718702
            // disable overflow, for the unlikely possibility that you
            // are compiling with overflow-checking enabled
            unchecked
            {
                int Hash = SomeLowValuePrimeA;

                Hash = (SomeLowValuePrimeB * Hash) + Radius.GetHashCode();
                Hash = (SomeLowValuePrimeB * Hash) + Height.GetHashCode();
                Hash = (SomeLowValuePrimeB * Hash) + Rotation.GetHashCode();

                return Hash;
            }
        }
        public override string ToString()
        {
            return
                string.Format("Radius = {0}, Height = {1}, Rotation = {2}°",
                Radius.ToString("0.0000"), Height.ToString("0.0000"), Rotation.ToString("0.00"));
        }

        // Implicit Operators
        public static implicit operator Cylindrical(Spherical p_coordinate)
        {
            return
                new Cylindrical(p_coordinate);
        }
        public static implicit operator Cylindrical(Cartesian3D p_coordinate)
        {
            return
                new Cylindrical(p_coordinate);
        }
    }
    public class Cartesian3D
    {
        // Public Fields
        /// <summary>
        /// Length-axis value (Graphically -Left, +Right)
        /// </summary>
        public readonly double Length;

        /// <summary>
        /// Height-axis value (Graphically -Down, +Up)
        /// </summary>
        public readonly double Height;

        /// <summary>
        /// Depth-axis value (Graphically -Toward camera, +Away from camera)
        /// </summary>
        public readonly double Depth;

        // Constructor
        public Cartesian3D(double p_length, double p_height, double p_depth)
        {
            // Init
            Length = p_length;
            Height = p_height;
            Depth = p_depth;
        }
        public Cartesian3D(Cylindrical p_coordinate)
        {
            // Init
            double Radius = p_coordinate.Radius;
            Height = p_coordinate.Height;
            double Rotation = p_coordinate.Rotation * Constants.RADIANS_PER_DEGREE;

            Length = Radius * Math.Cos(Rotation);
            Depth = Radius * Math.Sin(Rotation);
        }

        // Overrides
        public override bool Equals(object obj)
        {
            if (this == (obj as Cartesian3D))
                return true;

            return false;
        }
        public static bool operator ==(Cartesian3D A, Cartesian3D B)
        {
            // Both null
            if (Object.Equals(A, null) && Object.Equals(B, null))
                return true;

            // Either Null
            if (Object.Equals(A, null) || Object.Equals(B, null))
                return false;

            // Compare fields/properties
            if (Math.Abs(A.Length - B.Length) > Coordinates.Cartesian_Max_Resolution)
                return false;

            if (Math.Abs(A.Height - B.Height) > Coordinates.Cartesian_Max_Resolution)
                return false;

            if (Math.Abs(A.Depth - B.Depth) > Coordinates.Cartesian_Max_Resolution)
                return false;

            return true;
        }
        public static bool operator !=(Cartesian3D A, Cartesian3D B)
        {
            return !(A == B);
        }
        public override int GetHashCode()
        {
            // http://primes.utm.edu/lists/small/1000.txt
            int SomeLowValuePrimeA = 5;
            int SomeLowValuePrimeB = 353;

            // http://stackoverflow.com/a/720282/1718702
            // disable overflow, for the unlikely possibility that you
            // are compiling with overflow-checking enabled
            unchecked
            {
                int Hash = SomeLowValuePrimeA;

                Hash = (SomeLowValuePrimeB * Hash) + Length.GetHashCode();
                Hash = (SomeLowValuePrimeB * Hash) + Height.GetHashCode();
                Hash = (SomeLowValuePrimeB * Hash) + Depth.GetHashCode();

                return Hash;
            }
        }
        public override string ToString()
        {
            return
                string.Format("Length = {0}, Height = {1}, Depth = {2}",
                Length.ToString("0.0000"), Height.ToString("0.0000"), Depth.ToString("0.0000"));
        }

        // Implicit Operators
        public static implicit operator Cartesian3D(Cylindrical p_coordinate)
        {
            return
                new Cartesian3D(p_coordinate);
        }
        public static implicit operator Cartesian3D(Spherical p_coordinate)
        {
            Cylindrical Point = p_coordinate;

            return
                new Cartesian3D(Point);
        }
    }

    public class Spherical
    {
        // Public Fields
        /// <summary>
        /// Magnitude (rho) of vector.
        /// </summary>
        public readonly double Magnitude;

        /// <summary>
        /// Angle (phi) of elevation measured relative to the horizontal plane. Elevation of 90° is colinear with the positive polar-axis; -90° is colinear with the negative polar-axis.
        /// </summary>
        public readonly double Elevation;

        /// <summary>
        /// Angle (theta) of rotation about the polar axis. Rotation of 0° is colinear with positive length-axis; 90° is colinear with positive depth-axis.
        /// </summary>
        public readonly double Rotation;

        // Constructor
        public Spherical(double p_magnitude, double p_elevation, double p_rotation)
        {
            // Init
            Magnitude = p_magnitude;
            Elevation = p_elevation % 180;
            Rotation = p_rotation;

            // Correct Inverted Magnitude
            if (Magnitude < 0)
            {
                Magnitude = -Magnitude;
                Elevation = -Elevation;
                Rotation += 180;
            }

            // Limit Elevation from 90° to -90°
            if (Elevation > 90.0)
            {
                Elevation = 180 - Elevation;
                Rotation += 180;
            }
            else if (Elevation < -90.0)
            {
                Elevation = -180 - Elevation;
                Rotation += 180;
            }

            // Limit Rotation from 0° to 360°
            Rotation = Rotation % 360;
            Rotation = (Rotation < 0) ? Rotation + 360 : Rotation;
        }
        public Spherical(Cylindrical p_coordinate)
        {
            // Init
            double R = p_coordinate.Radius;
            double H = p_coordinate.Height;
            Rotation = p_coordinate.Rotation;

            Magnitude = Math.Sqrt(R * R + H * H);

            // Vector magnitude projected in the horizontal plane will always be positive
            // Therefore angle will always be between 0° and 180°
            Elevation = Math.Atan2(H, R) * Constants.DEGREES_PER_RADIAN;

            // Limit Rotation from 0° to 360°
            Rotation = Rotation % 360;
            Rotation = (Rotation < 0) ? Rotation + 360 : Rotation;
        }

        // Overrides
        public override bool Equals(object obj)
        {
            if (this == (obj as Spherical))
                return true;

            return false;
        }
        public static bool operator ==(Spherical A, Spherical B)
        {
            // Both null
            if (Object.Equals(A, null) && Object.Equals(B, null))
                return true;

            // Either Null
            if (Object.Equals(A, null) || Object.Equals(B, null))
                return false;

            // Compare fields/properties
            if (Math.Abs(A.Magnitude - B.Magnitude) > Coordinates.Cartesian_Max_Resolution)
                return false;

            if (Math.Abs(A.Elevation - B.Elevation) > Coordinates.Degree_Max_Resolution)
                return false;

            if (Math.Abs(A.Rotation - B.Rotation) > Coordinates.Degree_Max_Resolution)
                return false;

            return true;
        }
        public static bool operator !=(Spherical A, Spherical B)
        {
            return !(A == B);
        }
        public override int GetHashCode()
        {
            // http://primes.utm.edu/lists/small/1000.txt
            int SomeLowValuePrimeA = 503;
            int SomeLowValuePrimeB = 151;

            // http://stackoverflow.com/a/720282/1718702
            // disable overflow, for the unlikely possibility that you
            // are compiling with overflow-checking enabled
            unchecked
            {
                int Hash = SomeLowValuePrimeA;

                Hash = (SomeLowValuePrimeB * Hash) + Magnitude.GetHashCode();
                Hash = (SomeLowValuePrimeB * Hash) + Elevation.GetHashCode();
                Hash = (SomeLowValuePrimeB * Hash) + Rotation.GetHashCode();

                return Hash;
            }
        }
        public override string ToString()
        {
            return
                string.Format("Magnitude = {0}, Elevation = {1}°, Rotation = {2}°",
                Magnitude.ToString("0.0000"), Elevation.ToString("0.00"), Rotation.ToString("0.00"));
        }

        // Implicit Operators
        public static implicit operator Spherical(Cylindrical p_coordinate)
        {
            return
                new Spherical(p_coordinate);
        }
        public static implicit operator Spherical(Cartesian3D p_coordinate)
        {
            Cylindrical Point = p_coordinate;

            return
                new Spherical(Point);
        }
    }
}
\$\endgroup\$
  • \$\begingroup\$ You seem to be assuming there is only 1 way it can come to rest on the ground which seems like a bad assumption. And without that assumption being true I don't think your question is answerable. \$\endgroup\$ – ClassicThunder Sep 18 '14 at 20:20
  • \$\begingroup\$ hmm... valid point. I will rephrase the question then. Ultimately I need to determine maximum thickness of as disc (think CD shaped, but potentially "warped"). One moment and I will edit. I was trying to ask a bite-sized question instead of asking for everything at once. \$\endgroup\$ – HodlDwon Sep 18 '14 at 20:22
  • \$\begingroup\$ Ok, edits done. Hopefully that's more to the point. \$\endgroup\$ – HodlDwon Sep 18 '14 at 20:55
4
\$\begingroup\$

I'm not entirely sure I understand the question but, given a 3D point cloud with a 2D "shape" (like a disc) the problem is simply one of dimension reduction.

The algorithm you're looking for is called "Principal Component Analysis". How it works is not trivial, so I'll leave that for you to research if you're not already familiar with PCA. What it will give you is a series of eigenvectors and eigenvalues that will tell you how significant a dimension is. Here's a quick sketch of how this would look in 2D with inherently 1D data (but easily expanded to ND):

PCA Example

This information is usually used to simplify multi-dimensional data by reducing dimensions and removing ones that aren't that significant (e.g. in the above example the red; the data can be treated as linear along the blue without the loss of a lot of information). In your case however, you can find the smallest eigenvalue, and use that as your normal for measuring disk thickness in 3D (or alternatively, cross the two eigenvectors with the largest corresponding eigenvalues and the resulting vector is the normal to the disc pointing up relative to it).

Once you have that vector, it's trivial to find out which point is furthest along on the axis that vector creates. And there you have your thickness!

\$\endgroup\$
  • \$\begingroup\$ Hey thanks, I was not aware of the concept of PCA at all. It looks like exactly what I need though. This is actually very good as it looks like a very stable way to find the normal vector to the plane of the disc. That's important because if I scan the same disc twice, I want to get the same answer! Now to figure out the math... \$\endgroup\$ – HodlDwon Sep 19 '14 at 15:59
  • 1
    \$\begingroup\$ +1 for @Paraknight . Besides PCA, if the OP is looking for a robust algorithm for fitting planar faces to a point cloud, RANSAC may come in handy -> en.wikipedia.org/wiki/RANSAC \$\endgroup\$ – teodron Sep 19 '14 at 16:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.