DGtal 1.3.0
Loading...
Searching...
No Matches
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
43template <typename S, typename P, typename TSep, typename TImage>
44inline
45void
46DGtal::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
68template <typename S, typename P,typename TSep, typename TImage>
69inline
70void
71DGtal::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 (size_t i = 0; i < subRangePoints.size(); ++i)
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
114template <typename S,typename P, typename TSep, typename TImage>
115void
116DGtal::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 */
274template <typename S,typename P,typename TSep, typename TImage>
275inline
276DGtal::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
289template <typename S,typename P,typename TSep, typename TImage>
290inline
291DGtal::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
310template <typename S,typename P,typename TSep, typename TImage>
311inline
312typename DGtal::VoronoiMap<S, P, TSep, TImage>::Point
313DGtal::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
321template <typename S,typename P,typename TSep, typename TImage>
322inline
323typename DGtal::VoronoiMap<S, P, TSep, TImage>::Point::Coordinate
324DGtal::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
330template <typename S,typename P,typename TSep, typename TImage>
331inline
332void
333DGtal::VoronoiMap<S,P, TSep, TImage>::selfDisplay ( std::ostream & out ) const
334{
335 out << "[VoronoiMap] separable metric=" << *myMetricPtr ;
336}
337
338
339// // //
340// ///////////////////////////////////////////////////////////////////////////////
341
342template <typename S,typename P,typename TSep, typename TImage>
343inline
344std::ostream&
345DGtal::operator<< ( std::ostream & out,
346 const VoronoiMap<S,P,TSep, TImage> & object )
347{
348 object.selfDisplay( out );
349 return out;
350}