DGtal  0.9.2
OrderedAlphabet.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 OrderedAlphabet.ih
19  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20  * Laboratory of Mathematics (CNRS, UMR 5807), University of Savoie, France
21  * @author Laurent Provot (\c Laurent.Provot@loria.fr )
22  * LORIA (CNRS, UMR 7503), Nancy University, France
23  *
24  * @date 2010/07/01
25  *
26  * Implementation of inline methods defined in OrderedAlphabet.h
27  *
28  * This file is part of the DGtal library.
29  */
30 
31 ///////////////////////////////////////////////////////////////////////////////
32 // IMPLEMENTATION of inline methods.
33 ///////////////////////////////////////////////////////////////////////////////
34 
35 //////////////////////////////////////////////////////////////////////////////
36 #include <cstdlib>
37 //////////////////////////////////////////////////////////////////////////////
38 
39 
40 
41 ///////////////////////////////////////////////////////////////////////////////
42 // Implementation of inline methods //
43 
44 /**
45  * Constructor from letters
46  *
47  * @param first the first letter of the alphabet.
48  * @param nb the number of letters of the alphabet.
49  *
50  * Exemple: OrderedAlphabet( '0', 4 ) defines the alphabet for
51  * 4-connected freeman chains.
52  */
53 inline
54 DGtal::OrderedAlphabet::OrderedAlphabet( char first, unsigned int nb )
55  : myFirst( first ), myNb( nb )
56 {
57  myOrder = new unsigned int[ nb ];
58 
59  ASSERT( ( myOrder != 0 )
60  && "[DGtal::OrderedAlphabet::OrderedAlphabet( char first, int nb )] error in new: no memory left ?" );
61  for ( unsigned int i = 0; i < myNb; ++i )
62  myOrder[ i ] = i;
63 }
64 
65 /**
66  * @param c any valid letter in this alphabet.
67  *
68  * @return the index of the letter [c] in the order relation,
69  * starting from 0 to myNb-1.
70  */
71 inline
72 unsigned int
73 DGtal::OrderedAlphabet::order( char c ) const
74 {
75  ASSERT( ( c - myFirst ) < (int) myNb
76  && "[DGtal::OrderedAlphabet::order( char c )] invalid letter." );
77  return myOrder[ c - myFirst ];
78 }
79 
80 /**
81  * @param i the index of some letter in the order relation,
82  * between 0 and myNb-1.
83  *
84  * @return c the corresponding letter in this alphabet.
85  *
86  * NB: O(nb of letters in the alphabet).
87  */
88 inline
89 char
90 DGtal::OrderedAlphabet::letter( unsigned int i ) const
91 {
92  ASSERT( i < myNb );
93  for ( unsigned int j = 0; j < myNb; ++j )
94  if ( myOrder[ j ] == i )
95  return myFirst + j;
96  return myFirst;
97 }
98 
99 
100 /**
101  * @param c1 a letter in the alphabet
102  * @param c2 another letter in the same alphabet.
103  * @return 'true' iff c1 < c2
104  */
105 inline
106 bool
107 DGtal::OrderedAlphabet::less( char c1, char c2 ) const
108 {
109  return myOrder[ c1 - myFirst ] < myOrder[ c2 - myFirst ];
110 }
111 
112 /**
113  * @param c1 a letter in the alphabet
114  * @param c2 another letter in the same alphabet.
115  * @return 'true' iff c1 <= c2
116  */
117 inline
118 bool
119 DGtal::OrderedAlphabet::lessOrEqual( char c1, char c2 ) const
120 {
121  return myOrder[ c1 - myFirst ] <= myOrder[ c2 - myFirst ];
122 }
123 
124 /**
125  * @param c1 a letter in the alphabet
126  * @param c2 another letter in the same alphabet.
127  * @return 'true' iff c1 == c2
128  */
129 inline
130 bool
131 DGtal::OrderedAlphabet::equal( char c1, char c2 ) const
132 {
133  return c1 == c2;
134 }
135 
136 
137 ///////////////////////////////////////////////////////////////////////////////
138 // Implementation of inline functions and external operators //
139 
140 /**
141  * Overloads 'operator<<' for displaying objects of class 'OrderedAlphabet'.
142  * @param out the output stream where the object is written.
143  * @param object the object of class 'OrderedAlphabet' to write.
144  * @return the output stream after the writing.
145  */
146 inline
147 std::ostream&
148 DGtal::operator<< ( std::ostream & out,
149  const OrderedAlphabet & object )
150 {
151  object.selfDisplay ( out );
152  return out;
153 }
154 
155 
156 
157 // //
158 ///////////////////////////////////////////////////////////////////////////////
159 
160 
161 
162 /**
163  * Destructor.
164  */
165 inline
166 DGtal::OrderedAlphabet::~OrderedAlphabet()
167 {
168  if ( myOrder != 0 )
169  delete[ ] myOrder;
170 }
171 
172 /**
173  * @return the current ordered alphabet.
174  */
175 inline
176 std::string
177 DGtal::OrderedAlphabet::orderedAlphabet() const
178 {
179  char *tbl;
180  tbl = (char *)malloc((myNb + 1)*sizeof(char));
181  for ( unsigned int i = 0; i < myNb; ++i )
182  {
183  tbl[ myOrder[ i ] ] = myFirst + i;
184  }
185  tbl[ myNb ] = '\0';
186  std::string s( tbl );
187  free(tbl);
188  return s;
189 }
190 
191 /**
192  * Shift a0 < a1 < ... < an to a1 < ... < an < a0
193  */
194 inline
195 void
196 DGtal::OrderedAlphabet::shiftLeft()
197 {
198  unsigned int* p = myOrder;
199  unsigned int* q = myOrder + myNb;
200  while ( p != q )
201  {
202  unsigned int k = *p;
203  *p = ( k == 0 ) ? myNb - 1 : k - 1;
204  ++p;
205  }
206 }
207 
208 /**
209  * Shift a0 < a1 < ... < an to an < a0 < ... < an-1
210  */
211 inline
212 void
213 DGtal::OrderedAlphabet::shiftRight()
214 {
215  unsigned int* p = myOrder;
216  unsigned int* q = myOrder + myNb;
217  while ( p != q )
218  {
219  unsigned int k = *p + 1;
220  *p = ( k == myNb ) ? 0 : k;
221  ++p;
222  }
223 }
224 
225 /**
226  * Reverse the order a0 < a1 < ... < an to an < ... < a1 < a0
227  */
228 inline
229 void
230 DGtal::OrderedAlphabet::reverse()
231 {
232  unsigned int* p = myOrder;
233  unsigned int* q = myOrder + myNb;
234  while ( p != q )
235  {
236  *p = myNb - 1 - *p;
237  ++p;
238  }
239 }
240 
241 /**
242  * Reverse the order a0 < a1 < ... < an to a3 < a2 < a1 < a0 < an < ...
243  */
244 inline
245 void
246 DGtal::OrderedAlphabet::reverseAround12()
247 {
248  unsigned int* p = myOrder;
249  unsigned int* q = myOrder + myNb;
250  while ( p != q )
251  {
252  *p = ( myNb + 3 - *p ) % myNb;
253  ++p;
254  }
255 }
256 
257 
258 
259 
260 /**
261  * Gives the first lyndon factor of the word [w] starting at
262  * position [s] in this alphabet.
263  *
264  * @param len (returns) the length of the primitive Lyndon factor
265  * (which starts at position s).
266  *
267  * @param nb (returns) the number of times the Lyndon factor appears.
268  *
269  * @param w a word
270  * @param s the starting index in [w].
271  * @param e the index after the end in [w] (s<e).
272  */
273 inline
274 void
275 DGtal::OrderedAlphabet::firstLyndonFactor
276 ( size_t & len, size_t & nb,
277  const std::string & w,
278  index_t s, index_t e ) const
279 {
280  index_t i = s;
281  index_t j = s+1;
282  while ( ( j < e ) && ( lessOrEqual( w[ i ], w[ j ] ) ) )
283  {
284  if ( equal( w[ i ], w[ j ] ) )
285  ++i;
286  else
287  i = s;
288  ++j;
289  }
290  len = (size_t) j - i;
291  nb = ( (size_t) ( j - s ) ) / len;
292 }
293 
294 
295 /**
296  * Gives the first lyndon factor of the cyclic word [w]
297  * starting at position [s] in this alphabet.
298  *
299  * @param len (returns) the length of the primitive Lyndon factor
300  * (which starts at position s).
301  *
302  * @param nb (returns) the number of times the Lyndon factor appears.
303  *
304  * @param w a word
305  * @param s the starting index in [w].
306  * @param e the index after the end in [w] (s and e arbitrary).
307  */
308 inline
309 void
310 DGtal::OrderedAlphabet::firstLyndonFactorMod
311 ( size_t & len, size_t & nb,
312  const std::string & w,
313  index_t s, index_t e ) const
314 {
315  size_t modulo = (DGtal::OrderedAlphabet::size_t)w.size();
316  DGtal::ModuloComputer< Integer > mc( modulo );
317  index_t i = s;
318  index_t j = mc.next( s );
319  while ( ( j != e ) && ( lessOrEqual( w[ i ], w[ j ] ) ) )
320  {
321  if ( equal( w[ i ], w[ j ] ) )
322  mc.increment( i );
323  else
324  i = s;
325  mc.increment( j );
326  }
327  len = j >= i ? (size_t) ( j - i )
328  : (size_t) ( j + modulo - i );
329  if (len == 0)
330  nb = 0;
331  else
332  nb = ( (size_t) ( ( j + modulo - s ) % modulo ) ) / len;
333 }
334 
335 
336 /**
337  * @param len (returns) the length of the primitive Lyndon factor
338  * (which starts at position s).
339  *
340  * @param nb (returns) the number of times the Lyndon factor appears.
341  *
342  * @param w a word which starts with a1 or a2 at position s.
343  * @param s the starting index in [w].
344  * @param e the index after the end in [w] (s<e).
345  */
346 inline
347 bool DGtal::OrderedAlphabet::duvalPP( size_t & len, size_t & nb,
348  const std::string & w, index_t s, index_t e) const
349 {
350  ASSERT( ( order( w[ s ] ) == 1 ) || ( order( w[ s ] ) == 2 ) );
351  index_t i = s;
352  index_t j = s+1;
353  index_t p = s;
354  index_t q = s+1;
355  while ( ( j < e ) && ( lessOrEqual( w[ i ], w[ j ] ) ) )
356  {
357  if ( equal( w[ i ], w[ j ] ) )
358  {
359  if ( j == q )
360  q += (p-s+1);
361  ++i;
362  }
363  else
364  {
365  if ( ( j != q ) || ( order ( w[ j ] ) != 2 ) )
366  {
367  len = j-s; nb = 0;
368  return false;
369  }
370  index_t tmp = p;
371  p = q;
372  q += q - tmp;
373  i = s;
374  }
375  ++j;
376  }
377  len = (size_t) j - i;
378  nb = ( (size_t) (j-s) ) / len;
379  return true;
380 }
381 
382 /**
383  * @param len (returns) the length of the primitive Lyndon factor
384  * (which starts at position s).
385  *
386  * @param nb (returns) the number of times the Lyndon factor appears.
387  *
388  * @param n1 (returns) the number of occurrences of the letter a1
389  * in the Lyndon factor
390  *
391  * @param n2 (returns) the number of occurrences of the letter a2
392  * in the Lyndon factor
393  *
394  * @param Lf1 (returns) the number of occurrences of the letter a1
395  * from 's' to the first lower leaning point.
396  *
397  * @param Lf2 (returns) the number of occurrences of the letter a2
398  * from 's' to the first lower leaning point.
399  *
400  * @param w a word which starts with a1 or a2 at position s.
401  * @param s the starting index in [w].
402  * @param e the index after the end in [w] (s<e).
403  */
404 inline
405 bool
406 DGtal::OrderedAlphabet::duvalPPtoDSS
407 ( size_t & len, size_t & nb,
408  unsigned int & n1, unsigned int & n2,
409  unsigned int & lf1, unsigned int & lf2,
410  const std::string & w,
411  index_t s, index_t e
412  ) const
413 {
414  ASSERT( ( order( w[ s ] ) == 1 ) || ( order( w[ s ] ) == 2 ) );
415  index_t i = s;
416  index_t j = s+1;
417  index_t p = s;
418  index_t q = s+1;
419  unsigned int slope1 = (order( w[ i ] ) == 1) ? 1 : 0;
420  unsigned int slope2 = (order( w[ i ] ) == 2) ? 1 : 0;
421  lf1 = n1 = slope1;
422  lf2 = n2 = slope2;
423  nb = 1;
424 
425  while ( ( j < e ) && ( lessOrEqual( w[ i ], w[ j ] ) ) ) {
426 
427  //This 'if/else if' is added to compute the vector defined by
428  //the Christoffel word, this is usefull in order to compute the
429  //leaning points.
430  if (order( w[ j ] ) == 1)
431  ++slope1;
432  else if (order( w[ j ] ) == 2)
433  ++slope2;
434 
435  if ( equal( w[ i ], w[ j ] ) ) {
436  if ( j == q ) {
437  q += (p-s+1);
438 
439  //A repetition is read when j==s+(nb+1)*(p-s+1)-1
440  }
441  if ( j == nb * (p-s+1) + p ) {
442  ++nb;
443  }
444  ++i;
445  } else {
446  if ( ( j != q ) || ( order ( w[ j ] ) != 2 ) ) {
447  break;
448  }
449  index_t tmp = p;
450  p = q;
451  q += q - tmp;
452  i = s;
453 
454  // Extra information to compute the leaning points
455  lf1 = lf1 + (nb-1)*n1;
456  lf2 = lf2 + (nb-1)*n2;
457  n1 = slope1;
458  n2 = slope2;
459  nb = 1;
460  }
461  ++j;
462  }
463  len = (size_t) j - i;
464 
465  ASSERT( nb == (size_t) (j-s) / len);
466  return true;
467 }
468 
469 
470 /**
471  * Adaptation of Duval's algorithm to extract the first Lyndon factor
472  * (FLF). Whilst scanning the Lyndon factor, it also checks whether it
473  * is a Christoffel word or not. It returns 'true' if the FLF is
474  * indeed a Christoffel word, otherwise returns false. It starts the
475  * extraction at position [s] in the cyclic word [w].
476  *
477  * The alphabet takes the form a0 < a1 < a2 < ... < an-1. [w] starts
478  * with a1 or a2 at position s.
479  *
480  * See [Provencal, Lachaud 2009].
481  *
482  * @param len (returns) the length of the primitive Lyndon factor
483  * (which starts at position s).
484  *
485  * @param nb (returns) the number of times the Lyndon factor appears.
486  *
487  * @param w a (cyclic) word which starts with a1 or a2 at position s.
488  *
489  * @param s the starting index in [w].
490  * @param e the index after the end in [w] (s and e arbitrary).
491  */
492 inline
493 bool
494 DGtal::OrderedAlphabet::duvalPPMod( size_t & len, size_t & nb,
495  const std::string & w,
496  index_t s, index_t e ) const
497 {
498  ASSERT( ( order( w[ s ] ) == 1 )
499  || ( order( w[ s ] ) == 2 ) );
500  size_t modulo = (DGtal::OrderedAlphabet::size_t)w.size();
501  ModuloComputer< Integer > mc( modulo );
502  ModuloComputer< Integer >::UnsignedInteger i = s;
503  index_t j = mc.next( s );
504  unsigned int p = 1;
505  unsigned int q = 2;
506  while ( ( j != e ) && ( lessOrEqual( w[ i ], w[ j ] ) ) )
507  {
508  if ( equal( w[ i ], w[ j ] ) )
509  {
510  if ( j == mc.cast( s + q - 1 ) )
511  q += p;
512  mc.increment( i );
513  }
514  else
515  {
516  if ( ( j != mc.cast( s + q - 1 ) ) || ( order ( w[ j ] ) != 2 ) )
517  {
518  len = j; nb = 0;
519  return false;
520  }
521  unsigned int tmp = p;
522  p = q;
523  q += q - tmp;
524  i = s;
525  }
526  mc.increment( j );
527  }
528  len = j >= i ? (size_t) ( j - i )
529  : (size_t) ( j + modulo - i );
530  nb = ( (size_t) ( ( j + modulo - s ) % modulo ) ) / len;
531  return true;
532 }
533 
534 
535 ///////////////////////////////////////////////////////////////////////////////
536 // Interface - public :
537 
538 /**
539  * Writes/Displays the object on an output stream.
540  * @param out the output stream where the object is written.
541  */
542 inline
543 void
544 DGtal::OrderedAlphabet::selfDisplay ( std::ostream & out ) const
545 {
546  out << "[OrderedAlphabet]";
547  out << " " << orderedAlphabet() << std::endl;
548 }
549 
550 /**
551  * Checks the validity/consistency of the object.
552  * @return 'true' if the object is valid, 'false' otherwise.
553  */
554 inline
555 bool
556 DGtal::OrderedAlphabet::isValid() const
557 {
558  return true;
559 }
560 
561 
562 
563 ///////////////////////////////////////////////////////////////////////////////
564 // ----------------------- MLP services ------------------------------
565 
566 /**
567  * Extracts the next edge of the minimum length polygon starting from
568  * position [s] on the word [w]. The alphabet may be modified
569  * (reversed or shifted). The output alphabet is some a0 < a1 < a2 < ...
570  *
571  * @param nb_a1 (returns) the number of letters a1 in the extracted edge (a1
572  * in the output alphabet)
573  *
574  * @param nb_a2 (returns) the number of letters a2 in the extracted edge (a2
575  * in the output alphabet)
576  *
577  * @param w the input (cyclic) word (may be modified in the process).
578  *
579  * @param s the starting point in the word (updated).
580  *
581  * @param cvx (updates) this boolean is flipped only if a change of
582  * convexity is detected.
583  *
584  * @return the number of letters of the extracted edge.
585  */
586 inline
587 DGtal::OrderedAlphabet::size_t
588 DGtal::OrderedAlphabet::nextEdge( size_t & nb_a1,
589  size_t & nb_a2,
590  std::string & w,
591  index_t & s,
592  bool & cvx )
593 {
594  ModuloComputer< Integer > mc( (const unsigned int)w.size() );
595  size_t l;
596  size_t len;
597  size_t nb;
598  bool inC = duvalPPMod( len, nb, w, s, s );
599  if ( ! inC )
600  // case : change of convexity
601  {
602  // JOL : temporary change of letter w[ s ]
603  char tmp = w[ s ];
604  index_t tmp_s = s;
605  w[ s ] = letter( 2 ); // a3
606  this->reverseAround12();
607  cvx = ! cvx;
608  l = nextEdge( nb_a1, nb_a2, w, s, cvx );
609  // JOL : former letter is put back in string.
610  w[ tmp_s ] = tmp;
611  }
612  else if ( ( len == 1 ) && ( order( w[ s ] ) == 1 ) )
613  // case u=a1 => Quadrant change
614  {
615  this->shiftRight();
616  s = mc.cast( s + nb );
617  nb_a1 = 0;
618  nb_a2 = nb - 1;
619  l = nb;
620  }
621  else
622  { // standard (convex) case.
623  l = len * nb;
624  char a2 = letter( 2 );
625  nb_a1 = len;
626  nb_a2 = 0;
627  index_t ss = s;
628  s = mc.cast( s + l );
629  while ( len != 0 )
630  {
631  if ( w[ ss ] == a2 ) ++nb_a2;
632  mc.increment( ss );
633  --len;
634  }
635  nb_a1 -= nb_a2;
636  nb_a1 *= nb;
637  nb_a2 *= nb;
638  }
639  return l;
640 }
641 
642 
643 ///////////////////////////////////////////////////////////////////////////////
644 // Internals - private :
645 
646 // //
647 ///////////////////////////////////////////////////////////////////////////////