LambdaMST3D.ih

/**
 *  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 <http://www.gnu.org/licenses/>.
 *
 **/
 
/**
 * @file LambdaMST3D.ih
 * @author Kacper Pluta (\c kacper.pluta@esiee.fr )
 * Laboratoire d'Informatique Gaspard-Monge - LIGM, France
 *
 * @date 2014/10/06
 *
 * This file is part of the DGtal library.
 */
 
namespace DGtal
{
 
  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
  inline
  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::LambdaMST3DEstimator() : myBegin(), myEnd(), dssSegments ( 0 ), myFunctor() {}
 
 
  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
  inline
  void
  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::init ( ConstIterator itb, ConstIterator ite )
  {
    myBegin = itb;
    myEnd = ite;
  }
 
  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
  inline
  void
  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::attach ( Alias<TSegmentation> segmentComputer )
  {
    dssSegments = &segmentComputer;
  }
 
  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
  inline
  bool
  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::isValid () const
  {
    return ( dssSegments != 0 );
  }
 
 
  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
  inline
  DSSFilter &
  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::getDSSFilter ( )
  {
    return myDSSFilter;
  }
 
  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
  inline
  typename LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::Value
  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::treatOrphan ( OrphanDSSIterator begin, OrphanDSSIterator end,
                                                                                     const Point & p )
  {
    Value tangent, prev, partial;
    for ( auto it = begin ; it != end; ++it )
    {
      // the returned type is signed but dssLen should never be negative
      const auto dssLen = std::distance ( it->begin(), it->end() ) + 1;
      prev = partial;
      const auto pos = myDSSFilter. position ( *it, p );
      partial = myFunctor ( *it, pos, dssLen );
      if ( prev.first.cosineSimilarity ( partial.first ) > M_PI_2 )
        partial.first = -partial.first;
      tangent += partial;
    }
    return tangent;
  }
 
 
  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
  template < typename DSSesIterator, typename OrphanIterator >
  inline
  void
  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::treatOrphans ( DSSesIterator begin,
                                                                                    DSSesIterator end,
                                                                                    OrphanIterator obegin,
                                                                                    OrphanIterator oend,
                                                                                    std::multimap < Point, Value > & outValues )
  {
    for ( auto DSS = begin; DSS != end; ++DSS )
    {
      for ( auto orphan = obegin; orphan != oend; ++orphan )
      {
        if ( ! DSS->isInDSS ( *orphan ) && myDSSFilter.admissibility ( *DSS, *orphan ) )
        {
          // the returned type is signed but dssLen should never be negative
          const auto dssLen = std::distance ( DSS->begin ( ), DSS->end ( ) ) + 1;
          const auto pos = myDSSFilter. position ( *DSS, *orphan );
          outValues.insert ( std::make_pair ( *orphan,  myFunctor ( *DSS, pos, dssLen ) ) );
        }
      }
    }
  }
 
 
  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
  inline
  typename LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::RealVector
  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::eval( const Point & p )
  {
    assert ( isValid() );
    std::vector < SegmentComputer > admissibleDSS;
    bool isOrphan = true;
    typename TSegmentation::SegmentComputerIterator DSS = dssSegments->begin();
    typename TSegmentation::SegmentComputerIterator lastDSS = dssSegments->end();
    Value tangent, partial, prev;
  
    for ( ; DSS != lastDSS; ++DSS )
    {
      if ( myDSSFilter ( *DSS ) )
        continue;
 
      // collect DSSes that may need to be required to ensure coverage if the only covering DSS was filtered out
      if ( isOrphan && ! DSS->isInDSS ( p ) && myDSSFilter.admissibility ( *DSS, p ) )
        admissibleDSS.push_back ( *DSS );
 
      if ( DSS->isInDSS ( p ) )
      {
        // the returned type is signed but dssLen should never be negative
        const auto dssLen = std::distance ( DSS.begin(), DSS.end() ) + 1;
        prev = partial;
        // the returned type is signed but pos should never be negative
        const auto  pos = std::distance ( DSS.begin(), std::find ( DSS.begin ( ), DSS.end ( ), p ) ) + 1;
        partial = myFunctor ( *DSS, pos, dssLen );
        if ( partial.first.norm() > 0. && prev.first.norm() > 0. && prev.first.cosineSimilarity ( partial.first ) > M_PI_2 )
              partial.first = -partial.first;
        tangent += partial;
        isOrphan = false;
      }
    }
 
    if ( isOrphan && admissibleDSS.size ( ) > 0 )
      tangent = treatOrphan ( admissibleDSS.cbegin ( ), admissibleDSS.cend ( ), p );
    else if ( isOrphan && admissibleDSS.size ( ) == 0 )
      throw std::runtime_error ( "The DSS filter is not well constructed! The point is not DSS covered!" );
 
    if ( tangent.second != 0. )
      return tangent.first / tangent.second;
    else
      return tangent.first;
  }
 
  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
  template < typename OutputIterator >
  inline
  OutputIterator
  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::eval ( ConstIterator itb, ConstIterator ite,
                                                                            OutputIterator result )
  {
    assert ( myBegin != myEnd && isValid() && myBegin <= itb && ite <= myEnd && itb != ite );
    std::multimap < Point, Value > outValues;
    dssSegments->setSubRange ( itb, ite );
    typename TSegmentation::SegmentComputerIterator DSS = dssSegments->begin();
    typename TSegmentation::SegmentComputerIterator lastDSS = dssSegments->end();
    std::vector< Point > orphans;
 
    for(; DSS != lastDSS; ++DSS)
    {
      // coolect potential orphans
      if ( myDSSFilter ( *DSS ) )
      {
        for ( const auto & point : *DSS )
          if ( std::find ( orphans.cbegin (  ), orphans.cend ( ), point ) == orphans.cend ( ) )
            orphans.push_back ( point );
        continue;
      }
 
      auto dssLen = std::distance ( DSS.begin(), DSS.end() );
      for ( unsigned int indexOfPointInDSS = 0; indexOfPointInDSS < dssLen; indexOfPointInDSS++ )
      {
        outValues.insert ( std::make_pair ( *(DSS.begin ( ) + indexOfPointInDSS),
                           myFunctor ( *DSS, indexOfPointInDSS + 1, dssLen + 1 ) ) );
        // if a point is covered then it is not an orphan
        auto orphan = std::find ( orphans.begin ( ), orphans.end ( ), *(DSS.begin ( ) + indexOfPointInDSS) );
        if ( orphan != orphans.end ( ) )
          orphans.erase ( orphan );
      }
    }
    if ( ! orphans.empty ( ) )
      treatOrphans ( dssSegments->begin ( ),  dssSegments->end ( ), orphans.cbegin ( ), orphans.cend ( ), outValues );
    accumulate< OutputIterator >( outValues, itb, ite, result );
    return result;
  }
 
  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
  template <typename OutputIterator>
  inline
  void
  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::accumulate ( std::multimap < Point, Value > & outValues,
                                                                       ConstIterator itb, ConstIterator ite,
                                                                       OutputIterator & result )
  {
    Value prev = outValues.lower_bound( *itb )->second;
    Value accum_prev = outValues.lower_bound( *itb )->second;
    for ( ConstIterator itt = itb; itt != ite; ++itt )
    {
      typename std::multimap< Point, Value >::const_iterator it  = outValues.lower_bound ( *itt );
      typename std::multimap< Point, Value >::const_iterator it2 = outValues.upper_bound ( *itt );
      Value tangent;
      for (; it != it2; it++ )
      {
        Value partial = it->second;
        if ( partial.first.norm() > 0. && prev.first.norm() > 0. && prev.first.cosineSimilarity ( partial.first ) > M_PI_2 )
              partial.first = -partial.first;
        prev = partial;
        tangent += partial;
      }
      outValues.erase ( *itt );
      // avoid tangent flapping
      if ( accum_prev.first.norm() > 0. && tangent.first.norm() > 0. && accum_prev.first.cosineSimilarity ( tangent.first ) > M_PI_2 )
        tangent.first = -tangent.first;
      accum_prev = tangent;
      if ( tangent.second != 0 )
        *result++ = ( tangent.first / tangent.second );
      else
        *result++ = tangent.first;
    }
  }
}