DGtal  1.5.beta
ArithmeticalDSSConvexHull.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 ArithmeticalDSSConvexHull.ih
19  * @author Tristan Roussillon (\c tristan.roussillon@liris.cnrs.fr )
20  * Laboratoire d'InfoRmatique en Image et Systèmes d'information - LIRIS (CNRS, UMR 5205), CNRS, France
21  *
22  * @date 2013/10/28
23  *
24  * Implementation of inline methods defined in ArithmeticalDSSConvexHull.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 
33 #include "DGtal/geometry/curves/ArithmeticalDSSConvexHull.h"
34 //////////////////////////////////////////////////////////////////////////////
35 
36 //----------------------------------------------------------------------------
37 template<typename DSL, typename OutputIterator>
38 inline
39 typename DSL::Vector
40 DGtal::functions::
41 smartCH(const DSL& aDSL,
42  const typename DSL::Point& aFirstPoint,
43  const typename DSL::Position& aLength,
44  OutputIterator uIto, OutputIterator lIto)
45 {
46  typedef typename DSL::Vector Vector;
47  typedef typename DSL::Position Position;
48  typedef typename DSL::Coordinate Coordinate;
49  ASSERT( aDSL.isValid() );
50  ASSERT( aLength >= NumberTraits<Position>::ONE );
51 
52  //position functor
53  functors::PositionFunctorFrom2DPoint<Vector,Position> position(aDSL.shift());
54 
55  //steps
56  Vector step = aDSL.steps().first;
57  Coordinate rStep = DSL::toCoordinate( aDSL.remainder(step) );
58  Vector shift = -aDSL.shift();
59  Coordinate rShift = DSL::toCoordinate( aDSL.remainder(shift) );
60 
61  //intercept
62  Coordinate intercept = DSL::toCoordinate( aDSL.mu() - aDSL.remainder(aFirstPoint) );
63 
64  Position lastPosition = position(aFirstPoint) + aLength;
65 
66  return smartCH(aFirstPoint, intercept, lastPosition, step, rStep, shift, rShift, position, uIto, lIto);
67 }
68 
69 
70 //----------------------------------------------------------------------------
71 template<typename PointVector, typename Coordinate, typename Position,
72  typename PositionFunctor, typename OutputIterator>
73 inline
74 PointVector
75 DGtal::functions::
76 smartCH(const PointVector& aFirstPoint,
77  const Coordinate& aRemainderBound,
78  const Position& aPositionBound,
79  const PointVector& aStep,
80  const Coordinate& aRStep,
81  const PointVector& aShift,
82  const Coordinate& aRShift,
83  const PositionFunctor& aPositionFunctor,
84  OutputIterator uIto, OutputIterator lIto)
85 {
86  BOOST_STATIC_ASSERT(( concepts::ConceptUtils::SameType<Coordinate, typename PointVector::Coordinate>::value ));
87  BOOST_STATIC_ASSERT(( concepts::ConceptUtils::SameType<Position, typename PositionFunctor::Position>::value ));
88  BOOST_CONCEPT_ASSERT(( boost_concepts::IncrementableIteratorConcept<OutputIterator> ));
89  BOOST_CONCEPT_ASSERT(( boost_concepts::WritableIteratorConcept<OutputIterator,PointVector> ));
90 
91  //----------------- init --------------------//
92  //functors
93  functors::StrictTruncationFunctor<Coordinate> sTruncation;
94  functors::LargeTruncationFunctor<Coordinate> lTruncation;
95 
96  PointVector L, U, V; //point/vectors
97  Coordinate rL, rU, rV; //remainders of the above point/vectors
98  Coordinate q; //quotient of the divisions
99 
100  U = aFirstPoint;
101  rU = NumberTraits<Coordinate>::ZERO;
102  *uIto++ = U;
103 
104  L = U + aShift;
105  rL = rU + aRShift;
106  *lIto++ = L;
107 
108  q = (aRemainderBound - (rU + aRStep)) / aRShift;
109  V = aStep + aShift * q;
110  rV = aRStep + q * aRShift;
111 
112  //----------------- main --------------------//
113  bool stop = false;
114  while ( (rV != 0) && (!stop) )
115  {
116 
117  if (rV > 0)
118  {
119 
120  stop = smartCHNextVertex(aPositionBound,
121  aRemainderBound,
122  L, rL, U, rU, V, rV,
123  lIto, aPositionFunctor,
124  sTruncation, lTruncation);
125  }
126  else
127  { //if (rV < 0)
128 
129  stop = smartCHNextVertex(aPositionBound,
130  aRemainderBound,
131  U, rU, L, rL, V, rV,
132  uIto, aPositionFunctor,
133  lTruncation, sTruncation);
134  }
135 
136  ASSERT( (V[0]*(L-U)[1] - V[1]*(L-U)[0]) == 1 ); //eq. 6
137  ASSERT( (stop) || ((!stop)&&(rL + rV < aRemainderBound)) ); //eq. 5
138  ASSERT( (stop) || ((!stop)&&(rU + rV >= aRemainderBound)) ); //eq. 5
139  ASSERT( (stop) || ((!stop)&&((aPositionFunctor(L)-aPositionFunctor(aFirstPoint)) < aPositionFunctor(V))) ); //eq. 10
140  ASSERT( (stop) || ((!stop)&&((aPositionFunctor(U)-aPositionFunctor(aFirstPoint)) < aPositionFunctor(V))) ); //eq. 10
141  }
142 
143  return V;
144 }
145 
146 
147 //----------------------------------------------------------------------------
148 template<typename Position, typename Coordinate, typename PointVector,
149  typename OutputIterator,
150  typename PositionFunctor,
151  typename TruncationFunctor1, typename TruncationFunctor2>
152 inline
153 bool
154 DGtal::functions::
155 smartCHNextVertex(const Position& positionBound,
156  const Coordinate& remainderBound,
157  PointVector& X,
158  Coordinate& rX,
159  const PointVector& Y,
160  const Coordinate& rY,
161  PointVector& V,
162  Coordinate& rV,
163  OutputIterator ito,
164  const PositionFunctor& pos,
165  const TruncationFunctor1& f1,
166  const TruncationFunctor2& f2)
167 {
168  ASSERT( rV != NumberTraits<Coordinate>::ZERO );
169  ASSERT( ((V[0]*(X-Y)[1] - V[1]*(X-Y)[0]) == 1)||
170  ((V[0]*(X-Y)[1] - V[1]*(X-Y)[0]) == -1) );
171  BOOST_CONCEPT_ASSERT(( boost_concepts::IncrementableIteratorConcept<OutputIterator> ));
172  BOOST_CONCEPT_ASSERT(( boost_concepts::WritableIteratorConcept<OutputIterator,PointVector> ));
173 
174  bool stop = false;
175  Coordinate q = f1(remainderBound - rX, rV); //first ray casting
176  X += V * q;
177  if ( pos(X) <= positionBound )
178  {
179  rX += q * rV;
180  *ito++ = X;
181 
182  q = f2( remainderBound - (rY + rV), (rX - rY) ); //second ray casting
183  V += (X - Y) * q;
184  if ( (pos(V) + pos(Y)) <= positionBound )
185  {
186  rV += q * (rX - rY);
187  }
188  else
189  {
190  stop = true;
191  PointVector XY = (X - Y);
192  V -= (XY) * q;
193  q = ( (positionBound - (pos(Y) + pos(V))) / pos(XY) );
194  if (q > 0)
195  V += XY * q;
196  }
197  }
198  else
199  {
200  stop = true;
201  X -= V * q;
202  q = ( (positionBound - pos(X)) / pos(V) );
203  ASSERT(q >= 0);
204  if (q > 0)
205  {
206  X += V * q;
207  *ito++ = X;
208  }
209  PointVector XY = (X - Y);
210  q = ( (positionBound - (pos(Y) + pos(V))) / pos(XY) );
211  if (q > 0)
212  V += XY * q;
213  }
214 
215  return stop;
216 }
217 
218 
219 //----------------------------------------------------------------------------
220 template<typename DSS, typename OutputIterator>
221 inline
222 typename DSS::Vector
223 DGtal::functions::
224 reversedSmartCH(const DSS& aDSS,
225  const typename DSS::Position& aPositionBound,
226  OutputIterator uIto, OutputIterator lIto)
227 {
228  ASSERT( aDSS.position(aDSS.back()) <= aPositionBound );
229 
230  typedef typename DSS::Vector Vector;
231  typedef typename DSS::Position Position;
232 
233  //position functor
234  functors::PositionFunctorFrom2DPoint<Vector,Position> position(aDSS.dsl().shift());
235  //direction vector
236  Vector v(aDSS.b(), aDSS.a());
237 
238  return reversedSmartCH( aDSS.Uf(), aDSS.Lf()-aDSS.dsl().shift(), v,
239  aDSS.position(aDSS.back()), aPositionBound,
240  position, uIto, lIto );
241 }
242 
243 
244 //----------------------------------------------------------------------------
245 template<typename PointVector, typename Position,
246  typename PositionFunctor, typename OutputIterator>
247 inline
248 PointVector
249 DGtal::functions::
250 reversedSmartCH(PointVector U, PointVector L, PointVector V,
251  const Position& aFirstPosition, const Position& aLastPosition,
252  const PositionFunctor& aPositionFunctor,
253  OutputIterator uIto, OutputIterator lIto)
254 {
255  BOOST_STATIC_ASSERT(( concepts::ConceptUtils::SameType<Position, typename PositionFunctor::Position>::value ));
256  BOOST_CONCEPT_ASSERT(( boost_concepts::IncrementableIteratorConcept<OutputIterator> ));
257  BOOST_CONCEPT_ASSERT(( boost_concepts::WritableIteratorConcept<OutputIterator,PointVector> ));
258  typedef typename PointVector::Coordinate Coordinate;
259  BOOST_STATIC_ASSERT(( concepts::ConceptUtils::SameType<Coordinate, typename PointVector::Coordinate>::value ));
260  ASSERT( aFirstPosition <= aLastPosition );
261  ASSERT( (V[0]*(L-U)[1] - V[1]*(L-U)[0]) == NumberTraits<Coordinate>::ONE ); //eq. 6
262 
263  //functors
264  functors::StrictTruncationFunctor<Coordinate> sTruncation;
265  functors::LargeTruncationFunctor<Coordinate> lTruncation;
266 
267  //----------------- init --------------------//
268  *lIto++ = L;
269  *uIto++ = U;
270  bool stop = false;
271  if (aPositionFunctor(L) >= aPositionFunctor(U))
272  {
273  if ((aPositionFunctor(U)+aPositionFunctor(V)) <= aLastPosition)
274  stop = true;
275  }
276  else
277  {
278  if ((aPositionFunctor(L)+aPositionFunctor(V)) <= aLastPosition)
279  stop = true;
280  }
281 
282  //----------------- main --------------------//
283  while ( (aPositionFunctor(L) != aPositionFunctor(U)) && (!stop) )
284  {
285  if (aPositionFunctor(L) > aPositionFunctor(U))
286  {
287  stop = smartCHPreviousVertex(L, U, V, aFirstPosition, aLastPosition,
288  lIto, aPositionFunctor, sTruncation, lTruncation);
289  }
290  else
291  {
292  stop = smartCHPreviousVertex(U, L, V, aFirstPosition, aLastPosition,
293  uIto, aPositionFunctor, sTruncation, lTruncation);
294  }
295  }
296 
297  return V;
298 }
299 
300 //----------------------------------------------------------------------------
301 template<typename PointVector, typename Position,
302  typename OutputIterator,
303  typename TruncationFunctor1, typename TruncationFunctor2,
304  typename PositionFunctor>
305 inline
306 bool
307 DGtal::functions::
308 smartCHPreviousVertex(PointVector& X, const PointVector& Y, PointVector& V,
309  const Position& aFirstPosition, const Position& aLastPosition,
310  OutputIterator ito,
311  const PositionFunctor& pos,
312  const TruncationFunctor1& f1,
313  const TruncationFunctor2& f2)
314 {
315  BOOST_CONCEPT_ASSERT(( boost_concepts::IncrementableIteratorConcept<OutputIterator> ));
316  BOOST_CONCEPT_ASSERT(( boost_concepts::WritableIteratorConcept<OutputIterator,PointVector> ));
317  typedef PointVector Vector;
318 
319  //pre-conditions
320  ASSERT(pos(X) > pos(Y));
321  ASSERT((pos(X)-aFirstPosition) < pos(V)); //eq 10
322  ASSERT((pos(Y)-aFirstPosition) < pos(V)); //eq 10
323 
324  bool stop = false;
325 
326  //update V
327  Vector XY = X-Y;
328  Position q = f1( (aFirstPosition+pos(V)-pos(Y)), pos(XY) );
329  ASSERT(q > 0);
330  V -= q * XY;
331  ASSERT(pos(V) <= (pos(X)-aFirstPosition)); //eq 9 (NB: equality at position 0)
332  if (pos(Y) + pos(V) > aLastPosition)
333  {
334  //update X
335  q = f2((pos(X)-aFirstPosition), pos(V));
336  ASSERT(q >= 0);
337  X -= q * V;
338  ASSERT((pos(X)-aFirstPosition) < pos(V)); //eq 10
339 
340  if (pos(X) + pos(V) <= aLastPosition)
341  { //X can be increased
342  stop = true;
343  q = (aLastPosition - (pos(X) + pos(V))) / pos(V);
344  X += q * V;
345  }
346  *ito++ = X;
347  }
348  else
349  {
350  stop = true;
351  //V can be increased
352  q = (aLastPosition - (pos(Y) + pos(V)) ) / pos(XY);
353  ASSERT(q >= 0);
354  V += q * XY;
355  //update X
356  q = f2( (pos(X) - pos(Y)), pos(V) );
357  if (q >= 0)
358  {
359  X -= q * V;
360  *ito++ = X;
361  }
362  //there is nothing else to do
363  ASSERT(pos(X) <= aLastPosition);
364  ASSERT(pos(Y) <= aLastPosition);
365  ASSERT(pos(Y) + pos(V) <= aLastPosition);
366  }
367 
368  return stop;
369 }