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