1 |
|
/******************************************************************** |
2 |
|
* |
3 |
|
* Author: Manfredi Leto (Xemet) |
4 |
|
* License: GPL Version 2 |
5 |
|
* System: Linux |
6 |
|
* |
7 |
|
* Copyright (c) 2009 All rights reserved. |
8 |
|
* |
9 |
|
********************************************************************/ |
10 |
|
|
11 |
|
/* Those functions are needed to calculate NURBS points */ |
12 |
|
|
13 |
|
#include <math.h> |
14 |
|
#include <algorithm> |
15 |
|
#include "canon.hh" |
16 |
|
|
17 |
|
static void unit(PLANE_POINT &p) { |
18 |
|
double h = hypot(p.X, p.Y); |
19 |
|
if(h != 0) { p.X/=h; p.Y/=h; } |
20 |
|
} |
21 |
|
|
22 |
1 |
std::vector<unsigned int> knot_vector_creator(unsigned int n, unsigned int k) { |
23 |
|
|
24 |
|
unsigned int i; |
25 |
|
std::vector<unsigned int> knot_vector; |
26 |
9 |
for (i=0; i<=n+k; i++) { |
27 |
8 |
if (i < k) |
28 |
6 |
knot_vector.push_back(0); |
29 |
5 |
else if (i >= k && i <= n) |
30 |
4 |
knot_vector.push_back(i - k + 1); |
31 |
|
else |
32 |
6 |
knot_vector.push_back(n - k + 2); |
33 |
|
} |
34 |
1 |
return knot_vector; |
35 |
|
|
36 |
|
} |
37 |
|
|
38 |
4200 |
double Nmix(unsigned int i, unsigned int k, double u, |
39 |
|
std::vector<unsigned int> knot_vector) { |
40 |
|
|
41 |
4200 |
if (k == 1){ |
42 |
5376 |
if ((u >= knot_vector[i]) && (u <= knot_vector[i+1])) { |
43 |
|
return 1; |
44 |
|
} else { |
45 |
1296 |
return 0;} |
46 |
2184 |
} else if (k > 1) { |
47 |
7056 |
if ((knot_vector[i+k-1]-knot_vector[i] == 0) && |
48 |
1512 |
(knot_vector[i+k]-knot_vector[i+1] != 0)) { |
49 |
2016 |
return ((knot_vector[i+k] - u)*Nmix(i+1,k-1,u,knot_vector))/ |
50 |
2016 |
(knot_vector[i+k]-knot_vector[i+1]); |
51 |
5544 |
} else if ((knot_vector[i+k]-knot_vector[i+1] == 0) && |
52 |
504 |
(knot_vector[i+k-1]-knot_vector[i] != 0)) { |
53 |
1512 |
return ((u - knot_vector[i])*Nmix(i,k-1,u,knot_vector))/ |
54 |
2016 |
(knot_vector[i+k-1]-knot_vector[i]); |
55 |
1176 |
} else if ((knot_vector[i+k-1]-knot_vector[i] == 0) && |
56 |
|
(knot_vector[i+k]-knot_vector[i+1] == 0)) { |
57 |
|
return 0; |
58 |
|
} else { |
59 |
3528 |
return ((u - knot_vector[i])*Nmix(i,k-1,u,knot_vector))/ |
60 |
7056 |
(knot_vector[i+k-1]-knot_vector[i]) + ((knot_vector[i+k] - u)* |
61 |
5880 |
Nmix(i+1,k-1,u,knot_vector))/(knot_vector[i+k]-knot_vector[i+1]); |
62 |
|
} |
63 |
|
} |
64 |
|
else return -1; |
65 |
|
} |
66 |
|
|
67 |
|
|
68 |
|
|
69 |
140 |
double Rden(double u, unsigned int k, |
70 |
|
std::vector<CONTROL_POINT> nurbs_control_points, |
71 |
|
std::vector<unsigned int> knot_vector) { |
72 |
|
|
73 |
|
unsigned int i; |
74 |
140 |
double d = 0.0; |
75 |
1680 |
for (i=0; i<(nurbs_control_points.size()); i++) |
76 |
2100 |
d = d + Nmix(i,k,u,knot_vector)*nurbs_control_points[i].W; |
77 |
140 |
return d; |
78 |
|
} |
79 |
|
|
80 |
14 |
PLANE_POINT nurbs_point(double u, unsigned int k, |
81 |
|
std::vector<CONTROL_POINT> nurbs_control_points, |
82 |
|
std::vector<unsigned int> knot_vector) { |
83 |
|
|
84 |
|
unsigned int i; |
85 |
|
PLANE_POINT point; |
86 |
14 |
point.X = 0; |
87 |
14 |
point.Y = 0; |
88 |
168 |
for (i=0; i<(nurbs_control_points.size()); i++) { |
89 |
280 |
point.X = point.X + nurbs_control_points[i].X*Nmix(i,k,u,knot_vector) |
90 |
350 |
*nurbs_control_points[i].W/Rden(u,k,nurbs_control_points,knot_vector); |
91 |
280 |
point.Y = point.Y + nurbs_control_points[i].Y*Nmix(i,k,u,knot_vector) |
92 |
350 |
*nurbs_control_points[i].W/Rden(u,k,nurbs_control_points,knot_vector); |
93 |
|
} |
94 |
14 |
return point; |
95 |
|
} |
96 |
|
|
97 |
|
#define DU (1e-5) |
98 |
|
PLANE_POINT nurbs_tangent(double u, unsigned int k, |
99 |
|
std::vector<CONTROL_POINT> nurbs_control_points, |
100 |
|
std::vector<unsigned int> knot_vector) { |
101 |
|
unsigned int n = nurbs_control_points.size() - 1; |
102 |
|
double umax = n - k + 2; |
103 |
|
double ulo = std::max(0.0, u-DU), uhi = std::min(umax, u+DU); |
104 |
|
PLANE_POINT P1 = nurbs_point(ulo, k, nurbs_control_points, knot_vector); |
105 |
|
PLANE_POINT P3 = nurbs_point(uhi, k, nurbs_control_points, knot_vector); |
106 |
|
PLANE_POINT r = {(P3.X - P1.X) / (uhi-ulo), (P3.Y - P1.Y) / (uhi-ulo)}; |
107 |
|
unit(r); |
108 |
|
return r; |
109 |
|
} |