Point Cloud Library (PCL) 1.13.0
Loading...
Searching...
No Matches
opennurbs_math.h
1/* $NoKeywords: $ */
2/*
3//
4// Copyright (c) 1993-2012 Robert McNeel & Associates. All rights reserved.
5// OpenNURBS, Rhinoceros, and Rhino3D are registered trademarks of Robert
6// McNeel & Associates.
7//
8// THIS SOFTWARE IS PROVIDED "AS IS" WITHOUT EXPRESS OR IMPLIED WARRANTY.
9// ALL IMPLIED WARRANTIES OF FITNESS FOR ANY PARTICULAR PURPOSE AND OF
10// MERCHANTABILITY ARE HEREBY DISCLAIMED.
11//
12// For complete openNURBS copyright information see <http://www.opennurbs.org>.
13//
14////////////////////////////////////////////////////////////////
15*/
16
17#if !defined(ON_MATH_INC_)
18#define ON_MATH_INC_
19
20class ON_3dVector;
21class ON_Interval;
22class ON_Line;
23class ON_Arc;
24class ON_Plane;
25
26/*
27Description:
28 Class for carefully adding long list of numbers.
29*/
30class ON_CLASS ON_Sum
31{
32public:
33
34 /*
35 Description:
36 Calls ON_Sum::Begin(x)
37 */
38 void operator=(double x);
39
40 /*
41 Description:
42 Calls ON_Sum::Plus(x);
43 */
44 void operator+=(double x);
45
46 /*
47 Description:
48 Calls ON_Sum::Plus(-x);
49 */
50 void operator-=(double x);
51
52 /*
53 Description:
54 Creates a sum that is ready to be used.
55 */
57
58 /*
59 Description:
60 If a sum is being used more than once, call Begin()
61 before starting each sum.
62 Parameters:
63 starting_value - [in] Initial value of sum.
64 */
65 void Begin( double starting_value = 0.0 );
66
67 /*
68 Description:
69 Add x to the current sum.
70 Parameters:
71 x - [in] value to add to the current sum.
72 */
73 void Plus( double x );
74
75 /*
76 Description:
77 Calculates the total sum.
78 Parameters:
79 error_estimate - [out] if not NULL, the returned value of
80 *error_estimate is an estimate of the error in the sum.
81 Returns:
82 Total of the sum.
83 Remarks:
84 You can get subtotals by mixing calls to Plus() and Total().
85 In delicate sums, some precision may be lost in the final
86 total if you call Total() to calculate subtotals.
87 */
88 double Total( double* error_estimate = NULL );
89
90 /*
91 Returns:
92 Number of summands.
93 */
94 int SummandCount() const;
95
96private:
97 enum {
98 sum1_max_count=256,
99 sum2_max_count=512,
100 sum3_max_count=1024
101 };
102 double m_sum_err;
103 double m_pos_sum;
104 double m_neg_sum;
105
106 int m_zero_count; // number of zeros added
107 int m_pos_count; // number of positive numbers added
108 int m_neg_count; // number of negative numbers added
109
110 int m_pos_sum1_count;
111 int m_pos_sum2_count;
112 int m_pos_sum3_count;
113 double m_pos_sum1[sum1_max_count];
114 double m_pos_sum2[sum2_max_count];
115 double m_pos_sum3[sum3_max_count];
116
117 int m_neg_sum1_count;
118 int m_neg_sum2_count;
119 int m_neg_sum3_count;
120 double m_neg_sum1[sum1_max_count];
121 double m_neg_sum2[sum2_max_count];
122 double m_neg_sum3[sum3_max_count];
123
124 double SortAndSum( int, double* );
125};
126
127/*
128Description:
129 Abstract function with an arbitrary number of parameters
130 and values. ON_Evaluator is used to pass functions to
131 local solvers.
132*/
133class ON_CLASS ON_Evaluator
134{
135public:
136
137 /*
138 Description:
139 Construction of the class for a function that takes
140 parameter_count input functions and returns
141 value_count values. If the domain is infinite, pass
142 a NULL for the domain[] and periodic[] arrays. If
143 the domain is finite, pass a domain[] array with
144 parameter_count increasing intervals. If one or more of
145 the parameters is periodic, pass the fundamental domain
146 in the domain[] array and a true in the periodic[] array.
147 Parameters:
148 parameter_count - [in] >= 1. Number of input parameters
149 value_count - [in] >= 1. Number of output values.
150 domain - [in] If not NULL, then this is an array
151 of parameter_count increasing intervals
152 that defines the domain of the function.
153 periodic - [in] if not NULL, then this is an array of
154 parameter_count bools where b[i] is true if
155 the i-th parameter is periodic. Valid
156 increasing finite domains must be specificed
157 when this parameter is not NULL.
158 */
160 int parameter_count,
161 int value_count,
162 const ON_Interval* domain,
163 const bool* periodic
164 );
165
166 virtual ~ON_Evaluator();
167
168 /*
169 Description:
170 Evaluate the function that takes m_parameter_count parameters
171 and returns a m_value_count dimensional point.
172 Parameters:
173 parameters - [in] array of m_parameter_count evaluation parameters
174 values - [out] array of m_value_count function values
175 jacobian - [out] If NULL, simply evaluate the value of the function.
176 If not NULL, this is the jacobian of the function.
177 jacobian[i][j] = j-th partial of the i-th value
178 0 <= i < m_value_count,
179 0 <= j < m_parameter_count
180 If not NULL, then all the memory for the
181 jacobian is allocated, you just need to fill
182 in the answers.
183 Example:
184 If f(u,v) = square of the distance from a fixed point P to a
185 surface evaluated at (u,v), then
186
187 values[0] = (S-P)o(S-P)
188 jacobian[0] = ( 2*(Du o (S-P)), 2*(Dv o (S-P)) )
189
190 where S, Du, Dv = surface point and first partials evaluated
191 at u=parameters[0], v = parameters[1].
192
193 If the function takes 3 parameters, say (x,y,z), and returns
194 two values, say f(x,y,z) and g(z,y,z), then
195
196 values[0] = f(x,y,z)
197 values[1] = g(x,y,z)
198
199 jacobian[0] = (DfDx, DfDy, DfDz)
200 jacobian[1] = (DgDx, DgDy, DgDz)
201
202 where dfx denotes the first partial of f with respect to x.
203
204 Returns:
205 0 = unable to evaluate
206 1 = successful evaluation
207 2 = found answer, terminate search
208 */
209 virtual int Evaluate(
210 const double* parameters,
211 double* values,
212 double** jacobian
213 ) = 0;
214
215 /*
216 Description:
217 OPTIONAL ability to evaluate the hessian in the case when
218 m_value_count is one. If your function has more that
219 one value or it is not feasable to evaluate the hessian,
220 then do not override this function. The default implementation
221 returns -1.
222 Parameters:
223 parameters - [in] array of m_parameter_count evaluation parameters
224 value - [out] value of the function (one double)
225 gradient - [out] The gradient of the function. This is a vector
226 of length m_parameter_count; gradient[i] is
227 the first partial of the function with respect to
228 the i-th parameter.
229 hessian - [out] The hessian of the function. This is an
230 m_parameter_count x m_parameter_count
231 symmetric matrix: hessian[i][j] is the
232 second partial of the function with respect
233 to the i-th and j-th parameters. The evaluator
234 is responsible for filling in both the upper
235 and lower triangles. Since the matrix is
236 symmetrix, you should do something like evaluate
237 the upper triangle and copy the values to the
238 lower tiangle.
239 Returns:
240 -1 = Hessian evaluation not available.
241 0 = unable to evaluate
242 1 = successful evaluation
243 2 = found answer, terminate search
244 */
245 virtual int EvaluateHessian(
246 const double* parameters,
247 double* value,
248 double* gradient,
249 double** hessian
250 );
251
252 // Number of the function's input parameters. This number
253 // is >= 1 and is specified in the constructor.
255
256 // Number of the function's output values. This number
257 // is >= 1 and is specified in the constructor.
258 const int m_value_count;
259
260 /*
261 Description:
262 Functions can have finite or infinite domains. Finite domains
263 are specified by passing the domain[] array to the constructor
264 or filling in the m_domain[] member variable. If
265 m_domain.Count() == m_parameter_count > 0, then the function
266 has finite domains.
267 Returns:
268 True if the domain of the function is finite.
269 */
270 bool FiniteDomain() const;
271
272 /*
273 Description:
274 If a function has a periodic parameter, then the m_domain
275 interval for that parameter is the fundamental domain and
276 the m_bPeriodicParameter bool for that parameter is true.
277 A parameter is periodic if, and only if,
278 m_domain.Count() == m_parameter_count, and
279 m_bPeriodicParameter.Count() == m_parameter_count, and
280 m_bPeriodicParameter[parameter_index] is true.
281 Returns:
282 True if the function parameter is periodic.
283 */
285 int parameter_index
286 ) const;
287
288 /*
289 Description:
290 If a function has a periodic parameter, then the m_domain
291 interval for that parameter is the fundamental domain and
292 the m_bPeriodicParameter bool for that parameter is true.
293 A parameter is periodic if, and only if,
294 m_domain.Count() == m_parameter_count, and
295 m_bPeriodicParameter.Count() == m_parameter_count, and
296 m_bPeriodicParameter[parameter_index] is true.
297 Returns:
298 The domain of the parameter. If the domain is infinite,
299 the (-1.0e300, +1.0e300) is returned.
300 */
302 int parameter_index
303 ) const;
304
305
306 // If the function has a finite domain or periodic
307 // parameters, then m_domain[] is an array of
308 // m_parameter_count finite increasing intervals.
310
311 // If the function has periodic parameters, then
312 // m_bPeriodicParameter[] is an array of m_parameter_count
313 // bools. If m_bPeriodicParameter[i] is true, then
314 // the i-th parameter is periodic and m_domain[i] is
315 // the fundamental domain for that parameter.
317
318private:
319 ON_Evaluator(); // prohibit default constructor
320 ON_Evaluator& operator=(const ON_Evaluator&); // prohibit operator= (can't copy const members)
321};
322
323/*
324Description:
325 Test a double to make sure it is a valid number.
326Returns:
327 True if x != ON_UNSET_VALUE and _finite(x) is true.
328*/
329ON_DECL
330bool ON_IsValid( double x );
331
332ON_DECL
333bool ON_IsValidFloat( float x );
334
335/*
336class ON_CLASS ON_TimeLimit
337{
338 ON_TimeLimit();
339 ON_TimeLimit(ON__UINT64 time_limit_seconds);
340 void SetTimeLimit(ON__UINT64 time_limit_seconds);
341 bool Continue() const;
342 bool IsSet() const;
343private:
344 ON__UINT64 m_time_limit[2];
345};
346*/
347
348// The ON_IS_FINITE and ON_IS_VALID defines are much faster
349// than calling ON_IsValid(), but need to be used when
350// the macro expansion works.
351
352#if defined(ON_LITTLE_ENDIAN)
353
354// works on little endian CPUs with IEEE doubles
355#define ON_IS_FINITE(x) (0x7FF0 != (*((unsigned short*)(&x) + 3) & 0x7FF0))
356#define ON_IS_VALID(x) (x != ON_UNSET_VALUE && 0x7FF0 != (*((unsigned short*)(&x) + 3) & 0x7FF0))
357#define ON_IS_VALID_FLOAT(x) (x != ON_UNSET_FLOAT)
358//TODO - ADD FAST ugly bit check#define ON_IS_VALID_FLOAT(x) (x != ON_UNSET_FLOAT && 0x7FF0 != (*((unsigned short*)(&x) + 3) & 0x7FF0))
359
360#elif defined(ON_BIG_ENDIAN)
361
362// works on big endian CPUs with IEEE doubles
363#define ON_IS_FINITE(x) (0x7FF0 != (*((unsigned short*)(&x)) & 0x7FF0))
364#define ON_IS_VALID(x) (x != ON_UNSET_VALUE && 0x7FF0 != (*((unsigned short*)(&x)) & 0x7FF0))
365#define ON_IS_VALID_FLOAT(x) (x != ON_UNSET_FLOAT)
366//TODO - ADD FAST ugly bit check#define ON_IS_VALID_FLOAT(x) (x != ON_UNSET_FLOAT && 0x7FF0 != (*((unsigned short*)(&x) + 3) & 0x7FF0))
367
368#else
369
370// Returns true if x is a finite double. Specifically,
371// _finite returns a nonzero value (true) if its argument x
372// is not infinite, that is, if -INF < x < +INF.
373// It returns 0 (false) if the argument is infinite or a NaN.
374//
375// If you are trying to compile opennurbs on a platform
376// that does not support finite(), then see if you can
377// use _fpclass(), fpclass(), _isnan(), or isnan(). If
378// you can't find anything, then just set this
379// function to return true.
380
381#if defined(_GNU_SOURCE)
382// if you are using an older version of gcc, use finite()
383//#define ON_IS_FINITE(x) (finite(x)?true:false)
384#define ON_IS_FINITE(x) (isfinite(x)?true:false)
385#else
386#define ON_IS_FINITE(x) (_finite(x)?true:false)
387#endif
388
389#define ON_IS_VALID(x) (x != ON_UNSET_VALUE && ON_IS_FINITE(x))
390#define ON_IS_VALID_FLOAT(x) (x != ON_UNSET_FLOAT && ON_IS_FINITE(x))
391
392#endif
393
394
395ON_DECL
396float ON_ArrayDotProduct( // returns AoB
397 int, // size of arrays (can be zero)
398 const float*, // A[]
399 const float* // B[]
400 );
401
402ON_DECL
403void ON_ArrayScale(
404 int, // size of arrays (can be zero)
405 float, // a
406 const float*, // A[]
407 float* // returns a*A[]
408 );
409
410ON_DECL
411void ON_Array_aA_plus_B(
412 int, // size of arrays (can be zero)
413 float, // a
414 const float*, // A[]
415 const float*, // B[]
416 float* // returns a*A[] + B[]
417 );
418
419ON_DECL
420double ON_ArrayDotProduct( // returns AoB
421 int, // size of arrays (can be zero)
422 const double*, // A[]
423 const double* // B[]
424 );
425
426ON_DECL
427double ON_ArrayDotDifference( // returns A o ( B - C )
428 int, // size of arrays (can be zero)
429 const double*, // A[]
430 const double*, // B[]
431 const double* // C[]
432 );
433
434ON_DECL
435double ON_ArrayMagnitude( // returns sqrt(AoA)
436 int, // size of arrays (can be zero)
437 const double* // A[]
438 );
439
440ON_DECL
441double ON_ArrayMagnitudeSquared( // returns AoA
442 int, // size of arrays (can be zero)
443 const double* // A[]
444 );
445
446ON_DECL
447double ON_ArrayDistance( // returns sqrt((A-B)o(A-B))
448 int, // size of arrays (can be zero)
449 const double*, // A[]
450 const double* // B[]
451 );
452
453ON_DECL
454double ON_ArrayDistanceSquared( // returns (A-B)o(A-B)
455 int, // size of arrays (can be zero)
456 const double*, // A[]
457 const double* // B[]
458 );
459
460ON_DECL
461void ON_ArrayScale(
462 int, // size of arrays (can be zero)
463 double, // a
464 const double*, // A[]
465 double* // returns a*A[]
466 );
467
468ON_DECL
469void ON_Array_aA_plus_B(
470 int, // size of arrays (can be zero)
471 double, // a
472 const double*, // A[]
473 const double*, // B[]
474 double* // returns a*A[] + B[]
475 );
476
477ON_DECL
478int ON_SearchMonotoneArray( // find a value in an increasing array
479 // returns -1: t < array[0]
480 // i: array[i] <= t < array[i+1] ( 0 <= i < length-1 )
481 // length-1: t == array[length-1]
482 // length: t >= array[length-1]
483 const double*, // array[]
484 int, // length of array
485 double // t = value to search for
486 );
487
488
489/*
490Description:
491 Compute a binomial coefficient.
492Parameters:
493 i - [in]
494 j - [in]
495Returns:
496 (i+j)!/(i!j!), if 0 <= i and 0 <= j, and 0 otherwise.
497See Also:
498 ON_TrinomialCoefficient()
499Remarks:
500 If (i+j) <= 52, this function is fast and returns the exact
501 value of the binomial coefficient.
502
503 For (i+j) > 52, the coefficient is computed recursively using
504 the formula bc(i,j) = bc(i-1,j) + bc(i,j-1).
505 For (i+j) much larger than 60, this is inefficient.
506 If you need binomial coefficients for large i and j, then you
507 should probably be using something like Stirling's Formula.
508 (Look up "Stirling" or "Gamma function" in a calculus book.)
509*/
510ON_DECL
511double ON_BinomialCoefficient(
512 int i,
513 int j
514 );
515
516
517/*
518Description:
519 Compute a trinomial coefficient.
520Parameters:
521 i - [in]
522 j - [in]
523 k - [in]
524Returns:
525 (i+j+k)!/(i!j!k!), if 0 <= i, 0 <= j and 0<= k, and 0 otherwise.
526See Also:
527 ON_BinomialCoefficient()
528Remarks:
529 The trinomial coefficient is computed using the formula
530
531 (i+j+k)! (i+j+k)! (j+k)!
532 -------- = -------- * -------
533 i! j! k! i! (j+k)! j! k!
534
535 = ON_BinomialCoefficient(i,j+k)*ON_BinomialCoefficient(j,k)
536
537*/
538ON_DECL
539double ON_TrinomialCoefficient(
540 int i,
541 int j,
542 int k
543 );
544
545
546ON_DECL
547ON_BOOL32 ON_GetParameterTolerance(
548 double, double, // domain
549 double, // parameter in domain
550 double*, double* // parameter tolerance (tminus, tplus) returned here
551 );
552
553
554ON_DECL
555ON_BOOL32 ON_IsValidPointList(
556 int, // dim
557 ON_BOOL32, // true for homogeneous rational points
558 int, // count
559 int, // stride
560 const float*
561 );
562
563ON_DECL
564ON_BOOL32 ON_IsValidPointList(
565 int, // dim
566 ON_BOOL32, // true for homogeneous rational points
567 int, // count
568 int, // stride
569 const double*
570 );
571
572/*
573Description:
574 Determine if a list of points is planar.
575Parameters:
576 bRational - [in]
577 false if the points are euclidean (x,y,z)
578 true if the points are homogeneous rational (x,y,z,w)
579 point_count - [in]
580 number of points
581 point_stride - [in]
582 number of doubles between point x coordinates
583 first point's x coordinate = points[0],
584 second point's x coordinate = points[point_stride],...
585 points - [in]
586 point coordinates (3d or 4d homogeneous rational)
587 boxMin - [in]
588 boxMax - [in]
589 optional 3d bounding box - pass nulls if not readily available
590 tolerance - [in] >= 0.0
591 plane_equation0 - [in]
592 If you want to test for planarity in a specific plane,
593 pass the plane equation in here. If you want to find
594 a plane containing the points, pass null here.
595 plane_equation - [out]
596 If this point is not null, then the equation of the plane
597 containing the points is retuened here.
598Returns:
599 0 - points are not coplanar to the specified tolerance
600 1 - points are coplanar to the specified tolerance
601 2 - points are colinear to the specified tolerance
602 (in this case, plane_equation is not a unique answer)
603 3 - points are coincident to the specified tolerance
604 (in this case, plane_equation is not a unique answer)
605*/
606ON_DECL
607int ON_IsPointListPlanar(
608 bool bRational,
609 int count,
610 int stride,
611 const double* points,
612 const double* boxMin,
613 const double* boxMax,
614 double tolerance,
615 ON_PlaneEquation* plane_equation
616 );
617
618ON_DECL
619ON_BOOL32 ON_IsValidPointGrid(
620 int, // dim
621 ON_BOOL32, // true for homogeneous rational points
622 int, int, // point_count0, point_count1,
623 int, int, // point_stride0, point_stride1,
624 const double*
625 );
626
627ON_DECL
628bool ON_ReversePointList(
629 int, // dim
630 ON_BOOL32, // true for homogeneous rational points
631 int, // count
632 int, // stride
633 double*
634 );
635
636ON_DECL
637ON_BOOL32 ON_ReversePointGrid(
638 int, // dim
639 ON_BOOL32, // true for homogeneous rational points
640 int, int, // point_count0, point_count1,
641 int, int, // point_stride0, point_stride1,
642 double*,
643 int // dir = 0 or 1
644 );
645
646ON_DECL
647bool ON_SwapPointListCoordinates(
648 int, // count
649 int, // stride
650 float*,
651 int, int // coordinates to swap
652 );
653
654ON_DECL
655bool ON_SwapPointListCoordinates(
656 int, // count
657 int, // stride
658 double*,
659 int, int // coordinates to swap
660 );
661
662ON_DECL
663ON_BOOL32 ON_SwapPointGridCoordinates(
664 int, int, // point_count0, point_count1,
665 int, int, // point_stride0, point_stride1,
666 double*,
667 int, int // coordinates to swap
668 );
669
670ON_DECL
671bool ON_TransformPointList(
672 int, // dim
673 ON_BOOL32, // true for homogeneous rational points
674 int, // count
675 int, // stride
676 float*,
677 const ON_Xform&
678 );
679
680ON_DECL
681bool ON_TransformPointList(
682 int, // dim
683 ON_BOOL32, // true for homogeneous rational points
684 int, // count
685 int, // stride
686 double*,
687 const ON_Xform&
688 );
689
690ON_DECL
691ON_BOOL32 ON_TransformPointGrid(
692 int, // dim
693 ON_BOOL32, // true for homogeneous rational points
694 int, int, // point_count0, point_count1,
695 int, int, // point_stride0, point_stride1,
696 double*,
697 const ON_Xform&
698 );
699
700ON_DECL
701ON_BOOL32 ON_TransformVectorList(
702 int, // dim
703 int, // count
704 int, // stride
705 float*,
706 const ON_Xform&
707 );
708
709ON_DECL
710ON_BOOL32 ON_TransformVectorList(
711 int, // dim
712 int, // count
713 int, // stride
714 double*,
715 const ON_Xform&
716 );
717
718/*
719Parameters:
720 dim - [in]
721 >= 1
722 is_rat - [in]
723 true if the points are rational and points[dim] is the "weight"
724 pointA - [in]
725 pointB - [in]
726 point coordinates
727Returns:
728 True if the input is valid and for each coordinate pair,
729 |a-b| <= ON_ZERO_TOLERANCE
730 or |a-b| <= (fabs(a)+fabs(b))*ON_RELATIVE_TOLERANCE.
731 False otherwise.
732*/
733ON_DECL
734bool ON_PointsAreCoincident(
735 int dim,
736 int is_rat,
737 const double* pointA,
738 const double* pointB
739 );
740
741/*
742Description
743 See ON_PointsAreCoincident() for a description of when opennurbs
744 considers two points to be conincident.
745Parameters:
746 dim - [in]
747 >= 1
748 is_rat - [in]
749 true if the points are rational and points[dim] is the "weight"
750 point_count - [in]
751 number of points >= 2
752 point_stride - [in]
753 >= (0 != is_rat) ? (dim+1) : dim
754 points - [in]
755 point coordinates
756Returns:
757 True if the first and last points are coincident and all other
758 points in the list are coincident with the previous point.
759 False if there are points that are not coincident or
760 point_count < 2 or other input parameters are invalid.
761*/
762ON_DECL
763bool ON_PointsAreCoincident(
764 int dim,
765 int is_rat,
766 int point_count,
767 int point_stride,
768 const double* points
769 );
770
771ON_DECL
772int ON_ComparePoint( // returns
773 // -1: first < second
774 // 0: first == second
775 // +1: first > second
776 int dim, // dim (>=0)
777 ON_BOOL32 israt, // true for rational CVs
778 const double* cv0, // first CV
779 const double* cv1 // secont CV
780 );
781
782ON_DECL
783int ON_ComparePointList( // returns
784 // -1: first < second
785 // 0: first == second
786 // +1: first > second
787 int, // dim (>=0)
788 ON_BOOL32, // true for rational CVs
789 int, // count
790 // first point list
791 int, // stride
792 const double*, // point
793 // second point list
794 int, // stride
795 const double* // point
796 );
797
798ON_DECL
799ON_BOOL32 ON_IsPointListClosed(
800 int, // dim
801 int, // true for homogeneos rational points
802 int, // count
803 int, // stride
804 const double*
805 );
806
807ON_DECL
808ON_BOOL32 ON_IsPointGridClosed(
809 int, // dim
810 ON_BOOL32, // true for homogeneous rational points
811 int, int, // point_count0, point_count1,
812 int, int, // point_stride0, point_stride1,
813 const double*,
814 int // dir = 0 or 1
815 );
816
817ON_DECL
818int ON_SolveQuadraticEquation( // solve a*X^2 + b*X + c = 0
819 // returns 0: two distinct real roots (r0 < r1)
820 // 1: one real root (r0 = r1)
821 // 2: two complex conjugate roots (r0 +/- (r1)*sqrt(-1))
822 // -1: failure - a = 0, b != 0 (r0 = r1 = -c/b)
823 // -2: failure - a = 0, b = 0 c != 0 (r0 = r1 = 0.0)
824 // -3: failure - a = 0, b = 0 c = 0 (r0 = r1 = 0.0)
825 double, double, double, // a, b, c
826 double*, double* // roots r0 and r1 returned here
827 );
828
829ON_DECL
830ON_BOOL32 ON_SolveTriDiagonal( // solve TriDiagMatrix( a,b,c )*X = d
831 int, // dimension of d and X (>=1)
832 int, // number of equations (>=2)
833 double*, // a[n-1] = sub-diagonal (a is modified)
834 const double*, // b[n] = diagonal
835 double*, // c[n-1] = supra-diagonal
836 const double*, // d[n*dim]
837 double* // X[n*dim] = unknowns
838 );
839
840// returns rank - if rank != 2, system is under determined
841// If rank = 2, then solution to
842//
843// a00*x0 + a01*x1 = b0,
844// a10*x0 + a11*x1 = b1
845//
846// is returned
847ON_DECL
848int ON_Solve2x2(
849 double, double, // a00 a01 = first row of 2x2 matrix
850 double, double, // a10 a11 = second row of 2x2 matrix
851 double, double, // b0 b1
852 double*, double*, // x0, x1 if not NULL, then solution is returned here
853 double* // if not NULL, then pivot_ratio returned here
854 );
855
856// Description:
857// Solves a system of 3 linear equations and 2 unknowns.
858//
859// x*col0[0] + y*col1[0] = d0
860// x*col0[1] + y*col1[1] = d0
861// x*col0[2] + y*col1[2] = d0
862//
863// Parameters:
864// col0 - [in] coefficents for "x" unknown
865// col1 - [in] coefficents for "y" unknown
866// d0 - [in] constants
867// d1 - [in]
868// d2 - [in]
869// x - [out]
870// y - [out]
871// error - [out]
872// pivot_ratio - [out]
873//
874// Returns:
875// rank of the system.
876// If rank != 2, system is under determined
877// If rank = 2, then the solution is
878//
879// (*x)*[col0] + (*y)*[col1]
880// + (*error)*((col0 X col1)/|col0 X col1|)
881// = (d0,d1,d2).
882ON_DECL
883int ON_Solve3x2(
884 const double[3], // col0
885 const double[3], // col1
886 double, // d0
887 double, // d1
888 double, // d2
889 double*, // x
890 double*, // y
891 double*, // error
892 double* // pivot_ratio
893 );
894
895/*
896Description:
897 Use Gauss-Jordan elimination with full pivoting to solve
898 a system of 3 linear equations and 3 unknowns(x,y,z)
899
900 x*row0[0] + y*row0[1] + z*row0[2] = d0
901 x*row1[0] + y*row1[1] + z*row1[2] = d1
902 x*row2[0] + y*row2[1] + z*row2[2] = d2
903
904Parameters:
905 row0 - [in] first row of 3x3 matrix
906 row1 - [in] second row of 3x3 matrix
907 row2 - [in] third row of 3x3 matrix
908 d0 - [in]
909 d1 - [in]
910 d2 - [in] (d0,d1,d2) right hand column of system
911 x_addr - [in] first unknown
912 y_addr - [in] second unknown
913 z_addr - [in] third unknown
914 pivot_ratio - [out] if not NULL, the pivot ration is
915 returned here. If the pivot ratio is "small",
916 then the matrix may be singular or ill
917 conditioned. You should test the results
918 before you use them. "Small" depends on the
919 precision of the input coefficients and the
920 use of the solution. If you can't figure out
921 what "small" means in your case, then you
922 must check the solution before you use it.
923
924Returns:
925 The rank of the 3x3 matrix (0,1,2, or 3)
926 If ON_Solve3x3() is successful (returns 3), then
927 the solution is returned in
928 (*x_addr, *y_addr, *z_addr)
929 and *pivot_ratio = min(|pivots|)/max(|pivots|).
930 If the return code is < 3, then (0,0,0) is returned
931 as the "solution".
932
933See Also:
934 ON_Solve2x2
935 ON_Solve3x2
936 ON_Solve4x4
937*/
938ON_DECL
939int ON_Solve3x3(
940 const double row0[3],
941 const double row1[3],
942 const double row2[3],
943 double d0,
944 double d1,
945 double d2,
946 double* x_addr,
947 double* y_addr,
948 double* z_addr,
949 double* pivot_ratio
950 );
951
952/*
953Description:
954 Use Gauss-Jordan elimination with full pivoting to solve
955 a system of 4 linear equations and 4 unknowns(x,y,z,w)
956
957 x*row0[0] + y*row0[1] + z*row0[2] + w*row0[3] = d0
958 x*row1[0] + y*row1[1] + z*row1[2] + w*row1[3] = d1
959 x*row2[0] + y*row2[1] + z*row2[2] + w*row2[3] = d2
960 x*row3[0] + y*row3[1] + z*row3[2] + w*row3[2] = d3
961
962Parameters:
963 row0 - [in] first row of 4x4 matrix
964 row1 - [in] second row of 4x4 matrix
965 row2 - [in] third row of 4x4 matrix
966 row3 - [in] forth row of 4x4 matrix
967 d0 - [in]
968 d1 - [in]
969 d2 - [in]
970 d3 - [in] (d0,d1,d2,d3) right hand column of system
971 x_addr - [in] first unknown
972 y_addr - [in] second unknown
973 z_addr - [in] third unknown
974 w_addr - [in] forth unknown
975 pivot_ratio - [out] if not NULL, the pivot ration is
976 returned here. If the pivot ratio is "small",
977 then the matrix may be singular or ill
978 conditioned. You should test the results
979 before you use them. "Small" depends on the
980 precision of the input coefficients and the
981 use of the solution. If you can't figure out
982 what "small" means in your case, then you
983 must check the solution before you use it.
984
985Returns:
986 The rank of the 4x4 matrix (0,1,2,3, or 4)
987 If ON_Solve4x4() is successful (returns 4), then
988 the solution is returned in
989 (*x_addr, *y_addr, *z_addr, *w_addr)
990 and *pivot_ratio = min(|pivots|)/max(|pivots|).
991 If the return code is < 4, then, it a solution exists,
992 on is returned. However YOU MUST CHECK THE SOLUTION
993 IF THE RETURN CODE IS < 4.
994
995See Also:
996 ON_Solve2x2
997 ON_Solve3x2
998 ON_Solve3x3
999*/
1000ON_DECL
1001int
1002ON_Solve4x4(
1003 const double row0[4],
1004 const double row1[4],
1005 const double row2[4],
1006 const double row3[4],
1007 double d0,
1008 double d1,
1009 double d2,
1010 double d3,
1011 double* x_addr,
1012 double* y_addr,
1013 double* z_addr,
1014 double* w_addr,
1015 double* pivot_ratio
1016 );
1017
1018/*
1019Description:
1020 Use Gauss-Jordan elimination to find a numerical
1021 solution to M*X = B where M is a n x n matrix,
1022 B is a known n-dimensional vector and X is
1023 an unknown.
1024Paramters:
1025 bFullPivot - [in] if true, full pivoting is used,
1026 otherwise partial pivoting is used. In rare
1027 cases full pivoting can produce a more accurate
1028 answer and never produces a less accurate answer.
1029 However full pivoting is slower. If speed is an
1030 issue, then experiement with bFullPivot=false
1031 and see if it makes a difference. Otherwise,
1032 set it to true.
1033 bNormalize - [in]
1034 If bNormalize is true, then the rows of the
1035 matrix are scaled so the sum of their squares
1036 is one. This doesn't make the solution more
1037 accurate but in some cases it makes the pivot
1038 ratio more meaningful. Set bNormalize to
1039 false unless you have a reason for setting it
1040 to true.
1041 n - [in] size of the matrix and vectors.
1042 M - [in] n x n matrix. The values in M are
1043 changed as the solution is calculated.
1044 If you need to preserve M for future use,
1045 pass in a copy.
1046 B - [in] n-dimensional vector. The values in
1047 B are changed as the solution is calculated.
1048 If you need to preserve B for future use,
1049 pass in a copy.
1050 X - [out] solution to M*X = B.
1051Returns:
1052 If the returned value is <= 0.0, the input matrix
1053 has rank < n and no solution is returned in X.
1054 If the returned value is > 0.0, then a solution is
1055 returned in X and the returned value is the ratio
1056 (minimum pivot)/(maximum pivot). This value is
1057 called the pivot ratio and will be denoted "pr"
1058 the discussion below. If pr <= 1e-15, then
1059 M was nearly degenerate and the solution should be
1060 used with caution. If an accurate solution is
1061 critcial, then check the solution anytime pr <= 1e-10
1062 In general, the difference between M*X and B will be
1063 reasonably small. However, when the pr is small
1064 there tend to be vector E, substantually different
1065 from zero, such that M*(X+E) - B is also reasonably
1066 small.
1067See Also:
1068 ON_Solve2x2
1069 ON_Solve3x3
1070 ON_Solve4x4
1071 ON_Solve3x2
1072*/
1073ON_DECL
1074double ON_SolveNxN(bool bFullPivot, bool bNormalize, int n, double* M[], double B[], double X[]);
1075
1076
1077// return false if determinant is (nearly) singular
1078ON_DECL
1079ON_BOOL32 ON_EvJacobian(
1080 double, // ds o ds
1081 double, // ds o dt
1082 double, // dt o dt
1083 double* // jacobian = determinant ( ds_o_ds dt_o_dt / ds_o_dt ds_o_dt )
1084 );
1085
1086/*
1087Description:
1088 Finds scalars x and y so that the component of V in the plane
1089 of A and B is x*A + y*B.
1090Parameters:
1091 V - [in]
1092 A - [in] nonzero and not parallel to B
1093 B - [in] nonzero and not parallel to A
1094 x - [out]
1095 y - [out]
1096Returns:
1097 1 - The rank of the problem is 2. The decomposition is unique.
1098 0 - The rank less than 2. Either there is no solution or there
1099 are infinitely many solutions.
1100
1101See Also:
1102 ON_Solve2x2
1103*/
1104ON_DECL
1105int ON_DecomposeVector(
1106 const ON_3dVector& V,
1107 const ON_3dVector& A,
1108 const ON_3dVector& B,
1109 double* x, double* y
1110 );
1111
1112
1113/*
1114Description:
1115 Evaluate partial derivatives of surface unit normal
1116Parameters:
1117 ds - [in]
1118 dt - [in] surface first partial derivatives
1119 dss - [in]
1120 dst - [in]
1121 dtt - [in] surface second partial derivatives
1122 ns - [out]
1123 nt - [out] First partial derivatives of surface unit normal
1124 (If the Jacobian is degenerate, ns and nt are set to zero.)
1125Returns:
1126 true if Jacobian is nondegenerate
1127 false if Jacobian is degenerate
1128*/
1129ON_DECL
1130ON_BOOL32 ON_EvNormalPartials(
1131 const ON_3dVector& ds,
1132 const ON_3dVector& dt,
1133 const ON_3dVector& dss,
1134 const ON_3dVector& dst,
1135 const ON_3dVector& dtt,
1136 ON_3dVector& ns,
1137 ON_3dVector& nt
1138 );
1139
1140ON_DECL
1141ON_BOOL32
1142ON_Pullback3dVector( // use to pull 3d vector back to surface parameter space
1143 const ON_3dVector&, // 3d vector
1144 double, // signed distance from vector location to closet point on surface
1145 // < 0 if point is below with respect to Du x Dv
1146 const ON_3dVector&, // ds surface first partials
1147 const ON_3dVector&, // dt
1148 const ON_3dVector&, // dss surface 2nd partials
1149 const ON_3dVector&, // dst (used only when dist != 0)
1150 const ON_3dVector&, // dtt
1151 ON_2dVector& // pullback
1152 );
1153
1154ON_DECL
1155ON_BOOL32
1156ON_GetParameterTolerance(
1157 double, // t0 domain
1158 double, // t1
1159 double, // t parameter in domain
1160 double*, // tminus parameter tolerance (tminus, tplus) returned here
1161 double* // tplus
1162 );
1163
1164
1165ON_DECL
1166ON_BOOL32 ON_EvNormal(
1167 int, // limit_dir 0=default,1=from quadrant I, 2 = from quadrant II, ...
1168 const ON_3dVector&, const ON_3dVector&, // first partials (Du,Dv)
1169 const ON_3dVector&, const ON_3dVector&, const ON_3dVector&, // optional second partials (Duu, Duv, Dvv)
1170 ON_3dVector& // unit normal returned here
1171 );
1172
1173// returns false if first returned tangent is zero
1174ON_DECL
1175bool ON_EvTangent(
1176 const ON_3dVector&, // first derivative
1177 const ON_3dVector&, // second derivative
1178 ON_3dVector& // Unit tangent returned here
1179 );
1180
1181// returns false if first derivtive is zero
1182ON_DECL
1183ON_BOOL32 ON_EvCurvature(
1184 const ON_3dVector&, // first derivative
1185 const ON_3dVector&, // second derivative
1186 ON_3dVector&, // Unit tangent returned here
1187 ON_3dVector& // Curvature returned here
1188 );
1189
1190ON_DECL
1191ON_BOOL32 ON_EvPrincipalCurvatures(
1192 const ON_3dVector&, // Ds,
1193 const ON_3dVector&, // Dt,
1194 const ON_3dVector&, // Dss,
1195 const ON_3dVector&, // Dst,
1196 const ON_3dVector&, // Dtt,
1197 const ON_3dVector&, // N, // unit normal to surface (use ON_EvNormal())
1198 double*, // gauss, // = Gaussian curvature = kappa1*kappa2
1199 double*, // mean, // = mean curvature = (kappa1+kappa2)/2
1200 double*, // kappa1, // = largest principal curvature value (may be negative)
1201 double*, // kappa2, // = smallest principal curvature value (may be negative)
1202 ON_3dVector&, // K1, // kappa1 unit principal curvature direction
1203 ON_3dVector& // K2 // kappa2 unit principal curvature direction
1204 // output K1,K2,N is right handed frame
1205 );
1206
1207ON_DECL
1208ON_BOOL32 ON_EvPrincipalCurvatures(
1209 const ON_3dVector&, // Ds,
1210 const ON_3dVector&, // Dt,
1211 double l, // Dss*N Second fundamental form coefficients
1212 double m, // Dst*N,
1213 double n, // Dtt*N,
1214 const ON_3dVector&, // N, // unit normal to surface (use ON_EvNormal())
1215 double*, // gauss, // = Gaussian curvature = kappa1*kappa2
1216 double*, // mean, // = mean curvature = (kappa1+kappa2)/2
1217 double*, // kappa1, // = largest principal curvature value (may be negative)
1218 double*, // kappa2, // = smallest principal curvature value (may be negative)
1219 ON_3dVector&, // K1, // kappa1 unit principal curvature direction
1220 ON_3dVector& // K2 // kappa2 unit principal curvature direction
1221 // output K1,K2,N is right handed frame
1222 );
1223
1224/*
1225Description:
1226 Evaluate sectional curvature from surface derivatives and
1227 section plane normal.
1228Parameters:
1229 S10, S01 - [in]
1230 surface 1st partial derivatives
1231 S20, S11, S02 - [in]
1232 surface 2nd partial derivatives
1233 planeNormal - [in]
1234 unit normal to section plane
1235 K - [out] Sectional curvature
1236 Curvature of the intersection curve of the surface
1237 and plane through the surface point where the partial
1238 derivatives were evaluationed.
1239Returns:
1240 True if successful.
1241 False if first partials are not linearly independent, in
1242 which case the K is set to zero.
1243*/
1244ON_DECL
1245bool ON_EvSectionalCurvature(
1246 const ON_3dVector& S10,
1247 const ON_3dVector& S01,
1248 const ON_3dVector& S20,
1249 const ON_3dVector& S11,
1250 const ON_3dVector& S02,
1251 const ON_3dVector& planeNormal,
1252 ON_3dVector& K
1253 );
1254
1255
1256ON_DECL
1257ON_3dVector ON_NormalCurvature(
1258 const ON_3dVector&, // surface 1rst partial (Ds)
1259 const ON_3dVector&, // surface 1rst partial (Dt)
1260 const ON_3dVector&, // surface 1rst partial (Dss)
1261 const ON_3dVector&, // surface 1rst partial (Dst)
1262 const ON_3dVector&, // surface 1rst partial (Dtt)
1263 const ON_3dVector&, // surface unit normal
1264 const ON_3dVector& // unit tangent direction
1265 );
1266
1267/*
1268Description:
1269 Determing if two curvatrues are different enough
1270 to qualify as a curvature discontinuity.
1271Parameters:
1272 Km - [in]
1273 Kp - [in]
1274 Km and Kp should be curvatures evaluated at the same
1275 parameters using limits from below (minus) and above (plus).
1276 The assumption is that you have already compared the
1277 points and tangents and consider to curve to be G1 at the
1278 point in question.
1279 cos_angle_tolerance - [in]
1280 If the inut value of cos_angle_tolerance >= -1.0
1281 and cos_angle_tolerance <= 1.0 and
1282 Km o Kp < cos_angle_tolerance*|Km|*|Kp|, then
1283 true is returned. Otherwise it is assumed Km and Kp
1284 are parallel. If the curve being tested is nonplanar,
1285 then use something like cos(2*tangent angle tolerance)
1286 for this parameter. If the curve being tested is planar,
1287 then 0.0 will work fine.
1288 curvature_tolerance - [in]
1289 If |Kp-Km| <= curvature_tolerance,
1290 then false is returned, otherwise other tests are used
1291 to determing continuity.
1292 zero_curvature - [in] (ignored if < 2^-110 = 7.7037197787136e-34)
1293 If |K| <= zero_curvature, then K is treated as zero.
1294 When in doubt, use ON_ZERO_CURVATURE_TOLERANCE.
1295 radius_tolerance - [in]
1296 If radius_tolerance >= 0.0 and the difference between the
1297 radii of curvature is >= radius_tolerance, then true
1298 is returned.
1299 relative_tolerance - [in]
1300 If relative_tolerance > 0 and
1301 |(|Km| - |Kp|)|/max(|Km|,|Kp|) > relative_tolerance,
1302 then true is returned. Note that if the curvatures are
1303 nonzero and rm and rp are the radii of curvature, then
1304 |(|Km| - |Kp|)|/max(|Km|,|Kp|) = |rm-rp|/max(rm,rp).
1305 This means the relative_tolerance insures both the scalar
1306 curvature and the radii of curvature agree to the specified
1307 number of decimal places.
1308 When in double use ON_RELATIVE_CURVATURE_TOLERANCE, which
1309 is currently 0.05.
1310Returns:
1311 False if the curvatures should be considered G2.
1312 True if the curvatures are different enough that the curve should be
1313 considered not G2.
1314 In addition to the tests described under the curvature_tolerance and
1315 radius_tolerance checks, other hurestic tests are used.
1316*/
1317ON_DECL
1318bool ON_IsCurvatureDiscontinuity(
1319 const ON_3dVector Km,
1320 const ON_3dVector Kp,
1321 double cos_angle_tolerance,
1322 double curvature_tolerance,
1323 double zero_curvature,
1324 double radius_tolerance,
1325 double relative_tolerance
1326 );
1327
1328ON_DECL
1329bool ON_IsCurvatureDiscontinuity(
1330 const ON_3dVector Km,
1331 const ON_3dVector Kp,
1332 double cos_angle_tolerance,
1333 double curvature_tolerance,
1334 double zero_curvature,
1335 double radius_tolerance
1336 );
1337
1338
1339/*
1340Description:
1341 This function is used to test curvature continuity
1342 in IsContinuous and GetNextDiscontinuity functions
1343 when the continuity parameter is ON::G2_continuous.
1344Parameters:
1345 Km - [in]
1346 Curve's vector curvature evaluated from below
1347 Kp - [in]
1348 Curve's vector curvature evaluated from below
1349Returns:
1350 True if the change from Km to Kp should be considered
1351 G2 continuous.
1352*/
1353ON_DECL
1354bool ON_IsG2CurvatureContinuous(
1355 const ON_3dVector Km,
1356 const ON_3dVector Kp,
1357 double cos_angle_tolerance,
1358 double curvature_tolerance
1359 );
1360
1361/*
1362Description:
1363 This function is used to test curvature continuity
1364 in IsContinuous and GetNextDiscontinuity functions
1365 when the continuity parameter is ON::Gsmooth_continuous.
1366Parameters:
1367 Km - [in]
1368 Curve's vector curvature evaluated from below
1369 Kp - [in]
1370 Curve's vector curvature evaluated from below
1371Returns:
1372 True if the change from Km to Kp should be considered
1373 Gsmooth continuous.
1374*/
1375ON_DECL
1376bool ON_IsGsmoothCurvatureContinuous(
1377 const ON_3dVector Km,
1378 const ON_3dVector Kp,
1379 double cos_angle_tolerance,
1380 double curvature_tolerance
1381 );
1382
1383/*
1384Description:
1385 Test curve continuity from derivative values.
1386Parameters:
1387 c - [in] type of continuity to test for. Read ON::continuity
1388 comments for details.
1389 Pa - [in] point on curve A.
1390 D1a - [in] first derviative of curve A.
1391 D2a - [in] second derviative of curve A.
1392 Pb - [in] point on curve B.
1393 D1b - [in] first derviative of curve B.
1394 D3b - [in] second derviative of curve B.
1395 point_tolerance - [in] if the distance between two points is
1396 greater than point_tolerance, then the curve is not C0.
1397 d1_tolerance - [in] if the difference between two first derivatives is
1398 greater than d1_tolerance, then the curve is not C1.
1399 d2_tolerance - [in] if the difference between two second derivatives is
1400 greater than d2_tolerance, then the curve is not C2.
1401 cos_angle_tolerance - [in] default = cos(1 degree) Used only when
1402 c is ON::G1_continuous or ON::G2_continuous. If the cosine
1403 of the angle between two tangent vectors
1404 is <= cos_angle_tolerance, then a G1 discontinuity is reported.
1405 curvature_tolerance - [in] (default = ON_SQRT_EPSILON) Used only when
1406 c is ON::G2_continuous. If K0 and K1 are curvatures evaluated
1407 from above and below and |K0 - K1| > curvature_tolerance,
1408 then a curvature discontinuity is reported.
1409Returns:
1410 true if the curve has at least the c type continuity at
1411 the parameter t.
1412*/
1413ON_DECL
1414ON_BOOL32 ON_IsContinuous(
1415 ON::continuity c,
1416 ON_3dPoint Pa,
1417 ON_3dVector D1a,
1418 ON_3dVector D2a,
1419 ON_3dPoint Pb,
1420 ON_3dVector D1b,
1421 ON_3dVector D2b,
1422 double point_tolerance=ON_ZERO_TOLERANCE,
1423 double d1_tolerance=ON_ZERO_TOLERANCE,
1424 double d2_tolerance=ON_ZERO_TOLERANCE,
1425 double cos_angle_tolerance=ON_DEFAULT_ANGLE_TOLERANCE_COSINE,
1426 double curvature_tolerance=ON_SQRT_EPSILON
1427 );
1428
1429
1430ON_DECL
1431bool ON_TuneupEvaluationParameter(
1432 int side,
1433 double s0, double s1, // segment domain
1434 double *s // segment parameter
1435 );
1436
1437
1438ON_DECL
1439int ON_Compare2dex( const ON_2dex* a, const ON_2dex* b);
1440
1441ON_DECL
1442int ON_Compare3dex( const ON_3dex* a, const ON_3dex* b);
1443
1444ON_DECL
1445int ON_Compare4dex( const ON_4dex* a, const ON_4dex* b);
1446
1447ON_DECL
1448const ON_2dex* ON_BinarySearch2dexArray(
1449 int key_i,
1450 const ON_2dex* base,
1451 std::size_t nel
1452 );
1453
1454// These simple intersectors are fast and detect transverse intersections.
1455// If the intersection is not a simple transverse case, then they
1456// return false and you will have to use one of the slower but fancier
1457// models.
1458
1459// returns closest points between the two infinite lines
1460ON_DECL
1461bool ON_Intersect(
1462 const ON_Line&,
1463 const ON_Line&,
1464 double*, // parameter on first line
1465 double* // parameter on second line
1466 );
1467
1468// Returns false unless intersection is a single point
1469// If returned parameter is < 0 or > 1, then the line
1470// segment between line.m_point[0] and line.m_point[1]
1471// does not intersect the plane
1472ON_DECL
1473bool ON_Intersect(
1474 const ON_Line&,
1475 const ON_Plane&,
1476 double* // parameter on line
1477 );
1478
1479ON_DECL
1480bool ON_Intersect(
1481 const ON_Plane&,
1482 const ON_Plane&,
1483 ON_Line& // intersection line is returned here
1484 );
1485
1486ON_DECL
1487bool ON_Intersect(
1488 const ON_Plane&,
1489 const ON_Plane&,
1490 const ON_Plane&,
1491 ON_3dPoint& // intersection point is returned here
1492 );
1493
1494// returns 0 = no intersections,
1495// 1 = intersection = single point,
1496// 2 = intersection = circle
1497// If 0 is returned, returned circle has radius=0
1498// and center = point on sphere closest to plane.
1499// If 1 is returned, intersection is a single
1500// point and returned circle has radius=0
1501// and center = intersection point on sphere.
1502ON_DECL
1503int ON_Intersect(
1504 const ON_Plane&, const ON_Sphere&, ON_Circle&
1505 );
1506
1507// Intersects an infinte line and sphere and returns
1508// 0 = no intersections,
1509// 1 = one intersection,
1510// 2 = 2 intersections
1511// If 0 is returned, first point is point
1512// on line closest to sphere and 2nd point is the point
1513// on the sphere closest to the line.
1514// If 1 is returned, first point is obtained by evaluating
1515// the line and the second point is obtained by evaluating
1516// the sphere.
1517ON_DECL
1518int ON_Intersect(
1519 const ON_Line&,
1520 const ON_Sphere&,
1521 ON_3dPoint&,
1522 ON_3dPoint& // intersection point(s) returned here
1523 );
1524
1525
1526// Intersects an infinte line and cylinder and returns
1527// 0 = no intersections,
1528// 1 = one intersection,
1529// 2 = 2 intersections
1530// 3 = line lies on cylinder
1531//
1532// If 0 is returned, first point is point
1533// on line closest to cylinder and 2nd point is the point
1534// on the cylinder closest to the line.
1535// If 1 is returned, first point is obtained by evaluating
1536// the line and the second point is obtained by evaluating
1537// the cylinder.
1538//
1539// The value of cylinder.IsFinite() determines if the
1540// intersection is performed on the finite or infinite cylinder.
1541ON_DECL
1542int ON_Intersect(
1543 const ON_Line&, // [in]
1544 const ON_Cylinder&, // [in]
1545 ON_3dPoint&, // [out] first intersection point
1546 ON_3dPoint& // [out] second intersection point
1547 );
1548
1549// Description:
1550// Intersect an infinte line and circle.
1551// Parameters:
1552// line - [in]
1553// circle - [in]
1554// line_t0 - [out] line parameter of first intersection point
1555// circle_point0 - [out] first intersection point on circle
1556// line_t1 - [out] line parameter of second intersection point
1557// circle_point1 - [out] second intersection point on circle
1558// Returns:
1559// 0 No intersection
1560// 1 One intersection at line.PointAt(*line_t0)
1561// 2 Two intersections at line.PointAt(*line_t0)
1562// and line.PointAt(*line_t1).
1563ON_DECL
1564int ON_Intersect(
1565 const ON_Line& line,
1566 const ON_Circle& circle,
1567 double* line_t0,
1568 ON_3dPoint& circle_point0,
1569 double* line_t1,
1570 ON_3dPoint& circle_point1
1571 );
1572
1573
1574
1575// Description:
1576// Intersect a infinte line and arc.
1577// Parameters:
1578// line - [in]
1579// arc - [in]
1580// line_t0 - [out] line parameter of first intersection point
1581// arc_point0 - [out] first intersection point on arc
1582// line_t1 - [out] line parameter of second intersection point
1583// arc_point1 - [out] second intersection point on arc
1584// Returns:
1585// 0 No intersection
1586// 1 One intersection at line.PointAt(*line_t0)
1587// 2 Two intersections at line.PointAt(*line_t0)
1588// and line.PointAt(*line_t1).
1589ON_DECL
1590int ON_Intersect(
1591 const ON_Line& line,
1592 const ON_Arc& arc,
1593 double* line_t0,
1594 ON_3dPoint& arc_point0,
1595 double* line_t1,
1596 ON_3dPoint& arc_point1
1597 );
1598
1599// Description:
1600// Intersect a plane and a circle.
1601// Parameters:
1602// plane - [in]
1603// circle - [in]
1604// point0 - [out] first intersection point
1605// point1 - [out] second intersection point
1606// Returns:
1607// 0 No intersection
1608// 1 One intersection at point0
1609// 2 Two intersections at point0
1610// and point1.
1611// 3 Circle lies on plane
1612ON_DECL
1613int ON_Intersect(
1614 const ON_Plane& plane,
1615 const ON_Circle& circle,
1616 ON_3dPoint& point0,
1617 ON_3dPoint& point1
1618 );
1619
1620// Description:
1621// Intersect a plane and an arc.
1622// Parameters:
1623// plane - [in]
1624// arc - [in]
1625// point0 - [out] first intersection point
1626// point1 - [out] second intersection point
1627// Returns:
1628// 0 No intersection
1629// 1 One intersection at point0
1630// 2 Two intersections at point0
1631// and point1.
1632// 3 Arc lies on plane
1633ON_DECL
1634int ON_Intersect(
1635 const ON_Plane& plane,
1636 const ON_Arc& arc,
1637 ON_3dPoint& point0,
1638 ON_3dPoint& point1
1639 );
1640
1641
1642// returns 0 = no, 1 = yes, 2 = points are coincident and on line
1643ON_DECL
1644int ON_ArePointsOnLine(
1645 int, // dimension of points
1646 int, // is_rat = true if homogeneous rational
1647 int, // count = number of points
1648 int, // stride ( >= is_rat?(dim+1) :dim)
1649 const double*, // point array
1650 const ON_BoundingBox&, // if needed, use ON_GetBoundingBox(dim,is_rat,count,stride,point)
1651 const ON_Line&,
1652 double // tolerance (if 0.0, a tolerance based on bounding box size is used)
1653 );
1654
1655// returns 0 = no, 1 = yes, 2 = points are coincident and on line
1656ON_DECL
1657int ON_ArePointsOnPlane(
1658 int, // dimension of points
1659 int, // is_rat = true if homogeneous rational
1660 int, // count = number of points
1661 int, // stride ( >= is_rat?(dim+1) :dim)
1662 const double*, // point array
1663 const ON_BoundingBox&, // if needed, use ON_GetBoundingBox(dim,is_rat,count,stride,point)
1664 const ON_Plane&,
1665 double // tolerance (if 0.0, a tolerance based on bounding box size is used)
1666 );
1667
1668/*
1669Description:
1670 Use the quotient rule to compute derivatives of a one parameter
1671 rational function F(t) = X(t)/W(t), where W is a scalar
1672 and F and X are vectors of dimension dim.
1673Parameters:
1674 dim - [in]
1675 der_count - [in] number of derivative (>=0)
1676 v_stride - [in] (>= dim+1)
1677 v - [in/out]
1678 v[] is an array of length (der_count+1)*v_stride.
1679 The input v[] array contains derivatives of the numerator and
1680 denominator functions in the order (X, W), (Xt, Wt), (Xtt, Wtt), ...
1681 In general, the (dim+1) coordinates of the d-th derivative
1682 are in (v[n],...,v[n+dim]) where n = d*v_stride.
1683 In the output v[] array the derivatives of X are replaced with
1684 the derivatives of F and the derivatives of W are divided by
1685 w = v[dim].
1686Returns:
1687 True if input is valid; i.e., v[dim] != 0.
1688See Also:
1689 ON_EvaluateQuotientRule2
1690 ON_EvaluateQuotientRule3
1691*/
1692ON_DECL
1693bool ON_EvaluateQuotientRule(
1694 int dim,
1695 int der_count,
1696 int v_stride,
1697 double *v
1698 );
1699
1700/*
1701Description:
1702 Use the quotient rule to compute partial derivatives of a two parameter
1703 rational function F(s,t) = X(s,t)/W(s,t), where W is a scalar
1704 and F and X are vectors of dimension dim.
1705Parameters:
1706 dim - [in]
1707 der_count - [in] number of derivative (>=0)
1708 v_stride - [in] (>= dim+1)
1709 v - [in/out]
1710 v[] is an array of length (der_count+2)*(der_count+1)*v_stride.
1711 The input array contains derivatives of the numerator and denominator
1712 functions in the order X, W, Xs, Ws, Xt, Wt, Xss, Wss, Xst, Wst, Xtt, Wtt, ...
1713 In general, the (i,j)-th derivatives are in the (dim+1) entries of v[]
1714 v[k], ..., answer[k+dim], where k = ((i+j)*(i+j+1)/2 + j)*v_stride.
1715 In the output v[] array the derivatives of X are replaced with
1716 the derivatives of F and the derivatives of W are divided by
1717 w = v[dim].
1718Returns:
1719 True if input is valid; i.e., v[dim] != 0.
1720See Also:
1721 ON_EvaluateQuotientRule
1722 ON_EvaluateQuotientRule3
1723*/
1724ON_DECL
1725bool ON_EvaluateQuotientRule2(
1726 int dim,
1727 int der_count,
1728 int v_stride,
1729 double *v
1730 );
1731
1732/*
1733Description:
1734 Use the quotient rule to compute partial derivatives of a 3 parameter
1735 rational function F(r,s,t) = X(r,s,t)/W(r,s,t), where W is a scalar
1736 and F and X are vectors of dimension dim.
1737Parameters:
1738 dim - [in]
1739 der_count - [in] number of derivative (>=0)
1740 v_stride - [in] (>= dim+1)
1741 v - [in/out]
1742 v[] is an array of length
1743 v_stride*(der_count+1)*(der_count+2)*(der_count+3)/6.
1744 The input v[] array contains derivatives of the numerator and
1745 denominator functions in the order (X, W), (Xr, Wr), (Xs, Ws),
1746 (Xt, Wt), (Xrr, Wrr), (Xrs, Wrs), (Xrt, Wrt), (Xss, Wss),
1747 (Xst, Wst), (Xtt, Wtt), ...
1748 In general, the (dim+1) coordinates of the derivative
1749 (Dr^i Ds^j Dt^k, i+j+k=d) are at v[n], ..., v[n+dim] where
1750 n = v_stride*( d*(d+1)*(d+2)/6 + (d-i)*(d-i+1)/2 + k ).
1751 In the output v[] array the derivatives of X are replaced with
1752 the derivatives of F and the derivatives of W are divided by
1753 w = v[dim].
1754Returns:
1755 True if input is valid; i.e., v[dim] != 0.
1756See Also:
1757 ON_EvaluateQuotientRule
1758 ON_EvaluateQuotientRule2
1759*/
1760ON_DECL
1761bool ON_EvaluateQuotientRule3(
1762 int dim,
1763 int der_count,
1764 int v_stride,
1765 double *v
1766 );
1767
1768ON_DECL
1769bool ON_GetPolylineLength(
1770 int, // dimension of points
1771 ON_BOOL32, // bIsRational true if points are homogeneous rational
1772 int, // number of points
1773 int, // stride between points
1774 const double*, // points
1775 double* // length returned here
1776 );
1777
1778
1779/*
1780Description:
1781 Find the index of the point in the point_list that is closest to P.
1782Parameters:
1783 point_count - [in]
1784 point_list - [in]
1785 P - [in]
1786 closest_point_index - [out]
1787Returns:
1788 True if successful and *closest_point_index is set.
1789 False if input is not valid, in which case *closest_point_index
1790 is undefined.
1791*/
1792ON_DECL
1793bool ON_GetClosestPointInPointList(
1794 int point_count,
1795 const ON_3dPoint* point_list,
1796 ON_3dPoint P,
1797 int* closest_point_index
1798 );
1799
1800/*
1801Description:
1802 Test math library functions.
1803Parameters:
1804 function_index - [in] Determines which math library function is called.
1805
1806 1: z = x+y
1807 2: z = x-y
1808 3: z = x*y
1809 4: z = x/y
1810 5: z = fabs(x)
1811 6: z = exp(x)
1812 7: z = log(x)
1813 8: z = logb(x)
1814 9: z = log10(x)
1815 10: z = pow(x,y)
1816 11: z = sqrt(x)
1817 12: z = sin(x)
1818 13: z = cos(x)
1819 14: z = tan(x)
1820 15: z = sinh(x)
1821 16: z = cosh(x)
1822 17: z = tanh(x)
1823 18: z = asin(x)
1824 19: z = acos(x)
1825 20: z = atan(x)
1826 21: z = atan2(y,x)
1827 22: z = fmod(x,y)
1828 23: z = modf(x,&y)
1829 24: z = frexp(x,&y)
1830
1831 double x - [in]
1832 double y - [in]
1833Returns:
1834 Returns the "z" value listed in the function_index parameter
1835 description.
1836Remarks:
1837 This function is used to test the results of class floating
1838 point functions. It is primarily used to see what happens
1839 when opennurbs is used as a DLL and illegal operations are
1840 performed.
1841*/
1842ON_DECL
1843double ON_TestMathFunction(
1844 int function_index,
1845 double x,
1846 double y
1847 );
1848
1849// If performance is important, then
1850// you are better off using ((b<a)?a:b)
1851ON_DECL double ON_Max(double a, double b);
1852
1853// If performance is important, then
1854// you are better off using ((b<a)?a:b)
1855ON_DECL float ON_Max(float a, float b);
1856
1857// If performance is important, then
1858// you are better off using ((b<a)?a:b)
1859ON_DECL int ON_Max(int a, int b);
1860
1861// If performance is important, then
1862// you are better off using ((a<b)?a:b)
1863ON_DECL double ON_Min(double a, double b);
1864
1865// If performance is important, then
1866// you are better off using ((a<b)?a:b)
1867ON_DECL float ON_Min(float a, float b);
1868
1869// If performance is important, then
1870// you are better off using ((a<b)?a:b)
1871ON_DECL int ON_Min(int a, int b);
1872
1873// Do not call ON_Round() in any opennurbs code, tl code
1874// or any other code that does critical calculations or
1875// when there is any possibility that x is invalid or
1876// fabs(x)>2147483647. Use floor(x+0.5) instead.
1877ON_DECL int ON_Round(double x);
1878
1879
1880/*
1881Description:
1882 Find the equation of the parabola, ellipse or hyperbola
1883 (non-degenerate conic) that passes through six distinct points.
1884Parameters:
1885 stride - [in] (>=2)
1886 points array stride
1887 points2d - [in] (>=2)
1888 i-th point is (points[i*stride],points[i*stride+1])
1889 conic - [out]
1890 Coefficients of the conic equation.
1891 The points on the conic satisfy the equation
1892 0 = conic[0]*x^2 + conic[1]*xy + conic[2]*y^2
1893 + conic[3]*x + conic[4]*y + conic[5]
1894 max_pivot - [out] (can be null)
1895 min_pivot - [out] (can be null)
1896 zero_pivot - [out] (can be null)
1897 If there are some near duplicates in the input point set,
1898 the calculation is not stable. If you want to get an
1899 estimate of the validity of the solution, then inspect
1900 the returned values. max_pivot should around 1,
1901 min_pivot should be > 1e-4 or so, and zero_pivot should
1902 be < 1e-10 or so. If the returned pivots don't satisify
1903 these condtions, then exercise caution when using the
1904 returned solution.
1905Returns:
1906 True if a there is an ellipse, parabola or hyperbola through the
1907 six points.
1908 False if the input is invalid or the conic degenerate (the
1909 points lie on one or two lines).
1910 If false is returned, then conic[0]=...=conic[5] = 0 and
1911 *min_pivot = *max_pivot = *zero_pivot = 0.
1912*/
1913ON_DECL bool ON_GetConicEquationThrough6Points(
1914 int stride,
1915 const double* points2d,
1916 double conic[6],
1917 double* max_pivot,
1918 double* min_pivot,
1919 double* zero_pivot
1920 );
1921
1922/*
1923Description:
1924 Test a conic equation to see if it defines and ellipse. If so,
1925 return the center and axes of the ellipse.
1926Parameters:
1927 conic - [in]
1928 Coefficients of the conic equation.
1929 The points on the conic satisfy the equation
1930 0 = conic[0]*x^2 + conic[1]*xy + conic[2]*y^2
1931 + conic[3]*x + conic[4]*y + conic[5]
1932 center - [out]
1933 major_axis - [out]
1934 minor_axis - [out]
1935 major_radius - [out]
1936 minor_radius - [out]
1937Returns:
1938 True if the conic is an ellipse and the center and axes were found.
1939 False if the conic is not an ellipse, in which case the input values
1940 of center, major_axis, minor_axis, major_radius, and minor_radius
1941 are not changed.
1942*/
1943ON_DECL bool ON_IsConicEquationAnEllipse(
1944 const double conic[6],
1945 ON_2dPoint& center,
1946 ON_2dVector& major_axis,
1947 ON_2dVector& minor_axis,
1948 double* major_radius,
1949 double* minor_radius
1950 );
1951
1952/*
1953Description:
1954 Get the conic equation of an ellipse.
1955Parameters:
1956 a - [in] (a>0)
1957 b - [in] (b>0)
1958 a and b are the lengths of the axes. Either one
1959 may be largest and they can be equal.
1960 x0 - [in]
1961 y0 - [in]
1962 (x0,y0) is the enter of the ellipse.
1963 alpha - [in] (angle in radians)
1964 When alpha is 0, a corresponds to the x-axis and
1965 b corresponds to the y axis. If alpha is non-zero
1966 it specifies the rotation of these axes.
1967 conic - [out]
1968 Coefficients of the conic equation.
1969 The points on the conic satisfy the equation
1970 0 = conic[0]*x^2 + conic[1]*xy + conic[2]*y^2
1971 + conic[3]*x + conic[4]*y + conic[5]
1972 center - [out]
1973 major_axis - [out]
1974 minor_axis - [out]
1975 major_radius - [out]
1976 minor_radius - [out]
1977Remarks:
1978 Here is the way to evaluate a point on the ellipse:
1979
1980
1981 double t = ellipse paramter in radians;
1982 double x = a*cos(t);
1983 double y = b*sin(t);
1984 ON_2dPoint ellipse_point;
1985 ellipse_point.x = x0 + x*cos(alpha) + y*sin(alpha);
1986 ellipse_point.y = y0 - x*sin(alpha) + y*cos(alpha);
1987
1988Returns:
1989 True if the input is valid and conic[] was filled in.
1990 Falis if the input is not valid. In this case the values in conic[]
1991 are not changed.
1992*/
1993ON_DECL bool ON_GetEllipseConicEquation(
1994 double a, double b,
1995 double x0, double y0,
1996 double alpha,
1997 double conic[6]
1998 );
1999
2000/*
2001Descripton:
2002 Return the length of a 2d vector (x,y)
2003Returns:
2004 sqrt(x^2 + y^2) calculated in as precisely and safely as possible.
2005*/
2006ON_DECL double ON_Length2d( double x, double y );
2007
2008/*
2009Descripton:
2010 Return the length of a 3d vector (x,y,z)
2011Returns:
2012 sqrt(x^2 + y^2 + z^2) calculated in as precisely and safely as possible.
2013*/
2014ON_DECL double ON_Length3d( double x, double y, double z );
2015
2016
2017/*
2018Description:
2019 Convert a double x to the largest float f such that
2020 the mathematical value of f is at most the value of x.
2021Parameters:
2022 x - [in]
2023Returns
2024 The largest float f such that the mathematical value
2025 of f is at most the value of x.
2026*/
2027ON_DECL float ON_FloatFloor(double x);
2028
2029/*
2030Description:
2031 Convert a double x to the smallest float f such that
2032 the mathematical value of f is at least the value of x.
2033Parameters:
2034 x - [in]
2035Returns
2036 The smallest float f such that the mathematical value
2037 of f is at least the value of x.
2038*/
2039ON_DECL float ON_FloatCeil(double x);
2040
2041#endif
bool Periodic(int parameter_index) const
ON_SimpleArray< bool > m_bPeriodicParameter
virtual int EvaluateHessian(const double *parameters, double *value, double *gradient, double **hessian)
virtual int Evaluate(const double *parameters, double *values, double **jacobian)=0
const int m_value_count
ON_Interval Domain(int parameter_index) const
ON_SimpleArray< ON_Interval > m_domain
const int m_parameter_count
ON_Evaluator(int parameter_count, int value_count, const ON_Interval *domain, const bool *periodic)
virtual ~ON_Evaluator()
bool FiniteDomain() const
void operator=(double x)
void operator-=(double x)
void Plus(double x)
int SummandCount() const
void operator+=(double x)
void Begin(double starting_value=0.0)
double Total(double *error_estimate=NULL)