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