From fbf4cfbe8b925d7adccaa2ed5c44f52641d48996 Mon Sep 17 00:00:00 2001
From: pcain-usgs <pcain@usgs.gov>
Date: Tue, 9 Feb 2021 16:40:28 -0700
Subject: [PATCH] Refactor SVD transforms
---
etc/adjusted/synthetic.json | 2 +-
.../adjusted/transform/QRFactorization.py | 4 +-
.../transform/RotationTranslation3D.py | 57 ---------------
.../transform/RotationTranslationXY.py | 46 ++----------
geomagio/adjusted/transform/SVD.py | 71 ++++++++++++++++++-
geomagio/adjusted/transform/__init__.py | 3 +-
test/adjusted_test/adjusted_test.py | 8 +--
7 files changed, 84 insertions(+), 107 deletions(-)
delete mode 100644 geomagio/adjusted/transform/RotationTranslation3D.py
diff --git a/etc/adjusted/synthetic.json b/etc/adjusted/synthetic.json
index 015f75cd1..29f466f5c 100644
--- a/etc/adjusted/synthetic.json
+++ b/etc/adjusted/synthetic.json
@@ -238,7 +238,7 @@
1
]
],
- "RotationTranslation3D": [
+ "SVD": [
[
0.8100811050998874,
-0.3030939677547838,
diff --git a/geomagio/adjusted/transform/QRFactorization.py b/geomagio/adjusted/transform/QRFactorization.py
index 866744926..fa8187382 100644
--- a/geomagio/adjusted/transform/QRFactorization.py
+++ b/geomagio/adjusted/transform/QRFactorization.py
@@ -8,6 +8,8 @@ from .SVD import SVD
class QRFactorization(SVD):
"""Calculates affine using singular value decomposition with QR factorization"""
+ ndims = 2
+
def get_weighted_values_lstsq(
self,
values: Tuple[List[float], List[float], List[float]],
@@ -40,14 +42,12 @@ class QRFactorization(SVD):
abs_stacked = self.get_stacked_values(
values=absolutes,
weighted_values=weighted_absolutes,
- ndims=2,
)
# RHS, or independent variables
ord_stacked = self.get_stacked_values(
values=ordinates,
weighted_values=weighted_ordinates,
- ndims=2,
)
abs_stacked = self.get_weighted_values_lstsq(
diff --git a/geomagio/adjusted/transform/RotationTranslation3D.py b/geomagio/adjusted/transform/RotationTranslation3D.py
deleted file mode 100644
index d609b4f53..000000000
--- a/geomagio/adjusted/transform/RotationTranslation3D.py
+++ /dev/null
@@ -1,57 +0,0 @@
-import numpy as np
-import scipy.linalg as spl
-from typing import List, Tuple
-
-from .SVD import SVD
-
-
-class RotationTranslation3D(SVD):
- """Calculates affine using using singular value decomposition,
- constrained to 3D scaled rotation and translation(no shear)"""
-
- def calculate(
- self,
- ordinates: Tuple[List[float], List[float], List[float]],
- absolutes: Tuple[List[float], List[float], List[float]],
- weights: List[float],
- ) -> np.array:
- weighted_ordinates = self.get_weighted_values(values=ordinates, weights=weights)
- weighted_absolutes = self.get_weighted_values(values=absolutes, weights=weights)
- # generate weighted "covariance" matrix
- H = np.dot(
- self.get_stacked_values(
- values=ordinates,
- weighted_values=weighted_ordinates,
- ndims=3,
- ),
- np.dot(
- np.diag(weights),
- self.get_stacked_values(
- values=absolutes,
- weighted_values=weighted_absolutes,
- ndims=3,
- ).T,
- ),
- )
- # Singular value decomposition, then rotation matrix from L&R eigenvectors
- # (the determinant guarantees a rotation, and not a reflection)
- U, S, Vh = np.linalg.svd(H)
- if np.sum(S) < 3:
- print("Poorly conditioned or singular matrix, returning NaNs")
- return np.nan * np.ones((4, 4))
-
- R = np.dot(Vh.T, np.dot(np.diag([1, 1, np.linalg.det(np.dot(Vh.T, U.T))]), U.T))
-
- # now get translation using weighted centroids and R
- T = np.array(
- [weighted_absolutes[0], weighted_absolutes[1], weighted_absolutes[2]]
- ) - np.dot(
- R, [weighted_ordinates[0], weighted_ordinates[1], weighted_ordinates[2]]
- )
-
- return [
- [R[0, 0], R[0, 1], R[0, 2], T[0]],
- [R[1, 0], R[1, 1], R[1, 2], T[1]],
- [R[2, 0], R[2, 1], R[2, 2], T[2]],
- [0.0, 0.0, 0.0, 1.0],
- ]
diff --git a/geomagio/adjusted/transform/RotationTranslationXY.py b/geomagio/adjusted/transform/RotationTranslationXY.py
index b0cf427d9..ade36e195 100644
--- a/geomagio/adjusted/transform/RotationTranslationXY.py
+++ b/geomagio/adjusted/transform/RotationTranslationXY.py
@@ -9,47 +9,12 @@ class RotationTranslationXY(SVD):
constrained to rotation and translation in XY(no scale or shear),
and only translation in Z"""
- def calculate(
- self,
- ordinates: Tuple[List[float], List[float], List[float]],
- absolutes: Tuple[List[float], List[float], List[float]],
- weights: List[float],
- ) -> np.array:
- if weights is None:
- weights = np.ones_like(ordinates[0])
- weighted_ordinates = self.get_weighted_values(values=ordinates, weights=weights)
- weighted_absolutes = self.get_weighted_values(values=absolutes, weights=weights)
- # generate weighted "covariance" matrix
- H = np.dot(
- self.get_stacked_values(
- values=ordinates,
- weighted_values=weighted_ordinates,
- ndims=2,
- ),
- np.dot(
- np.diag(weights),
- self.get_stacked_values(
- values=absolutes,
- weighted_values=weighted_absolutes,
- ndims=2,
- ).T,
- ),
- )
+ ndims = 2
- # Singular value decomposition, then rotation matrix from L&R eigenvectors
- # (the determinant guarantees a rotation, and not a reflection)
- U, S, Vh = np.linalg.svd(H)
- if np.sum(S) < 2:
- print("Poorly conditioned or singular matrix, returning NaNs")
- return np.nan * np.ones((4, 4))
-
- R = np.dot(Vh.T, np.dot(np.diag([1, np.linalg.det(np.dot(Vh.T, U.T))]), U.T))
-
- # now get translation using weighted centroids and R
- T = np.array([weighted_absolutes[0], weighted_absolutes[1]]) - np.dot(
- R, [weighted_ordinates[0], weighted_ordinates[1]]
- )
+ def get_rotation_matrix(self, U, Vh):
+ return np.dot(Vh.T, np.dot(np.diag([1, np.linalg.det(np.dot(Vh.T, U.T))]), U.T))
+ def format_matrix(self, R, T, weighted_absolutes, weighted_ordinates):
return [
[R[0, 0], R[0, 1], 0.0, T[0]],
[R[1, 0], R[1, 1], 0.0, T[1]],
@@ -57,7 +22,8 @@ class RotationTranslationXY(SVD):
0.0,
0.0,
1.0,
- np.array(weighted_absolutes[2]) - np.array(weighted_ordinates[2]),
+ np.array(weighted_absolutes[self.ndims])
+ - np.array(weighted_ordinates[self.ndims]),
],
[0.0, 0.0, 0.0, 1.0],
]
diff --git a/geomagio/adjusted/transform/SVD.py b/geomagio/adjusted/transform/SVD.py
index 201dbc417..d1b8505dc 100644
--- a/geomagio/adjusted/transform/SVD.py
+++ b/geomagio/adjusted/transform/SVD.py
@@ -7,7 +7,28 @@ from .Transform import Transform
class SVD(Transform):
"""Instance of Transform. Applies singular value decomposition to generate matrices"""
- def get_stacked_values(self, values, weighted_values, ndims=3) -> np.array:
+ ndims = 3
+
+ def get_covariance_matrix(self, absolutes, ordinates, weights):
+ weighted_ordinates = self.get_weighted_values(values=ordinates, weights=weights)
+ weighted_absolutes = self.get_weighted_values(values=absolutes, weights=weights)
+ # generate weighted "covariance" matrix
+ H = np.dot(
+ self.get_stacked_values(
+ values=ordinates,
+ weighted_values=weighted_ordinates,
+ ),
+ np.dot(
+ np.diag(weights),
+ self.get_stacked_values(
+ values=absolutes,
+ weighted_values=weighted_absolutes,
+ ).T,
+ ),
+ )
+ return H
+
+ def get_stacked_values(self, values, weighted_values) -> np.array:
"""Supports intermediate mathematical steps by differencing and shaping values for SVD
Attributes
@@ -20,7 +41,7 @@ class SVD(Transform):
-------
Stacked and differenced values from their weighted counterparts
"""
- return np.vstack([[values[i] - weighted_values[i]] for i in range(ndims)])
+ return np.vstack([[values[i] - weighted_values[i]] for i in range(self.ndims)])
def get_weighted_values(
self,
@@ -35,3 +56,49 @@ class SVD(Transform):
np.average(values[1], weights=weights),
np.average(values[2], weights=weights),
)
+
+ def get_rotation_matrix(self, U, Vh):
+ return np.dot(
+ Vh.T, np.dot(np.diag([1, 1, np.linalg.det(np.dot(Vh.T, U.T))]), U.T)
+ )
+
+ def get_translation_matrix(self, R, weighted_absolutes, weighted_ordinates):
+ return np.array([weighted_absolutes[i] for i in range(self.ndims)]) - np.dot(
+ R, [weighted_ordinates[i] for i in range(self.ndims)]
+ )
+
+ def calculate(
+ self,
+ ordinates: Tuple[List[float], List[float], List[float]],
+ absolutes: Tuple[List[float], List[float], List[float]],
+ weights: List[float],
+ ) -> np.array:
+ weighted_ordinates = self.get_weighted_values(values=ordinates, weights=weights)
+ weighted_absolutes = self.get_weighted_values(values=absolutes, weights=weights)
+ # generate weighted "covariance" matrix
+ H = self.get_covariance_matrix(absolutes, ordinates, weights)
+ # Singular value decomposition, then rotation matrix from L&R eigenvectors
+ # (the determinant guarantees a rotation, and not a reflection)
+ U, S, Vh = np.linalg.svd(H)
+ if np.sum(S) < 3:
+ print("Poorly conditioned or singular matrix, returning NaNs")
+ return np.nan * np.ones((4, 4))
+ R = self.get_rotation_matrix(U, Vh)
+ # now get translation using weighted centroids and R
+ T = self.get_translation_matrix(R, weighted_absolutes, weighted_ordinates)
+
+ return self.format_matrix(R, T, weighted_absolutes, weighted_ordinates)
+
+ def format_matrix(
+ self,
+ R,
+ T,
+ weighted_absolutes,
+ weighted_ordinates,
+ ):
+ return [
+ [R[0, 0], R[0, 1], R[0, 2], T[0]],
+ [R[1, 0], R[1, 1], R[1, 2], T[1]],
+ [R[2, 0], R[2, 1], R[2, 2], T[2]],
+ [0.0, 0.0, 0.0, 1.0],
+ ]
diff --git a/geomagio/adjusted/transform/__init__.py b/geomagio/adjusted/transform/__init__.py
index 85156ca47..0fa99ab9a 100644
--- a/geomagio/adjusted/transform/__init__.py
+++ b/geomagio/adjusted/transform/__init__.py
@@ -1,11 +1,11 @@
from .LeastSq import LeastSq
from .QRFactorization import QRFactorization
from .Rescale3D import Rescale3D
-from .RotationTranslation3D import RotationTranslation3D
from .RotationTranslationXY import RotationTranslationXY
from .ShearYZ import ShearYZ
from .Transform import Transform
from .TranslateOrigins import TranslateOrigins
+from .SVD import SVD
from .ZRotationHScale import ZRotationHscale
from .ZRotationHScaleZBaseline import ZRotationHscaleZbaseline
from .ZRotationShear import ZRotationShear
@@ -19,6 +19,7 @@ __all__ = [
"ShearYZ",
"Transform",
"TranslateOrigins",
+ "SVD",
"ZRotationHscale",
"ZRotationHscaleZbaseline",
"ZRotationShear",
diff --git a/test/adjusted_test/adjusted_test.py b/test/adjusted_test/adjusted_test.py
index 276e626fb..6e4ef6368 100644
--- a/test/adjusted_test/adjusted_test.py
+++ b/test/adjusted_test/adjusted_test.py
@@ -10,7 +10,7 @@ from geomagio.adjusted.transform import (
LeastSq,
QRFactorization,
Rescale3D,
- RotationTranslation3D,
+ SVD,
RotationTranslationXY,
ShearYZ,
TranslateOrigins,
@@ -79,15 +79,15 @@ def test_Rescale3D_synthetic():
)
-def test_RotationTranslation3D_synthetic():
+def test_SVD_synthetic():
ordinates, absolutes, weights = get_sythetic_variables()
assert_array_almost_equal(
- RotationTranslation3D().calculate(
+ SVD().calculate(
ordinates=ordinates,
absolutes=absolutes,
weights=weights,
),
- get_expected_synthetic_result("RotationTranslation3D"),
+ get_expected_synthetic_result("SVD"),
decimal=3,
)
--
GitLab