1 /*=============================================================================
2 File: nurbsSub.cpp
3 Purpose:
4 Revision: $Id: nurbsSub.cpp,v 1.2 2002/05/13 21:07:46 philosophil Exp $
5 Created by: Philippe Lavoie (20 Januray, 1999)
6 Modified by:
7
8 Copyright notice:
9 Copyright (C) 1999 Philippe Lavoie
10
11 This library is free software; you can redistribute it and/or
12 modify it under the terms of the GNU Library General Public
13 License as published by the Free Software Foundation; either
14 version 2 of the License, or (at your option) any later version.
15
16 This library is distributed in the hope that it will be useful,
17 but WITHOUT ANY WARRANTY; without even the implied warranty of
18 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 Library General Public License for more details.
20
21 You should have received a copy of the GNU Library General Public
22 License along with this library; if not, write to the Free
23 Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
24 =============================================================================*/
25 #include "nurbsSub.h"
26 #include <iostream>
27 #include <cstdio>
28 #include <stdio.h>
29
30 /*!
31 */
32 namespace PLib {
33
34
35 const int FALSE_ = 0 ;
36 const int TRUE_ = 1 ;
37
38 struct BOOL{
39 int value ;
BOOLPLib::BOOL40 BOOL(): value(FALSE_) {;}
BOOLPLib::BOOL41 BOOL(int a) { if(a>0) value=TRUE_; else value=FALSE_;}
BOOLPLib::BOOL42 BOOL(const BOOL& b) { value = b.value ; }
operator !PLib::BOOL43 int operator!() { if(value>0) return FALSE_; return TRUE_ ;}
operator intPLib::BOOL44 operator int() const {return value;}
operator =PLib::BOOL45 int& operator=(int a) { value = a ; return value ; }
46 };
47
48 /*!
49 \brief Constructor from a NurbsSurface
50
51 Constructor from a NurbsSurface
52
53 \param s the NurbsSurface to construct from
54
55 \author Philippe Lavoie
56 \date 20 January 1999
57 */
58 template <class T>
NurbsSubSurface(const NurbsSurface<T,3> & s)59 NurbsSubSurface<T>::NurbsSubSurface(const NurbsSurface<T,3>& s): rsurf(s) {
60 surf = 0 ;
61 render = 0 ;
62 //initSurf();
63
64 }
65
66 #define MAXORDER 20 /* Maximum order allowed (for local array sizes) */
67
68 template <class T>
69 struct NurbSurface {
70 /* Number of Points in the U and V directions, respectivly */
71 int numU, numV;
72 /* Order of the surface in U and V (must be >= 2, < MAXORDER) */
73 int orderU, orderV;
74 /* Knot vectors, indexed as [0..numU+orderU-1] and [0..numV+orderV-1] */
75 T * kvU, * kvV;
76 /* Control points, indexed as points[0..numV-1][0..numU-1] */
77 /* Note the w values are *premultiplied* with the x, y and z values */
78 Matrix<HPoint_nD<T,3> >& points;
79
80 /* These fields are added to support subdivision */
81 BOOL strV0, strVn, /* Edge straightness flags for subdivision */
82 strU0, strUn;
83 BOOL flatV, flatU; /* Surface flatness flags for subdivision */
84 SurfSample<T> c00, c0n,
85 cn0, cnn; /* Corner data structures for subdivision */
86 RenderMesh<T> *render;
87
88 static T epsilon ;
89
NurbSurfacePLib::NurbSurface90 NurbSurface(): points(tmpPoints) { render = 0 ; }
91 NurbSurface(Matrix<HPoint_nD<T,3> >& s, int du, int dv);
92
93 protected:
94 Matrix<HPoint_nD<T,3> > tmpPoints ;
95
96 };
97
98 #define CHECK( n ) \
99 { if (!(n)) { fprintf( stderr, "Ran out of memory\n" ); exit(-1); } }
100
101 #define DIVW( rpt, pt ) \
102 { (pt)->x() = (rpt)->x() / (rpt)->w(); \
103 (pt)->y() = (rpt)->y() / (rpt)->w(); \
104 (pt)->z() = (rpt)->z() / (rpt)->w(); }
105
106
107 /* Function prototypes */
108
109 template <class T> void DrawSubdivision( NurbSurface<T> *, T tolerance );
110 template <class T> void DrawEvaluation( NurbSurface<T> * );
111
112 template <class T> int FindBreakPoint( T u, T * kv, int m, int k );
113 template <class T> void AllocNurb( NurbSurface<T> *, T *, T * );
114 template <class T> void CloneNurb( NurbSurface<T> *, NurbSurface<T> * );
115 template <class T> void FreeNurb( NurbSurface<T> * );
116 template <class T> void RefineSurface( NurbSurface<T> *, NurbSurface<T> *, BOOL );
117
118 template <class T> void CalcPoint( T, T, NurbSurface<T> *, Point_nD<T,3> *, Point_nD<T,3> *, Point_nD<T,3> * );
119
120
121 /*!
122 \brief Draw the subdivision of the NURBS surface
123
124 Draw the subdivision of the NURBS surface
125
126 \param tolerance the accepted tolerance
127
128 \author Philippe Lavoie
129 \date 20 January 1999
130 */
131 template <class T>
drawSubdivision(T tolerance)132 void NurbsSubSurface<T>::drawSubdivision(T tolerance)
133 {
134 initSurf();
135 surf->render = render ;
136 DrawSubdivision(surf,tolerance);
137 }
138
139 /*!
140 \brief perform the subdivision of the NURBS and write the result in a PS file.
141
142 \param f the file name to write to
143 \param tolerance the accepted tolerance
144
145 \author Philippe Lavoie
146 \date 20 January 1999
147 */
148 template <class T>
drawSubdivisionPS(const char * f,T tolerance)149 void NurbsSubSurface<T>::drawSubdivisionPS(const char* f, T tolerance)
150 {
151 ofstream fout ;
152 fout.open(f) ;
153 if(fout)
154 drawSubdivisionPS(fout,tolerance) ;
155 fout.close();
156 }
157
158 /*!
159 \brief perform the subdivision of the NURBS and write the result in a VRML file.
160
161
162 \param f the file name to write to
163 \param tolerance the accepted tolerance
164
165 \author Philippe Lavoie
166 \date 20 January 1999
167 */
168 template <class T>
drawSubdivisionVRML(const char * f,T tolerance,const Color & col)169 void NurbsSubSurface<T>::drawSubdivisionVRML(const char* f, T tolerance, const Color& col)
170 {
171 ofstream fout ;
172 fout.open(f) ;
173 if(fout)
174 drawSubdivisionVRML(fout,tolerance,col) ;
175 fout.close();
176 }
177
178 /*!
179 \brief perform the subdivision of the NURBS and write the result in a VRML file.
180
181
182 \param f the file name to write to
183 \param tolerance the accepted tolerance
184
185 \author Philippe Lavoie
186 \date 30 April 1999
187 */
188 template <class T>
drawSubdivisionVRML97(const char * f,T tolerance,const Color & col)189 void NurbsSubSurface<T>::drawSubdivisionVRML97(const char* f, T tolerance, const Color& col)
190 {
191 ofstream fout ;
192 fout.open(f) ;
193 if(fout)
194 drawSubdivisionVRML97(fout,tolerance,col) ;
195 fout.close();
196 }
197
198 /*!
199 \brief perform the subdivision of the NURBS and write the result in a PS file.
200
201
202 \param os the ostream to write to
203 \param tolerance the accepted tolerance
204
205 \author Philippe Lavoie
206 \date 20 January 1999
207 */
208 template <class T>
drawSubdivisionPS(ostream & os,T tolerance)209 void NurbsSubSurface<T>::drawSubdivisionPS(ostream &os, T tolerance)
210 {
211 if(render)
212 delete render ;
213 render = new RenderMeshPS<T>(os) ;
214 drawSubdivision(tolerance) ;
215 }
216
217 /*!
218 \brief perform the subdivision of the NURBS and write the result in a VRML file.
219
220 \param os the ostream to write to
221 \param tolerance the accepted tolerance
222
223
224 \author Philippe Lavoie
225 \date 20 January 1999
226 */
227 template <class T>
drawSubdivisionVRML(ostream & os,T tolerance,const Color & col)228 void NurbsSubSurface<T>::drawSubdivisionVRML(ostream &os, T tolerance, const Color& col)
229 {
230 if(render)
231 delete render ;
232 render = new RenderMeshVRML<T>(os,col) ;
233 drawSubdivision(tolerance) ;
234 }
235
236 /*!
237 \brief perform the subdivision of the NURBS and write the result in a VRML file.
238
239 \param os the ostream to write to
240 \param tolerance the accepted tolerance
241
242
243 \author Philippe Lavoie
244 \date 30 April 1999
245 */
246 template <class T>
drawSubdivisionVRML97(ostream & os,T tolerance,const Color & col)247 void NurbsSubSurface<T>::drawSubdivisionVRML97(ostream &os, T tolerance, const Color& col)
248 {
249 if(render)
250 delete render ;
251 render = new RenderMeshVRML97<T>(os,col) ;
252 drawSubdivision(tolerance) ;
253 }
254
255 /*!
256 \brief perform the subdivision of the NURBS and write the result in a VRML file.
257
258 \param os the ostream to write to
259 \param tolerance the accepted tolerance
260
261
262 \author Philippe Lavoie
263 \date 20 January 1999
264 */
265 template <class T>
266 //void NurbsSubSurface<T>::drawSubdivisionPoints(vector<Point_nD<T,3> > &pnts, T tolerance)
drawSubdivisionPoints(BasicArray<Point_nD<T,3>> & pnts,T tolerance)267 void NurbsSubSurface<T>::drawSubdivisionPoints(BasicArray<Point_nD<T,3> >&pnts, T tolerance)
268 {
269 if(render)
270 delete render ;
271 render = new RenderMeshPoints<T>(pnts) ;
272 drawSubdivision(tolerance) ;
273 }
274
275 /*!
276 \brief initialise the subdivision surface
277
278
279 \param tolerance the accepted tolerance
280
281 \author Philippe Lavoie
282 \date 20 January 1999
283 */
284 template <class T>
initSurf()285 void NurbsSubSurface<T>::initSurf()
286 {
287 if(surf)
288 delete surf ;
289 surf = new NurbSurface<T> ;
290 surf->numU = rsurf.ctrlPnts().rows() ;
291 surf->numV = rsurf.ctrlPnts().cols() ;
292 surf->orderU = rsurf.degreeU()+1 ;
293 surf->orderV = rsurf.degreeV()+1 ;
294 surf->kvU = new T[surf->numU + surf->orderU];
295 surf->kvV = new T[surf->numV + surf->orderV];
296 surf->points.resize(surf->numV,surf->numU);
297
298 surf->flatV = FALSE_;
299 surf->flatU = FALSE_;
300 surf->strU0 = FALSE_;
301 surf->strUn = FALSE_;
302 surf->strV0 = FALSE_;
303 surf->strVn = FALSE_;
304
305 surf->render = render ;
306
307 int i;
308 for(i=rsurf.knotU().n()-1;i>=0;--i){
309 surf->kvU[i] = rsurf.knotU()[i] ;
310 }
311 for(i=rsurf.knotV().n()-1;i>=0;--i){
312 surf->kvV[i] = rsurf.knotV()[i] ;
313 }
314 for(i=rsurf.ctrlPnts().rows()-1;i>=0;--i)
315 for(int j=rsurf.ctrlPnts().cols()-1;j>=0;--j)
316 surf->points(j,i) = rsurf.ctrlPnts()(i,j) ;
317 }
318
319 template <class T>
NurbSurface(Matrix<HPoint_nD<T,3>> & s,int du,int dv)320 NurbSurface<T>::NurbSurface(Matrix<HPoint_nD<T,3> >& s, int du, int dv): points(s), orderU(du+1), orderV(dv+1){
321 numU = points.rows() ;
322 numV = points.cols() ;
323 kvU = new T[numU + orderU] ;
324 kvV = new T[numV + orderV] ;
325 render = 0 ;
326 }
327
328 /*!
329 \brief Destructor
330
331 \author Philippe Lavoie
332 \date 20 January 1999
333 */
334 template <class T>
~NurbsSubSurface()335 NurbsSubSurface<T>::~NurbsSubSurface()
336 {
337 if(surf){
338 delete []surf->kvU ;
339 delete []surf->kvV ;
340 }
341 }
342
343 /*!
344 \brief projects from world to screen coordinates
345
346 A Post Script point is the projection of the point from the
347 homogenous space to the 2D paper surface with the axis multiplied
348 by 100 and with an offset of 200.
349
350 There is no perspective projection performed.
351
352 \param worldPt the point in world coordinate
353 \param screenPt the point in the VRML coordinate
354
355 \author Philippe Lavoie
356 \date 20 January 1999 */
357 template <class T>
screenProject(const HPoint_nD<T,3> & worldPt,Point_nD<T,3> & screenPt)358 void RenderMeshPS<T>::screenProject(const HPoint_nD<T,3> &worldPt, Point_nD<T,3> &screenPt )
359 {
360 screenPt.x() = worldPt.x() / worldPt.w() * 100 + 200;
361 screenPt.y() = worldPt.y() / worldPt.w() * 100 + 200;
362 screenPt.z() = worldPt.z() / worldPt.w() * 100 + 200;
363 }
364
365 /*!
366 \brief projects from world to screen coordinates
367
368 In the case of a VRML file, the world and screen coordinate are the same.
369 Except that one is in homogenous space and the other in normal space.
370
371 \param worldPt the point in world coordinate
372 \param screenPt the point in the VRML coordinate
373
374 \author Philippe Lavoie
375 \date 20 January 1999
376 */
377 template <class T>
screenProject(const HPoint_nD<T,3> & worldPt,Point_nD<T,3> & screenPt)378 void RenderMeshVRML<T>::screenProject(const HPoint_nD<T,3> &worldPt, Point_nD<T,3> &screenPt )
379 {
380 screenPt = project(worldPt) ;
381 }
382
383
384 /*!
385 \brief projects from world to screen coordinates
386
387 In the case of a VRML file, the world and screen coordinate are the same.
388 Except that one is in homogenous space and the other in normal space.
389
390 \param worldPt the point in world coordinate
391 \param screenPt the point in the VRML coordinate
392
393 \author Philippe Lavoie
394 \date 30 April 1999
395 */
396 template <class T>
screenProject(const HPoint_nD<T,3> & worldPt,Point_nD<T,3> & screenPt)397 void RenderMeshVRML97<T>::screenProject(const HPoint_nD<T,3> &worldPt, Point_nD<T,3> &screenPt )
398 {
399 screenPt = project(worldPt) ;
400 if(init){
401 p_min = p_max = screenPt ;
402 init = 0 ;
403 }
404 if(screenPt.x()<p_min.x()) p_min.x() = screenPt.x();
405 if(screenPt.y()<p_min.y()) p_min.y() = screenPt.y();
406 if(screenPt.z()<p_min.z()) p_min.z() = screenPt.z();
407 if(screenPt.x()>p_max.x()) p_max.x() = screenPt.x();
408 if(screenPt.y()>p_max.y()) p_max.y() = screenPt.y();
409 if(screenPt.z()>p_max.z()) p_max.z() = screenPt.z();
410 }
411
412
413 /*!
414 \brief projects from world to screen coordinates
415
416 The world and screen coordinate are the same. Except that one is
417 in homogenous space and the other in normal space.
418
419 \param worldPt the point in world coordinate
420 \param screenPt the point in the normal space
421
422 \author Philippe Lavoie
423 \date 20 January 1999
424 */
425 template <class T>
screenProject(const HPoint_nD<T,3> & worldPt,Point_nD<T,3> & screenPt)426 void RenderMeshPoints<T>::screenProject(const HPoint_nD<T,3> &worldPt, Point_nD<T,3> &screenPt )
427 {
428 screenPt = project(worldPt) ;
429 }
430
431
LineTo(ostream & out,short x,short y)432 inline void LineTo( ostream& out, short x, short y )
433 {
434 out << (int)x << " " << (int)y << " lineto\n" ;
435 }
436
MoveTo(ostream & out,short x,short y)437 inline void MoveTo( ostream& out, short x, short y )
438 {
439 out << (int)x << " " << (int)y << " moveto\n" ;
440 }
441
442 /*!
443 \brief write the header of a PS file
444
445 \author Philippe Lavoie
446 \date 20 January 1999
447 */
448 template <class T>
drawHeader()449 void RenderMeshPS<T>::drawHeader(){
450 out << "%!PS-Adobe-2.1\n";
451 out << "%%Title: code5_0.ps (20)\n" ;
452 out << "%%Creator: color_grid_generator\n" ;
453 out << "%%BoundingBox: 0 0 500 500\n" ;
454 out << "%%Pages: 0\n" ;
455 out << "%%EndComments\n" ;
456 out << "0 setlinewidth\n" ;
457 out << "0 0 0 setrgbcolor\n" ;
458 }
459
460 /*!
461 \brief write the footer of a PS file
462
463 \author Philippe Lavoie
464 \date 20 January 1999
465 */
466 template <class T>
drawFooter()467 void RenderMeshPS<T>::drawFooter(){
468 out << "showpage\n%%EOF\n" ;
469 }
470
471
472 /*!
473 \brief Draw a triangle
474
475 \param v0
476 \param v1
477 \param v2
478
479 \author Philippe Lavoie
480 \date 20 January 1999
481 */
482 template <class T>
drawTriangle(const SurfSample<T> & v0,const SurfSample<T> & v1,const SurfSample<T> & v2)483 void RenderMeshPS<T>::drawTriangle( const SurfSample<T> &v0, const SurfSample<T> &v1, const SurfSample<T> & v2 )
484 {
485 out << "newpath\n" ;
486 MoveTo( out, (short) (v0.point.x() * 100 + 200), (short) (v0.point.y() * 100 + 200) );
487 LineTo( out, (short) (v1.point.x() * 100 + 200), (short) (v1.point.y() * 100 + 200) );
488 LineTo( out, (short) (v2.point.x() * 100 + 200), (short) (v2.point.y() * 100 + 200) );
489 LineTo( out, (short) (v0.point.x() * 100 + 200), (short) (v0.point.y() * 100 + 200) );
490 out << "stroke\n" ;
491 }
492
493 /*!
494 \brief Draw a line
495
496 \param v0
497 \param v1
498
499 \author Philippe Lavoie
500 \date 18 may 1999
501 */
502 template <class T>
drawLine(const SurfSample<T> & v0,const SurfSample<T> & v1)503 void RenderMeshPS<T>::drawLine( const SurfSample<T> &v0, const SurfSample<T> &v1)
504 {
505 out << "newpath\n" ;
506 MoveTo(out , (short) (v0.point.x() * 100 + 200), (short) (v0.point.y() * 100 + 200) );
507 LineTo( out, (short) (v1.point.x() * 100 + 200), (short) (v1.point.y() * 100 + 200) );
508 out << "stroke\n" ;
509 }
510
511
512 /*!
513 \brief write the header information for a VRML file
514
515 \author Philippe Lavoie
516 \date 30 April 1999
517 */
518 template <class T>
drawHeader()519 void RenderMeshVRML97<T>::drawHeader()
520 {
521 size = 0 ;
522 out << "#VRML V2.0 utf8\n" ;
523 out << "# Generated by Phil's NURBS library\n" ;
524 out << "\nGroup {\n" ;
525 out << "\n children [\n" ;
526 out << "\tDEF T Transform {\n";
527 out << "\t children [\n" ;
528 out << "\t\tShape {\n" ;
529 out << "\t\t appearance Appearance {\n" ;
530 out << "\t\t\tmaterial Material{ diffuseColor " << float(color.r/255.0) << ' ' << float(color.g/255.0) << ' ' << float(color.b/255.0) << " }\n" ;
531 out << "\t\t }\n" ;
532 out << "\t\t geometry IndexedFaceSet {\n" ;
533 out << "\t\t\tsolid FALSE\n" ;
534 out << "\t\t\tcoord Coordinate {\n" ;
535 out << "\t\t\t point [\n" ;
536 }
537
538 /*!
539 \brief write the footer information for a VRML file
540
541 Write the footer information for a VRML file
542
543 \author Philippe Lavoie
544 \date 30 April 1999
545 */
546 template <class T>
drawFooter()547 void RenderMeshVRML97<T>::drawFooter(){
548
549 out << "\t\t\t ]\n" ; // point [
550 out << "\t\t\t}\n" ; // coord
551
552 out << "\t\t\t coordIndex [\n" ;
553
554 for(int i=0;i<size;++i){
555 out << "\t\t\t\t" << 3*i << ", " << 3*i+1 << ", " << 3*i+2 << ", -1,\n" ;
556 }
557 out << "\t\t\t ]\n" ;
558 out << "\t\t\t}\n" ; // IndexedFaceSet
559 out << "\t\t}\n" ;
560 out << "\t ]\n" ;
561 out << "\t}\n" ;
562 out << " ]\n" ;
563 out << "}\n" ;
564
565 Point_nD<T,3> p_mid((p_max.x()+p_min.x())/T(2),
566 (p_max.y()+p_min.y())/T(2),
567 (p_max.z()+p_min.z())/T(2));
568
569 T x_axis = p_max.x() - p_min.x() ;
570 T y_axis = p_max.y() - p_min.y() ;
571
572 T axis = (x_axis< y_axis) ? y_axis : x_axis ;
573 axis *= T(2) ;
574
575 out << "Viewpoint {\n\t position " << p_mid.x() << ' ' << p_mid.y() << ' ' << p_max.z()+axis << "\n\t description \"top\"\n}\n" ;
576
577
578 out << "NavigationInfo { type \"EXAMINE\" }\n" ;
579
580 }
581
582 /*!
583 \brief write the header information for a VRML file
584
585 \author Philippe Lavoie
586 \date 20 January 1999
587 */
588 template <class T>
drawHeader()589 void RenderMeshVRML<T>::drawHeader()
590 {
591 size = 0 ;
592 out << "#VRML V1.0 ascii\n" ;
593 out << "# Generated by Phil's NURBS library\n" ;
594 out << "\nSeparator {\n" << "\tMaterialBinding { value PER_VERTEX_INDEXED }\n" ;
595 out << "\tMaterial{ ambientColor 0.25 0.25 0.25\n\t\tdiffuseColor " << float(color.r/255.0) << ' ' << float(color.g/255.0) << ' ' << float(color.b/255.0) << " }\n" ;
596 out << "\tCoordinate3 {\n" ;
597 out << "\t\tpoint [\n" ;
598 }
599
600 /*!
601 \brief write the footer information for a VRML file
602
603 Write the footer information for a VRML file
604
605 \author Philippe Lavoie
606 \date 20 January 1999
607 */
608 template <class T>
drawFooter()609 void RenderMeshVRML<T>::drawFooter(){
610
611 out << "\t\t]\n" ; // point [
612 out << "\t}\n" ; // cordinate3
613
614 out << "\tIndexedFaceSet{\n" ;
615 out << "\t\tcoordIndex[\n" ;
616
617 for(int i=0;i<size;++i){
618 out << "\t\t\t" << 3*i << ", " << 3*i+1 << ", " << 3*i+2 << ", -1,\n" ;
619 }
620 out << "\t\t]\n" ;
621 out << "\t}\n" ; // IndexedFaceSet
622
623 out << "}\n" ;
624 }
625
626
627 /*!
628 \brief draws the triangle
629
630 This function draws the triangle points to the ostream.
631
632 \param v0 a corner point of the triangle
633 \param v1 a corner point of the triangle
634 \param v2 a corner point of the triangle
635
636 \author Philippe Lavoie
637 \date 20 January 1999
638 */
639 template <class T>
drawTriangle(const SurfSample<T> & v0,const SurfSample<T> & v1,const SurfSample<T> & v2)640 void RenderMeshVRML<T>::drawTriangle(const SurfSample<T> &v0, const SurfSample<T> &v1, const SurfSample<T> & v2 )
641 {
642 ++size ;
643 out << "\t\t\t" << v0.point.x() << " " << v0.point.y() << " " << v0.point.z() << ",\n" ;
644 out << "\t\t\t" << v1.point.x() << " " << v1.point.y() << " " << v1.point.z() << ",\n" ;
645 out << "\t\t\t" << v2.point.x() << " " << v2.point.y() << " " << v2.point.z() << ",\n" ;
646
647 }
648
649 /*!
650 \brief draws the triangle
651
652 This function draws the triangle points to the ostream.
653
654 \param v0 a corner point of the triangle
655 \param v1 a corner point of the triangle
656 \param v2 a corner point of the triangle
657
658 \author Philippe Lavoie
659 \date 30 April 1999
660 */
661 template <class T>
drawTriangle(const SurfSample<T> & v0,const SurfSample<T> & v1,const SurfSample<T> & v2)662 void RenderMeshVRML97<T>::drawTriangle(const SurfSample<T> &v0, const SurfSample<T> &v1, const SurfSample<T> & v2 )
663 {
664 ++size ;
665 out << "\t\t\t\t" << v0.point.x() << " " << v0.point.y() << " " << v0.point.z() << ",\n" ;
666 out << "\t\t\t\t" << v1.point.x() << " " << v1.point.y() << " " << v1.point.z() << ",\n" ;
667 out << "\t\t\t\t" << v2.point.x() << " " << v2.point.y() << " " << v2.point.z() << ",\n" ;
668
669 }
670
671 /*!
672 \brief write the header information for a mesh file
673
674 \author Philippe Lavoie
675 \date 20 January 1999
676 */
677 template <class T>
drawHeader()678 void RenderMeshPoints<T>::drawHeader()
679 {
680 points.clear();
681 }
682
683 /*!
684 \brief empty function
685
686
687 \author Philippe Lavoie
688 \date 20 January 1999
689 */
690 template <class T>
drawFooter()691 void RenderMeshPoints<T>::drawFooter(){
692 }
693
694
695 /*!
696 \brief draws the triangle
697
698 Adds the triangle points to the point vector.
699
700 \param v0 a corner point of the triangle
701 \param v1 a corner point of the triangle
702 \param v2 a corner point of the triangle
703
704 \author Philippe Lavoie
705 \date 20 January 1999
706 */
707 template <class T>
drawTriangle(const SurfSample<T> & v0,const SurfSample<T> & v1,const SurfSample<T> & v2)708 void RenderMeshPoints<T>::drawTriangle(const SurfSample<T> &v0, const SurfSample<T> &v1, const SurfSample<T> & v2 )
709 {
710 // naive method
711 points.push_back(v0.point);
712 points.push_back(v1.point);
713 points.push_back(v2.point);
714 }
715
716 template <class T>
717 int
FindBreakPoint(T u,T * kv,int m,int k)718 FindBreakPoint( T u, T * kv, int m, int k )
719 {
720 int i;
721
722 if (u == kv[m+1]) /* Special case for closed interval */
723 return m;
724
725 i = m + k;
726 while ((u < kv[i]) && (i > 0))
727 i--;
728 return( i );
729 }
730
731 /*
732 * Compute Bi,k(u), for i = 0..k.
733 * u is the parameter of the spline to find the basis functions for
734 * brkPoint is the start of the knot interval ("segment")
735 * kv is the knot vector
736 * k is the order of the curve
737 * bvals is the array of returned basis values.
738 *
739 * (From Bartels, Beatty & Barsky, p.387)
740 */
741
742 template <class T> void
BasisFunctions(T u,int brkPoint,T * kv,int k,T * bvals)743 BasisFunctions( T u, int brkPoint, T * kv, int k, T * bvals )
744 {
745 int r, s, i;
746 T omega;
747
748 bvals[0] = 1.0;
749 for (r = 2; r <= k; r++)
750 {
751 i = brkPoint - r + 1;
752 bvals[r - 1] = 0.0;
753 for (s = r-2; s >= 0; s--)
754 {
755 i++;
756 if (i < 0)
757 omega = 0;
758 else
759 omega = (u - kv[i]) / (kv[i + r - 1] - kv[i]);
760 bvals[s + 1] = bvals[s + 1] + (1 - omega) * bvals[s];
761 bvals[s] = omega * bvals[s];
762 }
763 }
764 }
765
766 /*
767 * Compute derivatives of the basis functions Bi,k(u)'
768 */
769 template <class T> void
BasisDerivatives(T u,int brkPoint,T * kv,int k,T * dvals)770 BasisDerivatives( T u, int brkPoint, T * kv, int k, T * dvals )
771 {
772 register int s, i;
773 T omega, knotScale;
774
775 BasisFunctions( u, brkPoint, kv, k - 1, dvals );
776
777 dvals[k-1] = 0.0; /* BasisFunctions misses this */
778
779 knotScale = kv[brkPoint + 1] - kv[brkPoint];
780
781 i = brkPoint - k + 1;
782 for (s = k - 2; s >= 0; s--)
783 {
784 i++;
785 omega = knotScale * ((T)(k-1)) / (kv[i+k-1] - kv[i]);
786 dvals[s + 1] += -omega * dvals[s];
787 dvals[s] *= omega;
788 }
789 }
790
791 /*
792 * Calculate a point p on NurbSurface<T> n at a specific u, v using the tensor product.
793 * If utan and vtan are not nil, compute the u and v tangents as well.
794 *
795 * Note the valid parameter range for u and v is
796 * (kvU[orderU] <= u < kvU[numU), (kvV[orderV] <= v < kvV[numV])
797 */
798
799 template <class T>
800 void
CalcPoint(T u,T v,NurbSurface<T> * n,Point_nD<T,3> * p,Point_nD<T,3> * utan,Point_nD<T,3> * vtan)801 CalcPoint(T u, T v, NurbSurface<T> * n,
802 Point_nD<T,3> * p, Point_nD<T,3> * utan, Point_nD<T,3> * vtan)
803 {
804 int i, j, ri, rj;
805 HPoint_nD<T,3> * cp;
806 T tmp;
807 T wsqrdiv;
808 int ubrkPoint, ufirst;
809 T bu[MAXORDER], buprime[MAXORDER];
810 int vbrkPoint, vfirst;
811 T bv[MAXORDER], bvprime[MAXORDER];
812 HPoint_nD<T,3> r, rutan, rvtan;
813
814 r.x() = 0.0;
815 r.y() = 0.0;
816 r.z() = 0.0;
817 r.w() = 0.0;
818
819 rutan = r;
820 rvtan = r;
821
822 /* Evaluate non-uniform basis functions (and derivatives) */
823
824 ubrkPoint = FindBreakPoint( u, n->kvU, n->numU-1, n->orderU );
825 ufirst = ubrkPoint - n->orderU + 1;
826 BasisFunctions( u, ubrkPoint, n->kvU, n->orderU, bu );
827 if (utan)
828 BasisDerivatives( u, ubrkPoint, n->kvU, n->orderU, buprime );
829
830 vbrkPoint = FindBreakPoint( v, n->kvV, n->numV-1, n->orderV );
831 vfirst = vbrkPoint - n->orderV + 1;
832 BasisFunctions( v, vbrkPoint, n->kvV, n->orderV, bv );
833 if (vtan)
834 BasisDerivatives( v, vbrkPoint, n->kvV, n->orderV, bvprime );
835
836 /* Weight control points against the basis functions */
837
838 for (i = 0; i < n->orderV; i++)
839 for (j = 0; j < n->orderU; j++)
840 {
841 ri = n->orderV - 1 - i;
842 rj = n->orderU - 1 - j;
843
844 tmp = bu[rj] * bv[ri];
845 cp = &( n->points(i+vfirst,j+ufirst) );
846 r.x() += cp->x() * tmp;
847 r.y() += cp->y() * tmp;
848 r.z() += cp->z() * tmp;
849 r.w() += cp->w() * tmp;
850
851 if (utan)
852 {
853 tmp = buprime[rj] * bv[ri];
854 rutan.x() += cp->x() * tmp;
855 rutan.y() += cp->y() * tmp;
856 rutan.z() += cp->z() * tmp;
857 rutan.w() += cp->w() * tmp;
858 }
859 if (vtan)
860 {
861 tmp = bu[rj] * bvprime[ri];
862 rvtan.x() += cp->x() * tmp;
863 rvtan.y() += cp->y() * tmp;
864 rvtan.z() += cp->z() * tmp;
865 rvtan.w() += cp->w() * tmp;
866 }
867 }
868
869 /* Project tangents, using the quotient rule for differentiation */
870
871 wsqrdiv = 1.0 / (r.w() * r.w());
872 if (utan)
873 {
874 utan->x() = (r.w() * rutan.x() - rutan.w() * r.x()) * wsqrdiv;
875 utan->y() = (r.w() * rutan.y() - rutan.w() * r.y()) * wsqrdiv;
876 utan->z() = (r.w() * rutan.z() - rutan.w() * r.z()) * wsqrdiv;
877 }
878 if (vtan)
879 {
880 vtan->x() = (r.w() * rvtan.x() - rvtan.w() * r.x()) * wsqrdiv;
881 vtan->y() = (r.w() * rvtan.y() - rvtan.w() * r.y()) * wsqrdiv;
882 vtan->z() = (r.w() * rvtan.z() - rvtan.w() * r.z()) * wsqrdiv;
883 }
884
885 p->x() = r.x() / r.w();
886 p->y() = r.y() / r.w();
887 p->z() = r.z() / r.w();
888 }
889
890 /*
891 * Draw a mesh of points by evaluating the surface at evenly spaced
892 * points.
893 */
894 template <class T>
895 void
DrawEvaluation(NurbSurface<T> * n)896 DrawEvaluation( NurbSurface<T> * n )
897 {
898 Point_nD<T,3> p, utan, vtan;
899 register int i, j;
900 register T u, v, d;
901 SurfSample<T> ** pts ;
902
903 int Granularity = 10; /* Controls the number of steps in u and v */
904
905 /* Allocate storage for the grid of points generated */
906
907 CHECK( pts = new SurfSample<T>* [Granularity+1]);
908 CHECK( pts[0] = new SurfSample<T>[(Granularity+1)*(Granularity+1)]);
909
910 for (i = 1; i <= Granularity; i++)
911 pts[i] = &(pts[0][(Granularity+1) * i]);
912
913 /* Compute points on curve */
914
915 for (i = 0; i <= Granularity; i++)
916 {
917 v = ((T) i / (T) Granularity)
918 * (n->kvV[n->numV] - n->kvV[n->orderV-1])
919 + n->kvV[n->orderV-1];
920
921 for (j = 0; j <= Granularity; j++)
922 {
923 u = ((T) j / (T) Granularity)
924 * (n->kvU[n->numU] - n->kvU[n->orderU-1])
925 + n->kvU[n->orderU-1];
926
927 CalcPoint( u, v, n, &(pts[i][j].point), &utan, &vtan );
928 p = crossProduct(utan,vtan) ; //(void) V3Cross( &utan, &vtan, &p );
929 d = norm(p) ; // d = V3Length( &p );
930 if (d != 0.0)
931 {
932 p.x() /= d;
933 p.y() /= d;
934 p.z() /= d;
935 }
936 else
937 {
938 p.x() = 0;
939 p.y() = 0;
940 p.z() = 0;
941 }
942 pts[i][j].normLen = d;
943 pts[i][j].normal = p;
944 pts[i][j].u = u;
945 pts[i][j].v = v;
946 }
947 }
948
949 /* Draw the grid */
950
951 for (i = 0; i < Granularity; i++)
952 for (j = 0; j < Granularity; j++)
953 {
954 n->render->drawTriangle( pts[i][j], pts[i+1][j+1], pts[i+1][j] );
955 n->render->drawTriangle( pts[i][j], pts[i][j+1], pts[i+1][j+1] );
956 }
957
958 delete [] pts[0];
959 delete [] pts ;
960 }
961
962 /*
963 * NurbRefine.c - Given a refined knot vector, add control points to a surface.
964 *
965 * John Peterson
966 */
967
968
969 /*
970 * Given the original knot vector ukv, and a new knotvector vkv, compute
971 * the "alpha matrix" used to generate the corresponding new control points.
972 * This routines allocates the alpha matrix if it isn't allocated already.
973 *
974 * This is from Bartels, Beatty & Barsky, p. 407
975 */
976 template <class T> void
CalcAlpha(T * ukv,T * wkv,int m,int n,int k,T *** alpha)977 CalcAlpha( T * ukv, T * wkv, int m, int n, int k, T *** alpha )
978 {
979 int i, j;
980 int brkPoint, r, rm1, last, s;
981 T omega;
982 T aval[MAXORDER];
983
984 if (! *alpha) /* Must allocate alpha */
985 {
986 CHECK( *alpha = new T* [k+1]);
987 for (i = 0; i <= k; i++)
988 CHECK( (*alpha)[i] = new T[(m + n + 1)]);
989 }
990
991 for (j = 0; j <= m + n; j++)
992 {
993 brkPoint = FindBreakPoint( wkv[j], ukv, m, k );
994 aval[0] = 1.0;
995 for (r = 2; r <= k; r++)
996 {
997 rm1 = r - 1;
998 last = minimum( rm1, brkPoint );
999 i = brkPoint - last;
1000 if (last < rm1)
1001 aval[last] = aval[last] * (wkv[j + r - 1] - ukv[i])
1002 / (ukv[i + r - 1] - ukv[i]);
1003 else
1004 aval[last] = 0.0;
1005
1006 for (s = last - 1; s >= 0; s-- )
1007 {
1008 i++;
1009 omega = (wkv[j + r - 1] - ukv[i]) / (ukv[i + r - 1] - ukv[i]);
1010 aval[s + 1] = aval[s+1] + (1 - omega) * aval[s];
1011 aval[s] = omega * aval[s];
1012 }
1013 }
1014 last = minimum( k - 1, brkPoint );
1015 for (i = 0; i <= k; i++)
1016 (*alpha)[i][j] = 0.0;
1017 for (s = 0; s <= last; s++)
1018 (*alpha)[last - s][j] = aval[s];
1019 }
1020 }
1021
1022 /*
1023 * Apply the alpha matrix computed above to the rows (or columns)
1024 * of the surface. If dirflag is TRUE_ do the U's (row), else do V's (col).
1025 */
1026 template <class T>
1027 void
RefineSurface(NurbSurface<T> * src,NurbSurface<T> * dest,BOOL dirflag)1028 RefineSurface( NurbSurface<T> * src, NurbSurface<T> * dest, BOOL dirflag )
1029 {
1030 int i, j, out;
1031 HPoint_nD<T,3> * dp, * sp;
1032 int i1, brkPoint, maxj, maxout;
1033 register T tmp;
1034 T ** alpha = 0;
1035
1036 // Compute the alpha matrix and indexing variables for the requested direction
1037
1038 if (dirflag)
1039 {
1040 CalcAlpha( src->kvU, dest->kvU, src->numU - 1, dest->numU - src->numU,
1041 src->orderU, &alpha );
1042 maxj = dest->numU;
1043 maxout = src->numV;
1044 }
1045 else
1046 {
1047 CalcAlpha( src->kvV, dest->kvV, src->numV - 1, dest->numV - src->numV,
1048 src->orderV, &alpha );
1049 maxj = dest->numV;
1050 maxout = dest->numU;
1051 }
1052
1053 /* Apply the alpha matrix to the original control points, generating new ones */
1054
1055 for (out = 0; out < maxout; out++)
1056 for (j = 0; j < maxj; j++)
1057 {
1058 if (dirflag)
1059 {
1060 dp = &(dest->points(out,j));
1061 brkPoint = FindBreakPoint( dest->kvU[j], src->kvU,
1062 src->numU-1, src->orderU );
1063 i1 = maximum( brkPoint - src->orderU + 1, 0 );
1064 sp = &(src->points(out,i1));
1065 } else {
1066 dp = &(dest->points(j,out));
1067 brkPoint = FindBreakPoint( dest->kvV[j], src->kvV,
1068 src->numV-1, src->orderV );
1069 i1 = maximum( brkPoint - src->orderV + 1, 0 );
1070 sp = &(src->points(i1,out));
1071 }
1072 dp->x() = 0.0;
1073 dp->y() = 0.0;
1074 dp->z() = 0.0;
1075 dp->w() = 0.0;
1076 for (i = i1; i <= brkPoint; i++)
1077 {
1078 tmp = alpha[i - i1][j];
1079 sp = (dirflag ? &(src->points(out,i)) : &(src->points(i,out)) );
1080 dp->x() += tmp * sp->x();
1081 dp->y() += tmp * sp->y();
1082 dp->z() += tmp * sp->z();
1083 dp->w() += tmp * sp->w();
1084 }
1085 }
1086
1087 /* Free up the alpha matrix */
1088 for (i = 0; i <= (dirflag ? src->orderU : src->orderV); i++)
1089 delete [] alpha[i] ;
1090 delete [] alpha ;
1091 }
1092
1093 /*
1094 * NurbUtils.c - Code for Allocating, freeing, & copying NURB surfaces.
1095 *
1096 * John Peterson
1097 */
1098
1099
1100 /*
1101 * Allocate memory for a NURB (assumes numU, numV, orderU
1102 * and orderV have been set). If ukv or vkv are not NIL, they
1103 * are assumed to be pointers to valid knot vectors.
1104 */
1105
1106 template <class T>
1107 void
AllocNurb(NurbSurface<T> * n,T * ukv,T * vkv)1108 AllocNurb( NurbSurface<T> * n, T * ukv, T * vkv )
1109 {
1110 int i;
1111
1112 if (! ukv)
1113 n->kvU = new T[n->numU + n->orderU] ;
1114 else
1115 n->kvU = ukv;
1116 if (! vkv)
1117 n->kvV = new T[n->numV + n->orderV];
1118 else
1119 n->kvV = vkv;
1120
1121 n->points.resize(n->numV,n->numU) ;
1122 }
1123
1124 /*
1125 * Release storage for a patch
1126 */
1127
1128 template <class T>
1129 void
FreeNurb(NurbSurface<T> * n)1130 FreeNurb( NurbSurface<T> * n )
1131 {
1132 int i;
1133
1134 if (n->kvU) delete [] n->kvU ;
1135 n->kvU = 0;
1136 if (n->kvV) delete [] n->kvV ;
1137 n->kvV = 0;
1138 delete n ;
1139 n = 0 ;
1140 // Don't touch render, it might still be used.
1141 }
1142
1143 /*
1144 * Clone a nurb (deep copy)
1145 */
1146
1147 template <class T>
1148 void
CloneNurb(NurbSurface<T> * src,NurbSurface<T> * dst)1149 CloneNurb( NurbSurface<T> * src, NurbSurface<T> * dst )
1150 {
1151 int i, j;
1152 T * srcp, *dstp;
1153
1154 //*dst = *src; /* Copy fields that don't change */
1155 dst->numU = src->numU ;
1156 dst->numV = src->numV ;
1157 dst->orderU = src->orderU ;
1158 dst->orderV = src->orderV ;
1159
1160 dst->strU0 = src->strU0 ;
1161 dst->strUn = src->strUn ;
1162 dst->strV0 = src->strV0 ;
1163 dst->strVn = src->strVn ;
1164
1165 dst->kvU = 0;
1166 dst->kvV = 0; /* So they get allocated */
1167 dst->points = 0;
1168
1169 AllocNurb( dst, (T*)0, (T*)0 );
1170
1171 /* Copy kv's */
1172 srcp = src->kvU;
1173 dstp = dst->kvU;
1174 for (i = 0; i < src->numU + src->orderU; i++)
1175 *dstp++ = *srcp++;
1176
1177 srcp = src->kvV;
1178 dstp = dst->kvV;
1179 for (i = 0; i < src->numV + src->orderV; i++)
1180 *dstp++ = *srcp++;
1181
1182 /* Copy control points */
1183 for (i = 0; i < src->numV; i++)
1184 for (j = 0; j < src->numU; j++)
1185 dst->points(i,j) = src->points(i,j);
1186 }
1187
1188 /*
1189 * NurbSubdiv.c - Perform adaptive subdivision on a NURB surface.
1190 *
1191 * John Peterson
1192 *
1193 */
1194
1195 #define DIVPT( p, dn ) { ((p).x()) /= (dn); ((p).y()) /= (dn); ((p).z()) /= (dn); }
1196
1197 #define maxV(surf) ((surf)->numV-1)
1198 #define maxU(surf) ((surf)->numU-1)
1199
1200 /*
1201 * Split a knot vector at the center, by adding multiplicity k knots near
1202 * the middle of the parameter range. Tries to start with an existing knot,
1203 * but will add a new knot value if there's nothing in "the middle" (e.g.,
1204 * a Bezier curve).
1205 */
1206 template <class T> int
SplitKV(T * srckv,T ** destkv,int * splitPt,int m,int k)1207 SplitKV( T * srckv,
1208 T ** destkv,
1209 int * splitPt, /* Where the knot interval is split */
1210 int m, int k )
1211 {
1212 int i, last;
1213 int middex, extra, same; /* "middex" ==> "middle index" */
1214 T midVal;
1215
1216 extra = 0;
1217 last = (m + k);
1218
1219 middex = last / 2;
1220 midVal = srckv[middex];
1221
1222 /* Search forward and backward to see if multiple knot is already there */
1223
1224 i = middex+1;
1225 same = 1;
1226 while ((i < last) && (srckv[i] == midVal)) {
1227 i++;
1228 same++;
1229 }
1230
1231 i = middex-1;
1232 while ((i > 0) && (srckv[i] == midVal)) {
1233 i--;
1234 middex--; /* middex is start of multiple knot */
1235 same++;
1236 }
1237
1238 if (i <= 0) /* No knot in middle, must create it */
1239 {
1240 midVal = (srckv[0] + srckv[last]) / 2.0;
1241 middex = last / 2;
1242 while (srckv[middex + 1] < midVal)
1243 middex++;
1244 same = 0;
1245 }
1246
1247 extra = k - same;
1248 *destkv = new T[m+k+extra+1];
1249
1250 if (same < k) /* Must add knots */
1251 {
1252 for (i = 0; i <= middex; i++)
1253 (*destkv)[i] = srckv[i];
1254
1255 for (i = middex+1; i <= middex+extra; i++)
1256 (*destkv)[i] = midVal;
1257
1258 for (i = middex + k - same + 1; i <= m + k + extra; i++)
1259 (*destkv)[i] = srckv[i - extra];
1260 }
1261 else
1262 {
1263 for (i = 0; i <= m + k; i++)
1264 (*destkv)[i] = srckv[i];
1265 }
1266
1267 *splitPt = (extra < k) ? middex - 1 : middex;
1268 return( extra );
1269 }
1270
1271 /*
1272 * Given a line defined by firstPt and lastPt, project midPt onto
1273 * that line. Used for fixing "cracks".
1274 */
1275 template <class T> void
ProjectToLine(Point_nD<T,3> * firstPt,Point_nD<T,3> * lastPt,Point_nD<T,3> * midPt)1276 ProjectToLine( Point_nD<T,3> * firstPt, Point_nD<T,3> * lastPt, Point_nD<T,3> * midPt )
1277 {
1278 Point_nD<T,3> base, v0, vm;
1279 T fraction, denom;
1280
1281 base = *firstPt;
1282
1283 v0 = *lastPt - base ; // (void) V3Sub( lastPt, &base, &v0 );
1284 vm = *midPt - base ; // (void) V3Sub( midPt, &base, &vm );
1285
1286 denom = norm2(v0) ; // V3SquaredLength( &v0 );
1287 //fraction = (denom == 0.0) ? 0.0 : (V3Dot( &v0, &vm ) / denom);
1288 fraction = (denom == 0.0) ? 0.0 : (v0*vm ) / denom;
1289
1290 midPt->x() = base.x() + fraction * v0.x();
1291 midPt->y() = base.y() + fraction * v0.y();
1292 midPt->z() = base.z() + fraction * v0.z();
1293 }
1294
1295 /*
1296 * If a normal has collapsed to zero (normLen == 0.0) then try
1297 * and fix it by looking at its neighbors. If all the neighbors
1298 * are sick, then re-compute them from the plane they form.
1299 * If that fails too, then we give up...
1300 */
1301 template <class T> void
FixNormals(SurfSample<T> * s0,SurfSample<T> * s1,SurfSample<T> * s2)1302 FixNormals( SurfSample<T> * s0, SurfSample<T> * s1, SurfSample<T> * s2 )
1303 {
1304 BOOL goodnorm;
1305 int i, j, ok;
1306 T dist;
1307 SurfSample<T> * V[3];
1308 Point_nD<T,3> normal;
1309
1310 V[0] = s0; V[1] = s1; V[2] = s2;
1311
1312 /* Find a reasonable normalal */
1313 for (ok = 0, goodnorm = FALSE_;
1314 (ok < 3L) && !(goodnorm = (V[ok]->normLen > 0.0)); ok++);
1315
1316 if (! goodnorm) /* All provided normals are zilch, try and invent one */
1317 {
1318 normal.x() = 0.0; normal.y() = 0.0; normal.z() = 0.0;
1319
1320 for (i = 0; i < 3L; i++)
1321 {
1322 j = (i + 1) % 3L;
1323 normal.x() += (V[i]->point.y() - V[j]->point.y()) * (V[i]->point.z() + V[j]->point.z());
1324 normal.y() += (V[i]->point.z() - V[j]->point.z()) * (V[i]->point.x() + V[j]->point.x());
1325 normal.z() += (V[i]->point.x() - V[j]->point.x()) * (V[i]->point.y() + V[j]->point.y());
1326 }
1327 //dist = V3Length( &norm );
1328 dist = norm(normal) ;
1329 if (dist == 0.0)
1330 return; /* This sucker's hopeless... */
1331
1332 DIVPT( normal, dist );
1333
1334 for (i = 0; i < 3; i++)
1335 {
1336 V[i]->normal = normal;
1337 V[i]->normLen = dist;
1338 }
1339 }
1340 else /* Replace a sick normal with a healthy one nearby */
1341 {
1342 for (i = 0; i < 3; i++)
1343 if ((i != ok) && (V[i]->normLen == 0.0))
1344 V[i]->normal = V[ok]->normal;
1345 }
1346 return;
1347 }
1348
1349 /*
1350 * Normalize the normal in a sample. If it's degenerate,
1351 * flag it as such by setting the normLen to 0.0
1352 */
1353 template <class T> void
AdjustNormal(SurfSample<T> * samp)1354 AdjustNormal( SurfSample<T> * samp )
1355 {
1356 // If it's not degenerate, do the normalization now */
1357 samp->normLen = norm(samp->normal) ; // V3Length( &(samp->normal) );
1358
1359 if (samp->normLen < samp->epsilon )
1360 samp->normLen = 0.0;
1361 else
1362 DIVPT( (samp->normal), samp->normLen );
1363 }
1364
1365 /*
1366 * Compute the normal of a corner point of a mesh. The
1367 * base is the value of the point at the corner, indU and indV
1368 * are the mesh indices of that point (either 0 or numU|numV).
1369 */
1370 template <class T> void
GetNormal(NurbSurface<T> * n,int indV,int indU)1371 GetNormal( NurbSurface<T> * n, int indV, int indU )
1372 {
1373 Point_nD<T,3> tmpL, tmpR; /* "Left" and "Right" of the base point */
1374 SurfSample<T> * crnr;
1375
1376 if ( (indU && indV) || ((! indU) && (!indV)) )
1377 {
1378 if (indU)
1379 crnr = &(n->cnn);
1380 else
1381 crnr = &(n->c00);
1382 DIVW( &(n->points(indV,(indU ? (indU-1) : 1))), &tmpL );
1383 DIVW( &(n->points((indV ? (indV-1) : 1),indU)), &tmpR );
1384 }
1385 else
1386 {
1387 if (indU)
1388 crnr = &(n->c0n);
1389 else
1390 crnr = &(n->cn0);
1391 DIVW( &(n->points(indV,(indU ? (indU-1) : 1))), &tmpR );
1392 DIVW( &(n->points((indV ? (indV-1) : 1),indU)), &tmpL );
1393 }
1394
1395 tmpL -= crnr->point ; //(void) V3Sub( &tmpL, &(crnr->point), &tmpL );
1396 tmpR -= crnr->point ;//(void) V3Sub( &tmpR, &(crnr->point), &tmpR );
1397 crnr->normal = crossProduct(tmpL,tmpR); //(void) V3Cross( &tmpL, &tmpR, &(crnr->normal) );
1398 AdjustNormal( crnr );
1399 }
1400
1401 /*
1402 * Build the new corners in the two new surfaces, computing both
1403 * point on the surface aint with the normal. Prevent cracks that may occur.
1404 */
1405 template <class T> void
MakeNewCorners(NurbSurface<T> * parent,NurbSurface<T> * kid0,NurbSurface<T> * kid1,BOOL dirflag)1406 MakeNewCorners( NurbSurface<T> * parent,
1407 NurbSurface<T> * kid0,
1408 NurbSurface<T> * kid1,
1409 BOOL dirflag )
1410 {
1411 DIVW( &(kid0->points(maxV(kid0),maxU(kid0))), &(kid0->cnn.point) );
1412 GetNormal( kid0, maxV(kid0), maxU(kid0) );
1413
1414 if (dirflag)
1415 {
1416 kid0->strUn = FALSE_; /* Must re-test new edge straightness */
1417
1418 DIVW( &(kid0->points(0,maxU(kid0))), &(kid0->c0n.point) );
1419 GetNormal( kid0, 0, maxU(kid0) );
1420 /*
1421 * Normals must be re-calculated for kid1 in case the surface
1422 * was split at a c1 (or c0!) discontinutiy
1423 */
1424 kid1->c00.point = kid0->c0n.point;
1425 GetNormal( kid1, 0, 0 );
1426 kid1->cn0.point = kid0->cnn.point;
1427 GetNormal( kid1, maxV(kid1), 0 );
1428
1429 /*
1430 * Prevent cracks from forming by forcing the points on the seam to
1431 * lie aint any straight edges. (Must do this BEFORE finding normals)
1432 */
1433 if (parent->strV0)
1434 ProjectToLine( &(parent->c00.point),
1435 &(parent->c0n.point),
1436 &(kid0->c0n.point) );
1437 if (parent->strVn)
1438 ProjectToLine( &(parent->cn0.point),
1439 &(parent->cnn.point),
1440 &(kid0->cnn.point) );
1441
1442 kid1->c00.point = kid0->c0n.point;
1443 kid1->cn0.point = kid0->cnn.point;
1444 kid1->strU0 = FALSE_;
1445 }
1446 else
1447 {
1448 kid0->strVn = FALSE_;
1449
1450 DIVW( &(kid0->points(maxV(kid0),0)), &(kid0->cn0.point) );
1451 GetNormal( kid0, maxV(kid0), 0 );
1452 kid1->c00.point = kid0->cn0.point;
1453 GetNormal( kid1, 0, 0 );
1454 kid1->c0n.point = kid0->cnn.point;
1455 GetNormal( kid1, 0, maxU(kid1) );
1456
1457 if (parent->strU0)
1458 ProjectToLine( &(parent->c00.point),
1459 &(parent->cn0.point),
1460 &(kid0->cn0.point) );
1461 if (parent->strUn)
1462 ProjectToLine( &(parent->c0n.point),
1463 &(parent->cnn.point),
1464 &(kid0->cnn.point) );
1465
1466 kid1->c00.point = kid0->cn0.point;
1467 kid1->c0n.point = kid0->cnn.point;
1468 kid1->strV0 = FALSE_;
1469 }
1470 }
1471
1472 /*
1473 * Split a surface into two halves. First inserts multiplicity k knots
1474 * in the center of the parametric range. After refinement, the two
1475 * resulting surfaces are copied into separate data structures. If the
1476 * parent surface had straight edges, the points of the children are
1477 * projected onto those edges.
1478 */
1479 template <class T> void
SplitSurface(NurbSurface<T> * parent,NurbSurface<T> * kid0,NurbSurface<T> * kid1,BOOL dirflag)1480 SplitSurface( NurbSurface<T> * parent,
1481 NurbSurface<T> * kid0, NurbSurface<T> * kid1,
1482 BOOL dirflag ) /* If TRUE_ subdivided in U, else in V */
1483 {
1484 NurbSurface<T> *tmp;
1485 T * newkv;
1486 int i, j, splitPt;
1487
1488 tmp = new NurbSurface<T> ;
1489
1490 //
1491 // Add a multiplicty k knot to the knot vector in the direction
1492 // specified by dirflag, and refine the surface. This creates two
1493 // adjacent surfaces with c0 discontinuity at the seam.
1494 //
1495
1496 //tmp = *parent; // Copy order, # of points, etc.
1497 tmp->numU = parent->numU ;
1498 tmp->numV = parent->numV ;
1499 tmp->orderU = parent->orderU ;
1500 tmp->orderV = parent->orderV ;
1501
1502 tmp->strU0 = parent->strU0 ;
1503 tmp->strUn = parent->strUn ;
1504 tmp->strV0 = parent->strV0 ;
1505 tmp->strVn = parent->strVn ;
1506
1507 tmp->render = parent->render ;
1508
1509 if (dirflag)
1510 {
1511 tmp->numU = parent->numU + SplitKV( parent->kvU,
1512 &newkv,
1513 &splitPt,
1514 maxU(parent),
1515 parent->orderU );
1516 AllocNurb( tmp, newkv, (T*)0 );
1517 for (i = 0; i < tmp->numV + tmp->orderV; i++)
1518 tmp->kvV[i] = parent->kvV[i];
1519 }
1520 else
1521 {
1522 tmp->numV = parent->numV + SplitKV( parent->kvV,
1523 &newkv,
1524 &splitPt,
1525 maxV(parent),
1526 parent->orderV );
1527 AllocNurb( tmp, (T*)0, newkv );
1528 for (i = 0; i < tmp->numU + tmp->orderU; i++)
1529 tmp->kvU[i] = parent->kvU[i];
1530 }
1531 RefineSurface( parent, tmp, dirflag );
1532
1533 //
1534 // Build the two child surfaces, and copy the data from the refined
1535 // version of the parent (tmp) into the two children
1536 //
1537
1538 // First half
1539
1540 // *kid0 = *parent; // copy various edge flags and orders
1541 kid0->orderU = parent->orderU ;
1542 kid0->orderV = parent->orderV ;
1543
1544 kid0->strU0 = parent->strU0 ;
1545 kid0->strUn = parent->strUn ;
1546 kid0->strV0 = parent->strV0 ;
1547 kid0->strVn = parent->strVn ;
1548
1549 kid0->c00 = parent->c00 ;
1550 kid0->c0n = parent->c0n ;
1551 kid0->cn0 = parent->cn0 ;
1552 kid0->cnn = parent->cnn ;
1553
1554 kid0->render = parent->render ;
1555
1556 kid0->numU = dirflag ? splitPt+1 : parent->numU;
1557 kid0->numV = dirflag ? parent->numV : splitPt+1;
1558 kid0->kvU = kid0->kvV = 0;
1559 AllocNurb( kid0, (T*)0, (T*)0 );
1560
1561 for (i = 0; i < kid0->numV; i++) // Copy the point and kv data
1562 for (j = 0; j < kid0->numU; j++)
1563 kid0->points(i,j) = tmp->points(i,j) ;
1564 for (i = 0; i < kid0->orderU + kid0->numU; i++)
1565 kid0->kvU[i] = tmp->kvU[i];
1566 for (i = 0; i < kid0->orderV + kid0->numV; i++)
1567 kid0->kvV[i] = tmp->kvV[i];
1568
1569 // Second half
1570
1571 splitPt++;
1572 //*kid1 = *parent;
1573 kid1->orderU = parent->orderU ;
1574 kid1->orderV = parent->orderV ;
1575
1576 kid1->strU0 = parent->strU0 ;
1577 kid1->strUn = parent->strUn ;
1578 kid1->strV0 = parent->strV0 ;
1579 kid1->strVn = parent->strVn ;
1580
1581 kid1->c00 = parent->c00 ;
1582 kid1->c0n = parent->c0n ;
1583 kid1->cn0 = parent->cn0 ;
1584 kid1->cnn = parent->cnn ;
1585
1586 kid1->render = parent->render ;
1587
1588 kid1->numU = dirflag ? tmp->numU - splitPt : parent->numU;
1589 kid1->numV = dirflag ? parent->numV : tmp->numV - splitPt;
1590 kid1->kvU = kid1->kvV = 0;
1591 AllocNurb( kid1, (T*)0, (T*)0 );
1592
1593 for (i = 0; i < kid1->numV; i++) // Copy the point and kv data
1594 for (j = 0; j < kid1->numU; j++)
1595 kid1->points(i,j)
1596 = tmp->points(dirflag ? i: (i + splitPt) ,dirflag ? (j + splitPt) : j);
1597 for (i = 0; i < kid1->orderU + kid1->numU; i++)
1598 kid1->kvU[i] = tmp->kvU[dirflag ? (i + splitPt) : i];
1599 for (i = 0; i < kid1->orderV + kid1->numV; i++)
1600 kid1->kvV[i] = tmp->kvV[dirflag ? i : (i + splitPt)];
1601
1602 // Construct new corners on the boundry between the two kids
1603 MakeNewCorners( parent, kid0, kid1, dirflag );
1604
1605 FreeNurb( tmp ); // Get rid of refined parent
1606
1607 }
1608
1609 /*
1610 * Test if a particular row or column of control points in a mesh
1611 * is "straight" with respect to a particular tolerance. Returns TRUE_
1612 * if it is.
1613 */
1614
1615 #define GETPT( i ) (( dirflag ? (n->points(crvInd,i)) : (n->points(i,crvInd)) ))
1616
1617 #define EPSILON n->epsilon
1618
1619 template <class T> BOOL
IsCurveStraight(NurbSurface<T> * n,T tolerance,int crvInd,BOOL dirflag)1620 IsCurveStraight( NurbSurface<T> * n,
1621 T tolerance,
1622 int crvInd,
1623 BOOL dirflag ) /* If TRUE_, test in U direction, else test in V */
1624 {
1625 Point_nD<T,3> p, vec, prod;
1626 Point_nD<T,3> cp, e0;
1627 int i, last;
1628 T linelen, dist;
1629
1630 /* Special case: lines are automatically straight. */
1631 if ((dirflag ? n->numU : n->numV) == 2)
1632 return( TRUE_ );
1633
1634 last = (dirflag ? n->numU : n->numV) - 1;
1635 n->render->screenProject( GETPT( 0 ), e0 );
1636
1637 /* Form an initial line to test the other points against (skiping degen lines) */
1638
1639 linelen = 0.0;
1640 for (i = last; (i > 0) && (linelen < EPSILON); i--)
1641 {
1642 n->render->screenProject( GETPT( i ), cp );
1643 vec = cp - e0 ;
1644 linelen = norm(vec) ;
1645 }
1646
1647 DIVPT( vec, linelen );
1648
1649 if (linelen > EPSILON) /* If no non-degenerate lines found, it's all degen */
1650 for (i = 1; i <= last; i++)
1651 {
1652 /* The cross product of the vector defining the
1653 * initial line with the vector of the current point
1654 * gives the distance to the line. */
1655 n->render->screenProject( GETPT( i ), cp );
1656 p = cp - e0 ;
1657
1658 prod = crossProduct(p,vec) ;
1659 dist = norm(prod) ;
1660
1661 if (dist > tolerance)
1662 return( FALSE_ );
1663 }
1664
1665 return( TRUE_ );
1666 }
1667
1668 /*
1669 * Check to see if a surface is flat. Tests are only performed on edges and
1670 * directions that aren't already straight. If an edge is flagged as straight
1671 * (from the parent surface) it is assumed it will stay that way.
1672 */
1673 template <class T> BOOL
TestFlat(NurbSurface<T> * n,T tolerance)1674 TestFlat( NurbSurface<T> * n, T tolerance )
1675 {
1676 int i;
1677 BOOL straight;
1678 Point_nD<T,3> cp00, cp0n, cpn0, cpnn, planeEqn;
1679 T dist,d ;
1680
1681 /* Check edge straightness */
1682
1683 if (! n->strU0)
1684 n->strU0 = IsCurveStraight( n, tolerance, 0, FALSE_ );
1685 if (! n->strUn)
1686 n->strUn = IsCurveStraight( n, tolerance, maxU(n), FALSE_ );
1687 if (! n->strV0)
1688 n->strV0 = IsCurveStraight( n, tolerance, 0, TRUE_ );
1689 if (! n->strVn)
1690 n->strVn = IsCurveStraight( n, tolerance, maxV(n), TRUE_ );
1691
1692 /* Test to make sure control points are straight in U and V */
1693
1694 straight = TRUE_;
1695 if ( (! n->flatU) && (n->strV0) && (n->strVn) )
1696 for (i = 1;
1697 (i < maxV(n)) && (straight = IsCurveStraight( n, tolerance, i, TRUE_ ));
1698 i++);
1699
1700 if (straight && n->strV0 && n->strVn)
1701 n->flatU = TRUE_;
1702
1703 straight = TRUE_;
1704 if ( (! n->flatV) && (n->strU0) && (n->strUn) )
1705 for (i = 1;
1706 (i < maxU(n)) && (straight = IsCurveStraight( n, tolerance, i, FALSE_ ));
1707 i++);
1708
1709 if (straight && n->strU0 && n->strUn)
1710 n->flatV = TRUE_;
1711
1712 if ( (! n->flatV) || (! n->flatU) )
1713 return( FALSE_ );
1714
1715 // The surface can pass the above tests but still be twisted.
1716
1717 n->render->screenProject( (n->points(0,0)), cp00 );
1718 n->render->screenProject( (n->points(0,maxU(n))), cp0n );
1719 n->render->screenProject( (n->points(maxV(n),0)), cpn0 );
1720 n->render->screenProject( (n->points(maxV(n),maxU(n))), cpnn );
1721
1722 cp0n -= cp00 ; // Make edges into vectors
1723 cpn0 -= cp00 ;
1724
1725
1726 // Compute the plane equation from two adjacent sides, and
1727 // measure the distance from the far point to the plane. If it's
1728 // larger than tolerance, the surface is twisted.
1729
1730 planeEqn = crossProduct(cpn0,cp0n) ;
1731 planeEqn = planeEqn.unitLength() ; // Normalize to keep adds in sync w/ mults
1732
1733 d = planeEqn * cp00 ;
1734 dist = fabs( ( planeEqn * cpnn ) - d );
1735
1736 if ( dist > tolerance ) // Surface is twisted
1737 return( FALSE_ );
1738 else
1739 return( TRUE_ );
1740 }
1741
1742 /*
1743 * Turn a sufficiently flat surface into triangles.
1744 */
1745 template <class T>
EmitTriangles(NurbSurface<T> * n)1746 void EmitTriangles( NurbSurface<T> * n )
1747 {
1748 Point_nD<T,3> vecnn, vec0n; // Diagonal vectors
1749 T len2nn, len20n; // Diagonal lengths squared
1750 T u0, un, v0, vn; // Texture coords;
1751
1752 //
1753 // Measure the distance aint the two diagonals to decide the best
1754 // way to cut the rectangle into triangles.
1755 //
1756
1757 vecnn = n->c00.point - n->cnn.point ;
1758 vec0n = n->c0n.point - n->cn0.point ;
1759
1760 len2nn = norm2(vecnn) ;
1761 len20n = norm2(vec0n) ;
1762
1763 if (maximum(len2nn, len20n) < n->epsilon)
1764 return; // Triangles are too small to render
1765
1766 //
1767 // Assign the texture coordinates
1768 //
1769 u0 = n->kvU[n->orderU-1];
1770 un = n->kvU[n->numU];
1771 v0 = n->kvV[n->orderV-1];
1772 vn = n->kvV[n->numV];
1773 n->c00.u = u0; n->c00.v = v0;
1774 n->c0n.u = un; n->c0n.v = v0;
1775 n->cn0.u = u0; n->cn0.v = vn;
1776 n->cnn.u = un; n->cnn.v = vn;
1777
1778 //
1779 // If any normals are sick, fix them now.
1780 //
1781 if ((n->c00.normLen == 0.0) || (n->cnn.normLen == 0.0) || (n->cn0.normLen == 0.0))
1782 FixNormals( &(n->c00), &(n->cnn), &(n->cn0) );
1783 if (n->c0n.normLen == 0.0)
1784 FixNormals( &(n->c00), &(n->c0n), &(n->cnn) );
1785
1786 if ( len2nn < len20n )
1787 {
1788 n->render->drawTriangle( n->c00, n->cnn, n->cn0 );
1789 n->render->drawTriangle( n->c00, n->c0n, n->cnn );
1790 }
1791 else
1792 {
1793 n->render->drawTriangle( n->c0n, n->cnn, n->cn0 );
1794 n->render->drawTriangle( n->c0n, n->cn0, n->c00 );
1795 }
1796 }
1797
1798 /*
1799 * The recursive subdivision algorithm. Test if the surface is flat.
1800 * If so, split it into triangles. Otherwise, split it into two halves,
1801 * and invoke the procedure on each half.
1802 */
1803 template <class T> void
DoSubdivision(NurbSurface<T> * n,T tolerance,BOOL dirflag,int level)1804 DoSubdivision( NurbSurface<T> * n, T tolerance, BOOL dirflag, int level )
1805 {
1806 NurbSurface<T> *left, *right;
1807
1808 left = new NurbSurface<T>;
1809 right = new NurbSurface<T>;
1810
1811 if (TestFlat( n, tolerance ))
1812 {
1813 EmitTriangles( n );
1814 }
1815 else
1816 {
1817 if ( ((! n->flatV) && (! n->flatU)) || ((n->flatV) && (n->flatU)) )
1818 dirflag = !dirflag; // If twisted or curved in both directions,
1819 else // then alternate subdivision direction
1820 {
1821 if (n->flatU) // Only split in directions that aren't flat
1822 dirflag = FALSE_;
1823 else
1824 dirflag = TRUE_;
1825 }
1826 SplitSurface( n, left, right, dirflag );
1827 DoSubdivision( left, tolerance, dirflag, level + 1 );
1828 DoSubdivision( right, tolerance, dirflag, level + 1 );
1829 FreeNurb( left );
1830 FreeNurb( right );
1831 }
1832 }
1833
1834 /*
1835 * Main entry point for subdivision */
1836 template <class T> void
DrawSubdivision(NurbSurface<T> * surf,T tolerance)1837 DrawSubdivision( NurbSurface<T> * surf, T tolerance )
1838 {
1839 surf->flatV = FALSE_;
1840 surf->flatU = FALSE_;
1841 surf->strU0 = FALSE_;
1842 surf->strUn = FALSE_;
1843 surf->strV0 = FALSE_;
1844 surf->strVn = FALSE_;
1845
1846 //
1847 // Initialize the projected corners of the surface
1848 // and the normals.
1849 //
1850 DIVW( &(surf->points(0,0)), &surf->c00.point );
1851 DIVW( &(surf->points(0,surf->numU-1)), &surf->c0n.point );
1852 DIVW( &(surf->points(surf->numV-1,0)), &surf->cn0.point );
1853 DIVW( &(surf->points(surf->numV-1,surf->numU-1)), &surf->cnn.point );
1854
1855 GetNormal( surf, 0, 0 );
1856 GetNormal( surf, 0, maxU(surf) );
1857 GetNormal( surf, maxV(surf), 0 );
1858 GetNormal( surf, maxV(surf), maxU(surf) );
1859
1860 RenderMesh<T> *render ;
1861
1862 render = surf->render ;
1863 render->drawHeader();
1864 DoSubdivision( surf, tolerance, TRUE_, 0 );
1865 // Note surf is deallocated by the subdivision process
1866 render->drawFooter();
1867 }
1868
1869 /*!
1870 \brief the copy operator
1871
1872 \param s the surface sample to copy
1873
1874 \author Philippe Lavoie
1875 \date 20 January 1999
1876 */
1877 template <class T>
operator =(const SurfSample<T> & s)1878 SurfSample<T>& SurfSample<T>::operator=(const SurfSample<T>& s) {
1879 point = s.point ;
1880 normal = s.normal ;
1881 normLen = s.normLen ;
1882 u = s.u ;
1883 v = s.v ;
1884 return *this ;
1885 }
1886
1887 /*
1888 template <class T>
1889 class NurbsCurveTess : public NurbsCurve<T,2> {
1890 public:
1891 NurbsCurveTess(const NurbsSurface<T,3>& rs) ;
1892
1893 void tesselate(T tol) ;
1894
1895 protected:
1896 int isStraight() ;
1897
1898 const NurbsSurface<T,3>& rsurf ;
1899 T tolerance ;
1900 }
1901
1902 template <class T>
1903 NurbsCurveTess::NurbsCurveTess(const NurbsSurface<T,3>& rs) : rsurf(rs) {
1904 tolerance = 0.01 ;
1905 }
1906
1907 template <class T>
1908 void NurbsCurveTess::tesselate(T tol){
1909 tolerance = tol ;
1910 if(knot.n() != 2*(deg_+1)){
1911
1912 }
1913 }
1914 */
1915
1916 }
1917