DGtal 1.4.0
Loading...
Searching...
No Matches
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//-----------------------------------------------------------------------------
43template <typename TRealPoint, typename TRealVector>
44typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::SurfaceMesh
45DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
46makeSphere( 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//-----------------------------------------------------------------------------
103template <typename TRealPoint, typename TRealVector>
104typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
105DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
106sphereMeanCurvatures( 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//-----------------------------------------------------------------------------
122template <typename TRealPoint, typename TRealVector>
123typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
124DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
125sphereGaussianCurvatures( 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//-----------------------------------------------------------------------------
141template <typename TRealPoint, typename TRealVector>
142typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
143DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
144sphereFirstPrincipalCurvatures( 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//-----------------------------------------------------------------------------
160template <typename TRealPoint, typename TRealVector>
161typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
162DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
163sphereSecondPrincipalCurvatures( 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//-----------------------------------------------------------------------------
179template <typename TRealPoint, typename TRealVector>
180typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
181DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
182sphereFirstPrincipalDirections( 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//-----------------------------------------------------------------------------
210template <typename TRealPoint, typename TRealVector>
211typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
212DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
213sphereSecondPrincipalDirections( 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//-----------------------------------------------------------------------------
243template <typename TRealPoint, typename TRealVector>
244typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::SurfaceMesh
245DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
246makeLantern( 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//-----------------------------------------------------------------------------
314template <typename TRealPoint, typename TRealVector>
315typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
316DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
317lanternMeanCurvatures( 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//-----------------------------------------------------------------------------
326template <typename TRealPoint, typename TRealVector>
327typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
328DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
329lanternGaussianCurvatures( 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//-----------------------------------------------------------------------------
339template <typename TRealPoint, typename TRealVector>
340typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
341DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
342lanternFirstPrincipalCurvatures( 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//-----------------------------------------------------------------------------
357template <typename TRealPoint, typename TRealVector>
358typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
359DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
360lanternSecondPrincipalCurvatures( 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//-----------------------------------------------------------------------------
375template <typename TRealPoint, typename TRealVector>
376typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
377DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
378lanternFirstPrincipalDirections( 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//-----------------------------------------------------------------------------
394template <typename TRealPoint, typename TRealVector>
395typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
396DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
397lanternSecondPrincipalDirections( 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//-----------------------------------------------------------------------------
422template <typename TRealPoint, typename TRealVector>
423typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::SurfaceMesh
424DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
425makeTorus( 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//-----------------------------------------------------------------------------
494template <typename TRealPoint, typename TRealVector>
495typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
496DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
497torusMeanCurvatures( 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//-----------------------------------------------------------------------------
524template <typename TRealPoint, typename TRealVector>
525typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
526DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
527torusGaussianCurvatures( 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//-----------------------------------------------------------------------------
555template <typename TRealPoint, typename TRealVector>
556typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
557DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
558torusFirstPrincipalCurvatures( 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//-----------------------------------------------------------------------------
572template <typename TRealPoint, typename TRealVector>
573typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
574DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
575torusSecondPrincipalCurvatures( 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//-----------------------------------------------------------------------------
589template <typename TRealPoint, typename TRealVector>
590typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
591DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
592torusFirstPrincipalDirections( 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//-----------------------------------------------------------------------------
619template <typename TRealPoint, typename TRealVector>
620typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
621DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
622torusSecondPrincipalDirections( 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///////////////////////////////////////////////////////////////////////////////