DGtal  1.4.beta
SurfaceMeshHelper.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 SurfaceMeshHelper.ih
19  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21  *
22  * @date 2020/02/18
23  *
24  * Implementation of inline methods defined in SurfaceMeshHelper.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 #include <limits>
33 //////////////////////////////////////////////////////////////////////////////
34 
35 ///////////////////////////////////////////////////////////////////////////////
36 // IMPLEMENTATION of inline methods.
37 ///////////////////////////////////////////////////////////////////////////////
38 
39 ///////////////////////////////////////////////////////////////////////////////
40 // ----------------------- Standard services ------------------------------
41 
42 //-----------------------------------------------------------------------------
43 template <typename TRealPoint, typename TRealVector>
44 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::SurfaceMesh
45 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
46 makeSphere( const Scalar radius, const RealPoint& center,
47  Size m, Size n, const NormalsType normals )
48 {
49  m = std::max( m, (Size)1 ); // nb latitudes (except poles)
50  n = std::max( n, (Size)3 ); // nb longitudes
51  // vertices are numbered as follows
52  // 0 : bottom pole
53  // 1 -- m: first row
54  // m+1 -- 2*m: second row
55  // ...
56  // (n-1)*m+1 -- n*m: before last row
57  // n*m+1 : top pole
58  const Size nbv = 2 + n * m;
59  std::vector< RealPoint > p( nbv ); // positions
60  std::vector< std::vector< Size > > faces; // faces = ( incident vertices )
61  p[ 0 ] = center + radius * RealPoint( 0.0, 0.0, -1.0 );
62  p[ m*n+1 ] = center + radius * RealPoint( 0.0, 0.0, 1.0 );
63  const auto fv = [n] (Size i, Size j) -> Size { return i * n + j + (Size)1; };
64  for ( Size i = 0; i < m; ++i )
65  for ( Size j = 0; j < n; ++j )
66  {
67  const Scalar theta = ( 2.0 * M_PI * (double) j ) / (double) n;
68  const Scalar phi = ( M_PI * (double) (i+1) ) / (double) (m+1) - 0.5 * M_PI;
69  p[ fv( i, j ) ] = center + radius * RealPoint( cos( theta ) * cos( phi ),
70  sin( theta ) * cos( phi ),
71  sin( phi ) );
72  }
73  for ( Size j = 0; j < n; ++j )
74  faces.push_back( std::vector< Size > { 0u, fv( 0, (j+1) % n ), fv( 0, j ) } );
75  for ( Size i = 0; i < m-1; ++i )
76  for ( Size j = 0; j < n; ++j )
77  {
78  faces.push_back( std::vector< Size >
79  { fv( i, j ), fv( i, (j+1) % n ),
80  fv( i+1, (j+1) % n ), fv( i+1, j ) } );
81  }
82  for ( Size j = 0; j < n; ++j )
83  faces.push_back( std::vector< Size > { fv( m-1, j ), fv( m-1, (j+1) % n ), m*n+(Size)1 } );
84  SurfaceMesh smesh( p.cbegin(), p.cend(), faces.cbegin(), faces.cend() );
85  if ( normals == NormalsType::VERTEX_NORMALS )
86  {
87  std::vector< RealVector > vnormals( nbv );
88  for ( Size k = 0; k < nbv; ++k )
89  vnormals[ k ] = ( p[ k ] - center ).getNormalized();
90  smesh.setVertexNormals( vnormals.cbegin(), vnormals.cend() );
91  }
92  else if ( normals == NormalsType::FACE_NORMALS )
93  {
94  std::vector< RealVector > fnormals( smesh.nbFaces() );
95  for ( Size k = 0; k < fnormals.size(); ++k )
96  fnormals[ k ] = ( smesh.faceCentroid( k ) - center ).getNormalized();
97  smesh.setFaceNormals( fnormals.cbegin(), fnormals.cend() );
98  }
99  return smesh;
100 }
101 
102 //-----------------------------------------------------------------------------
103 template <typename TRealPoint, typename TRealVector>
104 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
105 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
106 sphereMeanCurvatures( const Scalar radius, Size m, Size n )
107 {
108  n = std::max( n, (Size)1 ); // nb latitudes (except poles)
109  m = std::max( m, (Size)3 ); // nb longitudes
110  // vertices are numbered as follows
111  // 0 : bottom pole
112  // 1 -- m: first row
113  // m+1 -- 2*m: second row
114  // ...
115  // (n-1)*m+1 -- n*m: before last row
116  // n*m+1 : top pole
117  const Size nbv = 2 + n * m;
118  return Scalars( nbv, 1.0 / radius );
119 }
120 
121 //-----------------------------------------------------------------------------
122 template <typename TRealPoint, typename TRealVector>
123 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
124 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
125 sphereGaussianCurvatures( const Scalar radius, Size m, Size n )
126 {
127  n = std::max( n, (Size)1 ); // nb latitudes (except poles)
128  m = std::max( m, (Size)3 ); // nb longitudes
129  // vertices are numbered as follows
130  // 0 : bottom pole
131  // 1 -- m: first row
132  // m+1 -- 2*m: second row
133  // ...
134  // (n-1)*m+1 -- n*m: before last row
135  // n*m+1 : top pole
136  const Size nbv = 2 + n * m;
137  return Scalars( nbv, 1.0 / ( radius * radius ) );
138 }
139 
140 //-----------------------------------------------------------------------------
141 template <typename TRealPoint, typename TRealVector>
142 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
143 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
144 sphereFirstPrincipalCurvatures( const Scalar radius, Size m, Size n )
145 {
146  n = std::max( n, (Size)1 ); // nb latitudes (except poles)
147  m = std::max( m, (Size)3 ); // nb longitudes
148  // vertices are numbered as follows
149  // 0 : bottom pole
150  // 1 -- m: first row
151  // m+1 -- 2*m: second row
152  // ...
153  // (n-1)*m+1 -- n*m: before last row
154  // n*m+1 : top pole
155  const Size nbv = 2 + n * m;
156  return Scalars( nbv, 1.0 / radius );
157 }
158 
159 //-----------------------------------------------------------------------------
160 template <typename TRealPoint, typename TRealVector>
161 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
162 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
163 sphereSecondPrincipalCurvatures( const Scalar radius, Size m, Size n )
164 {
165  n = std::max( n, (Size)1 ); // nb latitudes (except poles)
166  m = std::max( m, (Size)3 ); // nb longitudes
167  // vertices are numbered as follows
168  // 0 : bottom pole
169  // 1 -- m: first row
170  // m+1 -- 2*m: second row
171  // ...
172  // (n-1)*m+1 -- n*m: before last row
173  // n*m+1 : top pole
174  const Size nbv = 2 + n * m;
175  return Scalars( nbv, 1.0 / radius );
176 }
177 
178 //-----------------------------------------------------------------------------
179 template <typename TRealPoint, typename TRealVector>
180 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
181 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
182 sphereFirstPrincipalDirections( const Scalar , Size m, Size n )
183 {
184  // since the sphere is made of umbilic points, returned directions
185  // are arbitrary tangent vectors. Here we return latitude vectors.
186  n = std::max( n, (Size)1 ); // nb latitudes (except poles)
187  m = std::max( m, (Size)3 ); // nb longitudes
188  // vertices are numbered as follows
189  // 0 : bottom pole
190  // 1 -- m: first row
191  // m+1 -- 2*m: second row
192  // ...
193  // (n-1)*m+1 -- n*m: before last row
194  // n*m+1 : top pole
195  const Size nbv = 2 + n * m;
196  std::vector< RealVector > d1( nbv ); // directions
197  d1[ 0 ] = RealVector( 0.0, 1.0, 0.0 );
198  d1[ m*n+1 ] = RealVector( 0.0, 1.0, 0.0 );
199  const auto fv = [n] (Size i, Size j) -> Size { return i * n + j + (Size)1; };
200  for ( Size i = 0; i < m; ++i )
201  for ( Size j = 0; j < n; ++j )
202  {
203  const Scalar theta = ( 2.0 * M_PI * (double) j ) / (double) n;
204  d1[ fv( i, j ) ] = RealVector( -sin( theta ), cos( theta ), 0.0 );
205  }
206  return d1;
207 }
208 
209 //-----------------------------------------------------------------------------
210 template <typename TRealPoint, typename TRealVector>
211 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
212 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
213 sphereSecondPrincipalDirections( const Scalar , Size m, Size n )
214 {
215  // since the sphere is made of umbilic points, returned directions
216  // are arbitrary tangent vectors. Here we return latitude vectors.
217  n = std::max( n, (Size)1 ); // nb latitudes (except poles)
218  m = std::max( m, (Size)3 ); // nb longitudes
219  // vertices are numbered as follows
220  // 0 : bottom pole
221  // 1 -- m: first row
222  // m+1 -- 2*m: second row
223  // ...
224  // (n-1)*m+1 -- n*m: before last row
225  // n*m+1 : top pole
226  const Size nbv = 2 + n * m;
227  std::vector< RealVector > d2( nbv ); // directions
228  d2[ 0 ] = RealVector( 1.0, 0.0, 0.0 );
229  d2[ m*n+1 ] = RealVector( -1.0, 0.0, 0.0 );
230  const auto fv = [n] (Size i, Size j) -> Size { return i * n + j + (Size)1; };
231  for ( Size i = 0; i < m; ++i )
232  for ( Size j = 0; j < n; ++j )
233  {
234  const Scalar theta = ( 2.0 * M_PI * (double) j ) / (double) n;
235  const Scalar phi = ( M_PI * (double) (i+1) ) / (double) (m+1) - 0.5 * M_PI;
236  d2[ fv( i, j ) ] = RealVector( -cos( theta ) * sin( phi ),
237  -sin( theta ) * sin( phi ), cos( phi ) );
238  }
239  return d2;
240 }
241 
242 //-----------------------------------------------------------------------------
243 template <typename TRealPoint, typename TRealVector>
244 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::SurfaceMesh
245 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
246 makeLantern( const Scalar radius, const Scalar height, const RealPoint& center,
247  Size m, Size n, const NormalsType normals )
248 {
249  m = std::max( m, (Size)2 ); // nb latitudes
250  n = std::max( n, (Size)3 ); // nb longitudes
251  // vertices are numbered as follows
252  // 0 -- m-1: first row
253  // m -- 2m-1: second row (shifted)
254  // ...
255  // (n-1)*m -- nm-1: last row
256  const Size nbv = n * m;
257  std::vector< RealPoint > p( nbv ); // positions
258  std::vector< std::vector< Size > > faces; // faces = ( incident vertices )
259  const auto fv = [n] (Size i, Size j) -> Size { return i * n + j; };
260  for ( Size i = 0; i < m; ++i )
261  for ( Size j = 0; j < n; ++j )
262  {
263  const Scalar theta = ( i % 2 == 0 )
264  ? ( 2.0 * M_PI * (double) j ) / (double) n
265  : ( 2.0 * M_PI * ( 0.5 + (double) j ) ) / (double) n;
266  const Scalar h = height * ( (double) i / (double) (m-1) - 0.5 );
267  p[ fv( i, j ) ] = center + RealPoint( radius * cos( theta ),
268  radius * sin( theta ),
269  h );
270  }
271  for ( Size i = 0; i < m-1; ++i )
272  for ( Size j = 0; j < n; ++j )
273  if ( i % 2 == 0 )
274  {
275  faces.push_back( std::vector< Size >
276  { fv( i, j ), fv( i, (j+1) % n ), fv( i+1, j ) } );
277  faces.push_back( std::vector< Size >
278  { fv( i, (j+1) % n ), fv( i+1, (j+1) % n ), fv( i+1, j )} );
279  }
280  else
281  {
282  faces.push_back( std::vector< Size >
283  { fv( i, j ), fv( i+1, (j+1)%n ), fv( i+1, j ) } );
284  faces.push_back( std::vector< Size >
285  { fv( i+1, (j+1)%n ), fv( i, j ), fv( i, (j+1) % n ) } );
286  }
287  SurfaceMesh smesh( p.cbegin(), p.cend(), faces.cbegin(), faces.cend() );
288  if ( normals == NormalsType::VERTEX_NORMALS )
289  {
290  std::vector< RealVector > vnormals( nbv );
291  for ( Size k = 0; k < nbv; ++k )
292  {
293  RealVector _n = p[ k ] - center;
294  _n[ 2 ] = 0.0;
295  vnormals[ k ] = _n.getNormalized();
296  }
297  smesh.setVertexNormals( vnormals.cbegin(), vnormals.cend() );
298  }
299  else if ( normals == NormalsType::FACE_NORMALS )
300  {
301  std::vector< RealVector > fnormals( smesh.nbFaces() );
302  for ( Size k = 0; k < fnormals.size(); ++k )
303  {
304  RealVector _n = smesh.faceCentroid( k ) - center;
305  _n[ 2 ] = 0.0;
306  fnormals[ k ] = _n.getNormalized();
307  }
308  smesh.setFaceNormals( fnormals.cbegin(), fnormals.cend() );
309  }
310  return smesh;
311 }
312 
313 //-----------------------------------------------------------------------------
314 template <typename TRealPoint, typename TRealVector>
315 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
316 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
317 lanternMeanCurvatures( const Scalar radius, Size m, Size n )
318 {
319  n = std::max( n, (Size)2 ); // nb latitudes
320  m = std::max( m, (Size)3 ); // nb longitudes
321  const Size nbv = n * m;
322  return Scalars( nbv, 1.0 / ( 2.0 * radius ) );
323 }
324 
325 //-----------------------------------------------------------------------------
326 template <typename TRealPoint, typename TRealVector>
327 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
328 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
329 lanternGaussianCurvatures( const Scalar radius, Size m, Size n )
330 {
331  (void)radius; //unused param.
332  n = std::max( n, (Size)2 ); // nb latitudes
333  m = std::max( m, (Size)3 ); // nb longitudes
334  const Size nbv = n * m;
335  return Scalars( nbv, 0.0 );
336 }
337 
338 //-----------------------------------------------------------------------------
339 template <typename TRealPoint, typename TRealVector>
340 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
341 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
342 lanternFirstPrincipalCurvatures( const Scalar radius, Size m, Size n )
343 {
344  (void)radius; //unused param.
345  m = std::max( m, (Size)2 ); // nb latitudes
346  n = std::max( n, (Size)3 ); // nb longitudes
347  // vertices are numbered as follows
348  // 0 -- m-1: first row
349  // m -- 2m-1: second row (shifted)
350  // ...
351  // (n-1)*m -- nm-1: last row
352  const Size nbv = n * m;
353  return Scalars( nbv, 0.0 );
354 }
355 
356 //-----------------------------------------------------------------------------
357 template <typename TRealPoint, typename TRealVector>
358 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
359 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
360 lanternSecondPrincipalCurvatures( const Scalar radius, Size m, Size n )
361 {
362  (void)radius; //unused param.
363  m = std::max( m, (Size)2 ); // nb latitudes
364  n = std::max( n, (Size)3 ); // nb longitudes
365  // vertices are numbered as follows
366  // 0 -- m-1: first row
367  // m -- 2m-1: second row (shifted)
368  // ...
369  // (n-1)*m -- nm-1: last row
370  const Size nbv = n * m;
371  return Scalars( nbv, 1.0 / radius );
372 }
373 
374 //-----------------------------------------------------------------------------
375 template <typename TRealPoint, typename TRealVector>
376 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
377 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
378 lanternFirstPrincipalDirections( const Scalar radius, Size m, Size n )
379 {
380  (void)radius; //unused param.
381  m = std::max( m, (Size)2 ); // nb latitudes
382  n = std::max( n, (Size)3 ); // nb longitudes
383  // vertices are numbered as follows
384  // 0 -- m-1: first row
385  // m -- 2m-1: second row (shifted)
386  // ...
387  // (n-1)*m -- nm-1: last row
388  const Size nbv = n * m;
389  std::vector< RealVector > d1( nbv, RealVector( 0.0, 0.0, -1.0 ) ); // directions
390  return d1;
391 }
392 
393 //-----------------------------------------------------------------------------
394 template <typename TRealPoint, typename TRealVector>
395 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
396 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
397 lanternSecondPrincipalDirections( const Scalar radius, Size m, Size n )
398 {
399  (void)radius; //unused param.
400  m = std::max( m, (Size)2 ); // nb latitudes
401  n = std::max( n, (Size)3 ); // nb longitudes
402  // vertices are numbered as follows
403  // 0 -- m-1: first row
404  // m -- 2m-1: second row (shifted)
405  // ...
406  // (n-1)*m -- nm-1: last row
407  const Size nbv = n * m;
408  std::vector< RealVector > d2( nbv ); // directions
409  const auto fv = [n] (Size i, Size j) -> Size { return i * n + j; };
410  for ( Size i = 0; i < m; ++i )
411  for ( Size j = 0; j < n; ++j )
412  {
413  const Scalar theta = ( i % 2 == 0 )
414  ? ( 2.0 * M_PI * (double) j ) / (double) n
415  : ( 2.0 * M_PI * ( 0.5 + (double) j ) ) / (double) n;
416  d2[ fv( i, j ) ] = RealVector( -sin( theta ), cos( theta ), 0.0 );
417  }
418  return d2;
419 }
420 
421 //-----------------------------------------------------------------------------
422 template <typename TRealPoint, typename TRealVector>
423 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::SurfaceMesh
424 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
425 makeTorus( const Scalar big_radius, const Scalar small_radius, const RealPoint& center,
426  Size m, Size n, const int twist, const NormalsType normals )
427 {
428  m = std::max( m, (Size)3 ); // nb latitudes
429  n = std::max( n, (Size)3 ); // nb longitudes
430  // vertices are numbered as follows
431  // 0 -- m-1: first row
432  // m -- 2m-1: second row
433  // ...
434  // (n-1)*m -- nm-1: last row
435  const Size nbv = n * m;
436  std::vector< RealPoint > p( nbv ); // positions
437  std::vector< std::vector< Size > > faces; // faces = ( incident vertices )
438  const auto fv = [n] (Size i, Size j) -> Size { return i * n + j; };
439  for ( Size i = 0; i < m; ++i )
440  for ( Size j = 0; j < n; ++j )
441  {
442  const Scalar theta = ( 2.0 * M_PI * (double) i ) / (double) m;
443  const Scalar phi = ( 2.0 * M_PI * (double) j ) / (double) n
444  + twist * 2.0 * M_PI * (double) i / ( double ) m;
445  p[ fv( i, j ) ] = center
446  + RealPoint( ( big_radius + small_radius * cos( phi ) ) * cos( theta ),
447  ( big_radius + small_radius * cos( phi ) ) * sin( theta ),
448  small_radius * sin( phi ) );
449  }
450  for ( Size i = 0; i < m; ++i )
451  for ( Size j = 0; j < n; ++j )
452  {
453  faces.push_back( std::vector< Size >
454  { fv( i, j ), fv( (i+1)%m, j ), fv( i, (j+1) % n ) } );
455  faces.push_back( std::vector< Size >
456  { fv( i, (j+1) % n ), fv( (i+1)%m, j ), fv( (i+1)%m, (j+1) % n ) } );
457  }
458  SurfaceMesh smesh( p.cbegin(), p.cend(), faces.cbegin(), faces.cend() );
459  if ( normals == NormalsType::VERTEX_NORMALS )
460  {
461  std::vector< RealVector > vnormals( nbv );
462  for ( Size i = 0; i < m; ++i )
463  {
464  const Scalar theta = ( 2.0 * M_PI * (double) i ) / (double) m;
465  const RealPoint c = center + RealPoint( big_radius * cos( theta ),
466  big_radius * sin( theta ), 0.0 );
467  for ( Size j = 0; j < n; ++j )
468  vnormals[ fv( i, j ) ] = ( p[ fv( i, j ) ] - c ).getNormalized();
469  }
470  smesh.setVertexNormals( vnormals.cbegin(), vnormals.cend() );
471  }
472  else if ( normals == NormalsType::FACE_NORMALS )
473  {
474  std::vector< RealVector > fnormals( smesh.nbFaces() );
475  Size k = 0;
476  for ( Size i = 0; i < m; ++i )
477  {
478  const Scalar theta = ( 2.0 * M_PI * (double) i ) / (double) m;
479  const RealPoint c = center + RealPoint( big_radius * cos( theta ),
480  big_radius * sin( theta ), 0.0 );
481  for ( Size j = 0; j < n; ++j )
482  {
483  fnormals[ k ] = ( smesh.faceCentroid( k ) - c ).getNormalized();
484  fnormals[ k+1 ] = ( smesh.faceCentroid( k+1 ) - c ).getNormalized();
485  k += 2;
486  }
487  }
488  smesh.setFaceNormals( fnormals.cbegin(), fnormals.cend() );
489  }
490  return smesh;
491 }
492 
493 //-----------------------------------------------------------------------------
494 template <typename TRealPoint, typename TRealVector>
495 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
496 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
497 torusMeanCurvatures( const Scalar big_radius, const Scalar small_radius,
498  Size m, Size n, const int twist )
499 {
500  m = std::max( m, (Size)3 ); // nb latitudes
501  n = std::max( n, (Size)3 ); // nb longitudes
502  const Size nbv = n * m;
503  Scalars hvalues( 2*nbv );
504  // H = - ( R + 2r cos phi ) / ( 2r(R + r cos phi) )
505  // G = cos phi / ( r(R + r cos phi) )
506  auto H = [&] ( int i, int j )
507  {
508  const Scalar phi = ( 2.0 * M_PI * (double) j ) / (double) n
509  + twist * 2.0 * M_PI * (double) i / ( double ) m;
510  return ( big_radius + 2.0 * small_radius * cos( phi ) )
511  / ( 2.0 * small_radius * ( big_radius + small_radius * cos( phi ) ) );
512  };
513  Index f = 0;
514  for ( Size i = 0; i < m; ++i )
515  for ( Size j = 0; j < n; ++j )
516  {
517  hvalues[ f++ ] = ( H( i, j ) + H( (i+1)%m, j ) + H( i, (j+1) % n ) ) / 3.0;
518  hvalues[ f++ ] = ( H( i, (j+1)%n ) + H( (i+1)%m, j ) + H( (i+1)%m, (j+1)%n ) ) / 3.0;
519  }
520  return hvalues;
521 }
522 
523 //-----------------------------------------------------------------------------
524 template <typename TRealPoint, typename TRealVector>
525 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
526 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
527 torusGaussianCurvatures( const Scalar big_radius, const Scalar small_radius,
528  Size m, Size n, const int twist )
529 {
530  m = std::max( m, (Size)3 ); // nb latitudes
531  n = std::max( n, (Size)3 ); // nb longitudes
532  const Size nbv = n * m;
533  Scalars gvalues( 2*nbv );
534  // H = - ( R + 2r cos phi ) / ( 2r(R + r cos phi) )
535  // G = cos phi / ( r(R + r cos phi) )
536  auto G = [&] ( int i, int j )
537  {
538  const Scalar phi = ( 2.0 * M_PI * (double) j ) / (double) n
539  + twist * 2.0 * M_PI * (double) i / ( double ) m;
540  return cos( phi )
541  / ( small_radius * ( big_radius + small_radius * cos( phi ) ) );
542  };
543  Index f = 0;
544  for ( Size i = 0; i < m; ++i )
545  for ( Size j = 0; j < n; ++j )
546  {
547  gvalues[ f++ ] = ( G( i, j ) + G( (i+1)%m, j ) + G( i, (j+1) % n ) ) / 3.0;
548  gvalues[ f++ ] = ( G( i, (j+1)%n ) + G( (i+1)%m, j ) + G( (i+1)%m, (j+1)%n ) ) / 3.0;
549  }
550  return gvalues;
551 
552 }
553 
554 //-----------------------------------------------------------------------------
555 template <typename TRealPoint, typename TRealVector>
556 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
557 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
558 torusFirstPrincipalCurvatures( const Scalar big_radius, const Scalar small_radius,
559  Size m, Size n, const int twist )
560 {
561  const auto H = torusMeanCurvatures( big_radius, small_radius, m, n, twist );
562  const auto G = torusGaussianCurvatures( big_radius, small_radius, m, n, twist );
563  // k1 = H - sqrt( fabs( H * H - G ))
564  Scalars k1( H.size() );
565  const auto compute_k1 =
566  [&] ( Size i ) { return H[ i ] - sqrt( fabs( H[ i ] * H[ i ] - G[ i ] ) ); };
567  for ( Size i = 0; i < k1.size(); ++i ) k1[ i ] = compute_k1( i );
568  return k1;
569 }
570 
571 //-----------------------------------------------------------------------------
572 template <typename TRealPoint, typename TRealVector>
573 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
574 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
575 torusSecondPrincipalCurvatures( const Scalar big_radius, const Scalar small_radius,
576  Size m, Size n, const int twist )
577 {
578  const auto H = torusMeanCurvatures( big_radius, small_radius, m, n, twist );
579  const auto G = torusGaussianCurvatures( big_radius, small_radius, m, n, twist );
580  // k2 = H + sqrt( fabs( H * H - G ))
581  Scalars k2( H.size() );
582  const auto compute_k2 =
583  [&] ( Size i ) { return H[ i ] + sqrt( fabs( H[ i ] * H[ i ] - G[ i ] ) ); };
584  for ( Size i = 0; i < k2.size(); ++i ) k2[ i ] = compute_k2( i );
585  return k2;
586 }
587 
588 //-----------------------------------------------------------------------------
589 template <typename TRealPoint, typename TRealVector>
590 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
591 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
592 torusFirstPrincipalDirections( const Scalar big_radius, const Scalar small_radius,
593  Size m, Size n, const int twist )
594 {
595  (void)big_radius;//not used
596  (void)small_radius;
597  (void)twist;
598 
599  m = std::max( m, (Size)3 ); // nb latitudes
600  n = std::max( n, (Size)3 ); // nb longitudes
601  // vertices are numbered as follows
602  // 0 -- m-1: first row
603  // m -- 2m-1: second row
604  // ...
605  // (n-1)*m -- nm-1: last row
606  const Size nbv = n * m;
607  std::vector< RealVector > d1( nbv ); // directions
608  const auto fv = [n] (Size i, Size j) -> Size { return i * n + j; };
609  for ( Size i = 0; i < m; ++i )
610  for ( Size j = 0; j < n; ++j )
611  {
612  const Scalar theta = ( 2.0 * M_PI * (double) i ) / (double) m;
613  d1[ fv( i, j ) ] = RealVector( sin( theta ), cos( theta ), 0.0 );
614  }
615  return d1;
616 }
617 
618 //-----------------------------------------------------------------------------
619 template <typename TRealPoint, typename TRealVector>
620 typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
621 DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
622 torusSecondPrincipalDirections( const Scalar big_radius, const Scalar small_radius,
623  Size m, Size n, const int twist )
624 {
625  (void)big_radius; //not used
626  (void)small_radius; //not used
627 
628  m = std::max( m, (Size)3 ); // nb latitudes
629  n = std::max( n, (Size)3 ); // nb longitudes
630  // vertices are numbered as follows
631  // 0 -- m-1: first row
632  // m -- 2m-1: second row
633  // ...
634  // (n-1)*m -- nm-1: last row
635  const Size nbv = n * m;
636  std::vector< RealVector > d2( nbv ); // directions
637  const auto fv = [n] (Size i, Size j) -> Size { return i * n + j; };
638  for ( Size i = 0; i < m; ++i )
639  for ( Size j = 0; j < n; ++j )
640  {
641  const Scalar theta = ( 2.0 * M_PI * (double) i ) / (double) m;
642  const Scalar phi = ( 2.0 * M_PI * (double) j ) / (double) n
643  + twist * 2.0 * M_PI * (double) i / ( double ) m;
644  d2[ fv( i, j ) ] = RealVector( -( sin( phi ) ) * cos( theta ),
645  -( sin( phi ) ) * sin( theta ),
646  cos( phi ) );
647  }
648  return d2;
649 }
650 
651 ///////////////////////////////////////////////////////////////////////////////
652 ///////////////////////////////////////////////////////////////////////////////