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