DGtal  1.5.beta
LambdaMST3D.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file LambdaMST3D.ih
19  * @author Kacper Pluta (\c kacper.pluta@esiee.fr )
20  * Laboratoire d'Informatique Gaspard-Monge - LIGM, France
21  *
22  * @date 2014/10/06
23  *
24  * This file is part of the DGtal library.
25  */
26 
27 namespace DGtal
28 {
29 
30  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
31  inline
32  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::LambdaMST3DEstimator() : myBegin(), myEnd(), dssSegments ( 0 ), myFunctor() {}
33 
34 
35  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
36  inline
37  void
38  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::init ( ConstIterator itb, ConstIterator ite )
39  {
40  myBegin = itb;
41  myEnd = ite;
42  }
43 
44  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
45  inline
46  void
47  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::attach ( Alias<TSegmentation> segmentComputer )
48  {
49  dssSegments = &segmentComputer;
50  }
51 
52  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
53  inline
54  bool
55  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::isValid () const
56  {
57  return ( dssSegments != 0 );
58  }
59 
60 
61  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
62  inline
63  DSSFilter &
64  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::getDSSFilter ( )
65  {
66  return myDSSFilter;
67  }
68 
69  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
70  inline
71  typename LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::Value
72  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::treatOrphan ( OrphanDSSIterator begin, OrphanDSSIterator end,
73  const Point & p )
74  {
75  Value tangent, prev, partial;
76  for ( auto it = begin ; it != end; ++it )
77  {
78  // the returned type is signed but dssLen should never be negative
79  const auto dssLen = std::distance ( it->begin(), it->end() ) + 1;
80  prev = partial;
81  const auto pos = myDSSFilter. position ( *it, p );
82  partial = myFunctor ( *it, pos, dssLen );
83  if ( prev.first.cosineSimilarity ( partial.first ) > M_PI_2 )
84  partial.first = -partial.first;
85  tangent += partial;
86  }
87  return tangent;
88  }
89 
90 
91  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
92  template < typename DSSesIterator, typename OrphanIterator >
93  inline
94  void
95  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::treatOrphans ( DSSesIterator begin,
96  DSSesIterator end,
97  OrphanIterator obegin,
98  OrphanIterator oend,
99  std::multimap < Point, Value > & outValues )
100  {
101  for ( auto DSS = begin; DSS != end; ++DSS )
102  {
103  for ( auto orphan = obegin; orphan != oend; ++orphan )
104  {
105  if ( ! DSS->isInDSS ( *orphan ) && myDSSFilter.admissibility ( *DSS, *orphan ) )
106  {
107  // the returned type is signed but dssLen should never be negative
108  const auto dssLen = std::distance ( DSS->begin ( ), DSS->end ( ) ) + 1;
109  const auto pos = myDSSFilter. position ( *DSS, *orphan );
110  outValues.insert ( std::make_pair ( *orphan, myFunctor ( *DSS, pos, dssLen ) ) );
111  }
112  }
113  }
114  }
115 
116 
117  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
118  inline
119  typename LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::RealVector
120  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::eval( const Point & p )
121  {
122  assert ( isValid() );
123  std::vector < SegmentComputer > admissibleDSS;
124  bool isOrphan = true;
125  typename TSegmentation::SegmentComputerIterator DSS = dssSegments->begin();
126  typename TSegmentation::SegmentComputerIterator lastDSS = dssSegments->end();
127  Value tangent, partial, prev;
128 
129  for ( ; DSS != lastDSS; ++DSS )
130  {
131  if ( myDSSFilter ( *DSS ) )
132  continue;
133 
134  // collect DSSes that may need to be required to ensure coverage if the only covering DSS was filtered out
135  if ( isOrphan && ! DSS->isInDSS ( p ) && myDSSFilter.admissibility ( *DSS, p ) )
136  admissibleDSS.push_back ( *DSS );
137 
138  if ( DSS->isInDSS ( p ) )
139  {
140  // the returned type is signed but dssLen should never be negative
141  const auto dssLen = std::distance ( DSS.begin(), DSS.end() ) + 1;
142  prev = partial;
143  // the returned type is signed but pos should never be negative
144  const auto pos = std::distance ( DSS.begin(), std::find ( DSS.begin ( ), DSS.end ( ), p ) ) + 1;
145  partial = myFunctor ( *DSS, pos, dssLen );
146  if ( partial.first.norm() > 0. && prev.first.norm() > 0. && prev.first.cosineSimilarity ( partial.first ) > M_PI_2 )
147  partial.first = -partial.first;
148  tangent += partial;
149  isOrphan = false;
150  }
151  }
152 
153  if ( isOrphan && admissibleDSS.size ( ) > 0 )
154  tangent = treatOrphan ( admissibleDSS.cbegin ( ), admissibleDSS.cend ( ), p );
155  else if ( isOrphan && admissibleDSS.size ( ) == 0 )
156  throw std::runtime_error ( "The DSS filter is not well constructed! The point is not DSS covered!" );
157 
158  if ( tangent.second != 0. )
159  return tangent.first / tangent.second;
160  else
161  return tangent.first;
162  }
163 
164  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
165  template < typename OutputIterator >
166  inline
167  OutputIterator
168  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::eval ( ConstIterator itb, ConstIterator ite,
169  OutputIterator result )
170  {
171  assert ( myBegin != myEnd && isValid() && myBegin <= itb && ite <= myEnd && itb != ite );
172  std::multimap < Point, Value > outValues;
173  dssSegments->setSubRange ( itb, ite );
174  typename TSegmentation::SegmentComputerIterator DSS = dssSegments->begin();
175  typename TSegmentation::SegmentComputerIterator lastDSS = dssSegments->end();
176  std::vector< Point > orphans;
177 
178  for(; DSS != lastDSS; ++DSS)
179  {
180  // coolect potential orphans
181  if ( myDSSFilter ( *DSS ) )
182  {
183  for ( const auto & point : *DSS )
184  if ( std::find ( orphans.cbegin ( ), orphans.cend ( ), point ) == orphans.cend ( ) )
185  orphans.push_back ( point );
186  continue;
187  }
188 
189  auto dssLen = std::distance ( DSS.begin(), DSS.end() );
190  for ( unsigned int indexOfPointInDSS = 0; indexOfPointInDSS < dssLen; indexOfPointInDSS++ )
191  {
192  outValues.insert ( std::make_pair ( *(DSS.begin ( ) + indexOfPointInDSS),
193  myFunctor ( *DSS, indexOfPointInDSS + 1, dssLen + 1 ) ) );
194  // if a point is covered then it is not an orphan
195  auto orphan = std::find ( orphans.begin ( ), orphans.end ( ), *(DSS.begin ( ) + indexOfPointInDSS) );
196  if ( orphan != orphans.end ( ) )
197  orphans.erase ( orphan );
198  }
199  }
200  if ( ! orphans.empty ( ) )
201  treatOrphans ( dssSegments->begin ( ), dssSegments->end ( ), orphans.cbegin ( ), orphans.cend ( ), outValues );
202  accumulate< OutputIterator >( outValues, itb, ite, result );
203  return result;
204  }
205 
206  template < typename TSpace, typename TSegmentation, typename Functor, typename DSSFilter >
207  template <typename OutputIterator>
208  inline
209  void
210  LambdaMST3DEstimator< TSpace, TSegmentation, Functor, DSSFilter >::accumulate ( std::multimap < Point, Value > & outValues,
211  ConstIterator itb, ConstIterator ite,
212  OutputIterator & result )
213  {
214  Value prev = outValues.lower_bound( *itb )->second;
215  Value accum_prev = outValues.lower_bound( *itb )->second;
216  for ( ConstIterator itt = itb; itt != ite; ++itt )
217  {
218  typename std::multimap< Point, Value >::const_iterator it = outValues.lower_bound ( *itt );
219  typename std::multimap< Point, Value >::const_iterator it2 = outValues.upper_bound ( *itt );
220  Value tangent;
221  for (; it != it2; it++ )
222  {
223  Value partial = it->second;
224  if ( partial.first.norm() > 0. && prev.first.norm() > 0. && prev.first.cosineSimilarity ( partial.first ) > M_PI_2 )
225  partial.first = -partial.first;
226  prev = partial;
227  tangent += partial;
228  }
229  outValues.erase ( *itt );
230  // avoid tangent flapping
231  if ( accum_prev.first.norm() > 0. && tangent.first.norm() > 0. && accum_prev.first.cosineSimilarity ( tangent.first ) > M_PI_2 )
232  tangent.first = -tangent.first;
233  accum_prev = tangent;
234  if ( tangent.second != 0 )
235  *result++ = ( tangent.first / tangent.second );
236  else
237  *result++ = tangent.first;
238  }
239  }
240 }