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