DGtal  0.9.2
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 #ifdef VERBOSE
33 #include <boost/lexical_cast.hpp>
34 #endif
35 //////////////////////////////////////////////////////////////////////////////
36 
37 ///////////////////////////////////////////////////////////////////////////////
38 // IMPLEMENTATION of inline methods.
39 ///////////////////////////////////////////////////////////////////////////////
40 
41 ///////////////////////////////////////////////////////////////////////////////
42 // ----------------------- Standard services ------------------------------
43 
44 
45 
46 /**
47  * Destructor.
48  */
49 template <typename S, typename P,typename TSep, typename TImage>
50 inline
51 DGtal::VoronoiMap<S, P,TSep, TImage>::~VoronoiMap()
52 {
53 }
54 
55 
56 template <typename S, typename P, typename TSep, typename TImage>
57 inline
58 typename DGtal::VoronoiMap<S,P, TSep, TImage>::Self &
59 DGtal::VoronoiMap<S,P, TSep, TImage>::operator=(const Self &aOtherVoronoiMap )
60 {
61  if (this != &aOtherVoronoiMap)
62  {
63  myMetricPtr = aOtherVoronoiMap.myMetricPtr;
64  myImagePtr = aOtherVoronoiMap.myImagePtr;
65  myPointPredicatePtr = aOtherVoronoiMap.myPointPredicatePtr;
66  myDomainPtr = aOtherVoronoiMap.myDomainPtr;
67  myInfinity = aOtherVoronoiMap.myInfinity;
68  myLowerBoundCopy = aOtherVoronoiMap.myLowerBoundCopy;
69  myUpperBoundCopy = aOtherVoronoiMap.myUpperBoundCopy;
70  }
71  return *this;
72 }
73 
74 template <typename S, typename P, typename TSep, typename TImage>
75 inline
76 void
77 DGtal::VoronoiMap<S,P, TSep, TImage>::compute( )
78 {
79  //We copy the image extent
80  myLowerBoundCopy = myDomainPtr->lowerBound();
81  myUpperBoundCopy = myDomainPtr->upperBound();
82 
83  //Point outside the domain
84  myInfinity = myDomainPtr->upperBound() + Point::diagonal(1);
85 
86  //Init
87  for(typename Domain::ConstIterator it = myDomainPtr->begin(), itend = myDomainPtr->end();
88  it != itend;
89  ++it)
90  if ( (*myPointPredicatePtr)( *it ))
91  myImagePtr->setValue ( *it, myInfinity );
92  else
93  myImagePtr->setValue ( *it, *it );
94 
95  //We process the remaining dimensions
96  for ( Dimension dim = 0; dim< S::dimension ; dim++ )
97  computeOtherSteps ( dim );
98 }
99 
100 template <typename S, typename P,typename TSep, typename TImage>
101 inline
102 void
103 DGtal::VoronoiMap<S,P, TSep, TImage>::computeOtherSteps ( const Dimension dim ) const
104 {
105 #ifdef VERBOSE
106  std::string title = "Voro dimension " + boost::lexical_cast<std::string>( dim ) ;
107  trace.beginBlock ( title );
108 #endif
109 
110  typedef typename Domain::ConstSubRange::ConstIterator ConstDomIt;
111 
112  //We setup the subdomain iterator
113  //the iterator will scan dimension using the order:
114  // {n-1, n-2, ... 1} (we skip the '0' dimension).
115  std::vector<Dimension> subdomain;
116  subdomain.reserve(S::dimension - 1);
117  for ( int k = 0; k < (int)S::dimension ; k++)
118  if ( static_cast<Dimension>(((int)S::dimension - 1 - k)) != dim)
119  subdomain.push_back( (int)S::dimension - 1 - k );
120 
121  Domain localDomain(myLowerBoundCopy, myUpperBoundCopy);
122 
123 #ifdef WITH_OPENMP
124  //Parallel loop
125  std::vector<Point> subRangePoints;
126  //Starting point precomputation
127  for (ConstDomIt it = localDomain.subRange( subdomain ).begin(),
128  itend = localDomain.subRange( subdomain ).end();
129  it != itend; ++it)
130  subRangePoints.push_back( *it );
131  //We run the 1D problems in //
132 #pragma omp parallel for schedule(dynamic)
133  for (size_t i = 0; i < subRangePoints.size(); ++i)
134  computeOtherStep1D ( subRangePoints[i], dim);
135 
136 #else
137  //We solve the 1D problems sequentially
138  for (ConstDomIt it = localDomain.subRange( subdomain ).begin(),
139  itend = localDomain.subRange( subdomain ).end();
140  it != itend; ++it)
141  computeOtherStep1D ( *it, dim);
142 #endif
143 
144 #ifdef VERBOSE
145  trace.endBlock();
146 #endif
147 }
148 
149 // //////////////////////////////////////////////////////////////////////:
150 // ////////////////////////// Other Phases
151 template <typename S,typename P, typename TSep, typename TImage>
152 void
153 DGtal::VoronoiMap<S,P,TSep, TImage>::computeOtherStep1D ( const Point &startingPoint,
154  const Dimension dim) const
155 {
156  Point point = startingPoint;
157  Point endpoint = startingPoint;
158  Point psite;
159  int nbSites = -1;
160  std::vector<Point> Sites;
161 
162  ASSERT(dim < S::dimension);
163 
164  //Reserve
165  Sites.reserve( myUpperBoundCopy[dim] - myLowerBoundCopy[dim] +1);
166 
167  //endpoint of the 1D row
168  endpoint[dim] = myUpperBoundCopy[dim];
169 
170  //Pruning the list of sites (dim=0 implies no hibben sites)
171  if (dim==0)
172  {
173  for(Abscissa i = myLowerBoundCopy[dim] ; i <= myUpperBoundCopy[dim] ; i++)
174  {
175  psite = myImagePtr->operator()(point);
176  if ( psite != myInfinity )
177  {
178  nbSites++;
179  Sites.push_back( psite );
180  }
181  point[dim] ++;
182  }
183  }
184  else
185  {
186  //Pruning the list of sites
187  for(Abscissa i = myLowerBoundCopy[dim] ; i <= myUpperBoundCopy[dim] ; i++)
188  {
189  psite = myImagePtr->operator()(point);
190  if ( psite != myInfinity )
191  {
192  while ((nbSites >= 1) &&
193  ( myMetricPtr->hiddenBy(Sites[nbSites-1], Sites[nbSites] ,
194  psite, startingPoint, endpoint, dim) ))
195  {
196  nbSites --;
197  Sites.pop_back();
198  }
199  nbSites++;
200  Sites.push_back( psite );
201  }
202  point[dim] ++;
203  }
204  }
205 
206  //No sites found
207  if (nbSites == -1)
208  return;
209 
210  int k = 0;
211 
212  //Rewriting
213  point[dim] = myLowerBoundCopy[dim];
214  for(Abscissa i = myLowerBoundCopy[dim] ; i <= myUpperBoundCopy[dim] ; i++)
215  {
216  while ( (k < nbSites) &&
217  ( myMetricPtr->closest(point, Sites[k], Sites[k+1])
218  != DGtal::ClosestFIRST ))
219  k++;
220 
221  myImagePtr->setValue(point, Sites[k]);
222  point[dim]++;
223  }
224 }
225 
226 
227 /**
228  * Constructor.
229  */
230 template <typename S,typename P,typename TSep, typename TImage>
231 inline
232 DGtal::VoronoiMap<S,P, TSep, TImage>::VoronoiMap( ConstAlias<Domain> aDomain,
233  ConstAlias<PointPredicate> aPredicate,
234  ConstAlias<SeparableMetric> aMetric):
235  myDomainPtr(&aDomain), myPointPredicatePtr(&aPredicate),
236  myMetricPtr(&aMetric)
237 {
238  myImagePtr = CountedPtr<OutputImage>( new OutputImage(aDomain) );
239  compute();
240 }
241 
242 template <typename S,typename P,typename TSep, typename TImage>
243 inline
244 void
245 DGtal::VoronoiMap<S,P, TSep, TImage>::selfDisplay ( std::ostream & out ) const
246 {
247  out << "[VoronoiMap] separable metric=" << *myMetricPtr ;
248 }
249 
250 
251 // // //
252 // ///////////////////////////////////////////////////////////////////////////////
253 
254 template <typename S,typename P,typename TSep, typename TImage>
255 inline
256 std::ostream&
257 DGtal::operator<< ( std::ostream & out,
258  const VoronoiMap<S,P,TSep, TImage> & object )
259 {
260  object.selfDisplay( out );
261  return out;
262 }