Exterior is also a numpy array now

Author: alexlib

openptv_python/calibration.py

@@ -6,20 +6,35 @@
 import numpy as np
 from numba import njit
 
-from openptv_python.vec_utils import vec_set
-
 
 @njit
-def rotation_matrix(phi: float, omega: float, kappa: float) -> np.ndarray:
-    """Calculate the necessary trigonometric functions to rotate the Dmatrix of Exterior ext_par."""
-    cp = np.cos(phi)
-    sp = np.sin(phi)
+def rotation_matrix(ext: np.ndarray) -> None:
+    """Calculate the necessary trigonometric functions to rotate the Dmatrix of Exterior ext_par.
+
+    Rotation is performed by multiplication of three rotation matrices,
+    rotation around X axis is performed first, then around Y axis, and finally around Z axis.
+    rotation around X axis does not change values along X axis, i.e. X_omega = X
+
+    Maas, H.G., Gruen, A. & Papantoniou, D. Particle tracking velocimetry in
+    three-dimensional flows. Experiments in Fluids 15, 133–146 (1993).
+    https://doi.org/10.1007/BF00190953
+
+    """
+    omega, phi, kappa = ext[0]['omega'], ext[0]['phi'], ext[0]['kappa']
+
+
     co = np.cos(omega)
     so = np.sin(omega)
+
+    cp = np.cos(phi)
+    sp = np.sin(phi)
+
     ck = np.cos(kappa)
     sk = np.sin(kappa)
 
-    dm = np.zeros((3, 3), dtype=np.float64)
+    # dm = np.zeros((3, 3), dtype=np.float64)
+    dm = ext[0]['dm'] # shortcut to the dm field of the first element of the array
+
     dm[0, 0] = cp * ck
     dm[0, 1] = -cp * sk
     dm[0, 2] = sp
@@ -29,59 +44,39 @@ def rotation_matrix(phi: float, omega: float, kappa: float) -> np.ndarray:
     dm[2, 0] = so * sk - co * sp * ck
     dm[2, 1] = so * ck + co * sp * sk
     dm[2, 2] = co * cp
-    return dm
 
+    # print(dm)
 
-class Exterior:
-    """Exterior orientation data structure."""
+    return None
 
-    def __init__(self, x0=0.0, y0=0.0, z0=0.0, omega=0.0, phi=0.0, kappa=0.0, dm=None):
-        self.x0 = x0
-        self.y0 = y0
-        self.z0 = z0
-        self.omega = omega
-        self.phi = phi
-        self.kappa = kappa
-        self.dm = dm if dm is not None else np.identity(3, dtype=np.float64)
+exterior_dtype = np.dtype([
+    ('x0', np.float64),
+    ('y0', np.float64),
+    ('z0', np.float64),
+    ('omega', np.float64),
+    ('phi', np.float64),
+    ('kappa', np.float64),
+    ('dm', np.float64, (3, 3))
+    ])
+Exterior = np.zeros(1, dtype=exterior_dtype).view(np.recarray) # initialize memory
+rotation_matrix(Exterior)             # rotation should be a unit matrix
+assert np.allclose(np.eye(3), Exterior[0]['dm'])
 
-    def update_rotation_matrix(self) -> None:
-        """Create rotation matrix using three angles of the camera."""
-        self.dm = rotation_matrix(self.phi, self.omega, self.kappa)
-
-    def set_rotation_matrix(self, dm: np.ndarray) -> None:
-        """Set the rotation matrix of the camera."""
-        self.dm = dm
-
-    def set_pos(self, pos: np.ndarray) -> None:
-        """Set the position of the camera."""
-        pos = np.array(pos, dtype = np.float64)
-
-        if pos.shape != (3,):
-            raise ValueError(
-                "Illegal array argument "
-                + str(pos)
-                + " for x, y, z. Expected array/list of 3 numbers"
-            )
-        self.x0, self.y0, self.z0 = pos
 
-    def set_angles(self, angles: np.ndarray) -> None:
-        """Set the angles of the camera."""
-        self.omega, self.phi, self.kappa = angles
+# rotation_matrix(exterior) # inplace update of the rotation matrix
 
-        # adjust rotation matrix
-        self.update_rotation_matrix()
 
-    def increment_attribute(self, attr_name, increment_value):
-        """Update the value of an attribute by increment_value."""
-        if hasattr(self, attr_name):
-            setattr(self, attr_name, getattr(
-                self, attr_name) + increment_value)
+    # def increment_attribute(self, attr_name, increment_value):
+    #     """Update the value of an attribute by increment_value."""
+    #     if hasattr(self, attr_name):
+    #         setattr(self, attr_name, getattr(
+    #             self, attr_name) + increment_value)
 
-    def __repr__(self) -> str:
-        """Return a string representation of the Exterior object."""
-        output = f"Exterior: x0={self.x0}, y0={self.y0}, z0={self.z0}\n"
-        output += f"omega={self.omega}, phi={self.phi}, kappa={self.kappa}\n"
-        return output
+    # def __repr__(self) -> str:
+    #     """Return a string representation of the Exterior object."""
+    #     output = f"Exterior: x0={self['x0']}, y0={self['y0']}, z0={self['z0']}\n"
+    #     output += f"omega={self['omega']}, phi={self['phi']}, kappa={self['kappa']}\n"
+    #     return output
 
 
 class Interior:
@@ -100,10 +95,6 @@ def set_back_focal_distance(self, cc: float) -> None:
         """Set the back focal distance of the camera."""
         self.cc = cc
 
-def default_glass_vec() -> np.ndarray:
-    """Return default glass vector."""
-    return vec_set(0.0, 0.0, 1.0)
-
 class ap_52:
     """Additional parameters for distortion correction."""
 
@@ -149,7 +140,7 @@ class Calibration:
 
     def __init__(self, ext_par=None, int_par=None, glass_par=None, added_par=None, mmlut=None):
         if ext_par is None:
-            ext_par = Exterior()
+            ext_par = Exterior.copy()
         if int_par is None:
             int_par = Interior()
         if glass_par is None:
@@ -165,6 +156,7 @@ def __init__(self, ext_par=None, int_par=None, glass_par=None, added_par=None, m
         self.added_par = added_par
         self.mmlut = mmlut
 
+
     @classmethod
     def from_file(cls, ori_file: str, add_file: str):
         """
@@ -197,7 +189,7 @@ def from_file(cls, ori_file: str, add_file: str):
             # skip line
             fp.readline()
 
-            # read 3 lines and set a rotation matrix
+            # read 3 lines of rotation matrix, but recalculate it from angles
             _ = [[float(x) for x in fp.readline().split()] for _ in range(3)]
 
             # I skip reading rotation matrix as it's set up by set_angles.
@@ -254,23 +246,42 @@ def write(self, ori_file: str, addpar_file: str):
         if not success:
             raise IOError("Failed to write ori file")
 
+
+    def increment_attribute(self, attr_name, increment_value):
+        """Update the value of an attribute by increment_value."""
+        if hasattr(self, attr_name):
+            setattr(self, attr_name, getattr(
+                self, attr_name) + increment_value)
+
+    def update_rotation_matrix(self) -> None:
+        """Update the rotation matrix of the exterior orientation."""
+        rotation_matrix(self.ext_par)
+
     def set_rotation_matrix(self, dm: np.ndarray) -> None:
         """Set exterior rotation matrix."""
         if dm.shape != (3, 3):
             raise ValueError("Illegal argument for exterior rotation matrix")
-        self.ext_par.set_rotation_matrix(dm)
+        self.ext_par[0]['dm'] = dm
 
-    def set_pos(self, x_y_z_np: np.ndarray) -> None:
+    def set_pos(self, pos: np.ndarray) -> None:
         """
         Set exterior position.
 
         Parameter: x_y_z_np - numpy array of 3 elements for x, y, z.
         """
-        self.ext_par.set_pos(x_y_z_np)
+        pos = np.array(pos, dtype = np.float64)
 
-    def get_pos(self):
+        if pos.shape != (3,):
+            raise ValueError(
+                "Illegal array argument "
+                + str(pos)
+                + " for x, y, z. Expected array/list of 3 numbers"
+            )
+        self.ext_par['x0'], self.ext_par['y0'], self.ext_par['z0'] = pos
+
+    def get_pos(self) -> np.ndarray:
         """Return array of 3 elements representing exterior's x, y, z."""
-        return np.r_[self.ext_par.x0, self.ext_par.y0, self.ext_par.z0]
+        return np.r_[self.ext_par['x0'], self.ext_par['y0'], self.ext_par['z0']]
 
     def set_angles(self, o_p_k_np: np.ndarray) -> None:
         """
@@ -285,15 +296,16 @@ def set_angles(self, o_p_k_np: np.ndarray) -> None:
                 f"Illegal array argument {o_p_k_np} for "
                 "omega, phi, kappa. Expected array or list of 3 float"
             )
-        self.ext_par.set_angles(o_p_k_np)
+        self.ext_par['omega'], self.ext_par['phi'], self.ext_par['kappa'] = o_p_k_np
+        self.update_rotation_matrix()
 
-    def get_angles(self):
+    def get_angles(self) -> np.ndarray:
         """Return an array of 3 elements representing omega, phi, kappa."""
-        return np.r_[self.ext_par.omega, self.ext_par.phi, self.ext_par.kappa]
+        return np.r_[self.ext_par['omega'], self.ext_par['phi'], self.ext_par['kappa']]
 
-    def get_rotation_matrix(self):
+    def get_rotation_matrix(self) -> np.ndarray:
         """Return a 3x3 numpy array that represents Exterior's rotation matrix."""
-        return self.ext_par.dm
+        return self.ext_par['dm'].copy()
 
     def set_primary_point(self, prim_point_pos: np.ndarray) -> None:
         """
@@ -418,7 +430,7 @@ def copy(self, new_copy):
 
 
 def write_ori(
-    ext_par: Exterior,
+    ext_par: np.ndarray,
     int_par: Interior,
     glass: np.ndarray,
     added_par: ap_52,
@@ -429,9 +441,9 @@ def write_ori(
     success = False
 
     with open(filename, "w", encoding="utf-8") as fp:
-        fp.write(f"{ext_par.x0:.8f} {ext_par.y0:.8f} {ext_par.z0:.8f}\n")
-        fp.write(f"{ext_par.omega:.8f} {ext_par.phi:.8f} {ext_par.kappa:.8f}\n\n")
-        for row in ext_par.dm:
+        fp.write(f"{ext_par[0]['x0']:.8f} {ext_par[0]['y0']:.8f} {ext_par[0]['z0']:.8f}\n")
+        fp.write(f"{ext_par[0]['omega']:.8f} {ext_par[0]['phi']:.8f} {ext_par[0]['kappa']:.8f}\n\n")
+        for row in ext_par[0]['dm']:
             fp.write(f"{row[0]:.7f} {row[1]:.7f} {row[2]:.7f}\n")
         fp.write(f"\n{int_par.xh:.4f} {int_par.yh:.4f}\n{int_par.cc:.4f}\n")
         fp.write(
@@ -470,16 +482,16 @@ def read_ori(ori_file: str, add_file: str) -> Calibration:
     return ret
 
 
-def compare_exterior(e1: Exterior, e2: Exterior) -> bool:
+def compare_exterior(e1: np.ndarray, e2: np.ndarray) -> bool:
     """Compare exterior orientation parameters."""
     return (
-        np.allclose(e1.dm, e2.dm, atol=1e-6)
-        and (e1.x0 == e2.x0)
-        and (e1.y0 == e2.y0)
-        and (e1.z0 == e2.z0)
-        and (e1.omega == e2.omega)
-        and (e1.phi == e2.phi)
-        and (e1.kappa == e2.kappa)
+        np.allclose(e1['dm'], e2['dm'], atol=1e-6)
+        and (e1['x0'] == e2['x0'])
+        and (e1['y0'] == e2['y0'])
+        and (e1['z0'] == e2['z0'])
+        and (e1['omega'] == e2['omega'])
+        and (e1['phi'] == e2['phi'])
+        and (e1['kappa'] == e2['kappa'])
     )
 
 

openptv_python/imgcoord.py

@@ -9,7 +9,6 @@
 from .multimed import back_trans_point, multimed_nlay, trans_cam_point
 from .parameters import MultimediaPar
 from .trafo import flat_to_dist
-from .vec_utils import vec_set
 
 
 def flat_image_coord(
@@ -46,43 +45,18 @@ def flat_image_coord(
     x_t, y_t = multimed_nlay(cal_t, mm, pos_t)
     # print(f"x_t {x_t}, y_t {y_t}")
 
-    pos_t = vec_set(x_t, y_t, pos_t[2])
+    pos_t = np.r_[x_t, y_t, pos_t[2]]
     pos = back_trans_point(pos_t, mm, cal.glass_par, cross_p, cross_c)
 
-    # print(f"pos {pos}")
-
-    deno = (
-        cal.ext_par.dm[0][2] * (pos[0] - cal.ext_par.x0)
-        + cal.ext_par.dm[1][2] * (pos[1] - cal.ext_par.y0)
-        + cal.ext_par.dm[2][2] * (pos[2] - cal.ext_par.z0)
-    )
-
-    # print(f"deno {deno}")
+    dm = cal.ext_par[0].dm
+    origin = np.r_[cal.ext_par.x0, cal.ext_par.y0, cal.ext_par.z0]
 
+    deno = np.dot(dm[:,2], (pos - origin))
     if deno == 0:
         deno = 1
 
-    x = (
-        -cal.int_par.cc
-        * (
-            cal.ext_par.dm[0][0] * (pos[0] - cal.ext_par.x0)
-            + cal.ext_par.dm[1][0] * (pos[1] - cal.ext_par.y0)
-            + cal.ext_par.dm[2][0] * (pos[2] - cal.ext_par.z0)
-        )
-        / deno
-    )
-
-    y = (
-        -cal.int_par.cc
-        * (
-            cal.ext_par.dm[0][1] * (pos[0] - cal.ext_par.x0)
-            + cal.ext_par.dm[1][1] * (pos[1] - cal.ext_par.y0)
-            + cal.ext_par.dm[2][1] * (pos[2] - cal.ext_par.z0)
-        )
-        / deno
-    )
-
-    # print(f"x {x}, y {y}")
+    x = (-cal.int_par.cc * np.dot(dm[:, 0], (pos - origin)) / deno)
+    y = (-cal.int_par.cc * np.dot(dm[:, 1], (pos - origin)) / deno)
 
     return x, y
 

openptv_python/multimed.py

@@ -4,7 +4,7 @@
 import numpy as np
 from numba import njit
 
-from .calibration import Calibration, Exterior
+from .calibration import Calibration
 from .parameters import (
     ControlPar,
     MultimediaPar,
@@ -116,7 +116,7 @@ def fast_multimed_r_nlay(
 
 
 def trans_cam_point(
-    ex: Exterior, mm: MultimediaPar, glass_dir: np.ndarray, pos: np.ndarray
+    ex: np.ndarray, mm: MultimediaPar, glass_dir: np.ndarray, pos: np.ndarray
 ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, float]:
     """Transform the camera and point coordinates to the glass coordinates.
 
@@ -127,7 +127,7 @@ def trans_cam_point(
 
     pos_t, cross_p, cross_c = trans_cam_point(ex, mm, glass, pos, ex_t)
     """
-    origin = np.array([ex.x0, ex.y0, ex.z0], dtype=np.float64)
+    origin = np.r_[ex.x0, ex.y0, ex.z0]
     pos = pos.astype(np.float64)
 
     return fast_trans_cam_point(
@@ -196,7 +196,7 @@ def back_trans_point(
     return fast_back_trans_point(glass, mm.d[0], cross_c, cross_p, pos_t)
 
 @njit
-def fast_back_trans_point(glass_direction: np.recarray, d: float, cross_c, cross_p, pos_t) -> np.ndarray:
+def fast_back_trans_point(glass_direction: np.ndarray, d: float, cross_c, cross_p, pos_t) -> np.ndarray:
     """Run numba faster version of back projection."""
     # Calculate the glass direction vector