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