GCC Code Coverage Report

Directory: ./
File: openvdb/openvdb/math/Maps.cc
Date: 2022-07-25 17:40:05
Exec Total Coverage
Lines: 67 79 84.8%
Functions: 10 12 83.3%
Branches: 60 132 45.5%
Line Branch Exec Source
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3
4 #include "Maps.h"
5 #include <mutex>
6
7 namespace openvdb {
8 OPENVDB_USE_VERSION_NAMESPACE
9 namespace OPENVDB_VERSION_NAME {
10 namespace math {
11
12 namespace {
13
14 // Declare this at file scope to ensure thread-safe initialization.
15 // NOTE: Do *NOT* move this into Maps.h or else we will need to pull in
16 // Windows.h with things like 'rad2' defined!
17 std::mutex sInitMapRegistryMutex;
18
19 } // unnamed namespace
20
21
22 ////////////////////////////////////////
23
24
25 // Caller is responsible for calling this function serially.
26 MapRegistry*
27 9324 MapRegistry::staticInstance()
28 {
29
3/4
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 9322 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
 9324 static MapRegistry registry;
30 9324 return &registry;
31 }
32
33
34 MapRegistry*
35 ✗ MapRegistry::instance()
36 {
37 std::lock_guard<std::mutex> lock(sInitMapRegistryMutex);
38 ✗ return staticInstance();
39 }
40
41
42 MapBase::Ptr
43 137 MapRegistry::createMap(const Name& name)
44 {
45 std::lock_guard<std::mutex> lock(sInitMapRegistryMutex);
46
2/4
✓ Branch 1 taken 137 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 137 times.
✗ Branch 5 not taken.
 274 MapDictionary::const_iterator iter = staticInstance()->mMap.find(name);
47
48
2/4
✓ Branch 1 taken 137 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 137 times.
 137 if (iter == staticInstance()->mMap.end()) {
49 ✗ OPENVDB_THROW(LookupError, "Cannot create map of unregistered type " << name);
50 }
51
52
1/2
✓ Branch 1 taken 137 times.
✗ Branch 2 not taken.
 274 return (iter->second)();
53 }
54
55
56 bool
57 157 MapRegistry::isRegistered(const Name& name)
58 {
59 std::lock_guard<std::mutex> lock(sInitMapRegistryMutex);
60
3/6
✓ Branch 1 taken 157 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 157 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 157 times.
✗ Branch 7 not taken.
 471 return (staticInstance()->mMap.find(name) != staticInstance()->mMap.end());
61 }
62
63
64 void
65 2681 MapRegistry::registerMap(const Name& name, MapBase::MapFactory factory)
66 {
67 std::lock_guard<std::mutex> lock(sInitMapRegistryMutex);
68
69
3/6
✓ Branch 1 taken 2681 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2681 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✓ Branch 7 taken 2681 times.
 5362 if (staticInstance()->mMap.find(name) != staticInstance()->mMap.end()) {
70 ✗ OPENVDB_THROW(KeyError, "Map type " << name << " is already registered");
71 }
72
73
3/6
✓ Branch 1 taken 2681 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 2681 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 2681 times.
✗ Branch 7 not taken.
 2681 staticInstance()->mMap[name] = factory;
74 2681 }
75
76
77 void
78 ✗ MapRegistry::unregisterMap(const Name& name)
79 {
80 std::lock_guard<std::mutex> lock(sInitMapRegistryMutex);
81 ✗ staticInstance()->mMap.erase(name);
82 }
83
84
85 void
86 693 MapRegistry::clear()
87 {
88 std::lock_guard<std::mutex> lock(sInitMapRegistryMutex);
89
1/2
✓ Branch 1 taken 693 times.
✗ Branch 2 not taken.
 693 staticInstance()->mMap.clear();
90 693 }
91
92
93 ////////////////////////////////////////
94
95 // Utility methods for decomposition
96
97
98 SymmetricMap::Ptr
99 1 createSymmetricMap(const Mat3d& m)
100 {
101 // test that the mat3 is a rotation || reflection
102
1/2
✗ Branch 1 not taken.
✓ Branch 2 taken 1 times.
 1 if (!isSymmetric(m)) {
103 ✗ OPENVDB_THROW(ArithmeticError,
104 "3x3 Matrix initializing symmetric map was not symmetric");
105 }
106 Vec3d eigenValues;
107 Mat3d Umatrix;
108
109 1 bool converged = math::diagonalizeSymmetricMatrix(m, Umatrix, eigenValues);
110
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (!converged) {
111 ✗ OPENVDB_THROW(ArithmeticError, "Diagonalization of the symmetric matrix failed");
112 }
113
114 1 UnitaryMap rotation(Umatrix);
115 1 ScaleMap diagonal(eigenValues);
116 1 CompoundMap<UnitaryMap, ScaleMap> first(rotation, diagonal);
117
118 1 UnitaryMap rotationInv(Umatrix.transpose());
119
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 2 return SymmetricMap::Ptr( new SymmetricMap(first, rotationInv));
120 }
121
122
123 PolarDecomposedMap::Ptr
124 1 createPolarDecomposedMap(const Mat3d& m)
125 {
126 // Because our internal libary left-multiplies vectors against matrices
127 // we are constructing M = Symmetric * Unitary instead of the more
128 // standard M = Unitary * Symmetric
129 Mat3d unitary, symmetric, mat3 = m.transpose();
130
131 // factor mat3 = U * S where U is unitary and S is symmetric
132 1 bool gotPolar = math::polarDecomposition(mat3, unitary, symmetric);
133
1/2
✗ Branch 0 not taken.
✓ Branch 1 taken 1 times.
 1 if (!gotPolar) {
134 ✗ OPENVDB_THROW(ArithmeticError, "Polar decomposition of transform failed");
135 }
136 // put the result in a polar map and then copy it into the output polar
137 1 UnitaryMap unitary_map(unitary.transpose());
138 1 SymmetricMap::Ptr symmetric_map = createSymmetricMap(symmetric);
139
140
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 return PolarDecomposedMap::Ptr(new PolarDecomposedMap(*symmetric_map, unitary_map));
141 }
142
143
144 FullyDecomposedMap::Ptr
145
1/2
✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
 1 createFullyDecomposedMap(const Mat4d& m)
146 {
147 if (!isAffine(m)) {
148 ✗ OPENVDB_THROW(ArithmeticError,
149 "4x4 Matrix initializing Decomposition map was not affine");
150 }
151
152 1 TranslationMap translate(m.getTranslation());
153
1/2
✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
 1 PolarDecomposedMap::Ptr polar = createPolarDecomposedMap(m.getMat3());
154
155
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 UnitaryAndTranslationMap rotationAndTranslate(polar->secondMap(), translate);
156
157
2/4
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
 2 return FullyDecomposedMap::Ptr(new FullyDecomposedMap(polar->firstMap(), rotationAndTranslate));
158 }
159
160
161 MapBase::Ptr
162
2/2
✓ Branch 0 taken 21 times.
✓ Branch 1 taken 101 times.
 122 simplify(AffineMap::Ptr affine)
163 {
164
2/2
✓ Branch 0 taken 21 times.
✓ Branch 1 taken 101 times.
 122 if (affine->isScale()) { // can be simplified into a ScaleMap
165
166
2/2
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 9 times.
 21 Vec3d scale = affine->applyMap(Vec3d(1,1,1));
167
4/4
✓ Branch 0 taken 12 times.
✓ Branch 1 taken 9 times.
✓ Branch 2 taken 11 times.
✓ Branch 3 taken 1 times.
 21 if (isApproxEqual(scale[0], scale[1]) && isApproxEqual(scale[0], scale[2])) {
168 11 return MapBase::Ptr(new UniformScaleMap(scale[0]));
169 } else {
170
1/2
✓ Branch 2 taken 10 times.
✗ Branch 3 not taken.
 10 return MapBase::Ptr(new ScaleMap(scale));
171 }
172
173
2/2
✓ Branch 1 taken 3 times.
✓ Branch 2 taken 98 times.
 101 } else if (affine->isScaleTranslate()) { // can be simplified into a ScaleTranslateMap
174
175 3 Vec3d translate = affine->applyMap(Vec3d(0,0,0));
176
1/2
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
 3 Vec3d scale = affine->applyMap(Vec3d(1,1,1)) - translate;
177
178
2/4
✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✓ Branch 3 taken 3 times.
 3 if (isApproxEqual(scale[0], scale[1]) && isApproxEqual(scale[0], scale[2])) {
179 3 return MapBase::Ptr(new UniformScaleTranslateMap(scale[0], translate));
180 } else {
181 ✗ return MapBase::Ptr(new ScaleTranslateMap(scale, translate));
182 }
183 }
184
185 // could not simplify the general Affine map.
186 return StaticPtrCast<MapBase, AffineMap>(affine);
187 }
188
189
190 Mat4d
191 4 approxInverse(const Mat4d& mat4d)
192 {
193
2/2
✓ Branch 1 taken 1 times.
✓ Branch 2 taken 3 times.
 4 if (std::abs(mat4d.det()) >= 3 * math::Tolerance<double>::value()) {
194 try {
195
1/2
✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
 1 return mat4d.inverse();
196 ✗ } catch (ArithmeticError& ) {
197 // Mat4 code couldn't invert.
198 }
199 }
200
201 const Mat3d mat3 = mat4d.getMat3();
202 const Mat3d mat3T = mat3.transpose();
203 const Vec3d trans = mat4d.getTranslation();
204
205 // absolute tolerance used for the symmetric test.
206 const double tol = 1.e-6;
207
208 // only create the pseudoInverse for symmetric
209 bool symmetric = true;
210
2/2
✓ Branch 0 taken 9 times.
✓ Branch 1 taken 3 times.
 12 for (int i = 0; i < 3; ++i ) {
211
2/2
✓ Branch 0 taken 27 times.
✓ Branch 1 taken 9 times.
 36 for (int j = 0; j < 3; ++j ) {
212
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 25 times.
 27 if (!isApproxEqual(mat3[i][j], mat3T[i][j], tol)) {
213 symmetric = false;
214 }
215 }
216 }
217
218
2/2
✓ Branch 0 taken 1 times.
✓ Branch 1 taken 2 times.
 3 if (!symmetric) {
219
220 // not symmetric, so just zero out the mat3 inverse and reverse the translation
221
222 1 Mat4d result = Mat4d::zero();
223 result.setTranslation(-trans);
224 1 result[3][3] = 1.f;
225 1 return result;
226
227 } else {
228
229 // compute the pseudo inverse
230
231 Mat3d eigenVectors;
232 Vec3d eigenValues;
233
234 2 diagonalizeSymmetricMatrix(mat3, eigenVectors, eigenValues);
235
236 2 Mat3d d = Mat3d::identity();
237
2/2
✓ Branch 0 taken 6 times.
✓ Branch 1 taken 2 times.
 8 for (int i = 0; i < 3; ++i ) {
238
2/2
✓ Branch 0 taken 2 times.
✓ Branch 1 taken 4 times.
 6 if (std::abs(eigenValues[i]) < 10.0 * math::Tolerance<double>::value()) {
239 2 d[i][i] = 0.f;
240 } else {
241 4 d[i][i] = 1.f/eigenValues[i];
242 }
243 }
244 // assemble the pseudo inverse
245
246 2 Mat3d pseudoInv = eigenVectors * d * eigenVectors.transpose();
247 2 Vec3d invTrans = -trans * pseudoInv;
248
249 2 Mat4d result = Mat4d::identity();
250 result.setMat3(pseudoInv);
251 result.setTranslation(invTrans);
252
253 2 return result;
254 }
255 }
256
257 } // namespace math
258 } // namespace OPENVDB_VERSION_NAME
259 } // namespace openvdb
260