Commit 570f0b11 authored by Andreas Schwarzkopf's avatar Andreas Schwarzkopf

[ADD ##371] Surface detection approach of Lari/Habib 2014

 * This approach describes the surface detection of Zahra Lari and Ayman Habib (2014).
 * Paper: "An adaptive approach for the segmentation and extraction of planar
 * andlinear/cylindrical features from laser scanning data", ISPRS Journal of
 * Photogrammetry and Remote Sensing, 2014
 *
 * Eventually the process covers also the linear and cylindrical feature segmentation
 * but in our requirements only the planar segmentation is required. Steps:
 *   - Original laser point cloud as input data
 *   - PCA-based classification of planar features
 *   - Selection of appropriate representation models for the detected planar features
 *   - Local point density estimation along the detected planar features
 *   - Precise estimation of the planar features segmentation attributes
 *   - Parameter-domain segmentation and extraction of planar features
 *   - Boundary tracking for resolving the parameter-domain segmentation ambiguities of
 *     planar features
 *   - Extracted planar features in spatial domain
 *
 * Estimations (Based on the paper because I haven't completely impelemented yet):
 *   - Fair segmentation results
 *   - Relatively slow
 *
 * Implementation progress:
 *   - Done to the step "Precise estimation of the planar features segmentation
 *     attributes".
 *
 *
 * Other updates:
 *   - Renamed module "Pooints - Crop" to "Points - Transform". It optained the feature
 *     to translate and to rotate the point set.
 *   - Wew code refactoring.
parent b755bfa0
......@@ -33,7 +33,7 @@
#include "buildingsDetection/WMBuildingsDetection.h"
#include "buildingsDetectionByPCA/WMBuildingsDetectionByPCA.h"
#include "elevationImageExport/WMElevationImageExport.h"
#include "pointsCrop/WMPointsCrop.h"
#include "pointsTransform/WMPointsTransform.h"
#include "pointsCutOutliers/WMPointsCutOutliers.h"
#include "pointsGroupSelector/WMPointsGroupSelector.h"
#include "readLAS/WMReadLAS.h"
......@@ -58,9 +58,9 @@ extern "C" void WLoadModule( WModuleList& m ) // NOLINT
m.push_back( boost::shared_ptr< WModule >( new WMBuildingsDetection ) );
m.push_back( boost::shared_ptr< WModule >( new WMBuildingsDetectionByPCA ) );
m.push_back( boost::shared_ptr< WModule >( new WMElevationImageExport ) );
m.push_back( boost::shared_ptr< WModule >( new WMPointsCrop ) );
m.push_back( boost::shared_ptr< WModule >( new WMPointsCutOutliers ) );
m.push_back( boost::shared_ptr< WModule >( new WMPointsGroupSelector ) );
m.push_back( boost::shared_ptr< WModule >( new WMPointsTransform ) );
m.push_back( boost::shared_ptr< WModule >( new WMReadLAS ) );
m.push_back( boost::shared_ptr< WModule >( new WMSurfaceDetectionByLari ) );
m.push_back( boost::shared_ptr< WModule >( new WMSurfaceDetectionByPCL ) );
......
......@@ -28,10 +28,10 @@
#include <vector>
#include "core/graphicsEngine/WTriangleMesh.h"
#include "core/dataHandler/WDataSetPoints.h"
#include "../datastructures/octree/WOctNode.h"
#include "../datastructures/quadtree/WQuadNode.h"
#include "../datastructures/quadtree/WQuadTree.h"
#include "../datastructures/octree/WOctree.h"
#include "../common/datastructures/octree/WOctNode.h"
#include "../common/datastructures/quadtree/WQuadNode.h"
#include "../common/datastructures/quadtree/WQuadTree.h"
#include "../common/datastructures/octree/WOctree.h"
/**
* Class that detects buildings using the WDataSetPoints
......
......@@ -42,7 +42,7 @@
#include "WMBuildingsDetection.xpm"
#include "WMBuildingsDetection.h"
#include "WBuildingDetector.h"
#include "../datastructures/octree/WOctree.h"
#include "../common/datastructures/octree/WOctree.h"
// This line is needed by the module loader to actually find your module.
//W_LOADABLE_MODULE( WMBuildingsDetection )
......
......@@ -42,10 +42,10 @@
#include <osg/ShapeDrawable>
#include <osg/Geode>
#include "core/dataHandler/WDataSetPoints.h"
#include "../datastructures/octree/WOctree.h"
#include "../datastructures/quadtree/WQuadTree.h"
#include "../common/datastructures/octree/WOctree.h"
#include "../common/datastructures/quadtree/WQuadTree.h"
#include "../datastructures/WDataSetPointsGrouped.h"
#include "../common/datastructures/WDataSetPointsGrouped.h"
......
......@@ -42,7 +42,7 @@
#include "WMBuildingsDetectionByPCA.xpm"
#include "WMBuildingsDetectionByPCA.h"
#include "WPCADetector.h"
#include "../datastructures/octree/WOctree.h"
#include "../common/datastructures/octree/WOctree.h"
// This line is needed by the module loader to actually find your module.
//W_LOADABLE_MODULE( WMBuildingsDetectionByPCA )
......
......@@ -42,10 +42,10 @@
#include <osg/ShapeDrawable>
#include <osg/Geode>
#include "core/dataHandler/WDataSetPoints.h"
#include "../datastructures/octree/WOctree.h"
#include "../datastructures/quadtree/WQuadTree.h"
#include "../common/datastructures/octree/WOctree.h"
#include "../common/datastructures/quadtree/WQuadTree.h"
#include "../datastructures/WDataSetPointsGrouped.h"
#include "../common/datastructures/WDataSetPointsGrouped.h"
......
......@@ -28,10 +28,10 @@
#include <vector>
#include "core/graphicsEngine/WTriangleMesh.h"
#include "core/dataHandler/WDataSetPoints.h"
#include "../datastructures/octree/WOctNode.h"
#include "../datastructures/quadtree/WQuadNode.h"
#include "../datastructures/quadtree/WQuadTree.h"
#include "../datastructures/octree/WOctree.h"
#include "../common/datastructures/octree/WOctNode.h"
#include "../common/datastructures/quadtree/WQuadNode.h"
#include "../common/datastructures/quadtree/WQuadTree.h"
#include "../common/datastructures/octree/WOctree.h"
#include "structure/WPcaDetectOctNode.h"
#include "core/common/math/principalComponentAnalysis/WPrincipalComponentAnalysis.h"
#include "core/common/WProgress.h"
......
......@@ -26,7 +26,7 @@
#ifndef WPCADETECTOCTNODE_H
#define WPCADETECTOCTNODE_H
#include <vector>
#include "../../datastructures/octree/WOctNode.h"
#include "../../common/datastructures/octree/WOctNode.h"
#include "core/common/math/linearAlgebra/WPosition.h"
#include "core/dataHandler/WDataSetPoints.h"
......
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV-Leipzig and CNCF-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#include <vector>
#include "WLeastSquares.h"
WLeastSquares::WLeastSquares()
{
m_positions = 0;
m_dimensions = 3;
m_verticalDimension = 2;
}
WLeastSquares::WLeastSquares( size_t dimensions )
{
m_positions = 0;
m_dimensions = dimensions;
m_verticalDimension = 2;
}
WLeastSquares::~WLeastSquares()
{
}
void WLeastSquares::analyzeData( vector<WPosition>* data )
{
m_positions = data;
calculatePerpendicularDimension();
clearMatrices();
calculateMatrices();
calculateHessescheNormalForm();
}
vector<double> WLeastSquares::getHessescheNormalForm()
{
return m_hessescheNormalForm;
}
void WLeastSquares::calculatePerpendicularDimension()
{
WPrincipalComponentAnalysis pca;
pca.analyzeData( m_positions );
vector<double> eigenValues = pca.getEigenValues();
double smallestEigenValue = eigenValues[0];
size_t perpendicularEigenVector = 0;
for( size_t index = 0; index < eigenValues.size(); index++ )
if( eigenValues[index] < smallestEigenValue )
{
smallestEigenValue = eigenValues[index];
perpendicularEigenVector = index;
}
m_verticalDimension = 0;
WVector3d eigenVector = pca.getDirections()[perpendicularEigenVector];
double biggestValue = abs( eigenVector[0] );
for( size_t index = 0; index < eigenVector.size(); index++ )
if( abs( eigenVector[index] ) > biggestValue )
{
biggestValue = abs( eigenVector[index] );
m_verticalDimension = index;
}
}
void WLeastSquares::clearMatrices()
{
MatrixXd newMatrixA( m_dimensions, m_dimensions );
m_matrixA = newMatrixA;
MatrixXd newMatrixB( m_dimensions, 1 );
m_matrixB = newMatrixB;
for( size_t row = 0; row < m_dimensions; row++ )
{
for( size_t col = 0; col < m_dimensions; col++ )
m_matrixA( row, col ) = 0.0;
m_matrixB( row, 0 ) = 0.0;
}
}
void WLeastSquares::calculateMatrices()
{
for( size_t index = 0; index < m_positions->size(); index++ )
{
WPosition pos = m_positions->at( index );
for( size_t row = 0; row < m_dimensions; row++ )
{
size_t rowMat = row <= m_verticalDimension ?row :row - 1;
if( row == m_verticalDimension ) rowMat = m_dimensions - 1;
for( size_t col = 0; col < m_dimensions; col++ )
{
size_t colMat = col <= m_verticalDimension ?col :col - 1;
if( col == m_verticalDimension )
colMat = m_dimensions - 1;
double resultA = 1.0;
if( row != m_verticalDimension )
resultA *= pos[row];
if( col != m_verticalDimension )
resultA *= pos[col];
m_matrixA( rowMat, colMat ) += resultA;
}
double resultB = pos[m_verticalDimension];
if( row != m_verticalDimension )
resultB *= pos[row];
m_matrixB( rowMat, 0 ) += resultB;
}
}
}
void WLeastSquares::calculateHessescheNormalForm()
{
MatrixXd matrixX = m_matrixA.inverse()*m_matrixB;
m_hessescheNormalForm.reserve( m_dimensions + 1 );
m_hessescheNormalForm.resize( m_dimensions + 1 );
for( size_t index = 0; index < m_hessescheNormalForm.size(); index++ )
{
size_t indexMat = index <= m_verticalDimension ?index : index - 1;
m_hessescheNormalForm[index] =
index != m_verticalDimension
?-matrixX( indexMat, 0 ) :1;
}
double sumSquared = 0;
for( size_t index = 0; index < m_dimensions; index++ )
sumSquared += pow( m_hessescheNormalForm[index], 2.0 );
sumSquared = pow( sumSquared, 0.5 );
for( size_t index = 0; index < m_dimensions; index++ )
m_hessescheNormalForm[index] /= sumSquared;
}
vector<double> WLeastSquares::getParametersXYZ0()
{
vector<double> parameters;
double squaredSum = 0;
for( size_t index = 0; index < m_dimensions; index++ )
squaredSum += pow( m_hessescheNormalForm[index], 2.0 );
for(size_t index = 0; index < m_dimensions; index++)
parameters.push_back( -m_hessescheNormalForm[index] *
m_hessescheNormalForm[ m_dimensions] / squaredSum );
return parameters;
}
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV-Leipzig and CNCF-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#ifndef WLEASTSQUARES_H
#define WLEASTSQUARES_H
#include <vector>
#include "core/common/math/linearAlgebra/WPosition.h"
#include "core/common/math/linearAlgebra/WVectorFixed.h"
#include "core/common/math/principalComponentAnalysis/WPrincipalComponentAnalysis.h"
#include <Eigen/Dense>
using std::abs;
using std::cout;
using std::endl;
using std::pow;
using std::vector;
/**
* Class that does the least sqhares adjustment on arbitrary points with n dimensions.
* The procedure finds a plane with the least errors sum. Each error is multiplied with
* itself (count by the power of two) to calculate the entire error amount. The class
* calculates a plane with the smallest total error.
*/
class WLeastSquares
{
public:
WLeastSquares();
/**
* Establishs a new least squares adjustment instance using an arbitraryly
* dimensional space.
* \param dimensions Dimension definition of the space for the calculation.
*/
explicit WLeastSquares( size_t dimensions );
virtual ~WLeastSquares();
/**
* Launchs the least squares adjustment.
* \param data Point data to analyze the best fitted plane for.
*/
void analyzeData( vector<WPosition>* data );
/**
* Returns the plane formula for the best fitted plane.
* \return The hessesche Normal Form of the best fitted plane. First n numbers (by
* the dimensions coordinate count) represent the normal vector of the best
* fitted plane. The last one is the perpendicular euclidian distance to the
* coordinate system orign.
*/
vector<double> getHessescheNormalForm();
/**
* Returns the parameters X_0, Y_0 and Z_0 etc. of the best fitted plane in an n
* dimensional space.
* \return Parameter space coordinates corresponding to Lari/Habib 2014.
*/
vector<double> getParametersXYZ0();
private:
/**
* Returns the dimension with the smallest extent
* \return Dimension with the biggest eigen vector coordonate extent of the smallest Eigen
* Value.
*/
void calculatePerpendicularDimension();
/**
* Erases the matrices required during the calculation.
*/
void clearMatrices();
/**
* The algorithm solves the matrix formula A*x=B. This method calculates A and B.
*/
void calculateMatrices();
/**
* The algorithm solves the matrix formula A*x=B. This method at first calculates
* the matrix x and then the Hessesche Normal Form of the best fitted plane.
*/
void calculateHessescheNormalForm();
/**
* Dimension definition of the coordinate system.
*/
size_t m_dimensions;
/**
* Points to calculate the best fitted plane for.
*/
vector<WPosition>* m_positions;
/**
* Space for the calculated Hessesche Normal Form.
*/
vector<double> m_hessescheNormalForm;
/**
* Dimension with the biggest eigen vector coordonate extent of the smallest Eigen Value.
*/
size_t m_verticalDimension;
/**
* Matrix A of the least squares adjustment,
*/
MatrixXd m_matrixA;
/**
* Matrix A of the least squares adjustment,
*/
MatrixXd m_matrixB;
};
#endif // WLEASTSQUARES_H
......@@ -43,7 +43,7 @@
#include <osg/Geode>
#include "core/dataHandler/WDataSetPoints.h"
#include "../datastructures/WDataSetPointsGrouped.h"
#include "../common/datastructures/WDataSetPointsGrouped.h"
......@@ -69,10 +69,10 @@
#include "bitmapImage/WBmpImage.h"
#include "bitmapImage/WBmpSaver.h"
#include "../datastructures/quadtree/WQuadTree.h"
#include "../datastructures/quadtree/WQuadNode.h"
#include "../datastructures/octree/WOctree.h"
#include "../datastructures/octree/WOctNode.h"
#include "../common/datastructures/quadtree/WQuadTree.h"
#include "../common/datastructures/quadtree/WQuadNode.h"
#include "../common/datastructures/octree/WOctree.h"
#include "../common/datastructures/octree/WOctNode.h"
/**
* Tool that outlines an elevation imate to a triangle mesh.
......
......@@ -42,7 +42,7 @@
#include <osg/ShapeDrawable>
#include <osg/Geode>
#include "core/dataHandler/WDataSetPoints.h"
#include "../datastructures/quadtree/WQuadTree.h"
#include "../common/datastructures/quadtree/WQuadTree.h"
......@@ -70,7 +70,7 @@
#include "bitmapImage/WBmpSaver.h"
#include "WElevationImageOutliner.h"
#include "../datastructures/WDataSetPointsGrouped.h"
#include "../common/datastructures/WDataSetPointsGrouped.h"
// forward declarations to reduce compile dependencies
template< class T > class WModuleInputData;
......
......@@ -28,10 +28,10 @@
#include <vector>
#include "../../datastructures/quadtree/WQuadNode.h"
#include "../../datastructures/quadtree/WQuadTree.h"
#include "../../datastructures/octree/WOctree.h"
#include "../../datastructures/WDataSetPointsGrouped.h"
#include "../../common/datastructures/quadtree/WQuadNode.h"
#include "../../common/datastructures/quadtree/WQuadTree.h"
#include "../../common/datastructures/octree/WOctree.h"
#include "../../common/datastructures/WDataSetPointsGrouped.h"
/**
* Image object. Currently it's used for saving bmp files.
......
......@@ -26,8 +26,8 @@
#include <vector>
#include "WCutOutliersDeamon.h"
#include "../datastructures/octree/WOctNode.h"
#include "../datastructures/octree/WOctree.h"
#include "../common/datastructures/octree/WOctNode.h"
#include "../common/datastructures/octree/WOctree.h"
WCutOutliersDeamon::WCutOutliersDeamon()
{
......
......@@ -28,7 +28,7 @@
#include <vector>
#include "core/graphicsEngine/WTriangleMesh.h"
#include "core/dataHandler/WDataSetPoints.h"
#include "../datastructures/octree/WOctNode.h"
#include "../common/datastructures/octree/WOctNode.h"
/**
* This is an outliers cut algorithm it simply groups all the points in cube groups.
......
......@@ -41,7 +41,7 @@
#include "core/kernel/WModuleInputData.h"
#include "WMPointsCutOutliers.xpm"
#include "WMPointsCutOutliers.h"
#include "../datastructures/octree/WOctree.h"
#include "../common/datastructures/octree/WOctree.h"
#include "WCutOutliersDeamon.h"
// This line is needed by the module loader to actually find your module.
......
......@@ -42,7 +42,7 @@
#include <osg/ShapeDrawable>
#include <osg/Geode>
#include "core/dataHandler/WDataSetPoints.h"
#include "../datastructures/octree/WOctree.h"
#include "../common/datastructures/octree/WOctree.h"
......
......@@ -42,7 +42,7 @@
#include "WMPointsGroupSelector.xpm"
#include "WMPointsGroupSelector.h"
#include "WVoxelOutliner.h"
#include "../datastructures/octree/WOctree.h"
#include "../common/datastructures/octree/WOctree.h"
// This line is needed by the module loader to actually find your module.
//W_LOADABLE_MODULE( WMPointsGroupSelector )
......
......@@ -41,9 +41,9 @@
#include <osg/ShapeDrawable>
#include <osg/Geode>
#include "../datastructures/octree/WOctree.h"
#include "../common/datastructures/octree/WOctree.h"
#include "../datastructures/WDataSetPointsGrouped.h"
#include "../common/datastructures/WDataSetPointsGrouped.h"
#include "core/dataHandler/WDataSetPoints.h"
......
......@@ -28,8 +28,8 @@
#include <vector>
#include "core/graphicsEngine/WTriangleMesh.h"
#include "core/dataHandler/WDataSetPoints.h"
#include "../datastructures/octree/WOctNode.h"
#include "../datastructures/octree/WOctree.h"
#include "../common/datastructures/octree/WOctNode.h"
#include "../common/datastructures/octree/WOctree.h"
/**
* Tool to draw an octree to a WTriangle mesh in order to e. g. display it using the plugin
......
......@@ -39,44 +39,44 @@
#include "core/graphicsEngine/WGEManagedGroupNode.h"
#include "core/kernel/WKernel.h"
#include "core/kernel/WModuleInputData.h"
#include "WMPointsCrop.xpm"
#include "WMPointsCrop.h"
#include "../datastructures/octree/WOctree.h"
#include "WMPointsTransform.xpm"
#include "WMPointsTransform.h"
#include "../common/datastructures/octree/WOctree.h"
// This line is needed by the module loader to actually find your module.
//W_LOADABLE_MODULE( WMPointsCrop )
//W_LOADABLE_MODULE( WMPointsTransform )
//TODO(aschwarzkopf): Reenable above after solving the toolbox problem
WMPointsCrop::WMPointsCrop():
WMPointsTransform::WMPointsTransform():
WModule(),
m_propCondition( new WCondition() )
{
}
WMPointsCrop::~WMPointsCrop()
WMPointsTransform::~WMPointsTransform()
{
}
boost::shared_ptr< WModule > WMPointsCrop::factory() const
boost::shared_ptr< WModule > WMPointsTransform::factory() const
{
return boost::shared_ptr< WModule >( new WMPointsCrop() );
return boost::shared_ptr< WModule >( new WMPointsTransform() );
}
const char** WMPointsCrop::getXPMIcon() const
const char** WMPointsTransform::getXPMIcon() const
{
return WMPointsCrop_xpm;
return WMPointsTransform_xpm;
}
const std::string WMPointsCrop::getName() const
const std::string WMPointsTransform::getName() const
{
return "Points - Crop";
return "Points - Transform";
}
const std::string WMPointsCrop::getDescription() const
const std::string WMPointsTransform::getDescription() const
{
return "Crops point data to a selection.";
}
void WMPointsCrop::connectors()
void WMPointsTransform::connectors()
{
m_input = WModuleInputData< WDataSetPoints >::createAndAdd( shared_from_this(), "input", "The mesh to display" );
......@@ -88,26 +88,58 @@ void WMPointsCrop::connectors()
WModule::connectors();
}