DGtal  1.5.beta
PowerMap.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 PowerMap.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 PowerMap.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 W, typename Sep, typename Im>
49 inline
50 void
51 DGtal::PowerMap<W, Sep,Im>::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 the map: the power map at point p is:
62  // - p if p is an input weighted point (with weight > 0);
63  // - myInfinity otherwise.
64  for( auto const & pt : *myDomainPtr )
65  if ( myWeightImagePtr->domain().isInside( pt ) &&
66  ( myWeightImagePtr->operator()( pt ) > 0 ) )
67  myImagePtr->setValue ( pt, pt );
68  else
69  myImagePtr->setValue ( pt, myInfinity );
70 
71  //We process the dimensions one by one
72  for ( Dimension dim = 0; dim < W::Domain::Space::dimension ; dim++ )
73  computeOtherSteps ( dim );
74 }
75 
76 template < typename W, typename Sep, typename Im>
77 inline
78 void
79 DGtal::PowerMap<W, Sep,Im>::computeOtherSteps ( const Dimension dim ) const
80 {
81 #ifdef VERBOSE
82  std::string title = "Powermap dimension " + boost::lexical_cast<std::string>( dim ) ;
83  trace.beginBlock ( title );
84 #endif
85 
86  //We setup the subdomain iterator
87  //the iterator will scan dimension using the order:
88  // {n-1, n-2, ... 1} (we skip the '0' dimension).
89  std::vector<Dimension> subdomain;
90  subdomain.reserve(W::Domain::Space::dimension - 1);
91  for (unsigned int k = 0; k < W::Domain::Space::dimension ; k++)
92  if ( ((int)W::Domain::Space::dimension - 1 - k) != dim)
93  subdomain.push_back( (int)W::Domain::Space::dimension - 1 - k );
94 
95  Domain localDomain(myLowerBoundCopy, myUpperBoundCopy);
96 
97 
98 #ifdef WITH_OPENMP
99  //Parallel loop
100  std::vector<Point> subRangePoints;
101  //Starting point precomputation
102  for ( auto const & pt : localDomain.subRange( subdomain ) )
103  subRangePoints.push_back( pt );
104 
105  //We run the 1D problems in //
106 #pragma omp parallel for schedule(dynamic)
107  for (int i = 0; i < (int)subRangePoints.size(); ++i) //MSCV needs signed
108  computeOtherStep1D ( subRangePoints[i], dim);
109 
110 #else
111  //We solve the 1D problems sequentially
112  for ( auto const & pt : localDomain.subRange( subdomain ) )
113  computeOtherStep1D ( pt, dim);
114 #endif
115 
116 #ifdef VERBOSE
117  trace.endBlock();
118 #endif
119 }
120 
121 // //////////////////////////////////////////////////////////////////////:
122 // ////////////////////////// Other Phases
123 template <typename W, typename Sep, typename Im>
124 void
125 DGtal::PowerMap<W,Sep,Im>::computeOtherStep1D ( const Point &startingPoint,
126  const Dimension dim) const
127 {
128  ASSERT(dim < Space::dimension);
129 
130  // Default starting and ending point for a cycle
131  Point startPoint = startingPoint;
132  Point endPoint = startingPoint;
133  startPoint[dim] = myLowerBoundCopy[dim];
134  endPoint[dim] = myUpperBoundCopy[dim];
135 
136  // Extent along current dimension.
137  const auto extent = myUpperBoundCopy[dim] - myLowerBoundCopy[dim] + 1;
138 
139  // Site storage.
140  std::vector<Point> Sites; // Site coordinates with unbounded coordinates (can be outside the domain along periodic dimensions).
141  std::vector<Point> boundedSites; // Site coordinates with bounded coordinates (always inside the domain).
142 
143  // Reserve sites storage.
144  // +1 along periodic dimension in order to store two times the site that is on break index.
145  Sites.reserve( extent + ( isPeriodic(dim) ? 1 : 0 ) );
146  boundedSites.reserve( extent + ( isPeriodic(dim) ? 1 : 0 ) );
147 
148  // Pruning the list of sites and defining cycle bounds.
149  // In the periodic case, the cycle bounds depend on the so-called break index
150  // that defines the start point.
151  if ( dim == 0 )
152  {
153  // For dim = 0, no sites are hidden.
154  for ( auto point = startPoint ; point[dim] <= myUpperBoundCopy[dim] ; ++point[dim] )
155  {
156  const Point psite = myImagePtr->operator()( point );
157  if ( psite != myInfinity )
158  {
159  Sites.push_back( psite );
160  boundedSites.push_back( psite );
161  }
162  }
163 
164  // If no sites are found, then there is nothing to do.
165  if ( Sites.size() == 0 )
166  return;
167 
168  // In the periodic case and along the first dimension, the break index
169  // is at the first site found.
170  if ( isPeriodic(dim) )
171  {
172  startPoint[dim] = Sites[0][dim];
173  endPoint[dim] = startPoint[dim] + extent - 1;
174 
175  // The first site is also the last site (with appropriate shift).
176  Sites.push_back( Sites[0] + Point::base(dim, extent) );
177  boundedSites.push_back( Sites[0] );
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 site found.
186  auto minPowerDist = DGtal::NumberTraits< typename PowerSeparableMetric::Weight >::max();
187 
188  for ( auto point = startPoint; point[dim] <= myUpperBoundCopy[dim]; ++point[dim] )
189  {
190  const Point psite = myImagePtr->operator()( point );
191 
192  if ( psite != myInfinity )
193  {
194  const auto powerDist = myMetricPtr->powerDistance( point, psite, myWeightImagePtr->operator()( projectPoint( psite, dim-1 ) ) );
195  if ( powerDist < minPowerDist )
196  {
197  minPowerDist = powerDist;
198  startPoint[dim] = point[dim];
199  }
200  }
201  }
202 
203  // If no sites are found, then there is nothing to do.
204  if ( minPowerDist == DGtal::NumberTraits< typename PowerSeparableMetric::Weight >::max() )
205  return;
206 
207  endPoint[dim] = startPoint[dim] + extent - 1;
208  }
209 
210  if ( myPeriodicityIndex.size() == 0 )
211  {
212  // Pruning the list of sites for both periodic and non-periodic cases.
213  for( auto point = startPoint ; point[dim] <= myUpperBoundCopy[dim] ; ++point[dim] )
214  {
215  const Point psite = myImagePtr->operator()(point);
216 
217  if ( psite != myInfinity )
218  {
219  while (( Sites.size() >= 2 ) &&
220  ( myMetricPtr->hiddenByPower(Sites[Sites.size()-2], myWeightImagePtr->operator()(Sites[Sites.size()-2]),
221  Sites[Sites.size()-1], myWeightImagePtr->operator()(Sites[Sites.size()-1]),
222  psite, myWeightImagePtr->operator()(psite),
223  startingPoint, endPoint, dim) ))
224  {
225  Sites.pop_back();
226  boundedSites.pop_back();
227  }
228 
229  Sites.push_back( psite );
230  boundedSites.push_back( psite );
231  }
232  }
233  }
234  else
235  {
236  // Pruning the list of sites for both periodic and non-periodic cases.
237  for( auto point = startPoint ; point[dim] <= myUpperBoundCopy[dim] ; ++point[dim] )
238  {
239  const Point psite = myImagePtr->operator()(point);
240 
241  if ( psite != myInfinity )
242  {
243  const Point boundedPSite = projectPoint( psite, dim-1 );
244 
245  while (( Sites.size() >= 2 ) &&
246  ( myMetricPtr->hiddenByPower(Sites[Sites.size()-2], myWeightImagePtr->operator()( projectPoint( Sites[Sites.size()-2], dim-1 ) ),
247  Sites[Sites.size()-1], myWeightImagePtr->operator()( projectPoint( Sites[Sites.size()-1], dim-1 ) ),
248  psite, myWeightImagePtr->operator()( boundedPSite ),
249  startingPoint, endPoint, dim) ))
250  {
251  Sites.pop_back();
252  boundedSites.pop_back();
253  }
254 
255  Sites.push_back( psite );
256  boundedSites.push_back( boundedPSite );
257  }
258  }
259  }
260 
261 
262  // Pruning the remaining list of sites in the periodic case.
263  if ( isPeriodic(dim) )
264  {
265  auto point = startPoint;
266  point[dim] = myLowerBoundCopy[dim];
267  for ( ; point[dim] <= endPoint[dim] - extent + 1; ++point[dim] ) // +1 in order to add the break-index site at the cycle's end.
268  {
269  Point psite = myImagePtr->operator()(point);
270 
271  if ( psite != myInfinity )
272  {
273  const Point boundedPSite = projectPoint( psite, dim-1 );
274 
275  // Site coordinates must be between startPoint and endPoint.
276  psite[dim] += extent;
277 
278  while (( Sites.size() >= 2 ) &&
279  ( myMetricPtr->hiddenByPower(Sites[Sites.size()-2], myWeightImagePtr->operator()(boundedSites[Sites.size()-2]),
280  Sites[Sites.size()-1], myWeightImagePtr->operator()(boundedSites[Sites.size()-1]),
281  psite, myWeightImagePtr->operator()(boundedPSite),
282  startingPoint, endPoint, dim) ))
283  {
284  Sites.pop_back();
285  boundedSites.pop_back();
286  }
287 
288  Sites.push_back( psite );
289  boundedSites.push_back( boundedPSite );
290  }
291  }
292  }
293  }
294 
295  // No sites found
296  if ( Sites.size() == 0 )
297  return;
298 
299  // Rewriting for both periodic and non-periodic cases.
300  std::size_t siteId = 0;
301  auto point = startPoint;
302 
303  for ( ; point[dim] <= myUpperBoundCopy[dim] ; ++point[dim] )
304  {
305  while ( ( siteId < Sites.size()-1 ) &&
306  ( myMetricPtr->closestPower(point,
307  Sites[siteId], myWeightImagePtr->operator()(boundedSites[siteId]),
308  Sites[siteId+1], myWeightImagePtr->operator()(boundedSites[siteId+1]))
309  != DGtal::ClosestFIRST ))
310  siteId++;
311 
312  myImagePtr->setValue(point, Sites[siteId]);
313  }
314 
315  // Continuing rewriting in the periodic case.
316  if ( isPeriodic(dim) )
317  {
318  for ( ; point[dim] <= endPoint[dim] ; ++point[dim] )
319  {
320  while ( ( siteId < Sites.size()-1 ) &&
321  ( myMetricPtr->closestPower(point,
322  Sites[siteId], myWeightImagePtr->operator()(boundedSites[siteId]),
323  Sites[siteId+1],myWeightImagePtr->operator()(boundedSites[siteId+1]))
324  != DGtal::ClosestFIRST ))
325  siteId++;
326 
327  myImagePtr->setValue(point - Point::base(dim, extent), Sites[siteId] - Point::base(dim, extent) );
328  }
329  }
330 
331 }
332 
333 
334 /**
335  * Constructor.
336  */
337 template <typename W,typename TSep,typename Im>
338 inline
339 DGtal::PowerMap<W,TSep,Im>::PowerMap( ConstAlias<Domain> aDomain,
340  ConstAlias<WeightImage> aWeightImage,
341  ConstAlias<PowerSeparableMetric> aMetric)
342  : myDomainPtr(&aDomain)
343  , myDomainExtent( aDomain->upperBound() - aDomain->lowerBound() + Point::diagonal(1) )
344  , myMetricPtr(&aMetric)
345  , myWeightImagePtr(&aWeightImage)
346 {
347  myPeriodicitySpec.fill( false );
348  myImagePtr = CountedPtr<OutputImage>(new OutputImage(aDomain));
349  compute();
350 }
351 
352 template <typename W,typename TSep,typename Im>
353 inline
354 DGtal::PowerMap<W,TSep,Im>::PowerMap( ConstAlias<Domain> aDomain,
355  ConstAlias<WeightImage> aWeightImage,
356  ConstAlias<PowerSeparableMetric> aMetric,
357  PeriodicitySpec const & aPeriodicitySpec )
358  : myDomainPtr(&aDomain)
359  , myDomainExtent( aDomain->upperBound() - aDomain->lowerBound() + Point::diagonal(1) )
360  , myMetricPtr(&aMetric)
361  , myWeightImagePtr(&aWeightImage)
362  , myPeriodicitySpec(aPeriodicitySpec)
363 {
364  // Finding periodic dimension index.
365  for ( Dimension i = 0; i < Space::dimension; ++i )
366  if ( isPeriodic(i) )
367  myPeriodicityIndex.push_back( i );
368 
369  myImagePtr = CountedPtr<OutputImage>(new OutputImage(aDomain));
370  compute();
371 }
372 
373 template <typename W,typename TSep,typename Im>
374 inline
375 typename DGtal::PowerMap<W,TSep,Im>::Point
376 DGtal::PowerMap<W,TSep,Im>::projectPoint( Point aPoint, const Dimension aMaxDim ) const
377 {
378  for ( std::size_t i = 0; i < myPeriodicityIndex.size() && myPeriodicityIndex[i] <= aMaxDim; ++i )
379  aPoint[ myPeriodicityIndex[i] ] = projectCoordinate( aPoint[ myPeriodicityIndex[i] ], myPeriodicityIndex[i] );
380 
381  return aPoint;
382 }
383 
384 template <typename W,typename TSep,typename Im>
385 inline
386 typename DGtal::PowerMap<W,TSep,Im>::Point
387 DGtal::PowerMap<W,TSep,Im>::projectPoint( Point aPoint ) const
388 {
389  for ( auto const & dim : myPeriodicityIndex )
390  aPoint[ dim ] = projectCoordinate( aPoint[ dim ], dim );
391 
392  return aPoint;
393 }
394 
395 template <typename W,typename TSep,typename Im>
396 inline
397 typename DGtal::PowerMap<W,TSep,Im>::Point::Coordinate
398 DGtal::PowerMap<W,TSep,Im>::projectCoordinate( typename Point::Coordinate aCoordinate, const Dimension aDim ) const
399 {
400  ASSERT( aCoordinate - myDomainPtr->lowerBound()[aDim] + myDomainExtent[aDim] >= 0 );
401  return ( aCoordinate - myDomainPtr->lowerBound()[aDim] + myDomainExtent[aDim] ) % myDomainExtent[aDim] + myDomainPtr->lowerBound()[aDim];
402 }
403 
404 template <typename W,typename TSep,typename Im>
405 inline
406 void
407 DGtal::PowerMap<W,TSep,Im>::selfDisplay ( std::ostream & out ) const
408 {
409  out << "[PowerMap] power separable metric=" << *myMetricPtr ;
410 }
411 
412 
413 // // //
414 // ///////////////////////////////////////////////////////////////////////////////
415 
416 template <typename W,typename TSep,typename Im>
417 inline
418 std::ostream&
419 DGtal::operator<< ( std::ostream & out,
420  const PowerMap<W,TSep,Im> & object )
421 {
422  object.selfDisplay( out );
423  return out;
424 }
425 
426 // // //
427 // ///////////////////////////////////////////////////////////////////////////////