CEllipsoid.cpp

Go to the documentation of this file.

/* +---------------------------------------------------------------------------+
   |                     Mobile Robot Programming Toolkit (MRPT)               |
   |                          http://www.mrpt.org/                             |
   |                                                                           |
   | Copyright (c) 2005-2017, 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 "opengl-precomp.h"  // Precompiled header


#include <mrpt/opengl/CEllipsoid.h>
#include <mrpt/math/CMatrix.h>
#include <mrpt/math/geometry.h>
#include <mrpt/math/matrix_serialization.h>
#include <mrpt/utils/CStream.h>

#include "opengl_internals.h"

using namespace mrpt;
using namespace mrpt::opengl;
using namespace mrpt::utils;
using namespace mrpt::math;
using namespace std;

IMPLEMENTS_SERIALIZABLE( CEllipsoid, CRenderizableDisplayList, mrpt::opengl )

/*---------------------------------------------------------------
                                                        render
  ---------------------------------------------------------------*/
void   CEllipsoid::render_dl() const
{
#if MRPT_HAS_OPENGL_GLUT
        MRPT_START

    const size_t dim = m_cov.getColCount();

        if(m_eigVal(0,0) != 0.0 && m_eigVal(1,1) != 0.0 && (dim==2 || m_eigVal(2,2) != 0.0) && m_quantiles!=0.0)
        {
                glEnable(GL_BLEND);
                checkOpenGLError();
                glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
                checkOpenGLError();
                glLineWidth(m_lineWidth);
                checkOpenGLError();

                if (dim==2)
                {
                        glDisable(GL_LIGHTING);  // Disable lights when drawing lines

                        // ---------------------
                        //     2D ellipse
                        // ---------------------

                        /* Equivalent MATLAB code:
                         *
                         * q=1;
                         * [vec val]=eig(C);
                         * M=(q*val*vec)';
                         * R=M*[x;y];
                         * xx=R(1,:);yy=R(2,:);
                         * plot(xx,yy), axis equal;
                         */

                        double                  ang;
                        unsigned int    i;

                        // Compute the new vectors for the ellipsoid:
                        CMatrixDouble   M;
                        M.noalias() = double(m_quantiles) * m_eigVal * m_eigVec.adjoint();

                        glBegin( GL_LINE_LOOP );

                        // Compute the points of the 2D ellipse:
                        for (i=0,ang=0;i<m_2D_segments;i++,ang+= (M_2PI/m_2D_segments))
                        {
                                double ccos = cos(ang);
                                double ssin = sin(ang);

                                const float x = ccos * M.get_unsafe(0,0) + ssin * M.get_unsafe(1,0);
                                const float y = ccos * M.get_unsafe(0,1) + ssin * M.get_unsafe(1,1);

                                glVertex2f( x,y );
                        } // end for points on ellipse

                        glEnd();

                        // 2D: Save bounding box:
                        const double max_radius = m_quantiles * std::max( m_eigVal(0,0), m_eigVal(1,1) );
                        m_bb_min = mrpt::math::TPoint3D(-max_radius,-max_radius, 0);
                        m_bb_max = mrpt::math::TPoint3D(max_radius,max_radius, 0);
                        // Convert to coordinates of my parent:
                        m_pose.composePoint(m_bb_min, m_bb_min);
                        m_pose.composePoint(m_bb_max, m_bb_max);

                        glEnable(GL_LIGHTING);
                }
                else
                {
                        // ---------------------
                        //    3D ellipsoid
                        // ---------------------
                        GLfloat         mat[16];

                        //  A homogeneous transformation matrix, in this order:
                        //
                        //     0  4  8  12
                        //     1  5  9  13
                        //     2  6  10 14
                        //     3  7  11 15
                        //
                        mat[3] = mat[7] = mat[11] = 0;
                        mat[15] = 1;
                        mat[12] = mat[13] = mat[14] = 0;

                        mat[0] = m_eigVec(0,0); mat[1] = m_eigVec(1,0); mat[2] = m_eigVec(2,0); // New X-axis
                        mat[4] = m_eigVec(0,1); mat[5] = m_eigVec(1,1); mat[6] = m_eigVec(2,1); // New X-axis
                        mat[8] = m_eigVec(0,2); mat[9] = m_eigVec(1,2); mat[10] = m_eigVec(2,2);        // New X-axis

                        GLUquadricObj   *obj = gluNewQuadric();
                        checkOpenGLError();

                        if (!m_drawSolid3D) glDisable(GL_LIGHTING);  // Disable lights when drawing lines

                        gluQuadricDrawStyle( obj, m_drawSolid3D ? GLU_FILL : GLU_LINE);

                        glPushMatrix();
                        glMultMatrixf( mat );
                        glScalef(m_eigVal(0,0)*m_quantiles,m_eigVal(1,1)*m_quantiles,m_eigVal(2,2)*m_quantiles);

                        gluSphere( obj, 1,m_3D_segments,m_3D_segments);
                        checkOpenGLError();

                        glPopMatrix();

                        gluDeleteQuadric(obj);
                        checkOpenGLError();

                        // 3D: Save bounding box:
                        const double max_radius = m_quantiles * std::max( m_eigVal(0,0), std::max(m_eigVal(1,1), m_eigVal(2,2) ) );
                        m_bb_min = mrpt::math::TPoint3D(-max_radius,-max_radius, 0);
                        m_bb_max = mrpt::math::TPoint3D(max_radius,max_radius, 0);
                        // Convert to coordinates of my parent:
                        m_pose.composePoint(m_bb_min, m_bb_min);
                        m_pose.composePoint(m_bb_max, m_bb_max);
                }


                glDisable(GL_BLEND);

                glEnable(GL_LIGHTING);
        }
        MRPT_END_WITH_CLEAN_UP( \
                cout << "Covariance matrix leading to error is:" << endl << m_cov << endl; \
        );
#endif
}

/*---------------------------------------------------------------
   Implements the writing to a CStream capability of
     CSerializable objects
  ---------------------------------------------------------------*/
void  CEllipsoid::writeToStream(mrpt::utils::CStream &out,int *version) const
{
        if (version)
                *version = 1;
        else
        {
                writeToStreamRender(out);
                out << m_cov << m_drawSolid3D << m_quantiles << (uint32_t)m_2D_segments << (uint32_t)m_3D_segments << m_lineWidth;
        }
}

/*---------------------------------------------------------------
        Implements the reading from a CStream capability of
                CSerializable objects
  ---------------------------------------------------------------*/
void  CEllipsoid::readFromStream(mrpt::utils::CStream &in,int version)
{
        switch(version)
        {
        case 0:
        case 1:
                {
                        uint32_t        i;
                        readFromStreamRender(in);
                        if (version==0)
                        {
                                CMatrix c;
                                in >> c; m_cov = c.cast<double>();
                        }
                        else
                        {
                                in >> m_cov;
                        }

                        in >> m_drawSolid3D >> m_quantiles;
                        in >> i; m_2D_segments = i;
                        in >> i; m_3D_segments = i;
                        in >> m_lineWidth;

                        // Update cov. matrix cache:
                        m_prevComputedCov = m_cov;
                        setCovMatrix(m_cov);

                } break;
        default:
                MRPT_THROW_UNKNOWN_SERIALIZATION_VERSION(version)

        };
        CRenderizableDisplayList::notifyChange();
}

bool quickSolveEqn(double a,double b_2,double c,double &t)      {
        double delta=square(b_2)-a*c;
        if (delta==0) return (t=-b_2/a)>=0;
        else if (delta>0)       {
                delta=sqrt(delta);
                if ((t=(-b_2-delta)/a)>=0) return true;
                else return (t=(-b_2+delta)/a)>=0;
        }       else return false;
}

bool CEllipsoid::traceRay(const mrpt::poses::CPose3D &o,double &dist) const     {
        if (m_cov.getRowCount()!=3) return false;
        TLine3D lin,lin2;
        createFromPoseX(o-this->m_pose,lin);
        lin.unitarize();        //By adding this line, distance from any point of the line to its base is exactly equal to the "t".
        for (size_t i=0;i<3;i++)        {
                lin2.pBase[i]=0;
                lin2.director[i]=0;
                for (size_t j=0;j<3;j++)        {
                        double vji=m_eigVec(j,i);
                        lin2.pBase[i]+=vji*lin.pBase[j];
                        lin2.director[i]+=vji*lin.director[j];
                }
        }
        double a=0,b_2=0,c=-square(m_quantiles);
        for (size_t i=0;i<3;i++)        {
                double ev=m_eigVal(i,i);
                a+=square(lin2.director[i]/ev);
                b_2+=lin2.director[i]*lin2.pBase[i]/square(ev);
                c+=square(lin2.pBase[i]/ev);
        }
        return quickSolveEqn(a,b_2,c,dist);
}

void CEllipsoid::setCovMatrix( const mrpt::math::CMatrixDouble &m, int resizeToSize)
{
        MRPT_START

        ASSERT_( m.getColCount() == m.getRowCount() );
        ASSERT_( size(m,1)==2 || size(m,1)==3 || (resizeToSize>0 && (resizeToSize==2 || resizeToSize==3)));

        m_cov = m;
        if (resizeToSize>0 && resizeToSize<(int)size(m,1))
                m_cov.setSize(resizeToSize,resizeToSize);

        if (m_cov==m_prevComputedCov)
                return; // Done.

        CRenderizableDisplayList::notifyChange();

        // Handle the special case of an ellipsoid of volume = 0
        const double d=m_cov.det();
        if (d==0 || d!=d) // Note: "d!=d" is a great test for invalid numbers, don't remove!
        {
                // All zeros:
                m_prevComputedCov = m_cov;
                m_eigVec.zeros(3,3);
                m_eigVal.zeros(3,3);
        }
        else
        {
                // Not null matrix: compute the eigen-vectors & values:
                m_prevComputedCov = m_cov;
                if (m_cov.eigenVectors(m_eigVec,m_eigVal)) {
                        m_eigVal = m_eigVal.array().sqrt().matrix();
                        // Do the scale at render to avoid recomputing the m_eigVal for different m_quantiles
                } else {
                        m_eigVec.zeros(3,3);
                        m_eigVal.zeros(3,3);
                }
        }


        MRPT_END
}

void CEllipsoid::setCovMatrix( const mrpt::math::CMatrixFloat &m, int resizeToSize)
{
        CRenderizableDisplayList::notifyChange();
        setCovMatrix( CMatrixDouble(m), resizeToSize);
}

/** Evaluates the bounding box of this object (including possible children) in the coordinate frame of the object parent. */
void CEllipsoid::getBoundingBox(mrpt::math::TPoint3D &bb_min, mrpt::math::TPoint3D &bb_max) const
{
        bb_min = m_bb_min;
        bb_max = m_bb_max;
}