DGtal 1.3.0
Loading...
Searching...
No Matches
ArithmeticalDSSFactory.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 ArithmeticalDSSFactory.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 ArithmeticalDSSFactory.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#include "DGtal/kernel/NumberTraits.h"
35#include "DGtal/base/OneItemOutputIterator.h"
36//////////////////////////////////////////////////////////////////////////////
37
38///////////////////////////////////////////////////////////////////////////////
39// IMPLEMENTATION of inline methods.
40///////////////////////////////////////////////////////////////////////////////
41
42//-----------------------------------------------------------------------------
43template <typename TCoordinate, typename TInteger, unsigned short adjacency>
44inline
45DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
46DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::
47createSubsegment(const DSL& aDSL, const Point& aF, const Point& aL)
48{
49 ASSERT( aDSL.isInDSL( aF ) );
50 ASSERT( aDSL.isInDSL( aL ) );
51 ASSERT( aDSL.beforeOrEqual(aF, aL) );
52
53 //specific case: DSS of one point
54 if (aF == aL)
55 {
56 return DSS( aF );
57 }
58 //otherwise
59 else
60 {
61 //running smartCH to compute the minimal characteristics AND the leaning points
62 OneItemOutputIterator<Point> lastUpperVertex, lastLowerVertex;
63 using namespace functions;
64
65 Vector v = smartCH( aDSL, aF, ( aDSL.position(aL) - aDSL.position(aF) ),
66 lastUpperVertex, lastLowerVertex );
67
68 //DSS construction
69 Point U = lastUpperVertex.get();
70 DSL dsl( v[1], v[0], DSL::remainder(v[1], v[0], U) );
71 Point L = lastLowerVertex.get() + dsl.shift();
72 //NB: U (resp. L) can be the first or the last upper (resp. lower) leaning point
73 // of the DSS according to the relative order of their position.
74 ASSERT( dsl.position(v) != 0 );
75 Position qUForward = ( dsl.position(aL) - dsl.position(U) ) / dsl.position(v);
76 Position qLForward = ( dsl.position(aL) - dsl.position(L) ) / dsl.position(v);
77 Position qUBackward = ( dsl.position(U) - dsl.position(aF) ) / dsl.position(v);
78 Position qLBackward = ( dsl.position(L) - dsl.position(aF) ) / dsl.position(v);
79 return DSS( v[1], v[0], dsl.mu(), DSL::remainder(v[1], v[0], L), aF, aL,
80 (U - qUBackward * v), (U + qUForward * v ),
81 (L - qLBackward * v), (L + qLForward * v ),
82 dsl.steps(), dsl.shift() );
83 }
84}
85
86//-----------------------------------------------------------------------------
87template <typename TCoordinate, typename TInteger, unsigned short adjacency>
88inline
89DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
90DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::
91createSubsegment(const DSS& aDSS, const Point& aF, const Point& aL)
92{
93 ASSERT( aDSS.dsl().isInDSL( aF ) );
94 ASSERT( aDSS.dsl().isInDSL( aL ) );
95 ASSERT( aDSS.dsl().beforeOrEqual(aF, aL) );
96 ASSERT( aDSS.dsl().beforeOrEqual(aDSS.back(), aL) );
97 ASSERT( aDSS.dsl().beforeOrEqual(aF, aDSS.front()) );
98
99 //specific case: DSS of one point
100 if (aF == aL)
101 {
102 return DSS( aF );
103 }
104 //otherwise
105 else
106 {
107 DSS leftSubsegment = createLeftSubsegment(aDSS, aL);
108 DSS tmp = leftSubsegment.negate();
109 DSS subsegment = createLeftSubsegment(tmp, aF);
110 return subsegment.negate();
111 }
112}
113
114//-----------------------------------------------------------------------------
115template <typename TCoordinate, typename TInteger, unsigned short adjacency>
116inline
117DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
118DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::
119createLeftSubsegment(const DSS& aDSS, const Point& aL)
120{
121 //running reversedSmartCH to compute the minimal characteristics AND the leaning points
122 OneItemOutputIterator<Point> lastUpperVertex, lastLowerVertex;
123 using namespace functions;
124 Vector v = reversedSmartCH( aDSS, aDSS.position(aL),
125 lastUpperVertex, lastLowerVertex );
126
127 //DSS construction
128 Point U = lastUpperVertex.get();
129 DSL dsl( v[1], v[0], DSL::remainder(v[1], v[0], U) );
130 Point L = lastLowerVertex.get() + dsl.shift();
131 //NB: U (resp. L) can be the first or the last upper (resp. lower) leaning point
132 // of the DSS according to the relative order of their position.
133 ASSERT( dsl.position(v) != 0 );
134 Position qUForward = ( dsl.position(aL) - dsl.position(U) ) / dsl.position(v);
135 Position qLForward = ( dsl.position(aL) - dsl.position(L) ) / dsl.position(v);
136 Position qUBackward = ( dsl.position(U) - dsl.position(aDSS.back()) ) / dsl.position(v);
137 Position qLBackward = ( dsl.position(L) - dsl.position(aDSS.back()) ) / dsl.position(v);
138 return DSS( v[1], v[0], dsl.mu(), DSL::remainder(v[1], v[0], L), aDSS.back(), aL,
139 (U - qUBackward * v), (U + qUForward * v ),
140 (L - qLBackward * v), (L + qLForward * v ),
141 dsl.steps(), dsl.shift() );
142}
143//-----------------------------------------------------------------------------
144template <typename TCoordinate, typename TInteger, unsigned short adjacency>
145inline
146DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
147DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::createPattern(const Point& aF, const Point& aL)
148{
149
150 //specific case: a and b are both null
151 if (aF == aL)
152 {
153 return DSS( aF );
154 }
155 //otherwise
156 else
157 {
158
159 //upper leaning points
160 Point Uf, Ul, Lf, Ll;
161 Uf = aF;
162 Ul = aL;
163
164 //direction vector
165 Vector u2 = aL - aF;
166 IntegerComputer<Coordinate> computer;
167 Coordinate gcd = computer.gcd(u2[0], u2[1]);
168
169 //irreductible direction vector
170 Vector u = Vector(u2[0]/gcd, u2[1]/gcd);
171
172 //intercept
173 Integer mu = DSL::remainder(u[1], u[0], Uf);
174
175 //DSL
176 DSL dsl(u[1], u[0], mu);
177
178 //lower leaning points
179 if ( dsl.steps().second == Vector(NumberTraits<TCoordinate>::ZERO,NumberTraits<TCoordinate>::ZERO) )
180 { //specific case: there is only one step,
181 //ie. either a == 0, or b == 0, or
182 //(a == 1 and b == 1) in the naive case
183 Lf = Uf;
184 Ll = Ul;
185
186 return DSS(dsl.a(), dsl.b(),
187 mu, mu, aF, aL,
188 Uf, Ul, Lf, Ll,
189 dsl.steps(), dsl.shift());
190 }
191 else
192 { //otherwise
193 Coordinate signedUnity = -1;
194 Vector v = bezoutVector(u[1], u[0], signedUnity);
195 Lf = Uf + v - dsl.shift()*signedUnity;
196 Ll = Uf + u*(gcd-1) + v - dsl.shift()*signedUnity;
197
198 return DSS(dsl.a(), dsl.b(),
199 mu, DSL::remainder(u[1], u[0], Lf),
200 aF, aL,
201 Uf, Ul, Lf, Ll,
202 dsl.steps(), dsl.shift());
203 }
204 }
205}
206
207//-----------------------------------------------------------------------------
208template <typename TCoordinate, typename TInteger, unsigned short adjacency>
209inline
210DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
211DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::createDSS(const Coordinate& aA, const Coordinate& aB, const Point& aFirst, const Point& aLast, const Point& aU)
212{
213 Point Uf, Ul, Lf, Ll;
214 DSL dsl(aA,aB,DSL::remainder(aA,aB,aU));
215
216 Vector u(aB,aA);
217
218
219 //Upper leaning points
220
221 Vector t = (aLast-aU)/u;
222 Integer k1 = (t[0]>t[1])?t[1]:t[0];
223
224 Ul = aU + u*k1;
225
226 t = (aU-aFirst)/u;
227 Integer k2 = (t[0]>t[1])?t[1]:t[0];
228 Uf = aU - u*k2;
229
230
231 //Lower leaning points
232 if ( dsl.steps().second == Vector(NumberTraits<TCoordinate>::ZERO,NumberTraits<TCoordinate>::ZERO) )
233 { //specific case: there is only one step,
234 //ie. either a == 0, or b == 0, or
235 //(a == 1 and b == 1) in the naive case
236 Lf = Uf;
237 Ll = Ul;
238 }
239 else
240 { //otherwise
241 Coordinate signedUnity = -1;
242 Vector v = bezoutVector(aA, aB, signedUnity);
243 v += dsl.shift(); // "forward" vector from an upper to a lower leaning point
244 Vector w = v - u; // "backward" vector from an upper to a lower leaning point
245 if(dsl.beforeOrEqual(aFirst,Uf+w)) // there is a lower leaning point before Uf
246 Lf = Uf + w;
247 else
248 Lf = Uf + v;
249
250 if(dsl.beforeOrEqual(Ul+v,aLast))
251 Ll = Ul + v;
252 else
253 Ll = Ul + w;
254 }
255
256 DSS res(dsl.a(), dsl.b(), aFirst, aLast, Uf, Ul, Lf, Ll);
257
258 assert(res.isValid());
259
260 return res;
261
262}
263
264
265
266//-----------------------------------------------------------------------------
267template <typename TCoordinate, typename TInteger, unsigned short adjacency>
268inline
269DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
270DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::createReversedPattern(const Point& aF, const Point& aL)
271{
272 DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency> pattern = createPattern(aL, aF);
273 return pattern.negate();
274}
275
276//-----------------------------------------------------------------------------
277template <typename TCoordinate, typename TInteger, unsigned short adjacency>
278inline
279typename DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::Vector
280DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::
281bezoutVector(const Coordinate& aA,
282 const Coordinate& aB,
283 const Coordinate& aR)
284{
285 ASSERT( (aR == 1)||(aR == -1) );
286
287 //compute one Bezout vector
288 IntegerComputer<Coordinate> computer;
289 Vector v = computer.extendedEuclid(aA, -aB, aR);
290 ASSERT( (aA*v[0]-aB*v[1]) == aR );
291
292 // compute the one whose component
293 // have a sign equal to the direction
294 // vector components
295 // modify v if and only if aB and v[0] or aA and v[1] have strictly different signs
296
297 if ( (aB > 0)&&(v[0] < 0) )
298 v += Vector(aB,aA);
299 else if ( (aB < 0)&&(v[0] > 0) )
300 v += Vector(aB,aA);
301 else if ( ( aA > 0 )&&(v[1] < 0) )
302 v += Vector(aB,aA);
303 else if ( ( aA < 0 )&&(v[1] > 0) )
304 v += Vector(aB,aA);
305
306 ASSERT( (aA*v[0]-aB*v[1]) == aR );
307
308 // Assert that for the pairs (aA,v[1]) and (aB,v[0]) either one member is null, or they have the same sign
309 ASSERT( (aA ==0 || v[1] == 0 || (aA > 0 && v[1] > 0) || (aA < 0 && v[1] <0)) && (aB ==0 || v[0] == 0 || (aB > 0 && v[0] > 0) || (aB < 0 && v[0] <0)));
310
311 return v;
312}
313
314
315