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