1 /*
2 * NurbsArc.cpp
3 *
4 * Copyright (C) 2003 Th. Rothermel
5 *
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program (see the file "COPYING" for details); if
18 * not, write to the Free Software Foundation, Inc., 675 Mass Ave,
19 * Cambridge, MA 02139, USA.
20 */
21
22 #include <math.h>
23 #include <stdio.h>
24 #include <stdlib.h>
25
26 #ifndef _STDAFX_H
27 #include "stdafx.h"
28 #endif
29
30 #ifndef _WIN32
31 # include "stdlib.h"
32 #endif
33
34 #ifndef _NURBSARC_H
35 #include "NurbsArc.h"
36 #endif
37 #ifndef _NURBS_CURVE_DEGREE_ELEVATE_H
38 #include "NurbsCurveDegreeElevate.h"
39 #endif
40
41 //definition of "numerical 0", change it as you like (and the program keeps working...)
42 #define NUMZERO 1e-6
43
NurbsArc(int narcs,float degree,Vec3f sPoint,Vec3f point,Vec3f vector,int pDegree)44 NurbsArc::NurbsArc(int narcs, float degree, Vec3f sPoint, Vec3f point, Vec3f vector, int pDegree)
45 {
46 m_valid = true;
47 /*Build an arc from NURBS somewhere in space, but orthogonal to a straight line
48 defined by point and vector. sPoint is the starting point of the arc.
49 The algorithm has been adapted from Piegl/Tiller "The NURBS Book",
50 2nd ed., p. 346,
51 */
52
53 /*
54 if(pDegree<2){
55 assert((pDegree<2) != 0);
56 pDegree = 2;
57 }
58
59 if(degree>=(180*narcs)){
60 assert((degree>=(180*narcs)) != 0);
61 narcs = narcs + 1; //Otherwise it may happen, that a arc is build with one point at
62 //infinity (-->weight=(cos(M_PI))).
63 }
64 */
65
66 Vec3f rPoint, X, Y, P0, T0, P1, P2, T2;
67 int index;
68 float r;
69 float angle = 0;
70 float rDegree = degree / 180 * M_PI;
71 float delta = (float) (rDegree/narcs);
72 float w1 = cos(delta / 2);
73 MyArray<float> cosines, sines;
74
75 for(int i=1; i<=narcs; i++){
76 angle = angle + delta;
77 cosines[i] = cos(angle);
78 sines[i] = sin(angle);
79 }
80
81 rPoint = PointToLine(point, vector, sPoint);
82 X = sPoint - rPoint;
83 r = X.length();
84
85 if(r<NUMZERO){
86 /* Points that are have a distance < NUMZERO from axis of rotation
87 will be considered to be on the axis
88 */
89 points[0] = rPoint;
90 weights[0] = 1;
91 index = 0;
92 for(int i=1; i<=narcs; i++){
93 points[index+1] = rPoint;
94 points[index+2] = rPoint;
95 weights[index+1] = w1;
96 weights[index+2] = 1;
97 index = index + 2;
98 }
99 }
100 else{
101 X.normalize();
102 Y = vector.cross(X);
103 P0 = sPoint;
104 T0 = Y;
105 points[0] = P0;
106 weights[0] = 1;
107 index = 0;
108 angle = 0;
109 for(int i=1; i<=narcs; i++){
110 P2 = rPoint + X * r * cosines[i] + Y * r * sines[i];
111 points[index+2] = P2;
112 weights[index+2] = 1;
113 T2 = X * (-sines[i]) + Y * cosines[i];
114 //Lines intersect?
115 if (Intersect3DLines(P0, T0, P2, T2) != 1) {
116 m_valid = false;
117 return;
118 }
119 points[index+1] = m_intersection;
120 weights[index+1] = w1;
121 P0 = P2;
122 T0 = T2;
123 index = index + 2;
124 }
125 }
126
127 //setup knotvector
128 makeKnotvector(narcs);
129
130 //if pDegree>2 degree elevate arc
131
132 if(pDegree>2)
133 {
134 int dimension = getPointSize();
135 int upDegree = pDegree - 2;
136 int deg = 2;
137 int k;
138
139 k = getPointSize();
140
141 Vec3f *tPoints = new Vec3f[k];
142 float *tWeights = new float[k];
143
144 for(int i=0; i<k; i++){
145 tPoints[i] = getControlPoints(i);
146 tWeights[i] = getWeights(i);
147 }
148
149 k = knots.size();
150
151 MyArray<float> tKnots(k);
152
153 for(int i=0; i<k; i++){
154 tKnots[i] = getKnots(i);
155 }
156
157 NurbsCurveDegreeElevate elevatedArc(tPoints, tWeights, tKnots, dimension,
158 deg, upDegree);
159
160 for(int i=0; i<elevatedArc.getPointSize(); i++){
161 points[i] = elevatedArc.getControlPoints(i);
162 weights[i] = elevatedArc.getWeights(i);
163 }
164
165 for(int i=0; i<elevatedArc.getKnotSize(); i++){
166 knots[i] = elevatedArc.getKnots(i);
167 }
168 }
169 }
170
171 int
Intersect3DLines(Vec3f & P1,Vec3f & T1,Vec3f & P2,Vec3f & T2)172 NurbsArc::Intersect3DLines(Vec3f &P1, Vec3f &T1, Vec3f &P2, Vec3f &T2)
173 {
174 Vec3f N;
175
176 float t, det;
177
178 if(T1==T2){
179 return 0;
180 }
181
182 det = T1.x * T2.y * (P2.z - P1.z) + T1.y * T2.z * (P2.x - P1.x) + T1.z * T2.x * (P2.y - P1.y) -
183 ((P2.x - P1.x) * T2.y * T1.z + (P2.y - P1.y) * T2.z * T1.x + (P2.z - P1.z) * T2.x *T1.y);
184
185 if(fabs(det)>(NUMZERO)){
186 return 0;
187 }
188
189
190 N = (T1.cross(T2)).cross(T2);
191 t = N.dot(P2-P1) / N.dot(T1);
192 m_intersection = P1 + T1 * t;
193
194 return 1 ;
195 }
196
197 Vec3f
PointToLine(Vec3f lPoint,Vec3f lVector,Vec3f point)198 NurbsArc::PointToLine(Vec3f lPoint, Vec3f lVector, Vec3f point)
199 {
200 Vec3f rPoint;
201 float t;
202
203 assert(lVector.length()>NUMZERO);
204 t = (lVector.dot(point - lPoint)) / (lVector.length()*lVector.length());
205
206 rPoint = lPoint + lVector * t;
207
208 return rPoint;
209 }
210
211
212 void
makeKnotvector(int narcs)213 NurbsArc::makeKnotvector(int narcs)
214 {
215
216 int n = 2 * narcs;
217 int j = n + 1;
218 for(int i=0; i<3; i++){
219 knots[i] = 0;
220 knots[i+j] = 1;
221 }
222 int index = 3;
223 float k = 1;
224 for(int i=1; i<narcs; i++){
225 knots[index] = k / narcs;
226 knots[index+1] = k / narcs;
227 index = index + 2;
228 k = k + 1;
229 }
230
231 }
232
233