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