Added support for omnidirectional cameras

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]
 

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

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 day="24" hour="16" minute="36" month="9" second="49.1429" utc="-4" year="2014"/>
+        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'))

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)