DGtal 1.3.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 radius, 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 radius, 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 n = std::max( n, (Size)2 ); // nb latitudes
332 m = std::max( m, (Size)3 ); // nb longitudes
333 const Size nbv = n * m;
334 return Scalars( nbv, 0.0 );
335}
336
337//-----------------------------------------------------------------------------
338template <typename TRealPoint, typename TRealVector>
339typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
340DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
341lanternFirstPrincipalCurvatures( const Scalar radius, Size m, Size n )
342{
343 m = std::max( m, (Size)2 ); // nb latitudes
344 n = std::max( n, (Size)3 ); // nb longitudes
345 // vertices are numbered as follows
346 // 0 -- m-1: first row
347 // m -- 2m-1: second row (shifted)
348 // ...
349 // (n-1)*m -- nm-1: last row
350 const Size nbv = n * m;
351 return Scalars( nbv, 0.0 );
352}
353
354//-----------------------------------------------------------------------------
355template <typename TRealPoint, typename TRealVector>
356typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
357DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
358lanternSecondPrincipalCurvatures( const Scalar radius, Size m, Size n )
359{
360 m = std::max( m, (Size)2 ); // nb latitudes
361 n = std::max( n, (Size)3 ); // nb longitudes
362 // vertices are numbered as follows
363 // 0 -- m-1: first row
364 // m -- 2m-1: second row (shifted)
365 // ...
366 // (n-1)*m -- nm-1: last row
367 const Size nbv = n * m;
368 return Scalars( nbv, 1.0 / radius );
369}
370
371//-----------------------------------------------------------------------------
372template <typename TRealPoint, typename TRealVector>
373typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
374DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
375lanternFirstPrincipalDirections( const Scalar radius, Size m, Size n )
376{
377 m = std::max( m, (Size)2 ); // nb latitudes
378 n = std::max( n, (Size)3 ); // nb longitudes
379 // vertices are numbered as follows
380 // 0 -- m-1: first row
381 // m -- 2m-1: second row (shifted)
382 // ...
383 // (n-1)*m -- nm-1: last row
384 const Size nbv = n * m;
385 std::vector< RealVector > d1( nbv, RealVector( 0.0, 0.0, -1.0 ) ); // directions
386 return d1;
387}
388
389//-----------------------------------------------------------------------------
390template <typename TRealPoint, typename TRealVector>
391typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
392DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
393lanternSecondPrincipalDirections( const Scalar radius, Size m, Size n )
394{
395 m = std::max( m, (Size)2 ); // nb latitudes
396 n = std::max( n, (Size)3 ); // nb longitudes
397 // vertices are numbered as follows
398 // 0 -- m-1: first row
399 // m -- 2m-1: second row (shifted)
400 // ...
401 // (n-1)*m -- nm-1: last row
402 const Size nbv = n * m;
403 std::vector< RealVector > d2( nbv ); // directions
404 const auto fv = [n] (Size i, Size j) -> Size { return i * n + j; };
405 for ( Size i = 0; i < m; ++i )
406 for ( Size j = 0; j < n; ++j )
407 {
408 const Scalar theta = ( i % 2 == 0 )
409 ? ( 2.0 * M_PI * (double) j ) / (double) n
410 : ( 2.0 * M_PI * ( 0.5 + (double) j ) ) / (double) n;
411 d2[ fv( i, j ) ] = RealVector( -sin( theta ), cos( theta ), 0.0 );
412 }
413 return d2;
414}
415
416//-----------------------------------------------------------------------------
417template <typename TRealPoint, typename TRealVector>
418typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::SurfaceMesh
419DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
420makeTorus( const Scalar big_radius, const Scalar small_radius, const RealPoint& center,
421 Size m, Size n, const int twist, const NormalsType normals )
422{
423 m = std::max( m, (Size)3 ); // nb latitudes
424 n = std::max( n, (Size)3 ); // nb longitudes
425 // vertices are numbered as follows
426 // 0 -- m-1: first row
427 // m -- 2m-1: second row
428 // ...
429 // (n-1)*m -- nm-1: last row
430 const Size nbv = n * m;
431 std::vector< RealPoint > p( nbv ); // positions
432 std::vector< std::vector< Size > > faces; // faces = ( incident vertices )
433 const auto fv = [n] (Size i, Size j) -> Size { return i * n + j; };
434 for ( Size i = 0; i < m; ++i )
435 for ( Size j = 0; j < n; ++j )
436 {
437 const Scalar theta = ( 2.0 * M_PI * (double) i ) / (double) m;
438 const Scalar phi = ( 2.0 * M_PI * (double) j ) / (double) n
439 + twist * 2.0 * M_PI * (double) i / ( double ) m;
440 p[ fv( i, j ) ] = center
441 + RealPoint( ( big_radius + small_radius * cos( phi ) ) * cos( theta ),
442 ( big_radius + small_radius * cos( phi ) ) * sin( theta ),
443 small_radius * sin( phi ) );
444 }
445 for ( Size i = 0; i < m; ++i )
446 for ( Size j = 0; j < n; ++j )
447 {
448 faces.push_back( std::vector< Size >
449 { fv( i, j ), fv( (i+1)%m, j ), fv( i, (j+1) % n ) } );
450 faces.push_back( std::vector< Size >
451 { fv( i, (j+1) % n ), fv( (i+1)%m, j ), fv( (i+1)%m, (j+1) % n ) } );
452 }
453 SurfaceMesh smesh( p.cbegin(), p.cend(), faces.cbegin(), faces.cend() );
454 if ( normals == NormalsType::VERTEX_NORMALS )
455 {
456 std::vector< RealVector > vnormals( nbv );
457 for ( Size i = 0; i < m; ++i )
458 {
459 const Scalar theta = ( 2.0 * M_PI * (double) i ) / (double) m;
460 const RealPoint c = center + RealPoint( big_radius * cos( theta ),
461 big_radius * sin( theta ), 0.0 );
462 for ( Size j = 0; j < n; ++j )
463 vnormals[ fv( i, j ) ] = ( p[ fv( i, j ) ] - c ).getNormalized();
464 }
465 smesh.setVertexNormals( vnormals.cbegin(), vnormals.cend() );
466 }
467 else if ( normals == NormalsType::FACE_NORMALS )
468 {
469 std::vector< RealVector > fnormals( smesh.nbFaces() );
470 Size k = 0;
471 for ( Size i = 0; i < m; ++i )
472 {
473 const Scalar theta = ( 2.0 * M_PI * (double) i ) / (double) m;
474 const RealPoint c = center + RealPoint( big_radius * cos( theta ),
475 big_radius * sin( theta ), 0.0 );
476 for ( Size j = 0; j < n; ++j )
477 {
478 fnormals[ k ] = ( smesh.faceCentroid( k ) - c ).getNormalized();
479 fnormals[ k+1 ] = ( smesh.faceCentroid( k+1 ) - c ).getNormalized();
480 k += 2;
481 }
482 }
483 smesh.setFaceNormals( fnormals.cbegin(), fnormals.cend() );
484 }
485 return smesh;
486}
487
488//-----------------------------------------------------------------------------
489template <typename TRealPoint, typename TRealVector>
490typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
491DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
492torusMeanCurvatures( const Scalar big_radius, const Scalar small_radius,
493 Size m, Size n, const int twist )
494{
495 m = std::max( m, (Size)3 ); // nb latitudes
496 n = std::max( n, (Size)3 ); // nb longitudes
497 const Size nbv = n * m;
498 Scalars hvalues( 2*nbv );
499 // H = - ( R + 2r cos phi ) / ( 2r(R + r cos phi) )
500 // G = cos phi / ( r(R + r cos phi) )
501 auto H = [&] ( int i, int j )
502 {
503 const Scalar phi = ( 2.0 * M_PI * (double) j ) / (double) n
504 + twist * 2.0 * M_PI * (double) i / ( double ) m;
505 return ( big_radius + 2.0 * small_radius * cos( phi ) )
506 / ( 2.0 * small_radius * ( big_radius + small_radius * cos( phi ) ) );
507 };
508 Index f = 0;
509 for ( Size i = 0; i < m; ++i )
510 for ( Size j = 0; j < n; ++j )
511 {
512 hvalues[ f++ ] = ( H( i, j ) + H( (i+1)%m, j ) + H( i, (j+1) % n ) ) / 3.0;
513 hvalues[ f++ ] = ( H( i, (j+1)%n ) + H( (i+1)%m, j ) + H( (i+1)%m, (j+1)%n ) ) / 3.0;
514 }
515 return hvalues;
516}
517
518//-----------------------------------------------------------------------------
519template <typename TRealPoint, typename TRealVector>
520typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
521DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
522torusGaussianCurvatures( const Scalar big_radius, const Scalar small_radius,
523 Size m, Size n, const int twist )
524{
525 m = std::max( m, (Size)3 ); // nb latitudes
526 n = std::max( n, (Size)3 ); // nb longitudes
527 const Size nbv = n * m;
528 Scalars gvalues( 2*nbv );
529 // H = - ( R + 2r cos phi ) / ( 2r(R + r cos phi) )
530 // G = cos phi / ( r(R + r cos phi) )
531 auto G = [&] ( int i, int j )
532 {
533 const Scalar phi = ( 2.0 * M_PI * (double) j ) / (double) n
534 + twist * 2.0 * M_PI * (double) i / ( double ) m;
535 return cos( phi )
536 / ( small_radius * ( big_radius + small_radius * cos( phi ) ) );
537 };
538 Index f = 0;
539 for ( Size i = 0; i < m; ++i )
540 for ( Size j = 0; j < n; ++j )
541 {
542 gvalues[ f++ ] = ( G( i, j ) + G( (i+1)%m, j ) + G( i, (j+1) % n ) ) / 3.0;
543 gvalues[ f++ ] = ( G( i, (j+1)%n ) + G( (i+1)%m, j ) + G( (i+1)%m, (j+1)%n ) ) / 3.0;
544 }
545 return gvalues;
546
547}
548
549//-----------------------------------------------------------------------------
550template <typename TRealPoint, typename TRealVector>
551typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
552DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
553torusFirstPrincipalCurvatures( const Scalar big_radius, const Scalar small_radius,
554 Size m, Size n, const int twist )
555{
556 const auto H = torusMeanCurvatures( big_radius, small_radius, m, n, twist );
557 const auto G = torusGaussianCurvatures( big_radius, small_radius, m, n, twist );
558 // k1 = H - sqrt( fabs( H * H - G ))
559 Scalars k1( H.size() );
560 const auto compute_k1 =
561 [&] ( Size i ) { return H[ i ] - sqrt( fabs( H[ i ] * H[ i ] - G[ i ] ) ); };
562 for ( Size i = 0; i < k1.size(); ++i ) k1[ i ] = compute_k1( i );
563 return k1;
564}
565
566//-----------------------------------------------------------------------------
567template <typename TRealPoint, typename TRealVector>
568typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::Scalars
569DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
570torusSecondPrincipalCurvatures( const Scalar big_radius, const Scalar small_radius,
571 Size m, Size n, const int twist )
572{
573 const auto H = torusMeanCurvatures( big_radius, small_radius, m, n, twist );
574 const auto G = torusGaussianCurvatures( big_radius, small_radius, m, n, twist );
575 // k2 = H + sqrt( fabs( H * H - G ))
576 Scalars k2( H.size() );
577 const auto compute_k2 =
578 [&] ( Size i ) { return H[ i ] + sqrt( fabs( H[ i ] * H[ i ] - G[ i ] ) ); };
579 for ( Size i = 0; i < k2.size(); ++i ) k2[ i ] = compute_k2( i );
580 return k2;
581}
582
583//-----------------------------------------------------------------------------
584template <typename TRealPoint, typename TRealVector>
585typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
586DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
587torusFirstPrincipalDirections( const Scalar big_radius, const Scalar small_radius,
588 Size m, Size n, const int twist )
589{
590 m = std::max( m, (Size)3 ); // nb latitudes
591 n = std::max( n, (Size)3 ); // nb longitudes
592 // vertices are numbered as follows
593 // 0 -- m-1: first row
594 // m -- 2m-1: second row
595 // ...
596 // (n-1)*m -- nm-1: last row
597 const Size nbv = n * m;
598 std::vector< RealVector > d1( nbv ); // directions
599 const auto fv = [n] (Size i, Size j) -> Size { return i * n + j; };
600 for ( Size i = 0; i < m; ++i )
601 for ( Size j = 0; j < n; ++j )
602 {
603 const Scalar theta = ( 2.0 * M_PI * (double) i ) / (double) m;
604 d1[ fv( i, j ) ] = RealVector( sin( theta ), cos( theta ), 0.0 );
605 }
606 return d1;
607}
608
609//-----------------------------------------------------------------------------
610template <typename TRealPoint, typename TRealVector>
611typename DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::RealVectors
612DGtal::SurfaceMeshHelper<TRealPoint, TRealVector>::
613torusSecondPrincipalDirections( const Scalar big_radius, const Scalar small_radius,
614 Size m, Size n, const int twist )
615{
616 m = std::max( m, (Size)3 ); // nb latitudes
617 n = std::max( n, (Size)3 ); // nb longitudes
618 // vertices are numbered as follows
619 // 0 -- m-1: first row
620 // m -- 2m-1: second row
621 // ...
622 // (n-1)*m -- nm-1: last row
623 const Size nbv = n * m;
624 std::vector< RealVector > d2( nbv ); // directions
625 const auto fv = [n] (Size i, Size j) -> Size { return i * n + j; };
626 for ( Size i = 0; i < m; ++i )
627 for ( Size j = 0; j < n; ++j )
628 {
629 const Scalar theta = ( 2.0 * M_PI * (double) i ) / (double) m;
630 const Scalar phi = ( 2.0 * M_PI * (double) j ) / (double) n
631 + twist * 2.0 * M_PI * (double) i / ( double ) m;
632 d2[ fv( i, j ) ] = RealVector( -( sin( phi ) ) * cos( theta ),
633 -( sin( phi ) ) * sin( theta ),
634 cos( phi ) );
635 }
636 return d2;
637}
638
639///////////////////////////////////////////////////////////////////////////////
640///////////////////////////////////////////////////////////////////////////////
Generated on Sun Nov 27 2022 15:15:51 for DGtal by doxygen 1.9.5