1 // libSNL - Simple Nurbs Library
2 // Copyright 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 Library General Public License for more details.
13 //
14 // You should have received a copy of the GNU Library General Public License
15 // along with this program; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.
17
18 // *** NURBS routines that don't belong in classes ***
19
20 // *** Legacy OpenCDSM routines are included ***
21
22 // *** NOTE: OpenCDSM functions still use t,u for surface parameters ***
23
24 #include "snlNurbsCommon.h"
25 #include "snlVector.h"
26
27 //#define PROJ_COMMENT
28
29 static int maxIterations = 0;
30 static unsigned maxPasses = 0;
31
newtonIterStepSurf(snlPoint * derivPts,snlPoint * projPt,knot * deltaU,knot * deltaV)32 bool newtonIterStepSurf ( snlPoint* derivPts, snlPoint* projPt, knot* deltaU, knot* deltaV )
33 {
34 // Evaluate a single step in a newton iteration over a surface
35 // -----------------------------------------------------------
36 // derivPts: Pre calculated rational 0th to 2nd derivatives at desired paramters.
37 // This functions expectes a 3 X 3 array.
38 // projPt: Point being projected onto surface.
39 // deltaU: Change in u parameter to return.
40 // deltaV: Change in v parameter to return.
41 //
42 // returns: True if processing was succesfull. False if deltas would be infinite.
43 //
44 // Notes: Theory from "The NURBS Book 2nd Ed" pages 232 - 234.
45
46 basis jcbn [ 2 ] [ 2 ]; // 2 X 2 Jacobian matrix.
47
48 snlVector dist ( *projPt, derivPts [ 0 ] ); // Distance from point to be projected to surface.
49
50 jcbn [ 0 ] [ 0 ] = derivPts [ 3 ].lengthSqrd() + dist.dot ( derivPts [ 6 ] );
51 jcbn [ 0 ] [ 1 ] = snlVector ( derivPts [ 3 ] ).dot ( derivPts [ 1 ] )
52 + dist.dot ( derivPts [ 4 ] );
53 jcbn [ 1 ] [ 0 ] = jcbn [ 0 ] [ 1 ];
54 jcbn [ 1 ] [ 1 ] = derivPts [ 1 ].lengthSqrd( ) + dist.dot ( derivPts [ 2 ] );
55
56 return solve2X2LinEqn ( jcbn [ 0 ] [ 0 ], jcbn [ 1 ] [ 0 ], jcbn [ 0 ] [ 1 ], jcbn [ 1 ] [ 1 ],
57 - ( dist.dot ( derivPts [ 3 ] ) ),
58 - ( dist.dot ( derivPts [ 1 ] ) ),
59 deltaU, deltaV );
60 }
61
newtonIterStepCurve(snlPoint * derivPts,snlPoint * projPt,knot * paramDelta)62 bool newtonIterStepCurve ( snlPoint* derivPts, snlPoint* projPt, knot* paramDelta )
63 {
64 // Evaluate a single step in a newton iteration over a curve
65 // ---------------------------------------------------------
66 // derivPts: Pre calculated rational 0th to 2nd derivatives at desired paramter. Size 3.
67 // projPt: Point being projected onto curve.
68 // paramDelta: Change in t parameter to return.
69 //
70 // returns: True if processing was succesfull. False if paramDelta would be infinite.
71 //
72 // Notes: Theory from "The NURBS Book 2nd Ed" pages 230 - 231.
73
74 snlVector dist ( *projPt, derivPts [ 0 ] ); // Distance from point to be projected to surface.
75
76 if ( ( dist.dot ( derivPts [ 2 ] ) + derivPts [ 1 ].lengthSqrd() ) == 0.0 ) return false;
77
78 *paramDelta = - ( dist.dot ( derivPts [ 1 ] ) /
79 ( dist.dot ( derivPts [ 2 ] ) + derivPts [ 1 ].lengthSqrd() ) );
80
81 return true;
82 }
83
lineIterStepCurve(snlPoint * derivPts,snlPoint * projPt,knot * paramDelta)84 bool lineIterStepCurve ( snlPoint* derivPts, snlPoint* projPt, knot* paramDelta )
85 {
86 // Evaluate a single step in a line iteration over a degree 1 curve
87 // ----------------------------------------------------------------
88 // derivPts: Pre calculated rational 0th to 2nd derivatives at desired paramter. Size 3.
89 // projPt: Point being projected onto curve.
90 // paramDelta: Change in t parameter to return.
91 //
92 // returns: True if processing was succesfull. False if paramDelta would be infinite.
93
94 basis length;
95
96 // Check for divide by zero error.
97 length = derivPts [ 1 ].lengthSqrd();
98
99 if ( length == 0.0 ) return false;
100
101 snlVector dist ( *projPt, derivPts [ 0 ] ); // Distance from point to be projected to surface.
102
103 // Vectors are nose to tail so have to negate. Look at definition of dot product.
104 *paramDelta = -( dist.dot ( derivPts [ 1 ] ) ) / length;
105
106 return true;
107 }
108
solve2X2LinEqn(double a1,double a2,double a3,double a4,double b1,double b2,double * x1,double * x2)109 int solve2X2LinEqn ( double a1, double a2, double a3, double a4,
110 double b1, double b2, double* x1, double* x2 )
111 {
112 /* Solve a 2 X 2 system of linear equations
113 // ----------------------------------------
114 //
115 // a1 to a4: Coefficients of the system.
116 // b1, b2: If the system is represented as Ax = b, these are the b's.
117 //
118 // _ _ _ _
119 // | | | |
120 // | a1.x1 a3.x2 | | b1 |
121 // | | = | |
122 // | a2.x1 a4.x2 | | b2 |
123 // |_ _| |_ _|
124 //
125 // x1, x2: Pointers to return solutions in.
126 //
127 // Returns: False if det is 0, otherwise True. ie Returned values would be infinite.
128 */
129
130 double det, det1, det2;
131
132 det = a1 * a4 - a2 * a3;
133
134 det1 = b1 * a4 - b2 * a3;
135
136 det2 = a1 * b2 - a2 * b1;
137
138 if ( det != 0 )
139 {
140 *x1 = det1 / det;
141 *x2 = det2 / det;
142
143 return 1;
144 }
145
146 return 0;
147 }
148
projPtSurf(snlSurface & surface,snlPoint * toProj,int toProjNum,sLocn * best,double iterTol,double cosTol,unsigned maxPass)149 ptrList < sLocn >* projPtSurf ( snlSurface& surface, snlPoint* toProj, int toProjNum, sLocn* best,
150 double iterTol, double cosTol, unsigned maxPass )
151 {
152 // Project points onto given surface
153 // ---------------------------------
154 // surface: Surface to project to.
155 // toProj: Point to project.
156 // toProjNum: Number of points to project
157 // best: Array to return best evaluted locations in. Must be toProjNum in size.
158 // iterTol: Iteration tolerance. Maximum distance between surface points after
159 // successive iterations. If below this value then processing stops.
160 // cosTol: Cosine Tolerance. Maxium allowable deviation of cosine from zero.
161 // i.e. How close to a surface normal it is.
162 // maxPass: Maximum number of passes to be made. Each subsequent pass doubles the number of
163 // sections each span is divided into. For ordinary cases this only needs to be 1.
164 //
165 //
166 // Notes: Evaluates more than one point at once.
167 // Mostly OpenCDSM legacy code and uses OpenCDSM parameter names.
168 //
169 // This is messy inefficient code! But it mostly works ;-)
170
171 knot paramT, paramU; // Evaluation parameters.
172 knot pTStart, pTEnd; // Param T start and end for span.
173 knot pUStart, pUEnd; // Param U start and end for span.
174 knot pTStep, pUStep; // Param T and U incremental steps.
175 knot pTStepDenom, pUStepDenom; // PxStep denominators.
176
177 snlPoint evalPt; // Current evaled non-homogeneous point.
178
179 snlPoint* evalDerivs = 0; // Derivatives of point needed for newton iterations.
180 snlPoint* newEvalDerivs = 0;
181
182 bool withinTol = false; // Found a point within tolerance.
183
184 bool noBest = true; // No best results have been initialised yet.
185 int index;
186
187 unsigned cPass;
188
189 bool* withinTols = new bool [ toProjNum ];
190
191 ptrList < projLocn >* bestLists = new ptrList < projLocn > [ toProjNum ];
192 ptrList < projLocn >* cList;
193
194 projLocn* cBest;
195
196 basis bestDist; // Used for point selection with multiple hits.
197
198 ptrList < sLocn >* ambig = new ptrList < sLocn >;
199
200 unsigned degT = surface.degreeU();
201 unsigned degU = surface.degreeV();
202
203 bool degT1 = ( degT == 1 ); // Degree T is 1.
204 bool degU1 = ( degU == 1 ); // Degree U is 1.
205
206 for ( index = 0; index < toProjNum; index ++ ) withinTols [ index ] = false;
207
208 #ifdef PROJ_COMMENT
209
210 cout << "\nNew Projection\n";
211 cout.precision ( 15 );
212
213 #endif
214
215 const snlKnotVector& kntsT = surface.knotVectorU();
216 const snlKnotVector& kntsU = surface.knotVectorV();
217
218 // Calculate bisection step multipliers used to account for surfaces with both
219 // high curvature and high discontinuity ( small number of non-zero spans, large number of knots ).
220
221 int stepMultiT = (int) ( surface.maxCurvatureU() * 2.0 ) + ( surface.sizeU() / kntsT.getNumSpans() );
222 int stepMultiU = (int) ( surface.maxCurvatureV() * 2.0 ) + ( surface.sizeV() / kntsU.getNumSpans() );
223
224 #ifdef PROJ_COMMENT
225 cout << "stepMultiT: " << stepMultiT << " stepMultiU: " << stepMultiU << "\n";
226 #endif
227
228 // Evaluate surface points that correspond to beginning and end of knot spans
229 // and store point with least distance from surface points.
230
231 for ( cPass = 1; cPass <= maxPass; cPass ++ )
232 {
233
234 #ifdef PROJ_COMMENT
235
236 cout << "Pass: " << cPass << "\n";
237
238 #endif
239
240 unsigned minSpanT = kntsT.findSpan ( kntsT.min() );
241 unsigned maxSpanT = kntsT.findSpan ( kntsT.max() );
242
243 unsigned minSpanU = kntsU.findSpan ( kntsU.min() );
244 unsigned maxSpanU = kntsU.findSpan ( kntsU.max() );
245
246 // Seed best point.
247 paramT = kntsT.val ( minSpanT );
248 paramU = kntsU.val ( minSpanU );
249
250 evalPt = surface.eval ( paramT, paramU );
251
252 // On the first loop the value found is the best.
253 if ( noBest )
254 {
255 // Initialise pointer lists
256 for ( index = 0; index < toProjNum; index ++ )
257 {
258 cBest = new projLocn;
259
260 cBest -> paramT = paramT;
261 cBest -> paramU = paramU;
262 cBest -> pt = evalPt;
263
264 bestLists [ index ].append ( cBest, true );
265 }
266
267 noBest = false;
268 }
269
270 for ( unsigned spanT = minSpanT; spanT <= maxSpanT; spanT ++ )
271 {
272 pTStart = kntsT.val ( spanT );
273 pTEnd = kntsT.val ( spanT + 1 );
274
275 // Only work with non-zero length knot spans.
276 if ( pTStart < pTEnd )
277 {
278 // Set T bisection step value.
279 pTStepDenom = ( ( degT + 1 ) * cPass * stepMultiT );
280 pTStep = ( pTEnd - pTStart ) / pTStepDenom;
281
282 for ( unsigned spanU = minSpanU; spanU <= maxSpanU; spanU ++ )
283 {
284 pUStart = kntsU.val ( spanU );
285 pUEnd = kntsU.val ( spanU + 1 );
286
287 // Only work with non-zero length knot spans.
288 if ( pUStart < pUEnd )
289 {
290 // Set U bisection step value.
291 pUStepDenom = ( ( degU + 1 ) * cPass * stepMultiU );
292 pUStep = ( pUEnd - pUStart ) / pUStepDenom;
293
294 // Bisect span "rectangle", spanT -> spanT + 1, spanU -> spanU + 1.
295 // Step through each section. Some points are checked twice for the sake
296 // of code simplicity.
297
298 for ( unsigned sectT = 0; sectT <= (unsigned) pTStepDenom; sectT ++ )
299 {
300 paramT = pTStart + pTStep * sectT;
301
302 for ( unsigned sectU = 0; sectU <= (unsigned) pUStepDenom; sectU ++ )
303 {
304 paramU = pUStart + pUStep * sectU;
305
306 // Evaluate
307 evalPt = surface.eval ( paramT, paramU );
308
309 // Iterate through each point to be projected.
310 for ( index = 0; index < toProjNum; index ++ )
311 {
312 if ( withinTols [ index ] ) continue;
313
314 cList = bestLists + index;
315
316 cBest = cList -> first();
317
318 // See if evalStart [ 0 ] is the best fit so far.
319 if ( snlVector ( toProj [ index ], cBest -> pt ).length()
320 > snlVector ( toProj [ index ], evalPt ).length() )
321 {
322 // If the evaluated point is closer to the point to be projected
323 // then it is the new best.
324 cBest -> paramT = paramT;
325 cBest -> paramU = paramU;
326 cBest -> pt = evalPt;
327
328 cList -> truncate();
329 }
330 else if ( ( evalPt == ( cBest -> pt ) ) &&
331 ( cBest -> paramT != paramT ||
332 cBest -> paramU != paramU ) )
333 {
334 // Coincident point found.
335
336 // Make sure point hasn't been found before.
337
338 bool exists = false;
339
340 while ( cBest )
341 {
342 if ( paramT == cBest -> paramT && paramU == cBest -> paramU )
343 exists = true;
344
345 cBest = cList -> next();
346 }
347
348 if ( ! exists )
349 {
350 cBest = new projLocn;
351 cBest -> paramT = paramT;
352 cBest -> paramU = paramU;
353 cBest -> pt = evalPt;
354
355 cList -> append ( cBest, true );
356 }
357 }
358 }
359 }
360 }
361 }
362 }
363 }
364 }
365
366 // Whoa! that is far too many nested for loops *grin*.
367
368 // Okey dokey. Now it is time to try a newton iteration or two ( thousand ) on the point of best fit.
369
370 for ( index = 0; index < toProjNum; index ++ )
371 {
372 if ( withinTols [ index] ) continue;
373
374 cBest = bestLists [ index ].first();
375
376 double guessDistance = snlVector ( cBest -> pt, toProj [ index ] ).length();
377
378 #ifdef PROJ_COMMENT
379 cout << "Processing Index: " << index << "\n";
380 cout << "Number in bestList: " << bestLists [ index ].count() << "\n";
381 cout << "Original Point: "; toProj [ index ].print(); cout << "\n";
382 #endif
383
384 while ( cBest )
385 {
386 #ifdef PROJ_COMMENT
387 cout << "Guess: ";cBest -> pt.print(); cout << "\n";
388 #endif
389
390 for ( int count = 0; count < MAX_PROJ_ITER; count ++ )
391 {
392 withinTol = false;
393
394 if ( maxIterations < count ) maxIterations = count;
395
396 #ifdef PROJ_COMMENT
397
398 cout << "Iteration Number: " << count << "\n";
399
400 #endif
401
402 // If proj point is within iteration tolerance from evaluated point then stop.
403
404 cBest -> dist = snlVector ( cBest -> pt, toProj [ index ] ).length();
405
406 if ( ( cBest -> dist ) < iterTol )
407 {
408 withinTol = true;
409 break;
410 }
411
412 // Get 1st and 2nd derivatives of point.
413 if ( ! evalDerivs )
414 evalDerivs = surface.evalDerivs ( cBest -> paramT, cBest -> paramU, 2, 2 );
415
416 // Calc projPt to surface vector.
417 snlVector projToSurf ( toProj [ index ], evalDerivs [ 0 ] );
418
419 // Check to see if cosine of determination function is under tolerance ( cosTol ).
420
421 snlVector tVect ( evalDerivs [ 3 ] );
422 basis cosT = projToSurf.calcAbsCos ( tVect );
423
424 snlVector uVect ( evalDerivs [ 1 ] );
425 basis cosU = projToSurf.calcAbsCos ( uVect );
426
427 #ifdef PROJ_COMMENT
428
429 cout << "Iteration: " << count << " DotT: " << projToSurf.dot ( tVect ) << " DotU: " << projToSurf.dot ( uVect ) << "\n";
430 cout << "Cosine T: " << cosT << "\n";
431 cout << "Cosine U: " << cosU << "\n";
432
433 #endif
434
435 if ( cosU <= cosTol && cosT <= cosTol )
436 {
437 withinTol = true;
438 break;
439 }
440
441 knot deltaT;
442 knot deltaU;
443
444 knot newT;
445 knot newU;
446
447 // Get new paramaters.
448 if ( !degT1 && !degU1 )
449 {
450 if ( ! newtonIterStepSurf ( evalDerivs, toProj + index, &deltaT, &deltaU ) )
451 {
452 break; // Param deltas would have gone to infinity.
453 }
454 }
455 else
456 {
457 // At least one of the degrees is 1.
458
459 snlPoint tDerivs [ 3 ];
460
461 tDerivs [ 0 ] = evalDerivs [ 0 ];
462 tDerivs [ 1 ] = evalDerivs [ 3 ];
463 tDerivs [ 2 ] = evalDerivs [ 6 ];
464
465 if ( degT1 )
466 {
467 if ( ! lineIterStepCurve ( tDerivs, toProj + index, &deltaT ) )
468 break;
469
470 }
471 else
472 {
473 if ( ! newtonIterStepCurve ( tDerivs, toProj + index, &deltaT ) )
474 break;
475 }
476
477 if ( degU1 )
478 {
479 if ( ! lineIterStepCurve ( evalDerivs, toProj + index, &deltaU ) )
480 break;
481 }
482 else
483 {
484 if ( ! newtonIterStepCurve ( evalDerivs, toProj + index, &deltaU ) )
485 break;
486 }
487 }
488
489 // Clamp params to knot bounds.
490 newT = cBest -> paramT + deltaT;
491 newU = cBest -> paramU + deltaU;
492
493 if ( newT > ( kntsT.max() ) )
494 {
495 newT = kntsT.max();
496 deltaT = newT - cBest -> paramT;
497 }
498
499 if ( newT < ( kntsT.min() ) )
500 {
501 newT = kntsT.min();
502 deltaT = newT - cBest -> paramT;
503 }
504
505 if ( newU > ( kntsU.max() ) )
506 {
507 newU = kntsU.max();
508 deltaU = newU - cBest -> paramU;
509 }
510
511 if ( newU < ( kntsU.min() ) )
512 {
513 newU = kntsU.min();
514 deltaU = newU - cBest -> paramU;
515 }
516
517 #ifdef PROJ_COMMENT
518 cout << "OldT: " << cBest -> paramT << "\n";
519 cout << "OldU: " << cBest -> paramU << "\n";
520 cout << "NewT: " << newT << "\n";
521 cout << "NewU: " << newU << "\n";
522 cout << "deltaT: " << deltaT << "\n";
523 cout << "deltaU: " << deltaT << "\n";
524 #endif
525
526 // If parameters haven't changed then break.
527 if ( deltaT == 0.0 && deltaU == 0.0 )
528 break;
529
530 // Evaluate new parameters.
531
532 newEvalDerivs = surface.evalDerivs ( newT, newU, 2, 2 );
533
534 // Generate new best.
535 cBest -> pt = newEvalDerivs [ 0 ];
536 cBest -> paramT = newT;
537 cBest -> paramU = newU;
538
539 // Calculate new distance to surface.
540 cBest -> dist = snlVector ( cBest -> pt, toProj [ index ] ).length();
541
542 // Check for below tolerance changes in surface point.
543 snlVector chgT ( evalDerivs [ 3 ] );
544 chgT *= deltaT;
545
546 snlVector chgU ( evalDerivs [ 1 ] );
547 chgU *= deltaU;
548
549 // Change vector.
550 snlVector chgV = chgT + chgU;
551
552 #ifdef PROJ_COMMENT
553 cout << "Change in T direction: " << chgT.length() << "\n";
554 cout << "Change in U direction: " << chgT.length() << "\n";
555 cout << "Change in point: " << chgV.length() << "\n";
556 #endif
557
558 if ( chgV.length() <= iterTol )
559 {
560 withinTol = true;
561 break;
562 }
563
564 // Reorganise evaluation arrays.
565 if ( evalDerivs ) delete[] evalDerivs;
566 evalDerivs = newEvalDerivs;
567 newEvalDerivs = 0;
568 }
569
570 // Detect whether the point to be projected is within tolerances
571 // Projected points that are further from surface than original
572 // guess are not within tolerance.
573
574 if ( withinTol && ( cBest -> dist <= guessDistance ) )
575 {
576 withinTols [ index ] = withinTol;
577 cBest -> flag = true;
578 }
579 else
580 cBest -> flag = false;
581
582 cBest = bestLists [ index ].next();
583 }
584
585 if ( withinTols [ index ] )
586 {
587 // ** Process best points **
588
589 // Get best distance.
590
591 cBest = bestLists [ index ].first();
592
593 bestDist = cBest -> dist;
594
595 cBest = bestLists [ index ].next();
596
597 while ( cBest )
598 {
599 if ( bestDist > ( cBest -> dist ) ) bestDist = ( cBest -> dist );
600
601 cBest = bestLists [ index ].next();
602 }
603
604 if ( bestDist < iterTol ) bestDist = iterTol; // Accounts for discrepencies in Newton convergence.
605
606 cBest = bestLists [ index ].first();
607
608 int numBest = 0;
609
610 while ( cBest )
611 {
612 // Get the point closest to the surface.
613
614 if ( ( cBest -> flag ) && ( ( cBest -> dist ) <= bestDist ) )
615 {
616 if ( ! numBest )
617 {
618 best [ index ].paramT = cBest -> paramT;
619 best [ index ].paramU = cBest -> paramU;
620 best [ index ].pt = cBest -> pt;
621 }
622 else
623 {
624 // Add item to ambiguity list.
625 sLocn* ambLocn = new sLocn;
626 ambLocn -> paramT = cBest -> paramT;
627 ambLocn -> paramU = cBest -> paramU;
628 ambLocn -> pt = cBest -> pt;
629 ambLocn -> flag = index;
630
631 ambig -> append ( ambLocn, true );
632 }
633
634 numBest ++;
635 }
636 else
637 cBest -> flag = false;
638
639 cBest = bestLists [ index ].next();
640 }
641
642 best [ index ].flag = numBest;
643 }
644 else
645 {
646 // Make sure at least something is in the best array.
647 cBest = bestLists [ index ].first();
648
649 if ( cBest )
650 {
651 best [ index ].paramT = cBest -> paramT;
652 best [ index ].paramU = cBest -> paramU;
653 best [ index ].pt = cBest -> pt;
654 best [ index ].flag = 1;
655 }
656 }
657
658 // Clean up.
659 if ( evalDerivs )
660 {
661 delete[] evalDerivs;
662 evalDerivs = 0;
663 }
664
665 if ( newEvalDerivs )
666 {
667 delete[] newEvalDerivs;
668 newEvalDerivs = 0;
669 }
670
671 #ifdef PROJ_COMMENT
672 cout << "Best Found: "; best [ index ].pt.print(); cout << "\n";
673 #endif
674 }
675
676 // All point tolerances must be within threshold to break out of pass loop.
677
678 withinTol = true;
679
680 for ( index = 0; index < toProjNum; index ++ )
681 {
682 if ( ! withinTols [ index ] ) withinTol = false;
683 }
684
685 if ( withinTol ) break;
686 }
687
688 delete[] withinTols;
689
690 if ( maxPasses < cPass ) maxPasses = cPass;
691
692 #ifdef PROJ_COMMENT
693
694 cout << "Max Passes: " << maxPasses << "\n";
695 cout << "Max Iterations: " << maxIterations << "\n";
696
697 #endif
698
699 if ( ambig -> count() )
700 resolveAmbig ( best, toProjNum, ambig );
701
702 return ambig;
703 }
704
resolveAmbig(sLocn * projns,int projSize,ptrList<sLocn> * ambig)705 void resolveAmbig ( sLocn* projns, int projSize, ptrList < sLocn >* ambig )
706 {
707 // Attempt to resolve projection ambiguities.
708 // ------------------------------------------
709 // projns: Projections.
710 // projSize: Size of projns array.
711 // ambig: Ambiguities.
712
713 // Generate array so that ambiguous points before and after
714 // current point being considered can be easily checked.
715
716 sLocn* locn;
717
718 if ( projSize < 3 ) return; // Can't get a trend with less than three points.
719
720 bool* ambigIndexes = new bool [ projSize ];
721
722 for ( int index = 0; index < projSize; index ++ )
723 ambigIndexes [ index ] = false;
724
725 for ( locn = ambig -> first(); locn; locn = ambig -> next() )
726 ambigIndexes [ locn -> flag ] = true;
727
728 bool reverse = false;
729
730 // Process projection ambiguities.
731
732 // If parametric distance of ambiguous point is closer to
733 // neighbours then this point is considered a better match.
734
735 // Process forwards.
736
737 for ( locn = ambig -> first(); locn; locn = ambig -> next() )
738 {
739 // Don't compare to ambiguous neighbours.
740
741 if ( ! ( locn -> flag ) )
742 {
743 reverse = true;
744 continue;
745 }
746
747 if ( ambigIndexes [ ( locn -> flag ) - 1 ] )
748 {
749 reverse = true;
750 continue;
751 }
752
753 ambigIndexes [ ( locn -> flag ) ] = false; // Don't process in reverse.
754
755 knot cDist = paramDistSqrd ( projns [ ( locn -> flag ) - 1 ].paramT,
756 projns [ ( locn -> flag ) - 1 ].paramU,
757 projns [ ( locn -> flag ) ].paramT,
758 projns [ ( locn -> flag ) ].paramU );
759
760 knot tDist = paramDistSqrd ( projns [ ( locn -> flag ) - 1 ].paramT,
761 projns [ ( locn -> flag ) - 1 ].paramU,
762 locn -> paramT,
763 locn -> paramU );
764
765 // If the sLocn being tested is a better match then swap data.
766
767 if ( tDist < cDist )
768 {
769 sLocn tmpLocn = *locn;
770 *locn = projns [ tmpLocn.flag ];
771 projns [ tmpLocn.flag ] = tmpLocn;
772
773 // Make sure flags aren't swapped.
774 projns [ tmpLocn.flag ].flag = locn -> flag;
775 locn -> flag = tmpLocn.flag;
776 }
777 }
778
779 // Process Backwards. Try and pick up points missed during forward processing.
780
781 if ( reverse )
782 {
783 for ( locn = ambig -> last(); locn; locn = ambig -> prev() )
784 {
785 if ( ! ambigIndexes [ ( locn -> flag ) ] )
786 continue; // Has already been processed.
787
788 // Don't compare to ambiguous neighbours.
789
790 if ( ( locn -> flag ) >= projSize - 1 ) // Is last point.
791 continue;
792
793 if ( ambigIndexes [ ( locn -> flag ) + 1 ] )
794 continue;
795
796 knot cDist = paramDistSqrd ( projns [ ( locn -> flag ) + 1 ].paramT,
797 projns [ ( locn -> flag ) + 1 ].paramU,
798 projns [ ( locn -> flag ) ].paramT,
799 projns [ ( locn -> flag ) ].paramU );
800
801 knot tDist = paramDistSqrd ( projns [ ( locn -> flag ) + 1 ].paramT,
802 projns [ ( locn -> flag ) + 1 ].paramU,
803 locn -> paramT,
804 locn -> paramU );
805
806 // If the sLocn being tested is a better match then swap data.
807
808 if ( tDist < cDist )
809 {
810 sLocn tmpLocn = *locn;
811 *locn = projns [ tmpLocn.flag ];
812 projns [ tmpLocn.flag ] = tmpLocn;
813
814 // Make sure flags aren't swapped.
815 projns [ tmpLocn.flag ].flag = locn -> flag;
816 locn -> flag = tmpLocn.flag;
817 }
818 }
819
820 }
821
822 delete[] ambigIndexes;
823 }
824
paramDistSqrd(knot t1,knot u1,knot t2,knot u2)825 knot paramDistSqrd ( knot t1, knot u1, knot t2, knot u2 )
826 {
827 return ( t2 - t1 ) * ( t2 - t1 ) + ( u2 - u1 ) * ( u2 - u1 );
828 }
829
830