exterior is a numpy array. works well with recarray view

Author: alexlib

openptv_python/calibration.py

@@ -20,7 +20,7 @@ def rotation_matrix(ext: np.ndarray) -> None:
     https://doi.org/10.1007/BF00190953
 
     """
-    omega, phi, kappa = ext[0]['omega'], ext[0]['phi'], ext[0]['kappa']
+    omega, phi, kappa = ext['omega'], ext['phi'], ext['kappa']
 
 
     co = np.cos(omega)
@@ -33,7 +33,7 @@ def rotation_matrix(ext: np.ndarray) -> None:
     sk = np.sin(kappa)
 
     # 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 = ext['dm'] # shortcut to the dm field of the first element of the array
 
     dm[0, 0] = cp * ck
     dm[0, 1] = -cp * sk
@@ -58,9 +58,10 @@ def rotation_matrix(ext: np.ndarray) -> None:
     ('kappa', np.float64),
     ('dm', np.float64, (3, 3))
     ])
-Exterior = np.zeros(1, dtype=exterior_dtype).view(np.recarray) # initialize memory
+# Exterior = np.zeros(1, dtype=exterior_dtype).view(np.recarray) # initialize memory
+Exterior = np.array((0, 0, 0, 0, 0, 0, np.eye(3)), dtype = exterior_dtype).view(np.recarray)
 rotation_matrix(Exterior)             # rotation should be a unit matrix
-assert np.allclose(np.eye(3), Exterior[0]['dm'])
+assert np.allclose(np.eye(3), Exterior['dm'])
 
 
 # rotation_matrix(exterior) # inplace update of the rotation matrix
@@ -249,9 +250,12 @@ def write(self, ori_file: str, addpar_file: str):
 
     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)
+        if hasattr(self.ext_par, attr_name):
+            setattr(self.ext_par, attr_name, getattr(
+                self.ext_par, attr_name) + increment_value)
+        if hasattr(self.int_par, attr_name):
+            setattr(self.int_par, attr_name, getattr(
+                self.int_par, attr_name) + increment_value)
 
     def update_rotation_matrix(self) -> None:
         """Update the rotation matrix of the exterior orientation."""
@@ -441,9 +445,9 @@ def write_ori(
     success = False
 
     with open(filename, "w", encoding="utf-8") as fp:
-        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"{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"{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(

openptv_python/imgcoord.py

@@ -48,7 +48,7 @@ def flat_image_coord(
     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)
 
-    dm = cal.ext_par[0].dm
+    dm = cal.ext_par.dm
     origin = np.r_[cal.ext_par.x0, cal.ext_par.y0, cal.ext_par.z0]
 
     deno = np.dot(dm[:,2], (pos - origin))

openptv_python/multimed.py

@@ -1,4 +1,3 @@
-import math
 from typing import List, Tuple
 
 import numpy as np
@@ -12,7 +11,7 @@
 )
 from .ray_tracing import ray_tracing
 from .trafo import correct_brown_affine, pixel_to_metric
-from .vec_utils import norm, vec_set
+from .vec_utils import vec_norm
 
 
 def multimed_nlay(
@@ -64,9 +63,9 @@ def fast_multimed_r_nlay(
     n2: np.ndarray,
     n3: float,
     d: np.ndarray,
-    x0,
-    y0,
-    z0,
+    x0: float,
+    y0: float,
+    z0: float,
     pos: np.ndarray
 ) -> float:
     """Calculate the radial shift for the multimedia model.
@@ -87,22 +86,23 @@ def fast_multimed_r_nlay(
     for i in range(1, nlay):
         zout += d[i]
 
-    r = norm(X-x0, Y-y0, 0)
+
+    r = vec_norm(np.array([X-x0, Y-y0, 0]))
     rq = r
 
     it = 0
     while (rdiff > 0.001 or rdiff < -0.001) and it < n_iter:
         zdiff = z0 - Z
         if zdiff == 0:
             zdiff = 1.0
-        beta1 = math.atan(rq / zdiff)
+        beta1 = np.arctan(rq / zdiff)
         for i in range(nlay):
-            beta2[i] = math.asin(math.sin(beta1) * n1 / n2[i])
-        beta3 = math.asin(math.sin(beta1) * n1 / n3)
+            beta2[i] = np.arcsin(np.sin(beta1) * n1 / n2[i])
+        beta3 = np.arcsin(np.sin(beta1) * n1 / n3)
 
-        rbeta = (z0 - d[0]) * math.tan(beta1) - zout * math.tan(beta3)
+        rbeta = (z0 - d[0]) * np.tan(beta1) - zout * np.tan(beta3)
         for i in range(nlay):
-            rbeta += d[i] * math.tan(beta2[i])
+            rbeta += d[i] * np.tan(beta2[i])
 
         rdiff = r - rbeta
         rq += rdiff
@@ -127,14 +127,14 @@ def trans_cam_point(
 
     pos_t, cross_p, cross_c = trans_cam_point(ex, mm, glass, pos, ex_t)
     """
-    origin = np.r_[ex.x0, ex.y0, ex.z0]
+    origin = np.r_[ex.x0, ex.y0, ex.z0] # type: ignore
     pos = pos.astype(np.float64)
 
     return fast_trans_cam_point(
         origin, mm.d[0], glass_dir, pos)
 
 
-# @njit(fastmath=True)
+@njit(fastmath=True)
 def fast_trans_cam_point(
     primary_point: np.ndarray,
     d: float,
@@ -266,8 +266,8 @@ def init_mmlut(vpar: VolumePar, cpar: ControlPar, cal: Calibration) -> Calibrati
     z_min = min(vpar.z_min_lay)
     z_max = max(vpar.z_max_lay)
 
-    z_min -= math.fmod(z_min, rw)
-    z_max += rw - math.fmod(z_max, rw)
+    z_min -= np.fmod(z_min, rw)
+    z_max += rw - np.fmod(z_max, rw)
 
     z_min_t = z_min
     z_max_t = z_max
@@ -292,8 +292,11 @@ def init_mmlut(vpar: VolumePar, cpar: ControlPar, cal: Calibration) -> Calibrati
             if xyz_t[2] > z_max_t:
                 z_max_t = xyz_t[2]
 
-            R = norm(xyz_t[0] - cal_t.ext_par.x0,
-                     xyz_t[1] - cal_t.ext_par.y0, 0)
+            R = vec_norm(
+                np.r_[xyz_t[0] - cal_t.ext_par.x0,
+                     xyz_t[1] - cal_t.ext_par.y0,
+                     0]
+                )
 
             if R > Rmax:
                 Rmax = R
@@ -308,14 +311,14 @@ def init_mmlut(vpar: VolumePar, cpar: ControlPar, cal: Calibration) -> Calibrati
             if xyz_t[2] > z_max_t:
                 z_max_t = xyz_t[2]
 
-            R = norm(xyz_t[0] - cal_t.ext_par.x0,
-                     xyz_t[1] - cal_t.ext_par.y0, 0)
+            R = vec_norm(np.r_[xyz_t[0] - cal_t.ext_par.x0,
+                     xyz_t[1] - cal_t.ext_par.y0, 0])
 
             if R > Rmax:
                 Rmax = R
 
     # round values (-> enlarge)
-    Rmax += rw - math.fmod(Rmax, rw)
+    Rmax += rw - np.fmod(Rmax, rw)
 
     # get # of rasterlines in r, z
     nr = int(Rmax / rw + 1)
@@ -334,11 +337,11 @@ def init_mmlut(vpar: VolumePar, cpar: ControlPar, cal: Calibration) -> Calibrati
 
         for i in range(nr):
             for j in range(nz):
-                xyz = vec_set(Ri[i] + cal_t.ext_par.x0,
-                              cal_t.ext_par.y0, Zi[j])
+                xyz = np.r_[Ri[i] + cal_t.ext_par.x0,
+                              cal_t.ext_par.y0, Zi[j]]
                 data.flat[i * nz + j] = multimed_r_nlay(cal_t, cpar.mm, xyz)
 
-        print("filled mmlut data with {data}")
+        # print(f"filled mmlut data with {data}")
         cal.mmlut.data = data
 
     return cal

openptv_python/orientation.py

@@ -162,15 +162,15 @@ def num_deriv_exterior(
     steps = [dpos, dpos, dpos, dang, dang, dang]
 
     # print(f"exterior = {cal.ext_par}")
-    cal.ext_par.update_rotation_matrix()
+    cal.update_rotation_matrix()
     xs, ys = img_coord(pos, cal, cpar.mm)
     # print(f"  xs = {xs}, ys = {ys}")
 
     for pd in range(6):
-        cal.ext_par.increment_attribute(var[pd], steps[pd])
+        cal.increment_attribute(var[pd], steps[pd])
         # print(f"exterior = {cal.ext_par}")
         if pd > 2:
-            cal.ext_par.update_rotation_matrix()
+            cal.update_rotation_matrix()
 
         xpd, ypd = img_coord(pos, cal, cpar.mm)
         # print(f" xpd = {xpd}, ypd = {ypd}")
@@ -179,10 +179,10 @@ def num_deriv_exterior(
 
         # print(f"   x_ders[{pd}] = {x_ders[pd]}, y_ders[{pd}] = {y_ders[pd]}")
 
-        cal.ext_par.increment_attribute(var[pd], -steps[pd])
+        cal.increment_attribute(var[pd], -steps[pd])
         # print(f"exterior = {cal.ext_par}")
 
-    cal.ext_par.update_rotation_matrix()
+    cal.update_rotation_matrix()
 
     return (x_ders, y_ders)
 
@@ -335,7 +335,7 @@ def orient(
             xc, yc = correct_brown_affine(xc, yc, cal.added_par)
 
             # Projected 2D position on sensor of corresponding known point
-            cal.ext_par.update_rotation_matrix()
+            cal.update_rotation_matrix()
             xp, yp = img_coord(fix[i], cal, cpar.mm)
 
             # derivatives of distortion parameters
@@ -382,7 +382,7 @@ def orient(
 
             # Num. deriv. of projection coords over sensor distance from PP
             cal.int_par.cc += dm
-            cal.ext_par.update_rotation_matrix()
+            cal.update_rotation_matrix()
             xpd, ypd = img_coord(fix[i], cal, cpar.mm)
             X[n][6] = (xpd - xp) / dm
             X[n + 1][6] = ypd - yp
@@ -571,7 +571,7 @@ def orient(
         sigmabeta[i] = sigmabeta[NPAR] * np.sqrt(XPX[i][i])
 
     if stopflag:
-        cal.ext_par.update_rotation_matrix()
+        cal.update_rotation_matrix()
         return resi
     else:
         return None
@@ -611,7 +611,7 @@ def raw_orient(
         for i in range(nfix):
             xc, yc = pixel_to_metric(pix[i].x, pix[i].y, cpar)
             pos = fix[i]
-            cal.ext_par.update_rotation_matrix()
+            cal.update_rotation_matrix()
             xp, yp = img_coord(pos, cal, cpar.mm)
             X[n], X[n + 1] = num_deriv_exterior(cal, cpar, dm, drad, pos)
             y[n] = xc - xp
@@ -651,7 +651,7 @@ def raw_orient(
         cal.ext_par.kappa += beta[5]
 
     if stopflag:
-        cal.ext_par.update_rotation_matrix()
+        cal.update_rotation_matrix()
 
     return stopflag
 

tests/test_corresp.py

@@ -468,40 +468,3 @@ def test_four_camera_matching(self):
 
         # Assert that 16 triplets were matched.
         self.assertEqual(matched, 16)
-
-
-class TestSafelyAllocateAdjacencyLists(unittest.TestCase):
-    """Test the safely_allocate_adjacency_lists function."""
-
-    def test_correct_list_size(self):
-        """Test that the adjacency lists are correctly sized."""
-        num_cams = 5
-        target_counts = [3, 5, 2, 4, 1]
-        lists = safely_allocate_adjacency_lists(num_cams, target_counts)
-        self.assertEqual(len(lists), num_cams)
-        for i in range(num_cams):
-            self.assertEqual(len(lists[i]), num_cams)
-            for j in range(num_cams):
-                if i < j:
-                    self.assertTrue(len(lists[i][j]) >= target_counts[i]) # recarray is one length
-
-    def test_memory_error(self):
-        """Memory stress test."""
-        # available_memory = 8GB = 8 * 1024 * 1024 * 1024 bytes
-        # overhead = 200MB = 200 * 1024 * 1024 bytes
-        # item_size = 4 bytes (for integers)
-
-        # max_items = (8 * 1024 * 1024 * 1024 - 200 * 1024 * 1024) // 4 = 1,995,116,800
-
-        num_cams = 4
-        target_counts = [1000, 1000, 1000, 1000]
-        # with self.assertRaises(MemoryError):
-        _ = safely_allocate_adjacency_lists(num_cams, target_counts)
-
-        # target_counts = [int(1e3), int(1e3), int(1e3), int(1e10)]
-        # with self.assertRaises(MemoryError):
-        #     safely_allocate_adjacency_lists(num_cams, target_counts)
-
-
-if __name__ == "__main__":
-    unittest.main()

tests/test_epipolar.py

@@ -86,7 +86,10 @@ def test_two_cameras(self):
 
     def test_epi_mm_2D(self):
         """Test the epi_mm_2D function."""
-        test_Ex = Exterior(0.0, 0.0, 100.0, 0.0, 0.0, 0.0)
+        test_Ex = Exterior.copy()
+        test_Ex.z0 = 100.0
+
+
         test_I = Interior(0.0, 0.0, 100.0)
         test_G = np.array((0.0, 0.0, 50.0))
         test_addp = ap_52(0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0)
@@ -124,13 +127,20 @@ def test_epi_mm_2D(self):
 
     def test_epi_mm(self):
         """Test the epi_mm function."""
-        test_Ex = Exterior(10.0, 0.0, 100.0, 0.0, -0.01, 0.0)
+        test_Ex = Exterior.copy()
+        test_Ex.x0 = 10.0
+        test_Ex.z0 = 100.0
+        test_Ex.phi = -.01
+
         test_I = Interior(0.0, 0.0, 100.0)
         test_G = np.array((0.0, 0.0, 50.0))
         test_addp = ap_52(0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0)
         test_cal_1 = Calibration(test_Ex, test_I, test_G, test_addp)
 
-        test_Ex_2 = Exterior(-10.0, 0.0, 100.0, 0.0, 0.01, 0.0)
+        test_Ex_2 = Exterior.copy()
+        test_Ex_2.x0 = -10.0
+        test_Ex_2.z0 = 100.0
+        test_Ex_2.phi = .01
         test_cal_2 = Calibration(test_Ex_2, test_I, test_G, test_addp)
 
         test_mm = MultimediaPar(1, 1.0, [1.49, 0.0, 0.0], [
@@ -155,13 +165,16 @@ def test_epi_mm(self):
 
     def test_epi_mm_perpendicular(self):
         """Test the epi_mm function."""
-        test_Ex = Exterior(0.0, 0.0, 100.0, 0.0, 0.0, 0.0)
+        test_Ex = Exterior.copy()
+        test_Ex.z0 = 100.0
         test_I = Interior(0.0, 0.0, 100.0)
         test_G = np.array((0.0, 0.0, 50.0))
         test_addp = ap_52(0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0)
         test_cal_1 = Calibration(test_Ex, test_I, test_G, test_addp)
 
-        test_Ex_2 = Exterior(100.0, 0.0, 0.0, 0.0, 1.57, 0.0)
+        test_Ex_2 = Exterior.copy()
+        test_Ex_2.x0 = 100.0
+        test_Ex_2.phi = 1.57
         test_cal_2 = Calibration(test_Ex_2, test_I, test_G, test_addp)
 
         test_mm = MultimediaPar(1, 1.0, [1.0, 0.0, 0.0], [1.0, 0.0, 0.0], 1.0)
@@ -232,7 +245,9 @@ def test_find_candidate(self):
         ny = 3
         sumg = 100
 
-        test_Ex = Exterior(0.0, 0.0, 100.0, 0.0, 0.0, 0.0)
+        test_Ex = Exterior.copy()
+        test_Ex.z0 = 100.
+
         test_I = Interior(0.0, 0.0, 100.0)
         test_G = np.array((0.0, 0.0, 50.0))
         test_addp = ap_52(0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0)

tests/test_imgcoord.py

@@ -20,7 +20,9 @@ def test_flat_centered_cam(self):
         # gleaned from simple geometry.
         pos = vec_set(10, 5, -20)
         cal = Calibration(
-            ext_par=Exterior(0, 0, 40, 0, 0, 0, [[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+            ext_par = np.array((0, 0, 40, 0, 0, 0,
+                                [[1, 0, 0], [0, 1, 0], [0, 0, 1]]),
+                               dtype = Exterior.dtype).view(np.recarray),
             int_par=Interior(0, 0, 10),
             glass_par=np.array((0., 0., 20.)),
             added_par=ap_52(0, 0, 0, 0, 0, 1, 0),