Fix half-pixel offsets when converting between cube maps and equirect…

…angular

There were multiple sources of half-pixel offsets, which caused small inaccuracies when converting between cube maps and equirectangular images. These are now fixed.

Fixes issue haruishi43#8

equilib/cube2equi/numpy.py

@@ -128,10 +128,13 @@ def _equirect_facetype(h: int, w: int) -> np.ndarray:
         np.arange(4).repeat(w // 4)[None, :].repeat(h, 0), 3 * w // 8, 1
     )
 
+    half_pixel_angular_width = np.pi / w
+    pixel_angular_height = np.pi / h
+
     # Prepare ceil mask
     mask = np.zeros((h, w // 4), np.bool)
-    idx = np.linspace(-np.pi, np.pi, w // 4) / 4
-    idx = h // 2 - np.around(np.arctan(np.cos(idx)) * h / np.pi)
+    idx = np.linspace(-np.pi + half_pixel_angular_width, np.pi - half_pixel_angular_width, num=w // 4) / 4
+    idx = h // 2 - np.around(np.arctan(np.cos(idx)) * h / (np.pi - pixel_angular_height))
     idx = idx.astype(int_dtype)
     for i, j in enumerate(idx):
         mask[:j, i] = 1
@@ -150,8 +153,10 @@ def create_equi_grid(
     batch: int,
     dtype: np.dtype = np.dtype(np.float32),
 ) -> np.ndarray:
-    theta = np.linspace(-np.pi, np.pi, num=w_out, dtype=dtype)
-    phi = np.linspace(np.pi, -np.pi, num=h_out, dtype=dtype) / 2
+    half_pixel_angular_width = np.pi / w_out
+    half_pixel_angular_height = np.pi / h_out / 2
+    theta = np.linspace(-np.pi + half_pixel_angular_width, np.pi - half_pixel_angular_width, num=w_out, dtype=dtype)
+    phi = np.linspace(np.pi / 2 - half_pixel_angular_height, -np.pi / 2 + half_pixel_angular_height, num=h_out, dtype=dtype)
     theta, phi = np.meshgrid(theta, phi)
 
     # Get face id to each pixel: 0F 1R 2B 3L 4U 5D
@@ -179,8 +184,8 @@ def create_equi_grid(
             coor_y[mask] = -c * np.cos(theta[mask])
 
     # Final renormalize
-    coor_x = np.clip(np.clip(coor_x + 0.5, 0, 1) * w_face, 0, w_face - 1)
-    coor_y = np.clip(np.clip(coor_y + 0.5, 0, 1) * w_face, 0, w_face - 1)
+    coor_x = np.clip(coor_x + 0.5, 0, 1) * w_face
+    coor_y = np.clip(coor_y + 0.5, 0, 1) * w_face
 
     # change x axis of the x coordinate map
     for i in range(6):
@@ -189,6 +194,7 @@ def create_equi_grid(
 
     grid = np.stack((coor_y, coor_x), axis=0)
     grid = np.concatenate([grid[np.newaxis, ...]] * batch)
+    grid = grid - 0.5 # Offset pixel center
     return grid
 
 
@@ -254,7 +260,7 @@ def run(
             img=horizon, grid=grid, out=out, mode=mode
         )
     else:
-        out = numpy_grid_sample(img=horizon, grid=grid, out=out, mode=mode)
+        out = numpy_grid_sample(img=horizon, grid=grid, out=out, mode=mode, clamp_x=True, clamp_y=True)
 
     out = (
         out.astype(horizon_dtype)

equilib/cube2equi/torch.py

@@ -172,14 +172,14 @@ def create_equi_grid(
     dtype: torch.dtype = torch.float32,
     device: torch.device = torch.device("cpu"),
 ) -> torch.Tensor:
+
+    half_pixel_angular_width = math.pi / w_out
+    half_pixel_angular_height = math.pi / h_out / 2
     theta = torch.linspace(
-        -math.pi, math.pi, steps=w_out, dtype=dtype, device=device
+        -math.pi + half_pixel_angular_width, math.pi - half_pixel_angular_width, steps=w_out, dtype=dtype, device=device
     )
-    phi = (
-        torch.linspace(
-            math.pi, -math.pi, steps=h_out, dtype=dtype, device=device
-        )
-        / 2
+    phi = torch.linspace(
+        math.pi / 2 - half_pixel_angular_height, -math.pi / 2 + half_pixel_angular_height, steps=h_out, dtype=dtype, device=device
     )
     phi, theta = torch.meshgrid([phi, theta])
 
@@ -213,12 +213,8 @@ def create_equi_grid(
             coor_y[mask] = -c * torch.cos(theta[mask])
 
     # Final renormalize
-    coor_x = torch.clamp(
-        torch.clamp(coor_x + 0.5, 0, 1) * w_face, 0, w_face - 1
-    )
-    coor_y = torch.clamp(
-        torch.clamp(coor_y + 0.5, 0, 1) * w_face, 0, w_face - 1
-    )
+    coor_x = torch.clamp(coor_x + 0.5, 0, 1) * w_face
+    coor_y = torch.clamp(coor_y + 0.5, 0, 1) * w_face
 
     # change x axis of the x coordinate map
     for i in range(6):
@@ -230,6 +226,7 @@ def create_equi_grid(
     coor_y = coor_y.repeat(batch, 1, 1)
 
     grid = torch.stack((coor_y, coor_x), dim=-3).to(device)
+    grid = grid - 0.5 # offset pixel center
     return grid
 
 

equilib/equi2cube/numpy.py

@@ -100,8 +100,6 @@ def convert_grid(
         # but it was faster separately
         ui = (theta - np.pi) * w_equi / (2 * np.pi)
         uj = (np.pi / 2 - phi) * h_equi / np.pi  # NOTE: fixed here
-        ui += 0.5
-        uj += 0.5
         ui %= w_equi
         uj %= h_equi
     elif method == "faster":
@@ -112,8 +110,6 @@ def convert_grid(
         # out of range, it will not work with `faster`
         ui = (theta - np.pi) * w_equi / (2 * np.pi)
         uj = (np.pi / 2 - phi) * h_equi / np.pi  # NOTE: fixed here
-        ui += 0.5
-        uj += 0.5
         ui = np.where(ui < 0, ui + w_equi, ui)
         ui = np.where(ui >= w_equi, ui - w_equi, ui)
         uj = np.where(uj < 0, uj + h_equi, uj)
@@ -123,7 +119,7 @@ def convert_grid(
 
     # stack the pixel maps into a grid
     grid = np.stack((uj, ui), axis=-3)
-
+    grid = grid - 0.5 # offset pixel center
     return grid
 
 

equilib/equi2cube/torch.py

@@ -94,15 +94,11 @@ def convert_grid(
     if method == "robust":
         ui = (theta - pi) * w_equi / (2 * pi)
         uj = (pi / 2 - phi) * h_equi / pi  # FIXME: fixed here
-        ui += 0.5
-        uj += 0.5
         ui %= w_equi
         uj %= h_equi
     elif method == "faster":
         ui = (theta - pi) * w_equi / (2 * pi)
         uj = (pi / 2 - phi) * h_equi / pi  # FIXME: fixed here
-        ui += 0.5
-        uj += 0.5
         ui = torch.where(ui < 0, ui + w_equi, ui)
         ui = torch.where(ui >= w_equi, ui - w_equi, ui)
         uj = torch.where(uj < 0, uj + h_equi, uj)
@@ -112,7 +108,7 @@ def convert_grid(
 
     # stack the pixel maps into a grid
     grid = torch.stack((uj, ui), dim=-3)
-
+    grid = grid - 0.5 # offset pixel center
     return grid
 
 

equilib/numpy_utils/grid.py

@@ -104,7 +104,8 @@ def create_xyz_grid(
 ) -> np.ndarray:
     """xyz coordinates of the faces of the cube"""
     out = np.zeros((w_face, w_face * 6, 3), dtype=dtype)
-    rng = np.linspace(-0.5, 0.5, num=w_face, dtype=dtype)
+    pixel_half_width = 0.5 / w_face
+    rng = np.linspace(-0.5 + pixel_half_width, 0.5 - pixel_half_width, num=w_face, dtype=dtype)
 
     # Front face (x = 0.5)
     out[:, 0 * w_face : 1 * w_face, [1, 2]] = np.stack(