1 /*=========================================================================
2
3 Program: Visualization Toolkit
4 Module: UnitTestMath.cxx
5
6 Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7 All rights reserved.
8 See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9
10 This software is distributed WITHOUT ANY WARRANTY; without even
11 the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12 PURPOSE. See the above copyright notice for more information.
13
14 =========================================================================*/
15
16 #include "vtkMath.h"
17 #include "vtkMathUtilities.h"
18 #include "vtkSmartPointer.h"
19 #include "vtkUnsignedCharArray.h"
20 #include "vtkUnsignedShortArray.h"
21 #include "vtkType.h"
22
23 #include <vector>
24
25 static int TestPi();
26 static int TestDegreesFromRadians();
27 #ifndef VTK_LEGACY_REMOVE
28 static int TestRound();
29 #endif
30 static int TestFloor();
31 static int TestCeil();
32 static int TestCeilLog2();
33 static int TestIsPowerOfTwo();
34 static int TestNearestPowerOfTwo();
35 static int TestFactorial();
36 static int TestBinomial();
37 static int TestRandom();
38 static int TestAddSubtract();
39 static int TestMultiplyScalar();
40 static int TestMultiplyScalar2D();
41 static int TestDot();
42 static int TestOuter();
43 static int TestCross();
44 static int TestNorm();
45 static int TestNormalize();
46 static int TestPerpendiculars();
47 static int TestProjectVector();
48 static int TestProjectVector2D();
49 static int TestDistance2BetweenPoints();
50 static int TestAngleBetweenVectors();
51 static int TestGaussianAmplitude();
52 static int TestGaussianWeight();
53 static int TestDot2D();
54 static int TestNorm2D();
55 static int TestNormalize2D();
56 static int TestDeterminant2x2();
57 static int TestDeterminant3x3();
58 static int TestLUFactor3x3();
59 static int TestLUSolve3x3();
60 static int TestLinearSolve3x3();
61 static int TestMultiply3x3();
62 static int TestMultiplyMatrix();
63 static int TestTranspose3x3();
64 static int TestInvert3x3();
65 static int TestInvertMatrix();
66 static int TestIdentity3x3();
67 static int TestQuaternionToMatrix3x3();
68 static int TestMatrix3x3ToQuaternion();
69 static int TestMultiplyQuaternion();
70 static int TestOrthogonalize3x3();
71 static int TestDiagonalize3x3();
72 static int TestSingularValueDecomposition3x3();
73 static int TestSolveLinearSystem();
74 static int TestSolveLeastSquares();
75 static int TestSolveHomogeneousLeastSquares();
76 static int TestLUSolveLinearSystemEstimateMatrixCondition();
77 static int TestJacobiN();
78 static int TestClampValue();
79 static int TestClampValues();
80 static int TestClampAndNormalizeValue();
81 static int TestTensorFromSymmetricTensor();
82 static int TestGetScalarTypeFittingRange();
83 static int TestGetAdjustedScalarRange();
84 static int TestExtentIsWithinOtherExtent();
85 static int TestBoundsIsWithinOtherBounds();
86 static int TestPointIsWithinBounds();
87 static int TestSolve3PointCircle();
88 static int TestRGBToHSV();
89 static int TestInf();
90 static int TestNegInf();
91 static int TestNan();
92
UnitTestMath(int,char * [])93 int UnitTestMath(int,char *[])
94 {
95 int status = 0;
96
97 status += TestPi();
98
99 status += TestDegreesFromRadians();
100 #ifndef VTK_LEGACY_REMOVE
101 status += TestRound();
102 #endif
103 status += TestFloor();
104 status += TestCeil();
105 status += TestCeilLog2();
106 status += TestIsPowerOfTwo();
107 status += TestNearestPowerOfTwo();
108 status += TestFactorial();
109 status += TestBinomial();
110 status += TestRandom();
111 status += TestAddSubtract();
112 status += TestMultiplyScalar();
113 status += TestMultiplyScalar2D();
114 status += TestDot();
115 status += TestOuter();
116 status += TestCross();
117 status += TestNorm();
118 status += TestNormalize();
119 status += TestPerpendiculars();
120 status += TestProjectVector();
121 status += TestProjectVector2D();
122 status += TestDistance2BetweenPoints();
123 status += TestAngleBetweenVectors();
124 status += TestGaussianAmplitude();
125 status += TestGaussianWeight();
126 status += TestDot2D();
127 status += TestNorm2D();
128 status += TestNormalize2D();
129 status += TestDeterminant2x2();
130 status += TestDeterminant3x3();
131 status += TestLUFactor3x3();
132 status += TestLUSolve3x3();
133 status += TestLinearSolve3x3();
134 status += TestMultiply3x3();
135 status += TestMultiplyMatrix();
136 status += TestTranspose3x3();
137 status += TestInvert3x3();
138 status += TestInvertMatrix();
139 status += TestIdentity3x3();
140 status += TestQuaternionToMatrix3x3();
141 status += TestMatrix3x3ToQuaternion();
142 status += TestMultiplyQuaternion();
143 status += TestOrthogonalize3x3();
144 status += TestDiagonalize3x3();
145 status += TestSingularValueDecomposition3x3();
146 status += TestSolveLinearSystem();
147 status += TestSolveLeastSquares();
148 status += TestSolveHomogeneousLeastSquares();
149 status += TestLUSolveLinearSystemEstimateMatrixCondition();
150 status += TestJacobiN();
151 status += TestClampValue();
152 status += TestClampValues();
153 status += TestClampAndNormalizeValue();
154 status += TestTensorFromSymmetricTensor();
155 status += TestGetScalarTypeFittingRange();
156 status += TestGetAdjustedScalarRange();
157 status += TestExtentIsWithinOtherExtent();
158 status += TestBoundsIsWithinOtherBounds();
159 status += TestPointIsWithinBounds();
160 status += TestSolve3PointCircle();
161 status += TestRGBToHSV();
162 status += TestInf();
163 status += TestNegInf();
164 status += TestNan();
165 if (status != 0)
166 {
167 return EXIT_FAILURE;
168 }
169
170 vtkSmartPointer<vtkMath> math =
171 vtkSmartPointer<vtkMath>::New();
172 math->Print(std::cout);
173
174 return EXIT_SUCCESS;
175 }
176
177 // Validate by comparing to atan/4
TestPi()178 int TestPi()
179 {
180 int status = 0;
181 std::cout << "Pi..";
182 if (vtkMath::Pi() != std::atan(1.0) * 4.0)
183 {
184 std::cout << "Expected " << vtkMath::Pi()
185 << " but got " << std::atan(1.0) * 4.0;
186 ++status;
187 }
188
189 if (status)
190 {
191 std::cout << "..FAILED" << std::endl;
192 }
193 else
194 {
195 std::cout << ".PASSED" << std::endl;
196 }
197 return status;
198 }
199
200 // Validate against RadiansFromDegress
TestDegreesFromRadians()201 int TestDegreesFromRadians()
202 {
203 int status = 0;
204 std::cout << "DegreesFromRadians..";
205
206 unsigned int numSamples = 1000;
207 for (unsigned int i = 0; i < numSamples; ++i)
208 {
209 float floatDegrees = vtkMath::Random(-180.0, 180.0);
210 float floatRadians = vtkMath::RadiansFromDegrees(floatDegrees);
211 float result = vtkMath::DegreesFromRadians(floatRadians);
212 if (!vtkMathUtilities::FuzzyCompare(
213 result, floatDegrees,
214 std::numeric_limits<float>::epsilon()*128.0f))
215 {
216 std::cout << "Float Expected " << floatDegrees
217 << " but got " << result
218 << " difference is " << result - floatDegrees << " ";
219 std::cout << "eps ratio is: " << (result - floatDegrees)
220 / std::numeric_limits<float>::epsilon() << std::endl;
221 ++status;
222 }
223 }
224 for (unsigned int i = 0; i < numSamples; ++i)
225 {
226 double doubleDegrees = vtkMath::Random(-180.0, 180.0);
227 double doubleRadians = vtkMath::RadiansFromDegrees(doubleDegrees);
228 double result = vtkMath::DegreesFromRadians(doubleRadians);
229 if (!vtkMathUtilities::FuzzyCompare(
230 result, doubleDegrees,
231 std::numeric_limits<double>::epsilon()*256.0))
232 {
233 std::cout << " Double Expected " << doubleDegrees
234 << " but got " << result
235 << " difference is " << result - doubleDegrees;
236 std::cout << " eps ratio is: " << (result - doubleDegrees)
237 / std::numeric_limits<double>::epsilon() << std::endl;
238 ++status;
239 }
240 }
241 if (status)
242 {
243 std::cout << "..FAILED" << std::endl;
244 }
245 else
246 {
247 std::cout << ".PASSED" << std::endl;
248 }
249 return status;
250 }
251
252 #ifndef VTK_LEGACY_REMOVE
253 // Validate with http://en.wikipedia.org/wiki/Rounding#Rounding_to_integer
TestRound()254 int TestRound()
255 {
256 int status = 0;
257 std::cout << "Round..";
258 int result;
259 {
260 std::vector<float> values;
261 std::vector<int> expecteds;
262
263 values.push_back(23.67f); expecteds.push_back(24);
264 values.push_back(23.50f); expecteds.push_back(24);
265 values.push_back(23.35f); expecteds.push_back(23);
266 values.push_back(23.00f); expecteds.push_back(23);
267 values.push_back(0.00f); expecteds.push_back( 0);
268 values.push_back(-23.00f); expecteds.push_back(-23);
269 values.push_back(-23.35f); expecteds.push_back(-23);
270 values.push_back(-23.50f); expecteds.push_back(-24);
271 values.push_back(-23.67f); expecteds.push_back(-24);
272 for ( size_t i = 0; i < values.size(); ++i)
273 {
274 result = vtkMath::Round(values[i]);
275 if (result != expecteds[i])
276 {
277 std::cout << " Float Round(" << values[i] << ") got " << result
278 << " but expected " << expecteds[i];
279 ++status;
280 }
281 }
282 }
283 {
284 std::vector<double> values;
285 std::vector<int> expecteds;
286
287 values.push_back(23.67); expecteds.push_back(24);
288 values.push_back(23.50); expecteds.push_back(24);
289 values.push_back(23.35); expecteds.push_back(23);
290 values.push_back(23.00); expecteds.push_back(23);
291 values.push_back(0.00); expecteds.push_back( 0);
292 values.push_back(-23.00); expecteds.push_back(-23);
293 values.push_back(-23.35); expecteds.push_back(-23);
294 values.push_back(-23.50); expecteds.push_back(-24);
295 values.push_back(-23.67); expecteds.push_back(-24);
296 for ( size_t i = 0; i < values.size(); ++i)
297 {
298 result = vtkMath::Round(values[i]);
299 if (result != expecteds[i])
300 {
301 std::cout << " Double Round(" << values[i] << ") got " << result
302 << " but expected " << expecteds[i];
303 ++status;
304 }
305 }
306 }
307 if (status)
308 {
309 std::cout << "..FAILED" << std::endl;
310 }
311 else
312 {
313 std::cout << ".PASSED" << std::endl;
314 }
315 return status;
316 }
317 #endif
318
319 // Validate with http://en.wikipedia.org/wiki/Floor_and_ceiling_functions
TestFloor()320 int TestFloor()
321 {
322 int status = 0;
323 std::cout << "Floor..";
324
325 int result;
326 std::vector<double> values;
327 std::vector<int> expecteds;
328
329 values.push_back(2.4); expecteds.push_back(2);
330 values.push_back(2.7); expecteds.push_back(2);
331 values.push_back(-2.7); expecteds.push_back(-3);
332 values.push_back(-2.0); expecteds.push_back(-2);
333 for ( size_t i = 0; i < values.size(); ++i)
334 {
335 result = vtkMath::Floor(values[i]);
336 if (result != expecteds[i])
337 {
338 std::cout << " Floor(" << values[i] << ") got " << result
339 << " but expected " << expecteds[i];
340 ++status;
341 }
342 }
343 if (status)
344 {
345 std::cout << "..FAILED" << std::endl;
346 }
347 else
348 {
349 std::cout << ".PASSED" << std::endl;
350 }
351 return status;
352 }
353
354 // Validate with http://en.wikipedia.org/wiki/Floor_and_ceiling_functions
TestCeil()355 int TestCeil()
356 {
357 int status = 0;
358 std::cout << "Ceil..";
359
360 int result;
361 std::vector<double> values;
362 std::vector<int> expecteds;
363
364 values.push_back(2.4); expecteds.push_back(3);
365 values.push_back(2.7); expecteds.push_back(3);
366 values.push_back(-2.7); expecteds.push_back(-2);
367 values.push_back(-2.0); expecteds.push_back(-2);
368 for ( size_t i = 0; i < values.size(); ++i)
369 {
370 result = vtkMath::Ceil(values[i]);
371 if (result != expecteds[i])
372 {
373 std::cout << " Ceil(" << values[i] << ") got " << result
374 << " but expected " << expecteds[i];
375 ++status;
376 }
377 }
378 if (status)
379 {
380 std::cout << "..FAILED" << std::endl;
381 }
382 else
383 {
384 std::cout << ".PASSED" << std::endl;
385 }
386 return status;
387 }
388
389 // Validate by powers of 2 perturbations
TestCeilLog2()390 int TestCeilLog2()
391 {
392 int status = 0;
393 std::cout << "CeilLog2..";
394
395 int result;
396 std::vector<vtkTypeUInt64> values;
397 std::vector<int> expecteds;
398
399 for ( unsigned int p = 0; p < 30; ++p)
400 {
401 vtkTypeUInt64 shifted = (2 << p) + 1;
402 values.push_back(shifted); expecteds.push_back(p + 2);
403 shifted = (2 << p);
404 values.push_back(shifted); expecteds.push_back(p + 1);
405 }
406 for ( size_t i = 0; i < values.size(); ++i)
407 {
408 result = vtkMath::CeilLog2(values[i]);
409 if (result != expecteds[i])
410 {
411 std::cout << " CeilLog2(" << values[i] << ") got " << result
412 << " but expected " << expecteds[i];
413 ++status;
414 }
415 }
416
417 if (status)
418 {
419 std::cout << "..FAILED" << std::endl;
420 }
421 else
422 {
423 std::cout << ".PASSED" << std::endl;
424 }
425 return status;
426 }
427
428 // Validate by powers of 2 perturbations
TestIsPowerOfTwo()429 int TestIsPowerOfTwo()
430 {
431 int status = 0;
432 std::cout << "IsPowerOfTwo..";
433 bool result;
434
435 std::vector<vtkTypeUInt64> values;
436 std::vector<bool> expecteds;
437 int largestPower = std::numeric_limits<vtkTypeUInt64>::digits;
438 vtkTypeUInt64 shifted = 1;
439 for ( int p = 1; p < largestPower - 1; ++p)
440 {
441 shifted *= 2;
442 values.push_back(shifted); expecteds.push_back(true);
443 if (shifted != 2)
444 {
445 values.push_back(shifted - 1); expecteds.push_back(false);
446 }
447 if (shifted < std::numeric_limits<vtkTypeUInt64>::max() - 1)
448 {
449 values.push_back(shifted + 1); expecteds.push_back(false);
450 }
451 }
452 for ( size_t i = 0; i < values.size(); ++i)
453 {
454 result = vtkMath::IsPowerOfTwo(values[i]);
455 if (result != expecteds[i])
456 {
457 std::cout << " IsPowerOfTwo(" << values[i] << ") got " << result
458 << " but expected " << expecteds[i];
459 ++status;
460 }
461 }
462
463 if (status)
464 {
465 std::cout << "..FAILED" << std::endl;
466 }
467 else
468 {
469 std::cout << ".PASSED" << std::endl;
470 }
471 return status;
472 }
473
474 // Validate by powers of 2 perturbations
TestNearestPowerOfTwo()475 int TestNearestPowerOfTwo()
476 {
477 int status = 0;
478 std::cout << "NearestPowerOfTwo..";
479
480 std::vector<vtkTypeUInt64> values;
481 std::vector<int> expecteds;
482
483 values.push_back(0);
484 expecteds.push_back(1);
485
486 int numDigits = std::numeric_limits<int>::digits;
487 vtkTypeUInt64 shifted = 1;
488 for ( int p = 0; p < numDigits; ++p)
489 {
490 values.push_back(shifted);
491 expecteds.push_back(shifted);
492 if (shifted <= INT_MAX/2 )
493 {
494 values.push_back(shifted + 1);
495 expecteds.push_back(shifted * 2);
496 }
497 if (shifted != 2 )
498 {
499 values.push_back(shifted - 1);
500 expecteds.push_back(shifted);
501 }
502
503 shifted *= 2;
504 }
505
506 values.push_back(INT_MAX);
507 expecteds.push_back(INT_MIN);
508
509 for ( size_t i = 0; i < values.size(); ++i)
510 {
511 int result = vtkMath::NearestPowerOfTwo(values[i]);
512 if (result != expecteds[i])
513 {
514 std::cout << " NearestPowerOfTwo(" << values[i] << ") got " << result
515 << " but expected " << expecteds[i];
516 ++status;
517 }
518 }
519
520 if (status)
521 {
522 std::cout << "..FAILED" << std::endl;
523 }
524 else
525 {
526 std::cout << ".PASSED" << std::endl;
527 }
528 return status;
529 }
530
531 // Validate by alternate computation
TestFactorial()532 int TestFactorial()
533 {
534 int status = 0;
535 std::cout << "Factorial..";
536
537 std::vector<int> values;
538 std::vector<vtkTypeInt64> expecteds;
539 vtkTypeInt64 expected = 1;
540 for ( int f = 2; f < 10; ++f)
541 {
542 expected *= f;
543 values.push_back(f); expecteds.push_back(expected);
544 }
545 for ( size_t i = 0; i < values.size(); ++i)
546 {
547 int result = vtkMath::Factorial(values[i]);
548 if (result != expecteds[i])
549 {
550 std::cout << " Factorial(" << values[i] << ") got " << result
551 << " but expected " << expecteds[i];
552 ++status;
553 }
554 }
555
556 if (status)
557 {
558 std::cout << "..FAILED" << std::endl;
559 }
560 else
561 {
562 std::cout << ".PASSED" << std::endl;
563 }
564 return status;
565 }
566
567 // Validate by alternative computations
TestBinomial()568 int TestBinomial()
569 {
570 int status = 0;
571 int m;
572 int n;
573 std::cout << "Binomial..";
574
575 std::vector<int> mvalues;
576 std::vector<int> nvalues;
577
578 std::vector<vtkTypeInt64> expecteds;
579 double expected;
580 for (m = 1; m < 31; ++m)
581 {
582 for (n = 1; n <= m; ++n)
583 {
584 mvalues.push_back(m);
585 nvalues.push_back(n);
586 expected = 1;
587 for (int i = 1; i <= n; ++i)
588 {
589 expected *= static_cast<double>(m - i + 1) / i;
590 }
591 expecteds.push_back(static_cast<vtkTypeInt64>(expected));
592 }
593 }
594
595 for ( size_t i = 0; i < mvalues.size(); ++i)
596 {
597 int result = vtkMath::Binomial(mvalues[i], nvalues[i]);
598 if (result != expecteds[i])
599 {
600 std::cout << " Binomial("
601 << mvalues[i] << ", " << nvalues[i] << ") got " << result
602 << " but expected " << expecteds[i];
603 ++status;
604 }
605 }
606
607 // Now test the combination iterator
608 m = 6;
609 n = 3;
610 int more = 1;
611 int count = 0;
612 int *comb;
613 // First, m < n should produce 0
614 comb = vtkMath::BeginCombination(n, m);
615 if (comb != nullptr)
616 {
617 ++status;
618 std::cout << " Combinations("
619 << n << ", " << m << ") should return 0 "
620 << " but got " << comb;
621 }
622 comb = vtkMath::BeginCombination(m, n);
623 while (more)
624 {
625 ++count;
626 more = vtkMath::NextCombination (m, n, comb);
627 }
628 vtkMath::FreeCombination(comb);
629 if (count != vtkMath::Binomial(m, n))
630 {
631 ++status;
632 std::cout << " Combinations("
633 << m << ", " << n << ") got " << count
634 << " but expected " << vtkMath::Binomial(m, n);
635 }
636 if (status)
637 {
638 std::cout << "..FAILED" << std::endl;
639 }
640 else
641 {
642 std::cout << ".PASSED" << std::endl;
643 }
644 return status;
645 }
646
647 // No validation
TestRandom()648 int TestRandom()
649 {
650 int status = 0;
651 std::cout << "Random..";
652 // Not really a test of randomness, just covering code
653 int n = 1000;
654 vtkMath::RandomSeed(8775070);
655 vtkMath::GetSeed(); // just for coverage
656 double accum = 0.0;
657 for (int i = 0; i < n; ++i)
658 {
659 float random = vtkMath::Random();
660 accum += random;
661 if ( random < 0.0 || random > 1.0)
662 {
663 std::cout << "Random(): " << random << " out of range" << std::endl;
664 ++status;
665 }
666 random = vtkMath::Gaussian ();
667 accum += random;
668
669 random = vtkMath::Gaussian (0.0, 1.0);
670 accum += random;
671
672 random = vtkMath::Random (-1000.0, 1000.0);
673 accum += random;
674 if (random < -1000.0 || random > 1000.0)
675 {
676 std::cout << "Random (-1000.0, 1000.0): " << random << " out of range" << std::endl;
677 ++status;
678 }
679 }
680 if (accum == 0.0)
681 {
682 ++status;
683 }
684 if (status)
685 {
686 std::cout << "..FAILED" << std::endl;
687 }
688 else
689 {
690 std::cout << ".PASSED" << std::endl;
691 }
692 return status;
693 }
694
695 template<typename T>
AddSubtract()696 int AddSubtract()
697 {
698 int status = 0;
699 T da[3], db[3], dc[3], dd[3];
700 for (int n = 0; n < 100000; ++n)
701 {
702 for (int i = 0; i < 3; ++i)
703 {
704 da[i] = vtkMath::Random(-10.0, 10.0);
705 db[i] = vtkMath::Random(-10.0, 10.0);
706 }
707 vtkMath::Add(da, db, dc);
708 vtkMath::Subtract(dc, db, dd);
709 for (int i = 0; i < 3; ++i)
710 {
711 if (!vtkMathUtilities::FuzzyCompare(
712 da[i], dd[i],
713 std::numeric_limits<T>::epsilon() * (T) 256.0))
714 {
715 std::cout << " Add/Subtract got " << dd[i]
716 << " but expected " << da[i];
717 }
718 }
719 }
720 return status;
721 }
722
723 // Validate by a + b - b = a
TestAddSubtract()724 int TestAddSubtract()
725 {
726 int status = 0;
727 std::cout << "AddSubtract..";
728
729 status += AddSubtract<double>();
730 status += AddSubtract<float>();
731
732 if (status)
733 {
734 std::cout << "..FAILED" << std::endl;
735 }
736 else
737 {
738 std::cout << ".PASSED" << std::endl;
739 }
740 return status;
741 }
742
743 template<typename T>
MultiplyScalar()744 int MultiplyScalar()
745 {
746 int status = 0;
747 // first T
748 T da[3], db[3];
749 for (int n = 0; n < 100000; ++n)
750 {
751 for (int i = 0; i < 3; ++i)
752 {
753 da[i] = vtkMath::Random(-10.0, 10.0);
754 db[i] = da[i];
755 }
756 T scale = vtkMath::Random();
757 vtkMath::MultiplyScalar(da, scale);
758
759 for (int i = 0; i < 3; ++i)
760 {
761 if (!vtkMathUtilities::FuzzyCompare(
762 da[i], db[i] * scale,
763 std::numeric_limits<T>::epsilon() * (T) 256.0))
764 {
765 std::cout << " MultiplyScalar got " << da[i]
766 << " but expected " << db[i] * scale;
767 }
768 }
769 }
770
771 return status;
772 }
773
TestMultiplyScalar()774 int TestMultiplyScalar()
775 {
776 int status = 0;
777 std::cout << "MultiplyScalar..";
778
779 status += MultiplyScalar<double>();
780 status += MultiplyScalar<float>();
781
782 if (status)
783 {
784 std::cout << "..FAILED" << std::endl;
785 }
786 else
787 {
788 std::cout << ".PASSED" << std::endl;
789 }
790 return status;
791 }
792
TestMultiplyScalar2D()793 int TestMultiplyScalar2D()
794 {
795 int status = 0;
796 std::cout << "MultiplyScalar2D..";
797
798 // now 2D
799 // first double
800 double da[2], db[2];
801 for (int n = 0; n < 100000; ++n)
802 {
803 for (int i = 0; i < 2; ++i)
804 {
805 da[i] = vtkMath::Random(-10.0, 10.0);
806 db[i] = da[i];
807 }
808 double scale = vtkMath::Random();
809 vtkMath::MultiplyScalar2D(da, scale);
810
811 for (int i = 0; i < 2; ++i)
812 {
813 if (!vtkMathUtilities::FuzzyCompare(
814 da[i], db[i] * scale,
815 std::numeric_limits<double>::epsilon() * 256.0))
816 {
817 std::cout << " MultiplyScalar2D got " << da[i]
818 << " but expected " << db[i] * scale;
819 }
820 }
821 }
822
823 // then float
824 float fa[2], fb[2];
825 for (int n = 0; n < 100000; ++n)
826 {
827 for (int i = 0; i < 2; ++i)
828 {
829 fa[i] = vtkMath::Random(-10.0, 10.0);
830 fb[i] = fa[i];
831 }
832 float scale = vtkMath::Random();
833 vtkMath::MultiplyScalar2D(fa, scale);
834
835 for (int i = 0; i < 2; ++i)
836 {
837 if (!vtkMathUtilities::FuzzyCompare(
838 fa[i], fb[i] * scale,
839 std::numeric_limits<float>::epsilon() * 256.0f))
840 {
841 std::cout << " MultiplyScalar2D got " << fa[i]
842 << " but expected " << fb[i] * scale;
843 }
844 }
845 }
846
847 if (status)
848 {
849 std::cout << "..FAILED" << std::endl;
850 }
851 else
852 {
853 std::cout << ".PASSED" << std::endl;
854 }
855 return status;
856 }
857
858 class valueDouble3D
859 {
860 public:
861 valueDouble3D() = default;
valueDouble3D(double aa[3],double bb[3])862 valueDouble3D(double aa[3], double bb[3])
863 {
864 for (int i = 0; i < 3; ++i)
865 {
866 a[i] = aa[i];
867 b[i] = bb[i];
868 }
869 }
870 double a[3];
871 double b[3];
872 };
873
874 class valueFloat3D
875 {
876 public:
877 valueFloat3D() = default;
valueFloat3D(float aa[3],float bb[3])878 valueFloat3D(float aa[3], float bb[3])
879 {
880 for (int i = 0; i < 3; ++i)
881 {
882 a[i] = aa[i];
883 b[i] = bb[i];
884 }
885 }
886 float a[3];
887 float b[3];
888 };
TestDot()889 int TestDot()
890 {
891 int status = 0;
892 std::cout << "Dot..";
893
894 {
895 std::vector<valueDouble3D> values;
896 std::vector<double> expecteds;
897 for (int n = 0; n < 100; ++n)
898 {
899 valueDouble3D v;
900 double dot = 0.0;
901 for (int i = 0; i < 3; ++i)
902 {
903 v.a[i] = vtkMath::Random();
904 v.b[i] = vtkMath::Random();
905 dot += (v.a[i] * v.b[i]);
906 }
907 values.push_back(v);
908 expecteds.push_back(dot);
909 }
910 valueDouble3D test;
911 test.a[0] = 0.0; test.a[1] = 0.0; test.a[2] = 1.0;
912 test.b[0] = 1.0; test.b[1] = 0.0; test.b[2] = 0.0;
913 values.push_back(test);
914 expecteds.push_back(0.0);
915 test.a[0] = 0.0; test.a[1] = 0.0; test.a[2] = 1.0;
916 test.b[0] = 0.0; test.b[1] = 1.0; test.b[2] = 0.0;
917 values.push_back(test);
918 expecteds.push_back(0.0);
919 test.a[0] = 1.0; test.a[1] = 0.0; test.a[2] = 0.0;
920 test.b[0] = 0.0; test.b[1] = 1.0; test.b[2] = 0.0;
921 values.push_back(test);
922 expecteds.push_back(0.0);
923
924 for ( size_t i = 0; i < values.size(); ++i)
925 {
926 double result = vtkMath::Dot(values[i].a, values[i].b);
927 if (!vtkMathUtilities::FuzzyCompare(
928 result, expecteds[i],
929 std::numeric_limits<double>::epsilon()*128.0))
930 {
931 std::cout << " Dot got " << result
932 << " but expected " << expecteds[i];
933 ++status;
934 }
935 }
936 }
937
938 // now float
939 {
940 std::vector<valueFloat3D> values;
941 std::vector<float> expecteds;
942 for (int n = 0; n < 100; ++n)
943 {
944 valueFloat3D v;
945 float dot = 0.0;
946 for (int i = 0; i < 3; ++i)
947 {
948 v.a[i] = vtkMath::Random();
949 v.b[i] = vtkMath::Random();
950 dot += (v.a[i] * v.b[i]);
951 }
952 values.push_back(v);
953 expecteds.push_back(dot);
954 }
955 valueFloat3D test;
956 test.a[0] = 0.0; test.a[1] = 0.0; test.a[2] = 1.0;
957 test.b[0] = 1.0; test.b[1] = 0.0; test.b[2] = 0.0;
958 values.push_back(test);
959 expecteds.push_back(0.0);
960 test.a[0] = 0.0; test.a[1] = 0.0; test.a[2] = 1.0;
961 test.b[0] = 0.0; test.b[1] = 1.0; test.b[2] = 0.0;
962 values.push_back(test);
963 expecteds.push_back(0.0);
964 test.a[0] = 1.0; test.a[1] = 0.0; test.a[2] = 0.0;
965 test.b[0] = 0.0; test.b[1] = 1.0; test.b[2] = 0.0;
966 values.push_back(test);
967 expecteds.push_back(0.0);
968
969 for ( size_t i = 0; i < values.size(); ++i)
970 {
971 float result = vtkMath::Dot(values[i].a, values[i].b);
972 if (!vtkMathUtilities::FuzzyCompare(
973 result, expecteds[i],
974 std::numeric_limits<float>::epsilon()*128.0f))
975 {
976 std::cout << " Dot got " << result
977 << " but expected " << expecteds[i];
978 ++status;
979 }
980 }
981 }
982
983 if (status)
984 {
985 std::cout << "..FAILED" << std::endl;
986 }
987 else
988 {
989 std::cout << ".PASSED" << std::endl;
990 }
991 return status;
992 }
993
TestOuter()994 int TestOuter()
995 {
996 int status = 0;
997 std::cout << "Outer..";
998
999 if (status)
1000 {
1001 std::cout << "..FAILED" << std::endl;
1002 }
1003 else
1004 {
1005 std::cout << ".PASSED" << std::endl;
1006 }
1007 return status;
1008 }
1009
1010 // Verify by anticommutative property
1011 template<typename T>
Cross()1012 int Cross()
1013 {
1014 int status = 0;
1015 T a[3];
1016 T b[3];
1017 T c[3];
1018 T d[3];
1019
1020 for (int n = 0; n < 1000; ++n)
1021 {
1022 for (int i = 0; i < 3; ++i)
1023 {
1024 a[i] = vtkMath::Random(-1.0, 1.0);
1025 b[i] = vtkMath::Random(-1.0, 1.0);
1026 }
1027 vtkMath::Cross(a, b, c);
1028 vtkMath::MultiplyScalar(b, (T) -1.0);
1029 vtkMath::Cross(b, a, d);
1030 // a x b = -b x a
1031 for (int i = 0; i < 3; ++i)
1032 {
1033 if (!vtkMathUtilities::FuzzyCompare(
1034 c[i], d[i],
1035 std::numeric_limits<T>::epsilon()*(T)128.0))
1036 {
1037 std::cout << " Cross expected " << c[i]
1038 << " but got " << d[i];
1039 std::cout << "eps ratio is: " << (c[i] - d[i])
1040 / std::numeric_limits<T>::epsilon() << std::endl;
1041 ++status;
1042 }
1043 }
1044 }
1045 return status;
1046 }
1047
TestCross()1048 int TestCross()
1049 {
1050 int status = 0;
1051 std::cout << "Cross..";
1052
1053 status += Cross<double>();
1054 status += Cross<float>();
1055
1056 if (status)
1057 {
1058 std::cout << "..FAILED" << std::endl;
1059 }
1060 else
1061 {
1062 std::cout << ".PASSED" << std::endl;
1063 }
1064 return status;
1065 }
1066
1067 template<typename T, int NDimension>
Norm()1068 int Norm()
1069 {
1070 int status = 0;
1071 T x[NDimension];
1072
1073 for (int n = 0; n < 1000; ++n)
1074 {
1075 for (int i = 0; i < NDimension; ++i)
1076 {
1077 x[i] = (T) vtkMath::Random(-10.0, 10.0);
1078 }
1079
1080 T norm = vtkMath::Norm(x, NDimension);
1081
1082 for (int i = 0; i < NDimension; ++i)
1083 {
1084 x[i] /= norm;
1085 }
1086
1087 T unitNorm = vtkMath::Norm(x, NDimension);
1088 if (!vtkMathUtilities::FuzzyCompare(
1089 unitNorm, (T) 1.0,
1090 std::numeric_limits<T>::epsilon()*(T)128.0))
1091 {
1092 std::cout << "Norm Expected " << 1.0
1093 << " but got " << unitNorm;
1094 std::cout << " eps ratio is: " << ((T) 1.0 - unitNorm)
1095 / std::numeric_limits<T>::epsilon() << std::endl;
1096 ++status;
1097 }
1098 }
1099
1100 return status;
1101 }
1102
TestNorm()1103 int TestNorm()
1104 {
1105 int status = 0;
1106 std::cout << "Norm..";
1107
1108 status += Norm<double,1>();
1109 status += Norm<double,3>();
1110 status += Norm<double,1000>();
1111 status += Norm<float,1>();
1112 status += Norm<float,3>();
1113 status += Norm<float,1000>();
1114
1115 if (status)
1116 {
1117 std::cout << "..FAILED" << std::endl;
1118 }
1119 else
1120 {
1121 std::cout << ".PASSED" << std::endl;
1122 }
1123 return status;
1124 }
1125
1126 // Validate compute Norm, should be 1.0
1127
1128 template<typename T>
Normalize()1129 int Normalize()
1130 {
1131 int status = 0;
1132 for (int n = 0; n < 1000; ++n)
1133 {
1134 T a[3];
1135 for (int i = 0; i < 3; ++i)
1136 {
1137 a[i] = vtkMath::Random(-10000.0,
1138 10000.0);
1139 }
1140 vtkMath::Normalize(a);
1141 T value= vtkMath::Norm(a);
1142 T expected = 1.0;
1143 if (!vtkMathUtilities::FuzzyCompare(
1144 value, expected,
1145 std::numeric_limits<T>::epsilon() * (T) 128.0))
1146 {
1147 std::cout << " Normalize expected " << expected
1148 << " but got " << value;
1149 std::cout << "eps ratio is: " << value - expected
1150 / std::numeric_limits<T>::epsilon() << std::endl;
1151 ++status;
1152 }
1153 }
1154 return status;
1155 }
1156
TestNormalize()1157 int TestNormalize()
1158 {
1159 int status = 0;
1160 std::cout << "Normalize..";
1161
1162 status += Normalize<double>();
1163 status += Normalize<float>();
1164
1165 if (status)
1166 {
1167 std::cout << "..FAILED" << std::endl;
1168 }
1169 else
1170 {
1171 std::cout << ".PASSED" << std::endl;
1172 }
1173 return status;
1174 }
1175
TestPerpendiculars()1176 int TestPerpendiculars()
1177 {
1178 int status = 0;
1179 std::cout << "Perpendiculars..";
1180 {
1181 // first double
1182 double x[3], y[3], z[3];
1183 std::vector<valueDouble3D> values;
1184 std::vector<double> expecteds;
1185 for (int n = 0; n < 100; ++n)
1186 {
1187 for (int i = 0; i < 3; ++i)
1188 {
1189 x[i] = vtkMath::Random(-10.0, 10.0);
1190 }
1191 vtkMath::Perpendiculars(x, y, z,
1192 vtkMath::Random(-vtkMath::Pi(), vtkMath::Pi()));
1193 {
1194 valueDouble3D value(x, y);
1195 values.push_back(value);
1196 expecteds.push_back(0.0);
1197 }
1198 {
1199 valueDouble3D value(x, z);
1200 values.push_back(value);
1201 expecteds.push_back(0.0);
1202 }
1203 {
1204 valueDouble3D value(y, z);
1205 values.push_back(value);
1206 expecteds.push_back(0.0);
1207 }
1208 vtkMath::Perpendiculars(x, y, z,
1209 0.0);
1210 {
1211 valueDouble3D value(x, y);
1212 values.push_back(value);
1213 expecteds.push_back(0.0);
1214 }
1215 }
1216 for (size_t i = 0; i < values.size(); ++i)
1217 {
1218 double test = vtkMath::Dot(values[i].a, values[i].b);
1219 if (!vtkMathUtilities::FuzzyCompare(
1220 expecteds[i], test,
1221 std::numeric_limits<double>::epsilon() * 256.0))
1222 {
1223 std::cout << " Perpendiculars got " << test
1224 << " but expected " << expecteds[i];
1225 }
1226 }
1227 }
1228 {
1229 // then floats
1230 float x[3], y[3], z[3];
1231 std::vector<valueFloat3D> values;
1232 std::vector<float> expecteds;
1233 for (int n = 0; n < 100; ++n)
1234 {
1235 for (int i = 0; i < 3; ++i)
1236 {
1237 x[i] = vtkMath::Random(-10.0, 10.0);
1238 }
1239 vtkMath::Perpendiculars(x, y, z,
1240 vtkMath::Random(-vtkMath::Pi(), vtkMath::Pi()));
1241 {
1242 valueFloat3D value(x, y);
1243 values.push_back(value);
1244 expecteds.push_back(0.0);
1245 }
1246 {
1247 valueFloat3D value(x, z);
1248 values.push_back(value);
1249 expecteds.push_back(0.0);
1250 }
1251 {
1252 valueFloat3D value(y, z);
1253 values.push_back(value);
1254 expecteds.push_back(0.0);
1255 }
1256 vtkMath::Perpendiculars(x, y, z,
1257 0.0);
1258 {
1259 valueFloat3D value(x, y);
1260 values.push_back(value);
1261 expecteds.push_back(0.0);
1262 }
1263 }
1264 for (size_t i = 0; i < values.size(); ++i)
1265 {
1266 float test = vtkMath::Dot(values[i].a, values[i].b);
1267 if (!vtkMathUtilities::FuzzyCompare(
1268 expecteds[i], test,
1269 std::numeric_limits<float>::epsilon() * 256.0f))
1270 {
1271 std::cout << " Perpendiculars got " << test
1272 << " but expected " << expecteds[i];
1273 }
1274 }
1275 }
1276
1277 if (status)
1278 {
1279 std::cout << "..FAILED" << std::endl;
1280 }
1281 else
1282 {
1283 std::cout << ".PASSED" << std::endl;
1284 }
1285 return status;
1286 }
1287
1288 template<typename T>
ProjectVector()1289 int ProjectVector()
1290 {
1291 int status = 0;
1292 T a[3], b[3], c[3];
1293 for (int i = 0; i < 3; ++i)
1294 {
1295 a[i] = 0.0;
1296 b[i] = 0.0;
1297 }
1298 if (vtkMath::ProjectVector(a, b, c))
1299 {
1300 std::cout << "ProjectVector of a 0 vector should return false ";
1301 ++status;
1302 }
1303 return status;
1304 }
1305
1306 // Validate 0 vector case. TestMath does the rest
TestProjectVector()1307 int TestProjectVector()
1308 {
1309 int status = 0;
1310 std::cout << "ProjectVector..";
1311
1312 status += ProjectVector<double>();
1313 status += ProjectVector<float>();
1314
1315 if (status)
1316 {
1317 std::cout << "..FAILED" << std::endl;
1318 }
1319 else
1320 {
1321 std::cout << ".PASSED" << std::endl;
1322 }
1323 return status;
1324 }
1325
1326 template<typename T>
ProjectVector2D()1327 int ProjectVector2D()
1328 {
1329 int status = 0;
1330 T a[2], b[2], c[2];
1331 for (int i = 0; i < 2; ++i)
1332 {
1333 a[i] = 0.0;
1334 b[i] = 0.0;
1335 }
1336 if (vtkMath::ProjectVector2D(a, b, c))
1337 {
1338 std::cout << "ProjectVector2D of a 0 vector should return false ";
1339 ++status;
1340 }
1341 return status;
1342 }
1343
1344 // Validate 0 vector case. TestMath does the rest
TestProjectVector2D()1345 int TestProjectVector2D()
1346 {
1347 int status = 0;
1348 std::cout << "ProjectVector2D..";
1349
1350 status += ProjectVector2D<double>();
1351 status += ProjectVector2D<float>();
1352
1353 if (status)
1354 {
1355 std::cout << "..FAILED" << std::endl;
1356 }
1357 else
1358 {
1359 std::cout << ".PASSED" << std::endl;
1360 }
1361 return status;
1362 }
1363
TestDistance2BetweenPoints()1364 int TestDistance2BetweenPoints()
1365 {
1366 int status = 0;
1367 std::cout << "Distance2BetweenPoints..";
1368
1369 if (status)
1370 {
1371 std::cout << "..FAILED" << std::endl;
1372 }
1373 else
1374 {
1375 std::cout << ".PASSED" << std::endl;
1376 }
1377 return status;
1378 }
1379
TestAngleBetweenVectors()1380 int TestAngleBetweenVectors()
1381 {
1382 int status = 0;
1383 std::cout << "AngleBetweenVectors..";
1384
1385 if (status)
1386 {
1387 std::cout << "..FAILED" << std::endl;
1388 }
1389 else
1390 {
1391 std::cout << ".PASSED" << std::endl;
1392 }
1393 return status;
1394 }
1395
TestGaussianAmplitude()1396 int TestGaussianAmplitude()
1397 {
1398 int status = 0;
1399 std::cout << "GaussianAmplitude..";
1400
1401 if (status)
1402 {
1403 std::cout << "..FAILED" << std::endl;
1404 }
1405 else
1406 {
1407 std::cout << ".PASSED" << std::endl;
1408 }
1409 return status;
1410 }
1411
TestGaussianWeight()1412 int TestGaussianWeight()
1413 {
1414 int status = 0;
1415 std::cout << "GaussianWeight..";
1416
1417 if (status)
1418 {
1419 std::cout << "..FAILED" << std::endl;
1420 }
1421 else
1422 {
1423 std::cout << ".PASSED" << std::endl;
1424 }
1425 return status;
1426 }
1427
1428 class valueDouble2D
1429 {
1430 public:
1431 double a[2];
1432 double b[2];
1433 };
1434 class valueFloat2D
1435 {
1436 public:
1437 float a[2];
1438 float b[2];
1439 };
TestDot2D()1440 int TestDot2D()
1441 {
1442 int status = 0;
1443 std::cout << "Dot2D..";
1444
1445 {
1446 std::vector<valueDouble2D> values;
1447 std::vector<double> expecteds;
1448 for (int n = 0; n < 100; ++n)
1449 {
1450 valueDouble2D v;
1451 double dot = 0.0;
1452 for (int i = 0; i < 2; ++i)
1453 {
1454 v.a[i] = vtkMath::Random();
1455 v.b[i] = vtkMath::Random();
1456 dot += (v.a[i] * v.b[i]);
1457 }
1458 values.push_back(v);
1459 expecteds.push_back(dot);
1460 }
1461 valueDouble2D test;
1462 test.a[0] = 1.0; test.a[1] = 0.0;
1463 test.b[0] = 0.0; test.b[1] = 1.0;
1464 values.push_back(test);
1465 expecteds.push_back(0.0);
1466
1467 for ( size_t i = 0; i < values.size(); ++i)
1468 {
1469 double result = vtkMath::Dot2D(values[i].a, values[i].b);
1470 if (!vtkMathUtilities::FuzzyCompare(
1471 result, expecteds[i],
1472 std::numeric_limits<double>::epsilon()*128.0))
1473 {
1474 std::cout << " Dot got " << result
1475 << " but expected " << expecteds[i];
1476 ++status;
1477 }
1478 }
1479 }
1480
1481 // now float
1482 {
1483 std::vector<valueFloat2D> values;
1484 std::vector<float> expecteds;
1485 for (int n = 0; n < 100; ++n)
1486 {
1487 valueFloat2D v;
1488 float dot = 0.0;
1489 for (int i = 0; i < 2; ++i)
1490 {
1491 v.a[i] = vtkMath::Random();
1492 v.b[i] = vtkMath::Random();
1493 dot += (v.a[i] * v.b[i]);
1494 }
1495 values.push_back(v);
1496 expecteds.push_back(dot);
1497 }
1498 valueFloat2D test;
1499 test.a[0] = 0.0; test.a[1] = 1.0;
1500 test.b[0] = 1.0; test.b[1] = 0.0;
1501 values.push_back(test);
1502 expecteds.push_back(0.0);
1503
1504 for ( size_t i = 0; i < values.size(); ++i)
1505 {
1506 float result = vtkMath::Dot2D(values[i].a, values[i].b);
1507 if (!vtkMathUtilities::FuzzyCompare(
1508 result, expecteds[i],
1509 std::numeric_limits<float>::epsilon()*128.0f))
1510 {
1511 std::cout << " Dot got " << result
1512 << " but expected " << expecteds[i];
1513 ++status;
1514 }
1515 }
1516 }
1517
1518 if (status)
1519 {
1520 std::cout << "..FAILED" << std::endl;
1521 }
1522 else
1523 {
1524 std::cout << ".PASSED" << std::endl;
1525 }
1526 return status;
1527 }
1528
TestNorm2D()1529 int TestNorm2D()
1530 {
1531 int status = 0;
1532 std::cout << "Norm2D..";
1533
1534 if (status)
1535 {
1536 std::cout << "..FAILED" << std::endl;
1537 }
1538 else
1539 {
1540 std::cout << ".PASSED" << std::endl;
1541 }
1542 return status;
1543 }
1544
TestNormalize2D()1545 int TestNormalize2D()
1546 {
1547 int status = 0;
1548 std::cout << "Normalize2D..";
1549
1550 if (status)
1551 {
1552 std::cout << "..FAILED" << std::endl;
1553 }
1554 else
1555 {
1556 std::cout << ".PASSED" << std::endl;
1557 }
1558 return status;
1559 }
1560
TestDeterminant2x2()1561 int TestDeterminant2x2()
1562 {
1563 int status = 0;
1564 std::cout << "Determinant2x2..";
1565 // Frank Matrix
1566 {
1567 double a[2][2];
1568 for (int i = 1; i <= 2; ++i)
1569 {
1570 for (int j = 1; j <= 2; ++j)
1571 {
1572 if (j < i - 2)
1573 {
1574 a[i-1][j-1] = 0.0;
1575 }
1576 else if (j == (i - 1))
1577 {
1578 a[i-1][j-1] = 2 + 1 - i;
1579 }
1580 else
1581 {
1582 a[i-1][j-1] = 2 + 1 - j;
1583 }
1584 }
1585 }
1586 if (vtkMath::Determinant2x2(a[0], a[1]) != 1.0)
1587 {
1588 std::cout << "Determinant2x2 expected " << 1.0
1589 << " but got " << vtkMath::Determinant2x2(a[0], a[1])
1590 << std::endl;
1591 ++status;
1592 };
1593 }
1594 {
1595 float a[2][2];
1596 for (int i = 1; i <= 2; ++i)
1597 {
1598 for (int j = 1; j <= 2; ++j)
1599 {
1600 if (j < i - 2)
1601 {
1602 a[i-1][j-1] = 0.0;
1603 }
1604 else if (j == (i - 1))
1605 {
1606 a[i-1][j-1] = 2 + 1 - i;
1607 }
1608 else
1609 {
1610 a[i-1][j-1] = 2 + 1 - j;
1611 }
1612 }
1613 }
1614 if (vtkMath::Determinant2x2(a[0], a[1]) != 1.0)
1615 {
1616 std::cout << "Determinant2x2 expected " << 1.0
1617 << " but got " << vtkMath::Determinant2x2(a[0], a[1])
1618 << std::endl;
1619 ++status;
1620 };
1621 }
1622 if (status)
1623 {
1624 std::cout << "..FAILED" << std::endl;
1625 }
1626 else
1627 {
1628 std::cout << ".PASSED" << std::endl;
1629 }
1630 return status;
1631 }
1632
TestDeterminant3x3()1633 int TestDeterminant3x3()
1634 {
1635 int status = 0;
1636 std::cout << "Determinant3x3..";
1637
1638 // Frank Matrix
1639 {
1640 double a[3][3];
1641 for (int i = 1; i <= 3; ++i)
1642 {
1643 for (int j = 1; j <= 3; ++j)
1644 {
1645 if (j < i - 3)
1646 {
1647 a[i-1][j-1] = 0.0;
1648 }
1649 else if (j == (i - 1))
1650 {
1651 a[i-1][j-1] = 3 + 1 - i;
1652 }
1653 else
1654 {
1655 a[i-1][j-1] = 3 + 1 - j;
1656 }
1657 }
1658 }
1659 if (vtkMath::Determinant3x3(a[0], a[1], a[2]) != 1.0)
1660 {
1661 std::cout << "Determinant3x3 expected " << 1.0
1662 << " but got " << vtkMath::Determinant3x3(a[0], a[1], a[2])
1663 << std::endl;
1664 ++status;
1665 };
1666 }
1667 {
1668 float a[3][3];
1669 for (int i = 1; i <= 3; ++i)
1670 {
1671 for (int j = 1; j <= 3; ++j)
1672 {
1673 if (j < i - 3)
1674 {
1675 a[i-1][j-1] = 0.0;
1676 }
1677 else if (j == (i - 1))
1678 {
1679 a[i-1][j-1] = 3 + 1 - i;
1680 }
1681 else
1682 {
1683 a[i-1][j-1] = 3 + 1 - j;
1684 }
1685 }
1686 }
1687 if (vtkMath::Determinant3x3(a[0], a[1], a[2]) != 1.0)
1688 {
1689 std::cout << "Determinant3x3 expected " << 1.0
1690 << " but got " << vtkMath::Determinant3x3(a[0], a[1], a[2])
1691 << std::endl;
1692 ++status;
1693 };
1694 }
1695
1696 if (status)
1697 {
1698 std::cout << "..FAILED" << std::endl;
1699 }
1700 else
1701 {
1702 std::cout << ".PASSED" << std::endl;
1703 }
1704 return status;
1705 }
1706
1707 template<typename T>
LUFactor3x3()1708 int LUFactor3x3()
1709 {
1710 int status = 0;
1711
1712 T A[3][3];
1713 int index[3];
1714
1715 for (int n = 0; n < 1000; ++n)
1716 {
1717 for (int i =0; i < 3; ++i)
1718 {
1719 for (int j = 0; j < 3; ++j)
1720 {
1721 A[i][j] = vtkMath::Random(-10.0, 10.0);
1722 }
1723 }
1724 vtkMath::LUFactor3x3(A, index);
1725 }
1726 return status;
1727 }
1728
1729 // Just for coverage, validated as part of TestLUSolve3x3
TestLUFactor3x3()1730 int TestLUFactor3x3()
1731 {
1732 int status = 0;
1733 std::cout << "LUFactor3x3..";
1734
1735 status += LUFactor3x3<double>();
1736
1737 if (status)
1738 {
1739 std::cout << "..FAILED" << std::endl;
1740 }
1741 else
1742 {
1743 std::cout << ".PASSED" << std::endl;
1744 }
1745 return status;
1746 }
1747 template<typename T>
LUSolve3x3()1748 int LUSolve3x3()
1749 {
1750 int status = 0;
1751
1752 // Generate a Hilbert Matrix
1753 T mat[3][3];
1754 int index[3];
1755 T lhs[3];
1756 T rhs[3];
1757
1758 for (int n = 0; n < 1000; ++n)
1759 {
1760 for (int i = 0; i < 3; ++i)
1761 {
1762 lhs[i] = vtkMath::Random(-1.0, 1.0);
1763 }
1764
1765 for (int i = 1; i <= 3; ++i)
1766 {
1767 rhs[i-1] = 0.0;
1768 for (int j = 1; j <= 3; ++j)
1769 {
1770 mat[i-1][j-1] = 1.0 / (i + j - 1);
1771 rhs[i-1] += (mat[i-1][j-1] * lhs[j-1]);
1772 }
1773 }
1774 vtkMath::LUFactor3x3(mat, index);
1775 vtkMath::LUSolve3x3(mat, index, rhs);
1776 for (int i = 0; i < 3; ++i)
1777 {
1778 if (!vtkMathUtilities::FuzzyCompare(
1779 lhs[i], rhs[i],
1780 std::numeric_limits<T>::epsilon() * (T) 256.0))
1781 {
1782 std::cout << " LUSolve3x3(T) expected " << lhs[i]
1783 << " but got " << rhs[i];
1784 ++status;
1785 }
1786 }
1787 }
1788 return status;
1789 }
1790
TestLUSolve3x3()1791 int TestLUSolve3x3()
1792 {
1793 int status = 0;
1794 std::cout << "LUSolve3x3..";
1795
1796 status += LUSolve3x3<double>();
1797 status += LUSolve3x3<float>();
1798
1799 if (status)
1800 {
1801 std::cout << "..FAILED" << std::endl;
1802 }
1803 else
1804 {
1805 std::cout << ".PASSED" << std::endl;
1806 }
1807 return status;
1808 }
1809
1810 template<typename T>
LinearSolve3x3()1811 int LinearSolve3x3()
1812 {
1813 int status = 0;
1814
1815 // Generate a Hilbert Matrix
1816 T mat[3][3];
1817 T lhs[3];
1818 T rhs[3];
1819 T solution[3];
1820
1821 for (int n = 0; n < 2; ++n)
1822 {
1823 for (int i = 0; i < 3; ++i)
1824 {
1825 lhs[i] = vtkMath::Random(-1.0, 1.0);
1826 }
1827
1828 for (int i = 1; i <= 3; ++i)
1829 {
1830 rhs[i-1] = 0.0;
1831 for (int j = 1; j <= 3; ++j)
1832 {
1833 mat[i-1][j-1] = 1.0 / (i + j - 1);
1834 rhs[i-1] += (mat[i-1][j-1] * lhs[j-1]);
1835 }
1836 }
1837 vtkMath::LinearSolve3x3(mat, rhs, solution);
1838
1839 for (int i = 0; i < 3; ++i)
1840 {
1841 if (!vtkMathUtilities::FuzzyCompare(
1842 lhs[i], solution[i],
1843 std::numeric_limits<T>::epsilon() * (T) 512.0))
1844 {
1845 std::cout << " LinearSolve3x3(T) expected " << lhs[i]
1846 << " but got " << solution[i];
1847 ++status;
1848 }
1849 }
1850 }
1851 return status;
1852 }
1853
TestLinearSolve3x3()1854 int TestLinearSolve3x3()
1855 {
1856 int status = 0;
1857 std::cout << "LinearSolve3x3..";
1858
1859 status += LinearSolve3x3<double>();
1860 status += LinearSolve3x3<float>();
1861
1862 if (status)
1863 {
1864 std::cout << "..FAILED" << std::endl;
1865 }
1866 else
1867 {
1868 std::cout << ".PASSED" << std::endl;
1869 }
1870 return status;
1871 }
1872
1873 template<typename T>
Multiply3x3()1874 int Multiply3x3()
1875 {
1876 int status = 0;
1877 T A[3][3];
1878 T V[3];
1879 T U[3];
1880
1881 for (int i = 0; i < 3; ++i)
1882 {
1883 for (int j = 0; j < 3; ++j)
1884 {
1885 A[i][j] = vtkMath::Random(-10.0, 10.0);
1886 }
1887 V[i] = vtkMath::Random(-10.0, 10.0);
1888 }
1889
1890 vtkMath::Multiply3x3(A, V, U);
1891
1892 return status;
1893 }
1894
TestMultiply3x3()1895 int TestMultiply3x3()
1896 {
1897 int status = 0;
1898 std::cout << "Multiply3x3..";
1899
1900 status += Multiply3x3<double>();
1901 status += Multiply3x3<float>();
1902
1903 if (status)
1904 {
1905 std::cout << "..FAILED" << std::endl;
1906 }
1907 else
1908 {
1909 std::cout << ".PASSED" << std::endl;
1910 }
1911 return status;
1912 }
1913
1914 // For coverage only. Validated as part of TestInvertMatrix
TestMultiplyMatrix()1915 int TestMultiplyMatrix()
1916 {
1917 int status = 0;
1918 std::cout << "MultiplyMatrix..";
1919
1920 double a[3][3] =
1921 {{1.0, 2.0, 3.0},
1922 {4.0, 5.0, 6.0},
1923 {7.0, 8.0, 9.0}};
1924 double b[3][3] =
1925 {{1.0, 1.0, 1.0},
1926 {1.0, 1.0, 1.0},
1927 {1.0, 1.0, 1.0}};
1928 double c[3][3];
1929
1930 double *aa[3];
1931 double *bb[3];
1932 double *cc[3];
1933 for (int i = 0; i < 3; ++i)
1934 {
1935 aa[i] = a[i];
1936 bb[i] = b[i];
1937 cc[i] = c[i];
1938 }
1939 vtkMath::MultiplyMatrix(aa, bb, 3, 3, 3, 3, cc);
1940
1941 // WARNING: Number of columns of A must match number of rows of B.
1942 vtkMath::MultiplyMatrix(aa, bb, 3, 2, 3, 3, cc);
1943
1944 if (status)
1945 {
1946 std::cout << "..FAILED" << std::endl;
1947 }
1948 else
1949 {
1950 std::cout << ".PASSED" << std::endl;
1951 }
1952 return status;
1953 }
1954
TestTranspose3x3()1955 int TestTranspose3x3()
1956 {
1957 int status = 0;
1958 std::cout << "Transpose3x3..";
1959
1960 if (status)
1961 {
1962 std::cout << "..FAILED" << std::endl;
1963 }
1964 else
1965 {
1966 std::cout << ".PASSED" << std::endl;
1967 }
1968 return status;
1969 }
1970
TestInvert3x3()1971 int TestInvert3x3()
1972 {
1973 int status = 0;
1974 std::cout << "Invert3x3..";
1975 {
1976 // Generate a Hilbert Matrix
1977 double mat[3][3];
1978 double matI[3][3];
1979 double expected[3][3] =
1980 {{9.0, -36.0, 30.0},
1981 {-36.0, 192.0, -180.0},
1982 {30.0, -180.0, 180.0}};
1983
1984 for (int i = 1; i <= 3; ++i)
1985 {
1986 for (int j = 1; j <= 3; ++j)
1987 {
1988 mat[i-1][j-1] = 1.0 / (i + j - 1);
1989 }
1990 }
1991 vtkMath::Invert3x3(mat, matI);
1992 for (int i = 0; i < 3; ++i)
1993 {
1994 for (int j = 0; j < 3; ++j)
1995 {
1996 if (!vtkMathUtilities::FuzzyCompare(
1997 matI[i][j], expected[i][j],
1998 std::numeric_limits<double>::epsilon()*16384.0))
1999 {
2000 std::cout << " Invert3x3(double) expected " << expected[i][j]
2001 << " but got " << matI[i][j];
2002 ++status;
2003 }
2004 }
2005 }
2006 }
2007 {
2008 // Generate a Hilbert Matrix
2009 float mat[3][3];
2010 float matI[3][3];
2011 float expected[3][3] =
2012 {{9.0f, -36.0f, 30.0f},
2013 {-36.0f, 192.0f, -180.0f},
2014 {30.0f, -180.0f, 180.0f}};
2015
2016 for (int i = 1; i <= 3; ++i)
2017 {
2018 for (int j = 1; j <= 3; ++j)
2019 {
2020 mat[i-1][j-1] = 1.0 / (i + j - 1);
2021 }
2022 }
2023 vtkMath::Invert3x3(mat, matI);
2024 for (int i = 0; i < 3; ++i)
2025 {
2026 for (int j = 0; j < 3; ++j)
2027 {
2028 if (!vtkMathUtilities::FuzzyCompare(
2029 matI[i][j], expected[i][j],
2030 std::numeric_limits<float>::epsilon()*8192.0f))
2031 {
2032 std::cout << " Invert3x3(single) expected " << expected[i][j]
2033 << " but got " << matI[i][j];
2034 ++status;
2035 }
2036 }
2037 }
2038 }
2039 if (status)
2040 {
2041 std::cout << "..FAILED" << std::endl;
2042 }
2043 else
2044 {
2045 std::cout << ".PASSED" << std::endl;
2046 }
2047 return status;
2048 }
2049
2050 template<typename T, int NDimension>
InvertMatrix()2051 int InvertMatrix()
2052 {
2053 int status = 0;
2054
2055 // Generate a Hilbert Matrix
2056 T **mat = new T *[NDimension];
2057 T **orig = new T *[NDimension];
2058 T **matI = new T *[NDimension];
2059 T **ident = new T *[NDimension];
2060 int *tmp1 = new int [NDimension];
2061 T *tmp2 = new T [NDimension];
2062 for (int i = 1; i <= NDimension; ++i)
2063 {
2064 mat[i-1] = new T[NDimension];
2065 matI[i-1] = new T[NDimension];
2066 orig[i-1] = new T[NDimension];
2067 ident[i-1] = new T[NDimension];
2068 for (int j = 1; j <= NDimension; ++j)
2069 {
2070 orig[i-1][j-1] = mat[i-1][j-1] = 1.0 / (i + j - 1);
2071 }
2072 }
2073 if (vtkMath::InvertMatrix(mat, matI, NDimension, tmp1, tmp2) == 0)
2074 {
2075 delete [] mat;
2076 delete [] orig;
2077 delete [] matI;
2078 delete [] ident;
2079 delete [] tmp1;
2080 delete [] tmp2;
2081 return status;
2082 }
2083 vtkMath::MultiplyMatrix(orig,
2084 matI,
2085 NDimension, NDimension, NDimension, NDimension,
2086 ident);
2087
2088 T expected;
2089 for (int i = 0; i < NDimension; ++i)
2090 {
2091 for (int j = 0; j < NDimension; ++j)
2092 {
2093 if (i == j)
2094 {
2095 expected = 1.0;
2096 }
2097 else
2098 {
2099 expected = 0.0;
2100 }
2101 if (!vtkMathUtilities::FuzzyCompare(
2102 *(ident[i] + j), expected,
2103 std::numeric_limits<T>::epsilon() * (T) 100000.0))
2104 {
2105 std::cout << " InvertMatrix(T) expected " << expected
2106 << " but got " << *(ident[i] + j);
2107 std::cout << "eps ratio is: " << (*(ident[i] + j) - expected)
2108 / std::numeric_limits<T>::epsilon() << std::endl;
2109 ++status;
2110 }
2111 }
2112 }
2113 for (int i = 0; i < NDimension; i++)
2114 {
2115 delete [] mat[i];
2116 delete [] orig[i];
2117 delete [] matI[i];
2118 delete [] ident[i];
2119 }
2120 delete [] mat;
2121 delete [] orig;
2122 delete [] matI;
2123 delete [] ident;
2124 delete [] tmp1;
2125 delete [] tmp2;
2126 return status;
2127 }
TestInvertMatrix()2128 int TestInvertMatrix()
2129 {
2130 int status = 0;
2131 status += InvertMatrix<double,3> ();
2132 status += InvertMatrix<double,4> ();
2133 status += InvertMatrix<double,5> ();
2134
2135 if (status)
2136 {
2137 std::cout << "..FAILED" << std::endl;
2138 }
2139 else
2140 {
2141 std::cout << ".PASSED" << std::endl;
2142 }
2143 return status;
2144 }
2145
TestIdentity3x3()2146 int TestIdentity3x3()
2147 {
2148 int status = 0;
2149 std::cout << "Identity3x3..";
2150
2151 double m[3][3];
2152 vtkMath::Identity3x3(m);
2153
2154 double expected;
2155 for (int i = 0; i < 3; i++)
2156 {
2157 for (int j = 0; j < 3; j++)
2158 {
2159 if (i == j)
2160 {
2161 expected = 1.0;
2162 }
2163 else
2164 {
2165 expected = 0.0;
2166 }
2167 if (expected != m[i][j])
2168 {
2169 std::cout << " Identity expected " << expected
2170 << " but got " << m[i][j] << std::endl;
2171 ++status;
2172 }
2173 }
2174 }
2175
2176 if (status)
2177 {
2178 std::cout << "..FAILED" << std::endl;
2179 }
2180 else
2181 {
2182 std::cout << ".PASSED" << std::endl;
2183 }
2184 return status;
2185 }
2186
2187 template<typename T>
Matrix3x3ToQuaternion()2188 int Matrix3x3ToQuaternion()
2189 {
2190 int status = 0;
2191
2192 T A[3][3];
2193 T quat[4];
2194
2195 for (int n = 0; n < 1000; ++n)
2196 {
2197 for (int i = 0; i < 3; ++i)
2198 {
2199 for (int j = 0; j < 3; ++j)
2200 {
2201 A[i][j] = vtkMath::Random(-1.0, 1.0);
2202 }
2203 }
2204 vtkMath::Matrix3x3ToQuaternion(A, quat);
2205 }
2206 return status;
2207
2208 }
TestMatrix3x3ToQuaternion()2209 int TestMatrix3x3ToQuaternion()
2210 {
2211 int status = 0;
2212 std::cout << "Matrix3x3ToQuaternion..";
2213
2214 status += Matrix3x3ToQuaternion<double>();
2215 status += Matrix3x3ToQuaternion<float>();
2216
2217 if (status)
2218 {
2219 std::cout << "..FAILED" << std::endl;
2220 }
2221 else
2222 {
2223 std::cout << ".PASSED" << std::endl;
2224 }
2225 return status;
2226 }
2227
2228 template<typename T>
QuaternionToMatrix3x3()2229 int QuaternionToMatrix3x3()
2230 {
2231 int status = 0;
2232
2233 T A[3][3];
2234 T quat[4];
2235
2236 for (int n = 0; n < 1000; ++n)
2237 {
2238 quat[0] = vtkMath::Random(-vtkMath::Pi(), vtkMath::Pi());
2239 for (int i = 1; i < 2; ++i)
2240 {
2241 quat[i] = vtkMath::Random(-10.0,10.0);
2242
2243 }
2244 vtkMath::QuaternionToMatrix3x3(quat, A);
2245 }
2246 return status;
2247
2248 }
2249
TestQuaternionToMatrix3x3()2250 int TestQuaternionToMatrix3x3()
2251 {
2252 int status = 0;
2253 std::cout << "QuaternionToMatrix3x3..";
2254
2255 status += QuaternionToMatrix3x3<double>();
2256 status += QuaternionToMatrix3x3<float>();
2257
2258 if (status)
2259 {
2260 std::cout << "..FAILED" << std::endl;
2261 }
2262 else
2263 {
2264 std::cout << ".PASSED" << std::endl;
2265 }
2266 return status;
2267 }
2268
2269 template<typename T>
MultiplyQuaternion()2270 int MultiplyQuaternion()
2271 {
2272 int status = 0;
2273
2274 T q1[4];
2275 T q2[4];
2276 T q3[4];
2277 for (int n = 0; n < 1000; ++n)
2278 {
2279 q1[0] = vtkMath::Random(-vtkMath::Pi(), vtkMath::Pi());
2280 q2[0] = vtkMath::Random(-vtkMath::Pi(), vtkMath::Pi());
2281 vtkMath::MultiplyQuaternion(q1, q2, q3);
2282 }
2283
2284 return status;
2285 }
TestMultiplyQuaternion()2286 int TestMultiplyQuaternion()
2287 {
2288 int status = 0;
2289 std::cout << "MultiplyQuaternion..";
2290
2291 status += MultiplyQuaternion<double>();
2292 status += MultiplyQuaternion<float>();
2293
2294 if (status)
2295 {
2296 std::cout << "..FAILED" << std::endl;
2297 }
2298 else
2299 {
2300 std::cout << ".PASSED" << std::endl;
2301 }
2302 return status;
2303 }
2304
TestOrthogonalize3x3()2305 int TestOrthogonalize3x3()
2306 {
2307 int status = 0;
2308 std::cout << "Orthogonalize3x3..";
2309
2310 {
2311 // Generate a random matrix
2312 double mat[3][3];
2313 double matO[3][3];
2314 double matI[3][3];
2315
2316 for (int n = 0; n < 1000; ++n)
2317 {
2318 for (int i = 0; i < 3; ++i)
2319 {
2320 for (int j = 0; j < 3; ++j)
2321 {
2322 mat[i][j] = vtkMath::Random();
2323 }
2324 }
2325 vtkMath::Orthogonalize3x3(mat, matO);
2326 vtkMath::Transpose3x3(matO, mat);
2327 vtkMath::Multiply3x3(mat, matO, matI);
2328
2329 double identity[3][3];
2330 vtkMath::Identity3x3(identity);
2331 for (int i = 0; i < 3; ++i)
2332 {
2333 for (int j = 0; j < 3; ++j)
2334 {
2335 if (!vtkMathUtilities::FuzzyCompare(
2336 matI[i][j], identity[i][j],
2337 std::numeric_limits<double>::epsilon()*128.0))
2338 {
2339 std::cout << " Orthogonalize3x3 expected " << identity[i][j]
2340 << " but got " << matI[i][j];
2341 ++status;
2342 }
2343 }
2344 }
2345 }
2346 }
2347
2348 {
2349 // Generate a random matrix
2350 float mat[3][3];
2351 float matO[3][3];
2352 float matI[3][3];
2353
2354 for (int n = 0; n < 1000; ++n)
2355 {
2356 for (int i = 0; i < 3; ++i)
2357 {
2358 for (int j = 0; j < 3; ++j)
2359 {
2360 mat[i][j] = vtkMath::Random();
2361 }
2362 }
2363 vtkMath::Orthogonalize3x3(mat, matO);
2364 vtkMath::Transpose3x3(matO, mat);
2365 vtkMath::Multiply3x3(mat, matO, matI);
2366
2367 float identity[3][3];
2368 vtkMath::Identity3x3(identity);
2369 for (int i = 0; i < 3; ++i)
2370 {
2371 for (int j = 0; j < 3; ++j)
2372 {
2373 if (!vtkMathUtilities::FuzzyCompare(
2374 matI[i][j], identity[i][j],
2375 std::numeric_limits<float>::epsilon()*128.0f))
2376 {
2377 std::cout << " Orthogonalize3x3 expected " << identity[i][j]
2378 << " but got " << matI[i][j];
2379 ++status;
2380 }
2381 }
2382 }
2383 }
2384 }
2385
2386 if (status)
2387 {
2388 std::cout << "..FAILED" << std::endl;
2389 }
2390 else
2391 {
2392 std::cout << ".PASSED" << std::endl;
2393 }
2394 return status;
2395 }
2396
2397 template<typename T>
Diagonalize3x3()2398 int Diagonalize3x3()
2399 {
2400 int status = 0;
2401
2402 T mat[3][3];
2403 T eigenVector[3][3], eigenVectorT[3][3];
2404 T temp[3][3];
2405 T result[3][3];
2406 T eigen[3];
2407
2408 for (int n = 0; n < 0; ++n)
2409 {
2410 for (int i = 0; i < 3; ++i)
2411 {
2412 for (int j = i; j < 3; ++j)
2413 {
2414 mat[i][j] = mat[j][i] = vtkMath::Random(-1.0, 1.0);
2415 }
2416 }
2417
2418 vtkMath::Diagonalize3x3(mat, eigen, eigenVector);
2419
2420 // Pt * A * P = diagonal matrix with eigenvalues on diagonal
2421 vtkMath::Multiply3x3(mat, eigenVector, temp);
2422 vtkMath::Invert3x3(eigenVector, eigenVectorT);
2423 vtkMath::Multiply3x3(eigenVectorT, temp, result);
2424 T expected;
2425 for (int i = 0; i < 3; ++i)
2426 {
2427 for (int j = 0; j < 3; ++j)
2428 {
2429 if (i == j)
2430 {
2431 expected = eigen[i];
2432 }
2433 else
2434 {
2435 expected = 0.0;
2436 }
2437 if (!vtkMathUtilities::FuzzyCompare(
2438 result[i][j], expected,
2439 std::numeric_limits<T>::epsilon()*(T)128.0))
2440 {
2441 std::cout << " Diagonalize3x3 expected " << expected
2442 << " but got " << result[i][j];
2443 ++status;
2444 }
2445 }
2446 }
2447 }
2448
2449 // Now test for 2 and 3 equal eigenvalues
2450 vtkMath::Identity3x3(mat);
2451 mat[0][0] = 5.0;
2452 mat[1][1] = 5.0;
2453 mat[2][2] = 1.0;
2454
2455 vtkMath::Diagonalize3x3(mat, eigen, eigenVector);
2456 std::cout << "eigen: "
2457 << eigen[0] << ","
2458 << eigen[1] << ","
2459 << eigen[2] << std::endl;
2460
2461 vtkMath::Identity3x3(mat);
2462 mat[0][0] = 2.0;
2463 mat[1][1] = 2.0;
2464 mat[2][2] = 2.0;
2465
2466 vtkMath::Diagonalize3x3(mat, eigen, eigenVector);
2467 std::cout << "eigen: "
2468 << eigen[0] << ","
2469 << eigen[1] << ","
2470 << eigen[2] << std::endl;
2471 return status;
2472 }
2473
2474 // Validate Pt * A * P = diagonal matrix with eigenvalues on diagonal
TestDiagonalize3x3()2475 int TestDiagonalize3x3()
2476 {
2477 int status = 0;
2478 std::cout << "Diagonalize3x3..";
2479
2480 status += Diagonalize3x3<double>();
2481 status += Diagonalize3x3<float>();
2482
2483 if (status)
2484 {
2485 std::cout << "..FAILED" << std::endl;
2486 }
2487 else
2488 {
2489 std::cout << ".PASSED" << std::endl;
2490 }
2491 return status;
2492 }
2493
2494 template<typename T>
SingularValueDecomposition3x3()2495 int SingularValueDecomposition3x3()
2496 {
2497 int status = 0;
2498
2499 T a[3][3];
2500 T orig[3][3];
2501 T u[3][3];
2502 T w[3];
2503 T vt[3][3];
2504
2505 for (int n = 0; n < 1000; ++n)
2506 {
2507 for (int i = 0; i < 3; ++i)
2508 {
2509 for (int j = 0; j < 3; ++j)
2510 {
2511 orig[i][j] = a[i][j] = vtkMath::Random(-10.0, 10.0);
2512 }
2513 }
2514 vtkMath::SingularValueDecomposition3x3(a, u, w, vt);
2515
2516 T m[3][3];
2517 T W[3][3];
2518 vtkMath::Identity3x3(W);
2519 W[0][0] = w[0]; W[1][1] = w[1]; W[2][2] = w[2];
2520 vtkMath::Multiply3x3(u, W, m);
2521 vtkMath::Multiply3x3(m, vt, m);
2522
2523 for (int i = 0; i < 3; ++i)
2524 {
2525 for (int j = 0; j < 3; ++j)
2526 {
2527 if (!vtkMathUtilities::FuzzyCompare(
2528 m[i][j], orig[i][j],
2529 std::numeric_limits<T>::epsilon() * (T) 128.0))
2530 {
2531 std::cout << " SingularValueDecomposition3x3 expected " << orig[i][j]
2532 << " but got " << m[i][j];
2533 std::cout << " eps ratio is: " << (m[i][j] - orig[i][j])
2534 / std::numeric_limits<T>::epsilon() << std::endl;
2535 ++status;
2536 }
2537 }
2538 }
2539 }
2540 return status;
2541 }
2542 // Validate u * w * vt = m
TestSingularValueDecomposition3x3()2543 int TestSingularValueDecomposition3x3()
2544 {
2545 int status = 0;
2546 std::cout << "SingularValueDecomposition3x3..";
2547
2548 status += SingularValueDecomposition3x3<double>();
2549 status += SingularValueDecomposition3x3<float>();
2550
2551 if (status)
2552 {
2553 std::cout << "..FAILED" << std::endl;
2554 }
2555 else
2556 {
2557 std::cout << ".PASSED" << std::endl;
2558 }
2559 return status;
2560 }
2561
2562
2563 template<typename T, int NDimension>
SolveLinearSystem()2564 int SolveLinearSystem()
2565 {
2566 int status = 0;
2567
2568 for (int n = 0; n < 100; ++n)
2569 {
2570 // Generate a Random Matrix
2571 T **mat = new T *[NDimension];
2572 T *lhs = new T [NDimension];
2573 T *rhs = new T [NDimension];
2574
2575 for (int i = 0; i < NDimension; ++i)
2576 {
2577 mat[i] = new T[NDimension];
2578 lhs[i] =vtkMath::Random(-1.0, 1.0);
2579 for (int j = 0; j < NDimension; ++j)
2580 {
2581 *(mat[i] + j) = vtkMath::Random(-1.0, 1.0);
2582 }
2583 }
2584
2585 for (int i = 0; i < NDimension; ++i)
2586 {
2587 rhs[i] = 0.0;
2588 for (int j = 0; j < NDimension; ++j)
2589 {
2590 rhs[i] += (*(mat[i] + j) * lhs[j]);
2591 }
2592 }
2593 vtkMath::SolveLinearSystem(mat, rhs, NDimension);
2594
2595 for (int i = 0; i < NDimension; ++i)
2596 {
2597 if (!vtkMathUtilities::FuzzyCompare(
2598 lhs[i], rhs[i],
2599 std::numeric_limits<double>::epsilon()*32768.0))
2600 {
2601 std::cout << " SolveLinearSystem(double) expected " << lhs[i]
2602 << " but got " << rhs[i];
2603 std::cout << " eps ratio is: " << (lhs[i] - rhs[i])
2604 / std::numeric_limits<T>::epsilon() << std::endl;
2605 ++status;
2606 }
2607 }
2608
2609 if (NDimension == 1 || NDimension == 2)
2610 {
2611 for (int i = 0; i < NDimension; ++i)
2612 {
2613 for (int j = 0; j < NDimension; ++j)
2614 {
2615 *(mat[i] + j) = 0.0;
2616 }
2617 }
2618 if (vtkMath::SolveLinearSystem(mat, rhs, NDimension) != 0.0)
2619 {
2620 std::cout << " SolveLinearSystem for a zero matrix expected " << 0
2621 << " but got 1" ;
2622 }
2623 }
2624 for (int i = 0; i < NDimension; i++)
2625 {
2626 delete [] mat[i];
2627 }
2628 delete [] mat;
2629 delete [] rhs;
2630 delete [] lhs;
2631 }
2632 return status;
2633 }
2634
2635 // Validate with a known left hand side
TestSolveLinearSystem()2636 int TestSolveLinearSystem()
2637 {
2638 int status = 0;
2639 std::cout << "SolveLinearSystem..";
2640
2641 status += SolveLinearSystem<double,1> ();
2642 status += SolveLinearSystem<double,2> ();
2643 status += SolveLinearSystem<double,3> ();
2644 status += SolveLinearSystem<double,50> ();
2645
2646 if (status)
2647 {
2648 std::cout << "..FAILED" << std::endl;
2649 }
2650 else
2651 {
2652 std::cout << ".PASSED" << std::endl;
2653 }
2654 return status;
2655 }
2656
2657 // Validate with a known solution
TestSolveLeastSquares()2658 int TestSolveLeastSquares()
2659 {
2660 int status = 0;
2661 std::cout << "SolveLeastSquares..";
2662
2663 double **m = new double *[2];
2664 double **x = new double *[3];
2665 double **y = new double *[3];
2666
2667 for (int i = 0; i < 2; ++i)
2668 {
2669 m[i] = new double[1];
2670 }
2671 for (int i = 0; i < 3; ++i)
2672 {
2673 x[i] = new double[2];
2674 }
2675 x[0][0] = 1; x[0][1] = 4;
2676 x[1][0] = 1; x[1][1] = 2;
2677 x[2][0] = 2; x[2][1] = 3;
2678
2679 for (int i = 0; i < 3; ++i)
2680 {
2681 y[i] = new double[1];
2682 }
2683 y[0][0] = -2;
2684 y[1][0] = 6;
2685 y[2][0] = 1;
2686
2687 vtkMath::SolveLeastSquares(3, x, 2, y, 1, m);
2688
2689 std::vector<double> results;
2690 std::vector<double> expecteds;
2691 expecteds.push_back(3.0);
2692 results.push_back(m[0][0]);
2693 expecteds.push_back(-1.0);
2694 results.push_back(m[1][0]);
2695
2696 for ( size_t i = 0; i < results.size(); ++i)
2697 {
2698 if (!vtkMathUtilities::FuzzyCompare(
2699 results[i], expecteds[i],
2700 std::numeric_limits<double>::epsilon()*128.0))
2701 {
2702 std::cout << " Solve Least Squares got " << results[i]
2703 << " but expected " << expecteds[i];
2704 ++status;
2705 }
2706 }
2707
2708 // Now make one solution homogeous
2709 y[0][0] = 0.0;
2710 vtkMath::SolveLeastSquares(3, x, 2, y, 1, m);
2711
2712 // Now make all homogeous
2713 y[0][0] = 0.0;
2714 y[1][0] = 0.0;
2715 y[2][0] = 0.0;
2716 vtkMath::SolveLeastSquares(3, x, 2, y, 1, m);
2717
2718 // Insufficient number of samples. Underdetermined.
2719 if (vtkMath::SolveLeastSquares(1, x, 2, y, 1, m) != 0)
2720 {
2721 std::cout << " Solve Least Squares got " << 1
2722 << " but expected " << 0;
2723 ++status;
2724 }
2725
2726 if (status)
2727 {
2728 std::cout << "..FAILED" << std::endl;
2729 }
2730 else
2731 {
2732 std::cout << ".PASSED" << std::endl;
2733 }
2734
2735 for (int i = 0; i < 3; i++)
2736 {
2737 delete [] x[i];
2738 delete [] y[i];
2739 }
2740 for (int i = 0; i < 2; i++)
2741 {
2742 delete [] m[i];
2743 }
2744 delete [] x;
2745 delete [] y;
2746 delete [] m;
2747
2748 return status;
2749 }
2750
2751 // Only warning cases validate
2752 // No validation, just coverage
TestSolveHomogeneousLeastSquares()2753 int TestSolveHomogeneousLeastSquares()
2754 {
2755 int status = 0;
2756 std::cout << "SolveHomogenousLeastSquares..";
2757
2758 double **m = new double *[2];
2759 double **x = new double *[3];
2760 double **y = new double *[3];
2761
2762 for (int i = 0; i < 2; ++i)
2763 {
2764 m[i] = new double[1];
2765 }
2766 for (int i = 0; i < 3; ++i)
2767 {
2768 x[i] = new double[2];
2769 y[i] = new double[1];
2770 }
2771 x[0][0] = 1; x[0][1] = 2;
2772 x[1][0] = 2; x[1][1] = 4;
2773 x[2][0] = 3; x[2][1] = 6;
2774
2775 vtkMath::SolveHomogeneousLeastSquares(3, x, 1, m);
2776
2777 vtkMath::MultiplyMatrix(x, m, 3, 2, 2, 1, y);
2778
2779 std::vector<double> results;
2780 std::vector<double> expecteds;
2781
2782 for ( size_t i = 0; i < results.size(); ++i)
2783 {
2784 if (!vtkMathUtilities::FuzzyCompare(
2785 results[i], expecteds[i],
2786 std::numeric_limits<double>::epsilon()*128.0))
2787 {
2788 std::cout << " SolveHomogeneousLeastSquares got " << results[i]
2789 << " but expected " << expecteds[i];
2790 ++status;
2791 }
2792 }
2793
2794 // Insufficient number of samples. Underdetermined.
2795 if (vtkMath::SolveHomogeneousLeastSquares(3, x, 5, m) != 0)
2796 {
2797 std::cout << " SolveHomogeneousLeastSquares got " << 1
2798 << " but expected " << 0;
2799 ++status;
2800 }
2801
2802 if (status)
2803 {
2804 std::cout << "..FAILED" << std::endl;
2805 }
2806 else
2807 {
2808 std::cout << ".PASSED" << std::endl;
2809 }
2810
2811 for (int i = 0; i < 3; i++)
2812 {
2813 delete [] x[i];
2814 delete [] y[i];
2815 }
2816 for (int i = 0; i < 2; i++)
2817 {
2818 delete [] m[i];
2819 }
2820 delete [] x;
2821 delete [] y;
2822 delete [] m;
2823
2824 return status;
2825 }
2826
2827 template<typename T, int NDimension>
LUSolveLinearSystemEstimateMatrixCondition()2828 int LUSolveLinearSystemEstimateMatrixCondition()
2829 {
2830 int status = 0;
2831
2832 // Generate a Hilbert Matrix
2833 T **mat = new T *[NDimension];
2834 int index[NDimension];
2835
2836 for (int i = 1; i <= NDimension; ++i)
2837 {
2838 mat[i-1] = new T[NDimension];
2839 for (int j = 1; j <= NDimension; ++j)
2840 {
2841 mat[i-1][j-1] = 1.0 / (i + j - 1);
2842 }
2843 }
2844 vtkMath::LUFactorLinearSystem(mat, index, NDimension);
2845 T condition = vtkMath::EstimateMatrixCondition(mat, NDimension);
2846 std::cout << "Condition is: " << condition << std::endl;
2847
2848 T expected = condition;
2849 if (!vtkMathUtilities::FuzzyCompare(
2850 condition, expected,
2851 std::numeric_limits<T>::epsilon() * (T) 128.0))
2852 {
2853 std::cout << " EstimateMatrixCondition(T) expected " << expected
2854 << " but got " << condition;
2855 std::cout << "eps ratio is: " << condition - expected
2856 / std::numeric_limits<T>::epsilon() << std::endl;
2857 ++status;
2858 }
2859 for (int i = 0; i < NDimension; i++)
2860 {
2861 delete [] mat[i];
2862 }
2863 delete [] mat;
2864 return status;
2865 }
2866
2867 // Validate by obervation that the condition of a hilbert matrix
2868 // increases with dimension
TestLUSolveLinearSystemEstimateMatrixCondition()2869 int TestLUSolveLinearSystemEstimateMatrixCondition()
2870 {
2871 int status = 0;
2872 std::cout << "LUSolveLinearSystemEstimateMatrixCondition..";
2873 status += LUSolveLinearSystemEstimateMatrixCondition<double,10>();
2874 status += LUSolveLinearSystemEstimateMatrixCondition<double,8>();
2875 status += LUSolveLinearSystemEstimateMatrixCondition<double,6>();
2876 status += LUSolveLinearSystemEstimateMatrixCondition<double,4>();
2877 status += LUSolveLinearSystemEstimateMatrixCondition<double,3>();
2878 if (status)
2879 {
2880 std::cout << "..FAILED" << std::endl;
2881 }
2882 else
2883 {
2884 std::cout << ".PASSED" << std::endl;
2885 }
2886 return status;
2887 }
2888
2889 template<typename T, int NDimension>
JacobiN()2890 int JacobiN()
2891 {
2892 int status = 0;
2893
2894 T mat[NDimension][NDimension];
2895 T orig[NDimension][NDimension];
2896 T eigenVector[NDimension][NDimension], eigenVectorT[NDimension][NDimension];
2897 T temp[NDimension][NDimension];
2898 T result[NDimension][NDimension];
2899 T eigen[NDimension];
2900
2901 for (int n = 0; n < 10; ++n)
2902 {
2903 for (int i = 0; i < NDimension; ++i)
2904 {
2905 for (int j = i; j < NDimension; ++j)
2906 {
2907 mat[i][j] = mat[j][i] = vtkMath::Random(0.0, 1.0);
2908 orig[i][j] = orig[j][i] = mat[i][j];
2909 }
2910 }
2911
2912 // convert to jacobiN format
2913 T *origJ[NDimension];
2914 T *matJ[NDimension];
2915 T *eigenVectorJ[NDimension];
2916 T *eigenVectorTJ[NDimension];
2917 T *resultJ[NDimension];
2918 T *tempJ[NDimension];
2919 for (int i = 0; i < NDimension; ++i)
2920 {
2921 matJ[i] = mat[i];
2922 origJ[i] = orig[i];
2923 eigenVectorJ[i] = eigenVector[i];
2924 eigenVectorTJ[i] = eigenVectorT[i];
2925 tempJ[i] = temp[i];
2926 resultJ[i] = result[i];
2927 }
2928
2929 if (NDimension == 3)
2930 {
2931 vtkMath::Jacobi(matJ,
2932 eigen,
2933 eigenVectorJ);
2934 }
2935 else
2936 {
2937 vtkMath::JacobiN(matJ,
2938 NDimension,
2939 eigen,
2940 eigenVectorJ);
2941 }
2942
2943 // P^-1 * A * P = diagonal matrix with eigenvalues on diagonal
2944 vtkMath::MultiplyMatrix(origJ,
2945 eigenVectorJ,
2946 NDimension, NDimension, NDimension, NDimension,
2947 tempJ);
2948 vtkMath::InvertMatrix(eigenVectorJ, eigenVectorTJ, NDimension);
2949 vtkMath::MultiplyMatrix(eigenVectorTJ,
2950 tempJ,
2951 NDimension, NDimension, NDimension, NDimension,
2952 resultJ);
2953 T expected;
2954 for (int i = 0; i < NDimension; ++i)
2955 {
2956 for (int j = 0; j < NDimension; ++j)
2957 {
2958 if (i == j)
2959 {
2960 expected = eigen[i];
2961 }
2962 else
2963 {
2964 expected = 0.0;
2965 }
2966 if (!vtkMathUtilities::FuzzyCompare(
2967 result[i][j], expected,
2968 std::numeric_limits<T>::epsilon() * (T) 256.0))
2969 {
2970 std::cout << " JacobiN expected " << expected
2971 << " but got " << result[i][j];
2972 std::cout << "eps ratio is: " << (result[i][j] - expected)
2973 / std::numeric_limits<T>::epsilon() << std::endl;
2974 ++status;
2975 }
2976 }
2977 }
2978 }
2979 return status;
2980 }
2981
2982 // Validate P^-1 * A * P = diagonal matrix with eigenvalues on diagonal
TestJacobiN()2983 int TestJacobiN()
2984 {
2985 int status = 0;
2986 std::cout << "JacobiN..";
2987 status += JacobiN<double,3>();
2988 status += JacobiN<double,10>();
2989 status += JacobiN<double,50>();
2990
2991 if (status)
2992 {
2993 std::cout << "..FAILED" << std::endl;
2994 }
2995 else
2996 {
2997 std::cout << ".PASSED" << std::endl;
2998 }
2999 return status;
3000 }
3001
3002 template<typename T>
RGBToHSV()3003 int RGBToHSV()
3004 {
3005 int status = 0;
3006 T R, G, B;
3007 T H, S, V;
3008 T CR, CG, CB;
3009 for (int n = 0; n < 1000; ++n)
3010 {
3011 std::vector<T> values;
3012 std::vector<T> expecteds;
3013 R = (T) vtkMath::Random(0.0, 1.0);
3014 G = (T) vtkMath::Random(0.0, 1.0);
3015 B = (T) vtkMath::Random(0.0, 1.0);
3016
3017 vtkMath::RGBToHSV(R, G, B, &H, &S, &V);
3018 vtkMath::HSVToRGB(H, S, V, &CR, &CG, &CB);
3019 values.push_back(CR); values.push_back(CG); values.push_back(CB);
3020 expecteds.push_back(R); expecteds.push_back(G); expecteds.push_back(B);
3021
3022 for ( size_t i = 0; i < values.size(); ++i)
3023 {
3024 if (!vtkMathUtilities::FuzzyCompare(
3025 values[i], expecteds[i],
3026 std::numeric_limits<T>::epsilon() * (T) 128.0))
3027 {
3028 std::cout << " RGBToHSV got " << values[i]
3029 << " but expected " << expecteds[i];
3030 std::cout << " eps ratio is: " << (values[i] - expecteds[i])
3031 / std::numeric_limits<T>::epsilon() << std::endl;
3032 ++status;
3033 }
3034 }
3035 }
3036 return status;
3037 }
3038
3039 // Validate by rgb->hsv->rgb
TestRGBToHSV()3040 int TestRGBToHSV()
3041 {
3042 int status = 0;
3043 std::cout << "RGBToHSV..";
3044
3045 status += RGBToHSV<double>();
3046 status += RGBToHSV<float>();
3047 if (status)
3048 {
3049 std::cout << "..FAILED" << std::endl;
3050 }
3051 else
3052 {
3053 std::cout << ".PASSED" << std::endl;
3054 }
3055 return status;
3056 }
3057
3058 // Validate with known solutions
TestClampValue()3059 int TestClampValue()
3060 {
3061 int status = 0;
3062 std::cout << "ClampValue..";
3063
3064 double value;
3065 double clampedValue;
3066 double range[2] = {-1.0, 1.0};
3067
3068 value = -800.0;
3069 clampedValue = vtkMath::ClampValue(value, range[0], range[1]);
3070 if (clampedValue != range[0])
3071 {
3072 std::cout << " ClampValue expected " << range[0]
3073 << " but got " << value;
3074 ++status;
3075 }
3076
3077 value = 900.0;
3078 clampedValue = vtkMath::ClampValue(value, range[0], range[1]);
3079 if (clampedValue != range[1])
3080 {
3081 std::cout << " ClampValue expected " << range[1]
3082 << " but got " << value;
3083 ++status;
3084 }
3085
3086 value = 0.0;
3087 clampedValue = vtkMath::ClampValue(value, range[0], range[1]);
3088 if (clampedValue != 0.0)
3089 {
3090 std::cout << " ClampValue expected " << 0.0
3091 << " but got " << value;
3092 ++status;
3093 }
3094
3095
3096 value = -100.0;
3097 vtkMath::ClampValue(&value, range);
3098 if (value != range[0])
3099 {
3100 std::cout << " ClampValue expected " << range[0]
3101 << " but got " << value;
3102 ++status;
3103 }
3104 value = 100.0;
3105 vtkMath::ClampValue(&value, range);
3106 if (value != range[1])
3107 {
3108 std::cout << " ClampValue expected " << range[1]
3109 << " but got " << value;
3110 ++status;
3111 }
3112 value = -100.0;
3113 vtkMath::ClampValue(value, range, &clampedValue);
3114 if (clampedValue != range[0])
3115 {
3116 std::cout << " ClampValue expected " << range[0]
3117 << " but got " << clampedValue;
3118 ++status;
3119 }
3120
3121 value = 100.0;
3122 vtkMath::ClampValue(value, range, &clampedValue);
3123 if (clampedValue != range[1])
3124 {
3125 std::cout << " ClampValue expected " << range[1]
3126 << " but got " << clampedValue;
3127 ++status;
3128 }
3129
3130 value = 0.0;
3131 vtkMath::ClampValue(value, range, &clampedValue);
3132 if (clampedValue != value)
3133 {
3134 std::cout << " ClampValue expected " << value
3135 << " but got " << clampedValue;
3136 ++status;
3137 }
3138
3139 if (status)
3140 {
3141 std::cout << "..FAILED" << std::endl;
3142 }
3143 else
3144 {
3145 std::cout << ".PASSED" << std::endl;
3146 }
3147 return status;
3148 }
3149
3150 // Validate with known solutions
TestClampValues()3151 int TestClampValues()
3152 {
3153 int status = 0;
3154 std::cout << "ClampValues..";
3155
3156 double values[1000];
3157 double clampedValues[1000];
3158 for (int n = 0; n < 1000; ++n)
3159 {
3160 values[n] = vtkMath::Random(-2.0, 2.0);
3161 }
3162 double range[2] = {-1.0, 1.0};
3163 vtkMath::ClampValues( values, 1000, range, clampedValues);
3164 vtkMath::ClampValues( values, 1000, range);
3165
3166 for (int n = 0; n < 1000; ++n)
3167 {
3168 if (values[n] != clampedValues[n])
3169 {
3170 ++status;
3171 }
3172 }
3173
3174 vtkMath::ClampValues( nullptr, 1000, nullptr);
3175 vtkMath::ClampValues( nullptr, 1000, nullptr, nullptr);
3176
3177 if (status)
3178 {
3179 std::cout << "..FAILED" << std::endl;
3180 }
3181 else
3182 {
3183 std::cout << ".PASSED" << std::endl;
3184 }
3185 return status;
3186 }
3187
3188 // Validate with known solutions
TestClampAndNormalizeValue()3189 int TestClampAndNormalizeValue()
3190 {
3191 int status = 0;
3192 std::cout << "ClampAndNormalizeValue..";
3193
3194 double value;
3195 double result;
3196 double range[3] = {-1.0, 1.0};
3197
3198 value = -100.0;
3199 result = vtkMath::ClampAndNormalizeValue(value, range);
3200 if (result != 0.0)
3201 {
3202 std::cout << " ClampAndNormalizeValue expected " << 0.0
3203 << " but got " << result;
3204 ++status;
3205 }
3206 value = 100.0;
3207 result = vtkMath::ClampAndNormalizeValue(value, range);
3208 if (result != 1.0)
3209 {
3210 std::cout << " ClampAndNormalizeValue expected " << 1.0
3211 << " but got " << result;
3212 ++status;
3213 }
3214
3215 range[0] = 0.0; range[1] = 1.0;
3216 value = 0.5;
3217 result = vtkMath::ClampAndNormalizeValue(value, range);
3218 if (result != 0.5)
3219 {
3220 std::cout << " ClampValue expected " << 0.5
3221 << " but got " << result;
3222 ++status;
3223 }
3224
3225 range[0] = 1.0; range[1] = 1.0;
3226 value = 1.0;
3227 result = vtkMath::ClampAndNormalizeValue(value, range);
3228 if (result != 0.0)
3229 {
3230 std::cout << " ClampValue expected " << 0.0
3231 << " but got " << result;
3232 ++status;
3233 }
3234
3235 if (status)
3236 {
3237 std::cout << "..FAILED" << std::endl;
3238 }
3239 else
3240 {
3241 std::cout << ".PASSED" << std::endl;
3242 }
3243 return status;
3244 }
3245
3246 // Validate by checking symmetric tensor values
3247 // are in correct places
TestTensorFromSymmetricTensor()3248 int TestTensorFromSymmetricTensor()
3249 {
3250 int status = 0;
3251 std::cout << "TensorFromSymmetricTensor..";
3252 double symmTensor[9];
3253 for (int i = 0; i < 6; i++)
3254 {
3255 symmTensor[i] = vtkMath::Random();
3256 }
3257 double tensor[9];
3258 vtkMath::TensorFromSymmetricTensor(symmTensor, tensor);
3259 if (tensor[0] != symmTensor[0] ||
3260 tensor[1] != symmTensor[3] ||
3261 tensor[2] != symmTensor[5] ||
3262 tensor[3] != symmTensor[3] ||
3263 tensor[4] != symmTensor[1] ||
3264 tensor[5] != symmTensor[4] ||
3265 tensor[6] != symmTensor[5] ||
3266 tensor[7] != symmTensor[4] ||
3267 tensor[8] != symmTensor[2])
3268 {
3269 std::cout << " Unexpected results from TensorFromSymmetricTensor " << std::endl;
3270 ++status;
3271 }
3272
3273 vtkMath::TensorFromSymmetricTensor(symmTensor);
3274 for (int i = 0; i < 9; i++)
3275 {
3276 if (symmTensor[i] != tensor[i])
3277 {
3278 std::cout << " Unexpected results from in place TensorFromSymmetricTensor " << std::endl;
3279 ++status;
3280 break;
3281 }
3282 }
3283 if (status)
3284 {
3285 std::cout << "..FAILED" << std::endl;
3286 }
3287 else
3288 {
3289 std::cout << ".PASSED" << std::endl;
3290 }
3291 return status;
3292 }
3293
3294 // Validate by checking ranges with numeric_limits
TestGetScalarTypeFittingRange()3295 int TestGetScalarTypeFittingRange()
3296 {
3297 int status = 0;
3298 std::cout << "GetScalarTypeFittingRange..";
3299
3300 double rangeMin;
3301 double rangeMax;
3302
3303 rangeMin = (double) std::numeric_limits<char>::min();
3304 rangeMax = (double) std::numeric_limits<char>::max();
3305 if (vtkMath::GetScalarTypeFittingRange(
3306 rangeMin, rangeMax, 1.0, 0.0) != VTK_CHAR)
3307 {
3308 std::cout << " Bad fitting range for VTK_CHAR" << std::endl;
3309 ++status;
3310 }
3311
3312 rangeMin = (double) std::numeric_limits<unsigned char>::min();
3313 rangeMax = (double) std::numeric_limits<unsigned char>::max();
3314 if (vtkMath::GetScalarTypeFittingRange(
3315 rangeMin, rangeMax, 1.0, 0.0) != VTK_UNSIGNED_CHAR)
3316 {
3317 std::cout << " Bad fitting range for VTK_UNSIGNED_CHAR " << std::endl;
3318 ++status;
3319 }
3320
3321 rangeMin = (double) std::numeric_limits<short>::min();
3322 rangeMax = (double) std::numeric_limits<short>::max();
3323 if (vtkMath::GetScalarTypeFittingRange(
3324 rangeMin, rangeMax, 1.0, 0.0) != VTK_SHORT)
3325 {
3326 std::cout << " Bad fitting range for VTK_SHORT" << std::endl;
3327 ++status;
3328 }
3329
3330 rangeMin = (double) std::numeric_limits<unsigned short>::min();
3331 rangeMax = (double) std::numeric_limits<unsigned short>::max();
3332 if (vtkMath::GetScalarTypeFittingRange(
3333 rangeMin, rangeMax, 1.0, 0.0) != VTK_UNSIGNED_SHORT)
3334 {
3335 std::cout << " Bad fitting range for VTK_UNSIGNED_SHORT" << std::endl;
3336 ++status;
3337 }
3338
3339 rangeMin = (double) std::numeric_limits<int>::min();
3340 rangeMax = (double) std::numeric_limits<int>::max();
3341 if (vtkMath::GetScalarTypeFittingRange(
3342 rangeMin, rangeMax, 1.0, 0.0) != VTK_INT)
3343 {
3344 std::cout << " Bad fitting range for VTK_INT" << std::endl;
3345 ++status;
3346 }
3347
3348 rangeMin = (double) std::numeric_limits<unsigned int>::min();
3349 rangeMax = (double) std::numeric_limits<unsigned int>::max();
3350 if (vtkMath::GetScalarTypeFittingRange(
3351 rangeMin, rangeMax, 1.0, 0.0) != VTK_UNSIGNED_INT)
3352 {
3353 std::cout << " Bad fitting range for VTK_UNSIGNED_INT" << std::endl;
3354 ++status;
3355 }
3356
3357 rangeMin = (double) std::numeric_limits<long>::min();
3358 rangeMax = (double) std::numeric_limits<long>::max();
3359 int scalarType = vtkMath::GetScalarTypeFittingRange(
3360 rangeMin, rangeMax, 1.0, 0.0);
3361 if (sizeof(long) == sizeof(int))
3362 {
3363 if (scalarType != VTK_INT)
3364 {
3365 std::cout << " Bad fitting range for VTK_LONG" << std::endl;
3366 std::cout << " Expected " << VTK_INT << " but got " << scalarType;
3367 ++status;
3368 }
3369 }
3370 else
3371 {
3372 if (scalarType != VTK_LONG)
3373 {
3374 std::cout << " Bad fitting range for VTK_LONG" << std::endl;
3375 std::cout << " Expected " << VTK_LONG << " but got " << scalarType;
3376 ++status;
3377 }
3378 }
3379
3380 rangeMin = (double) std::numeric_limits<unsigned long>::min();
3381 rangeMax = (double) std::numeric_limits<unsigned long>::max();
3382 scalarType = vtkMath::GetScalarTypeFittingRange(
3383 rangeMin, rangeMax, 1.0, 0.0);
3384 if (sizeof(unsigned long) == sizeof(unsigned int))
3385 {
3386 if (scalarType != VTK_UNSIGNED_INT)
3387 {
3388 std::cout << " Bad fitting range for VTK_UNSIGNED_LONG" << std::endl;
3389 std::cout << " Expected " << VTK_UNSIGNED_INT << " but got " << scalarType;
3390 ++status;
3391 }
3392 }
3393 else
3394 {
3395 if (scalarType != VTK_UNSIGNED_LONG)
3396 {
3397 std::cout << " Bad fitting range for VTK_UNSIGNED_LONG" << std::endl;
3398 std::cout << " Expected " << VTK_UNSIGNED_LONG << " but got " << scalarType;
3399 ++status;
3400 }
3401 }
3402
3403 rangeMin = (double) std::numeric_limits<short>::min();
3404 rangeMax = (double) std::numeric_limits<short>::max();
3405 if (vtkMath::GetScalarTypeFittingRange(
3406 rangeMin, rangeMax, 1.1, 0.0) != VTK_FLOAT)
3407 {
3408 std::cout << " Bad fitting range for VTK_FLOAT" << std::endl;
3409 ++status;
3410 }
3411
3412 if (status)
3413 {
3414 std::cout << "..FAILED" << std::endl;
3415 }
3416 else
3417 {
3418 std::cout << ".PASSED" << std::endl;
3419 }
3420 return status;
3421 }
3422
3423 // Validate with known solutions
TestGetAdjustedScalarRange()3424 int TestGetAdjustedScalarRange()
3425 {
3426 int status = 0;
3427 std::cout << "GetAdjustedScalarRange..";
3428
3429 vtkSmartPointer<vtkUnsignedCharArray> uc =
3430 vtkSmartPointer<vtkUnsignedCharArray>::New();
3431 uc->SetNumberOfComponents(3);
3432 uc->SetNumberOfTuples(100);
3433 for (int i = 0; i < 100; ++i)
3434 {
3435 for (int j = 0; j < 3; ++j)
3436 {
3437 uc->SetComponent(i, j, vtkMath::Random(
3438 std::numeric_limits<unsigned char>::min(),
3439 std::numeric_limits<unsigned char>::max()));
3440 }
3441 }
3442 double range[2];
3443 vtkMath::GetAdjustedScalarRange(uc, 1, range);
3444 if (range[0] != uc->GetDataTypeMin() ||
3445 range[1] != uc->GetDataTypeMax())
3446 {
3447 std::cout << " GetAdjustedScalarRange(unsigned char) expected "
3448 << uc->GetDataTypeMin() << ", " << uc->GetDataTypeMax()
3449 << " but got " << range[0] << ", " << range[1]
3450 << std::endl;
3451 ++status;
3452 }
3453
3454 vtkSmartPointer<vtkUnsignedShortArray> us =
3455 vtkSmartPointer<vtkUnsignedShortArray>::New();
3456 us->SetNumberOfComponents(3);
3457 us->SetNumberOfTuples(10000);
3458 for (int i = 0; i < 10000; ++i)
3459 {
3460 us->SetComponent(i, 0, vtkMath::Random(
3461 std::numeric_limits<unsigned short>::min(),
3462 std::numeric_limits<unsigned short>::max()));
3463 us->SetComponent(i, 1, vtkMath::Random(
3464 std::numeric_limits<unsigned short>::min(),
3465 std::numeric_limits<unsigned char>::max() + 100));
3466 us->SetComponent(i, 2, vtkMath::Random(
3467 std::numeric_limits<unsigned short>::min(),
3468 std::numeric_limits<unsigned char>::max()));
3469 }
3470 vtkMath::GetAdjustedScalarRange(us, 0, range);
3471 if (range[0] != us->GetDataTypeMin() ||
3472 range[1] != us->GetDataTypeMax())
3473 {
3474 std::cout << " GetAdjustedScalarRange(unsigned short) expected "
3475 << us->GetDataTypeMin() << ", " << us->GetDataTypeMax()
3476 << " but got " << range[0] << ", " << range[1]
3477 << std::endl;
3478 ++status;
3479 }
3480
3481 vtkMath::GetAdjustedScalarRange(us, 1, range);
3482 if (range[0] != us->GetDataTypeMin() ||
3483 range[1] != 4095.0)
3484 {
3485 std::cout << " GetAdjustedScalarRange(unsigned short) expected "
3486 << us->GetDataTypeMin() << ", " << 4095.0
3487 << " but got " << range[0] << ", " << range[1]
3488 << std::endl;
3489 ++status;
3490 }
3491
3492 vtkMath::GetAdjustedScalarRange(us, 2, range);
3493 if (range[0] != us->GetDataTypeMin() ||
3494 range[1] >= uc->GetDataTypeMax())
3495 {
3496 std::cout << " GetAdjustedScalarRange(unsigned short) expected "
3497 << us->GetDataTypeMin() << ", " << ">= " << uc->GetDataTypeMax()
3498 << " but got " << range[0] << ", " << range[1]
3499 << std::endl;
3500 ++status;
3501 }
3502
3503 // Test nullptr array
3504 if (vtkMath::GetAdjustedScalarRange(nullptr, 1000, nullptr))
3505 {
3506 std::cout << " GetAdjustedScalarRange with a nullptr array expected "
3507 << 0
3508 << " but got " << 1
3509 << std::endl;
3510 ++status;
3511 }
3512 if (status)
3513 {
3514 std::cout << "..FAILED" << std::endl;
3515 }
3516 else
3517 {
3518 std::cout << ".PASSED" << std::endl;
3519 }
3520 return status;
3521 }
3522
3523 // Validate with known solutions
TestExtentIsWithinOtherExtent()3524 int TestExtentIsWithinOtherExtent()
3525 {
3526 int status = 0;
3527 std::cout << "ExtentIsWithinOtherExtent..";
3528
3529 if (vtkMath::ExtentIsWithinOtherExtent(nullptr, nullptr))
3530 {
3531 std::cout << " ExtentIsWithinOtherExtent expected 0 but got 1"
3532 << std::endl;
3533 ++status;
3534 }
3535
3536 int extent1[6];
3537 int extent2[6];
3538 extent1[0] = 100; extent1[1] = 101;
3539 extent1[2] = 100; extent1[3] = 101;
3540 extent1[4] = 100; extent1[5] = 101;
3541
3542 extent2[0] = 100; extent2[1] = 101;
3543 extent2[2] = 100; extent2[3] = 101;
3544 extent2[4] = 100; extent2[5] = 101;
3545
3546 if (!vtkMath::ExtentIsWithinOtherExtent(extent1, extent2))
3547 {
3548 std::cout << " ExtentIsWithinOtherExtent expected 1 but got 0"
3549 << std::endl;
3550 ++status;
3551 }
3552
3553 extent1[0] = 99; extent1[1] = 101;
3554 if (vtkMath::ExtentIsWithinOtherExtent(extent1, extent2))
3555 {
3556 std::cout << " ExtentIsWithinOtherExtent expected 0 but got 1"
3557 << std::endl;
3558 ++status;
3559 }
3560
3561 extent1[0] = 98; extent1[1] = 99;
3562 if (vtkMath::ExtentIsWithinOtherExtent(extent1, extent2))
3563 {
3564 std::cout << " ExtentIsWithinOtherExtent expected 0 but got 1"
3565 << std::endl;
3566 ++status;
3567 }
3568
3569 if (status)
3570 {
3571 std::cout << "..FAILED" << std::endl;
3572 }
3573 else
3574 {
3575 std::cout << ".PASSED" << std::endl;
3576 }
3577 return status;
3578 }
3579
3580 // Validate with known solutions
TestBoundsIsWithinOtherBounds()3581 int TestBoundsIsWithinOtherBounds()
3582 {
3583 int status = 0;
3584 std::cout << "BoundsIsWithinOtherBounds..";
3585
3586 if (vtkMath::BoundsIsWithinOtherBounds(nullptr, nullptr, nullptr))
3587 {
3588 std::cout << " BoundsIsWithinOtherBounds expected 0 but got 1"
3589 << std::endl;
3590 ++status;
3591 }
3592
3593 double delta[3];
3594 delta[0] = delta[1] = delta[2] = std::numeric_limits<double>::epsilon();
3595
3596 double bounds1[6];
3597 double bounds2[6];
3598 bounds1[0] = 1 - delta[0]; bounds1[1] = 2 + delta[0];
3599 bounds1[2] = 1; bounds1[3] = 2;
3600 bounds1[4] = 1; bounds1[5] = 2;
3601
3602 bounds2[0] = 1; bounds2[1] = 2;
3603 bounds2[2] = 1; bounds2[3] = 2;
3604 bounds2[4] = 1; bounds2[5] = 2;
3605
3606 if (!vtkMath::BoundsIsWithinOtherBounds(bounds1, bounds2, delta))
3607 {
3608 std::cout << " BoundsIsWithinOtherBounds expected 1 but got 0"
3609 << std::endl;
3610 ++status;
3611 }
3612
3613 bounds1[0] = 1 - 2.0 * delta[0]; bounds1[1] = 2 + 2.0 * delta[0];
3614 if (vtkMath::BoundsIsWithinOtherBounds(bounds1, bounds2, delta))
3615 {
3616 std::cout << " BoundsIsWithinOtherBounds expected 0 but got 1"
3617 << std::endl;
3618 ++status;
3619 }
3620
3621 bounds1[0] = 1 - 4.0 * delta[0]; bounds1[1] = 1 - 2.0 * delta[0];
3622 if (vtkMath::BoundsIsWithinOtherBounds(bounds1, bounds2, delta))
3623 {
3624 std::cout << " BoundsIsWithinOtherBounds expected 0 but got 1"
3625 << std::endl;
3626 ++status;
3627 }
3628
3629 if (status)
3630 {
3631 std::cout << "..FAILED" << std::endl;
3632 }
3633 else
3634 {
3635 std::cout << ".PASSED" << std::endl;
3636 }
3637 return status;
3638 }
3639
3640 // Validate with known solutions
TestPointIsWithinBounds()3641 int TestPointIsWithinBounds()
3642 {
3643 int status = 0;
3644 std::cout << "PointIsWithinBounds..";
3645
3646 if (vtkMath::PointIsWithinBounds(nullptr, nullptr, nullptr))
3647 {
3648 std::cout << " PointIsWithinBounds expected 0 but got 1"
3649 << std::endl;
3650 ++status;
3651 }
3652
3653 double delta[3];
3654 delta[0] = std::numeric_limits<double>::epsilon();
3655 delta[1] = std::numeric_limits<double>::epsilon() * 2.0;
3656 delta[2] = std::numeric_limits<double>::epsilon() * 256.0;
3657
3658 double bounds1[6];
3659 bounds1[0] = 1.0; bounds1[1] = 2.0;
3660 bounds1[2] = 1.0; bounds1[3] = 2.0;
3661 bounds1[4] = 1.0; bounds1[5] = 2.0;
3662
3663 double point[3];
3664 point[0] = bounds1[0] - delta[0];
3665 point[1] = bounds1[2] - delta[1];
3666 point[2] = bounds1[4];
3667
3668 if (!vtkMath::PointIsWithinBounds(point, bounds1, delta))
3669 {
3670 std::cout << " PointIsWithinBounds expected 1 but got 0"
3671 << std::endl;
3672 ++status;
3673 }
3674
3675 point[0] = bounds1[0] - delta[0];
3676 point[1] = bounds1[2] - delta[1];
3677 point[2] = bounds1[4] - 2.0 * delta[2];
3678
3679 if (vtkMath::PointIsWithinBounds(point, bounds1, delta))
3680 {
3681 std::cout << " PointIsWithinBounds expected 0 but got 1"
3682 << std::endl;
3683 ++status;
3684 }
3685
3686 point[0] = bounds1[1] + delta[0];
3687 point[1] = bounds1[3] + delta[1];
3688 point[2] = bounds1[5] + 2.0 * delta[2];
3689
3690 if (vtkMath::PointIsWithinBounds(point, bounds1, delta))
3691 {
3692 std::cout << " PointIsOtherBounds expected 0 but got 1"
3693 << std::endl;
3694 ++status;
3695 }
3696
3697 if (status)
3698 {
3699 std::cout << "..FAILED" << std::endl;
3700 }
3701 else
3702 {
3703 std::cout << ".PASSED" << std::endl;
3704 }
3705 return status;
3706 }
3707
3708 // Validate with with alternative solution
TestSolve3PointCircle()3709 int TestSolve3PointCircle()
3710 {
3711 int status = 0;
3712 std::cout << "Solve3PointCircle..";
3713
3714 for (int n = 0; n < 1000; ++n)
3715 {
3716 double A[3], B[3], C[3];
3717 double center[3];
3718 double a[3], b[3], aMinusb[3], aCrossb[3];
3719
3720 for (int i = 0; i < 3; ++i)
3721 {
3722 A[i] = vtkMath::Random(-1.0, 1.0);
3723 B[i] = vtkMath::Random(-1.0, 1.0);
3724 C[i] = vtkMath::Random(-1.0, 1.0);
3725 }
3726
3727 vtkMath::Subtract(A, C, a);
3728 vtkMath::Subtract(B, C, b);
3729 vtkMath::Subtract(a, b, aMinusb);
3730 vtkMath::Cross(a, b, aCrossb);
3731
3732 double expectedRadius;
3733 expectedRadius =
3734 (vtkMath::Norm(a) * vtkMath::Norm(b) * vtkMath::Norm(aMinusb)) /
3735 (2.0 * vtkMath::Norm(aCrossb));
3736
3737 double radius;
3738 radius = vtkMath::Solve3PointCircle(A, B, C, center);
3739 if (!vtkMathUtilities::FuzzyCompare(
3740 radius, expectedRadius,
3741 std::numeric_limits<double>::epsilon()*1024.0))
3742 {
3743 std::cout << " Solve3PointCircle radius expected " << expectedRadius
3744 << " but got " << radius;
3745 std::cout << "eps ratio is: " << (expectedRadius - radius)
3746 / std::numeric_limits<double>::epsilon() << std::endl;
3747 ++status;
3748 }
3749
3750 double ab[3], ba[3];
3751 double abMinusba[3];
3752 double abMinusbaCrossaCrossb[3];
3753
3754 vtkMath::Subtract(B, C, ab);
3755 vtkMath::Subtract(A, C, ba);
3756 vtkMath::MultiplyScalar(ab, vtkMath::Norm(a) * vtkMath::Norm(a));
3757 vtkMath::MultiplyScalar(ba, vtkMath::Norm(b) * vtkMath::Norm(b));
3758 vtkMath::Subtract(ab, ba, abMinusba);
3759 vtkMath::Cross(abMinusba, aCrossb, abMinusbaCrossaCrossb);
3760
3761 double expectedCenter[3];
3762 vtkMath::MultiplyScalar(
3763 abMinusbaCrossaCrossb,
3764 1.0 / (2.0 * vtkMath::Norm(aCrossb) * vtkMath::Norm(aCrossb)));
3765 vtkMath::Add(abMinusbaCrossaCrossb, C, expectedCenter);
3766
3767 for (int i = 0; i < 3; ++i)
3768 {
3769 if (!vtkMathUtilities::FuzzyCompare(
3770 center[i], expectedCenter[i],
3771 std::numeric_limits<double>::epsilon()*1024.0))
3772 {
3773 std::cout << " Solve3PointCircle center expected " << expectedCenter[i]
3774 << " but got " << center[i];
3775 std::cout << "eps ratio is: " << (expectedCenter[i] - center[i])
3776 / std::numeric_limits<double>::epsilon() << std::endl;
3777 ++status;
3778 }
3779 }
3780 }
3781
3782 if (status)
3783 {
3784 std::cout << "..FAILED" << std::endl;
3785 }
3786 else
3787 {
3788 std::cout << ".PASSED" << std::endl;
3789 }
3790 return status;
3791 }
3792
3793 // Tested by TestMath
TestInf()3794 int TestInf()
3795 {
3796 int status = 0;
3797 std::cout << "Inf..";
3798
3799 if (status)
3800 {
3801 std::cout << "..FAILED" << std::endl;
3802 }
3803 else
3804 {
3805 std::cout << ".PASSED" << std::endl;
3806 }
3807 return status;
3808 }
3809
3810 // Tested by TestMath
TestNegInf()3811 int TestNegInf()
3812 {
3813 int status = 0;
3814 std::cout << "NegInf..";
3815
3816 if (status)
3817 {
3818 std::cout << "..FAILED" << std::endl;
3819 }
3820 else
3821 {
3822 std::cout << ".PASSED" << std::endl;
3823 }
3824 return status;
3825 }
3826
3827 // Tested by TestMath
TestNan()3828 int TestNan()
3829 {
3830 int status = 0;
3831 std::cout << "Nan..";
3832
3833 if (status)
3834 {
3835 std::cout << "..FAILED" << std::endl;
3836 }
3837 else
3838 {
3839 std::cout << ".PASSED" << std::endl;
3840 }
3841 return status;
3842 }
3843