1+
2+ #include < cmath>
3+ #include < fstream>
4+ #include < iostream>
5+ #include < vector>
6+
7+ #include < CGAL/Simple_cartesian.h>
8+ #include < CGAL/IO/polygon_soup_io.h>
9+
10+ using Kernel = CGAL ::Simple_cartesian<double >;
11+ using Point = Kernel::Point_3;
12+ using Vector = Kernel::Vector_3;
13+ using Triangle = std::array<int , 3 >;
14+
15+ int main (int argc, char * argv[])
16+ {
17+ const std::string filename = argc > 1 ? argv[1 ] : " knot.off"
18+
19+ const int p = argc > 2 ? std::atoi (argv[2 ]) : 1 ; // Nb of major turns
20+ const int q = argc > 3 ? std::atoi (argv[3 ]) : 6 ; // Nb of knot turns
21+
22+ const double R = argc > 4 ? std::atof (argv[4 ]) : 2.0 ; // major radius
23+ const double r = argc > 5 ? std::atof (argv[5 ]) : 0.9 ; // knot radius
24+ const double tube_radius = argc > 6 ? std::atof (argv[6 ]) : 0.4 ;
25+
26+ const int tubular_segments = argc > 7 ? std::atoi (argv[7 ]) : 120 ;
27+ const int radial_segments = argc > 8 ? std::atoi (argv[8 ]) : 60 ;
28+
29+ std::vector<Point> vertices;
30+ std::vector<Triangle> faces;
31+
32+ // Generate knot polyline
33+ auto knot = [&](double t){
34+ double cq = std::cos (q * t);
35+ double sq = std::sin (q * t);
36+ double factor = R + r * cq;
37+ return Point (factor * std::cos (p * t),
38+ factor * std::sin (p * t),
39+ r * sq);
40+ };
41+ auto knot_derivative = [&](double t){
42+ double cq = std::cos (q * t);
43+ double sq = std::sin (q * t);
44+ double cp = std::cos (p * t);
45+ double sp = std::sin (p * t);
46+
47+ double factor = R + r*cq;
48+ double dfactor = -r*q*sq;
49+
50+ return Vector (dfactor*cp - factor*p*sp,
51+ dfactor*sp + factor*p*cp,
52+ r*q*cq);
53+ };
54+
55+ const double dt = 2.0 * M_PI / tubular_segments;
56+ const double dr = 2.0 * M_PI / radial_segments;
57+
58+ // A normal to the knot axis
59+ Vector N (1 ,0 ,0 );
60+ for (int i=0 ; i<tubular_segments; ++i)
61+ {
62+ double t = i * dt;
63+
64+ // knot point and its tangent
65+ Point C = knot (t);
66+ Vector T = knot_derivative (t);
67+ // Compute bases of the plane on C
68+ Vector B = CGAL::cross_product (N, T);
69+ N = -CGAL::cross_product (B, T);
70+ N /= CGAL::approximate_sqrt (N.squared_length ());
71+ B /= CGAL::approximate_sqrt (B.squared_length ());
72+
73+ for (int j=0 ; j<radial_segments; ++j){
74+ double phi = j * dr;
75+ vertices.push_back (C + tube_radius*(std::cos (phi)*N + std::sin (phi)*B));
76+ }
77+ }
78+
79+ auto idx = [&](int i, int j){
80+ i %= tubular_segments;
81+ j %= radial_segments;
82+ return i * radial_segments + j;
83+ };
84+
85+ for (int i=0 ; i<tubular_segments-1 ; ++i)
86+ for (int j=0 ; j<radial_segments; ++j){
87+ int v00 = idx (i , j );
88+ int v10 = idx (i+1 , j );
89+ int v01 = idx (i , j+1 );
90+ int v11 = idx (i+1 , j+1 );
91+
92+ faces.push_back ({v00,v10,v11});
93+ faces.push_back ({v00,v11,v01});
94+ }
95+
96+ // To avoid a twist, we compute the offset that minimize the distance between a point at tubular index N-1 and tubular index 0
97+ int n = tubular_segments-1 ;
98+ int offset = 0 ;
99+ double min = CGAL::squared_distance (vertices[idx (n, 0 )], vertices[idx (0 ,0 )]);
100+ for (int j=1 ; j<radial_segments; ++j){
101+ double sq = CGAL::squared_distance (vertices[idx (n, 0 )], vertices[idx (0 ,j)]);
102+ if (min >sq){
103+ min = sq;
104+ offset = j;
105+ }
106+ }
107+ for (int j=0 ; j<radial_segments; ++j){
108+ int v00 = idx (n, j );
109+ int v10 = idx (0 , j+offset);
110+ int v01 = idx (n, j+1 );
111+ int v11 = idx (0 , j+offset+1 );
112+
113+ faces.push_back ({v00,v10,v11});
114+ faces.push_back ({v00,v11,v01});
115+ }
116+
117+ CGAL::IO::write_polygon_soup (filename, vertices, faces);
118+
119+ return 0 ;
120+ }
0 commit comments