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.
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.
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/>.
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
24 * Implementation of inline methods defined in VoronoiMap.h
26 * This file is part of the DGtal library.
30 //////////////////////////////////////////////////////////////////////////////
34 #include <boost/lexical_cast.hpp>
37 #include "DGtal/kernel/NumberTraits.h"
39 //////////////////////////////////////////////////////////////////////////////
41 ///////////////////////////////////////////////////////////////////////////////
42 // IMPLEMENTATION of inline methods.
43 ///////////////////////////////////////////////////////////////////////////////
45 ///////////////////////////////////////////////////////////////////////////////
46 // ----------------------- Standard services ------------------------------
48 template <typename S, typename P, typename TSep, typename TImage>
51 DGtal::VoronoiMap<S,P, TSep, TImage>::compute( )
53 //We copy the image extent
54 myLowerBoundCopy = myDomainPtr->lowerBound();
55 myUpperBoundCopy = myDomainPtr->upperBound();
57 //Point outside the domain
58 for ( auto & coord : myInfinity )
59 coord = DGtal::NumberTraits< typename Point::Coordinate >::max();
62 for ( auto const & pt : *myDomainPtr )
63 if ( (*myPointPredicatePtr)( pt ))
64 myImagePtr->setValue ( pt, myInfinity );
66 myImagePtr->setValue ( pt, pt );
68 //We process the remaining dimensions
69 for ( Dimension dim = 0; dim< S::dimension ; dim++ )
70 computeOtherSteps ( dim );
73 template <typename S, typename P,typename TSep, typename TImage>
76 DGtal::VoronoiMap<S,P, TSep, TImage>::computeOtherSteps ( const Dimension dim ) const
79 std::string title = "VoronoiMap dimension " + boost::lexical_cast<std::string>( dim ) ;
80 trace.beginBlock ( title );
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 );
92 Domain localDomain(myLowerBoundCopy, myUpperBoundCopy);
96 std::vector<Point> subRangePoints;
97 //Starting point precomputation
98 for ( auto const & pt : localDomain.subRange( subdomain ) )
99 subRangePoints.push_back( pt );
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);
107 //We solve the 1D problems sequentially
108 for ( auto const & pt : localDomain.subRange( subdomain ) )
109 computeOtherStep1D ( pt, dim);
117 // //////////////////////////////////////////////////////////////////////:
118 // ////////////////////////// Other Phases
119 template <typename S,typename P, typename TSep, typename TImage>
121 DGtal::VoronoiMap<S,P,TSep, TImage>::computeOtherStep1D ( const Point &startingPoint,
122 const Dimension dim) const
124 ASSERT(dim < S::dimension);
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];
132 // Extent along current dimension.
133 const auto extent = myUpperBoundCopy[dim] - myLowerBoundCopy[dim] + 1;
136 std::vector<Point> Sites;
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 ) );
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.
147 // For dim = 0, no sites are hidden.
148 for ( auto point = startPoint ; point[dim] <= myUpperBoundCopy[dim] ; ++point[dim] )
150 const Point psite = myImagePtr->operator()( point );
151 if ( psite != myInfinity )
152 Sites.push_back( psite );
155 // If no sites are found, then there is nothing to do.
156 if ( Sites.size() == 0 )
159 // In the periodic case and along the first dimension, the break index
160 // is at the first site found.
161 if ( isPeriodic(dim) )
163 startPoint[dim] = Sites[0][dim];
164 endPoint[dim] = startPoint[dim] + extent - 1;
166 // The first site is also the last site (with appropriate shift).
167 Sites.push_back( Sites[0] + Point::base(dim, extent) );
172 // In the periodic case, the cycle depends on break index
173 if ( isPeriodic(dim) )
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();
178 for ( auto point = startPoint; point[dim] <= myUpperBoundCopy[dim]; ++point[dim] )
180 const Point psite = myImagePtr->operator()( point );
182 if ( psite != myInfinity )
184 const auto rawDist = myMetricPtr->rawDistance( point, psite );
185 if ( rawDist < minRawDist )
187 minRawDist = rawDist;
188 startPoint[dim] = point[dim];
193 // If no sites are found, then there is nothing to do.
194 if ( minRawDist == DGtal::NumberTraits< typename SeparableMetric::RawValue >::max() )
197 endPoint[dim] = startPoint[dim] + extent - 1;
200 // Pruning the list of sites for both periodic and non-periodic cases.
201 for( auto point = startPoint ; point[dim] <= myUpperBoundCopy[dim] ; ++point[dim] )
203 const Point psite = myImagePtr->operator()(point);
205 if ( psite != myInfinity )
207 while (( Sites.size() >= 2 ) &&
208 ( myMetricPtr->hiddenBy(Sites[Sites.size()-2], Sites[Sites.size()-1] ,
209 psite, startingPoint, endPoint, dim) ))
212 Sites.push_back( psite );
216 // Pruning the remaining list of sites in the periodic case.
217 if ( isPeriodic(dim) )
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.
223 Point psite = myImagePtr->operator()(point);
225 if ( psite != myInfinity )
227 // Site coordinates must be between startPoint and endPoint.
228 psite[dim] += extent;
230 while (( Sites.size() >= 2 ) &&
231 ( myMetricPtr->hiddenBy(Sites[Sites.size()-2], Sites[Sites.size()-1] ,
232 psite, startingPoint, endPoint, dim) ))
235 Sites.push_back( psite );
242 if ( Sites.size() == 0 )
245 // Rewriting for both periodic and non-periodic cases.
246 std::size_t siteId = 0;
247 auto point = startPoint;
249 for ( ; point[dim] <= myUpperBoundCopy[dim] ; ++point[dim] )
251 while ( ( siteId < Sites.size()-1 ) &&
252 ( myMetricPtr->closest(point, Sites[siteId], Sites[siteId+1])
253 != DGtal::ClosestFIRST ))
256 myImagePtr->setValue(point, Sites[siteId]);
259 // Continuing rewriting in the periodic case.
260 if ( isPeriodic(dim) )
262 for ( ; point[dim] <= endPoint[dim] ; ++point[dim] )
264 while ( ( siteId < Sites.size()-1 ) &&
265 ( myMetricPtr->closest(point, Sites[siteId], Sites[siteId+1])
266 != DGtal::ClosestFIRST ))
269 myImagePtr->setValue(point - Point::base(dim, extent), Sites[siteId] - Point::base(dim, extent) );
279 template <typename S,typename P,typename TSep, typename TImage>
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)
289 myPeriodicitySpec.fill( false );
290 myImagePtr = CountedPtr<OutputImage>( new OutputImage(aDomain) );
294 template <typename S,typename P,typename TSep, typename TImage>
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)
306 // Finding periodic dimension index.
307 for ( Dimension i = 0; i < Space::dimension; ++i )
309 myPeriodicityIndex.push_back( i );
311 myImagePtr = CountedPtr<OutputImage>( new OutputImage(aDomain) );
315 template <typename S,typename P,typename TSep, typename TImage>
317 typename DGtal::VoronoiMap<S, P, TSep, TImage>::Point
318 DGtal::VoronoiMap<S, P, TSep, TImage>::projectPoint( Point aPoint ) const
320 for ( auto const & dim : myPeriodicityIndex )
321 aPoint[ dim ] = projectCoordinate( aPoint[ dim ], dim );
326 template <typename S,typename P,typename TSep, typename TImage>
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
331 ASSERT( aCoordinate - myDomainPtr->lowerBound()[aDim] + myDomainExtent[aDim] >= 0 );
332 return ( aCoordinate - myDomainPtr->lowerBound()[aDim] + myDomainExtent[aDim] ) % myDomainExtent[aDim] + myDomainPtr->lowerBound()[aDim];
335 template <typename S,typename P,typename TSep, typename TImage>
338 DGtal::VoronoiMap<S,P, TSep, TImage>::selfDisplay ( std::ostream & out ) const
340 out << "[VoronoiMap] separable metric=" << *myMetricPtr ;
345 // ///////////////////////////////////////////////////////////////////////////////
347 template <typename S,typename P,typename TSep, typename TImage>
350 DGtal::operator<< ( std::ostream & out,
351 const VoronoiMap<S,P,TSep, TImage> & object )
353 object.selfDisplay( out );