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