From 1837e8f0fd5f70ca3ba5e0919631480fb8734134 Mon Sep 17 00:00:00 2001 From: Ian Maquignaz Date: Sat, 27 May 2023 11:28:34 -0400 Subject: [PATCH] Added support for omnidirectional cameras --- envmap/environmentmap.py | 77 ++++++++++++++++++++++++++- envmap/projections.py | 112 +++++++++++++++++++++++++++++++++++++++ envmap/xmlhelper.py | 81 +++++++++++++++++++++++++++- test/__init__.py | 0 test/test_projections.py | 29 +++++++++- 5 files changed, 294 insertions(+), 5 deletions(-) create mode 100644 test/__init__.py diff --git a/envmap/environmentmap.py b/envmap/environmentmap.py index 862ab47..9597d1a 100644 --- a/envmap/environmentmap.py +++ b/envmap/environmentmap.py @@ -144,6 +144,79 @@ def fromSkybox(cls, top, bottom, left, right, front, back): return cls(cube, format_="cube") + + @classmethod + def from_omnicam(cls, im, targetDim, targetFormat='skyangular', OcamCalib_=None, copy=True, order=1): + """ + Creates an EnvironmentMap from Omnidirectional Camera (OcamCalib) capture. + + :param im: Image path (str, pathlib.Path) or data (np.ndarray) representing + an ocam image. + :param OcamCalib: OCamCalib calibration dictionary. If not provided and + param im is an image path, then calibration will be loaded directly + from matching ".meta.xml" file. + :param targetFormat: Target format. + :param targetDim: Target dimension. + :param order: Interpolation order (0: nearest, 1: linear, ..., 5). + :param copy: When a numpy array is given, should it be copied. + + :type im: str, Path, np.ndarray + :type OcamCalib: dict + :type targetFormat: string + :type targetDim: integer + :type order: integer (0,1,...,5) + :type copy: bool + """ + + if not OcamCalib_ and isPath(im): + filename = os.path.splitext(str(im))[0] + metadata = EnvmapXMLParser("{}.meta.xml".format(filename)) + OcamCalib_ = metadata.get_calibration() + + assert OcamCalib_ is not None, ( + "Please provide OCam (metadata file not found).") + + if isPath(im): + # We received the filename + data = imread(str(im)) + elif type(im).__module__ == np.__name__: + # We received a numpy array + data = np.asarray(im, dtype='double') + if copy: + data = data.copy() + else: + raise Exception('Could not understand input. Please provide a ' + 'filename (str) or an image (np.ndarray).') + + e = EnvironmentMap(targetDim, targetFormat) + dx, dy, dz, valid = e.worldCoordinates() + u,v = world2ocam(dx, dy, dz, OcamCalib_) + + # Interpolate + # Repeat the first and last rows/columns for interpolation purposes + h, w, d = data.shape + source = np.empty((h + 2, w + 2, d)) + + source[1:-1, 1:-1] = data + source[0,1:-1] = data[0,:]; source[0,0] = data[0,0]; source[0,-1] = data[0,-1] + source[-1,1:-1] = data[-1,:]; source[-1,0] = data[-1,0]; source[-1,-1] = data[-1,-1] + source[1:-1,0] = data[:,0] + source[1:-1,-1] = data[:,-1] + + # To avoid displacement due to the padding + u += 0.5/data.shape[1] + v += 0.5/data.shape[0] + target = np.vstack((u.flatten()*data.shape[0], v.flatten()*data.shape[1])) + + data = np.zeros((u.shape[0], u.shape[1], d)) + for c in range(d): + map_coordinates(source[:,:,c], target, output=data[:,:,c].reshape(-1), cval=np.nan, order=order, prefilter=filter) + e.data = data + + # Done + return e + + def __hash__(self): """Provide a hash of the environment map type and size. Warning: doesn't take into account the data, just the type, @@ -228,7 +301,7 @@ def world2image(self, x, y, z): def world2pixel(self, x, y, z): """Returns the (u, v) coordinates (in the interval defined by the MxN image).""" - # Get (u,v) in [-1, 1] interval + # Get (u,v) in [0, 1] interval u,v = self.world2image(x, y, z) # de-Normalize coordinates to interval defined by the MxN image @@ -241,7 +314,7 @@ def world2pixel(self, x, y, z): def pixel2world(self, u, v): """Returns the (x, y, z) coordinates for pixel cordinates (u,v)(in the interval defined by the MxN image).""" - # Normalize coordinates to [-1, 1] interval + # Normalize coordinates to [0, 1] interval u = (u+0.5) / self.data.shape[1] v = (v+0.5) / self.data.shape[0] diff --git a/envmap/projections.py b/envmap/projections.py index 94f1336..fbcf07a 100644 --- a/envmap/projections.py +++ b/envmap/projections.py @@ -1,6 +1,7 @@ import numpy as np from numpy import logical_and as land, logical_or as lor +from .rotations import * def world2latlong(x, y, z): """Get the (u, v) coordinates of the point defined by (x, y, z) for @@ -295,3 +296,114 @@ def cube2world(u, v): return x.item(), y.item(), z.item(), valid.item() return x, y, z, valid + + +def ocam2world(u, v, ocam_calibration): + """ Project a point (u, v) in omnidirectional camera space to + (x, y, z) point in world coordinate space.""" + + # Step 0. De-Normalize coordinates to interval defined by the MxN image + # Where u=cols(x), v=rows(y) + width_cols = ocam_calibration['width'] + height_rows = ocam_calibration['height'] + + v = np.floor(v*height_rows).astype(int) + u = np.floor(u*width_cols).astype(int) + + # Step 1. Center & Skew correction + # M = affine matrix [c,d,xc; e,1,yc; 0,0,1] + # p' = M^-1 * (p) + M = ocam_calibration['affine_3x3'] + F = ocam_calibration['F'] + + # Inverse: M_ = M^-1 + M_ = np.linalg.inv(M) + + # Affine transform + w = np.ones_like(u) + assert u.shape == v.shape + save_original_shape = u.shape + p_uvw = np.array((v.reshape(-1),u.reshape(-1),w.reshape(-1))) + p_xyz = np.matmul(M_,p_uvw) + + # Add epsilon to mitigate NAN + p_xyz[ p_xyz==0 ] = np.finfo(p_xyz.dtype).eps + + # Step 2. Get unit-sphere world coordinate z + # Distance to center of image: p = sqrt(X^2 + Y^2) + p_z = np.linalg.norm(p_xyz[0:2], axis=0) + # Convert to z-coordinate with p_z = F(p) + p_xyz[2] = F(p_z) + + # Step 3. Normalize x,y,z to unit length of 1 (unit sphere) + p_xyz = p_xyz / np.linalg.norm(p_xyz, axis=0) + + # Step 4. Fix coordinate system alignment + # (rotate -90 degrees around z-axis) + # (x,y,z) -> (x,-z,y) as +y is up (not +z) + p_xyz = np.matmul(rotz(np.deg2rad(-90)),p_xyz) + x,y,z = ( + p_xyz[0].reshape(save_original_shape), + -p_xyz[2].reshape(save_original_shape), + -p_xyz[1].reshape(save_original_shape) + ) + + valid = np.ones(x.shape, dtype='bool') + return x,y,z, valid + + +def world2ocam(x, y, z, ocam_calibration): + """ Project a point (x, y, z) in world coordinate space to + a (u, v) point in omnidirectional camera space.""" + + # Step 1. Center & Skew correction + # M = affine matrix [c,d; e,1] + # T = translation vector + # p' = M^-1 * (p - T) + M = ocam_calibration['affine_3x3'] + F = ocam_calibration['F'] + + assert x.shape == y.shape and x.shape == z.shape, f'{x.shape} == {y.shape} == {z.shape}' + save_original_shape = x.shape + + # Step 2. Fix coordinate system alignment + # (x,y,z) -> (x,z,-y) as +z is up (not +y) + p_xyz = np.array((x.reshape(-1),y.reshape(-1),-z.reshape(-1))) + + # Add epsilon to mitigate NAN + p_xyz[ p_xyz==0 ] = np.finfo(p_xyz.dtype).eps + + # (rotate 90 degrees around z-axis) + p_xyz = np.array((p_xyz[0], p_xyz[2], -p_xyz[1])) + p_xyz = np.matmul(rotz(np.deg2rad(90)),p_xyz) + + # Step 3. 3D to 2D + m = p_xyz[2] / np.linalg.norm(p_xyz[0:2], axis=0) + + def poly_inverse(y): + F_ = F.copy() + F_.coef[1] -=y + F_r = F_.roots() + F_r = F_r[ (F_r >= 0) & (F_r.imag == 0) ] + if len(F_r)>0: + return F_r[0].real + else: + return np.nan + m_ = np.vectorize(poly_inverse)(m) + + uvw = p_xyz / np.linalg.norm(p_xyz[0:2], axis=0) * m_ + + # Step 4. Affine transform for center and skew + uvw[2,:] = 1 + uvw = np.nan_to_num(uvw) + uvw = np.matmul(M, uvw) + u,v = uvw[1].reshape(save_original_shape), uvw[0].reshape(save_original_shape) + + # Step 3. Normalize coordinates to interval [0,1] + # Where u=cols(x), v=rows(y) + width_cols = ocam_calibration['width'] + height_rows = ocam_calibration['height'] + u = (u+0.5) / width_cols + v = (v+0.5) / height_rows + + return u,v diff --git a/envmap/xmlhelper.py b/envmap/xmlhelper.py index cf59751..867deb3 100644 --- a/envmap/xmlhelper.py +++ b/envmap/xmlhelper.py @@ -1,3 +1,5 @@ +import numpy as np +import datetime as dt import xml.etree.ElementTree as ET @@ -9,26 +11,101 @@ def __init__(self, filename): self.tree = ET.parse(filename) self.root = self.tree.getroot() + def _getFirstChildTag(self, tag): for elem in self.root: if elem.tag == tag: return elem.attrib + def _getAttrib(self, node, attribute, default=None): if node: return node.get(attribute, default) return default + def getFormat(self): """Returns the format of the environment map.""" node = self._getFirstChildTag('data') return self._getAttrib(node, 'format', 'Unknown') + def getDate(self): - """Returns the date of the environment mapin dict format.""" + """Returns the date of the environment map in dict format.""" return self._getFirstChildTag('date') + + def get_datetime(self): + """Returns the date of the environment map in datetime format.""" + # Example: + date = self.root.find('date') + year = date.get('year') + month = date.get('month').zfill(2) + day = date.get('day').zfill(2) + hour = date.get('hour').zfill(2) + minute = date.get('minute').zfill(2) + second = str(int(float(date.get('second')))).zfill(2) + utc_offset = int(date.get('utc')) + if utc_offset > 0: + utc_offset = f'+{str(utc_offset).zfill(2)}' + else: + utc_offset = f'-{str(np.abs(utc_offset)).zfill(2)}' + return dt.datetime.fromisoformat(f"{year}-{month}-{day} {hour}:{minute}:{second}{utc_offset}:00") + + + def get_calibration(self): + """Returns the OCamCalib calibration metadata.""" + + calibration = {} + # Calibration Model + node = self.root.find('calibrationModel') + # Affine 2D = [c,d; + # e,1]; + c = float(node.get('c')) + d = float(node.get('d')) + e = float(node.get('e')) + affine_2x2 = np.array([[c,d],[e, 1]]) + calibration['c'] = c + calibration['d'] = d + calibration['e'] = e + calibration['affine_2x2'] = affine_2x2 + + # shape = [height, width] + height = int(node.get('height')) + width = int(node.get('width')) + calibration['height'] = height + calibration['width'] = width + calibration['shape'] = (height, width) + + # center = [xc, yc] + xc = float(node.get('xc')) + yc = float(node.get('yc')) + calibration['xc'] = xc + calibration['yc'] = yc + calibration['center'] = (xc, yc) + + # Affine 3D = [c,d,xc; + # e,1,yc; + # 0,0,1]; + affine_3x3 = np.array([[c,d,xc],[e, 1, yc],[0,0,1]]) + calibration['affine_3x3'] = affine_3x3 + + # ss = [a_0, a_1, ..., a_n] + ss = [ float(s.get('s')) for s in node.findall('ss') ] + calibration['ss'] = np.array(ss) + polydomain = xc if xc > yc else yc + polydomain = [-polydomain,polydomain] + calibration['F'] = \ + np.polynomial.polynomial.Polynomial( + ss, + domain=polydomain, + window=polydomain + ) + + return calibration + + def getExposure(self): """Returns the exposure of the environment map in EV.""" node = self._getFirstChildTag('exposure') - return self._getAttrib(node, 'EV') + return float(self._getAttrib(node, 'EV')) diff --git a/test/__init__.py b/test/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/test/test_projections.py b/test/test_projections.py index 92e0b7d..e8103a8 100644 --- a/test/test_projections.py +++ b/test/test_projections.py @@ -1,5 +1,4 @@ import pytest -import math import numpy as np from envmap import projections as t @@ -113,3 +112,31 @@ def test_projections_image(format_): np.testing.assert_array_almost_equal(u[valid], u_[valid], decimal=5) np.testing.assert_array_almost_equal(v[valid], v_[valid], decimal=5) + + +@pytest.mark.parametrize("calibration", + [ + # ( ocam calibration dict, from 20150909_144054_stack.meta.xml ) + ({ + 'height': 3840, + 'width': 5760, + 'F' : np.polynomial.polynomial.Polynomial( + np.array([-1283.8735,0,0.00035359,-1.2974e-07,6.7764e-11]), + domain=[-2895.3481,2895.3481], + window=[-2895.3481,2895.3481] + ), + 'affine_3x3':np.array([[0.99981,0.00024371,1904.2826],[3.7633e-06, 1, 2895.3481],[0,0,1]]), + }), + ] +) +def test_ocam(calibration): + + e = env.EnvironmentMap(64, 'skyangular', channels=2) + x,y,z, valid = e.worldCoordinates() + + u_, v_ = t.world2ocam(x, y, z, calibration) + x_, y_, z_, V = t.ocam2world(u_, v_, calibration) + + np.testing.assert_array_almost_equal(x[valid], x_[valid], decimal=3) + np.testing.assert_array_almost_equal(y[valid], y_[valid], decimal=3) + np.testing.assert_array_almost_equal(z[valid], z_[valid], decimal=3)