DGtal  1.5.beta
VoronoiMapComplete.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
19  *
20  * @author Robin Lamy
21  * @author David Coeurjolly (\c david.coeurjolly@liris.cnrs.fr )
22  * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS,
23  UMR 5205), CNRS, France
24  *
25  * initial 2012/08/14
26  * @date 2021/07/20
27 
28  * Implementation of inline methods defined in VoronoiMap.h
29  *
30  * This file is part of the DGtal library.
31  */
32 
33 //////////////////////////////////////////////////////////////////////////////
34 #include <cstdlib>
35 #include "DGtal/kernel/NumberTraits.h"
36 
37 ///////////////////////////////////////////////////////////////////////////////
38 // IMPLEMENTATION of inline methods.
39 ///////////////////////////////////////////////////////////////////////////////
40 // ----------------------- Standard services ------------------------------
41 
42 template <typename S, typename P, typename TSep, typename TImage>
43 inline void DGtal::VoronoiMapComplete<S, P, TSep, TImage>::compute()
44 {
45  // We copy the image extent
46  myLowerBoundCopy = myDomainPtr->lowerBound();
47  myUpperBoundCopy = myDomainPtr->upperBound();
48 
49  // Point outside the domain
50  for ( auto & coord : myInfinity )
51  coord = DGtal::NumberTraits<typename Point::Coordinate>::max();
52 
53  // Init
54  for ( auto const & pt : *myDomainPtr )
55  if ( ( *myPointPredicatePtr )( pt ) )
56  {
57  Value myVectorInfinity;
58  myVectorInfinity.insert( myInfinity );
59  myImagePtr->setValue( pt, myVectorInfinity );
60  }
61  else
62  {
63  Value external_point;
64  external_point.insert( pt );
65  myImagePtr->setValue( pt, external_point );
66  }
67 
68  // We process the remaining dimensions
69  for ( Dimension dim = 0; dim < S::dimension; dim++ )
70  {
71  computeOtherSteps( dim );
72  }
73 }
74 
75 template <typename S, typename P, typename TSep, typename TImage>
76 inline void DGtal::VoronoiMapComplete<S, P, TSep, TImage>::computeOtherSteps(
77 const Dimension dim ) const
78 {
79 #ifdef VERBOSE
80  std::string title = "VoronoiMapComplete dimension " + std::to_string( dim );
81  trace.beginBlock( title );
82 #endif
83 
84  // We setup the subdomain iterator
85  // the iterator will scan dimension using the order:
86  // {n-1, n-2, ... 1} (we skip the '0' dimension).
87  std::vector<Dimension> subdomain;
88  subdomain.reserve( S::dimension - 1 );
89  for ( int k = 0; k < (int)S::dimension; k++ )
90  if ( static_cast<Dimension>( ( (int)S::dimension - 1 - k ) ) != dim )
91  subdomain.push_back( (int)S::dimension - 1 - k );
92 
93  Domain localDomain( myLowerBoundCopy, myUpperBoundCopy );
94 
95 #ifdef WITH_OPENMP
96  // Parallel loop
97  std::vector<Point> subRangePoints;
98  // Starting point precomputation
99  for ( auto const & pt : localDomain.subRange( subdomain ) )
100  subRangePoints.push_back( pt );
101 
102  // We run the 1D problems in //
103 #pragma omp parallel for schedule( dynamic )
104  for ( int i = 0; i < static_cast<int>(subRangePoints.size()); ++i ) //MSVC requires signed type for openmp
105  computeOtherStep1D( subRangePoints[ i ], dim );
106 
107 #else
108  // We solve the 1D problems sequentially
109  for ( auto const & pt : localDomain.subRange( subdomain ) )
110  computeOtherStep1D( pt, dim );
111 #endif
112 
113 #ifdef VERBOSE
114  trace.endBlock();
115 #endif
116 }
117 
118 // //////////////////////////////////////////////////////////////////////:
119 // ////////////////////////// Other Phases
120 template <typename S, typename P, typename TSep, typename TImage>
121 void DGtal::VoronoiMapComplete<S, P, TSep, TImage>::computeOtherStep1D(
122 const Point & startingPoint, const Dimension dim ) const
123 {
124  ASSERT( dim < S::dimension );
125 
126  // Default starting and ending point for a cycle
127  Point startPoint = startingPoint;
128  Point endPoint = startingPoint;
129  startPoint[ dim ] = myLowerBoundCopy[ dim ];
130  endPoint[ dim ] = myUpperBoundCopy[ dim ];
131 
132  Value myVectorInfinity;
133  myVectorInfinity.emplace( myInfinity );
134 
135  // Extent along current dimension.
136  const auto extent = myUpperBoundCopy[ dim ] - myLowerBoundCopy[ dim ] + 1;
137 
138  // Site storage.
139  std::vector<Point> Sites;
140 
141  // Reserve sites storage.
142  // +1 along periodic dimension in order to store two times the site that is on
143  // break index.
144  Sites.reserve( extent + ( isPeriodic( dim ) ? 1 : 0 ) );
145 
146  // Pruning the list of sites and defining cycle bounds.
147  // In the periodic case, the cycle bounds depend on the so-called break index
148  // that defines the start point.
149  if ( dim == 0 )
150  {
151  // For dim = 0, no sites are hidden.
152  for ( auto point = startPoint; point[ dim ] <= myUpperBoundCopy[ dim ];
153  ++point[ dim ] )
154  {
155  const Value psite = myImagePtr->operator()( point );
156  if ( psite != myVectorInfinity )
157  {
158  for ( Point voronoiPoint : psite )
159  {
160  Sites.push_back( voronoiPoint );
161  }
162  }
163  }
164 
165  // If no sites are found, then there is nothing to do.
166  if ( Sites.size() == 0 )
167  return;
168 
169  // In the periodic case and along the first dimension, the break index
170  // is at the first site found.
171  if ( isPeriodic( dim ) )
172  {
173  startPoint[ dim ] = Sites[ 0 ][ dim ];
174  endPoint[ dim ] = startPoint[ dim ] + extent - 1;
175 
176  // The first site is also the last site (with appropriate shift).
177  Sites.push_back( Sites[ 0 ] + Point::base( dim, extent ) );
178  }
179  }
180  else
181  {
182  // In the periodic case, the cycle depends on break index
183  if ( isPeriodic( dim ) )
184  {
185  // Along other than the first dimension, the break index is at the lowest
186  // site found.
187  auto minRawDist =
188  DGtal::NumberTraits<typename SeparableMetric::RawValue>::max();
189 
190  for ( auto point = startPoint; point[ dim ] <= myUpperBoundCopy[ dim ];
191  ++point[ dim ] )
192  {
193  const Value psite = myImagePtr->operator()( point );
194 
195  if ( psite != myVectorInfinity )
196  {
197  const auto rawDist = myMetricPtr->rawDistance(
198  point, *( psite.begin() ) ); // All distances to Voronoi sites are
199  // equal, we take the first one
200  if ( rawDist < minRawDist )
201  {
202  minRawDist = rawDist;
203  startPoint[ dim ] = point[ dim ];
204  }
205  }
206  }
207 
208  // If no sites are found, then there is nothing to do.
209  if ( minRawDist ==
210  DGtal::NumberTraits<typename SeparableMetric::RawValue>::max() )
211  return;
212 
213  endPoint[ dim ] = startPoint[ dim ] + extent - 1;
214  }
215 
216  // Pruning the list of sites for both periodic and non-periodic cases.
217  for ( auto point = startPoint; point[ dim ] <= myUpperBoundCopy[ dim ];
218  ++point[ dim ] )
219  {
220  const Value psite = myImagePtr->operator()( point );
221 
222  if ( psite != myVectorInfinity )
223  {
224  for ( Point site : psite )
225  {
226  while ( ( Sites.size() >= 2 ) &&
227  ( myMetricPtr->hiddenBy( Sites[ Sites.size() - 2 ],
228  Sites[ Sites.size() - 1 ], site,
229  startingPoint, endPoint, dim ) ) )
230  {
231 
232  Sites.pop_back();
233  }
234 
235  Sites.push_back( site );
236  }
237  }
238  }
239 
240  // Pruning the remaining list of sites in the periodic case.
241  if ( isPeriodic( dim ) )
242  {
243  auto point = startPoint;
244  point[ dim ] = myLowerBoundCopy[ dim ];
245  for ( ; point[ dim ] <= endPoint[ dim ] - extent + 1;
246  ++point[ dim ] ) // +1 in order to add the break-index site at the
247  // cycle's end.
248  {
249  const Value psite = myImagePtr->operator()( point );
250 
251  if ( psite != myVectorInfinity )
252  {
253  for ( Point site : psite )
254  {
255  // Site coordinates must be between startPoint and endPoint.
256  site[ dim ] += extent;
257 
258  while ( ( Sites.size() >= 2 ) &&
259  ( myMetricPtr->hiddenBy( Sites[ Sites.size() - 2 ],
260  Sites[ Sites.size() - 1 ], site,
261  startingPoint, endPoint, dim ) ) )
262  Sites.pop_back();
263 
264  Sites.push_back( site );
265  }
266  }
267  }
268  }
269  }
270 
271  // No sites found
272  if ( Sites.size() == 0 )
273  return;
274 
275  // Rewriting for both periodic and non-periodic cases.
276  size_t siteId = 0;
277  auto point = startPoint;
278  IntegerLong minimal_distance;
279 
280  for ( ; point[ dim ] <= myUpperBoundCopy[ dim ]; ++point[ dim ] )
281  {
282  Value voronoiSites = myImagePtr->operator()( point );
283  minimal_distance =
284  myMetricPtr->rawDistance( *( voronoiSites.begin() ), point );
285 
286  for ( Point site : Sites )
287  {
288  if ( myMetricPtr->rawDistance( site, point ) < minimal_distance )
289  {
290  voronoiSites.clear();
291  voronoiSites.emplace( site );
292  minimal_distance = myMetricPtr->rawDistance( site, point );
293  }
294  else if ( ( myMetricPtr->rawDistance( site, point ) ==
295  minimal_distance ) &&
296  ( std::find( voronoiSites.begin(), voronoiSites.end(), site ) ==
297  voronoiSites.end() ) )
298  {
299  voronoiSites.emplace( site );
300  }
301  }
302  myImagePtr->setValue( point, voronoiSites );
303  }
304 
305  // Continuing rewriting in the periodic case.
306  if ( isPeriodic( dim ) )
307  {
308  for ( ; point[ dim ] <= endPoint[ dim ]; ++point[ dim ] )
309  {
310  Value voronoiSites;
311  siteId = 0;
312 
313  while ( siteId < Sites.size() - 1 )
314  { // Idem ici en prenant en compte la périodicité
315  if ( myMetricPtr->closest( point, Sites[ siteId ],
316  Sites[ siteId + 1 ] ) ==
317  DGtal::ClosestFIRST )
318  {
319  voronoiSites.insert( Sites[ siteId ] - Point::base( dim, extent ) );
320  }
321  siteId++;
322  }
323  myImagePtr->setValue( point - Point::base( dim, extent ), voronoiSites );
324  }
325  }
326 }
327 
328 /**
329  * Constructor.
330  */
331 template <typename S, typename P, typename TSep, typename TImage>
332 inline DGtal::VoronoiMapComplete<S, P, TSep, TImage>::VoronoiMapComplete(ConstAlias<Domain> aDomain,
333  ConstAlias<PointPredicate> aPredicate,
334  ConstAlias<SeparableMetric> aMetric )
335  : myDomainPtr( &aDomain )
336  , myPointPredicatePtr( &aPredicate )
337  , myDomainExtent( aDomain->upperBound() - aDomain->lowerBound() +
338  Point::diagonal( 1 ) )
339  , myMetricPtr( &aMetric )
340 {
341  myPeriodicitySpec.fill( false );
342  myImagePtr = CountedPtr<OutputImage>( new OutputImage( aDomain ) );
343  compute();
344 }
345 
346 template <typename S, typename P, typename TSep, typename TImage>
347 inline DGtal::VoronoiMapComplete<S, P, TSep, TImage>::VoronoiMapComplete(ConstAlias<Domain> aDomain,
348  ConstAlias<PointPredicate> aPredicate,
349  ConstAlias<SeparableMetric> aMetric,
350  PeriodicitySpec const & aPeriodicitySpec )
351  : myDomainPtr( &aDomain )
352  , myPointPredicatePtr( &aPredicate )
353  , myDomainExtent( aDomain->upperBound() - aDomain->lowerBound() +
354  Point::diagonal( 1 ) )
355  , myMetricPtr( &aMetric )
356  , myPeriodicitySpec( aPeriodicitySpec )
357 {
358  // Finding periodic dimension index.
359  for ( Dimension i = 0; i < Space::dimension; ++i )
360  if ( isPeriodic( i ) )
361  myPeriodicityIndex.push_back( i );
362 
363  myImagePtr = CountedPtr<OutputImage>( new OutputImage( aDomain ) );
364  compute();
365 }
366 
367 template <typename S, typename P, typename TSep, typename TImage>
368 inline typename DGtal::VoronoiMapComplete<S, P, TSep, TImage>::Point
369 DGtal::VoronoiMapComplete<S, P, TSep, TImage>::projectPoint(
370 Point aPoint ) const
371 {
372  for ( auto const & dim : myPeriodicityIndex )
373  aPoint[ dim ] = projectCoordinate( aPoint[ dim ], dim );
374 
375  return aPoint;
376 }
377 
378 template <typename S, typename P, typename TSep, typename TImage>
379 inline typename DGtal::VoronoiMapComplete<S, P, TSep, TImage>::Point::Coordinate
380 DGtal::VoronoiMapComplete<S, P, TSep, TImage>::projectCoordinate(typename Point::Coordinate aCoordinate,
381  const Dimension aDim ) const
382 {
383  ASSERT( aCoordinate - myDomainPtr->lowerBound()[ aDim ] +
384  myDomainExtent[ aDim ] >=
385  0 );
386  return ( aCoordinate - myDomainPtr->lowerBound()[ aDim ] +
387  myDomainExtent[ aDim ] ) %
388  myDomainExtent[ aDim ] +
389  myDomainPtr->lowerBound()[ aDim ];
390 }
391 
392 template <typename S, typename P, typename TSep, typename TImage>
393 inline void DGtal::VoronoiMapComplete<S, P, TSep, TImage>::selfDisplay(
394 std::ostream & out ) const
395 {
396  out << "[VoronoiMapComplete] separable metric=" << *myMetricPtr;
397 }
398 
399 // // //
400 // ///////////////////////////////////////////////////////////////////////////////
401 
402 template <typename S, typename P, typename TSep, typename TImage>
403 inline std::ostream &
404 DGtal::operator<<( std::ostream & out,
405  const VoronoiMapComplete<S, P, TSep, TImage> & object )
406 {
407  object.selfDisplay( out );
408  return out;
409 }