DGtal  1.0.0
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 std::pair<const Bel, Value> BelValue;
193  typedef std::map<Bel, Value> Buffer;
194 
195  typedef typename functors::SCellToIncidentPoints<KSpace>::Output Pair;
196  typedef std::vector<Pair> IncidentPoints;
197  functors::SCellToIncidentPoints<KSpace> getIncidentPoints( aK );
198 
199  /// 1) store a value for each bel
200  /// (== distance of the inner point to the interface)
201  IncidentPoints incidentPoints;
202  Buffer buffer;
203  for (TIteratorOnBels it = itb; it != ite; ++it)
204  {
205  //getting incident points
206  Pair points = getIncidentPoints( *it );
207  incidentPoints.push_back( points );
208 
209  //getting values
210  Value vin = aF( points.first );
211  Value vout = aF( points.second );
212 
213  //computing/storing the new value
214  Value e = std::max(vin, vout) - std::min(vin, vout);
215  Value v = (std::abs(vin)/e);
216  ASSERT( v >= 0 );
217  buffer.insert( BelValue( aK.unsigns( *it ), v ) );
218  }
219 
220  aSet.clear();
221  /// 2) for each inner/outer point
222  /// computes its distance from the values of its incident bels
223  typedef L2FirstOrderLocalDistanceFromCells<KSpace, Buffer, false> Computer4InnerPts;
224  typedef L2FirstOrderLocalDistanceFromCells<KSpace, Buffer, true> Computer4OuterPts;
225  Computer4InnerPts computerIn(aK, buffer);
226  Computer4OuterPts computerOut(aK, buffer);
227  for (typename IncidentPoints::const_iterator it = incidentPoints.begin();
228  it != incidentPoints.end(); ++it)
229  {
230  //computing the values
231  Value vin = computerIn( it->first );
232  Value vout = computerOut( it->second );
233 
234  //assignement
235  insertAndAlwaysSetValue( aImg, aSet, it->first, -k*vin );
236  insertAndAlwaysSetValue( aImg, aSet, it->second, k*vout );
237  }
238 }
239 
240 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
241 template <typename TIteratorOnPairs>
242 void
243 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
244 ::initFromIncidentPointsRange(const TIteratorOnPairs& itb, const TIteratorOnPairs& ite,
245  Image& aImg, AcceptedPointSet& aSet,
246  const Value& aValue,
247  bool aFlagIsPositive)
248 {
249  Value k = -1;
250  if (aFlagIsPositive) k = 1;
251 
252  aSet.clear();
253  for (TIteratorOnPairs it = itb; it != ite; ++it)
254  {
255  insertAndAlwaysSetValue( aImg, aSet, it->first, -k*aValue );
256  insertAndAlwaysSetValue( aImg, aSet, it->second, k*aValue );
257  }
258 }
259 
260 ///////////////////////////////////////////////////////////////////////////////
261 // Interface - public :
262 
263 
264 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
265 inline
266 void
267 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::compute()
268 {
269  Point p = Point::diagonal(0);
270  Value d = 0;
271  while ( addNewAcceptedPoint( p, d ) )
272  { }
273 }
274 
275 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
276 inline
277 bool
278 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
279 ::computeOneStep(Point& aPoint, Value& aValue)
280 {
281  return addNewAcceptedPoint(aPoint, aValue);
282 }
283 
284 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
285 inline
286 typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
287 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::min() const
288 {
289  return myMinValue;
290 }
291 
292 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
293 inline
294 typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
295 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::max() const
296 {
297  return myMaxValue;
298 }
299 
300 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
301 inline
302 typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
303 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::getMin() const
304 {
305  const AcceptedPointSet& set = myAcceptedPoints;
306  ASSERT( set.size() >= 1 );
307 
308  typename AcceptedPointSet::ConstIterator it = set.begin();
309  typename AcceptedPointSet::ConstIterator itEnd = set.end();
310  Value vmin = myImage( *it );
311  for (++it; it != itEnd; ++it)
312  {
313  Value v = myImage( *it );
314  if (v < vmin) vmin = v;
315  }
316  return vmin;
317 }
318 
319 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
320 inline
321 typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
322 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::getMax() const
323 {
324  const AcceptedPointSet& set = myAcceptedPoints;
325  ASSERT( set.size() >= 1 );
326 
327  typename AcceptedPointSet::ConstIterator it = set.begin();
328  typename AcceptedPointSet::ConstIterator itEnd = set.end();
329  Value vmax = myImage( *it );
330  for (++it; it != itEnd; ++it)
331  {
332  Value v = myImage( *it );
333  if (v > vmax) vmax = v;
334  }
335  return vmax;
336 }
337 
338 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
339 inline
340 bool
341 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::isValid() const
342 {
343  //area threshold
344  if ( (myAcceptedPoints.size() <= 0)
345  || (myAcceptedPoints.size() >= myAreaThreshold) ) return false;
346 
347  //distance threshold
348  if ( ( getMin() != min() ) || ( getMax() != max() ) ) return false;
349  if ( (std::abs(getMin()) >= myValueThreshold)
350  || (getMax() >= myValueThreshold) ) return false;
351 
352  //point predicate
353  bool flagIsOk = true;
354  const AcceptedPointSet& set = myAcceptedPoints;
355  typename AcceptedPointSet::ConstIterator it = set.begin();
356  typename AcceptedPointSet::ConstIterator itEnd = set.end();
357  for ( ; ( (it != itEnd)&&(flagIsOk == true) ); ++it)
358  {
359  if (myPointPredicate( *it ) == false) flagIsOk = false;
360  }
361  if (!flagIsOk) return false;
362 
363  return true;
364 }
365 
366 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
367 inline
368 void
369 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::selfDisplay ( std::ostream & out ) const
370 {
371  out << "[FMM " << dimension << "d] ";
372  out << myAcceptedPoints.size() << " accepted points (< " << myAreaThreshold << ")";
373  out << " and " << myCandidatePoints.size() << " candidates. ";
374  out << "dmin: " << min() << ", dmax: " << max();
375  out << " (abs < " << myValueThreshold << ")";
376 }
377 
378 
379 ///////////////////////////////////////////////////////////////////////////////
380 // Internals
381 
382 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
383 inline
384 void
385 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::init()
386 {
387 
388  myCandidatePoints.clear();
389 
390  typename AcceptedPointSet::Iterator it = myAcceptedPoints.begin();
391  typename AcceptedPointSet::Iterator itEnd = myAcceptedPoints.end();
392  for ( ; it != itEnd; ++it)
393  {
394  update( *it );
395  }
396 
397  myMinValue = getMin();
398  myMaxValue = getMax();
399 
400 }
401 
402 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
403 inline
404 bool
405 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
406 ::addNewAcceptedPoint(Point& aPoint, Value& aValue)
407 {
408 
409  if ( (myAcceptedPoints.size()+1) < myAreaThreshold )
410  {//if a new point can be accepted
411 
412  bool flagStop = false;
413  typename CandidatePointSet::iterator it = myCandidatePoints.begin();
414  typename CandidatePointSet::iterator itEnd = myCandidatePoints.end();
415  while ( (it != itEnd) && (!flagStop) )
416  { //while there are candidates and no point has been accepted
417 
418  //pair of min distance
419  PointValue minPair = *it;
420 
421  if ( std::abs(minPair.second) < myValueThreshold )
422  { //if distance below a given threshold
423 
424  //the point of min distance is removed from the set of candidates
425  myCandidatePoints.erase(*it);
426  //it can be inserted into the set of accepted points
427  if ( insertAndSetValue( myImage, myAcceptedPoints,
428  minPair.first, minPair.second ) )
429  { //if it does not belong to the set
430  //the set of candidates is updated with
431  //the neighbors of the new accepted point
432  aPoint = minPair.first;
433  aValue = minPair.second;
434  if (aValue > myMaxValue) myMaxValue = aValue;
435  if (aValue < myMinValue) myMinValue = aValue;
436  update( aPoint );
437  flagStop = true;
438  }
439  else
440  { //otherwise it has already been accepted
441  //with a smaller distance and the next candidate
442  //should be considered
443  it = myCandidatePoints.begin();
444  }
445 
446  }//end if distance below a given threshold
447  else return false;
448 
449  } //end while there are candidates
450 
451  return flagStop; //true if a point has been accepted
452 
453  } //end if a new point can be accepted
454  else return false;
455 }
456 
457 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
458 inline
459 void
460 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::update(const Point& aPoint)
461 {
462 
463  //neigbors
464  Point neighbor = aPoint;
465  for (Dimension k = 0; k < dimension; ++k)
466  {
467  typename Point::Coordinate c = neighbor[k];
468  neighbor[k] = (c+1);
469  addNewCandidate(neighbor);
470  neighbor[k] = (c-1);
471  addNewCandidate(neighbor);
472  neighbor[k] = c;
473  }
474 }
475 
476 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
477 inline
478 bool
479 DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::addNewCandidate(const Point& aPoint)
480 {
481 
482  //if it lies within the computation domain
483  //and if it is not already accepted
484  if ( (myPointPredicate(aPoint) )
485  && ( myAcceptedPoints.find(aPoint) == myAcceptedPoints.end() ) )
486  {
487  ASSERT( myPointFunctorPtr );
488  Value d = myPointFunctorPtr->operator()( aPoint );
489  PointValue newPair( aPoint, d );
490  //insert the new candidate with its distance
491  myCandidatePoints.insert(newPair);
492  return true;
493  }
494  else return false;
495 }
496 
497 
498 ///////////////////////////////////////////////////////////////////////////////
499 // Implementation of inline functions //
500 
501 template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
502 inline
503 std::ostream&
504 DGtal::operator<< ( std::ostream & out,
505  const FMM<TImage, TSet, TPointPredicate, TPointFunctor> & object )
506 {
507  object.selfDisplay( out );
508  return out;
509 }
510 
511 // //
512 ///////////////////////////////////////////////////////////////////////////////
513 
514