exampleMaximalSegmentSliceEstimation.cpp

Go to the documentation of this file.

#include <iostream>
#include "ConfigExamples.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/base/Common.h"
#include "DGtal/geometry/surfaces/estimation/MaximalSegmentSliceEstimation.h"
#include "DGtal/helpers/Shortcuts.h"
#include "DGtal/io/viewers/PolyscopeViewer.h"
 
using namespace std;
using namespace DGtal;
 
using KSpace    = Z3i::KSpace;
using SH3       = Shortcuts<KSpace>;
using Surface   = SH3::DigitalSurface;
using Surfel    = Surface::Surfel;
using RealPoint = SH3::RealPoint;
using Point     = SH3::Point;
 
RealPoint centerSurfel (KSpace const& K, SH3::SCell const& s)
{
    auto pointels = SH3::getPrimalVertices(K, s, true);
 
    Point p = K.uCoords(pointels[0]),
          u = K.uCoords(pointels[1]) - p,
          v = K.uCoords(pointels[3]) - p;
 
    static const RealPoint shift(-0.5, -0.5, -0.5);
 
    return p + 0.5 * u + 0.5 * v + shift;
}

 
int main()
{
    auto params = SH3::defaultParameters();
    params("polynomial", "ellipsoid")("gridstep", "1.0");
    auto implicit_shape  = SH3::makeImplicitShape3D  ( params );
    auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
    auto K               = SH3::getKSpace( params );
    auto binary_image    = SH3::makeBinaryImage( digitized_shape, params );
    auto surface         = SH3::makeDigitalSurface( binary_image, K, params );
    auto surfels         = SH3::getSurfelRange(surface);
 
    double gridstep = params["gridstep"].as<double>();
 
    PolyscopeViewer<> viewer(K);

    // Instantiation
    using Estimator = MaximalSegmentSliceEstimation<Surface>;
    Estimator estimator;
    estimator.init(gridstep, surfels.begin(), surfels.end());
    estimator.attach(surface);
 
    // Usage
    std::vector<Estimator::Quantity> quantities;
    estimator.eval(surfels.begin(), surfels.end(), std::back_inserter(quantities));
 
    Color fillColor = Color::White;
 
    for (std::size_t i = 0; i < surfels.size(); ++i)
    {
        const Surfel& s = surfels[i];
        const Estimator::Quantity& normal = quantities[i];
 
        const RealPoint& n = normal.getNormalized();
        RealPoint origin = centerSurfel(K, s);
 
        viewer.drawColor(fillColor);
        viewer << s;
 
        viewer.drawColor(Color::Green);
        viewer.drawLine(origin, origin + 1.5 * n);
    }
 
    viewer.show();
    return 0;
}
//                                                                           //