DGtal  0.9.2
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 //-----------------------------------------------------------------------------
43 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
44 inline
45 DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
46 DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::
47 createSubsegment(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 //-----------------------------------------------------------------------------
87 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
88 inline
89 DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
90 DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::
91 createSubsegment(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 //-----------------------------------------------------------------------------
115 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
116 inline
117 DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
118 DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::
119 createLeftSubsegment(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 //-----------------------------------------------------------------------------
144 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
145 inline
146 DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
147 DGtal::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 //-----------------------------------------------------------------------------
208 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
209 inline
210 DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
211 DGtal::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 //-----------------------------------------------------------------------------
267 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
268 inline
269 DGtal::ArithmeticalDSS<TCoordinate,TInteger,adjacency>
270 DGtal::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 //-----------------------------------------------------------------------------
277 template <typename TCoordinate, typename TInteger, unsigned short adjacency>
278 inline
279 typename DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::Vector
280 DGtal::ArithmeticalDSSFactory<TCoordinate,TInteger,adjacency>::
281 bezoutVector(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