testIntegralInvariantShortcuts.cpp

Go to the documentation of this file.

 
#include <iostream>
#include <vector>
 
#include "DGtal/base/Common.h"
#include "ConfigTest.h"
#include "DGtalCatch.h"
#include "DGtal/helpers/StdDefs.h"
#include "DGtal/helpers/Shortcuts.h"
#include "DGtal/helpers/ShortcutsGeometry.h"
 
 
using namespace std;
using namespace DGtal;
 
typedef Shortcuts<Z3i::KSpace>         SH3;
typedef ShortcutsGeometry<Z3i::KSpace> SHG3;
 
// Functions for testing class IntegralInvariantShortcuts.
 
TEST_CASE( "Testing IntegralInvariant Shortcuts API" )
{
  auto params = SH3::defaultParameters() | SHG3::defaultParameters() |  SHG3::parametersGeometryEstimation();
  params( "polynomial", "goursat" )( "gridstep", 1. );
  auto implicit_shape  = SH3::makeImplicitShape3D  ( params );
  auto digitized_shape = SH3::makeDigitizedImplicitShape3D( implicit_shape, params );
  auto binary_image    = SH3::makeBinaryImage( digitized_shape, params );
  auto K               = SH3::getKSpace( params );
  auto embedder        = SH3::getCellEmbedder( K );
  auto surface         = SH3::makeLightDigitalSurface( binary_image, K, params );
  auto surfels         = SH3::getSurfelRange( surface, params );
 
  trace.info() << "Nb surfels= " << surfels.size() << std::endl;
 
  //Computing some differential quantities
  params("r-radius", 3.0);
 
  //We compute the curvature tensor, the mean and the Gaussian curvature
  auto Hcurv     = SHG3::getIIMeanCurvatures( binary_image, surfels, params );
  auto Tcurv     = SHG3::getIIPrincipalCurvaturesAndDirections(binary_image, surfels, params);
  auto Kcurv     = SHG3::getIIGaussianCurvatures( binary_image, surfels, params);
  auto Ncurv     = SHG3::getIINormalVectors( binary_image, surfels, params );
 
  auto params_parallel = params;
  params_parallel( "ii-thread-number", 4 )( "ii-split-axis", 2 );
  auto HcurvParallel = SHG3::getIIMeanCurvatures( binary_image, surfels, params_parallel );
  auto TcurvParallel = SHG3::getIIPrincipalCurvaturesAndDirections( binary_image, surfels, params_parallel );
  auto KcurvParallel = SHG3::getIIGaussianCurvatures( binary_image, surfels, params_parallel );
  auto NcurvParallel = SHG3::getIINormalVectors( binary_image, surfels, params_parallel );
 
  std::vector<double> k1,k2,G;
  for(auto &result: Tcurv)
  {
    k1.push_back( std::get<0>(result) );
    k2.push_back( std::get<1>(result) );
    G.push_back( ( std::get<0>(result) * std::get<1>(result)) );
  }
 
  SECTION("Testing that mean/Gaussian/tensor curvature shortucut values match")
  {
    for(std::size_t i = 0; i < G.size(); ++i)
      REQUIRE( Kcurv[i] == Approx( G[i] ) );
  }
 
  SECTION("Testing that requesting the parallel II shortcut preserves curvature values")
  {
    REQUIRE( HcurvParallel.size() == Hcurv.size() );
    REQUIRE( KcurvParallel.size() == Kcurv.size() );
    REQUIRE( TcurvParallel.size() == Tcurv.size() );
 
    for ( std::size_t i = 0; i < Hcurv.size(); ++i )
      REQUIRE( HcurvParallel[ i ] == Approx( Hcurv[ i ] ) );
 
    for ( std::size_t i = 0; i < Kcurv.size(); ++i )
      REQUIRE( KcurvParallel[ i ] == Approx( Kcurv[ i ] ) );
 
    for ( std::size_t i = 0; i < Tcurv.size(); ++i )
    {
      REQUIRE( std::get<0>( TcurvParallel[ i ] ) == Approx( std::get<0>( Tcurv[ i ] ) ) );
      REQUIRE( std::get<1>( TcurvParallel[ i ] ) == Approx( std::get<1>( Tcurv[ i ] ) ) );
    }
  }
 
  SECTION("Testing that requesting the parallel II shortcut preserves normal vectors")
  {
    REQUIRE( NcurvParallel.size() == Ncurv.size() );
 
    for ( std::size_t i = 0; i < Ncurv.size(); ++i )
      for ( std::size_t d = 0; d < 3; ++d )
        REQUIRE( NcurvParallel[ i ][ d ] == Approx( Ncurv[ i ][ d ] ) );
  }
 
  SECTION("Testing that ii-split-axis accepts out-of-range values by clamping to a valid axis")
  {
    auto params_parallel_clamped = params;
    params_parallel_clamped( "ii-thread-number", 4 )( "ii-split-axis", 9 );
    auto HcurvParallelClamped = SHG3::getIIMeanCurvatures( binary_image, surfels, params_parallel_clamped );
 
    REQUIRE( HcurvParallelClamped.size() == Hcurv.size() );
    for ( std::size_t i = 0; i < Hcurv.size(); ++i )
      REQUIRE( HcurvParallelClamped[ i ] == Approx( Hcurv[ i ] ) );
  }
 
  SECTION("Testing on shifted domains")
  {
    auto SHIFT=512;
    auto domain = binary_image->domain();
    Z3i::Domain shifted(domain.lowerBound() + Z3i::Point::diagonal(SHIFT), domain.upperBound() + Z3i::Point::diagonal(SHIFT) );
 
    SH3::KSpace Ks;
    Ks.init( shifted.lowerBound(), shifted.upperBound(), true );
 
    CountedPtr<SH3::BinaryImage>  binary_image_shifted(new SH3::BinaryImage(shifted));
    for(auto p : binary_image->domain())
      binary_image_shifted->setValue(p+Z3i::Point::diagonal(SHIFT), binary_image->operator()(p));
 
    auto surfaceShifted  = SH3::makeLightDigitalSurface( binary_image_shifted, Ks, params );
    auto surfelsShifted  = SH3::getSurfelRange( surfaceShifted, params );
    trace.info() << "Nb surfels= " << surfels.size() << " "<<shifted<<" "<<Ks<<std::endl;
 
    //Computing some differential quantities
    auto KcurvShifted     = SHG3::getIIGaussianCurvatures( binary_image_shifted, surfelsShifted, params);
 
 
  }
}