DGtal 1.3.0
Loading...
Searching...
No Matches
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//----------------------------------------------------------------------------
37template<typename DSL, typename OutputIterator>
38inline
39typename DSL::Vector
40DGtal::functions::
41smartCH(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//----------------------------------------------------------------------------
71template<typename PointVector, typename Coordinate, typename Position,
72 typename PositionFunctor, typename OutputIterator>
73inline
74PointVector
75DGtal::functions::
76smartCH(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//----------------------------------------------------------------------------
148template<typename Position, typename Coordinate, typename PointVector,
149 typename OutputIterator,
150 typename PositionFunctor,
151 typename TruncationFunctor1, typename TruncationFunctor2>
152inline
153bool
154DGtal::functions::
155smartCHNextVertex(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//----------------------------------------------------------------------------
220template<typename DSS, typename OutputIterator>
221inline
222typename DSS::Vector
223DGtal::functions::
224reversedSmartCH(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//----------------------------------------------------------------------------
245template<typename PointVector, typename Position,
246 typename PositionFunctor, typename OutputIterator>
247inline
248PointVector
249DGtal::functions::
250reversedSmartCH(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//----------------------------------------------------------------------------
301template<typename PointVector, typename Position,
302 typename OutputIterator,
303 typename TruncationFunctor1, typename TruncationFunctor2,
304 typename PositionFunctor>
305inline
306bool
307DGtal::functions::
308smartCHPreviousVertex(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}
Generated on Sun Nov 27 2022 15:15:35 for DGtal by doxygen 1.9.5