lightweight_geom_data.cpp

Go to the documentation of this file.

/* +------------------------------------------------------------------------+
   |                     Mobile Robot Programming Toolkit (MRPT)            |
   |                          http://www.mrpt.org/                          |
   |                                                                        |
   | Copyright (c) 2005-2018, Individual contributors, see AUTHORS file     |
   | See: http://www.mrpt.org/Authors - All rights reserved.                |
   | Released under BSD License. See details in http://www.mrpt.org/License |
   +------------------------------------------------------------------------+ */
 
#include "math-precomp.h"  // Precompiled headers
 
#include <mrpt/core/exceptions.h>
#include <mrpt/math/lightweight_geom_data.h>
#include <mrpt/math/CQuaternion.h>
#include <mrpt/math/CMatrixFixedNumeric.h>
#include <mrpt/math/geometry.h>  // distance()
#include <mrpt/math/ops_containers.h>
#include <mrpt/math/homog_matrices.h>
#include <mrpt/serialization/stl_serialization.h>
#include <mrpt/serialization/CArchive.h>  // impl of << operator
 
using namespace std;  // For min/max, etc...
using namespace mrpt::serialization;  // CArchive, << operator for STL
 
using mrpt::DEG2RAD;
using mrpt::RAD2DEG;
 
namespace mrpt
{
namespace math
{
TPoint2D::TPoint2D(const TPose2D& p) : x(p.x), y(p.y) {}
TPoint2D::TPoint2D(const TPoint3D& p) : x(p.x), y(p.y) {}
TPoint2D::TPoint2D(const TPose3D& p) : x(p.x), y(p.y) {}
bool TPoint2D::operator<(const TPoint2D& p) const
{
        if (x < p.x)
                return true;
        else if (x > p.x)
                return false;
        else
                return y < p.y;
}
void TPoint2D::fromString(const std::string& s)
{
        CMatrixDouble m;
        if (!m.fromMatlabStringFormat(s))
                THROW_EXCEPTION("Malformed expression in ::fromString");
        ASSERTMSG_(
                m.rows() == 1 && m.cols() == 2, "Wrong size of vector in ::fromString");
        x = m.get_unsafe(0, 0);
        y = m.get_unsafe(0, 1);
}
 
TPose2D::TPose2D(const TPoint2D& p) : x(p.x), y(p.y), phi(0.0) {}
TPose2D::TPose2D(const TPoint3D& p) : x(p.x), y(p.y), phi(0.0) {}
TPose2D::TPose2D(const TPose3D& p) : x(p.x), y(p.y), phi(p.yaw) {}
void TPose2D::asString(std::string& s) const
{
        s = mrpt::format("[%f %f %f]", x, y, RAD2DEG(phi));
}
void TPose2D::fromString(const std::string& s)
{
        CMatrixDouble m;
        if (!m.fromMatlabStringFormat(s))
                THROW_EXCEPTION("Malformed expression in ::fromString");
        ASSERTMSG_(
                m.rows() == 1 && m.cols() == 3, "Wrong size of vector in ::fromString");
        x = m.get_unsafe(0, 0);
        y = m.get_unsafe(0, 1);
        phi = DEG2RAD(m.get_unsafe(0, 2));
}
mrpt::math::TPose2D mrpt::math::TPose2D::operator+(
        const mrpt::math::TPose2D& b) const
{
        const double A_cosphi = cos(this->phi), A_sinphi = sin(this->phi);
        // Use temporary variables for the cases (A==this) or (B==this)
        const double new_x = this->x + b.x * A_cosphi - b.y * A_sinphi;
        const double new_y = this->y + b.x * A_sinphi + b.y * A_cosphi;
        const double new_phi = mrpt::math::wrapToPi(this->phi + b.phi);
 
        return mrpt::math::TPose2D(new_x, new_y, new_phi);
}
 
mrpt::math::TPose2D mrpt::math::TPose2D::operator-(
        const mrpt::math::TPose2D& b) const
{
        const double B_cosphi = cos(b.phi), B_sinphi = sin(b.phi);
 
        const double new_x =
                (this->x - b.x) * B_cosphi + (this->y - b.y) * B_sinphi;
        const double new_y =
                -(this->x - b.x) * B_sinphi + (this->y - b.y) * B_cosphi;
        const double new_phi = mrpt::math::wrapToPi(this->phi - b.phi);
 
        return mrpt::math::TPose2D(new_x, new_y, new_phi);
}
 
// ----
void TTwist2D::asString(std::string& s) const
{
        s = mrpt::format("[%f %f %f]", vx, vy, RAD2DEG(omega));
}
void TTwist2D::fromString(const std::string& s)
{
        CMatrixDouble m;
        if (!m.fromMatlabStringFormat(s))
                THROW_EXCEPTION("Malformed expression in ::fromString");
        ASSERTMSG_(
                m.rows() == 1 && m.cols() == 3, "Wrong size of vector in ::fromString");
        vx = m.get_unsafe(0, 0);
        vy = m.get_unsafe(0, 1);
        omega = DEG2RAD(m.get_unsafe(0, 2));
}
// Transform the (vx,vy) components for a counterclockwise rotation of `ang`
// radians
void TTwist2D::rotate(const double ang)
{
        const double nvx = vx * cos(ang) - vy * sin(ang);
        const double nvy = vx * sin(ang) + vy * cos(ang);
        vx = nvx;
        vy = nvy;
}
bool TTwist2D::operator==(const TTwist2D& o) const
{
        return vx == o.vx && vy == o.vy && omega == o.omega;
}
bool TTwist2D::operator!=(const TTwist2D& o) const { return !(*this == o); }
mrpt::math::TPose2D TTwist2D::operator*(const double dt) const
{
        return mrpt::math::TPose2D(vx * dt, vy * dt, omega * dt);
}
 
// ----
void TTwist3D::asString(std::string& s) const
{
        s = mrpt::format(
                "[%f %f %f  %f %f %f]", vx, vy, vz, RAD2DEG(wx), RAD2DEG(wy),
                RAD2DEG(wz));
}
void TTwist3D::fromString(const std::string& s)
{
        CMatrixDouble m;
        if (!m.fromMatlabStringFormat(s))
                THROW_EXCEPTION("Malformed expression in ::fromString");
        ASSERTMSG_(
                m.rows() == 1 && m.cols() == 6, "Wrong size of vector in ::fromString");
        for (int i = 0; i < 3; i++) (*this)[i] = m.get_unsafe(0, i);
        for (int i = 0; i < 3; i++)
                (*this)[3 + i] = DEG2RAD(m.get_unsafe(0, 3 + i));
}
// Transform all 6 components for a change of reference frame from "A" to
// another frame "B" whose rotation with respect to "A" is given by `rot`. The
// translational part of the pose is ignored
void TTwist3D::rotate(const TPose3D& rot)
{
        const TTwist3D t = *this;
        mrpt::math::CMatrixDouble33 R;
        rot.getRotationMatrix(R);
        vx = R(0, 0) * t.vx + R(0, 1) * t.vy + R(0, 2) * t.vz;
        vy = R(1, 0) * t.vx + R(1, 1) * t.vy + R(1, 2) * t.vz;
        vz = R(2, 0) * t.vx + R(2, 1) * t.vy + R(2, 2) * t.vz;
 
        wx = R(0, 0) * t.wx + R(0, 1) * t.wy + R(0, 2) * t.wz;
        wy = R(1, 0) * t.wx + R(1, 1) * t.wy + R(1, 2) * t.wz;
        wz = R(2, 0) * t.wx + R(2, 1) * t.wy + R(2, 2) * t.wz;
}
bool TTwist3D::operator==(const TTwist3D& o) const
{
        return vx == o.vx && vy == o.vy && vz == o.vz && wx == o.wx && wy == o.wy &&
                   wz == o.wz;
}
bool TTwist3D::operator!=(const TTwist3D& o) const { return !(*this == o); }
TPoint3D::TPoint3D(const TPoint2D& p) : x(p.x), y(p.y), z(0.0) {}
TPoint3D::TPoint3D(const TPose2D& p) : x(p.x), y(p.y), z(0.0) {}
TPoint3D::TPoint3D(const TPose3D& p) : x(p.x), y(p.y), z(p.z) {}
bool TPoint3D::operator<(const TPoint3D& p) const
{
        if (x < p.x)
                return true;
        else if (x > p.x)
                return false;
        else if (y < p.y)
                return true;
        else if (y > p.y)
                return false;
        else
                return z < p.z;
}
void TPoint3D::fromString(const std::string& s)
{
        CMatrixDouble m;
        if (!m.fromMatlabStringFormat(s))
                THROW_EXCEPTION("Malformed expression in ::fromString");
        ASSERTMSG_(
                m.rows() == 1 && m.cols() == 3, "Wrong size of vector in ::fromString");
        x = m.get_unsafe(0, 0);
        y = m.get_unsafe(0, 1);
        z = m.get_unsafe(0, 2);
}
 
TPose3D::TPose3D(const TPoint2D& p)
        : x(p.x), y(p.y), z(0.0), yaw(0.0), pitch(0.0), roll(0.0)
{
}
TPose3D::TPose3D(const TPose2D& p)
        : x(p.x), y(p.y), z(0.0), yaw(p.phi), pitch(0.0), roll(0.0)
{
}
TPose3D::TPose3D(const TPoint3D& p)
        : x(p.x), y(p.y), z(p.z), yaw(0.0), pitch(0.0), roll(0.0)
{
}
void TPose3D::asString(std::string& s) const
{
        s = mrpt::format(
                "[%f %f %f %f %f %f]", x, y, z, RAD2DEG(yaw), RAD2DEG(pitch),
                RAD2DEG(roll));
}
void TPose3D::getAsQuaternion(
        mrpt::math::CQuaternion<double>& q,
        mrpt::math::CMatrixFixedNumeric<double, 4, 3>* out_dq_dr) const
{
        // See:
        // http://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
        const double cy = cos(yaw * 0.5), sy = sin(yaw * 0.5);
        const double cp = cos(pitch * 0.5), sp = sin(pitch * 0.5);
        const double cr = cos(roll * 0.5), sr = sin(roll * 0.5);
 
        const double ccc = cr * cp * cy;
        const double ccs = cr * cp * sy;
        const double css = cr * sp * sy;
        const double sss = sr * sp * sy;
        const double scc = sr * cp * cy;
        const double ssc = sr * sp * cy;
        const double csc = cr * sp * cy;
        const double scs = sr * cp * sy;
 
        q[0] = ccc + sss;
        q[1] = scc - css;
        q[2] = csc + scs;
        q[3] = ccs - ssc;
 
        // Compute 4x3 Jacobian: for details, see technical report:
        //   Parameterizations of SE(3) transformations: equivalences, compositions
        //   and uncertainty, J.L. Blanco (2010).
        //   http://www.mrpt.org/6D_poses:equivalences_compositions_and_uncertainty
        if (out_dq_dr)
        {
                alignas(MRPT_MAX_ALIGN_BYTES) const double nums[4 * 3] = {
                        -0.5 * q[3], 0.5 * (-csc + scs), -0.5 * q[1],
                        -0.5 * q[2], 0.5 * (-ssc - ccs), 0.5 * q[0],
                        0.5 * q[1],  0.5 * (ccc - sss),  0.5 * q[3],
                        0.5 * q[0],  0.5 * (-css - scc), -0.5 * q[2]};
                out_dq_dr->loadFromArray(nums);
        }
}
void TPose3D::composePoint(const TPoint3D l, TPoint3D& g) const
{
        CMatrixDouble33 R;
        this->getRotationMatrix(R);
        TPoint3D res;
        res.x = R(0, 0) * l.x + R(0, 1) * l.y + R(0, 2) * l.z + this->x;
        res.y = R(1, 0) * l.x + R(1, 1) * l.y + R(1, 2) * l.z + this->y;
        res.z = R(2, 0) * l.x + R(2, 1) * l.y + R(2, 2) * l.z + this->z;
 
        g = res;
}
void TPose3D::inverseComposePoint(const TPoint3D g, TPoint3D& l) const
{
        CMatrixDouble44 H;
        this->getInverseHomogeneousMatrix(H);
        TPoint3D res;
        res.x = H(0, 0) * g.x + H(0, 1) * g.y + H(0, 2) * g.z + H(0, 3);
        res.y = H(1, 0) * g.x + H(1, 1) * g.y + H(1, 2) * g.z + H(1, 3);
        res.z = H(2, 0) * g.x + H(2, 1) * g.y + H(2, 2) * g.z + H(2, 3);
 
        l = res;
}
void TPose3D::getRotationMatrix(mrpt::math::CMatrixDouble33& R) const
{
        const double cy = cos(yaw);
        const double sy = sin(yaw);
        const double cp = cos(pitch);
        const double sp = sin(pitch);
        const double cr = cos(roll);
        const double sr = sin(roll);
 
        alignas(MRPT_MAX_ALIGN_BYTES) const double rot_vals[] = {cy * cp,
                                                                                   cy * sp * sr - sy * cr,
                                                                                   cy * sp * cr + sy * sr,
                                                                                   sy * cp,
                                                                                   sy * sp * sr + cy * cr,
                                                                                   sy * sp * cr - cy * sr,
                                                                                   -sp,
                                                                                   cp * sr,
                                                                                   cp * cr};
        R.loadFromArray(rot_vals);
}
void TPose3D::SO3_to_yaw_pitch_roll(
        const mrpt::math::CMatrixDouble33& R, double& yaw, double& pitch,
        double& roll)
{
        ASSERTDEBMSG_(
                std::abs(
                        sqrt(square(R(0, 0)) + square(R(1, 0)) + square(R(2, 0))) - 1) <
                        3e-3,
                "Homogeneous matrix is not orthogonal & normalized!: " +
                        R.inMatlabFormat());
        ASSERTDEBMSG_(
                std::abs(
                        sqrt(square(R(0, 1)) + square(R(1, 1)) + square(R(2, 1))) - 1) <
                        3e-3,
                "Homogeneous matrix is not orthogonal & normalized!: " +
                        R.inMatlabFormat());
        ASSERTDEBMSG_(
                std::abs(
                        sqrt(square(R(0, 2)) + square(R(1, 2)) + square(R(2, 2))) - 1) <
                        3e-3,
                "Homogeneous matrix is not orthogonal & normalized!: " +
                        R.inMatlabFormat());
 
        // Pitch is in the range [-pi/2, pi/2 ], so this calculation is enough:
        pitch = atan2(-R(2, 0), hypot(R(0, 0), R(1, 0)));
 
        // Roll:
        if ((fabs(R(2, 1)) + fabs(R(2, 2))) <
                10 * std::numeric_limits<double>::epsilon())
        {
                // Gimbal lock between yaw and roll. This one is arbitrarily forced to
                // be zero.
                // Check
                // http://reference.mrpt.org/devel/classmrpt_1_1poses_1_1_c_pose3_d.html.
                // If cos(pitch)==0, the homogeneous matrix is:
                // When sin(pitch)==1:
                //  /0  cysr-sycr cycr+sysr x\   /0  sin(r-y) cos(r-y)  x\.
                //  |0  sysr+cycr sycr-cysr y| = |0  cos(r-y) -sin(r-y) y|
                //  |-1     0         0     z|   |-1    0         0     z|
                //  \0      0         0     1/   \0     0         0     1/
                //
                // And when sin(pitch)=-1:
                //  /0 -cysr-sycr -cycr+sysr x\   /0 -sin(r+y) -cos(r+y) x\.
                //  |0 -sysr+cycr -sycr-cysr y| = |0 cos(r+y)  -sin(r+y) y|
                //  |1      0          0     z|   |1    0          0     z|
                //  \0      0          0     1/   \0    0          0     1/
                //
                // Both cases are in a "gimbal lock" status. This happens because pitch
                // is vertical.
 
                roll = 0.0;
                if (pitch > 0)
                        yaw = atan2(R(1, 2), R(0, 2));
                else
                        yaw = atan2(-R(1, 2), -R(0, 2));
        }
        else
        {
                roll = atan2(R(2, 1), R(2, 2));
                // Yaw:
                yaw = atan2(R(1, 0), R(0, 0));
        }
}
 
void TPose3D::fromHomogeneousMatrix(const mrpt::math::CMatrixDouble44& HG)
{
        SO3_to_yaw_pitch_roll(HG.block<3, 3>(0, 0), yaw, pitch, roll);
        x = HG(0, 3);
        y = HG(1, 3);
        z = HG(2, 3);
}
void TPose3D::composePose(const TPose3D other, TPose3D& result) const
{
        CMatrixDouble44 me_H, o_H;
        this->getHomogeneousMatrix(me_H);
        other.getHomogeneousMatrix(o_H);
        result.fromHomogeneousMatrix(me_H * o_H);
}
void TPose3D::getHomogeneousMatrix(mrpt::math::CMatrixDouble44& HG) const
{
        CMatrixDouble33 R;
        getRotationMatrix(R);
        HG.block<3, 3>(0, 0) = R;
        HG(0, 3) = x;
        HG(1, 3) = y;
        HG(2, 3) = z;
        HG(3, 0) = HG(3, 1) = HG(3, 2) = 0.;
        HG(3, 3) = 1.;
}
void TPose3D::getInverseHomogeneousMatrix(mrpt::math::CMatrixDouble44& HG) const
{  // Get current HM & inverse in-place:
        this->getHomogeneousMatrix(HG);
        mrpt::math::homogeneousMatrixInverse(HG);
}
 
void TPose3D::fromString(const std::string& s)
{
        CMatrixDouble m;
        if (!m.fromMatlabStringFormat(s))
                THROW_EXCEPTION("Malformed expression in ::fromString");
        ASSERTMSG_(
                m.rows() == 1 && m.cols() == 6, "Wrong size of vector in ::fromString");
        x = m.get_unsafe(0, 0);
        y = m.get_unsafe(0, 1);
        z = m.get_unsafe(0, 2);
        yaw = DEG2RAD(m.get_unsafe(0, 3));
        pitch = DEG2RAD(m.get_unsafe(0, 4));
        roll = DEG2RAD(m.get_unsafe(0, 5));
}
 
void TPose3DQuat::fromString(const std::string& s)
{
        CMatrixDouble m;
        if (!m.fromMatlabStringFormat(s))
                THROW_EXCEPTION("Malformed expression in ::fromString");
        ASSERTMSG_(
                m.rows() == 1 && m.cols() == 7, "Wrong size of vector in ::fromString");
        for (int i = 0; i < m.cols(); i++) (*this)[i] = m.get_unsafe(0, i);
}
TPose3D operator-(const TPose3D& p)
{
        CMatrixDouble44 H;
        p.getInverseHomogeneousMatrix(H);
        TPose3D ret;
        ret.fromHomogeneousMatrix(H);
        return ret;
}
TPose3D operator-(const TPose3D& b, const TPose3D& a)
{
        // b - a = A^{-1} * B
        CMatrixDouble44 Hainv, Hb;
        a.getInverseHomogeneousMatrix(Hainv);
        b.getHomogeneousMatrix(Hb);
        TPose3D ret;
        ret.fromHomogeneousMatrix(Hainv * Hb);
        return ret;
}
 
// Text streaming:
std::ostream& operator<<(std::ostream& o, const TPoint2D& p)
{
        return (o << p.asString());
}
std::ostream& operator<<(std::ostream& o, const TPoint3D& p)
{
        return (o << p.asString());
}
std::ostream& operator<<(std::ostream& o, const TPose2D& p)
{
        return (o << p.asString());
}
std::ostream& operator<<(std::ostream& o, const TPose3D& p)
{
        return (o << p.asString());
}
std::ostream& operator<<(std::ostream& o, const TPose3DQuat& p)
{
        return (o << p.asString());
}
 
CArchive& operator>>(CArchive& in, mrpt::math::TSegment2D& s)
{
        return in >> s.point1 >> s.point2;
}
CArchive& operator<<(CArchive& out, const mrpt::math::TSegment2D& s)
{
        return out << s.point1 << s.point2;
}
CArchive& operator>>(CArchive& in, mrpt::math::TLine2D& l)
{
        return in >> l.coefs[0] >> l.coefs[1] >> l.coefs[2];
}
CArchive& operator<<(CArchive& out, const mrpt::math::TLine2D& l)
{
        return out << l.coefs[0] << l.coefs[1] << l.coefs[2];
}
 
CArchive& operator>>(CArchive& in, mrpt::math::TSegment3D& s)
{
        return in >> s.point1 >> s.point2;
}
CArchive& operator<<(CArchive& out, const mrpt::math::TSegment3D& s)
{
        return out << s.point1 << s.point2;
}
CArchive& operator>>(CArchive& in, mrpt::math::TLine3D& l)
{
        return in >> l.pBase >> l.director[0] >> l.director[1] >> l.director[2];
}
CArchive& operator<<(CArchive& out, const mrpt::math::TLine3D& l)
{
        return out << l.pBase << l.director[0] << l.director[1] << l.director[2];
}
CArchive& operator>>(CArchive& in, mrpt::math::TPlane& p)
{
        return in >> p.coefs[0] >> p.coefs[1] >> p.coefs[2] >> p.coefs[3];
}
CArchive& operator<<(CArchive& out, const mrpt::math::TPlane& p)
{
        return out << p.coefs[0] << p.coefs[1] << p.coefs[2] << p.coefs[3];
}
 
double TSegment2D::length() const { return math::distance(point1, point2); }
double TSegment2D::distance(const TPoint2D& point) const
{
        return std::abs(signedDistance(point));
}
double TSegment2D::signedDistance(const TPoint2D& point) const
{
        // It is reckoned whether the perpendicular line to the TSegment2D which
        // passes through point crosses or not the referred segment,
        // or what is the same, whether point makes an obtuse triangle with the
        // segment or not (being the longest segment one between the point and
        // either end of TSegment2D).
        const double d1 = math::distance(point, point1);
        if (point1 == point2) return d1;
 
        const double d2 = math::distance(point, point2);
        const double d3 = length();
        const double ds1 = square(d1);
        const double ds2 = square(d2);
        const double ds3 = square(d3);
        if (ds1 > (ds2 + ds3) || ds2 > (ds1 + ds3))
                // Fix sign:
                return min(d1, d2) *
                           (TLine2D(*this).signedDistance(point) < 0 ? -1 : 1);
        else
                return TLine2D(*this).signedDistance(point);
}
bool TSegment2D::contains(const TPoint2D& point) const
{
        return abs(math::distance(point1, point) + math::distance(point2, point) -
                           math::distance(point1, point2)) < getEpsilon();
}
void TSegment2D::generate3DObject(TSegment3D& s) const
{
        s = TSegment3D(*this);
}
TSegment2D::TSegment2D(const TSegment3D& s)
{
        point1 = TPoint2D(s.point1);
        point2 = TPoint2D(s.point2);
        if (point1 == point2)
                throw std::logic_error("Segment is normal to projection plane");
}
 
bool TSegment2D::operator<(const TSegment2D& s) const
{
        if (point1 < s.point1)
                return true;
        else if (s.point1 < point1)
                return false;
        else
                return point2 < s.point2;
}
 
double TSegment3D::length() const { return math::distance(point1, point2); }
double TSegment3D::distance(const TPoint3D& point) const
{
        return min(
                min(math::distance(point, point1), math::distance(point, point2)),
                TLine3D(*this).distance(point));
}
double TSegment3D::distance(const TSegment3D& segment) const
{
        Eigen::Vector3d u, v, w;
        TPoint3D diff_vect = point2 - point1;
        diff_vect.getAsVector(u);
        diff_vect = segment.point2 - segment.point1;
        diff_vect.getAsVector(v);
        diff_vect = point1 - segment.point1;
        diff_vect.getAsVector(w);
        double a = u.dot(u);  // always >= 0
        double b = u.dot(v);
        double c = v.dot(v);  // always >= 0
        double d = u.dot(w);
        double e = v.dot(w);
        double D = a * c - b * b;  // always >= 0
        double sc, sN, sD = D;  // sc = sN / sD, default sD = D >= 0
        double tc, tN, tD = D;  // tc = tN / tD, default tD = D >= 0
 
        // compute the line parameters of the two closest points
        if (D < 0.00000001)
        {  // the lines are almost parallel
                sN = 0.0;  // force using point P0 on segment S1
                sD = 1.0;  // to prevent possible division by 0.0 later
                tN = e;
                tD = c;
        }
        else
        {  // get the closest points on the infinite lines
                sN = (b * e - c * d);
                tN = (a * e - b * d);
                if (sN < 0.0)
                {  // sc < 0 => the s=0 edge is visible
                        sN = 0.0;
                        tN = e;
                        tD = c;
                }
                else if (sN > sD)
                {  // sc > 1 => the s=1 edge is visible
                        sN = sD;
                        tN = e + b;
                        tD = c;
                }
        }
 
        if (tN < 0.0)
        {  // tc < 0 => the t=0 edge is visible
                tN = 0.0;
                // recompute sc for this edge
                if (-d < 0.0)
                        sN = 0.0;
                else if (-d > a)
                        sN = sD;
                else
                {
                        sN = -d;
                        sD = a;
                }
        }
        else if (tN > tD)
        {  // tc > 1 => the t=1 edge is visible
                tN = tD;
                // recompute sc for this edge
                if ((-d + b) < 0.0)
                        sN = 0;
                else if ((-d + b) > a)
                        sN = sD;
                else
                {
                        sN = (-d + b);
                        sD = a;
                }
        }
        // finally do the division to get sc and tc
        sc = (fabs(sN) < 0.00000001 ? 0.0 : sN / sD);
        tc = (fabs(tN) < 0.00000001 ? 0.0 : tN / tD);
 
        // get the difference of the two closest points
        CVectorDouble dP = w + (sc * u) - (tc * v);  // = S1(sc) - S2(tc)
 
        return dP.norm();  // return the closest distance
}
bool TSegment3D::contains(const TPoint3D& point) const
{
        // Not very intuitive, but very fast, method.
        return abs(math::distance(point1, point) + math::distance(point2, point) -
                           math::distance(point1, point2)) < getEpsilon();
}
 
bool TSegment3D::operator<(const TSegment3D& s) const
{
        if (point1 < s.point1)
                return true;
        else if (s.point1 < point1)
                return false;
        else
                return point2 < s.point2;
}
 
double TLine2D::evaluatePoint(const TPoint2D& point) const
{
        return coefs[0] * point.x + coefs[1] * point.y + coefs[2];
}
bool TLine2D::contains(const TPoint2D& point) const
{
        return abs(distance(point)) < getEpsilon();
}
double TLine2D::distance(const TPoint2D& point) const
{
        return abs(evaluatePoint(point)) /
                   sqrt(coefs[0] * coefs[0] + coefs[1] * coefs[1]);
}
double TLine2D::signedDistance(const TPoint2D& point) const
{
        return evaluatePoint(point) /
                   sqrt(coefs[0] * coefs[0] + coefs[1] * coefs[1]);
}
void TLine2D::getNormalVector(double (&vector)[2]) const
{
        vector[0] = coefs[0];
        vector[1] = coefs[1];
}
void TLine2D::unitarize()
{
        double s = sqrt(coefs[0] * coefs[0] + coefs[1] * coefs[1]);
        for (size_t i = 0; i < 3; i++) coefs[i] /= s;
}
void TLine2D::getDirectorVector(double (&vector)[2]) const
{
        vector[0] = -coefs[1];
        vector[1] = coefs[0];
}
void TLine2D::generate3DObject(TLine3D& l) const { l = TLine3D(*this); }
void TLine2D::getAsPose2D(TPose2D& outPose) const
{
        // Line's director vector is (-coefs[1],coefs[0]).
        // If line is horizontal, force x=0. Else, force y=0. In both cases, we'll
        // find a suitable point.
        outPose.phi = atan2(coefs[0], -coefs[1]);
        if (abs(coefs[0]) < getEpsilon())
        {
                outPose.x = 0;
                outPose.y = -coefs[2] / coefs[1];
        }
        else
        {
                outPose.x = -coefs[2] / coefs[0];
                outPose.y = 0;
        }
}
void TLine2D::getAsPose2DForcingOrigin(
        const TPoint2D& origin, TPose2D& outPose) const
{
        if (!contains(origin))
                throw std::logic_error("Base point is not contained in the line");
        outPose = origin;
        // Line's director vector is (-coefs[1],coefs[0]).
        outPose.phi = atan2(coefs[0], -coefs[1]);
}
TLine2D::TLine2D(const TPoint2D& p1, const TPoint2D& p2)
{
        if (p1 == p2) throw logic_error("Both points are the same");
        coefs[0] = p2.y - p1.y;
        coefs[1] = p1.x - p2.x;
        coefs[2] = p2.x * p1.y - p2.y * p1.x;
}
TLine2D::TLine2D(const TSegment2D& s)
{
        coefs[0] = s.point2.y - s.point1.y;
        coefs[1] = s.point1.x - s.point2.x;
        coefs[2] = s.point2.x * s.point1.y - s.point2.y * s.point1.x;
        // unitarize(); //¿?
}
TLine2D::TLine2D(const TLine3D& l)
{
        // Line's projection to Z plane may be a point.
        if (hypot(l.director[0], l.director[1]) < getEpsilon())
                throw std::logic_error("Line is normal to projection plane");
        coefs[0] = -l.director[1];
        coefs[1] = l.director[0];
        coefs[2] = l.pBase.x * l.director[1] - l.pBase.y * l.director[0];
}
 
bool TLine3D::contains(const TPoint3D& point) const
{
        double dx = point.x - pBase.x;
        double dy = point.y - pBase.y;
        double dz = point.z - pBase.z;
        if (abs(dx) < getEpsilon() && abs(dy) < getEpsilon() &&
                abs(dz) < getEpsilon())
                return true;
        //       dx          dy          dz
        // if -----------=-----------=-----------, point is inside the line.
        //   director[0] director[1] director[2]
        return (abs(dx * director[1] - dy * director[0]) < getEpsilon()) &&
                   (abs(dx * director[2] - dz * director[0]) < getEpsilon()) &&
                   (abs(dy * director[2] - dz * director[1]) < getEpsilon());
}
double TLine3D::distance(const TPoint3D& point) const
{
        // Let d be line's base point minus the argument. Then,
        // d·director/(|d|·|director|) equals both vector's cosine.
        // So, d·director/|director| equals d's projection over director. Then,
        // distance is sqrt(|d|²-(d·director/|director|)²).
        double d[3] = {point.x - pBase.x, point.y - pBase.y, point.z - pBase.z};
        double dv = 0, d2 = 0, v2 = 0;
        for (size_t i = 0; i < 3; i++)
        {
                dv += d[i] * director[i];
                d2 += d[i] * d[i];
                v2 += director[i] * director[i];
        }
        return sqrt(d2 - (dv * dv) / v2);
}
void TLine3D::unitarize()
{
        double s = sqrt(squareNorm<3, double>(director));
        for (size_t i = 0; i < 3; i++) director[i] /= s;
}
TLine3D::TLine3D(const TPoint3D& p1, const TPoint3D& p2)
{
        if (abs(math::distance(p1, p2)) < getEpsilon())
                throw logic_error("Both points are the same");
        pBase = p1;
        director[0] = p2.x - p1.x;
        director[1] = p2.y - p1.y;
        director[2] = p2.z - p1.z;
}
TLine3D::TLine3D(const TSegment3D& s)
{
        pBase = s.point1;
        director[0] = s.point2.x - s.point1.x;
        director[1] = s.point2.y - s.point1.y;
        director[2] = s.point2.z - s.point1.z;
}
TLine3D::TLine3D(const TLine2D& l)
{
        director[0] = -l.coefs[1];
        director[1] = l.coefs[0];
        director[2] = 0;
        // We assume that either l.coefs[0] or l.coefs[1] is not null. Respectively,
        // either y or x can be used as free cordinate.
        if (abs(l.coefs[0]) >= getEpsilon())
        {
                pBase.x = -l.coefs[2] / l.coefs[0];
                pBase.y = 0;
        }
        else
        {
                pBase.x = 0;
                pBase.y = -l.coefs[1] / l.coefs[0];
        }
        pBase.z = 0;
}
 
double TPlane::evaluatePoint(const TPoint3D& point) const
{
        return dotProduct<3, double>(coefs, point) + coefs[3];
}
bool TPlane::contains(const TPoint3D& point) const
{
        return distance(point) < getEpsilon();
}
bool TPlane::contains(const TLine3D& line) const
{
        if (!contains(line.pBase)) return false;  // Base point must be contained
        return abs(getAngle(*this, line)) <
                   getEpsilon();  // Plane's normal must be normal to director vector
}
double TPlane::distance(const TPoint3D& point) const
{
        return abs(evaluatePoint(point)) / sqrt(squareNorm<3, double>(coefs));
}
double TPlane::distance(const TLine3D& line) const
{
        if (abs(getAngle(*this, line)) >= getEpsilon())
                return 0;  // Plane crosses with line
        else
                return distance(line.pBase);
}
void TPlane::getNormalVector(double (&vector)[3]) const
{
        vector[0] = coefs[0];
        vector[1] = coefs[1];
        vector[2] = coefs[2];
}
void TPlane::unitarize()
{
        double s = sqrt(squareNorm<3, double>(coefs));
        for (size_t i = 0; i < 4; i++) coefs[i] /= s;
}
 
// Returns a 6D pose such as its XY plane coincides with the plane
void TPlane::getAsPose3D(mrpt::math::TPose3D& outPose)
{
        double normal[3];
        getUnitaryNormalVector(normal);
        CMatrixDouble44 AXIS;
        generateAxisBaseFromDirectionAndAxis(normal, 2, AXIS);
        for (size_t i = 0; i < 3; i++)
                if (abs(coefs[i]) >= getEpsilon())
                {
                        AXIS.set_unsafe(i, 3, -coefs[3] / coefs[i]);
                        break;
                }
        outPose.fromHomogeneousMatrix(AXIS);
}
void TPlane::getAsPose3DForcingOrigin(const TPoint3D& newOrigin, TPose3D& pose)
{
        if (!contains(newOrigin))
                throw std::logic_error("Base point is not in the plane.");
        double normal[3];
        getUnitaryNormalVector(normal);
        CMatrixDouble44 AXIS;
        generateAxisBaseFromDirectionAndAxis(normal, 2, AXIS);
        for (size_t i = 0; i < 3; i++) AXIS.set_unsafe(i, 3, newOrigin[i]);
        pose.fromHomogeneousMatrix(AXIS);
}
TPlane::TPlane(const TPoint3D& p1, const TPoint3D& p2, const TPoint3D& p3)
{
        double dx1 = p2.x - p1.x;
        double dy1 = p2.y - p1.y;
        double dz1 = p2.z - p1.z;
        double dx2 = p3.x - p1.x;
        double dy2 = p3.y - p1.y;
        double dz2 = p3.z - p1.z;
        coefs[0] = dy1 * dz2 - dy2 * dz1;
        coefs[1] = dz1 * dx2 - dz2 * dx1;
        coefs[2] = dx1 * dy2 - dx2 * dy1;
        if (abs(coefs[0]) < getEpsilon() && abs(coefs[1]) < getEpsilon() &&
                abs(coefs[2]) < getEpsilon())
                throw logic_error("Points are linearly dependent");
        coefs[3] = -coefs[0] * p1.x - coefs[1] * p1.y - coefs[2] * p1.z;
}
TPlane::TPlane(const TPoint3D& p1, const TLine3D& r2)
{
        double dx1 = p1.x - r2.pBase.x;
        double dy1 = p1.y - r2.pBase.y;
        double dz1 = p1.z - r2.pBase.z;
        coefs[0] = dy1 * r2.director[2] - dz1 * r2.director[1];
        coefs[1] = dz1 * r2.director[0] - dx1 * r2.director[2];
        coefs[2] = dx1 * r2.director[1] - dy1 * r2.director[0];
        if (abs(coefs[0]) < getEpsilon() && abs(coefs[1]) < getEpsilon() &&
                abs(coefs[2]) < getEpsilon())
                throw logic_error("Point is contained in the line");
        coefs[3] = -coefs[0] * p1.x - coefs[1] * p1.y - coefs[2] * p1.z;
}
TPlane::TPlane(const TLine3D& r1, const TLine3D& r2)
{
        crossProduct3D(r1.director, r2.director, coefs);
        coefs[3] =
                -coefs[0] * r1.pBase.x - coefs[1] * r1.pBase.y - coefs[2] * r1.pBase.z;
        if (abs(coefs[0]) < getEpsilon() && abs(coefs[1]) < getEpsilon() &&
                abs(coefs[2]) < getEpsilon())
        {
                // Lines are parallel
                if (r1.contains(r2.pBase)) throw std::logic_error("Lines are the same");
                // Use a line's director vector and both pBase's difference to create
                // the plane.
                double d[3];
                for (size_t i = 0; i < 3; i++) d[i] = r1.pBase[i] - r2.pBase[i];
                crossProduct3D(r1.director, d, coefs);
                coefs[3] = -coefs[0] * r1.pBase.x - coefs[1] * r1.pBase.y -
                                   coefs[2] * r1.pBase.z;
        }
        else if (abs(evaluatePoint(r2.pBase)) >= getEpsilon())
                throw logic_error("Lines do not intersect");
}
 
template <class T>
inline void removeUnusedVertices(T& poly)
{
        size_t N = poly.size();
        if (N < 3) return;
        std::vector<size_t> unused;
        if (abs(mrpt::math::distance(poly[N - 1], poly[0]) +
                        mrpt::math::distance(poly[0], poly[1]) -
                        mrpt::math::distance(poly[N - 1], poly[1])) <
                mrpt::math::getEpsilon())
                unused.push_back(0);
        for (size_t i = 1; i < N - 1; i++)
                if (abs(mrpt::math::distance(poly[i - 1], poly[i]) +
                                mrpt::math::distance(poly[i], poly[i + 1]) -
                                mrpt::math::distance(poly[i - 1], poly[i + 1])) <
                        mrpt::math::getEpsilon())
                        unused.push_back(i);
        if (abs(mrpt::math::distance(poly[N - 2], poly[N - 1]) +
                        mrpt::math::distance(poly[N - 1], poly[0]) -
                        mrpt::math::distance(poly[N - 2], poly[0])) <
                mrpt::math::getEpsilon())
                unused.push_back(N - 1);
        unused.push_back(N);
        size_t diff = 1;
        for (size_t i = 0; i < unused.size() - 1; i++)
        {
                size_t last = unused[i + 1];
                for (size_t j = unused[i] + 1 - diff; j < last - diff; j++)
                        poly[j] = poly[j + diff];
        }
        poly.resize(N + 1 - unused.size());
}
template <class T>
inline void removeRepVertices(T& poly)
{
        size_t N = poly.size();
        if (N < 3) return;
        std::vector<size_t> rep;
        for (size_t i = 0; i < N - 1; i++)
                if (mrpt::math::distance(poly[i], poly[i + 1]) < getEpsilon())
                        rep.push_back(i);
        if (mrpt::math::distance(poly[N - 1], poly[0]) < getEpsilon())
                rep.push_back(N - 1);
        rep.push_back(N);
        size_t diff = 1;
        for (size_t i = 0; i < rep.size() - 1; i++)
        {
                size_t last = rep[i + 1];
                for (size_t j = rep[i] + 1 - diff; j < last - diff; j++)
                        poly[j] = poly[j + diff];
        }
        poly.resize(N + 1 - rep.size());
}
 
double TPolygon2D::distance(const TPoint2D& point) const
{
        if (contains(point)) return 0;
        std::vector<TSegment2D> sgs;
        getAsSegmentList(sgs);
 
        if (sgs.empty())
                THROW_EXCEPTION("Cannot compute distance to an empty polygon.");
 
        double distance = std::numeric_limits<double>::max();
 
        for (vector<TSegment2D>::const_iterator it = sgs.begin(); it != sgs.end();
                 ++it)
        {
                double d = (*it).distance(point);
                if (d < distance) distance = d;
        }
        return distance;
}
 
void TPolygon2D::getBoundingBox(
        TPoint2D& min_coords, TPoint2D& max_coords) const
{
        ASSERTMSG_(!this->empty(), "getBoundingBox() called on an empty polygon!");
        min_coords.x = min_coords.y = std::numeric_limits<double>::max();
        max_coords.x = max_coords.y = -std::numeric_limits<double>::max();
        for (size_t i = 0; i < size(); i++)
        {
                mrpt::keep_min(min_coords.x, (*this)[i].x);
                mrpt::keep_min(min_coords.y, (*this)[i].y);
                mrpt::keep_max(max_coords.x, (*this)[i].x);
                mrpt::keep_max(max_coords.y, (*this)[i].y);
        }
}
 
// isLeft(): tests if a point is Left|On|Right of an infinite line.
//    Input:  three points P0, P1, and P2
//    Return: >0 for P2 left of the line through P0 and P1
//            =0 for P2  on the line
//            <0 for P2  right of the line
//    See: Algorithm 1 "Area of Triangles and Polygons"
inline double isLeft(
        const mrpt::math::TPoint2D& P0, const mrpt::math::TPoint2D& P1,
        const mrpt::math::TPoint2D& P2)
{
        return ((P1.x - P0.x) * (P2.y - P0.y) - (P2.x - P0.x) * (P1.y - P0.y));
}
 
bool TPolygon2D::contains(const TPoint2D& P) const
{
        int wn = 0;  // the  winding number counter
 
        // loop through all edges of the polygon
        const size_t n = this->size();
        for (size_t i = 0; i < n; i++)  // edge from V[i] to  V[i+1]
        {
                if ((*this)[i].y <= P.y)
                {
                        // start y <= P.y
                        if ((*this)[(i + 1) % n].y > P.y)  // an upward crossing
                                if (isLeft((*this)[i], (*this)[(i + 1) % n], P) >
                                        0)  // P left of  edge
                                        ++wn;  // have  a valid up intersect
                }
                else
                {
                        // start y > P.y (no test needed)
                        if ((*this)[(i + 1) % n].y <= P.y)  // a downward crossing
                                if (isLeft((*this)[i], (*this)[(i + 1) % n], P) <
                                        0)  // P right of  edge
                                        --wn;  // have  a valid down intersect
                }
        }
 
        return wn != 0;
}
void TPolygon2D::getAsSegmentList(vector<TSegment2D>& v) const
{
        size_t N = size();
        v.resize(N);
        for (size_t i = 0; i < N - 1; i++)
                v[i] = TSegment2D(operator[](i), operator[](i + 1));
        v[N - 1] = TSegment2D(operator[](N - 1), operator[](0));
}
 
void TPolygon2D::generate3DObject(TPolygon3D& p) const
{
        p = TPolygon3D(*this);
}
// Auxiliar functor class to compute polygon's center
template <class T, int N>
class FAddPoint
{
   public:
        T& object;
        FAddPoint(T& o) : object(o)
        {
                for (size_t i = 0; i < N; i++) object[i] = 0.0;
        }
        void operator()(const T& o)
        {
                for (size_t i = 0; i < N; i++) object[i] += o[i];
        }
};
void TPolygon2D::getCenter(TPoint2D& p) const
{
        for_each(begin(), end(), FAddPoint<TPoint2D, 2>(p));
        size_t N = size();
        p.x /= N;
        p.y /= N;
}
bool TPolygon2D::isConvex() const
{
        size_t N = size();
        if (N <= 3) return true;
        vector<TSegment2D> sgms;
        getAsSegmentList(sgms);
        for (size_t i = 0; i < N; i++)
        {
                char s = 0;
                TLine2D l = TLine2D(sgms[i]);
                for (size_t j = 0; j < N; j++)
                {
                        double d = l.evaluatePoint(operator[](j));
                        if (abs(d) < getEpsilon())
                                continue;
                        else if (!s)
                                s = (d > 0) ? 1 : -1;
                        else if (s != ((d > 0) ? 1 : -1))
                                return false;
                }
        }
        return true;
}
void TPolygon2D::removeRepeatedVertices() { removeRepVertices(*this); }
void TPolygon2D::removeRedundantVertices()
{
        removeRepeatedVertices();
        removeUnusedVertices(*this);
}
void TPolygon2D::getPlotData(
        std::vector<double>& x, std::vector<double>& y) const
{
        size_t N = size();
        x.resize(N + 1);
        y.resize(N + 1);
        for (size_t i = 0; i < N; i++)
        {
                x[i] = operator[](i).x;
                y[i] = operator[](i).y;
        }
        x[N] = operator[](0).x;
        y[N] = operator[](0).y;
}
TPolygon2D::TPolygon2D(const TPolygon3D& p) : std::vector<TPoint2D>()
{
        size_t N = p.size();
        resize(N);
        for (size_t i = 0; i < N; i++) operator[](i) = TPoint2D(p[i]);
}
void TPolygon2D::createRegularPolygon(
        size_t numEdges, double radius, TPolygon2D& poly)
{
        if (numEdges < 3 || abs(radius) < getEpsilon())
                throw std::logic_error(
                        "Invalid arguments for regular polygon creations");
        poly.resize(numEdges);
        for (size_t i = 0; i < numEdges; i++)
        {
                double angle = i * M_PI * 2 / numEdges;
                poly[i] = TPoint2D(radius * cos(angle), radius * sin(angle));
        }
}
inline void TPolygon2D::createRegularPolygon(
        size_t numEdges, double radius, TPolygon2D& poly, const TPose2D& pose)
{
        createRegularPolygon(numEdges, radius, poly);
        for (size_t i = 0; i < numEdges; i++) poly[i] = pose.composePoint(poly[i]);
}
 
double TPolygon3D::distance(const TPoint3D& point) const
{
        TPlane pl;
        if (!getPlane(pl))
                throw std::logic_error("Polygon does not conform a plane");
        TPoint3D newPoint;
        TPolygon3D newPoly;
        TPose3D pose;
        pl.getAsPose3DForcingOrigin(operator[](0), pose);
        project3D(point, pose, newPoint);
        project3D(*this, pose, newPoly);
        double distance2D = TPolygon2D(newPoly).distance(TPoint2D(newPoint));
        return sqrt(newPoint.z * newPoint.z + distance2D * distance2D);
}
bool TPolygon3D::contains(const TPoint3D& point) const
{
        TPoint3D pMin, pMax;
        getPrismBounds(*this, pMin, pMax);
        if (point.x + getEpsilon() < pMin.x || point.y + getEpsilon() < pMin.y ||
                point.z + getEpsilon() < pMin.z || point.x > pMax.x + getEpsilon() ||
                point.y > pMax.y + getEpsilon() || point.z > pMax.z + getEpsilon())
                return false;
        TPlane plane;
        if (!getPlane(plane))
                throw std::logic_error("Polygon does not conform a plane");
        TPolygon3D projectedPoly;
        TPoint3D projectedPoint;
        TPose3D pose;
        // plane.getAsPose3DForcingOrigin(operator[](0),pose);
        plane.getAsPose3D(pose);
        CMatrixDouble44 P_inv;
        pose.getInverseHomogeneousMatrix(P_inv);
        pose.fromHomogeneousMatrix(P_inv);
        project3D(point, pose, projectedPoint);
        if (abs(projectedPoint.z) >= getEpsilon())
                return false;  // Point is not inside the polygon's plane.
        project3D(*this, pose, projectedPoly);
        return TPolygon2D(projectedPoly).contains(TPoint2D(projectedPoint));
}
void TPolygon3D::getAsSegmentList(vector<TSegment3D>& v) const
{
        size_t N = size();
        v.resize(N);
        for (size_t i = 0; i < N - 1; i++)
                v[i] = TSegment3D(operator[](i), operator[](i + 1));
        v[N - 1] = TSegment3D(operator[](N - 1), operator[](0));
}
bool TPolygon3D::getPlane(TPlane& p) const { return conformAPlane(*this, p); }
void TPolygon3D::getBestFittingPlane(TPlane& p) const
{
        getRegressionPlane(*this, p);
}
void TPolygon3D::getCenter(TPoint3D& p) const
{
        for_each(begin(), end(), FAddPoint<TPoint3D, 3>(p));
        size_t N = size();
        p.x /= N;
        p.y /= N;
        p.z /= N;
}
bool TPolygon3D::isSkew() const { return !mrpt::math::conformAPlane(*this); }
void TPolygon3D::removeRepeatedVertices() { removeRepVertices(*this); }
void TPolygon3D::removeRedundantVertices()
{
        removeRepeatedVertices();
        removeUnusedVertices(*this);
}
TPolygon3D::TPolygon3D(const TPolygon2D& p) : std::vector<TPoint3D>()
{
        size_t N = p.size();
        resize(N);
        for (size_t i = 0; i < N; i++) operator[](i) = p[i];
}
void TPolygon3D::createRegularPolygon(
        size_t numEdges, double radius, TPolygon3D& poly)
{
        if (numEdges < 3 || abs(radius) < getEpsilon())
                throw std::logic_error(
                        "Invalid arguments for regular polygon creations");
        poly.resize(numEdges);
        for (size_t i = 0; i < numEdges; i++)
        {
                double angle = i * 2 * M_PI / numEdges;
                poly[i] = TPoint3D(radius * cos(angle), radius * sin(angle), 0);
        }
}
inline void TPolygon3D::createRegularPolygon(
        size_t numEdges, double radius, TPolygon3D& poly, const TPose3D& pose)
{
        createRegularPolygon(numEdges, radius, poly);
        for (size_t i = 0; i < numEdges; i++) pose.composePoint(poly[i], poly[i]);
}
 
void TObject2D::generate3DObject(TObject3D& obj) const
{
        switch (type)
        {
                case GEOMETRIC_TYPE_POINT:
                        obj = TPoint3D(data.point);
                        break;
                case GEOMETRIC_TYPE_SEGMENT:
                        obj = TSegment3D(data.segment);
                        break;
                case GEOMETRIC_TYPE_LINE:
                        obj = TLine3D(data.line);
                        break;
                case GEOMETRIC_TYPE_POLYGON:
                        obj = TPolygon3D(*(data.polygon));
                        break;
                default:
                        obj = TObject3D();
                        break;
        }
}
void TObject2D::getPoints(
        const std::vector<TObject2D>& objs, std::vector<TPoint2D>& pnts)
{
        for (vector<TObject2D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isPoint()) pnts.push_back(it->data.point);
}
void TObject2D::getSegments(
        const std::vector<TObject2D>& objs, std::vector<TSegment2D>& sgms)
{
        for (vector<TObject2D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isSegment()) sgms.push_back(it->data.segment);
}
void TObject2D::getLines(
        const std::vector<TObject2D>& objs, std::vector<TLine2D>& lins)
{
        for (vector<TObject2D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isLine()) lins.push_back(it->data.line);
}
void TObject2D::getPolygons(
        const std::vector<TObject2D>& objs, std::vector<TPolygon2D>& polys)
{
        for (vector<TObject2D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isPolygon()) polys.push_back(*(it->data.polygon));
}
void TObject2D::getPoints(
        const std::vector<TObject2D>& objs, std::vector<TPoint2D>& pnts,
        std::vector<TObject2D>& remainder)
{
        for (vector<TObject2D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isPoint())
                        pnts.push_back(it->data.point);
                else
                        remainder.push_back(*it);
}
void TObject2D::getSegments(
        const std::vector<TObject2D>& objs, std::vector<TSegment2D>& sgms,
        std::vector<TObject2D>& remainder)
{
        for (vector<TObject2D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isSegment())
                        sgms.push_back(it->data.segment);
                else
                        remainder.push_back(*it);
}
void TObject2D::getLines(
        const std::vector<TObject2D>& objs, std::vector<TLine2D>& lins,
        std::vector<TObject2D>& remainder)
{
        for (vector<TObject2D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isLine())
                        lins.push_back(it->data.line);
                else
                        remainder.push_back(*it);
}
void TObject2D::getPolygons(
        const std::vector<TObject2D>& objs, std::vector<TPolygon2D>& polys,
        vector<TObject2D>& remainder)
{
        for (vector<TObject2D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isPolygon())
                        polys.push_back(*(it->data.polygon));
                else
                        remainder.push_back(*it);
}
void TObject3D::getPoints(
        const std::vector<TObject3D>& objs, std::vector<TPoint3D>& pnts)
{
        for (vector<TObject3D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isPoint()) pnts.push_back(it->data.point);
}
void TObject3D::getSegments(
        const std::vector<TObject3D>& objs, std::vector<TSegment3D>& sgms)
{
        for (vector<TObject3D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isSegment()) sgms.push_back(it->data.segment);
}
void TObject3D::getLines(
        const std::vector<TObject3D>& objs, std::vector<TLine3D>& lins)
{
        for (vector<TObject3D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isLine()) lins.push_back(it->data.line);
}
void TObject3D::getPlanes(
        const std::vector<TObject3D>& objs, std::vector<TPlane>& plns)
{
        for (vector<TObject3D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isPlane()) plns.push_back(it->data.plane);
}
void TObject3D::getPolygons(
        const std::vector<TObject3D>& objs, std::vector<TPolygon3D>& polys)
{
        for (vector<TObject3D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isPolygon()) polys.push_back(*(it->data.polygon));
}
void TObject3D::getPoints(
        const std::vector<TObject3D>& objs, std::vector<TPoint3D>& pnts,
        std::vector<TObject3D>& remainder)
{
        for (vector<TObject3D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isPoint())
                        pnts.push_back(it->data.point);
                else
                        remainder.push_back(*it);
}
void TObject3D::getSegments(
        const std::vector<TObject3D>& objs, std::vector<TSegment3D>& sgms,
        std::vector<TObject3D>& remainder)
{
        for (vector<TObject3D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isSegment())
                        sgms.push_back(it->data.segment);
                else
                        remainder.push_back(*it);
}
void TObject3D::getLines(
        const std::vector<TObject3D>& objs, std::vector<TLine3D>& lins,
        std::vector<TObject3D>& remainder)
{
        for (vector<TObject3D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isLine())
                        lins.push_back(it->data.line);
                else
                        remainder.push_back(*it);
}
void TObject3D::getPlanes(
        const std::vector<TObject3D>& objs, std::vector<TPlane>& plns,
        std::vector<TObject3D>& remainder)
{
        for (vector<TObject3D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isPlane())
                        plns.push_back(it->data.plane);
                else
                        remainder.push_back(*it);
}
void TObject3D::getPolygons(
        const std::vector<TObject3D>& objs, std::vector<TPolygon3D>& polys,
        vector<TObject3D>& remainder)
{
        for (vector<TObject3D>::const_iterator it = objs.begin(); it != objs.end();
                 ++it)
                if (it->isPolygon())
                        polys.push_back(*(it->data.polygon));
                else
                        remainder.push_back(*it);
}
 
CArchive& operator>>(CArchive& in, mrpt::math::TTwist2D& o)
{
        for (unsigned int i = 0; i < o.size(); i++) in >> o[i];
        return in;
}
CArchive& operator<<(CArchive& out, const mrpt::math::TTwist2D& o)
{
        for (unsigned int i = 0; i < o.size(); i++) out << o[i];
        return out;
}
 
CArchive& operator>>(CArchive& in, mrpt::math::TTwist3D& o)
{
        for (unsigned int i = 0; i < o.size(); i++) in >> o[i];
        return in;
}
CArchive& operator<<(CArchive& out, const mrpt::math::TTwist3D& o)
{
        for (unsigned int i = 0; i < o.size(); i++) out << o[i];
        return out;
}
 
CArchive& operator>>(CArchive& in, mrpt::math::TObject2D& o)
{
        uint16_t type;
        in >> type;
        switch (static_cast<unsigned char>(type))
        {
                case GEOMETRIC_TYPE_POINT:
                {
                        TPoint2D p;
                        in >> p;
                        o = p;
                }
                break;
                case GEOMETRIC_TYPE_SEGMENT:
                {
                        TSegment2D s;
                        in >> s;
                        o = s;
                }
                break;
                case GEOMETRIC_TYPE_LINE:
                {
                        TLine2D l;
                        in >> l;
                        o = l;
                }
                break;
                case GEOMETRIC_TYPE_POLYGON:
                {
                        TPolygon2D p;
                        in >> p;
                        o = p;
                }
                break;
                case GEOMETRIC_TYPE_UNDEFINED:
                {
                        o = TObject2D();
                }
                break;
                default:
                        throw std::logic_error(
                                "Unknown TObject2D type found while reading stream");
        }
        return in;
}
CArchive& operator<<(CArchive& out, const mrpt::math::TObject2D& o)
{
        out << static_cast<uint16_t>(o.getType());
        switch (o.getType())
        {
                case GEOMETRIC_TYPE_POINT:
                {
                        TPoint2D p;
                        o.getPoint(p);
                        return out << p;
                };
                case GEOMETRIC_TYPE_SEGMENT:
                {
                        TSegment2D s;
                        o.getSegment(s);
                        return out << s;
                };
                case GEOMETRIC_TYPE_LINE:
                {
                        TLine2D l;
                        o.getLine(l);
                        return out << l;
                };
                case GEOMETRIC_TYPE_POLYGON:
                {
                        TPolygon2D p;
                        o.getPolygon(p);
                        return out << p;
                };
        }
        return out;
}
 
CArchive& operator>>(CArchive& in, mrpt::math::TObject3D& o)
{
        uint16_t type;
        in >> type;
        switch (static_cast<unsigned char>(type))
        {
                case GEOMETRIC_TYPE_POINT:
                {
                        TPoint3D p;
                        in >> p;
                        o = p;
                }
                break;
                case GEOMETRIC_TYPE_SEGMENT:
                {
                        TSegment3D s;
                        in >> s;
                        o = s;
                }
                break;
                case GEOMETRIC_TYPE_LINE:
                {
                        TLine3D l;
                        in >> l;
                        o = l;
                }
                break;
                case GEOMETRIC_TYPE_PLANE:
                {
                        TPlane p;
                        in >> p;
                        o = p;
                }
                break;
                case GEOMETRIC_TYPE_POLYGON:
                {
                        TPolygon3D p;
                        in >> p;
                        o = p;
                }
                break;
                case GEOMETRIC_TYPE_UNDEFINED:
                {
                        o = TObject3D();
                }
                break;
                default:
                        throw std::logic_error(
                                "Unknown TObject3D type found while reading stream");
        }
        return in;
}
CArchive& operator<<(CArchive& out, const mrpt::math::TObject3D& o)
{
        out << static_cast<uint16_t>(o.getType());
        switch (o.getType())
        {
                case GEOMETRIC_TYPE_POINT:
                {
                        TPoint3D p;
                        o.getPoint(p);
                        return out << p;
                };
                case GEOMETRIC_TYPE_SEGMENT:
                {
                        TSegment3D s;
                        o.getSegment(s);
                        return out << s;
                };
                case GEOMETRIC_TYPE_LINE:
                {
                        TLine3D l;
                        o.getLine(l);
                        return out << l;
                };
                case GEOMETRIC_TYPE_PLANE:
                {
                        TPlane p;
                        o.getPlane(p);
                        return out << p;
                };
                case GEOMETRIC_TYPE_POLYGON:
                {
                        TPolygon3D p;
                        o.getPolygon(p);
                        return out << p;
                };
        }
        return out;
}
}  // namespace math
}  // namespace mrpt