diff --git a/cmake/CheckDGtalOptionalDependencies.cmake b/cmake/CheckDGtalOptionalDependencies.cmake index a972a38f84..eaa1ee2fd7 100644 --- a/cmake/CheckDGtalOptionalDependencies.cmake +++ b/cmake/CheckDGtalOptionalDependencies.cmake @@ -230,6 +230,8 @@ endif() set(OPENMP_FOUND_DGTAL 0) if(DGTAL_WITH_OPENMP) include(openmp) + target_link_libraries(DGtal PUBLIC OpenMP::OpenMP_CXX) + target_compile_definitions(DGtal PUBLIC -DDGTAL_WITH_OPENMP) set(DGtalLibDependencies ${DGtalLibDependencies} OpenMP::OpenMP_CXX) set(OPENMP_FOUND_DGTAL 1) set(DGTAL_WITH_OPENMP 1) diff --git a/src/DGtal/geometry/surfaces/estimation/DomainSplitter.h b/src/DGtal/geometry/surfaces/estimation/DomainSplitter.h new file mode 100644 index 0000000000..db723b0351 --- /dev/null +++ b/src/DGtal/geometry/surfaces/estimation/DomainSplitter.h @@ -0,0 +1,135 @@ +#pragma once + +#include "DGtal/kernel/domains/CDomain.h" +#include + +namespace DGtal +{ + /** + * @brief Data structure returned by Domain splitters + * + * This class offers the ability for the splitter to give information + * about the surfaces or object within the splits. + * For example, the hintVoxelCount field may allow to pre-allocate buffers. + */ + template + struct SplitInfo + { + BOOST_CONCEPT_ASSERT(( concepts::CDomain< Domain > )); + + Domain domain; //< The actual split domain + uint32_t hintVoxelCount = 0; //< An expected guess for the number of voxel + }; + + /** + * @brief Splits a domain evenly along all dimensions + * + * @tparam The model of domain to split + */ + template + struct RegularDomainSplitter + { + BOOST_CONCEPT_ASSERT(( concepts::CDomain< Domain > )); + + //Output spllitted domain type + typedef std::vector> SplitDomainsInfo; + + /** + * @brief Splits a domain + * + * This functions may split the domain in fewer part than wanted to ensure + * even subdomains. + * + * @param d The domain to split + * @param splitHint The targeted number of splits + */ + SplitDomainsInfo operator()(const Domain& d, uint32_t splitHint) const + { + // Find best match possible for even splitting + const uint32_t splitCount = std::floor(std::log(splitHint) / std::log(Domain::dimension)); + const uint32_t totalSplits = std::pow(splitCount, Domain::dimension); + + if (splitCount == 0) + return { SplitInfo{d, 0} }; + + auto splitSize = (d.upperBound() - d.lowerBound()) / (int32_t)splitCount; + SplitDomainsInfo result; + result.reserve(totalSplits); + + for (uint32_t i = 0; i < totalSplits; ++i) + { + auto start = d.lowerBound(); + auto idx = i; + for (uint32_t j = 0; j < Domain::dimension; ++j) + { + auto k = idx % splitCount; + start[j] += k * splitSize[j] + k; // +k ensure no overlap between domains + idx /= splitCount; + } + + // Make correction to ensure it remains within the domain + auto end = start + splitSize; + for (uint32_t j = 0; j < Domain::dimension; ++j) + end[j] = std::clamp(end[j], d.lowerBound()[j], d.upperBound()[j]); + + result.emplace_back(Domain(start, end), 0); + } + + return result; + }; + }; + + + /** + * @brief Splits a domain along one of the domain grid axis. + * + * @tparam The model of domain to split + */ + template + struct AxisDomainSplitter + { + + BOOST_CONCEPT_ASSERT(( concepts::CDomain< Domain > )); + + //Output spllitted domain type + typedef std::vector> SplitDomainsInfo; + + /** + * @brief Regularly splits a domain along one axis + * + * @param d The domain to split + * @param splitHint The targeted number of splits (clamped to the width of the domain) + * @param dim the split axis + */ + SplitDomainsInfo operator()(const Domain& d, uint32_t splitHint, typename Domain::Dimension dim) const + { + SplitDomainsInfo result; + if (splitHint == 0) + return result; + + auto lower = d.lowerBound(); + auto upper = d.upperBound(); + auto length = upper[dim] - lower[dim] + 1; + uint32_t splitCount = splitHint; + if (splitCount > length) + splitCount = length; + + result.reserve(splitCount); + auto base = length / splitCount; + auto rem = length % splitCount; + + auto start = lower; + for (uint32_t i = 0; i < splitCount; ++i) + { + auto size = base + (i < rem ? 1 : 0); + auto end = upper; + end[dim] = start[dim] + size - 1; + result.emplace_back(Domain(start, end), 0); + start[dim] = end[dim] + 1; + } + + return result; + }; + }; + +} diff --git a/src/DGtal/geometry/surfaces/estimation/ParallelIIEstimator.h b/src/DGtal/geometry/surfaces/estimation/ParallelIIEstimator.h new file mode 100644 index 0000000000..67469e686d --- /dev/null +++ b/src/DGtal/geometry/surfaces/estimation/ParallelIIEstimator.h @@ -0,0 +1,176 @@ +/** + * This program 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. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + **/ + +#include +#include + +#ifndef DGTAL_WITH_OPENMP +#error You need to have activated OpenMP (DGTAL_WITH_OPENMP) to include this file. +#endif + +// We require openmp for this estimator, even though it could run without +#include + +namespace DGtal +{ + /** + * @brief Run an Integral Invariant estimator in parallel + * + * This class is meant as an almost perfect replacement of + * other IIEstimator. The only difference is the constructor + * which needs the number of threads. + * + * TODO: Value to control output (ie. value + loc, just value or ordered value) + * TODO: Surface construction on subdomain fonctor ? + * + * @tparam TEstimator The model of Estimator + * @tparam TSplitter The function to split the domain + */ + template + class ParallelIIEstimator + { + public: + using Estimator = TEstimator; + using Splitter = TSplitter; + using Domain = typename TEstimator::Domain; + using Scalar = typename TEstimator::Scalar; + + using KSpace = typename TEstimator::KSpace; + using PointPredicate = typename TEstimator::PointPredicate; + using Surfel = typename KSpace::Surfel; + using SurfelSet = typename KSpace::SurfelSet; + using EstimatorQuantity = typename TEstimator::Quantity; + using Quantity = EstimatorQuantity; + + // Building the surface + using Boundary = LightImplicitDigitalSurface; + using Surface = DigitalSurface; + using Visitor = DepthFirstVisitor; + using VisitorRange = GraphVisitorRange; + + + /** + * @brief Constructor + * + * @tparam Args The estimator constructor arguments + * + * @param nbThread The number of thread to run in parallel. -1 means as many as possible + * @param args Constructor arguments to underlying estimators + */ + template + ParallelIIEstimator(int32_t nbThread, Args&&... args); + + /** + * Clears the object. It is now invalid. + */ + void clear(); + + /// @return the grid step. + Scalar h() const; + + /** + * Attach a shape, defined as a functor spel -> boolean + * + * @param[in] K the cellular grid space in which the shape is defined. + * @param aPointPredicate the shape of interest. The alias can be secured + * if a some counted pointer is handed. + */ + void attach(ConstAlias K, + ConstAlias pp); + + /** + * Set specific parameters: the radius of the ball. + * + * @param[in] dRadius the "digital" radius of the kernel (buy may be non integer). + */ + void setParams(double dRadius); + + + /** + * Checks the validity/consistency of the object. + * @return 'true' if the object is valid, 'false' otherwise. + */ + bool isValid() const; + + /** + * Model of CDigitalSurfaceLocalEstimator. Initialisation. + * + * @tparam SurfelConstIterator any model of forward readable iterator on Surfel. + * @param[in] _h grid size (must be >0). + * @param[in] ite iterator on the first surfel of the surface. + * @param[in] itb iterator after the last surfel of the surface. + */ + template + void init(double h_, ItA, ItB); + + /** + * -- Estimation -- + * + * Compute the integral invariant volume at surfel *it of + * a shape, then apply the VolumeFunctor to extract some + * geometric information. + * + * @tparam SurfelConstIterator type of Iterator on a Surfel + * + * @param[in] it iterator pointing on the surfel of the shape where + * we wish to evaluate some geometric information. + * + * @return quantity (normal vector) at surfel *it + */ + template + Quantity eval(It it); + /** + * -- Estimation -- + * + * Compute the integral invariant volume for a range of + * surfels [itb,ite) on a shape, then apply the + * VolumeFunctor to extract some geometric information. + * Return the result on an OutputIterator (param). + * + * @tparam OutputIterator type of Iterator of an array of Quantity + * @tparam SurfelConstIterator type of Iterator on a Surfel + * + * @param[in] itb iterator defining the start of the range of surfels + * where we wish to compute some geometric information. + * + * @param[in] ite iterator defining the end of the range of surfels + * where we wish to compute some geometric information. + * + * @param[in] result output iterator of results of the computation. + * @return the updated output iterator after all outputs. + */ + template + Oit eval(It itb, It ite, Oit rslt); + + /** + * Writes/Displays the object on an output stream. + * @param out the output stream where the object is written. + */ + void selfDisplay ( std::ostream & out ) const; + + private: + std::vector myEstimators; + Splitter mySplitter; + + CountedConstPtrOrConstPtr myPointPredicate; + CountedConstPtrOrConstPtr myKSpace; + + double myH; + double myRadius; + }; +} + +#include "ParallelIIEstimator.ih" diff --git a/src/DGtal/geometry/surfaces/estimation/ParallelIIEstimator.ih b/src/DGtal/geometry/surfaces/estimation/ParallelIIEstimator.ih new file mode 100644 index 0000000000..a592aadbb8 --- /dev/null +++ b/src/DGtal/geometry/surfaces/estimation/ParallelIIEstimator.ih @@ -0,0 +1,150 @@ +#include + +namespace DGtal +{ + template + template + ParallelIIEstimator:: + ParallelIIEstimator(int32_t nbThread, Args&&... args) + { + if (nbThread <= 0) + nbThread = std::max(1, omp_get_num_threads()); + + // Note: We can not use std::vector constructor or resize + // because it will copy the estimator but it can have shared + // state. + // We rather build a new one in a loop by calling the constructor + // implictely with emplace_back + myEstimators.reserve(nbThread); + for (int i = 0; i < nbThread; ++i) + myEstimators.emplace_back(std::forward(args)...); + } + + template + void ParallelIIEstimator::clear() + { + for (auto& estim : myEstimators) + estim.clear(); + } + + template + typename ParallelIIEstimator::Scalar + ParallelIIEstimator::h() const { return myEstimators[0].h(); } + + template + void ParallelIIEstimator::attach( + ConstAlias K, + ConstAlias pp) + { + clear(); + + myKSpace = K; + myPointPredicate = pp; + } + + template + void ParallelIIEstimator::setParams(double dRadius) + { + myRadius = dRadius; + } + + template + bool ParallelIIEstimator::isValid() const + { + bool valid = true; + for (const auto& estim : myEstimators) + valid = valid && estim.isValid(); + + return valid; + } + + template + template + void ParallelIIEstimator::init(double h_, ItA, ItB) + { + myH = h_; + } + + template + template + typename ParallelIIEstimator::Quantity + ParallelIIEstimator::eval(It it) + { + // TODO: Initialize properly the estimator (attaching shape, ...) + if (!myEstimators[0].isValid()) + { + auto ite = it; ite++; + myEstimators[0].init(myH, it, ite); + } + + return { + .location = *it, + .value = myEstimators[0].eval(it) + }; + } + + template + template + Oit ParallelIIEstimator::eval(It itb, It ite, Oit rslt) + { + Domain mainDomain(myKSpace->lowerBound(), myKSpace->upperBound()); + std::vector> domains = mySplitter(mainDomain, myEstimators.size()); + + // use std::map to keep insertion order so no lookup is needed afterward + std::map result; + std::vector*>> resultLocation(domains.size()); + for (auto it = itb; it != ite; it++) + { + result[*it] = 0.; // Safe, it will be overwritten (if domain really is a partition of the domain) + + // [23.2.4 Associative containers]: insersion does not invalidate pointer + for (unsigned int j = 0; j < domains.size(); ++j) + { + // Should be true for only one of the domain. The break ensure + // each surfel belongs to one and only one domain, meaning there + // will be no concurrent write afterward + if (domains[j].domain.isInside(myKSpace->interiorVoxel(*it))) + { + resultLocation[j].push_back(&(*result.find(*it))); + break; + } + } + } + + #pragma omp parallel for + for (uint32_t i = 0; i < domains.size(); ++i) + { + auto surfels = resultLocation[i] | std::views::transform([](const auto& kv) -> const Surfel& + { + return kv->first; + }); + auto values = resultLocation[i] | std::views::transform([&](auto& kv) -> EstimatorQuantity& + { + return kv->second; + }); + + auto& estim = myEstimators[i]; + + estim.attach(*myKSpace, *myPointPredicate); + estim.setParams(myRadius); + + estim.init(myH, surfels.begin(), surfels.end()); + estim.eval(surfels.begin(), surfels.end(), values.begin()); + } + + // Fuse result into the iterator + for (const auto& it : result) + *++rslt = it.second; + + return rslt; + } + + template + void + ParallelIIEstimator::selfDisplay( std::ostream & out ) const + { + out << "[ParallelIIEstimator Estim="; + myEstimators[0].selfDisplay(out); + out << "]"; + } +} diff --git a/tests/geometry/surfaces/CMakeLists.txt b/tests/geometry/surfaces/CMakeLists.txt index 50d555053f..786bfca07a 100644 --- a/tests/geometry/surfaces/CMakeLists.txt +++ b/tests/geometry/surfaces/CMakeLists.txt @@ -1,16 +1,9 @@ -set(TESTS_SRC - testArithmeticalDSSComputerOnSurfels +set(TESTS_SURFACES_SRC + testArithmeticalDSSComputerOnSurfels testChordGenericStandardPlaneComputer testDigitalPlanePredicate testPlaneProbingTetrahedronEstimator testPlaneProbingParallelepipedEstimator - ) - -foreach(FILE ${TESTS_SRC}) - DGtal_add_test(${FILE}) -endforeach() - -set(TESTS_SURFACES_SRC testIntegralInvariantShortcuts testNormalVectorEstimatorEmbedder testIntegralInvariantVolumeEstimator @@ -22,14 +15,19 @@ set(TESTS_SURFACES_SRC testSphericalHoughNormalVectorEstimator testDigitalSurfaceRegularization testShroudsRegularization + testDomainSplitter ) foreach(FILE ${TESTS_SURFACES_SRC}) DGtal_add_test(${FILE}) endforeach() +if ( DGTAL_WITH_OPENMP ) + DGtal_add_test(testParallelIntegralInvariantEstimator) +endif() + -if ( DGTAL_WITH_CGAL ) +if ( DGTAL_WITH_CGAL ) set(CGAL_TESTS_SRC testMonge ) foreach(FILE ${CGAL_TESTS_SRC}) @@ -38,7 +36,7 @@ if ( DGTAL_WITH_CGAL ) endif() -if ( DGTAL_WITH_POLYSCOPE_VIEWER ) +if ( DGTAL_WITH_POLYSCOPE_VIEWER ) set(POLYSCOPE_VIEWER_TESTS_SRC testLocalConvolutionNormalVectorEstimator testTensorVotingViewer) @@ -48,8 +46,6 @@ if ( DGTAL_WITH_POLYSCOPE_VIEWER ) endforeach() endif() - - if ( DGTAL_WITH_PONCA ) set(PONCA_TESTS_SRC testSphereFitting ) diff --git a/tests/geometry/surfaces/testDomainSplitter.cpp b/tests/geometry/surfaces/testDomainSplitter.cpp new file mode 100644 index 0000000000..83e9a5022b --- /dev/null +++ b/tests/geometry/surfaces/testDomainSplitter.cpp @@ -0,0 +1,109 @@ +/** + * This program 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. + * + * This program 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 General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program. If not, see . + * + **/ + +/** + * @file testDomainSplitter.cpp + * @ingroup Tests + * @author David Coeurjolly (\c david.coeurjolly@cnrs.fr) + * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), INSA-Lyon, France + * + * @date 2026/04/03 + * + * Functions for testing class DomainSplitter + * + * This file is part of the DGtal library. + */ + +/////////////////////////////////////////////////////////////////////////////// +#include +#include "DGtal/base/Common.h" +#include "DGtalCatch.h" + +#include "DGtal/geometry/surfaces/estimation/DomainSplitter.h" + +#ifdef DGTAL_WITH_POLYSCOPE_VIEWER +#include "DGtal/io/colormaps/HueShadeColorMap.h" +#include "DGtal/io/viewers/PolyscopeViewer.h" +#endif + +/////////////////////////////////////////////////////////////////////////////// + +using namespace DGtal; +using namespace Z3i; + +TEST_CASE( "Domain Regular Grid Splitter tests" ) +{ + Domain domain(Point(0,0,0), Point(16,32,64)); + + RegularDomainSplitter splitter; + + RegularDomainSplitter::SplitDomainsInfo output = splitter(domain,12); + REQUIRE( output.size() <= 12); + + trace.info() << "Original domain: "< cmap(0,(unsigned int)output.size()); + for(auto i=0; i< output.size(); ++i) + { + viewer << cmap(i); + viewer << output[i].domain; + } + viewer.show(); +#endif +} + +TEST_CASE( "Domain Axis Splitter tests" ) +{ + Domain domain(Point(0,0,0), Point(16,32,64)); + + AxisDomainSplitter splitter; + + AxisDomainSplitter::SplitDomainsInfo output = splitter(domain,3,0); + REQUIRE( output.size() == 3); + + trace.info() << "Original domain: "< cmap(0,(unsigned int)output.size()); + for(auto i=0; i< output.size(); ++i) + { + viewer << cmap(i); + viewer << output[i].domain; + } + viewer.show(); +#endif +} + +TEST_CASE( "Domain Axis Splitter tests (another direction)" ) +{ + Domain domain(Point(10,10,10), Point(16,32,64)); + + AxisDomainSplitter splitter; + + AxisDomainSplitter::SplitDomainsInfo output = splitter(domain,2,1); + REQUIRE( output.size() == 2); + + trace.info() << "Original domain: "<. + * + **/ + +/** + * @file testIntegralInvariantVolumeEstimator.cpp + * @ingroup Tests + * @author Bastien DOIGNIES (\c bastien.doignies@liris.cnrs.fr ) + * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), INSA-Lyon, France + * LAboratoire de MAthématiques - LAMA (CNRS, UMR 5127), Université de Savoie, France + * + * @date 2014/06/26 + * + * Functions for testing class IntegralInvariantVolumeEstimator and IIGeometricFunctor. + * + * This file is part of the DGtal library. + */ + +/////////////////////////////////////////////////////////////////////////////// +#include +#include "DGtal/base/Common.h" + +/// Shape +#include "DGtal/shapes/implicit/ImplicitBall.h" + +/// Digitization +#include "DGtal/shapes/GaussDigitizer.h" +#include "DGtal/topology/LightImplicitDigitalSurface.h" +#include "DGtal/topology/DigitalSurface.h" +#include "DGtal/graph/DepthFirstVisitor.h" +#include "DGtal/graph/GraphVisitorRange.h" + +/// Estimator +#include "DGtal/geometry/surfaces/estimation/IIGeometricFunctors.h" +#include "DGtal/geometry/surfaces/estimation/IntegralInvariantVolumeEstimator.h" +#include "DGtal/geometry/surfaces/estimation/ParallelIIEstimator.h" +#include "DGtal/geometry/surfaces/estimation/DomainSplitter.h" + + +/////////////////////////////////////////////////////////////////////////////// + + +using namespace DGtal; + +/////////////////////////////////////////////////////////////////////////////// +// Functions for testing class IntegralInvariantVolumeEstimator and IIGeometricFunctor. +/////////////////////////////////////////////////////////////////////////////// + +bool testCurvature2dP ( double h, double delta ) +{ + typedef ImplicitBall ImplicitShape; + typedef GaussDigitizer DigitalShape; + typedef LightImplicitDigitalSurface Boundary; + typedef DigitalSurface< Boundary > MyDigitalSurface; + typedef DepthFirstVisitor< MyDigitalSurface > Visitor; + typedef GraphVisitorRange< Visitor > VisitorRange; + typedef VisitorRange::ConstIterator VisitorConstIterator; + + typedef functors::IICurvatureFunctor MyIICurvatureFunctor; + typedef IntegralInvariantVolumeEstimator< Z2i::KSpace, DigitalShape, MyIICurvatureFunctor > MyIICurvatureEstimator; + typedef RegularDomainSplitter> Splitter; + typedef ParallelIIEstimator MyIICurvatureEstimatorA; + // typedef MyIICurvatureFunctor::Value Value; + typedef MyIICurvatureEstimatorA::Quantity Value; + + double re = 10; + double radius = 15; + double realValue = 1.0/radius; + + trace.beginBlock( "[PARALLEL] Shape initialisation ..." ); + + ImplicitShape ishape( Z2i::RealPoint( 0, 0 ), radius ); + DigitalShape dshape; + dshape.attach( ishape ); + dshape.init( Z2i::RealPoint( -20.0, -20.0 ), Z2i::RealPoint( 20.0, 20.0 ), h ); + + Z2i::KSpace K; + if ( !K.init( dshape.getLowerBound(), dshape.getUpperBound(), true ) ) + { + trace.error() << "Problem with Khalimsky space" << std::endl; + return false; + } + + Z2i::KSpace::Surfel bel = Surfaces::findABel( K, dshape, 10000 ); + Boundary boundary( K, dshape, SurfelAdjacency( true ), bel ); + MyDigitalSurface surf ( boundary ); + + trace.endBlock(); + + trace.beginBlock( "Curvature estimator initialisation ..."); + + VisitorRange range( new Visitor( surf, *surf.begin() )); + VisitorConstIterator ibegin = range.begin(); + VisitorConstIterator iend = range.end(); + + MyIICurvatureFunctor curvatureFunctor; + curvatureFunctor.init( h, re ); + + MyIICurvatureEstimatorA curvatureEstimator( 4, curvatureFunctor ); + curvatureEstimator.attach( K, dshape ); + curvatureEstimator.setParams( re/h ); + curvatureEstimator.init( h, ibegin, iend ); + + trace.endBlock(); + + trace.beginBlock( "Curvature estimator evaluation ..."); + + std::vector< Value > results; + std::back_insert_iterator< std::vector< Value > > resultsIt( results ); + curvatureEstimator.eval( ibegin, iend, resultsIt ); + + trace.endBlock(); + + trace.beginBlock ( "Comparing results of integral invariant 2D curvature ..." ); + + double mean = 0.0; + double rsize = static_cast(results.size()); + + if( rsize == 0 ) + { + trace.error() << "ERROR: surface is empty" << std::endl; + trace.endBlock(); + return false; + } + + for ( unsigned int i = 0; i < rsize; ++i ) + { + mean += results[ i ]; + } + mean /= rsize; + + if( mean != mean ) //NaN + { + trace.error() << "ERROR: result is NaN" << std::endl; + trace.endBlock(); + return false; + } + + double v = std::abs ( realValue - mean ); + + trace.warning() << "True value: " << realValue << std::endl; + trace.warning() << "Mean value: " << mean << std::endl; + trace.warning() << "Count value: " << rsize << std::endl; + trace.warning() << "Delta: " << delta << " |true - mean|: " << v << std::endl; + + if( v > delta ) + { + trace.endBlock(); + return false; + } + trace.endBlock(); + return true; +} + +bool testCurvature2d ( double h, double delta ) +{ + typedef ImplicitBall ImplicitShape; + typedef GaussDigitizer DigitalShape; + typedef LightImplicitDigitalSurface Boundary; + typedef DigitalSurface< Boundary > MyDigitalSurface; + typedef DepthFirstVisitor< MyDigitalSurface > Visitor; + typedef GraphVisitorRange< Visitor > VisitorRange; + typedef VisitorRange::ConstIterator VisitorConstIterator; + + typedef functors::IICurvatureFunctor MyIICurvatureFunctor; + typedef IntegralInvariantVolumeEstimator< Z2i::KSpace, DigitalShape, MyIICurvatureFunctor > MyIICurvatureEstimator; + typedef MyIICurvatureFunctor::Value Value; + + double re = 10; + double radius = 15; + double realValue = 1.0/radius; + + trace.beginBlock( "Shape initialisation ..." ); + + ImplicitShape ishape( Z2i::RealPoint( 0, 0 ), radius ); + DigitalShape dshape; + dshape.attach( ishape ); + dshape.init( Z2i::RealPoint( -20.0, -20.0 ), Z2i::RealPoint( 20.0, 20.0 ), h ); + + Z2i::KSpace K; + + if ( !K.init( dshape.getLowerBound(), dshape.getUpperBound(), true ) ) + { + trace.error() << "Problem with Khalimsky space" << std::endl; + return false; + } + + Z2i::KSpace::Surfel bel = Surfaces::findABel( K, dshape, 10000 ); + Boundary boundary( K, dshape, SurfelAdjacency( true ), bel ); + MyDigitalSurface surf ( boundary ); + + trace.endBlock(); + + trace.beginBlock( "Curvature estimator initialisation ..."); + + VisitorRange range( new Visitor( surf, *surf.begin() )); + VisitorConstIterator ibegin = range.begin(); + VisitorConstIterator iend = range.end(); + + MyIICurvatureFunctor curvatureFunctor; + curvatureFunctor.init( h, re ); + + MyIICurvatureEstimator curvatureEstimator( curvatureFunctor ); + curvatureEstimator.attach( K, dshape ); + curvatureEstimator.setParams( re/h ); + curvatureEstimator.init( h, ibegin, iend ); + std::cout << curvatureEstimator << std::endl; + + trace.endBlock(); + + trace.beginBlock( "Curvature estimator evaluation ..."); + + std::vector< Value > results; + std::back_insert_iterator< std::vector< Value > > resultsIt( results ); + curvatureEstimator.eval( ibegin, iend, resultsIt ); + + trace.endBlock(); + + trace.beginBlock ( "Comparing results of integral invariant 2D curvature ..." ); + + double mean = 0.0; + double rsize = static_cast(results.size()); + + if( rsize == 0 ) + { + trace.error() << "ERROR: surface is empty" << std::endl; + trace.endBlock(); + return false; + } + + for ( unsigned int i = 0; i < rsize; ++i ) + { + mean += results[ i ]; + } + mean /= rsize; + + if( mean != mean ) //NaN + { + trace.error() << "ERROR: result is NaN" << std::endl; + trace.endBlock(); + return false; + } + + double v = std::abs ( realValue - mean ); + + trace.warning() << "True value: " << realValue << std::endl; + trace.warning() << "Mean value: " << mean << std::endl; + trace.warning() << "Count value: " << rsize << std::endl; + trace.warning() << "Delta: " << delta << " |true - mean|: " << v << std::endl; + if( v > delta ) + { + trace.endBlock(); + return false; + } + trace.endBlock(); + return true; +} + + +/////////////////////////////////////////////////////////////////////////////// +// Standard services - public : +int main( int /*argc*/, char** /*argv*/ ) +{ + trace.beginBlock ( "Testing class ParrallelIIEstimator with IntegralInvariantVolumeEstimator in 2d" ); + + bool res = testCurvature2d( 0.05, 0.002 ) && testCurvature2dP(0.05, 0.002) ; + trace.emphase() << ( res ? "Passed." : "Error." ) << std::endl; + trace.endBlock(); + return res ? 0 : 1; +}