1 // libSNL - Simple Nurbs Library
2 // Copyright 2003 Scott A.E. Lanham, Australia.
3 // --------------------------------------------
4 // This program is free software; you can redistribute it and/or modify
5 // it under the terms of the GNU General Public License as published by
6 // the Free Software Foundation; either version 2 of the License, or
7 // (at your option) any later version.
8 //
9 // This program is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 // GNU General Public License for more details.
13
14 // *** General NURBS Surface ***
15
16 #include "snlSurface.h"
17 #include "snlUtil.h"
18
19 #include "snlNurbsCommon.h"
20
~snlSurface()21 snlSurface::~snlSurface()
22 {
23 if ( ctrlPtNet ) delete ctrlPtNet;
24
25 if ( knotVectU ) delete knotVectU;
26 if ( knotVectV ) delete knotVectV;
27
28 if ( trim_curves ) delete trim_curves;
29 }
30
snlSurface()31 snlSurface::snlSurface()
32 {
33 init();
34 }
35
init()36 void snlSurface::init()
37 {
38 //! Standard Initialisation.
39 // ------------------------
40
41 ctrlPtNet = 0;
42
43 knotVectU = 0;
44 knotVectV = 0;
45
46 trim_curves = new ptrList< snlCurve >;
47 }
48
snlSurface(const snlSurface & surfaceToCopy)49 snlSurface::snlSurface ( const snlSurface& surfaceToCopy )
50 {
51 //! Copy constructor.
52 // -----------------
53
54 init();
55
56 copyFrom ( surfaceToCopy );
57 }
58
snlSurface(int degreeU,int degreeV,unsigned sizeU,unsigned sizeV,snlPoint & origin,snlPoint & cornerMaxU,snlPoint & cornerMaxV)59 snlSurface::snlSurface ( int degreeU, int degreeV, unsigned sizeU, unsigned sizeV,
60 snlPoint& origin, snlPoint& cornerMaxU, snlPoint& cornerMaxV )
61 {
62 //! Construct a new NURBS surface.
63 // ------------------------------
64 //! @param degreeU Degree of surface in U direction.
65 //! @param degreeV Degree of surface in V direction.
66 //! @param sizeU Number of control points in the U dimension.
67 //! @param sizeV Number of control points in the V dimension.
68 //! @param origin Point at (u,v) = (0,0).
69 //! @param cornerMaxU Point at (u,v) = (MaxU,0).
70 //! @param cornerMaxV Point at (u,v) = (0,MaxV).
71
72 init();
73
74 ctrlPtNet = new snlCtrlPointNetSurface ( sizeU, sizeV, origin, cornerMaxU, cornerMaxV );
75
76 knotVectU = new snlKnotVector ( 0.0, 1.0, sizeU + degreeU + 1, degreeU );
77 knotVectV = new snlKnotVector ( 0.0, 1.0, sizeV + degreeV + 1, degreeV );
78
79 degU = degreeU;
80 degV = degreeV;
81 }
82
snlSurface(int degreeU,int degreeV,unsigned sizeU,unsigned sizeV,snlCtrlPoint * points,knot * knotsU,knot * knotsV)83 snlSurface::snlSurface ( int degreeU, int degreeV, unsigned sizeU, unsigned sizeV, snlCtrlPoint* points,
84 knot* knotsU, knot* knotsV )
85 {
86 //! Create surface from existing data.
87 // ----------------------------------
88 //! @param degreeU Degree in U direction.
89 //! @param degreeV Degree in V direction.
90 //! @param sizeU Size of U dimension.
91 //! @param sizeV Size of V dimension.
92 //! @param points Control points to use.
93 //! @param knotsU Array of knots for U direction.
94 //! @param knotsV Array of knots for V direction.
95 //!
96 //! @par Notes:
97 //! Assumes a clamped (open) knot vector.
98 //! @attention Does NOT COPY point and knot data. So don't delete them elsewhere.
99
100 init();
101
102 degU = degreeU;
103 degV = degreeV;
104
105 ctrlPtNet = new snlCtrlPointNetSurface ( points, sizeU, sizeV, false );
106
107 if ( knotsU )
108 knotVectU = new snlKnotVector ( knotsU, sizeU + degU + 1, degreeU );
109 else
110 knotVectU = new snlKnotVector ( 0.0, 1.0, sizeU + degU + 1, degreeU );
111
112 if ( knotsV )
113 knotVectV = new snlKnotVector ( knotsV, sizeV + degV + 1, degreeV );
114 else
115 knotVectV = new snlKnotVector ( 0.0, 1.0, sizeV + degV + 1, degreeV );
116 }
117
snlSurface(snlCurve & curve1,snlCurve & curve2,int direction)118 snlSurface::snlSurface ( snlCurve& curve1, snlCurve& curve2, int direction )
119 {
120 //! Generate ruled surface.
121 // -----------------------
122 //! @param curve1 First side of surface.
123 //! @param curve2 Second side of surface.
124 //! @param direction Parmetric direction defining curves lay in.
125
126 init();
127
128 // Curves may be modified so copy them.
129
130 snlCurve* curve_1 = new snlCurve ( curve1 );
131 snlCurve* curve_2 = new snlCurve ( curve2 );
132
133 // Make sure curves are compatible.
134
135 curve_1 -> makeCompatible ( curve_2 );
136
137 // Generate knot vectors.
138
139 if ( direction == SNL_U_DIR )
140 {
141 knotVectU = new snlKnotVector ( curve_1 -> knotVector() );
142 knotVectV = new snlKnotVector ( 0.0, 1.0, 4, 1 );
143 }
144 else
145 {
146 knotVectU = new snlKnotVector ( 0.0, 1.0, 4, 1 );
147 knotVectV = new snlKnotVector ( curve_1 -> knotVector() );
148 }
149
150 // Generate control points.
151
152 int size = curve_1 -> size();
153 int arraySize = size * 2;
154
155 snlCtrlPoint* ctrlPts = new snlCtrlPoint [ arraySize ];
156
157 const snlCtrlPoint* copyNet1 = curve_1 -> controlPointNet().getCtrlPts();
158 const snlCtrlPoint* copyNet2 = curve_2 -> controlPointNet().getCtrlPts();
159
160 int ctrlPtsIndex = 0;
161
162 if ( direction == SNL_U_DIR )
163 {
164 for ( int index = 0; index < size; index ++ )
165 {
166 ctrlPts [ ctrlPtsIndex ++ ] = copyNet1 [ index ];
167 ctrlPts [ ctrlPtsIndex ++ ] = copyNet2 [ index ];
168 }
169
170 ctrlPtNet = new snlCtrlPointNetSurface ( ctrlPts, size, 2, false );
171
172 degU = curve_1 -> degree();
173 degV = 1;
174 }
175 else
176 {
177 int ctrlPtsIndex2 = size;
178
179 for ( int index = 0; index < size; index ++ )
180 {
181 ctrlPts [ ctrlPtsIndex ++ ] = copyNet1 [ index ];
182 ctrlPts [ ctrlPtsIndex2 ++ ] = copyNet2 [ index ];
183 }
184
185 ctrlPtNet = new snlCtrlPointNetSurface ( ctrlPts, 2, size, false );
186
187 degU = 1;
188 degV = curve_1 -> degree();
189 }
190
191 // Clean up.
192
193 delete curve_1;
194 delete curve_2;
195 }
196
snlSurface(snlCurve & generator,snlPoint & axisStart,snlPoint & axisEnd,double angle)197 snlSurface::snlSurface ( snlCurve& generator, snlPoint& axisStart, snlPoint& axisEnd, double angle )
198 {
199 //! Construct surface of revolution.
200 // --------------------------------
201 //! @param generator Generating curve.
202 //! @param axisStart Starting point of axis generator is revolved about.
203 //! @param axisEnd Ending point of axis.
204 //! @param angle Angle in degrees to revolve generator about axis. Angle is in degrees so that
205 //! different precisions of PI do not affect surface closure.
206 //!
207 //! @par Notes:
208 //! Rotation is counter clockwise about axis vector ie Right hand rule. Curve defines V
209 //! direction.
210
211 init();
212
213 genSurfRevolution ( generator, axisStart, axisEnd, angle );
214 }
215
snlSurface(snlCurve ** curves,int numCurves,int dir)216 snlSurface::snlSurface ( snlCurve** curves, int numCurves, int dir )
217 {
218 //! Construct skinned surface.
219 // --------------------------
220 //! @param curves Curves to skin with.
221 //! @param numCurves Number of curves in curves array.
222 //! @param dir Parametric direction to skin over.
223 //!
224 //! @par Notes:
225 //! The skinned direction is degree 2. \n
226 //! It is assumed the curves are exactly the same size with the same knot vector.
227
228 init();
229
230 // Assemble points to be interpolated.
231
232 int numCurvePts = curves [ 0 ] -> size();
233 int numPts = numCurvePts * numCurves;
234
235 int sizeU = ( dir == SNL_U_DIR ? numCurves : numCurvePts );
236 int sizeV = ( dir == SNL_U_DIR ? numCurvePts : numCurves );
237
238 snlPoint* points = new snlPoint [ numPts ];
239
240 for ( int curveIndex = 0; curveIndex < numCurves; curveIndex ++ )
241 {
242 int index, step;
243
244 if ( dir == SNL_U_DIR )
245 {
246 index = curveIndex * numCurvePts;
247 step = 1;
248 }
249 else
250 {
251 index = curveIndex;
252 step = numCurves;
253 }
254
255 const snlCtrlPoint* curvePts = ( curves [ curveIndex ] -> controlPointNet() ).getCtrlPts();
256
257 for ( int ptIndex = 0; ptIndex < numCurvePts; ptIndex ++ )
258 {
259 points [ index ] = curvePts [ ptIndex ];
260
261 index += step;
262 }
263 }
264
265 // Generate parameters
266
267 knot* params = globalInterpGenParams ( SNL_GLOBAL_INTERP_CENTRIFUGAL, points, sizeU, sizeV, dir );
268
269 // Generate knot vectors.
270
271 if ( dir == SNL_U_DIR )
272 {
273 knotVectU = new snlKnotVector ( numCurves, 2, params );
274 knotVectV = new snlKnotVector ( curves [ 0 ] -> knotVector() );
275 }
276 else
277 {
278 knotVectU = new snlKnotVector ( curves [ 0 ] -> knotVector() );
279 knotVectV = new snlKnotVector ( numCurves, 2, params );
280 }
281
282 // Generate control points.
283
284 snlPoint* linePts = new snlPoint [ numCurves ];
285
286 snlCtrlPoint* ctrlPts = new snlCtrlPoint [ numPts ];
287
288 for ( int lineIndex = 0; lineIndex < numCurvePts; lineIndex ++ )
289 {
290 // Assemble points.
291
292 int index, step;
293
294 if ( dir == SNL_U_DIR )
295 {
296 index = lineIndex;
297 step = numCurvePts;
298 }
299 else
300 {
301 index = lineIndex * numCurves;
302 step = 1;
303 }
304
305 for ( int ptIndex = 0; ptIndex < numCurves; ptIndex ++ )
306 {
307 linePts [ ptIndex ] = points [ index ];
308 index += step;
309 }
310
311 // Interpolate between assembled points.
312
313 snlCtrlPoint* finalPts = snlCurve::genGlobalInterpPoints ( linePts, numCurves, params,
314 dir == SNL_U_DIR ? knotVectU : knotVectV );
315
316 // Copy points into surface control point array.
317
318 if ( dir == SNL_U_DIR )
319 {
320 index = lineIndex;
321 step = numCurvePts;
322 }
323 else
324 {
325 index = lineIndex * numCurves;
326 step = 1;
327 }
328
329 for ( int ptIndex = 0; ptIndex < numCurves; ptIndex ++ )
330 {
331 ctrlPts [ index ] = finalPts [ ptIndex ];
332 index += step;
333 }
334
335 delete[] finalPts;
336 }
337
338 ctrlPtNet = new snlCtrlPointNetSurface ( ctrlPts, sizeU, sizeV );
339
340 // Set other variables.
341
342 if ( dir == SNL_U_DIR )
343 {
344 degU = 2;
345 degV = curves [ 0 ] -> degree();
346 }
347 else
348 {
349 degU = curves [ 0 ] -> degree();
350 degV = 2;
351 }
352
353 // Clean up.
354
355 delete[] points;
356 delete[] params;
357 delete[] linePts;
358 }
359
snlSurface(int interpType,snlPoint * pointsInterp,int sizeU,int sizeV,int degreeU,int degreeV)360 snlSurface::snlSurface ( int interpType, snlPoint* pointsInterp, int sizeU, int sizeV, int degreeU, int degreeV )
361 {
362 //! Generate surface by interpolating point data.
363 // ---------------------------------------------
364 //! @param interpType Type of interpolation. Comes from enum SNL_INTERP_TYPES.
365 //! @param pointsInterp Points to interpolate. Point array.
366 //! @param sizeU Size of point array dimension in U direction.
367 //! @param sizeV Size of point array dimension in V direction.
368 //! @param degreeU Degree of U direction.
369 //! @param degreeV Degree of V direction.
370
371 init();
372
373 if ( interpType == SNL_GLOBAL_INTERP_CHORDLENGTH || interpType == SNL_GLOBAL_INTERP_CENTRIFUGAL )
374 genGlobalInterpSurf ( interpType, pointsInterp, sizeU, sizeV, degreeU, degreeV );
375 }
376
snlSurface(snlCurve * U1,snlCurve * U2,snlCurve * V1,snlCurve * V2)377 snlSurface::snlSurface ( snlCurve* U1, snlCurve* U2, snlCurve* V1, snlCurve* V2 )
378 {
379 //! Generate Coons Patch from given curves.
380 // ---------------------------------------
381 //! @param U1 First curve in U direction.
382 //! @param U2 Second curve in U direction.
383 //! @param V1 First curve in V direction.
384 //! @param V2 Second curve in V direction.
385
386 init();
387
388 genBilinearCoons( U1, U2, V1, V2 );
389 }
390
copyFrom(const snlSurface & surfaceToCopy)391 void snlSurface::copyFrom ( const snlSurface& surfaceToCopy )
392 {
393 //! Copy contents of another surface into this.
394 // -------------------------------------------
395
396 ctrlPtNet = new snlCtrlPointNetSurface ( *(surfaceToCopy.ctrlPtNet) );
397
398 knotVectU = new snlKnotVector ( *(surfaceToCopy.knotVectU) );
399 knotVectV = new snlKnotVector ( *(surfaceToCopy.knotVectV) );
400
401 degU = surfaceToCopy.degU;
402 degV = surfaceToCopy.degV;
403
404 snlCurve* trimCurve = ( surfaceToCopy.trim_curves ) -> first();
405
406 snlCurve* curveCopy;
407
408 while ( trimCurve )
409 {
410 curveCopy = new snlCurve ( *trimCurve );
411
412 addTrimCurve ( curveCopy );
413
414 trimCurve = ( surfaceToCopy.trim_curves ) -> next();
415 }
416 }
417
operator =(const snlSurface & surfaceToCopy)418 snlSurface& snlSurface::operator= ( const snlSurface& surfaceToCopy )
419 {
420 //! Assignment Operator.
421 // --------------------
422
423 if ( this != &surfaceToCopy )
424 {
425 // If surface already contains data then will have to delete it.
426
427 if ( ctrlPtNet ) delete ctrlPtNet;
428 if ( knotVectU ) delete knotVectU;
429 if ( knotVectV ) delete knotVectV;
430
431 trim_curves -> clear();
432
433 copyFrom ( surfaceToCopy );
434 }
435
436 return *this;
437 }
438
genBilinearCoons(snlCurve * curve_U1,snlCurve * curve_U2,snlCurve * curve_V1,snlCurve * curve_V2)439 void snlSurface::genBilinearCoons ( snlCurve* curve_U1, snlCurve* curve_U2, snlCurve* curve_V1, snlCurve* curve_V2 )
440 {
441 //! Generate Bilinear Coons Patch.
442 // ------------------------------
443 //! @param curve_U1 First curve in U direction.
444 //! @param curve_U2 Second curve in U direction.
445 //! @param curve_V1 First curve in V direction.
446 //! @param curve_V2 Second curve in V direction.
447 //!
448 //! @par Notes:
449 //! The curves must form a closed but non-sequential loop.
450 //! @verbatim
451 //! V1
452 //! -------->
453 //! | |
454 //! U1 | | U2
455 //! | |
456 //! V V
457 //! -------->
458 //! V2 @endverbatim
459
460 // Create Bilinear surface as base of this surface.
461
462 knotVectV = new snlKnotVector ( 0.0, 1.0, 4, 1 );
463 knotVectU = new snlKnotVector ( 0.0, 1.0, 4, 1 );
464
465 snlCtrlPoint* ctrlPts = new snlCtrlPoint [ 4 ];
466
467 ctrlPts [ 0 ] = ( curve_V1 -> controlPointNet() ).getPoint ( 0 );
468 ctrlPts [ 1 ] = ( curve_V1 -> controlPointNet() ).getPoint ( ( curve_V1 -> controlPointNet() ).size() - 1 );
469 ctrlPts [ 2 ] = ( curve_V2 -> controlPointNet() ).getPoint ( 0 );
470 ctrlPts [ 3 ] = ( curve_V2 -> controlPointNet() ).getPoint ( ( curve_V2 -> controlPointNet() ).size() - 1 );
471
472 degU = 1;
473 degV = 1;
474
475 ctrlPtNet = new snlCtrlPointNetSurface ( ctrlPts, 2, 2, false );
476
477 // Make sure curve pairs are compatible.
478
479 curve_U1 -> makeCompatible ( curve_U2 );
480 curve_V1 -> makeCompatible ( curve_V2 );
481
482 // Create U and V direction surfaces.
483
484 snlSurface* surf_U = new snlSurface ( *curve_U1, *curve_U2, (int) SNL_U_DIR );
485 snlSurface* surf_V = new snlSurface ( *curve_V1, *curve_V2, (int) SNL_V_DIR );
486
487 // Surfaces need to be compatible with each other.
488
489 surf_U -> makeCompatible ( surf_V, SNL_U_DIR );
490 surf_U -> makeCompatible ( surf_V, SNL_V_DIR );
491
492 makeCompatible ( surf_V, SNL_U_DIR );
493 makeCompatible ( surf_V, SNL_V_DIR );
494
495 makeCompatible ( surf_U, SNL_U_DIR );
496 makeCompatible ( surf_U, SNL_V_DIR );
497
498 // Generate new surface as addition of control points; surf_U + surf_V - bilinear ( this ).
499
500 *( surf_U -> ctrlPtNet ) += *( surf_V -> ctrlPtNet );
501 *( surf_U -> ctrlPtNet ) -= *ctrlPtNet;
502
503 snlCtrlPointNetSurface* newCtrlPtNet = surf_U -> ctrlPtNet; // Swap control point net into this.
504 surf_U -> ctrlPtNet = ctrlPtNet;
505 ctrlPtNet = newCtrlPtNet;
506
507 // Clean Up.
508
509 delete surf_U;
510 delete surf_V;
511
512 }
513
globalInterpGenParams(int type,snlPoint * points,int sizeU,int sizeV,int dir)514 knot* snlSurface::globalInterpGenParams ( int type, snlPoint* points, int sizeU, int sizeV,
515 int dir )
516 {
517 //! Generate parameters for global interpolation.
518 // ---------------------------------------------
519 //! @param type Type of global interpolation.
520 //! @param points Array points to interpolate. [ U ][ V ].
521 //! @param sizeU Array size in U direction.
522 //! @param sizeV Array size in V direction.
523 //! @param dir Point array direction to generate parameters for. See enum parametricDirections.
524 //!
525 //! @return 1D Array of knots that holds parametric positions in chosen direction. All
526 //! lines in dir have exactly the same parametric sequence. Hence this array is 1D.
527
528 int size, oSize;
529
530 if ( dir == SNL_U_DIR )
531 {
532 size = sizeU;
533 oSize = sizeV;
534 }
535 else
536 {
537 size = sizeV;
538 oSize = sizeU;
539 }
540
541 knot* chordLengths = new knot [ ( size - 1 ) * oSize ];
542
543 knot* totalChordLengths = new knot [ oSize ];
544
545 for ( int index = 0; index < oSize; index ++ )
546 totalChordLengths [ index ] = 0.0;
547
548 snlVector chord;
549
550 chord.homogeneous = true;
551
552 // Intermediate step. Calculate (square root of) chord length.
553
554 int chordIndex = 0;
555
556 for ( int cIndex = 1; cIndex < size; cIndex ++ )
557 {
558 int index;
559
560 if ( dir == SNL_U_DIR )
561 index = cIndex * sizeV;
562 else
563 index = cIndex;
564
565 for ( int oIndex = 0; oIndex < oSize; oIndex ++ )
566 {
567 if ( dir == SNL_U_DIR )
568 {
569 chord.calc ( points [ index - sizeV ], points [ index ] );
570 index ++;
571 }
572 else
573 {
574 chord.calc ( points [ index - 1 ], points [ index ] );
575 index += sizeV;
576 }
577
578 knot chordLength;
579
580 if ( type == SNL_GLOBAL_INTERP_CENTRIFUGAL )
581 chordLength = sqrt ( chord.length() );
582 else
583 chordLength = chord.length();
584
585 totalChordLengths [ oIndex ] += chordLength;
586
587 if ( chordLength == 0.0 )
588 {
589 cout << "Zero Chord Length @: " << cIndex << ", " << oIndex << "\n";
590 cout << "start pt: "; points [ index - sizeV - 1 ].print(); cout << "\n";
591 cout << "end pt: "; points [ index - 1 ].print(); cout << "\n";
592 }
593
594 chordLengths [ chordIndex ++ ] = chordLength;
595 }
596 }
597
598 // Calculate final parameter values.
599
600 knot* params = new knot [ size ];
601
602 params [ 0 ] = 0.0;
603 params [ size - 1 ] = 1.0;
604
605 chordIndex = 0;
606
607 for ( int cIndex = 1; cIndex < size - 1; cIndex ++ )
608 {
609 params [ cIndex ] = 0.0;
610
611 for ( int oIndex = 0; oIndex < oSize; oIndex ++ )
612 {
613 // Sum across all chords.
614 params [ cIndex ] += chordLengths [ chordIndex ++ ] / totalChordLengths [ oIndex ];
615 }
616
617 params [ cIndex ] /= (double) oSize;
618
619 params [ cIndex ] += params [ cIndex - 1 ];
620 }
621
622 delete[] totalChordLengths;
623
624 delete[] chordLengths;
625
626 return params;
627 }
628
genGlobalInterpSurf(int interpType,snlPoint * pointsInterp,int sizeU,int sizeV,int degreeU,int degreeV)629 void snlSurface::genGlobalInterpSurf ( int interpType, snlPoint* pointsInterp, int sizeU, int sizeV,
630 int degreeU, int degreeV )
631 {
632 //! Generate surface by globally interpolating point data.
633 // ------------------------------------------------------
634 //! @param interpType Type of interpolation. Comes from enum SNL_INTERP_TYPES.
635 //! @param pointsInterp Points to interpolate. Point array.
636 //! @param sizeU Size of point array dimension in U direction.
637 //! @param sizeV Size of point array dimension in V direction.
638 //! @param degreeU Degree of U direction.
639 //! @param degreeV Degree of V direction.
640
641 degU = degreeU;
642 degV = degreeV;
643
644 // Generate parameters that correspond to given points.
645
646 knot* paramsU = globalInterpGenParams ( interpType, pointsInterp, sizeU, sizeV, SNL_U_DIR );
647 knot* paramsV = globalInterpGenParams ( interpType, pointsInterp, sizeU, sizeV, SNL_V_DIR );
648
649 // Generate knot vectors.
650
651 knotVectU = new snlKnotVector ( sizeU, degreeU, paramsU );
652 knotVectV = new snlKnotVector ( sizeV, degreeV, paramsV );
653
654 // Generate intermediate control points using U direction.
655
656 int numPts = sizeU * sizeV;
657
658 snlPoint* intermPts = new snlPoint [ numPts ]; // Intermediate points.
659
660 snlPoint* curvePtsU = new snlPoint [ sizeU ];
661
662 for ( int indexV = 0; indexV < sizeV; indexV ++ )
663 {
664 // Find points for intermediate curve.
665
666 int curvePtsIndex = 0;
667
668 for ( int index = indexV; index < numPts; index += sizeV )
669 curvePtsU [ curvePtsIndex ++ ] = pointsInterp [ index ];
670
671 snlCtrlPoint* interpPoints = snlCurve::genGlobalInterpPoints ( curvePtsU, sizeU, paramsU, knotVectU );
672
673 // Place interpolted points into intermediate array.
674
675 curvePtsIndex = 0;
676
677 for ( int index = indexV; index < numPts; index += sizeV )
678 intermPts [ index ] = interpPoints [ curvePtsIndex ++ ];
679
680 delete[] interpPoints;
681 }
682
683 // Generate final control points.
684
685 snlCtrlPoint* finalPts = new snlCtrlPoint [ numPts ];
686
687 snlPoint* curvePtsV = new snlPoint [ sizeV ];
688
689 for ( int baseIndex = 0; baseIndex < numPts; baseIndex += sizeV )
690 {
691 // Each line in V direction is based at baseIndex.
692
693 int maxVLineIndex = baseIndex + sizeV;
694
695 // Find isoparametric curve points to evaluate.
696
697 int curvePtsIndex = 0;
698
699 for ( int vLineIndex = baseIndex; vLineIndex < maxVLineIndex; vLineIndex ++ )
700 curvePtsV [ curvePtsIndex ++ ] = intermPts [ vLineIndex ];
701
702 snlCtrlPoint* interpPoints = snlCurve::genGlobalInterpPoints ( curvePtsV, sizeV, paramsV, knotVectV );
703
704 // Place interpolated points into final control point array.
705
706 curvePtsIndex = 0;
707
708 for ( int vLineIndex = baseIndex; vLineIndex < maxVLineIndex; vLineIndex ++ )
709 finalPts [ vLineIndex ] = interpPoints [ curvePtsIndex ++ ];
710
711 delete[] interpPoints;
712 }
713
714 // Generate Control Point Net.
715
716 ctrlPtNet = new snlCtrlPointNetSurface ( finalPts, sizeU, sizeV, false ); // CtrlPtNet _owns_ finalPts array.
717
718 // Clean up.
719
720 delete[] intermPts;
721 delete[] curvePtsU;
722 delete[] curvePtsV;
723 delete[] paramsU;
724 delete[] paramsV;
725 }
726
degreeU() const727 int snlSurface::degreeU() const
728 {
729 //! Return degree of surface in U direction.
730 // ----------------------------------------
731
732 return degU;
733 }
734
degreeV() const735 int snlSurface::degreeV() const
736 {
737 //! Return degree of surface in V direction.
738 // ----------------------------------------
739
740 return degV;
741 }
742
sizeU() const743 unsigned snlSurface::sizeU() const
744 {
745 //! Return size of control point array in U direction.
746 // --------------------------------------------------
747
748 return ctrlPtNet -> getSizeU();
749 }
750
sizeV() const751 unsigned snlSurface::sizeV() const
752 {
753 //! Return size of control point array in V direction.
754 // --------------------------------------------------
755
756 return ctrlPtNet -> getSizeV();
757 }
758
759
minU()760 knot snlSurface::minU()
761 {
762 //! Return minimum U parameter.
763 // ---------------------------
764
765 return knotVectU -> min();
766 }
767
maxU()768 knot snlSurface::maxU()
769 {
770 //! Return maximum U parameter.
771 // ---------------------------
772
773 return knotVectU -> max();
774 }
775
minV()776 knot snlSurface::minV()
777 {
778 //! Return minimum V parameter.
779 // ---------------------------
780
781 return knotVectV -> min();
782 }
783
maxV()784 knot snlSurface::maxV()
785 {
786 //! Return maximum V parameter.
787 // ---------------------------
788
789 return knotVectV -> max();
790 }
791
controlPoints()792 const snlCtrlPoint* snlSurface::controlPoints()
793 {
794 //! Return pointer to array of surfaces' control points.
795 // ----------------------------------------------------
796
797 return ctrlPtNet -> getCtrlPts();
798 }
799
knotsU()800 const knot* snlSurface::knotsU()
801 {
802 //! Return pointer to array of knots in U direction.
803 // ------------------------------------------------
804
805 return knotVectU -> array();
806 }
807
knotsV()808 const knot* snlSurface::knotsV()
809 {
810 //! Return pointer to array of knots in V direction.
811 // ------------------------------------------------
812
813 return knotVectV -> array();
814 }
815
controlPointNet()816 snlCtrlPointNetSurface& snlSurface::controlPointNet()
817 {
818 //! Return reference to control point network object for surface.
819 // -------------------------------------------------------------
820
821 return *ctrlPtNet;
822 }
823
knotVectorU()824 const snlKnotVector& snlSurface::knotVectorU()
825 {
826 //! Return pointer to Knot vector in U direction.
827 // ---------------------------------------------
828
829 return *knotVectU;
830 }
831
knotVectorV()832 const snlKnotVector& snlSurface::knotVectorV()
833 {
834 //! Return pointer to Knot vector in U direction.
835 // ---------------------------------------------
836
837 return *knotVectV;
838 }
839
evalHmg(knot paramU,knot paramV,basis * basisU,basis * basisV) const840 snlPoint snlSurface::evalHmg ( knot paramU, knot paramV, basis* basisU, basis* basisV ) const
841 {
842 //! Evaluate Non Rational Homogeneous Surface Point.
843 // ------------------------------------------------
844 //! @param paramU Parameter in U direction to evaluate.
845 //! @param paramV Parameter in V direction to evaluate.
846 //! @param basisU Supplied basis function values. Must be 0 if not supplied.
847 //! @param basisV Supplied basis function values. Must be 0 if not supplied.
848 //!
849 //! @return Homogeneous point on surface.
850
851
852 snlPoint rPoint; // Return point.
853 snlPoint iPoint; // Intermediate point.
854
855 unsigned vStart; // The index where the V "line" starts in the control point array.
856
857 unsigned spanU = knotVectU -> findSpan ( paramU );
858 unsigned spanV = knotVectV -> findSpan ( paramV );
859
860 // Evaluate basis functions.
861
862 basis* bValsU;
863
864 if ( basisU )
865 bValsU = basisU;
866 else
867 bValsU = knotVectU -> evalBasis ( paramU );
868
869 basis* bValsV;
870
871 if ( basisV )
872 bValsV = basisV;
873 else
874 bValsV = knotVectV -> evalBasis ( paramV );
875
876 unsigned baseVIndex = spanV - (unsigned) degV; // Where in the V dimension the processing starts.
877
878 rPoint.w ( 0 );
879
880 // Get control point array.
881 const snlCtrlPoint* ctrlPts = ctrlPtNet -> getCtrlPts();
882
883 for ( int indexU = 0; indexU <= degU; indexU ++ )
884 {
885 iPoint.null(); // Set everything to zero.
886
887 vStart = ( spanU - (unsigned) degU + indexU ) * sizeV();
888
889 for ( int indexV = 0; indexV <= degV; indexV ++ )
890 {
891 snlPoint tmpPoint = ( ctrlPts [ vStart + baseVIndex + indexV ] ) * bValsV [ indexV ];
892 iPoint += tmpPoint;
893 }
894
895 snlPoint tmpPoint = iPoint * bValsU [ indexU ];
896 rPoint += tmpPoint;
897 }
898
899 if ( ! basisU ) delete[] bValsU;
900 if ( ! basisV ) delete[] bValsV;
901
902 return rPoint;
903 }
904
eval(knot paramU,knot paramV,basis * basisU,basis * basisV) const905 snlPoint snlSurface::eval ( knot paramU, knot paramV, basis* basisU, basis* basisV ) const
906 {
907 //! Evaluate rational non-homogeneous surface point.
908 // ------------------------------------------------
909 //! @param paramU Parameter in U direction to evaluate.
910 //! @param paramV Parameter in V direction to evaluate.
911 //! @param basisU Supplied basis function values. Must be 0 if not supplied.
912 //! @param basisV Supplied basis function values. Must be 0 if not supplied.
913 //!
914 //! @return Non-homogeneous point on surface.
915
916 knot minU = knotVectU -> min();
917 knot maxU = knotVectU -> max();
918 knot minV = knotVectV -> min();
919 knot maxV = knotVectV -> max();
920
921 // Clamp parameters.
922
923 if ( paramU > maxU )
924 {
925 cout.precision ( 16 );
926 //debug cout << "Surface Eval. Out of bounds U: " << paramU << " Min: " << minU << " Max: " << maxU << "\n";
927 paramU = maxU;
928 }
929
930 if ( paramU < minU )
931 {
932 cout.precision ( 16 );
933 //debug cout << "Surface Eval. Out of bounds U: " << paramU << " Min: " << minU << " Max: " << maxU << "\n";
934 paramU = minU;
935 }
936
937 if ( paramV > maxV )
938 {
939 cout.precision ( 16 );
940 //debug cout << "Surface Eval. Out of bounds V: " << paramV << " Min: " << minV << " Max: " << maxV << "\n";
941 paramV = maxV;
942 }
943
944 if ( paramV < minV )
945 {
946 cout.precision ( 16 );
947 //debug cout << "Surface Eval. Out of bounds V: " << paramV << " Min: " << minV << " Max: " << maxV << "\n";
948 paramV = minV;
949 }
950
951 snlPoint retPoint = evalHmg ( paramU, paramV, basisU, basisV );
952 retPoint.normalise();
953
954 return retPoint;
955 }
956
eval(knot paramU,knot paramV) const957 snlPoint snlSurface::eval ( knot paramU, knot paramV ) const
958 {
959 //! Evaluate rational non-homogeneous surface point.
960 // ------------------------------------------------
961 //! @param paramU Parameter in U direction to evaluate.
962 //! @param paramV Parameter in V direction to evaluate.
963 //!
964 //! @return Non-homogeneous point on surface.
965 //!
966 //! @par Notes:
967 //! Function duplication necessary to satisfy parent abstract class.
968
969 knot minU = knotVectU -> min();
970 knot maxU = knotVectU -> max();
971 knot minV = knotVectV -> min();
972 knot maxV = knotVectV -> max();
973
974 // Clamp parameters.
975
976 if ( paramU > maxU )
977 {
978 cout.precision ( 16 );
979 //debug cout << "Surface Eval. Out of bounds U: " << paramU << " Min: " << minU << " Max: " << maxU << "\n";
980 paramU = maxU;
981 }
982
983 if ( paramU < minU )
984 {
985 cout.precision ( 16 );
986 //debug cout << "Surface Eval. Out of bounds U: " << paramU << " Min: " << minU << " Max: " << maxU << "\n";
987 paramU = minU;
988 }
989
990 if ( paramV > maxV )
991 {
992 cout.precision ( 16 );
993 //debug cout << "Surface Eval. Out of bounds V: " << paramV << " Min: " << minV << " Max: " << maxV << "\n";
994 paramV = maxV;
995 }
996
997 if ( paramV < minV )
998 {
999 cout.precision ( 16 );
1000 //debug cout << "Surface Eval. Out of bounds V: " << paramV << " Min: " << minV << " Max: " << maxV << "\n";
1001 paramV = minV;
1002 }
1003
1004 snlPoint retPoint = evalHmg ( paramU, paramV );
1005 retPoint.normalise();
1006
1007 return retPoint;
1008 }
1009
evalDerivsHmg(knot paramU,knot paramV,unsigned derivU,unsigned derivV,basis * basisU,basis * basisV)1010 snlPoint* snlSurface::evalDerivsHmg ( knot paramU, knot paramV, unsigned derivU, unsigned derivV,
1011 basis* basisU, basis* basisV )
1012 {
1013 //! Evaluate Non Rational Homogeneous Surface Derivatives.
1014 // ------------------------------------------------------
1015 //! @param paramU Parameter in U direction to evaluate at.
1016 //! @param paramV Parameter in V direction to evaluate at.
1017 //! @param derivU Derivative order to evaluate in U direction.
1018 //! @param derivV Derivative order to evaluate in V direction.
1019 //! @param basisU Pre-computed basis functions values.
1020 //! @param basisV Pre-computed basis functions values.
1021 //!
1022 //! @return Array of snlPoint [ derivU + 1 ] [ derivV + 1 ]. Calling function
1023 //! must delete[] this array.
1024
1025 const snlPoint* cPnt;
1026 snlPoint* vPnts = new snlPoint [ derivV + 1 ];
1027
1028 // Find spans
1029 unsigned spanU = knotVectU -> findSpan ( paramU );
1030 unsigned spanV = knotVectV -> findSpan ( paramV );
1031
1032 // Evaluate basis functions.
1033
1034 basis* bValsU;
1035 basis* bValsV;
1036
1037 if ( basisU )
1038 bValsU = basisU;
1039 else
1040 bValsU = knotVectU -> evalBasisDeriv ( paramU, derivU );
1041
1042 if ( basisV )
1043 bValsV = basisV;
1044 else
1045 bValsV = knotVectV -> evalBasisDeriv ( paramV, derivV );
1046
1047 unsigned baseVIndex = spanV - degV; // Where in the V dimension the processing starts.
1048
1049 unsigned evalPtsSize = ( derivU + 1 ) * ( derivV + 1 );
1050
1051 snlPoint* evalPts = new snlPoint [ evalPtsSize ];
1052
1053 // Set eval points net to 0.
1054 for ( unsigned index = 0; index < evalPtsSize; index ++ )
1055 evalPts [ index ].null(); // Set x, y, z and w to zero.
1056
1057 // Get control point array.
1058 const snlCtrlPoint* ctrlPts = ctrlPtNet -> getCtrlPts();
1059
1060 // Just loop through control points that match non-zero basis funtions.
1061
1062 for ( unsigned indexU = 0; indexU <= (unsigned) degU; indexU ++ )
1063 {
1064 unsigned vStart = ( spanU - degU + indexU ) * sizeV();
1065
1066 // Zero vPnts array.
1067 for ( unsigned index = 0; index <= derivV; index ++ )
1068 vPnts [ index ].null();
1069
1070 // Utilise V direction basis values.
1071 for ( unsigned indexV = 0; indexV <= (unsigned) degV; indexV ++ )
1072 {
1073 // Find control point only once.
1074 cPnt = ctrlPts + vStart + baseVIndex + indexV;
1075
1076 // Calculate all V deriv values based on this point.
1077 for ( unsigned dIndexV = 0; dIndexV <= derivV; dIndexV ++ )
1078 {
1079 // Calc index into basis derivs for V direction.
1080 unsigned vDerivStart = ( dIndexV * ( degV + 1 ) ) + indexV;
1081
1082 snlPoint tmpPoint = ( *cPnt ) * bValsV [ vDerivStart ];
1083 vPnts [ dIndexV ] += tmpPoint;
1084 }
1085 }
1086
1087 // Multiply U deriv basis values and add to eval array.
1088
1089 for ( unsigned dIndexU = 0; dIndexU <= derivU; dIndexU ++ )
1090 {
1091 // Get basis value for current u deriv level.
1092 basis cBValU = bValsU [ dIndexU * ( degU + 1 ) + indexU ];
1093
1094 for ( unsigned dIndexV = 0; dIndexV <= derivV; dIndexV ++ )
1095 {
1096 unsigned evalIndex = dIndexU * ( derivV + 1 ) + dIndexV;
1097
1098 snlPoint tmpPoint = vPnts [ dIndexV ] * cBValU;
1099 evalPts [ evalIndex ] += tmpPoint;
1100 }
1101 }
1102 }
1103
1104 delete[] vPnts;
1105
1106 if ( ! basisU ) delete[] bValsU;
1107 if ( ! basisV ) delete[] bValsV;
1108
1109 return evalPts;
1110 }
1111
evalDerivs(knot paramU,knot paramV,unsigned derivU,unsigned derivV)1112 snlPoint* snlSurface::evalDerivs ( knot paramU, knot paramV, unsigned derivU, unsigned derivV )
1113 {
1114 //! Evaluate Rational Non-Homogeneous Surface Derivatives
1115 // -----------------------------------------------------
1116 //! @param paramU Parameter in U direction to evaluate at.
1117 //! @param paramV Parameter in V direction to evaluate at.
1118 //! @param derivU Derivative order to evaluate in U direction.
1119 //! @param derivV Derivative order to evaluate in V direction.
1120 //!
1121 //! @return Array of snlPoint [ derivU + 1 ] [ derivV + 1 ]. Calling function must
1122 //! delete[] array.
1123 //!
1124 //! @par Notes:
1125 //! Follows theory, "The NURBS Book 2nd ed", page 136 equation 4.20.
1126
1127 unsigned kIndex, lIndex; // Partial derivs in t and u directions respectively.
1128 unsigned iIndex, jIndex; // Current partial deriv levels for t and u directions respectively.
1129 unsigned cDIndex; // Current derivPts index.
1130 unsigned kBaseIndex; // Base k index into evalPts array.
1131 unsigned kBaseIndex2;
1132
1133 snlPoint iPoint; // Intermediate point. Holds x, y and z of homogeneous point.
1134 snlPoint iPoint2; // Another intermediate point.
1135
1136 // Get homogeneous derivatives.
1137 snlPoint* derivPts = evalDerivsHmg ( paramU, paramV, derivU, derivV );
1138
1139 // Array for returning points in.
1140 snlPoint* evalPts = new snlPoint [ ( derivU + 1 ) * ( derivV + 1 ) ];
1141
1142 unsigned dSizeV = derivV + 1;
1143
1144 for ( kIndex = 0; kIndex <= derivU; kIndex ++ )
1145 {
1146 kBaseIndex = kIndex * ( derivV + 1 );
1147
1148 for ( lIndex = 0; lIndex <= derivV; lIndex ++ )
1149 {
1150 cDIndex = kIndex * dSizeV + lIndex;
1151
1152 iPoint = derivPts [ cDIndex ];
1153
1154 // 0th k order derivatives.
1155 for ( jIndex = 1; jIndex <= lIndex; jIndex ++ )
1156 {
1157 iPoint.x ( iPoint.x() - binCoefs::binCoefArray [ lIndex ] [ jIndex ] * derivPts [ jIndex ].w() *
1158 evalPts [ kBaseIndex + lIndex - jIndex ].x() );
1159
1160 iPoint.y ( iPoint.y() - binCoefs::binCoefArray [ lIndex ] [ jIndex ] * derivPts [ jIndex ].w() *
1161 evalPts [ kBaseIndex + lIndex - jIndex ].y() );
1162
1163 iPoint.z ( iPoint.z() - binCoefs::binCoefArray [ lIndex ] [ jIndex ] * derivPts [ jIndex ].w() *
1164 evalPts [ kBaseIndex + lIndex - jIndex ].z() );
1165 }
1166
1167 // 0th j order and ( k, j ) order derivatives.
1168 for ( iIndex = 1; iIndex <= kIndex; iIndex ++ )
1169 {
1170 kBaseIndex2 = ( kIndex - iIndex ) * ( derivV + 1 );
1171
1172 cDIndex = iIndex * dSizeV;
1173
1174 iPoint.x ( iPoint.x() - binCoefs::binCoefArray [ kIndex ] [ iIndex ] * derivPts [ cDIndex ].w() *
1175 evalPts [ kBaseIndex2 + lIndex ].x() );
1176
1177 iPoint.y ( iPoint.y() - binCoefs::binCoefArray [ kIndex ] [ iIndex ] * derivPts [ cDIndex ].w() *
1178 evalPts [ kBaseIndex2 + lIndex ].y() );
1179
1180 iPoint.z ( iPoint.z() - binCoefs::binCoefArray [ kIndex ] [ iIndex ] * derivPts [ cDIndex ].w() *
1181 evalPts [ kBaseIndex2 + lIndex ].z() );
1182
1183 iPoint2.null();
1184
1185 for ( jIndex = 1; jIndex <= lIndex; jIndex ++ )
1186 {
1187 iPoint2.x ( iPoint2.x() + binCoefs::binCoefArray [ lIndex ] [ jIndex ] *
1188 derivPts [ cDIndex + jIndex ].w() *
1189 evalPts [ kBaseIndex2 + lIndex - jIndex ].x() );
1190
1191 iPoint2.y ( iPoint2.y() + binCoefs::binCoefArray [ lIndex ] [ jIndex ] *
1192 derivPts [ cDIndex + jIndex ].w() *
1193 evalPts [ kBaseIndex2 + lIndex - jIndex ].y() );
1194
1195 iPoint2.z ( iPoint2.z() + binCoefs::binCoefArray [ lIndex ] [ jIndex ] *
1196 derivPts [ cDIndex + jIndex ].w() *
1197 evalPts [ kBaseIndex2 + lIndex - jIndex ].z() );
1198 }
1199
1200 snlPoint tmpPoint = iPoint2 * binCoefs::binCoefArray [ kIndex ] [ iIndex ];
1201 iPoint -= tmpPoint;
1202 }
1203
1204 evalPts [ kBaseIndex + lIndex ] = iPoint / ( derivPts [ 0 ].w() );
1205 evalPts [ kBaseIndex + lIndex ].w ( 0 ); // No delta in w component at any time.
1206 }
1207 }
1208
1209 delete[] derivPts;
1210
1211 evalPts [ 0 ].w ( 1.0 ); // First entry in array is not a derivative.
1212
1213 return evalPts;
1214 }
1215
velocities(knot paramU,knot paramV,snlPoint & evalPoint,snlVector & velocityU,snlVector & velocityV,basis * basisU,basis * basisV)1216 void snlSurface::velocities ( knot paramU, knot paramV, snlPoint& evalPoint, snlVector& velocityU, snlVector& velocityV,
1217 basis* basisU, basis* basisV )
1218 {
1219 //! Compute velocities in U and V directions.
1220 // -----------------------------------------
1221 //! @param paramU U parameter to evaluate at.
1222 //! @param paramV V parameter to evaluate at.
1223 //! @param evalPoint Variable to return surface point at which velocity was evaluated.
1224 //! @param velocityU Varibale to return U directon velocity in.
1225 //! @param velocityV Varibale to return V directon velocity in.
1226 //! @param basisU Pre-computed basis function values.
1227 //! @param basisV Pre-computed basis funciton values.
1228 //!
1229 //! @par Notes:
1230 //! This function is designed for speed not correctness.
1231
1232 knot minU = knotVectU -> min();
1233 knot maxU = knotVectU -> max();
1234 knot minV = knotVectV -> min();
1235 knot maxV = knotVectV -> max();
1236
1237 // Clamp parameters.
1238
1239 if ( paramU > maxU )
1240 {
1241 cout.precision ( 16 );
1242 cout << "Surface Velocity Eval. Out of bounds U: " << paramU << " Min: "
1243 << minU << " Max: " << maxU << "\n";
1244 paramU = maxU;
1245 }
1246
1247 if ( paramU < minU )
1248 {
1249 cout.precision ( 16 );
1250 cout << "Surface Velocity Eval. Out of bounds U: " << paramU << " Min: " << minU
1251 << " Max: " << maxU << "\n";
1252 paramU = minU;
1253 }
1254
1255 if ( paramV > maxV )
1256 {
1257 cout.precision ( 16 );
1258 cout << "Surface Velocity Eval. Out of bounds U: " << paramV << " Min: " << minU
1259 << " Max: " << maxU << "\n";
1260 paramV = maxV;
1261 }
1262
1263 if ( paramV < minV )
1264 {
1265 cout.precision ( 16 );
1266 cout << "Surface Velocity Eval. Out of bounds U: " << paramV << " Min: " << minU
1267 << " Max: " << maxU << "\n";
1268 paramV = minV;
1269 }
1270
1271 // Get homogeneous velocities.
1272
1273 snlPoint* derivPts = evalDerivsHmg ( paramU, paramV, 1, 1, basisU, basisV );
1274
1275 basis wVal = derivPts [ 0 ].elements [ 3 ];
1276
1277 derivPts [ 0 ].normalise();
1278
1279 evalPoint = derivPts [ 0 ];
1280
1281 basis uWVal = derivPts [ 2 ].elements [ 3 ];
1282
1283 velocityU.elements [ 0 ] = ( derivPts [ 2 ].elements [ 0 ] - uWVal
1284 * derivPts [ 0 ].elements [ 0 ] ) / wVal;
1285 velocityU.elements [ 1 ] = ( derivPts [ 2 ].elements [ 1 ] - uWVal
1286 * derivPts [ 0 ].elements [ 1 ] ) / wVal;
1287 velocityU.elements [ 2 ] = ( derivPts [ 2 ].elements [ 2 ] - uWVal
1288 * derivPts [ 0 ].elements [ 2 ] ) / wVal;
1289
1290 basis vWVal = derivPts [ 1 ].elements [ 3 ];
1291
1292 velocityV.elements [ 0 ] = ( derivPts [ 1 ].elements [ 0 ] - vWVal
1293 * derivPts [ 0 ].elements [ 0 ] ) / wVal;
1294 velocityV.elements [ 1 ] = ( derivPts [ 1 ].elements [ 1 ] - vWVal
1295 * derivPts [ 0 ].elements [ 1 ] ) / wVal;
1296 velocityV.elements [ 2 ] = ( derivPts [ 1 ].elements [ 2 ] - vWVal
1297 * derivPts [ 0 ].elements [ 2 ] ) / wVal;
1298
1299 // Clean up.
1300
1301 delete[] derivPts;
1302 }
1303
insertKnot(knot iParam,int dir,bool reallocate)1304 void snlSurface::insertKnot ( knot iParam, int dir, bool reallocate )
1305 {
1306 //! For a Surface: Insert a Knot into Knot Vector and Calculate new Control Points
1307 // ------------------------------------------------------------------------------
1308 //! @param iParam Parameter value to insert.
1309 //! @param dir Direction to evaluate in. 0 = u, 1 = v.
1310 //! @param reallocate Reallocate memory for control points.
1311 //!
1312 //! @par Notes:
1313 //! ctrlPts MUST have an additional point space allocated at the end of
1314 //! each line in the array for the chosen direction.
1315
1316 unsigned count, index, lineIndex, offset;
1317 unsigned cDeg, oDeg; // Degree to be processed.
1318 snlKnotVector* cKnts; // Current knots
1319 snlKnotVector* oKnts; // Other knots.
1320 snlCtrlPoint pt1, pt2;
1321
1322 if ( dir == SNL_V_DIR )
1323 {
1324 cDeg = degV;
1325 oDeg = degU;
1326 cKnts = knotVectV;
1327 oKnts = knotVectU;
1328
1329 ctrlPtNet -> growV ( 1, reallocate );
1330 }
1331 else
1332 {
1333 cDeg = degU;
1334 oDeg = degV;
1335 cKnts = knotVectU;
1336 oKnts = knotVectV;
1337
1338 ctrlPtNet -> growU ( 1, reallocate );
1339 }
1340
1341 // Span new knot belongs to.
1342 unsigned span = cKnts -> findSpan ( iParam );
1343
1344 // Pre calculate alphas.
1345 double* alpha = new double [ cDeg ];
1346
1347 for ( count = 0; count < cDeg; count ++ )
1348 {
1349 index = span - cDeg + 1 + count;
1350 alpha [ count ] = ( iParam - ( cKnts -> val ( index ) ) )
1351 / ( cKnts -> val ( index + cDeg ) - cKnts -> val ( index ) );
1352 }
1353
1354 // Build temp array to store new array of control points in.
1355 snlCtrlPoint* tmpPts = new snlCtrlPoint [ cDeg ];
1356
1357 // Get pointer to control points.
1358 snlCtrlPoint* ctrlPts = ctrlPtNet -> getCtrlPtsPtr();
1359
1360 // Do for each "line" in direction of insertion.
1361 for ( lineIndex = 0; lineIndex < ( oKnts -> size() - oDeg - 1 ); lineIndex ++ )
1362 {
1363 // Calculate new control points.
1364 for ( count = 0; count < cDeg; count ++ )
1365 {
1366 index = span - cDeg + 1 + count;
1367
1368 // Get first and second ctrl points to process with
1369 if ( dir == SNL_V_DIR )
1370 {
1371 // V direction.
1372 pt1 = ctrlPts [ ( lineIndex * sizeV() ) + index ];
1373 pt2 = ctrlPts [ ( lineIndex * sizeV() ) + index - 1 ];
1374 }
1375 else
1376 {
1377 // U direction.
1378 pt1 = ctrlPts [ ( index * sizeV() ) + lineIndex ];
1379 pt2 = ctrlPts [ ( ( index - 1 ) * sizeV() ) + lineIndex ];
1380 }
1381
1382 tmpPts [ count ].x ( ( alpha[count] * pt1.x() ) + ( ( 1.0 - alpha[count] ) * pt2.x() ) );
1383 tmpPts [ count ].y ( ( alpha[count] * pt1.y() ) + ( ( 1.0 - alpha[count] ) * pt2.y() ) );
1384 tmpPts [ count ].z ( ( alpha[count] * pt1.z() ) + ( ( 1.0 - alpha[count] ) * pt2.z() ) );
1385 tmpPts [ count ].w ( ( alpha[count] * pt1.w() ) + ( ( 1.0 - alpha[count] ) * pt2.w() ) );
1386 }
1387
1388 // Place new points into array.
1389
1390 if ( dir == SNL_V_DIR )
1391 {
1392 // V direction insert
1393 offset = lineIndex * sizeV();
1394
1395 // Copy non-altered control points forward one position at the end of the array.
1396 for ( count = sizeV() - 1; count > span; count -- )
1397 ctrlPts [ offset + count ] = ctrlPts [ offset + count - 1 ];
1398
1399 // Copy new control points into array.
1400 for ( count = 0; count < cDeg; count ++ )
1401 {
1402 index = span - cDeg + 1 + count + offset;
1403 ctrlPts [ index ] = tmpPts [ count ];
1404 }
1405 }
1406 else
1407 {
1408 // U direction insert.
1409
1410 // Copy non-altered control points forward one position at the end of the array.
1411 for ( count = sizeU() - 1; count > span; count -- )
1412 {
1413 index = count * sizeV() + lineIndex;
1414 ctrlPts [ index ] = ctrlPts [ index - sizeV() ];
1415 }
1416
1417 // Copy new control points into array.
1418 for ( count = 0; count < cDeg; count ++ )
1419 {
1420 index = ( span - cDeg + 1 + count ) * sizeV() + lineIndex;
1421 ctrlPts [ index ] = tmpPts [ count ];
1422 }
1423 }
1424 }
1425
1426 // Insert new knot into knot vector
1427 cKnts -> insertKnot ( iParam );
1428
1429 delete[] tmpPts;
1430 delete[] alpha;
1431 }
1432
insertKnot(knot iParam,int dir,int numToInsert,bool reallocate)1433 void snlSurface::insertKnot ( knot iParam, int dir, int numToInsert, bool reallocate )
1434 {
1435 //! For a Surface: Insert a Knot into Knot Vector and Calculate new Control Points
1436 // ------------------------------------------------------------------------------
1437 //! @param iParam Parameter value to insert.
1438 //! @param dir Direction to evaluate in. 0 = u, 1 = v.
1439 //! @param numToInsert Number of knots to insert into location.
1440 //! @param reallocate Reallocate memory for control points.
1441 //!
1442 //! @par Notes:
1443 //! ctrlPts MUST have an additional point space allocated at the end of
1444 //! each line in the array for the chosen direction.
1445
1446 if ( numToInsert <= 0 )
1447 {
1448 cout << "Bad use of function: snlSurface::insertKnot. Can't insert " << numToInsert << " knots.\n";
1449 return;
1450 }
1451
1452 int count, index, lineIndex, offset;
1453 unsigned cDeg, oDeg; // Degree to be processed.
1454 snlKnotVector* cKnts; // Current knots
1455 snlKnotVector* oKnts; // Other knots.
1456 snlCtrlPoint pt1, pt2;
1457
1458 if ( dir == SNL_V_DIR )
1459 {
1460 cDeg = degV;
1461 oDeg = degU;
1462 cKnts = knotVectV;
1463 oKnts = knotVectU;
1464
1465 ctrlPtNet -> growV ( numToInsert, reallocate );
1466 }
1467 else
1468 {
1469 cDeg = degU;
1470 oDeg = degV;
1471 cKnts = knotVectU;
1472 oKnts = knotVectV;
1473
1474 ctrlPtNet -> growU ( numToInsert, reallocate );
1475 }
1476
1477 // Span new knot belongs to.
1478 unsigned span = cKnts -> findSpan ( iParam );
1479
1480 int multi = cKnts -> findMultiplicity ( iParam ); // Multiplicity of knot vector at param.
1481
1482 // Pre calculate alphas.
1483
1484 int numAlphas = cDeg - multi;
1485
1486 double* alpha = new double [ numAlphas * numToInsert ];
1487
1488 for ( int insertCount = 1; insertCount <= numToInsert; insertCount ++ )
1489 {
1490 for ( count = 0; count < numAlphas - insertCount + 1; count ++ )
1491 {
1492 index = span - cDeg + insertCount + count;
1493
1494 int alphaIndex = count + ( ( insertCount - 1 ) * numAlphas );
1495
1496 alpha [ alphaIndex ] = ( iParam - ( cKnts -> val ( index ) ) )
1497 / ( cKnts -> val ( index + cDeg - insertCount + 1 )
1498 - cKnts -> val ( index ) );
1499 }
1500 }
1501
1502 // Build temp array to store new array of control points in.
1503
1504 int numNewPts = cDeg - multi + numToInsert - 1;
1505
1506 snlCtrlPoint* tmpPts = new snlCtrlPoint [ numNewPts ];
1507
1508 // Get pointer to control points.
1509 snlCtrlPoint* ctrlPts = ctrlPtNet -> getCtrlPtsPtr();
1510
1511 // Do for each "line" in direction of insertion.
1512 for ( lineIndex = 0; lineIndex < (int)( oKnts -> size() - oDeg - 1 ); lineIndex ++ )
1513 {
1514 // Calculate new control points for first insertion.
1515 for ( count = 0; count < numAlphas; count ++ )
1516 {
1517 index = span - cDeg + 1 + count;
1518
1519 // Get first and second ctrl points to process with
1520 if ( dir == SNL_V_DIR )
1521 {
1522 // V direction.
1523 pt1 = ctrlPts [ ( lineIndex * sizeV() ) + index ];
1524 pt2 = ctrlPts [ ( lineIndex * sizeV() ) + index - 1 ];
1525 }
1526 else
1527 {
1528 // U direction.
1529 pt1 = ctrlPts [ ( index * sizeV() ) + lineIndex ];
1530 pt2 = ctrlPts [ ( ( index - 1 ) * sizeV() ) + lineIndex ];
1531 }
1532
1533 tmpPts [ count ].x ( ( alpha[count] * pt1.x() ) + ( ( 1.0 - alpha[count] ) * pt2.x() ) );
1534 tmpPts [ count ].y ( ( alpha[count] * pt1.y() ) + ( ( 1.0 - alpha[count] ) * pt2.y() ) );
1535 tmpPts [ count ].z ( ( alpha[count] * pt1.z() ) + ( ( 1.0 - alpha[count] ) * pt2.z() ) );
1536 tmpPts [ count ].w ( ( alpha[count] * pt1.w() ) + ( ( 1.0 - alpha[count] ) * pt2.w() ) );
1537 }
1538
1539 // Calculate subsequent points for each insertion.
1540
1541 for ( int insertCount = 0; insertCount < numToInsert - 1; insertCount ++ )
1542 {
1543 // Copy some last points calculated forward one position.
1544
1545 for ( int copyCount = 0; copyCount < insertCount + 1; copyCount ++ )
1546 tmpPts [ numAlphas + insertCount - copyCount ] = tmpPts [ numAlphas + insertCount - copyCount - 1 ];
1547
1548 // Calculate new points
1549
1550 int alphaOffset = ( insertCount + 1 ) * numAlphas;
1551
1552 index = numAlphas - 1;
1553
1554 for ( count = numAlphas - insertCount - 2; count >= 0 ; count -- )
1555 {
1556 tmpPts [ index ].x ( ( alpha [ alphaOffset + count ] * tmpPts [ index ].x() )
1557 + ( ( 1.0 - alpha [ alphaOffset + count ] ) * tmpPts [ index - 1 ].x() ) );
1558
1559 tmpPts [ index ].y ( ( alpha [ alphaOffset + count ] * tmpPts [ index ].y() )
1560 + ( ( 1.0 - alpha [ alphaOffset + count ] ) * tmpPts [ index - 1 ].y() ) );
1561
1562 tmpPts [ index ].z ( ( alpha [ alphaOffset + count ] * tmpPts [ index ].z() )
1563 + ( ( 1.0 - alpha [ alphaOffset + count ] ) * tmpPts [ index - 1 ].z() ) );
1564
1565 tmpPts [ index ].w ( ( alpha [ alphaOffset + count ] * tmpPts [ index ].w() )
1566 + ( ( 1.0 - alpha [ alphaOffset + count ] ) * tmpPts [ index - 1 ].w() ) );
1567
1568 -- index;
1569 }
1570 }
1571
1572 // Place new points into array.
1573
1574 if ( dir == SNL_V_DIR )
1575 {
1576 // V direction insert
1577 offset = lineIndex * sizeV();
1578
1579 // Copy non-altered control points forward.
1580 for ( count = sizeV() - 1; count > (int) ( span - multi ); count -- )
1581 ctrlPts [ offset + count ] = ctrlPts [ offset + count - numToInsert ];
1582
1583 index = span - cDeg + 1 + offset;
1584
1585 // Copy new control points into array.
1586 for ( count = 0; count < numNewPts; count ++ )
1587 {
1588 ctrlPts [ index ] = tmpPts [ count ];
1589 ++ index;
1590 }
1591 }
1592 else
1593 {
1594 // U direction insert.
1595
1596 int backIndexOffset = sizeV() * numToInsert;
1597
1598 index = ( sizeU() - 1 ) * sizeV() + lineIndex;
1599
1600 // Copy non-altered control points forward one position at the end of the array.
1601 for ( count = sizeU() - 1; count > (int) ( span - multi ); count -- )
1602 {
1603 ctrlPts [ index ] = ctrlPts [ index - backIndexOffset ];
1604 index -= sizeV();
1605 }
1606
1607 // Copy new control points into array.
1608
1609 index = ( span - cDeg + 1 ) * sizeV() + lineIndex;
1610
1611 for ( count = 0; count < numNewPts; count ++ )
1612 {
1613 ctrlPts [ index ] = tmpPts [ count ];
1614 index += sizeV();
1615 }
1616 }
1617 }
1618
1619 // Insert new knot into knot vector
1620 cKnts -> insertKnot ( iParam, numToInsert );
1621
1622 delete[] tmpPts;
1623 delete[] alpha;
1624 }
1625
removeKnots(int numKnots,unsigned removalIndex,int direction,double tolerance,bool reallocate)1626 double snlSurface::removeKnots ( int numKnots, unsigned removalIndex, int direction,
1627 double tolerance, bool reallocate )
1628 {
1629 //! For a Surface: Remove multiple Knots from Knot Vector and Calculate new Control Points
1630 // --------------------------------------------------------------------------------------
1631 //! @param numKnots Number of knots to remove.
1632 //! @param removalIndex Knot at index to remove.
1633 //! @param direction Direction to evaluate in. 0 = u, 1 = v.
1634 //! @param tolerance Maximum allowable error. A value of 0.0 means ignore.
1635 //! @param reallocate Reallocate memory for control points. If false then de-allocation of
1636 //! memory is left up to caller.
1637 //!
1638 //! @return Maximum error encountered.
1639
1640 if ( numKnots < 1 ) return 0.0;
1641
1642 double maxTol = 0.0;
1643
1644 snlKnotVector* knotVect;
1645
1646 if ( direction == SNL_U_DIR )
1647 {
1648 knotVect = knotVectU;
1649 }
1650 else
1651 {
1652 knotVect = knotVectV;
1653 }
1654
1655 double param = knotVect -> val ( removalIndex );
1656
1657 int multi = knotVect -> findMultiplicity ( removalIndex );
1658
1659 int numToRemove = numKnots > multi ? multi : numKnots;
1660
1661 for ( int count = 0; count < numToRemove; count ++ )
1662 {
1663 double tol = removeKnot ( removalIndex, direction, tolerance, reallocate );
1664
1665 if ( tol > maxTol ) maxTol = tol;
1666
1667 removalIndex = knotVect -> findSpan ( param );
1668 }
1669
1670 return maxTol;
1671 }
1672
removeKnot(unsigned rIndex,int dir,double tolerance,bool reallocate)1673 double snlSurface::removeKnot ( unsigned rIndex, int dir, double tolerance, bool reallocate )
1674 {
1675 //! For a Surface: Remove a Knot from Knot Vector and Calculate new Control Points
1676 // ------------------------------------------------------------------------------
1677 //! @param rIndex Knot at index to remove.
1678 //! @param dir Direction to evaluate in. 0 = u, 1 = v.
1679 //! @param tolerance Maximum allowable error. A value of 0.0 means ignore.
1680 //! @param reallocate Reallocate memory for control points.
1681 //!
1682 //! @return Maximum error encountered.
1683
1684 unsigned count, index, lineIndex, offset;
1685 unsigned cDeg, oDeg; // Degree to be processed.
1686 snlKnotVector* cKnts; // Current knots
1687 snlKnotVector* oKnts; // Other knots.
1688 snlCtrlPoint pt1;
1689
1690 if ( dir == SNL_V_DIR )
1691 {
1692 cDeg = degV;
1693 oDeg = degU;
1694 cKnts = knotVectV;
1695 oKnts = knotVectU;
1696 }
1697 else
1698 {
1699 cDeg = degU;
1700 oDeg = degV;
1701 cKnts = knotVectU;
1702 oKnts = knotVectV;
1703 }
1704
1705 knot rParam = cKnts -> val ( rIndex );
1706
1707 // Span knot to be removed belongs to. This will always adjust the removal index to
1708 // point to a non-zero span. ie Multiplicity check.
1709 unsigned rSpan = cKnts -> findSpan ( rParam );
1710
1711 // Find multiplicity of knot at index.
1712 unsigned multi = cKnts -> findMultiplicity ( rSpan );
1713
1714 // Calculate the number of equations.
1715 unsigned numEqns = cDeg - multi + 1;
1716
1717 // Pre calculate alphas.
1718 double* alpha = cKnts -> calcRemovalAlphas ( rSpan );
1719
1720 // Maximum error variable.
1721 double maxError = 0;
1722
1723 // Build temp array to store new array of control points in.
1724 // First and last points are not new.
1725 snlCtrlPoint* tmpPts = new snlCtrlPoint [ numEqns + 1 ];
1726
1727 // Copy control points. Only copy back if the tolerance is ok.
1728
1729 const snlCtrlPoint* ctrlPtsToCopy = controlPoints();
1730
1731 snlCtrlPoint* ctrlPts = new snlCtrlPoint [ sizeU() * sizeV() ];
1732
1733 for ( unsigned index = 0; index < ( sizeU() * sizeV() ); index ++ )
1734 ctrlPts [ index ] = ctrlPtsToCopy [ index ];
1735
1736 // Do for each "line" in direction of removal.
1737 for ( lineIndex = 0; lineIndex < ( oKnts -> size() - oDeg - 1 ); lineIndex ++ )
1738 {
1739 // Seed new points array.
1740 if ( dir == SNL_V_DIR )
1741 tmpPts [ 0 ] = ctrlPts [ ( lineIndex * sizeV() ) + ( rSpan - cDeg - 1 ) ];
1742 else
1743 tmpPts [ 0 ] = ctrlPts [ ( ( rSpan - cDeg - 1 ) * sizeV() ) + lineIndex ];
1744
1745 // Calculate new control points.
1746 for ( count = 1; count <= numEqns; count ++ )
1747 {
1748 index = rSpan - cDeg + count - 1;
1749
1750 // Get ctrl point to process with.
1751 if ( dir == SNL_V_DIR )
1752 {
1753 // U direction.
1754 pt1 = ctrlPts [ ( lineIndex * ( sizeV() ) ) + index ];
1755 }
1756 else
1757 {
1758 // T direction.
1759 pt1 = ctrlPts [ ( index * sizeV() ) + lineIndex ];
1760 }
1761
1762 // Calculate new control point.
1763
1764 tmpPts [ count ] = pt1;
1765 tmpPts [ count ] -= tmpPts [ count - 1 ] * ( 1.0 - alpha [ count - 1 ] );
1766 tmpPts [ count ] /= alpha [ count - 1 ];
1767 }
1768
1769 // Place new points into array.
1770
1771 if ( dir == SNL_V_DIR )
1772 {
1773 // V direction removal.
1774 offset = lineIndex * sizeV();
1775
1776 // Calculate maximum error.
1777
1778 snlCtrlPoint original = ctrlPts [ rSpan - cDeg + offset + numEqns - 1 ];
1779
1780 original.normalise();
1781
1782 tmpPts [ numEqns ].normalise();
1783
1784 snlVector errorVect ( original, tmpPts [ numEqns ] );
1785
1786 double error = errorVect.length();
1787
1788 if ( error > maxError ) maxError = error;
1789
1790 // Copy non-altered control points backward one position at the end of the array.
1791 for ( count = rSpan - multi; count < sizeV(); count ++ )
1792 ctrlPts [ offset + count ] = ctrlPts [ offset + count + 1 ];
1793
1794 // Copy new control points into array.
1795 for ( count = 0; count < ( numEqns - 1 ); count ++ )
1796 {
1797 index = rSpan - cDeg + count + offset;
1798 ctrlPts [ index ] = tmpPts [ count + 1 ];
1799 }
1800 }
1801 else
1802 {
1803 // U direction removal.
1804
1805 // Calculate maximum error.
1806
1807 snlCtrlPoint original = ctrlPts [ ( rSpan - cDeg + numEqns - 1 ) * sizeV() + lineIndex ];
1808
1809 original.normalise();
1810
1811 tmpPts [ numEqns ].normalise();
1812
1813 snlVector errorVect ( original, tmpPts [ numEqns ] );
1814
1815 double error = errorVect.length();
1816
1817 if ( error > maxError ) maxError = error;
1818
1819 // Copy non-altered control points backwards one position at the end of the line.
1820 for ( count = rSpan - multi; count > sizeU(); count ++ )
1821 {
1822 index = count * sizeV() + lineIndex;
1823 ctrlPts [ index ] = ctrlPts [ index + sizeV() ];
1824 }
1825
1826 // Copy new control points into array.
1827 for ( count = 0; count < ( numEqns - 1 ); count ++ )
1828 {
1829 index = ( rSpan - cDeg + count ) * sizeV() + lineIndex;
1830 ctrlPts [ index ] = tmpPts [ count + 1 ];
1831 }
1832 }
1833 }
1834
1835 // If maximum error was under tolerance then go ahead with removal.
1836 if ( maxError < tolerance || tolerance == 0.0 )
1837 {
1838 // Remove knot from knot vector
1839 cKnts -> removeKnot ( rSpan );
1840
1841 ctrlPtNet -> replacePoints ( ctrlPts );
1842
1843 if ( dir == SNL_V_DIR )
1844 ctrlPtNet -> shrinkV();
1845 else
1846 ctrlPtNet -> shrinkU();
1847 }
1848
1849 delete[] tmpPts;
1850 delete[] alpha;
1851
1852 return maxError;
1853 }
1854
createBezierSegments(int dir,int ** numKnotsAdded)1855 unsigned snlSurface::createBezierSegments ( int dir, int** numKnotsAdded )
1856 {
1857 //! Convert surface into Bezier segments.
1858 // -------------------------------------
1859 //! @param dir Direction to process in. u = 0, v = 1.
1860 //! @param numKnotsAdded Pointer to pointer that should hold the number of knots added to
1861 //! each span.
1862 //!
1863 //! @return Number of Bezier segments present.
1864
1865 unsigned cDeg; // Degree to be processed. Current degree, other degree.
1866 unsigned cSize, oSize; // Number of control points in current and other direction.
1867 snlKnotVector* cKnts; // Current knots
1868
1869 if ( dir == SNL_V_DIR )
1870 {
1871 cDeg = degV;
1872 cKnts = knotVectV;
1873 cSize = sizeV();
1874 oSize = sizeU();
1875 }
1876 else
1877 {
1878 cDeg = degU;
1879 cKnts = knotVectU;
1880 cSize = sizeU();
1881 oSize = sizeV();
1882 }
1883
1884 // Find number of non-zero spans.
1885 unsigned numSpans = cKnts -> getNumSpans();
1886
1887 if ( cDeg <= 1 ) return numSpans; // 1st degree curve sections are already Bezier segments.
1888
1889 int* knotsAdded = new int [ numSpans ]; // Last array element is unused and is left here so that a zero length array doesn't occur.
1890
1891 // Find first spans knot index.
1892 unsigned knotIndex = cKnts -> getFirstSpan();
1893
1894 // Resize control points array just once for all knot insertions.
1895
1896 unsigned nextSpan = knotIndex;
1897 unsigned extraKnots = 0;
1898
1899 // Find amount to resize by.
1900 for ( unsigned spanIndex = 1; spanIndex < numSpans; spanIndex ++ )
1901 {
1902 nextSpan = cKnts -> getNextSpan ( nextSpan );
1903
1904 extraKnots += cDeg - ( cKnts -> findMultiplicity ( nextSpan ) );
1905 }
1906
1907 // Append extra control point space to end of current control points.
1908 ctrlPtNet -> appendPointSpace ( oSize * extraKnots );
1909
1910 // *** Create Bezier segments ***
1911
1912 // Find knot index of second knot span.
1913 nextSpan = cKnts -> getNextSpan ( knotIndex );
1914
1915 for ( unsigned spanIndex = 0; spanIndex < numSpans - 1; spanIndex ++ )
1916 {
1917 // Increase multiplicity of span to degree.
1918
1919 unsigned multi = cKnts -> findMultiplicity ( nextSpan );
1920
1921 if ( cDeg - multi > 0 )
1922 {
1923 knot insertParam = cKnts -> val ( nextSpan );
1924
1925 insertKnot ( insertParam, dir, cDeg - multi, false );
1926
1927 // Re-adjust current span index to account for inserted knots.
1928 nextSpan = cKnts -> getNextSpan ( cKnts -> findSpan ( insertParam ) );
1929
1930 // Populate number of knots added array elements.
1931 knotsAdded [ spanIndex ] = cDeg - multi;
1932 }
1933 else
1934 nextSpan = cKnts -> getNextSpan ( nextSpan );
1935 }
1936
1937 if ( numKnotsAdded )
1938 *numKnotsAdded = knotsAdded;
1939 else
1940 delete[] knotsAdded;
1941
1942 return numSpans;
1943 }
1944
createBezierSegments(int * numU,int * numV)1945 void snlSurface::createBezierSegments ( int* numU, int* numV )
1946 {
1947 //! Convert surface into Bezier segments.
1948 // -------------------------------------
1949
1950 int numUSegments = createBezierSegments ( SNL_U_DIR );
1951 int numVSegments = createBezierSegments ( SNL_V_DIR );
1952
1953 if ( numU ) *numU = numUSegments;
1954 if ( numV ) *numV = numVSegments;
1955 }
1956
createConvexBezierSegments(int * numU,int * numV,double sensitivity)1957 void snlSurface::createConvexBezierSegments ( int* numU, int* numV, double sensitivity )
1958 {
1959 //! Convert surface into convex Bezier segments.
1960 // --------------------------------------------
1961 //! @param numU Pointer to variable to return number of U segments in.
1962 //! @param numV Pointer to variable to return number of V segments in.
1963 //! @param sensitivity Maximum concave angle allowed to be considered convex. Used to account
1964 //! for noise in the curvature of a relatively flat section.
1965
1966 int numUSegments = createBezierSegments ( SNL_U_DIR );
1967 int numVSegments = createBezierSegments ( SNL_V_DIR );
1968
1969 // Subdivide U direction until all segments are convex.
1970
1971 int maxDeg = degU > degV ? degU : degV;
1972
1973 snlCtrlPoint** testPoints = new snlCtrlPoint* [ maxDeg + 1 ];
1974
1975 bool keepChecking = true;
1976
1977 while ( keepChecking )
1978 {
1979 keepChecking = false;
1980
1981 int sizeV = ctrlPtNet -> getSizeV();
1982
1983 for ( int segment = 0; segment < numUSegments; segment ++ )
1984 {
1985 int indexU = segment * degU;
1986
1987 bool convex = true;
1988
1989 // Check all constant V lines for concave control point combinations.
1990
1991 for ( int indexV = 0; indexV < sizeV; indexV ++ )
1992 {
1993 ctrlPtNet -> locatePointsU ( indexU, indexV, degU + 1, testPoints );
1994
1995 if ( ! ( ctrlPtNet -> isConvex ( (snlPoint**) testPoints, degU + 1, sensitivity ) ) )
1996 {
1997 convex = false;
1998 break;
1999 }
2000 }
2001
2002 // If any part in U direction is concave then subdivide all patches in range.
2003
2004 if ( ! convex )
2005 {
2006 knot startKnot = knotVectU -> val ( ( segment + 1 ) * degU );
2007 knot endKnot = knotVectU -> val ( ( segment + 2 ) * degU );
2008
2009 // Insert multiple knots into middle of patch to create two new Bezier segments.
2010
2011 insertKnot ( ( ( endKnot - startKnot ) / 2.0 ) + startKnot, SNL_U_DIR, degU, true );
2012
2013 numUSegments ++;
2014 segment ++;
2015
2016 keepChecking = true;
2017 }
2018 }
2019 }
2020
2021 // Subdivide V direction until all segments are convex.
2022
2023 keepChecking = true;
2024
2025 while ( keepChecking )
2026 {
2027 keepChecking = false;
2028
2029 int sizeU = ctrlPtNet -> getSizeU();
2030
2031 for ( int segment = 0; segment < numVSegments; segment ++ )
2032 {
2033 int indexV = segment * degV;
2034
2035 bool convex = true;
2036
2037 // Check all U lines for concave control point combinations.
2038
2039 for ( int indexU = 0; indexU < sizeU; indexU ++ )
2040 {
2041 ctrlPtNet -> locatePointsV ( indexU, indexV, degV + 1, testPoints );
2042
2043 if ( ! ( ctrlPtNet -> isConvex ( (snlPoint**) testPoints, degV + 1, sensitivity ) ) )
2044 {
2045 convex = false;
2046 break;
2047 }
2048 }
2049
2050 // If any part in V direction is concave then subdivide all patches in range.
2051
2052 if ( ! convex )
2053 {
2054 knot startKnot = knotVectV -> val ( ( segment + 1 ) * degV );
2055 knot endKnot = knotVectV -> val ( ( segment + 2 ) * degV );
2056
2057 // Insert multiple knots into middle of patch to create two new Bezier segments.
2058
2059 insertKnot ( ( ( endKnot - startKnot ) / 2.0 ) + startKnot, SNL_V_DIR, degV, true );
2060
2061 numVSegments ++;
2062 segment ++;
2063
2064 keepChecking = true;
2065 }
2066 }
2067 }
2068
2069 if ( numU ) *numU = numUSegments;
2070 if ( numV ) *numV = numVSegments;
2071
2072 delete[] testPoints;
2073 }
2074
elevateDegree(int direction,int byDegree)2075 void snlSurface::elevateDegree ( int direction, int byDegree )
2076 {
2077 //! Elevate degree of surface.
2078 // --------------------------
2079 //! @param direction Parametric direction to elevate degree in. ( enum parametricDirections ).
2080 //! @param byDegree Number of degrees to elevate direction by.
2081 //! Convert curve into Bezier segments.
2082
2083 int* addedKnots = 0;
2084
2085 unsigned numSegments = createBezierSegments ( direction, &addedKnots );
2086
2087 // Grow control point net.
2088
2089 if ( direction == SNL_U_DIR )
2090 ctrlPtNet -> growU ( numSegments * byDegree, true );
2091 else
2092 ctrlPtNet -> growV ( numSegments * byDegree, true );
2093
2094 // Setup direction specific variables.
2095
2096 int cDeg, numLines, lineSize;
2097 snlKnotVector* cKnotVect;
2098
2099 if ( direction == SNL_U_DIR )
2100 {
2101 cDeg = degU;
2102 numLines = ctrlPtNet -> getSizeV();
2103 lineSize = ctrlPtNet -> getSizeU();
2104 cKnotVect = knotVectU;
2105 }
2106 else
2107 {
2108 cDeg = degV;
2109 numLines = ctrlPtNet -> getSizeU();
2110 lineSize = ctrlPtNet -> getSizeV();
2111 cKnotVect = knotVectV;
2112 }
2113
2114 // Elevate degree of Bezier segments.
2115
2116 int segmentSize = cDeg + 1;
2117 int newSegmentSize = segmentSize + byDegree;
2118
2119 snlCtrlPoint* tmpPts = new snlCtrlPoint [ newSegmentSize ];
2120
2121 snlCtrlPoint** linePtPtrs = new snlCtrlPoint* [ lineSize ];
2122
2123 snlCtrlPoint* linePts = new snlCtrlPoint [ lineSize ];
2124
2125 for ( int lineIndex = 0; lineIndex < numLines; lineIndex ++ )
2126 {
2127 // Generate new points per segment.
2128
2129 int ptsIndex = 0;
2130
2131 unsigned spanIndex = cDeg * 2;
2132
2133 if ( direction == SNL_U_DIR )
2134 {
2135 // Get line in U direction.
2136 ctrlPtNet -> locatePointsU ( 0, lineIndex, lineSize, linePtPtrs );
2137 }
2138 else
2139 {
2140 // Get line in V direction.
2141 ctrlPtNet -> locatePointsV ( lineIndex, 0, lineSize, linePtPtrs );
2142 }
2143
2144 // Elevate segments.
2145
2146 for ( unsigned segment = 0; segment < numSegments; segment ++ )
2147 {
2148 // Populate control points to process for line.
2149
2150 int segmentStartIndex = segment * ( segmentSize - 1 );
2151
2152 for ( int index = 0; index < segmentSize; index ++ )
2153 {
2154 linePts [ ptsIndex + index ] = *( linePtPtrs [ segmentStartIndex + index ] );
2155 }
2156
2157 // Process segment.
2158
2159 snlCurve::elevateBezierSegmentPointsDegree ( cDeg, byDegree, linePts + ptsIndex, tmpPts );
2160
2161 // Replace points in temp control point array. Index 0 remains unchanged.
2162 for ( int index = 1; index < newSegmentSize; index ++ )
2163 linePts [ ptsIndex + index ] = tmpPts [ index ];
2164
2165 ptsIndex += cDeg + byDegree;
2166
2167 if ( lineIndex == 0 )
2168 {
2169 // Add knots to knot vector. This is only done _once_.
2170 cKnotVect -> increaseMultiplicity ( spanIndex, byDegree );
2171 }
2172
2173 spanIndex += cDeg + byDegree;
2174 }
2175
2176 // Copy temp points into control point net.
2177
2178 for ( int index = 0; index < lineSize; index ++ )
2179 *( linePtPtrs [ index ] ) = linePts [ index ];
2180 }
2181
2182 // Make sure start clamp is of degree + 1 multiplicity.
2183
2184 cKnotVect -> increaseMultiplicity ( cDeg, byDegree );
2185
2186 // Increase degree indicator variables
2187
2188 cDeg += byDegree;
2189
2190 cKnotVect -> degree ( cDeg );
2191
2192 if ( direction == SNL_U_DIR )
2193 degU = cDeg;
2194 else
2195 degV = cDeg;
2196
2197 // Remove number of knots that were added during knot insertion.
2198 // No knots were added to degree one sections.
2199
2200 if ( cDeg - byDegree > 1 )
2201 {
2202 unsigned spanIndex = cKnotVect -> getFirstSpan();
2203
2204 spanIndex += cDeg;
2205
2206 for ( unsigned segment = 0; segment < numSegments - 1; segment ++ )
2207 {
2208 removeKnots ( addedKnots [ segment ], spanIndex, direction, 0.0, true );
2209
2210 spanIndex += cDeg - addedKnots [ segment ];
2211 }
2212 }
2213
2214 // Clean up.
2215
2216 if ( addedKnots ) delete[] addedKnots;
2217 delete[] tmpPts;
2218 delete[] linePtPtrs;
2219 delete[] linePts;
2220 }
2221
reduceDegree(int dir,unsigned numDeg,double tolerance)2222 double snlSurface::reduceDegree ( int dir, unsigned numDeg, double tolerance )
2223 {
2224 //! Reduce Surface Degree
2225 // ---------------------
2226 //! @param dir Direction to evaluate in. 0 = u, 1 = v.
2227 //! @param numDeg Number of degrees to reduce by.
2228 //! @param tolerance Maximum error. A value of 0.0 means no tolerance specified.
2229 //!
2230 //! @return Maximum error encountered.
2231 //!
2232 //! @par Notes:
2233 //! This function has not been optimised.
2234
2235 unsigned cDeg, oDeg; // Degree to be processed. Current degree, other degree.
2236 unsigned cSize, oSize; // Number of control points in current and other direction.
2237 snlKnotVector* cKnts; // Current knots
2238 snlKnotVector* oKnts; // Other knots.
2239
2240 double maxError = 0.0;
2241
2242 // Save knots and control points of original surface.
2243
2244 snlCtrlPointNetSurface* ctrlPtNetCopy = new snlCtrlPointNetSurface ( *ctrlPtNet );
2245
2246 snlKnotVector* knotVectUCopy = new snlKnotVector ( *knotVectU );
2247 snlKnotVector* knotVectVCopy = new snlKnotVector ( *knotVectV );
2248
2249 // Convert into Bezier segments.
2250 unsigned numSpans = createBezierSegments ( dir );
2251
2252 if ( dir == SNL_V_DIR )
2253 {
2254 cDeg = degV;
2255 oDeg = degU;
2256 cKnts = knotVectV;
2257 oKnts = knotVectU;
2258 cSize = sizeV();
2259 oSize = sizeU();
2260 }
2261 else
2262 {
2263 cDeg = degU;
2264 oDeg = degV;
2265 cKnts = knotVectU;
2266 oKnts = knotVectV;
2267 cSize = sizeU();
2268 oSize = sizeV();
2269 }
2270
2271 // *** Reduce degree of Bezier segments ***
2272
2273 snlCtrlPoint* ctrlPts = ctrlPtNet -> getCtrlPtsPtr();
2274
2275 for ( unsigned count = 0; count < numDeg; count ++ )
2276 {
2277 unsigned rDeg = cDeg - count; // Current degree during reduction.
2278
2279 // Pre-calculate alphas.
2280 double* alpha = new double [ rDeg - 1 ];
2281
2282 for ( unsigned index = 1; index < rDeg; index ++ )
2283 alpha [ index - 1 ] = (double) index / (double) rDeg;
2284
2285 for ( unsigned patchIndex = 0; patchIndex < numSpans; patchIndex ++ )
2286 {
2287 // Process up each "line" in the direction to be reduced.
2288 for ( unsigned lineIndex = 0; lineIndex < oSize; lineIndex ++ )
2289 {
2290 // Reduce patch.
2291
2292 snlCtrlPoint* cPnt; // Current point.
2293 snlCtrlPoint* pPnt; // Previous point.
2294
2295 // Get starting and previous point's pointer.
2296
2297 if ( dir == SNL_V_DIR )
2298 {
2299 // V Direction.
2300 pPnt = ctrlPts + (unsigned)( ( lineIndex * cSize ) + ( patchIndex * rDeg ) );
2301 cPnt = pPnt + 1;
2302 }
2303 else
2304 {
2305 // U Direction.
2306 pPnt = ctrlPts + (unsigned)( ( patchIndex * rDeg * oSize ) + lineIndex );
2307 cPnt = pPnt + sizeV();
2308 }
2309
2310 snlCtrlPoint newPoint;
2311
2312 // Caculate new internal patch points. ie First and last points aren't modified.
2313
2314 for ( unsigned pIndex = 1; pIndex < rDeg; pIndex ++ )
2315 {
2316 if ( pIndex > 1 )
2317 *pPnt = newPoint;
2318
2319 newPoint.x ( ( cPnt -> x() - ( ( pPnt -> x() ) * alpha [ pIndex - 1 ] ) )
2320 / ( 1 - alpha [ pIndex - 1 ] ) );
2321 newPoint.y ( ( cPnt -> y() - ( ( pPnt -> y() ) * alpha [ pIndex - 1 ] ) )
2322 / ( 1 - alpha [ pIndex - 1 ] ) );
2323 newPoint.z ( ( cPnt -> z() - ( ( pPnt -> z() ) * alpha [ pIndex - 1 ] ) )
2324 / ( 1 - alpha [ pIndex - 1 ] ) );
2325 newPoint.w ( ( cPnt -> w() - ( ( pPnt -> w() ) * alpha [ pIndex - 1 ] ) )
2326 / ( 1 - alpha [ pIndex - 1 ] ) );
2327
2328 // Find next points.
2329 pPnt = cPnt;
2330
2331 if ( dir == SNL_V_DIR )
2332 // V direction.
2333 cPnt ++;
2334 else
2335 // U direction
2336 cPnt += sizeV();
2337 }
2338
2339 // Calculate error.
2340
2341 snlCtrlPoint cPntCopy ( *cPnt );
2342 cPntCopy.normalise();
2343
2344 newPoint.normalise();
2345
2346 double cError = ( snlVector ( cPntCopy, newPoint ).length() );
2347
2348 if ( cError > maxError ) maxError = cError;
2349 }
2350
2351 // Adjust knot vector.
2352 cKnts -> removeKnot ( ( patchIndex * ( rDeg - 1 ) ) + 1 );
2353 }
2354
2355 // Reduce end clamp multiplicity.
2356 cKnts -> removeKnot ( ( cKnts -> size() ) - 1 );
2357
2358 // Rearrange array to facilitate removal of redundant control points.
2359
2360 snlCtrlPoint* copyFrom = ctrlPts;
2361 snlCtrlPoint* copyTo = ctrlPts;
2362
2363 if ( dir == SNL_V_DIR )
2364 {
2365 for ( unsigned lineIndex = 0; lineIndex < oSize; lineIndex ++ )
2366 {
2367 for ( unsigned patchIndex = 0; patchIndex < numSpans; patchIndex ++ )
2368 {
2369 for ( unsigned count = 0; count < rDeg - 1; count ++ )
2370 {
2371 *copyTo = *copyFrom;
2372 copyTo++;
2373 copyFrom++;
2374 }
2375
2376 copyFrom++;
2377 }
2378
2379 *copyTo = *copyFrom;
2380 copyTo++;
2381 copyFrom++;
2382 }
2383 }
2384 else
2385 {
2386 // U direction.
2387
2388 for ( unsigned patchIndex = 0; patchIndex < numSpans; patchIndex ++ )
2389 {
2390 for ( unsigned lineIndex = 0; lineIndex < rDeg - 1; lineIndex ++ )
2391 {
2392 for ( unsigned count = 0; count < oSize; count ++ )
2393 {
2394 *copyTo = *copyFrom;
2395 copyTo++;
2396 copyFrom++;
2397 }
2398 }
2399
2400 copyFrom += oSize; // Skip a line.
2401 }
2402
2403 // Copy last line.
2404
2405 for ( unsigned count = 0; count < oSize; count ++ )
2406 {
2407 *copyTo = *copyFrom;
2408 copyTo++;
2409 copyFrom++;
2410 }
2411 }
2412
2413 // Adjust current new size to reflect control point removals.
2414 cSize -= numSpans;
2415
2416 // Clean up.
2417 delete[] alpha;
2418 }
2419
2420 // Set control point net to correct size.
2421 ctrlPtNet -> truncatePointSpace ( numDeg * numSpans );
2422
2423 if ( dir == SNL_V_DIR )
2424 {
2425 ctrlPtNet -> setSizeV ( cSize );
2426 degV -= numDeg;
2427 }
2428 else
2429 {
2430 ctrlPtNet -> setSizeU ( cSize );
2431 degU -= numDeg;
2432 }
2433
2434 // Remove additional knots only if under tolerance.
2435
2436 if ( ( tolerance != 0.0 ) && ( tolerance > maxError ) )
2437 {
2438 // !@#$ Simplify surface here.
2439 }
2440
2441 if ( ( maxError < tolerance ) || ( tolerance == 0.0 ) )
2442 {
2443 // Can keep new knots and control points so delete original copies.
2444
2445 delete ctrlPtNetCopy;
2446 delete knotVectUCopy;
2447 delete knotVectVCopy;
2448 }
2449 else
2450 {
2451 // Reinstate old knots and control points.
2452
2453 delete ctrlPtNet;
2454 delete knotVectV;
2455 delete knotVectU;
2456
2457 ctrlPtNet = ctrlPtNetCopy;
2458 knotVectV = knotVectVCopy;
2459 knotVectU = knotVectUCopy;
2460
2461 if ( dir == SNL_V_DIR )
2462 degV += numDeg;
2463 else
2464 degU += numDeg;
2465 }
2466
2467 return maxError;
2468 }
2469
refine(double tolerance)2470 void snlSurface::refine ( double tolerance )
2471 {
2472 //! Refine surface control point net until tolerance is achieved.
2473 // -------------------------------------------------------------
2474 //! @par Notes:
2475 //! Guarantees that if a span is approximated by it's end points then the error will be no
2476 //! greater than tolerance.
2477
2478 // Refine in individual parametric directions on pass in each direction at a time.
2479
2480 bool tolOk = false;
2481
2482 while ( ! tolOk )
2483 {
2484 tolOk = true;
2485
2486 if ( ! refineHull_U ( tolerance, true ) ) tolOk = false;
2487 if ( ! refineHull_V ( tolerance, true ) ) tolOk = false;
2488 }
2489 }
2490
refineHull_UV(double tolerance)2491 void snlSurface::refineHull_UV ( double tolerance )
2492 {
2493 //! Refine control point net until it is guaranteed be within tolerance to surface.
2494 // -------------------------------------------------------------------------------
2495 //! @par Notes:
2496 //! Processes rectangular sections which is more accurate than doing each parametric
2497 //! direction in turn. This is also results in this funtion being VERY slow and extremely
2498 //! memory hungry.
2499
2500 if ( degU < 2 && degV < 2 ) return;
2501
2502 if ( degU < 2 )
2503 {
2504 refineHull_V ( tolerance );
2505 return;
2506 }
2507
2508 if ( degV < 2 )
2509 {
2510 refineHull_U ( tolerance );
2511 return;
2512 }
2513
2514 bool tolOk = false; // Tolerance is okay.
2515
2516 while ( ! tolOk )
2517 {
2518 tolOk = true;
2519
2520 // Control point net sizes must be obtained at each loop as they may change during the loop.
2521
2522 // Only process non-zero knot vector spans.
2523
2524 int numSpansU = knotVectU -> getNumSpans();
2525 int numSpansV = knotVectV -> getNumSpans();
2526
2527 bool* spanSubU = new bool [ numSpansU ]; // Span subdivision indicators.
2528 bool* spanSubV = new bool [ numSpansV ]; // If true then subdivide span.
2529
2530 knot* spanSubUVal = new double [ numSpansU ]; // Knot insertion value for span.
2531 knot* spanSubVVal = new double [ numSpansV ];
2532
2533 for ( int index = 0; index < numSpansU; index ++ )
2534 spanSubU [ index ] = false;
2535
2536 for ( int index = 0; index < numSpansV; index ++ )
2537 spanSubV [ index ] = false;
2538
2539 int spanU = 0;
2540
2541 int indexU = knotVectU -> getFirstSpan();
2542
2543 while ( indexU )
2544 {
2545 int indexV = knotVectV -> getFirstSpan();
2546
2547 int spanV = 0;
2548
2549 while ( indexV )
2550 {
2551 // Don't check rectangular span if both and U and V are already marked for subdivision.
2552
2553 if ( ! ( spanSubU [ spanU ] && spanSubV [ spanV ] ) )
2554 {
2555 double flatness = ctrlPtNet -> calcFlatness ( indexU - degU, indexV - degV, degU + 1, degV + 1 );
2556
2557 if ( flatness > tolerance )
2558 {
2559 if ( ! spanSubU [ spanU ] )
2560 {
2561 spanSubUVal [ spanU ] = ( ( knotVectU -> val ( indexU + 1 )
2562 - knotVectU -> val ( indexU ) ) / 2 )
2563 + knotVectU -> val ( indexU );
2564 }
2565
2566 if ( ! spanSubV [ spanV ] )
2567 {
2568 spanSubVVal [ spanV ] = ( ( knotVectV -> val ( indexV + 1 )
2569 - knotVectV -> val ( indexV ) ) / 2 )
2570 + knotVectV -> val ( indexV );
2571 }
2572
2573 spanSubU [ spanU ] = true;
2574 spanSubV [ spanV ] = true;
2575
2576 tolOk = false;
2577 }
2578 }
2579
2580 indexV = knotVectV -> getNextSpan ( indexV );
2581
2582 spanV ++;
2583 }
2584
2585 indexU = knotVectU -> getNextSpan ( indexU );
2586
2587 spanU ++;
2588 }
2589
2590 // Insert Knots where needed.
2591
2592 // Insert U knots.
2593
2594 for ( int spanIndex = 0; spanIndex < numSpansU; spanIndex ++ )
2595 {
2596 if ( spanSubU [ spanIndex ] )
2597 insertKnot ( spanSubUVal [ spanIndex ], SNL_U_DIR );
2598 }
2599
2600 // Insert V knots.
2601
2602 for ( int spanIndex = 0; spanIndex < numSpansV; spanIndex ++ )
2603 {
2604 if ( spanSubV [ spanIndex ] )
2605 insertKnot ( spanSubVVal [ spanIndex ], SNL_V_DIR );
2606 }
2607
2608 delete[] spanSubU;
2609 delete[] spanSubV;
2610
2611 delete[] spanSubUVal;
2612 delete[] spanSubVVal;
2613 }
2614 }
2615
refineHull_U(double tolerance,bool singlePass)2616 bool snlSurface::refineHull_U ( double tolerance, bool singlePass )
2617 {
2618 //! Refine control point net in U direction until it is
2619 //! guaranteed be within tolerance to surface.
2620 // ---------------------------------------------------
2621 //! @param tolerance Convex Hull of surface in U direction must be within tolerance of surface.
2622 //! @param singlePass Only do a single refinement pass.
2623 //!
2624 //! @return If tolerance is was okay. Should be true unless single pass is specified.
2625
2626 if ( degU < 2 ) return true;
2627
2628 bool tolOk = false; // Tolerance is okay.
2629
2630 while ( ! tolOk )
2631 {
2632 tolOk = true;
2633
2634 // Only process non-zero knot vector spans.
2635
2636 int numSpans = knotVectU -> getNumSpans();
2637
2638 bool* spanSub = new bool [ numSpans ]; // Span subdivision indicators. If true then subdivide span.
2639
2640 knot* spanSubVal = new double [ numSpans ]; // Knot insertion value for span.
2641
2642 for ( int index = 0; index < numSpans; index ++ )
2643 spanSub [ index ] = false;
2644
2645 for ( int indexV = 0; (unsigned) indexV < ( ctrlPtNet -> getSizeV() ); indexV ++ )
2646 {
2647 int indexU = knotVectU -> getFirstSpan();
2648
2649 int spanNum = 0;
2650
2651 while ( indexU )
2652 {
2653 // Test for flatness
2654
2655 double flatness;
2656
2657 flatness = ctrlPtNet -> calcFlatnessU ( indexU - degU, indexV, degU + 1, false );
2658
2659 if ( flatness > tolerance )
2660 {
2661 // Insert knot into surface. Half way between existing knots.
2662
2663 int insertIndex = indexU;
2664
2665 knot insertParam = ( ( knotVectU -> val ( insertIndex + 1 )
2666 - knotVectU -> val ( insertIndex ) ) / 2 )
2667 + knotVectU -> val ( insertIndex );
2668
2669 spanSubVal [ spanNum ] = insertParam;
2670
2671 spanSub [ spanNum ] = true;
2672
2673 tolOk = false;
2674 }
2675
2676 indexU = knotVectU -> getNextSpan ( indexU );
2677
2678 spanNum ++;
2679 }
2680 }
2681
2682 // Insert knots where needed to subdivide span.
2683
2684 for ( int spanNum = 0; spanNum < numSpans; spanNum ++ )
2685 {
2686 if ( spanSub [ spanNum ] )
2687 insertKnot ( spanSubVal [ spanNum ], SNL_U_DIR );
2688 }
2689
2690 delete[] spanSub;
2691 delete[] spanSubVal;
2692
2693 if ( singlePass ) break;
2694 }
2695
2696 return tolOk;
2697 }
2698
refineHull_V(double tolerance,bool singlePass)2699 bool snlSurface::refineHull_V ( double tolerance, bool singlePass )
2700 {
2701 //! Refine control point net in V direction until it is
2702 //! guaranteed be within tolerance to surface.
2703 // ---------------------------------------------------
2704
2705 if ( degV < 2 ) return true;
2706
2707 bool tolOk = false;
2708
2709 while ( ! tolOk )
2710 {
2711 tolOk = true;
2712
2713 // Only process non-zero knot vector spans.
2714
2715 int numSpans = knotVectV -> getNumSpans();
2716
2717 bool* spanSub = new bool [ numSpans ]; // Span subdivision indicators. If true then subdivide span.
2718
2719 knot* spanSubVal = new double [ numSpans ]; // Knot insertion value for span.
2720
2721 for ( int index = 0; index < numSpans; index ++ )
2722 spanSub [ index ] = false;
2723
2724 for ( int indexU = 0; (unsigned) indexU < ( ctrlPtNet -> getSizeU() ); indexU ++ )
2725 {
2726 int indexV = knotVectV -> getFirstSpan();
2727
2728 int spanNum = 0;
2729
2730 while ( indexV )
2731 {
2732 // Test for flatness
2733
2734 double flatness;
2735
2736 flatness = ctrlPtNet -> calcFlatnessV ( indexU, indexV - degV, degV + 1, false );
2737
2738 if ( flatness > tolerance )
2739 {
2740 // Insert knot into surface. Half way between existing knots.
2741
2742 int insertIndex = indexV;
2743
2744 knot insertParam = ( ( knotVectV -> val ( insertIndex + 1 )
2745 - knotVectV -> val ( insertIndex ) ) / 2 )
2746 + knotVectV -> val ( insertIndex );
2747
2748 spanSubVal [ spanNum ] = insertParam;
2749
2750 spanSub [ spanNum ] = true;
2751
2752 tolOk = false;
2753 }
2754
2755 indexV = knotVectV -> getNextSpan ( indexV );
2756
2757 spanNum ++;
2758 }
2759 }
2760
2761 // Insert knots where needed to subdivide span.
2762
2763 for ( int spanNum = 0; spanNum < numSpans; spanNum ++ )
2764 {
2765 if ( spanSub [ spanNum ] )
2766 insertKnot ( spanSubVal [ spanNum ], SNL_V_DIR );
2767 }
2768
2769 delete[] spanSub;
2770 delete[] spanSubVal;
2771
2772 if ( singlePass ) break;
2773 }
2774
2775 return tolOk;
2776 }
2777
refineHullBezier(double tolerance)2778 void snlSurface::refineHullBezier ( double tolerance )
2779 {
2780 //! Refine control point net using Bezier patch subdivision.
2781 // --------------------------------------------------------
2782 //! @param tolerance Control point net should be within this tolerance of surface.
2783
2784 // Break surface into Bezier segments ( patches ).
2785
2786 int numUSegments;
2787 int numVSegments;
2788
2789 createBezierSegments ( &numUSegments, &numVSegments );
2790
2791 // Walk through each segment and calculate it's flatness.
2792
2793 bool tolOk = false; // Tolerance has been achieved.
2794
2795 while ( ! tolOk )
2796 {
2797 tolOk = true;
2798
2799 numUSegments = knotVectU -> getNumSpans();
2800 numVSegments = knotVectV -> getNumSpans();
2801
2802 bool* splitSpanU = new bool [ numUSegments ]; // Split knot span corresponding to array index if true.
2803 bool* splitSpanV = new bool [ numVSegments ];
2804
2805 knot* splitKnotU = new knot [ numUSegments ]; // Knot value within span that span will be split at.
2806 knot* splitKnotV = new knot [ numVSegments ];
2807
2808 for ( int index = 0; index < numUSegments; index ++ )
2809 splitSpanU [ index ] = false;
2810
2811 for ( int index = 0; index < numVSegments; index ++ )
2812 splitSpanV [ index ] = false;
2813
2814 // Step through each segment and calculate control point distance to surface.
2815
2816 for ( int segU = 0; segU < numUSegments; segU ++ )
2817 {
2818 for ( int segV = 0; segV < numVSegments; segV ++ )
2819 {
2820
2821 if ( segU == 159 && segV == 2 )
2822 {
2823 cout << "stop\n";
2824 }
2825 // Don't test a segment if it has already being split.
2826
2827 if ( ! splitSpanU [ segU ] || ! splitSpanV [ segV ] )
2828 {
2829 double flatness = ctrlPtNet -> calcFlatness ( segU * degU, segV * degV, degU + 1, degV + 1 );
2830
2831 if ( flatness > tolerance )
2832 {
2833 tolOk = false;
2834 splitSpanU [ segU ] = true;
2835 splitSpanV [ segV ] = true;
2836 cout << "SegU: " << segU << " SegV: " << segV << " Flatness: " << flatness << "\n";
2837 }
2838 }
2839 }
2840 }
2841
2842 if ( ! tolOk )
2843 {
2844 // Split spans. First find parameters to insert.
2845
2846 int spanIndex = knotVectU -> getFirstSpan(); // Current span index.
2847
2848 for ( int spanNum = 0; spanNum < numUSegments; spanNum ++ )
2849 {
2850 if ( splitSpanU [ spanNum ] )
2851 {
2852 splitKnotU [ spanNum ] = ( ( knotVectU -> val ( spanIndex + 1 )
2853 - knotVectU -> val ( spanIndex ) ) / 2 )
2854 + knotVectU -> val ( spanIndex );
2855 }
2856
2857 spanIndex = knotVectU -> getNextSpan ( spanIndex );
2858 }
2859
2860 spanIndex = knotVectV -> getFirstSpan();
2861
2862 for ( int spanNum = 0; spanNum < numVSegments; spanNum ++ )
2863 {
2864 if ( splitSpanV [ spanNum ] )
2865 {
2866 splitKnotV [ spanNum ] = ( ( knotVectV -> val ( spanIndex + 1 )
2867 - knotVectV -> val ( spanIndex ) ) / 2 )
2868 + knotVectV -> val ( spanIndex );
2869 }
2870
2871 spanIndex = knotVectV -> getNextSpan ( spanIndex );
2872 }
2873
2874 // Insert knots into each vector.
2875
2876 for ( int spanNum = 0; spanNum < numUSegments; spanNum ++ )
2877 {
2878 if ( splitSpanU [ spanNum ] )
2879 insertKnot ( splitKnotU [ spanNum ], SNL_U_DIR, degU );
2880 }
2881
2882 for ( int spanNum = 0; spanNum < numVSegments; spanNum ++ )
2883 {
2884 if ( splitSpanV [ spanNum ] )
2885 insertKnot ( splitKnotV [ spanNum ], SNL_V_DIR, degV );
2886 }
2887 }
2888
2889 // Number of Bezier segments may have increased. If arrays are to be used again their sizes
2890 // will be invalid.
2891
2892 delete[] splitSpanU;
2893 delete[] splitSpanV;
2894
2895 delete[] splitKnotU;
2896 delete[] splitKnotV;
2897 }
2898 }
2899
maxCurvatureU()2900 double snlSurface::maxCurvatureU()
2901 {
2902 //! Return maximum curvature of surface in U direction.
2903 // ---------------------------------------------------
2904 //! @return Value between 0 and PI.
2905
2906 return ctrlPtNet -> maxCurvatureU();
2907 }
2908
maxCurvatureV()2909 double snlSurface::maxCurvatureV()
2910 {
2911 //! Return maximum curvature of surface in V direction.
2912 // ---------------------------------------------------
2913 //! @return Value between 0 and PI.
2914
2915 return ctrlPtNet -> maxCurvatureV();
2916 }
2917
calcNormal(knot paramU,knot paramV,snlPoint * evalPt)2918 snlVector snlSurface::calcNormal ( knot paramU, knot paramV, snlPoint* evalPt )
2919 {
2920 //! Calculate surface normal.
2921 // -------------------------
2922 //! @param paramU U parameter of location to generate normal from.
2923 //! @param paramV V parameter of location to generate normal from.
2924 //! @param evalPt Pointer to point that holds point evaluted at paramter. Optional.
2925 //!
2926 //! @return Pointer to normal to surface. Caller owns pointer.
2927
2928 snlVector velocityU;
2929 snlVector velocityV;
2930 snlPoint evalPoint;
2931
2932 if ( ! evalPt )
2933 velocities ( paramU, paramV, evalPoint, velocityU, velocityV );
2934 else
2935 velocities ( paramU, paramV, *evalPt, velocityU, velocityV );
2936
2937 snlVector normal;
2938
2939 normal.crossProduct ( velocityU, velocityV );
2940
2941 return normal;
2942 }
2943
findClosestCtrlPt(snlPoint * points,int numPoints)2944 snlSCtrlPtLocn* snlSurface::findClosestCtrlPt ( snlPoint* points, int numPoints )
2945 {
2946 //! Find closest control point to a given point.
2947 // --------------------------------------------
2948 //! @param points Array of points to process.
2949 //! @param numPoints Number of points in array to process.
2950 //!
2951 //! @return Array of found control points that best match given points.
2952 //! Returned array indexes correspond to given array indexes.
2953 //! Caller owns array and should delete it once no longer needed.
2954
2955 int index;
2956
2957 // Initialise return array.
2958
2959 snlSCtrlPtLocn* retArray = new snlSCtrlPtLocn [ numPoints ];
2960
2961 double maxDouble = DBL_MAX;
2962
2963 for ( index = 0; index < numPoints; index ++ )
2964 retArray [ index ].dist = maxDouble; // Make this very large.
2965
2966 const snlCtrlPoint* ctrlPts = controlPoints();
2967
2968 // Check test points against control points.
2969
2970 int uSize = sizeU();
2971 int vSize = sizeV();
2972
2973 index = 0;
2974
2975 double dist;
2976
2977 for ( int indexU = 0; indexU < uSize; indexU ++ )
2978 {
2979 for ( int indexV = 0; indexV < vSize; indexV ++ )
2980 {
2981 for ( int ptIndex = 0; ptIndex < numPoints; ptIndex ++ )
2982 {
2983 // Calculate distance to point.
2984
2985 dist = ctrlPts [ index ].distSqrd ( points [ ptIndex ] );
2986
2987 if ( retArray [ ptIndex ].dist > dist )
2988 {
2989 retArray [ ptIndex ].dist = dist;
2990 retArray [ ptIndex ].uIndex = indexU;
2991 retArray [ ptIndex ].vIndex = indexV;
2992 }
2993 }
2994
2995 index ++;
2996 }
2997 }
2998
2999 return retArray;
3000 }
3001
extractEdge(int edge)3002 snlCurve* snlSurface::extractEdge ( int edge )
3003 {
3004 //! Extract surface edge as curve.
3005 // ------------------------------
3006 //! @param edge enum surfaceEdges value that indicates which edge to return.
3007
3008 snlCtrlPoint** ctrlPtPtrs;
3009 snlKnotVector* knotVect;
3010
3011 int size;
3012
3013 switch ( edge )
3014 {
3015 case SNL_EDGE_UMIN:
3016
3017 size = ctrlPtNet -> getSizeV();
3018 ctrlPtPtrs = new snlCtrlPoint* [ size ];
3019 ctrlPtNet -> locatePointsV ( 0, 0, size, ctrlPtPtrs );
3020 knotVect = new snlKnotVector ( *knotVectV );
3021
3022 break;
3023
3024 case SNL_EDGE_VMIN:
3025
3026 size = ctrlPtNet -> getSizeU();
3027 ctrlPtPtrs = new snlCtrlPoint* [ size ];
3028 ctrlPtNet -> locatePointsU ( 0, 0, size, ctrlPtPtrs );
3029 knotVect = new snlKnotVector ( *knotVectU );
3030
3031 break;
3032
3033 case SNL_EDGE_UMAX:
3034
3035 size = ctrlPtNet -> getSizeV();
3036 ctrlPtPtrs = new snlCtrlPoint* [ size ];
3037 ctrlPtNet -> locatePointsV ( ctrlPtNet -> getSizeU() - 1, 0, size, ctrlPtPtrs );
3038 knotVect = new snlKnotVector ( *knotVectV );
3039
3040 break;
3041
3042 case SNL_EDGE_VMAX:
3043
3044 size = ctrlPtNet -> getSizeU();
3045 ctrlPtPtrs = new snlCtrlPoint* [ size ];
3046 ctrlPtNet -> locatePointsU ( 0, ctrlPtNet -> getSizeV() - 1, size, ctrlPtPtrs );
3047 knotVect = new snlKnotVector ( *knotVectU );
3048
3049 break;
3050 };
3051
3052 // Generate new control point array.
3053
3054 snlCtrlPoint* ctrlPts = new snlCtrlPoint [ size ];
3055
3056 for ( int index = 0; index < size; index ++ )
3057 ctrlPts [ index ] = * ctrlPtPtrs [ index ];
3058
3059 delete[] ctrlPtPtrs;
3060
3061 // Generate curve to return.
3062
3063 return new snlCurve ( size, ctrlPts, knotVect );
3064 }
3065
fillet(int edge,snlVector & frontFaceNormal,snlSurface & surface2,snlVector & frontFaceNormal2,double tolerance,double radius,bool trim1,bool trim2)3066 snlSurface* snlSurface::fillet ( int edge, snlVector& frontFaceNormal,
3067 snlSurface& surface2, snlVector& frontFaceNormal2,
3068 double tolerance, double radius, bool trim1, bool trim2 )
3069 {
3070 //! Create a fillet between this and another surface.
3071 // -------------------------------------------------
3072 //! @param edge Edge of this surface to fillet.
3073 //! @param frontFaceNormal Normal that defines side of this surface to place fillet on.
3074 //! @param surface2 surface to fillet.
3075 //! @param frontFaceNormal2 Normal that defines side of surface to place fillet on.
3076 //! @param tolerance Tolerance to surfaces fillet should comply with.
3077 //! @param radius Radius of fillet.
3078 //! @param trim1 Trim surface1 ( controlling surface ).
3079 //! @param trim2 Trim surface2 ( matching surface ).
3080 //!
3081 //! @par Notes:
3082 //! The normals used for orientation should always be relative
3083 //! to the first corner of the respective surface.
3084 //! @verbatim
3085 //! Edge orientation -
3086 //!
3087 //! uMin
3088 //! --------------------> V
3089 //! | |
3090 //! | |
3091 //! | |
3092 //! vMin | | vMax
3093 //! | |
3094 //! | |
3095 //! |________________|
3096 //! | uMax
3097 //! |
3098 //! v
3099 //!
3100 //! U @endverbatim
3101
3102 cout.precision ( 16 );
3103
3104 // Calculate normal orientation.
3105
3106 snlVector refNorm = calcNormal ( minU(), minV() );
3107
3108 double orientation1 = 1.0;
3109
3110 if ( refNorm.dot ( frontFaceNormal ) < 0.0 )
3111 orientation1 = - 1.0;
3112
3113 refNorm = surface2.calcNormal ( surface2.minU(), surface2.minV() );
3114
3115 double orientation2 = 1.0;
3116
3117 if ( refNorm.dot ( frontFaceNormal2 ) < 0.0 )
3118 orientation2 = - 1.0;
3119
3120 // Generate list of starting fillet locations.
3121
3122 ptrList < arcLocn > arcLocnList;
3123
3124 knot min_u = minU();
3125 knot max_u = maxU();
3126 knot min_v = minV();
3127 knot max_v = maxV();
3128
3129 knot paramU;
3130 knot paramV;
3131
3132 knot constParam;
3133
3134 bool stepU = false;
3135
3136 switch ( edge )
3137 {
3138 case SNL_EDGE_UMIN:
3139
3140 paramU = min_u;
3141 constParam = paramU;
3142 paramV = min_v;
3143 break;
3144
3145 case SNL_EDGE_VMIN:
3146
3147 paramU = min_u;
3148 paramV = min_v;
3149 constParam = paramV;
3150 stepU = true;
3151 break;
3152
3153 case SNL_EDGE_UMAX:
3154
3155 paramU = max_u;
3156 constParam = paramU;
3157 paramV = min_v;
3158 break;
3159
3160 case SNL_EDGE_VMAX:
3161
3162 paramU = min_u;
3163 paramV = max_v;
3164 constParam = paramV;
3165 stepU = true;
3166 break;
3167 };
3168
3169 snlCurve* edgeCurve = extractEdge ( edge );
3170
3171 edgeCurve -> refine ( tolerance );
3172
3173 snlVector normal;
3174
3175 // Generate seed arc locations.
3176
3177 snlCtrlPoint* edgePoints = ( edgeCurve -> controlPointNet() ).getCtrlPtsPtr();
3178
3179 int edgeCurveSize = edgeCurve -> size();
3180
3181 snlPoint* ptsToProj = new snlPoint [ edgeCurveSize ];
3182
3183 for ( int index = 0; index < edgeCurveSize; index ++ )
3184 ptsToProj [ index ] = edgePoints [ index ];
3185
3186 int numProjLocns;
3187
3188 snlSurfLocn* initialLocns = fastProject ( ptsToProj, edgeCurveSize, &numProjLocns, tolerance, 1.0e-8, 10, 1, 1 );
3189
3190 // Make sure first point is min param.
3191
3192 int arcIndex = 0;
3193
3194 knot lastParam;
3195
3196 for ( int index = -1; index < edgeCurveSize; index ++ )
3197 {
3198 if ( index == edgeCurveSize - 1 )
3199 {
3200 // Make sure max param is obtained for last point.
3201
3202 if ( stepU )
3203 {
3204 paramU = max_u;
3205 paramV = constParam;
3206 }
3207 else
3208 {
3209 paramU = constParam;
3210 paramV = max_v;
3211 }
3212 }
3213 else if ( index > - 1 )
3214 {
3215 paramU = initialLocns [ index ].paramU;
3216 paramV = initialLocns [ index ].paramV;
3217 }
3218
3219 // Check for duplicate parameter.
3220
3221 if ( index > -1 )
3222 {
3223 if ( ( stepU && paramU == lastParam ) || ( ! stepU && paramV == lastParam ) )
3224 continue;
3225 }
3226
3227 if ( stepU )
3228 lastParam = paramU;
3229 else
3230 lastParam = paramV;
3231
3232 arcLocn* locn = new arcLocn;
3233
3234 locn -> cParamU = paramU;
3235 locn -> cParamV = paramV;
3236
3237 velocities ( paramU, paramV, locn -> cPt, locn -> velU, locn -> velV );
3238
3239 // Controlling surface
3240
3241 normal.crossProduct ( locn -> velU, locn -> velV );
3242
3243 normal *= orientation1;
3244
3245 normal.length ( radius );
3246
3247 locn -> normPt = locn -> cPt + normal;
3248
3249 locn -> converged = false;
3250 locn -> stalled = false;
3251
3252 locn -> arcIndex = arcIndex ++;
3253
3254 arcLocnList.append ( locn, true );
3255
3256 }
3257
3258 delete[] initialLocns;
3259
3260 cout << "Finished stepping out. Num Locations: " << arcLocnList.count() << "\n";
3261 // Project control points to matching surface then calculate new control points
3262
3263 double normTol = tolerance / sqrt ( tolerance * tolerance + radius * radius ); // Normal tolerance for projection.
3264
3265 int arraySize = arcLocnList.count();
3266
3267 snlPoint* points = new snlPoint [ arraySize ];
3268
3269 while ( 1 )
3270 {
3271 // Assemble points to send to the projection function.
3272
3273 arcLocn* locn = arcLocnList.first();
3274
3275 int index = 0;
3276
3277 while ( locn )
3278 {
3279 if ( ! locn -> converged && ! locn -> stalled )
3280 points [ index ++ ] = locn -> normPt;
3281
3282 locn = arcLocnList.next();
3283 }
3284
3285 // Project to matching surface.
3286
3287 int numToSend = index;
3288
3289 snlSurfLocn* surfLocns = surface2.fastProject ( points, numToSend, 0, tolerance, normTol, 10, 1, 1 );
3290
3291 // Test for convergence and create next arc location to try if needed.
3292
3293 locn = arcLocnList.first();
3294
3295 while ( locn && ( locn -> converged || locn -> stalled ) )
3296 locn = arcLocnList.next();
3297
3298 bool complete = true;
3299
3300
3301 for ( index = 0; index < numToSend; index ++ )
3302 {
3303 // If distance to matching surface is within tolerance of arc radius then
3304 // arc location has converged to an answer.
3305
3306 if ( surfLocns [ index ].dist > radius - tolerance && surfLocns [ index ].dist < radius + tolerance )
3307 {
3308 locn -> converged = true;
3309 locn -> mParamU = surfLocns [ index ].paramU;
3310 locn -> mParamV = surfLocns [ index ].paramV;
3311 locn -> mPt = surfLocns [ index ].pt;
3312 }
3313 else
3314 {
3315 // Calculate new guess for arc location on controlling surface.
3316
3317 double distDelta;
3318
3319 snlVector mRadius ( locn -> normPt, surfLocns [ index ].pt );
3320
3321 // Calculate normal to matching surface at projected point.
3322
3323 snlVector mNormal = surface2.calcNormal ( surfLocns [ index ].paramU, surfLocns [ index ].paramV );
3324
3325 mNormal *= orientation2;
3326
3327 if ( mRadius.dot ( mNormal ) > 0.0 )
3328 distDelta = radius + mRadius.length();
3329 else
3330 distDelta = mRadius.length() - radius;
3331
3332 mRadius.length ( distDelta );
3333
3334 // Get component of delta in direction we are moving.
3335
3336 if ( stepU )
3337 {
3338 // U is fixed.
3339
3340 double length = locn -> velV.length();
3341 knot deltaV = mRadius.dot ( locn -> velV ) / ( length * length );
3342
3343 knot newV = locn -> cParamV + deltaV;
3344
3345 if ( newV > max_v ) newV = max_v;
3346 if ( newV < min_v ) newV = min_v;
3347
3348 if ( newV < locn -> cParamV + SNL_NUMERIC_NOISE && newV > locn -> cParamV - SNL_NUMERIC_NOISE )
3349 locn -> stalled = true;
3350 else
3351 {
3352 locn -> cParamV = newV;
3353 complete = false;
3354 }
3355 }
3356 else
3357 {
3358 // V is fixed.
3359
3360 double length = locn -> velU.length();
3361 knot deltaU = mRadius.dot ( locn -> velU ) / ( length * length );
3362
3363 knot newU = locn -> cParamU + deltaU;
3364
3365 if ( newU > max_u ) newU = max_u;
3366 if ( newU < min_u ) newU = min_u;
3367
3368 if ( newU < locn -> cParamU + SNL_NUMERIC_NOISE && newU > locn -> cParamU - SNL_NUMERIC_NOISE )
3369 locn -> stalled = true;
3370 else
3371 {
3372 locn -> cParamU = newU;
3373 complete = false;
3374 }
3375 }
3376
3377 // Recalculate controlling surface data.
3378
3379 velocities ( locn -> cParamU, locn -> cParamV, locn -> cPt, locn -> velU, locn -> velV );
3380
3381 normal.crossProduct ( locn -> velU, locn -> velV );
3382
3383 normal *= orientation1;
3384
3385 normal.length ( radius );
3386
3387 locn -> normPt = locn -> cPt + normal;
3388 }
3389
3390 locn = arcLocnList.next();
3391
3392 while ( locn && ( locn -> converged || locn -> stalled ) )
3393 locn = arcLocnList.next();
3394 }
3395
3396 delete[] surfLocns;
3397
3398 if ( complete ) break;
3399 }
3400
3401 // Generate surface from found arcs.
3402
3403 int arcCount = 0;
3404
3405 arcLocn* locn = arcLocnList.first();
3406
3407 // Get number of arcs to create surface from and calculate
3408 // the maximum arc angle.
3409
3410 double maxAngle = 0.0;
3411
3412 while ( locn )
3413 {
3414 if ( locn -> converged )
3415 {
3416 arcCount ++;
3417
3418 snlVector startArc ( locn -> normPt, locn -> cPt );
3419 snlVector endArc ( locn -> normPt, locn -> mPt );
3420
3421 double angle = startArc.angle ( endArc );
3422
3423 if ( angle > maxAngle ) maxAngle = angle;
3424
3425 locn -> arcAngle = angle;
3426 }
3427
3428 locn = arcLocnList.next();
3429 }
3430
3431 if ( ! arcCount ) return 0; // A surface couldn't be generated.
3432
3433 // Calculate number of sections. Must be common to all arcs so that knot vector is the same across them all.
3434
3435 int numSections = (int) ( ( maxAngle / ( M_PI / 2.0 ) ) + 1 );
3436
3437 // Generate array of arcs as curves.
3438
3439 snlCurve** curves = new snlCurve* [ arcCount ];
3440
3441 locn = arcLocnList.first();
3442
3443 int index = 0;
3444
3445 while ( locn )
3446 {
3447 if ( locn -> converged )
3448 curves [ index ++ ] = new snlCurve ( locn -> cPt, locn -> mPt, locn -> normPt, numSections );
3449
3450 locn = arcLocnList.next();
3451 }
3452
3453 // Generate skinned surface from arcs.
3454
3455 snlSurface* retSurf = new snlSurface ( curves, arcCount );
3456
3457 // Clean up and return.
3458
3459 for ( index = 0; index < arcCount; index ++ )
3460 delete curves [ index ];
3461
3462 delete[] curves;
3463
3464 delete edgeCurve;
3465
3466 return retSurf;
3467 }
3468
transform(snlTransform & transf)3469 void snlSurface::transform ( snlTransform& transf )
3470 {
3471 //! Apply transformation to surface.
3472 // --------------------------------
3473
3474 ctrlPtNet -> transform ( transf );
3475 }
3476
makeCompatible(snlSurface * surfaceToMatch,int direction)3477 void snlSurface::makeCompatible ( snlSurface* surfaceToMatch, int direction )
3478 {
3479 //! Make surfaces compatible in one parametric direction.
3480 // -----------------------------------------------------
3481 //! @param surfaceToMatch Surface to make this surface compatible with.
3482 //! @param direction Parametric direction. ( enum parametricDirections ).
3483
3484 // Synchronise degree.
3485
3486 if ( direction == SNL_U_DIR )
3487 {
3488 if ( degU > ( surfaceToMatch -> degU ) )
3489 surfaceToMatch -> elevateDegree ( direction, degU - ( surfaceToMatch -> degU ) );
3490
3491 if ( degU < ( surfaceToMatch -> degU ) )
3492 elevateDegree ( direction, ( surfaceToMatch -> degU ) - degU );
3493 }
3494 else
3495 {
3496 // SNL_V_DIR.
3497
3498 if ( degV > ( surfaceToMatch -> degV ) )
3499 surfaceToMatch -> elevateDegree ( direction, degV - ( surfaceToMatch -> degV ) );
3500
3501 if ( degV < ( surfaceToMatch -> degV ) )
3502 elevateDegree ( direction, ( surfaceToMatch -> degV ) - degV );
3503 }
3504
3505 // Sync knot vectors.
3506
3507 synchronise ( *surfaceToMatch, direction );
3508 surfaceToMatch -> synchronise ( *this, direction );
3509
3510 }
3511
synchronise(snlSurface & surface,int direction)3512 void snlSurface::synchronise ( snlSurface& surface, int direction )
3513 {
3514 //! In the specified direction, synchronise this surfaces knot vector to given surfaces knot
3515 //! vector.
3516 // ----------------------------------------------------------------------------------------
3517 //! @param surface Surface to snync to.
3518 //! @param direction Parametric direction. ( enum parametricDirections ).
3519 //!
3520 //! @par Notes:
3521 //! Knots are only ever added NOT removed.
3522 //! So if surface has less multiplicity at a particular span index
3523 //! then true synchronisation will not occur and the caller
3524 //! should call the synchronise function on surface with this object
3525 //! as it's argument.
3526
3527 int tDeg, oDeg; // This degree, other surfaces degree.
3528
3529 snlKnotVector* tKnotVect; // This surfaces knot vector in direction.
3530 snlKnotVector* oKnotVect; // Other surfaces knot vector in direction.
3531
3532 if ( direction == SNL_U_DIR )
3533 {
3534 tDeg = degU;
3535 oDeg = surface.degU;
3536 tKnotVect = knotVectU;
3537 oKnotVect = surface.knotVectU;
3538 }
3539 else
3540 {
3541 tDeg = degV;
3542 oDeg = surface.degV;
3543 tKnotVect = knotVectV;
3544 oKnotVect = surface.knotVectV;
3545 }
3546
3547 if ( tDeg != oDeg ) return; // Surfaces to sync to must have same degree.
3548
3549 // Make parametric ranges equal.
3550
3551 knot thisMin = tKnotVect -> min();
3552 knot thisMax = tKnotVect -> max();
3553
3554 knot compMin = oKnotVect -> min();
3555 knot compMax = oKnotVect -> max();
3556
3557 if ( thisMin != compMin || thisMax != compMax )
3558 {
3559 // Reparameterise both curves.
3560
3561 knot newMin = thisMin > compMin ? compMin : thisMin;
3562 knot newMax = thisMax > compMax ? thisMax : compMax;
3563
3564 tKnotVect -> reparameterise ( newMin, newMax );
3565 oKnotVect -> reparameterise ( newMin, newMax );
3566 }
3567
3568 // Sync knots.
3569
3570 unsigned numSpans = oKnotVect -> getNumSpans();
3571
3572 unsigned spanIndex = oKnotVect -> getFirstSpan();
3573
3574 for ( unsigned index = 0; index < numSpans; index ++ )
3575 {
3576 knot param = oKnotVect -> val ( spanIndex );
3577
3578 int multi = oKnotVect -> findMultiplicity ( spanIndex );
3579
3580 unsigned insertSpan = tKnotVect -> findSpan ( param ); // Where param would be inserted in this object.
3581
3582 // If knot already exists in this surface then reduce multiplicity to add.
3583
3584 if ( tKnotVect -> val ( insertSpan ) == param )
3585 multi -= tKnotVect -> findMultiplicity ( insertSpan );
3586
3587 if ( multi > 0 )
3588 insertKnot ( param, direction, multi, true );
3589
3590 // Get next span.
3591
3592 spanIndex = oKnotVect -> getNextSpan ( spanIndex );
3593 }
3594 }
3595
addTrimCurve(snlCurve * curve)3596 void snlSurface::addTrimCurve ( snlCurve* curve )
3597 {
3598 //! Add trimming curve to this surface.
3599 // -----------------------------------
3600 //! @param curve Trimming curve to add. This object owns the curve.
3601
3602 trim_curves -> append ( curve, true );
3603 }
3604
removeTrimCurve(snlCurve * curve)3605 bool snlSurface::removeTrimCurve ( snlCurve* curve )
3606 {
3607 //! Delete trimming curve from this surface.
3608 // ----------------------------------------
3609 //! @param curve Trimming curve to remove.
3610
3611 if ( ! trim_curves -> hasItem ( curve ) ) return false;
3612
3613 trim_curves -> remove ( curve );
3614
3615 return true;
3616 }
3617
print()3618 void snlSurface::print()
3619 {
3620 //! Print surfaces control points and knot vectors to std out.
3621 // ----------------------------------------------------------
3622
3623 ctrlPtNet -> print();
3624
3625 cout << "Knotvector U - ";
3626 knotVectU -> print();
3627
3628 cout << "Knotvector V - ";
3629 knotVectV -> print();
3630 }
3631
print_cpp()3632 void snlSurface::print_cpp()
3633 {
3634 //! Print surfaces control points and knot vectors to std out.
3635 // ----------------------------------------------------------
3636 //! @par Notes:
3637 //! Prints data to be used for direct inclusion into a c++ program.
3638
3639 cout << "int degreeU = " << degU << ";\n";
3640 cout << "int degreeV = " << degV << ";\n\n";
3641
3642 cout << "int sizeU = " << sizeU() << ";\n";
3643 cout << "int sizeV = " << sizeV() << ";\n\n";
3644
3645 ctrlPtNet -> print_cpp();
3646
3647 cout << "\n";
3648
3649 cout << "knot knotVectorU [ " << knotVectU -> size() << " ] = ";
3650 knotVectU -> print_cpp();
3651 cout << "\n\n";
3652
3653 cout << "knot knotVectorV [ " << knotVectV -> size() << " ] = ";
3654 knotVectV -> print_cpp();
3655 cout << "\n";
3656 }
3657
vertexNet(snlVertexNet * vNet,double tolerance,bool parametric)3658 void snlSurface::vertexNet ( snlVertexNet* vNet, double tolerance, bool parametric )
3659 {
3660 //! Return approximation to surface.
3661 // --------------------------------
3662 //! @param tolerance Tolerance to surface of approximation.
3663 //! @param parametric Do a parametric analysis of the surface as opposed to knot refinement.
3664 //! @param vNet Vertex net to fill with data.
3665
3666
3667 const snlCtrlPoint* ctrlPts;
3668 int sizeU;
3669 int sizeV;
3670
3671 snlSurface* tmpSurf = 0;
3672
3673 if ( tolerance > 0.0 )
3674 {
3675 tmpSurf = new snlSurface ( *this );
3676 tmpSurf -> refine ( tolerance );
3677 ctrlPts = ( tmpSurf -> ctrlPtNet ) -> getCtrlPts();
3678 sizeU = ( tmpSurf -> ctrlPtNet ) -> getSizeU();
3679 sizeV = ( tmpSurf -> ctrlPtNet ) -> getSizeV();
3680 }
3681 else
3682 {
3683 ctrlPts = ctrlPtNet -> getCtrlPts();
3684 sizeU = ctrlPtNet -> getSizeU();
3685 sizeV = ctrlPtNet -> getSizeV();
3686 }
3687
3688 vNet -> vertexNet ( ctrlPts, sizeU, sizeV );
3689
3690 if ( tmpSurf ) delete tmpSurf;
3691 }
3692
triangleMesh(snlTriangleMesh * triMesh,int toleranceType,double tolerance)3693 void snlSurface::triangleMesh ( snlTriangleMesh* triMesh, int toleranceType, double tolerance )
3694 {
3695 //! Return approximation to surface as triangle mesh.
3696 // -------------------------------------------------
3697 //! @param triMesh Triangle mesh to populate with data.
3698 //! @param toleranceType How the tolerance is applied. See snlSurfaceBase::meshToleranceType.
3699 //! @param tolerance Mesh must be within tolerance to surface.
3700
3701 // !@#$ TEMP CODE - Just triangulate a vertexNet so that higher level stuff can be built and tested first.
3702
3703 snlVertexNet* vNet = new snlVertexNet;
3704
3705 vertexNet ( vNet, tolerance, false );
3706
3707 const snlVertex* vertexes = vNet -> vertexes();
3708
3709 int sizeU = vNet -> sizeU();
3710 int sizeV = vNet -> sizeV();
3711
3712 triMesh -> addVertexes ( vertexes, sizeU * sizeV );
3713
3714 int leftEdge, rightEdge, diagonalEdge, bottomEdge, topEdge;
3715
3716 int* nextLeftEdges = new int [ sizeV - 1 ];
3717
3718 for ( int indexU = 0; indexU < sizeU - 1; indexU ++ )
3719 {
3720 int vertIndex = indexU * sizeV;
3721
3722 for ( int indexV = 0; indexV < sizeV - 1; indexV ++ )
3723 {
3724 // Create edges and triangles.
3725
3726 if ( indexU == 0 )
3727 leftEdge = triMesh -> addEdge ( vertIndex, vertIndex + 1 ); // Left Edge.
3728 else
3729 leftEdge = nextLeftEdges [ indexV ];
3730
3731 if ( indexV == 0 )
3732 bottomEdge = triMesh -> addEdge ( vertIndex, vertIndex + sizeV ); // Bottom Edge.
3733
3734 diagonalEdge = triMesh -> addEdge ( vertIndex + 1, vertIndex + sizeV ); // Diagonal.
3735
3736 topEdge = triMesh -> addEdge ( vertIndex + 1, vertIndex + sizeV + 1 ); // Top Edge.
3737
3738 rightEdge = triMesh -> addEdge ( vertIndex + sizeV, vertIndex + sizeV + 1 ); // Right Edge.
3739
3740 nextLeftEdges [ indexV ] = rightEdge;
3741
3742 // Create two triangles per quad.
3743
3744 triMesh -> addTriangle ( leftEdge, bottomEdge, diagonalEdge );
3745 triMesh -> addTriangle ( rightEdge, topEdge, diagonalEdge );
3746
3747 // Set edge for next quad.
3748
3749 bottomEdge = topEdge;
3750 }
3751 }
3752
3753 // Clean up.
3754
3755 delete[] nextLeftEdges;
3756 }
3757
genSurfRevolution(snlCurve & generator,snlPoint & axisStart,snlPoint & axisEnd,double angle)3758 void snlSurface::genSurfRevolution ( snlCurve& generator, snlPoint& axisStart, snlPoint& axisEnd, double angle )
3759 {
3760 //! Construct surface of revolution.
3761 // --------------------------------
3762 //! @param generator Generating curve.
3763 //! @param axisStart Starting point of axis generator is revolved about.
3764 //! @param axisEnd Ending point of axis.
3765 //! @param angle Angle in degrees to revolve generator about axis. Angle is in degrees so that
3766 //! different precisions of PI do not affect surface closure.
3767 //!
3768 //! @par Notes:
3769 //! Rotation is counter clockwise about axis vector. Right hand rule. Curve defines V
3770 //! direction.
3771
3772 // Clamp angle to 360 degrees.
3773
3774 if ( angle > 360.0 ) angle = 360.0;
3775
3776 double radAngle = ( angle / 180.0 ) * M_PI;
3777
3778 // Calculate number of sections and section angle.
3779
3780 int numSections = (int) ( ( angle / 90.0 ) + 1 );
3781
3782 double sectionAngle = radAngle / (double) numSections;
3783
3784 double stepAngle = sectionAngle / 2.0;
3785
3786 // Calculate mid point weight.
3787
3788 double midPtWeight = cos ( stepAngle );
3789
3790 // Setup rotation transforms.
3791
3792 int numRotations = numSections * 2;
3793
3794 snlTransform* rotations = new snlTransform [ numRotations ];
3795
3796 for ( int index = 0; index < numRotations; index ++ )
3797 rotations [ index ].rotate ( stepAngle * (double) ( index + 1 ), axisStart, axisEnd );
3798
3799 // Generate non rotated mid point control points.
3800
3801 int sizeV = generator.size();
3802
3803 const snlCtrlPoint* curvePts = ( generator.controlPointNet() ).getCtrlPts();
3804
3805 snlCtrlPoint* midPoints = new snlCtrlPoint [ sizeV ];
3806
3807 for ( int index = 0; index < sizeV; index ++ )
3808 {
3809 // Get vector from point to axis.
3810
3811 snlVector projection = projectToLine ( axisStart, axisEnd, curvePts [ index ] );
3812
3813 projection *= - 1.0; // Vector is pointing into the axis, we want it pointing out from it.
3814
3815 double projDist = projection.length();
3816
3817 if ( projDist != 0.0 )
3818 {
3819 // Calculate mid point to axis distance.
3820
3821 double dist = projDist / midPtWeight; // Mid point weight is the cosine of the step angle.
3822
3823 // Generate new control point.
3824
3825 projection.length ( dist - projDist );
3826 }
3827
3828 midPoints [ index ] = curvePts [ index ] + projection;
3829
3830 // Multiply by mid point weight.
3831
3832 midPoints [ index ].multiplyWeight ( midPtWeight );
3833 }
3834
3835 // Generate surface control points.
3836
3837 int sizeU = ( numSections * 2 ) + 1;
3838
3839 int totalSize = sizeU * sizeV;
3840
3841 snlCtrlPoint* surfPts = new snlCtrlPoint [ totalSize ];
3842
3843 int index = 0;
3844 int stepIndex = -1;
3845
3846 for ( int indexU = 0; indexU < sizeU; indexU ++ )
3847 {
3848 for ( int indexV = 0; indexV < sizeV; indexV ++ )
3849 {
3850 if ( ! indexU || ( indexU == ( sizeU - 1 ) && angle == 360.0 ) )
3851 // If this is the first U index then curve points are taken as is.
3852 surfPts [ index ] = curvePts [ indexV ];
3853
3854 else
3855 {
3856 // Calculate rotated control point.
3857
3858 if ( stepIndex % 2 )
3859 {
3860 // Not a mid point.
3861 surfPts [ index ] = curvePts [ indexV ];
3862 }
3863 else
3864 {
3865 // Is a mid point.
3866 surfPts [ index ] = midPoints [ indexV ];
3867 }
3868
3869 rotations [ stepIndex ].transform ( surfPts [ index ] );
3870 }
3871
3872 index ++;
3873 }
3874
3875 stepIndex ++;
3876 }
3877
3878 ctrlPtNet = new snlCtrlPointNetSurface ( surfPts, sizeU, sizeV );
3879
3880 // Generate Knot Vectors.
3881
3882 // U Knot Vector.
3883
3884 knot* uKnots = new knot [ sizeU + 3 ]; // Degree 2 knot vector.
3885
3886 for ( index = 0; index < 3; index ++ )
3887 {
3888 // End clamps.
3889
3890 uKnots [ index ] = 0.0;
3891 uKnots [ sizeU + index ] = 1.0;
3892 }
3893
3894 // Internal knots.
3895
3896 index = 3;
3897
3898 knot knotStep = 1.0 / (double) ( numSections );
3899 knot knotVal = knotStep;
3900
3901 for ( int step = 0; step < numSections - 1; step ++ )
3902 {
3903 uKnots [ index ++ ] = knotVal; // Multiplicity 2.
3904 uKnots [ index ++ ] = knotVal;
3905
3906 knotVal += knotStep;
3907 }
3908
3909 knotVectU = new snlKnotVector ( uKnots, sizeU + 3, 2 );
3910
3911 // V Knot Vector.
3912
3913 knotVectV = new snlKnotVector ( generator.knotVector() );
3914
3915 // Setup remaining variables.
3916
3917 degU = 2;
3918 degV = generator.degree();
3919
3920 // Clean up.
3921
3922 delete[] rotations;
3923 delete[] midPoints;
3924 }
3925
3926