DGtal 1.3.0
Loading...
Searching...
No Matches
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
39template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
40const 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
50template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
51inline
52DGtal::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
67template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
68inline
69DGtal::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
86template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
87inline
88DGtal::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
104template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
105inline
106DGtal::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
124template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
125inline
126DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::~FMM()
127{
128 if (myFlagIsOwning)
129 delete myPointFunctorPtr;
130}
131
132///////////////////////////////////////////////////////////////////////////////
133// Static functions :
134
135
136template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
137template <typename TIteratorOnPoints>
138void
139DGtal::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
152template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
153template <typename KSpace, typename TIteratorOnBels>
154void
155DGtal::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
177template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
178template <typename KSpace, typename TIteratorOnBels, typename TImplicitFunction>
179void
180DGtal::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
240template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
241template <typename TIteratorOnPairs>
242void
243DGtal::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
264template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
265inline
266void
267DGtal::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
275template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
276inline
277bool
278DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>
279::computeOneStep(Point& aPoint, Value& aValue)
280{
281 return addNewAcceptedPoint(aPoint, aValue);
282}
283
284template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
285inline
286typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
287DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::min() const
288{
289 return myMinValue;
290}
291
292template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
293inline
294typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
295DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::max() const
296{
297 return myMaxValue;
298}
299
300template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
301inline
302typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
303DGtal::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
319template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
320inline
321typename DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::Value
322DGtal::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
338template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
339inline
340bool
341DGtal::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
366template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
367inline
368void
369DGtal::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
382template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
383inline
384void
385DGtal::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
402template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
403inline
404bool
405DGtal::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#ifdef _MSC_VER
427 // After erasing, the iterator "it" becomes invalid and the test with "itEnd" won't work
428 it = itEnd;
429#endif //_MSC_VER
430 //it can be inserted into the set of accepted points
431 if ( insertAndSetValue( myImage, myAcceptedPoints,
432 minPair.first, minPair.second ) )
433 { //if it does not belong to the set
434 //the set of candidates is updated with
435 //the neighbors of the new accepted point
436 aPoint = minPair.first;
437 aValue = minPair.second;
438 if (aValue > myMaxValue) myMaxValue = aValue;
439 if (aValue < myMinValue) myMinValue = aValue;
440 update( aPoint );
441 flagStop = true;
442 }
443 else
444 { //otherwise it has already been accepted
445 //with a smaller distance and the next candidate
446 //should be considered
447 it = myCandidatePoints.begin();
448 }
449
450 }//end if distance below a given threshold
451 else return false;
452
453 } //end while there are candidates
454
455 return flagStop; //true if a point has been accepted
456
457 } //end if a new point can be accepted
458 else return false;
459}
460
461template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
462inline
463void
464DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::update(const Point& aPoint)
465{
466
467 //neigbors
468 Point neighbor = aPoint;
469 for (Dimension k = 0; k < dimension; ++k)
470 {
471 typename Point::Coordinate c = neighbor[k];
472 neighbor[k] = (c+1);
473 addNewCandidate(neighbor);
474 neighbor[k] = (c-1);
475 addNewCandidate(neighbor);
476 neighbor[k] = c;
477 }
478}
479
480template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
481inline
482bool
483DGtal::FMM<TImage, TSet, TPointPredicate, TPointFunctor>::addNewCandidate(const Point& aPoint)
484{
485
486 //if it lies within the computation domain
487 //and if it is not already accepted
488 if ( (myPointPredicate(aPoint) )
489 && ( myAcceptedPoints.find(aPoint) == myAcceptedPoints.end() ) )
490 {
491 ASSERT( myPointFunctorPtr );
492 Value d = myPointFunctorPtr->operator()( aPoint );
493 PointValue newPair( aPoint, d );
494 //insert the new candidate with its distance
495 myCandidatePoints.insert(newPair);
496 return true;
497 }
498 else return false;
499}
500
501
502///////////////////////////////////////////////////////////////////////////////
503// Implementation of inline functions //
504
505template <typename TImage, typename TSet, typename TPointPredicate, typename TPointFunctor >
506inline
507std::ostream&
508DGtal::operator<< ( std::ostream & out,
509 const FMM<TImage, TSet, TPointPredicate, TPointFunctor> & object )
510{
511 object.selfDisplay( out );
512 return out;
513}
514
515// //
516///////////////////////////////////////////////////////////////////////////////
517
518
Generated on Sun Nov 27 2022 15:15:42 for DGtal by doxygen 1.9.5