DGtal  0.9.2
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 #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 W, typename Sep, typename Im>
50 inline
51 DGtal::PowerMap< W,Sep,Im>::~PowerMap()
52 {
53 }
54 
55 
56 template < typename W, typename Sep, typename Im>
57 inline
58 typename DGtal::PowerMap<W, Sep,Im>::Self &
59 DGtal::PowerMap<W, Sep,Im>::operator=(const Self &aOtherPowerMap )
60 {
61  if (this != &aOtherPowerMap)
62  {
63  myMetricPtr = aOtherPowerMap.myMetricPtr;
64  myImagePtr = aOtherPowerMap.myImagePtr;
65  myWeightImagePtr = aOtherPowerMap.myWeightImagePtr;
66  myDomainPtr = aOtherPowerMap.myDomainPtr;
67  myInfinity = aOtherPowerMap.myInfinity;
68  myLowerBoundCopy = aOtherPowerMap.myLowerBoundCopy;
69  myUpperBoundCopy = aOtherPowerMap.myUpperBoundCopy;
70  }
71  return *this;
72 }
73 
74 template < typename W, typename Sep, typename Im>
75 inline
76 void
77 DGtal::PowerMap<W, Sep,Im>::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 the map: the power map at point p is:
87  // - p if p is an input weighted point (with weight > 0);
88  // - myInfinity otherwise.
89  for(typename Domain::ConstIterator it = myDomainPtr->begin(),
90  itend = myDomainPtr->end();
91  it != itend;
92  ++it)
93  if (myWeightImagePtr->domain().isInside( *it ) &&
94  (myWeightImagePtr->operator()( *it ) >0 ))
95  myImagePtr->setValue ( *it, *it );
96  else
97  myImagePtr->setValue ( *it, myInfinity );
98 
99  //We process the dimensions one by one
100  for ( Dimension dim = 0; dim < W::Domain::Space::dimension ; dim++ )
101  computeOtherSteps ( dim );
102 }
103 
104 template < typename W, typename Sep, typename Im>
105 inline
106 void
107 DGtal::PowerMap<W, Sep,Im>::computeOtherSteps ( const Dimension dim ) const
108 {
109 #ifdef VERBOSE
110  std::string title = "Powermap dimension " + boost::lexical_cast<std::string>( dim ) ;
111  trace.beginBlock ( title );
112 #endif
113 
114  typedef typename Domain::ConstSubRange::ConstIterator ConstDomIt;
115 
116  //We setup the subdomain iterator
117  //the iterator will scan dimension using the order:
118  // {n-1, n-2, ... 1} (we skip the '0' dimension).
119  std::vector<Dimension> subdomain;
120  subdomain.reserve(W::Domain::Space::dimension - 1);
121  for (unsigned int k = 0; k < W::Domain::Space::dimension ; k++)
122  if ( ((int)W::Domain::Space::dimension - 1 - k) != dim)
123  subdomain.push_back( (int)W::Domain::Space::dimension - 1 - k );
124 
125  Domain localDomain(myLowerBoundCopy, myUpperBoundCopy);
126 
127 
128 #ifdef WITH_OPENMP
129  //Parallel loop
130  std::vector<Point> subRangePoints;
131  //Starting point precomputation
132  for (ConstDomIt it = localDomain.subRange( subdomain ).begin(),
133  itend = localDomain.subRange( subdomain ).end();
134  it != itend; ++it)
135  subRangePoints.push_back( *it );
136  //We run the 1D problems in //
137 #pragma omp parallel for schedule(dynamic)
138  for (size_t i = 0; i < subRangePoints.size(); ++i)
139  computeOtherStep1D ( subRangePoints[i], dim);
140 
141 #else
142  //We solve the 1D problems sequentially
143  for (ConstDomIt it = localDomain.subRange( subdomain ).begin(),
144  itend = localDomain.subRange( subdomain ).end();
145  it != itend; ++it)
146  computeOtherStep1D ( *it, dim);
147 #endif
148 
149 #ifdef VERBOSE
150  trace.endBlock();
151 #endif
152 }
153 
154 // //////////////////////////////////////////////////////////////////////:
155 // ////////////////////////// Other Phases
156 template <typename W, typename Sep, typename Im>
157 void
158 DGtal::PowerMap<W,Sep,Im>::computeOtherStep1D ( const Point &startingPoint,
159  const Dimension dim) const
160 {
161  Point point = startingPoint;
162  Point endpoint = startingPoint;
163  Point psite;
164  int nbSites = -1;
165  std::vector<Point> Sites;
166 
167  //Reserve
168  Sites.reserve( myUpperBoundCopy[dim] - myLowerBoundCopy[dim] +1);
169 
170  //endpoint of the 1D row
171  endpoint[dim] = myUpperBoundCopy[dim];
172 
173  //Pruning the list of sites (dim=0 implies no hidden sites)
174  if (dim==0)
175  {
176  for(Abscissa i = myLowerBoundCopy[dim] ; i <= myUpperBoundCopy[dim] ; i++)
177  {
178  psite = myImagePtr->operator()( point );
179  if ( psite != myInfinity )
180  {
181  nbSites++;
182  Sites.push_back( psite );
183  }
184  point[dim] ++;
185  }
186  }
187  else
188  {
189  //Pruning the list of sites
190  for(Abscissa i = myLowerBoundCopy[dim] ; i <= myUpperBoundCopy[dim] ; i++)
191  {
192  psite = myImagePtr->operator()(point);
193  if ( psite != myInfinity )
194  {
195  while ((nbSites >= 1) &&
196  ( myMetricPtr->hiddenByPower(Sites[nbSites-1], myWeightImagePtr->operator()(Sites[nbSites-1]),
197  Sites[nbSites] , myWeightImagePtr->operator()(Sites[nbSites]),
198  psite, myWeightImagePtr->operator()(psite),
199  startingPoint, endpoint, dim) ))
200  {
201  nbSites --;
202  Sites.pop_back();
203  }
204  nbSites++;
205  Sites.push_back( psite );
206  }
207  point[dim] ++;
208  }
209  }
210 
211  //No sites found
212  if (nbSites == -1)
213  return;
214 
215  int k = 0;
216 
217  //Rewriting
218  point[dim] = myLowerBoundCopy[dim];
219  for(Abscissa i = myLowerBoundCopy[dim] ; i <= myUpperBoundCopy[dim] ; i++)
220  {
221  while ( (k < nbSites) &&
222  ( myMetricPtr->closestPower(point,
223  Sites[k], myWeightImagePtr->operator()(Sites[k]),
224  Sites[k+1],myWeightImagePtr->operator()(Sites[k+1]))
225  != DGtal::ClosestFIRST ))
226  k++;
227  myImagePtr->setValue(point, Sites[k]);
228  point[dim]++;
229  }
230 }
231 
232 
233 /**
234  * Constructor.
235  */
236 template <typename W,typename TSep,typename Im>
237 inline
238 DGtal::PowerMap<W,TSep,Im>::PowerMap( ConstAlias<Domain> aDomain,
239  ConstAlias<WeightImage> aWeightImage,
240  ConstAlias<PowerSeparableMetric> aMetric):
241  myDomainPtr(&aDomain),
242  myMetricPtr(&aMetric),
243  myWeightImagePtr(&aWeightImage)
244 
245 {
246  myImagePtr = CountedPtr<OutputImage>(new OutputImage(aDomain));
247  compute();
248 }
249 
250 template <typename W,typename TSep,typename Im>
251 inline
252 void
253 DGtal::PowerMap<W,TSep,Im>::selfDisplay ( std::ostream & out ) const
254 {
255  out << "[PowerMap] power separable metric=" << *myMetricPtr ;
256 }
257 
258 
259 // // //
260 // ///////////////////////////////////////////////////////////////////////////////
261 
262 template <typename W,typename TSep,typename Im>
263 inline
264 std::ostream&
265 DGtal::operator<< ( std::ostream & out,
266  const PowerMap<W,TSep,Im> & object )
267 {
268  object.selfDisplay( out );
269  return out;
270 }
271 
272 // // //
273 // ///////////////////////////////////////////////////////////////////////////////