demoSegmentation.cpp

/// example showing how to use opalsSegmentation
/// and access its results based on the segment manager
 
#include "opals/Exception.hpp"
#include "opals/ParamList.hpp"
#include "opals/ModuleDeleter.hpp"
#include "opals/ISegmentation.hpp"
#include "DM/IDatamanager.hpp"
#include "DM/IPolyline.hpp"
 
#include <iostream>
#include <iomanip>
 
#include <memory>
 
#include "opals/IImport.hpp"
#include "opals/INormals.hpp"
 
 
void output_addinfo(const DM::IAddInfo &info)
{
  //output all attribute values within an addinfo object
 
  std::cout << "attributes: ";
  if (!info.columns())
    std::cout << "NONE";
  else
  {
    for (unsigned i = 0; i < info.columns(); i++)
    {
      if (i > 0)
        std::cout << ", ";
      std::cout << info.name(i) << "=";
      if (info.isNull(i))
        std::cout << "NULL";
      else
      {
        switch (info.type(i))
        {
        case DM::ColumnType::int32:
          std::cout << info.getInt32(i);
          break;
        case DM::ColumnType::uint32:
          std::cout << info.getUInt32(i);
          break;
        case DM::ColumnType::int8:
          std::cout << info.getInt8(i);
          break;
        case DM::ColumnType::uint8:
          std::cout << info.getUInt8(i);
          break;
        case DM::ColumnType::int16:
          std::cout << info.getInt16(i);
          break;
        case DM::ColumnType::uint16:
          std::cout << info.getUInt16(i);
          break;
        case DM::ColumnType::float_:
          std::cout << info.getFloat(i);
          break;
        case DM::ColumnType::double_:
          std::cout << info.getDouble(i);
          break;
        case DM::ColumnType::int64:
          std::cout << info.getInt64(i);
          break;
        case DM::ColumnType::string:
          std::cout << info.getCStr(i);
          break;
        }
      }
    }
  }
  std::cout << std::endl;
}
 
 
void preprocessing(opals::Path name)
{
  std::shared_ptr<opals::IImport> imp(opals::IImport::New(), opals::ModuleDeleter());
  opals::Vector<opals::VectorOrRasterFile> infiles{1, name + ".laz"};
  imp->opts().inFile().set(infiles);
  imp->opts().tileSize().set(100);
  imp->opts().filter().set("generic[z>265]");
  imp->run();
  imp.reset();
 
  std::shared_ptr<opals::INormals> normals(opals::INormals::New(), opals::ModuleDeleter());
  
  normals->opts().inFile().set(name + ".odm");
  normals->opts().normalsAlg().set(opals::NormalsAlgorithm::simplePlane);
  normals->opts().neighbours().set(8);
  normals->opts().searchMode().set(opals::SearchMode::d2);
  normals->opts().searchRadius().set(1);
  normals->run();
}
 
int main(int argc, char** argv) 
{
  int errorCode = 0;
 
  try {
    //ATTENTION: programm needs to be excecuted within the opals demo directory or 
    //           copy strip21.laz in your current directory
 
    // necessary preprocessing steps
    preprocessing("strip21");
 
    //perform segmentation using plane extraction mode and alpha shape computation
    std::shared_ptr<opals::ISegmentation> module( opals::ISegmentation::New(), opals::ModuleDeleter() );
    module->opts().inFile().set("strip21.odm");
    opals::Vector< float > radius(1, 1.0);
    module->opts().searchRadius().set(radius);
    module->opts().minSegSize().set(50);
    module->opts().method().set(opals::SegmentationMethod::planeExtraction);
    module->opts().planeExtraction().maxDist().set(0.2);
    module->opts().planeExtraction().maxSigma().set(0.15);
    module->opts().planeExtraction().seedCalculator().set("NormalSigma0<0.02 AND Z > 275 ? NormalSigma0 : invalid");
    opals::Vector<opals::String> filters(1);
    filters[0] = "Region[529568.8 5338773.5 529863 5338773.9 529862.6 5338642.4 529706 5338685.4 529569.1 5338753.8]";
    module->opts().filter().set(filters);
    module->opts().alphaRadius().set(1);
    module->opts().sort().set(opals::SegmentSort::descendingPSize);
    module->run();
 
    // get segmentation manager(gives access to the segment objects)
    opals::SharedPtr<opals::ISegmentManager> segManager = module->opts().segments().get();
 
    // DM::IDatamanager::load(const char *filename, bool readOnly = false, bool threadSafety = true);
    DM::DatamanagerHandle odm = DM::IDatamanager::load("strip21.odm", true, false);
    DM::PointIndexHandle  pi = odm->getPointIndex();  // get point index object of odm
 
    std::cout << "The segmentation process has found " << segManager->sizeSegments() << " segment(s)" << std::endl;
 
    unsigned max_pt_output = 5;
 
    // iterate over all segments and output relevant segments properties
    auto segs = segManager->getSegments();
    for (auto it = segs.begin(); it != segs.end(); ++it)
    {
      const opals::ISegment &seg = **it;
 
      //print basic segment information
      std::cout << std::fixed << std::setprecision(3);
      std::cout << "Segment-Nr " << seg.getID() << "\t " << seg.sizePoints() << " points; "
                << "COG (" << seg.getCoG()->x() << " / " << seg.getCoG()->y() << " / " << seg.getCoG()->z() << ")" << std::endl;
 
      // get plane parameters
      opals::Vector<double> planeParams = seg.getPlaneParams();
      assert(planeParams.size() == 5);
      std::cout << std::fixed << std::setprecision(4);
      std::cout << "\tplane parameters: a=" << planeParams[0] << " b=" << planeParams[1] << " c=" << planeParams[2]
                << " d=" << planeParams[3] << std::setprecision(3) << " (sigma0=" << planeParams[4] << ")" << std::endl;
 
      // get alpha shape(can be empty in case the computation failed)
      DM::PolylineHandle alphaShape = seg.getAlphaShape();
      if (alphaShape)
      {
        std::cout << "\talpha shape consists of " << alphaShape->sizePart() << " parts:" << std::endl;
        unsigned idx = 0;
        for (auto itPart = alphaShape->beginPart(); itPart != alphaShape->endPart(); ++itPart, ++idx)
        {
          const DM::IPolyline::IPart &part = *itPart;
          std::cout << "\t\t" << idx << " part has " << part.sizePoint() << " points";
          output_addinfo(part.info());
        }
      }
 
      // get segment point ids and print the first "max_pt_output" points
      std::cout << "\touput the first " << max_pt_output << " points: " << std::endl;
      auto ptIds = seg.getPoints();
      unsigned idx = 0;
      for (auto it = ptIds.begin(); it != ptIds.end(); ++it, ++idx)
      {
        DM::PointHandle pt = pi->getPoint(*it);
        std::cout << "\t\t" << idx << std::setprecision(3) << ".point (" << pt->x() << " / " << pt->y() << " / " << pt->z() << ")" << std::endl;
        if (idx + 1 >= max_pt_output)
          break;
      }
    }
 
  }
  catch (opals::Exceptions::Base &e)
  {
    //in this case an error occurred, output error message
    errorCode = e.errorCode();
    std::cout << e.errorMessage() << std::endl;
  }
 
  return errorCode;
}