Refactor SVD transforms

etc/adjusted/synthetic.json

@@ -238,7 +238,7 @@
                 1
             ]
         ],
-        "RotationTranslation3D": [
+        "SVD": [
             [
                 0.8100811050998874,
                 -0.3030939677547838,

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(

geomagio/adjusted/transform/RotationTranslation3D.py

-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],
-        ]

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],
         ]

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],
+        ]

geomagio/adjusted/transform/__init__.py

 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",

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,
     )