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