DGtal  0.9.2
FMM.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 FMM.ih
19  * @author Tristan Roussillon (\c
20  * tristan.roussillon@liris.cnrs.fr ) Laboratoire d'InfoRmatique en
21  * Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS,
22  * France
23  *
24  *
25  * @date 2012/01/17
26  *
27  * @brief Implementation of inline methods defined in FMM.h
28  *
29  * This file is part of the DGtal library.
30  */
31 
32 
33 //////////////////////////////////////////////////////////////////////////////
34 #include <cstdlib>
35 //////////////////////////////////////////////////////////////////////////////
36 
37 #include "DGtal/topology/SCellsFunctors.h"
38 
39 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
40 const typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Dimension DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::dimension = Point::dimension;
41 
42 
43 ///////////////////////////////////////////////////////////////////////////////
44 // IMPLEMENTATION of inline methods.
45 ///////////////////////////////////////////////////////////////////////////////
46 
47 ///////////////////////////////////////////////////////////////////////////////
48 // ----------------------- Standard services ------------------------------
49 
50 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
51 inline
52 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
53 ::FMM(Image& aImg, AcceptedPointSet& aSet,
54  ConstAlias<PointPredicate> aPointPredicate)
55  : myImage( aImg ), myAcceptedPoints( aSet ),
56  myPointFunctorPtr( new PointFunctor(aImg, aSet) ),
57  myFlagIsOwning( true ),
58  myPointPredicate( aPointPredicate ),
59  myAreaThreshold( std::numeric_limits<Area>::max() ),
60  myValueThreshold( std::numeric_limits<Value>::max() )
61 {
62  if (myAcceptedPoints.size() == 0) throw InputException();
63  init();
64 }
65 
66 
67 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
68 inline
69 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
70 ::FMM(Image& aImg, AcceptedPointSet& aSet,
71  ConstAlias<PointPredicate> aPointPredicate,
72  const Area& aAreaThreshold,
73  const Value& aValueThreshold)
74  : myImage( aImg ), myAcceptedPoints( aSet ),
75  myPointFunctorPtr( new PointFunctor(aImg, aSet) ),
76  myFlagIsOwning( true ),
77  myPointPredicate( aPointPredicate ),
78  myAreaThreshold( aAreaThreshold ),
79  myValueThreshold( aValueThreshold )
80 {
81  if (myAcceptedPoints.size() == 0) throw InputException();
82  init();
83 }
84 
85 
86 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
87 inline
88 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
89 ::FMM(Image& aImg, AcceptedPointSet& aSet,
90  ConstAlias<PointPredicate> aPointPredicate,
91  PointFunctor& aPointFunctor)
92  : myImage( aImg ), myAcceptedPoints( aSet ),
93  myPointFunctorPtr( &aPointFunctor ),
94  myFlagIsOwning( false ),
95  myPointPredicate( aPointPredicate ),
96  myAreaThreshold( std::numeric_limits<Area>::max() ),
97  myValueThreshold( std::numeric_limits<Value>::max() )
98 {
99  if (myAcceptedPoints.size() == 0) throw InputException();
100  init();
101 }
102 
103 
104 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
105 inline
106 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
107 ::FMM(Image& aImg, AcceptedPointSet& aSet,
108  ConstAlias<PointPredicate> aPointPredicate,
109  const Area& aAreaThreshold,
110  const Value& aValueThreshold,
111  PointFunctor& aPointFunctor)
112  : myImage( aImg ), myAcceptedPoints( aSet ),
113  myPointFunctorPtr( &aPointFunctor ),
114  myFlagIsOwning( false ),
115  myPointPredicate( aPointPredicate ),
116  myAreaThreshold( aAreaThreshold ),
117  myValueThreshold( aValueThreshold )
118 {
119  if (myAcceptedPoints.size() == 0) throw InputException();
120  init();
121 }
122 
123 
124 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
125 inline
126 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::~FMM()
127 {
128  if (myFlagIsOwning)
129  delete myPointFunctorPtr;
130 }
131 
132 ///////////////////////////////////////////////////////////////////////////////
133 // Static functions :
134 
135 
136 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
137 template <typename TIteratorOnPoints>
138 void
139 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
140 ::initFromPointsRange(const TIteratorOnPoints& itb, const TIteratorOnPoints& ite,
141  Image& aImg, AcceptedPointSet& aSet,
142  const Value& aValue)
143 {
144 
145  aSet.clear();
146  for (TIteratorOnPoints it = itb; it != ite; ++it)
147  {
148  insertAndAlwaysSetValue( aImg, aSet, *it, aValue );
149  }
150 }
151 
152 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
153 template <typename KSpace, typename TIteratorOnBels>
154 void
155 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
156 ::initFromBelsRange(const KSpace& aK,
157  const TIteratorOnBels& itb, const TIteratorOnBels& ite,
158  Image& aImg, AcceptedPointSet& aSet,
159  const Value& aValue,
160  bool aFlagIsPositive)
161 {
162  Value k = -1;
163  if (aFlagIsPositive) k = 1;
164 
165  aSet.clear();
166  for (TIteratorOnBels it = itb; it != ite; ++it)
167  {
168  //getting incident points
169  functors::SCellToIncidentPoints<KSpace> func( aK );
170  typename functors::SCellToIncidentPoints<KSpace>::Output points = func( *it );
171  //assignement
172  insertAndAlwaysSetValue( aImg, aSet, points.first, -k*aValue );
173  insertAndAlwaysSetValue( aImg, aSet, points.second, k*aValue );
174  }
175 }
176 
177 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
178 template <typename KSpace, typename TIteratorOnBels, typename TImplicitFunction>
179 void
180 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
181 ::initFromBelsRange(const KSpace& aK,
182  const TIteratorOnBels& itb, const TIteratorOnBels& ite,
183  const TImplicitFunction& aF,
184  Image& aImg, AcceptedPointSet& aSet,
185  bool aFlagIsPositive)
186 {
187  Value k = -1;
188  if (aFlagIsPositive) k = 1;
189 
190  /// types
191  typedef typename KSpace::Cell Bel;
192  typedef typename TImplicitFunction::Value Value;
193  typedef std::pair<const Bel, Value> BelValue;
194  typedef std::map<Bel, Value> Buffer;
195 
196  typedef typename functors::SCellToIncidentPoints<KSpace>::Output Pair;
197  typedef std::vector<Pair> IncidentPoints;
198  functors::SCellToIncidentPoints<KSpace> getIncidentPoints( aK );
199 
200  /// 1) store a value for each bel
201  /// (== distance of the inner point to the interface)
202  IncidentPoints incidentPoints;
203  Buffer buffer;
204  for (TIteratorOnBels it = itb; it != ite; ++it)
205  {
206  //getting incident points
207  Pair points = getIncidentPoints( *it );
208  incidentPoints.push_back( points );
209 
210  //getting values
211  Value vin = aF( points.first );
212  Value vout = aF( points.second );
213 
214  //computing/storing the new value
215  Value e = std::max(vin, vout) - std::min(vin, vout);
216  Value v = (std::abs(vin)/e);
217  ASSERT( v >= 0 );
218  buffer.insert( BelValue( aK.unsigns( *it ), v ) );
219  }
220 
221  aSet.clear();
222  /// 2) for each inner/outer point
223  /// computes its distance from the values of its incident bels
224  typedef L2FirstOrderLocalDistanceFromCells<KSpace, Buffer, false> Computer4InnerPts;
225  typedef L2FirstOrderLocalDistanceFromCells<KSpace, Buffer, true> Computer4OuterPts;
226  Computer4InnerPts computerIn(aK, buffer);
227  Computer4OuterPts computerOut(aK, buffer);
228  for (typename IncidentPoints::const_iterator it = incidentPoints.begin();
229  it != incidentPoints.end(); ++it)
230  {
231  //computing the values
232  Value vin = computerIn( it->first );
233  Value vout = computerOut( it->second );
234 
235  //assignement
236  insertAndAlwaysSetValue( aImg, aSet, it->first, -k*vin );
237  insertAndAlwaysSetValue( aImg, aSet, it->second, k*vout );
238  }
239 }
240 
241 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
242 template <typename TIteratorOnPairs>
243 void
244 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
245 ::initFromIncidentPointsRange(const TIteratorOnPairs& itb, const TIteratorOnPairs& ite,
246  Image& aImg, AcceptedPointSet& aSet,
247  const Value& aValue,
248  bool aFlagIsPositive)
249 {
250  Value k = -1;
251  if (aFlagIsPositive) k = 1;
252 
253  aSet.clear();
254  for (TIteratorOnPairs it = itb; it != ite; ++it)
255  {
256  insertAndAlwaysSetValue( aImg, aSet, it->first, -k*aValue );
257  insertAndAlwaysSetValue( aImg, aSet, it->second, k*aValue );
258  }
259 }
260 
261 ///////////////////////////////////////////////////////////////////////////////
262 // Interface - public :
263 
264 
265 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
266 inline
267 void
268 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::compute()
269 {
270  Point p = Point::diagonal(0);
271  Value d = 0;
272  while ( addNewAcceptedPoint( p, d ) )
273  { }
274 }
275 
276 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
277 inline
278 bool
279 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
280 ::computeOneStep(Point& aPoint, Value& aValue)
281 {
282  return addNewAcceptedPoint(aPoint, aValue);
283 }
284 
285 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
286 inline
287 typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
288 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::min() const
289 {
290  return myMinValue;
291 }
292 
293 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
294 inline
295 typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
296 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::max() const
297 {
298  return myMaxValue;
299 }
300 
301 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
302 inline
303 typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
304 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::getMin() const
305 {
306  const AcceptedPointSet& set = myAcceptedPoints;
307  ASSERT( set.size() >= 1 );
308 
309  typename AcceptedPointSet::ConstIterator it = set.begin();
310  typename AcceptedPointSet::ConstIterator itEnd = set.end();
311  Value vmin = myImage( *it );
312  for (++it; it != itEnd; ++it)
313  {
314  Value v = myImage( *it );
315  if (v < vmin) vmin = v;
316  }
317  return vmin;
318 }
319 
320 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
321 inline
322 typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
323 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::getMax() const
324 {
325  const AcceptedPointSet& set = myAcceptedPoints;
326  ASSERT( set.size() >= 1 );
327 
328  typename AcceptedPointSet::ConstIterator it = set.begin();
329  typename AcceptedPointSet::ConstIterator itEnd = set.end();
330  Value vmax = myImage( *it );
331  for (++it; it != itEnd; ++it)
332  {
333  Value v = myImage( *it );
334  if (v > vmax) vmax = v;
335  }
336  return vmax;
337 }
338 
339 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
340 inline
341 bool
342 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::isValid() const
343 {
344  //area threshold
345  if ( (myAcceptedPoints.size() <= 0)
346  || (myAcceptedPoints.size() >= myAreaThreshold) ) return false;
347 
348  //distance threshold
349  if ( ( getMin() != min() ) || ( getMax() != max() ) ) return false;
350  if ( (std::abs(getMin()) >= myValueThreshold)
351  || (getMax() >= myValueThreshold) ) return false;
352 
353  //point predicate
354  bool flagIsOk = true;
355  const AcceptedPointSet& set = myAcceptedPoints;
356  typename AcceptedPointSet::ConstIterator it = set.begin();
357  typename AcceptedPointSet::ConstIterator itEnd = set.end();
358  for ( ; ( (it != itEnd)&&(flagIsOk == true) ); ++it)
359  {
360  if (myPointPredicate( *it ) == false) flagIsOk = false;
361  }
362  if (!flagIsOk) return false;
363 
364  return true;
365 }
366 
367 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
368 inline
369 void
370 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::selfDisplay ( std::ostream & out ) const
371 {
372  out << "[FMM " << dimension << "d] ";
373  out << myAcceptedPoints.size() << " accepted points (< " << myAreaThreshold << ")";
374  out << " and " << myCandidatePoints.size() << " candidates. ";
375  out << "dmin: " << min() << ", dmax: " << max();
376  out << " (abs < " << myValueThreshold << ")";
377 }
378 
379 
380 ///////////////////////////////////////////////////////////////////////////////
381 // Internals
382 
383 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
384 inline
385 void
386 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::init()
387 {
388 
389  myCandidatePoints.clear();
390 
391  typename AcceptedPointSet::Iterator it = myAcceptedPoints.begin();
392  typename AcceptedPointSet::Iterator itEnd = myAcceptedPoints.end();
393  for ( ; it != itEnd; ++it)
394  {
395  update( *it );
396  }
397 
398  myMinValue = getMin();
399  myMaxValue = getMax();
400 
401 }
402 
403 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
404 inline
405 bool
406 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
407 ::addNewAcceptedPoint(Point& aPoint, Value& aValue)
408 {
409 
410  if ( (myAcceptedPoints.size()+1) < myAreaThreshold )
411  {//if a new point can be accepted
412 
413  bool flagStop = false;
414  typename CandidatePointSet::iterator it = myCandidatePoints.begin();
415  typename CandidatePointSet::iterator itEnd = myCandidatePoints.end();
416  while ( (it != itEnd) && (!flagStop) )
417  { //while there are candidates and no point has been accepted
418 
419  //pair of min distance
420  PointValue minPair = *it;
421 
422  if ( std::abs(minPair.second) < myValueThreshold )
423  { //if distance below a given threshold
424 
425  //the point of min distance is removed from the set of candidates
426  myCandidatePoints.erase(*it);
427  //it can be inserted into the set of accepted points
428  if ( insertAndSetValue( myImage, myAcceptedPoints,
429  minPair.first, minPair.second ) )
430  { //if it does not belong to the set
431  //the set of candidates is updated with
432  //the neighbors of the new accepted point
433  aPoint = minPair.first;
434  aValue = minPair.second;
435  if (aValue > myMaxValue) myMaxValue = aValue;
436  if (aValue < myMinValue) myMinValue = aValue;
437  update( aPoint );
438  flagStop = true;
439  }
440  else
441  { //otherwise it has already been accepted
442  //with a smaller distance and the next candidate
443  //should be considered
444  it = myCandidatePoints.begin();
445  }
446 
447  }//end if distance below a given threshold
448  else return false;
449 
450  } //end while there are candidates
451 
452  return flagStop; //true if a point has been accepted
453 
454  } //end if a new point can be accepted
455  else return false;
456 }
457 
458 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
459 inline
460 void
461 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::update(const Point& aPoint)
462 {
463 
464  //neigbors
465  Point neighbor = aPoint;
466  for (Dimension k = 0; k < dimension; ++k)
467  {
468  typename Point::Coordinate c = neighbor[k];
469  neighbor[k] = (c+1);
470  addNewCandidate(neighbor);
471  neighbor[k] = (c-1);
472  addNewCandidate(neighbor);
473  neighbor[k] = c;
474  }
475 }
476 
477 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
478 inline
479 bool
480 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::addNewCandidate(const Point& aPoint)
481 {
482 
483  //if it lies within the computation domain
484  //and if it is not already accepted
485  if ( (myPointPredicate(aPoint) )
486  && ( myAcceptedPoints.find(aPoint) == myAcceptedPoints.end() ) )
487  {
488  ASSERT( myPointFunctorPtr );
489  Value d = myPointFunctorPtr->operator()( aPoint );
490  PointValue newPair( aPoint, d );
491  //insert the new candidate with its distance
492  myCandidatePoints.insert(newPair);
493  return true;
494  }
495  else return false;
496 }
497 
498 
499 ///////////////////////////////////////////////////////////////////////////////
500 // Implementation of inline functions //
501 
502 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
503 inline
504 std::ostream&
505 DGtal::operator<< ( std::ostream & out,
506  const FMM<TImage, TSet, TPointPredicate, TPointFunctor> & object )
507 {
508  object.selfDisplay( out );
509  return out;
510 }
511 
512 // //
513 ///////////////////////////////////////////////////////////////////////////////
514 
515