From 229b0d15233cbc305a141f054cbe5d1fd4a10340 Mon Sep 17 00:00:00 2001
From: "lodenthal@uni-bielefeld.de" <lodenthal@uni-bielefeld.de>
Date: Fri, 7 Sep 2018 17:11:38 +0200
Subject: [PATCH] repaired test after changes

---
 navipy/maths/euler.py                       |  98 ++++++++++--
 navipy/maths/homogeneous_transformations.py |  15 --
 navipy/maths/quaternion.py                  | 165 +++++++++++++-------
 navipy/processing/mcode.py                  |   4 +-
 navipy/processing/test_OpticFlow.py         |  48 +++---
 5 files changed, 224 insertions(+), 106 deletions(-)

diff --git a/navipy/maths/euler.py b/navipy/maths/euler.py
index 9c3e651..3b6b242 100644
--- a/navipy/maths/euler.py
+++ b/navipy/maths/euler.py
@@ -17,6 +17,16 @@ _AXES2TUPLE = {
 
 
 def R1(a):
+    """ rotation matrix around the x- axis
+
+    :param a: angle in degrees to be rotated
+    :returns: a matrix
+    :rtype: (np.ndarray)
+    ..ref: James Diebel
+           "Representing Attitude: Euler Angles, Unit Quaternions, and Rotation
+           Vectors."
+           (2006): p. 5 (2.4).
+    """
     if not isinstance(a, float) and not isinstance(a, int):
         raise TypeError("angle must be numeric value")
     R1 = np.array([[1, 0, 0],
@@ -26,6 +36,16 @@ def R1(a):
 
 
 def R2(a):
+    """ rotation matrix around the y- axis
+
+    :param a: angle in degrees to be rotated
+    :returns: a matrix
+    :rtype: (np.ndarray)
+    ..ref: James Diebel
+           "Representing Attitude: Euler Angles, Unit Quaternions, and Rotation
+           Vectors."
+           (2006): p. 5 (2.4).
+    """
     if not isinstance(a, float) and not isinstance(a, int):
         raise TypeError("angle must be numeric value")
     R2 = np.array([[c(a), 0, -s(a)],
@@ -35,6 +55,16 @@ def R2(a):
 
 
 def R3(a):
+    """ rotation matrix around the z- axis
+
+    :param a: angle in degrees to be rotated
+    :returns: a matrix
+    :rtype: (np.ndarray)
+    ..ref: James Diebel
+           "Representing Attitude: Euler Angles, Unit Quaternions, and Rotation
+           Vectors."
+           (2006): p. 5 (2.4).
+    """
     if not isinstance(a, float) and not isinstance(a, int):
         raise TypeError("angle must be numeric value")
     R3 = np.array([[c(a), s(a), 0],
@@ -43,7 +73,20 @@ def R3(a):
     return R3
 
 
-def matrix(ai, aj, ak, axes='rxyz'):
+def matrix(ai, aj, ak, axes='sxyz'):
+    """ rotation matrix around the three axis with the
+        order given by the axes parameter
+
+    :param ai: angle in degrees to be rotated about the first axis
+    :param aj: angle in degrees to be rotated about the second axis
+    :param ak: angle in degrees to be rotated about the third axis
+    :returns: a matrix
+    :rtype: (np.ndarray)
+    ..ref: James Diebel
+           "Representing Attitude: Euler Angles, Unit Quaternions, and Rotation
+           Vectors."
+           (2006): p. 9
+    """
     if axes not in list(_AXES2TUPLE.keys()):
         raise ValueError("the chosen convention is not supported")
     r, i, j, k = _AXES2TUPLE[axes]
@@ -52,7 +95,7 @@ def matrix(ai, aj, ak, axes='rxyz'):
                   np.dot(matrixes[j](aj),
                   matrixes[k](ak)))
     ID = np.identity(4)
-    ID[1:4,1:4] = Rijk
+    ID[:3, :3] = Rijk
     Rijk = ID
     if not r:
         return Rijk
@@ -66,6 +109,10 @@ def from_matrix(matrix, axes='sxyz'):
     axes : One of 24 axis sequences as string or encoded tuple
 
     Note that many Euler angle triplets can describe one matrix.
+    ..ref: James Diebel
+           "Representing Attitude: Euler Angles, Unit Quaternions, and Rotation
+           Vectors."
+           (2006): p. 23 - 31.
     """
     if not isinstance(matrix, np.ndarray) and not isinstance(matrix, list):
         raise TypeError("matrix must be np.array or list")
@@ -76,7 +123,7 @@ def from_matrix(matrix, axes='sxyz'):
     if not is_numeric_array(matrix):
         raise ValueError("matrix must contain numeric values")
     if matrix.shape[0] > 3:
-        matrix = matrix[1:4,1:4]
+        matrix = matrix[:3, :3]
     u = None
     rot, i, j, k = _AXES2TUPLE[axes]
     if rot:
@@ -132,10 +179,6 @@ def from_matrix(matrix, axes='sxyz'):
         u = [np.arctan2(matrix[1, 2], -matrix[0, 2]),
              np.arccos(matrix[2, 2]),
              np.arctan2(matrix[2, 1], matrix[2, 0])]
-
-    if u is None:
-        print("conv", axes, matrix)
-
     return u
 
 
@@ -149,16 +192,30 @@ def from_quaternion(quaternion, axes='sxyz'):
     #    raise ValueError('posorient must not contain nan')
     if axes not in list(_AXES2TUPLE.keys()):
         raise ValueError("the chosen convention is not supported")
-    print(quat.matrix(quaternion))
-    return from_matrix(quat.matrix(quaternion)[1:4,1:4], axes)
+    return from_matrix(quat.matrix(quaternion)[:3,:3], axes)
 
 
-def angle_rate_matrix(ai, aj, ak, axes='rxyz'):
+def angle_rate_matrix(ai, aj, ak, axes='sxyz'):
     """
     Return the Euler Angle Rates Matrix
 
     from Diebels Representing Attitude: Euler Angles,
     Unit Quaternions, and Rotation, 2006
+    rotation matrix around the three axis with the
+        order given by the axes parameter
+
+    :param ai: angle in degrees to be rotated about the first axis
+    :param aj: angle in degrees to be rotated about the second axis
+    :param ak: angle in degrees to be rotated about the third axis
+    :param axes: string representation for the order of axes to be rotated
+                 around and whether stationary or rotational
+                 axes should be used
+    :returns: a matrix
+    :rtype: (np.ndarray)
+    ..ref: James Diebel
+           "Representing Attitude: Euler Angles, Unit Quaternions, and Rotation
+           Vectors."
+           (2006): p. 9 (5.2)
     """
     if not isinstance(ai, float) and not isinstance(ai, int):
         raise TypeError("euler angle must be of type float")
@@ -197,9 +254,28 @@ def angle_rate_matrix(ai, aj, ak, axes='rxyz'):
     return rotM
 
 
-def angular_velocity(ai, aj, ak, dai, daj, dak, axes='rxyz'):
+def angular_velocity(ai, aj, ak, dai, daj, dak, axes='sxyz'):
     """
     Return the angular velocity
+
+    :param ai: angle in degrees to be rotated about the first axis
+    :param aj: angle in degrees to be rotated about the second axis
+    :param ak: angle in degrees to be rotated about the third axis
+    :param dai: time derivative in degrees/time of the angle to be rotated
+                about the first axis
+    :param daj: time derivative in degrees/time of the angle to be rotated
+                about the second axis
+    :param dak: time derivative in degrees/time of the angle to be rotated
+                about the third axis
+    :param axes: string representation for the order of axes to be rotated
+                 around and whether stationary or rotational axes should
+                 be used
+    :returns: a matrix
+    :rtype: (np.ndarray)
+    ..ref: James Diebel
+           "Representing Attitude: Euler Angles, Unit Quaternions, and Rotation
+           Vectors."
+           (2006): p. 9 (5.2)
     """
     if not isinstance(ai, float) and not isinstance(ai, int):
         raise TypeError("euler angle must be of type float")
diff --git a/navipy/maths/homogeneous_transformations.py b/navipy/maths/homogeneous_transformations.py
index 28e6bee..2ea6ee2 100644
--- a/navipy/maths/homogeneous_transformations.py
+++ b/navipy/maths/homogeneous_transformations.py
@@ -7,19 +7,6 @@ import navipy.maths.quaternion as quaternion
 from navipy.maths.tools import vector_norm
 from navipy.maths.tools import unit_vector
 
-def old2new(matrix):
-    new = np.identity(4)
-    new[1:4, 1:4] = matrix[:3, :3]
-    new[0, :] = np.flip(matrix[3, :], 0)
-    new[:, 0] = np.flip(matrix[:, 3], 0)
-    return new
-
-def new2old(matrix):
-    old = np.identity(4)
-    old[:3, :3] = matrix[1:4, 1:4]
-    old[3, :] = np.flip(matrix[0, :], 0)
-    old[:, 3] = np.flip(matrix[:, 0], 0)
-    return old
 
 def translation_matrix(direction):
     """
@@ -415,7 +402,6 @@ def compose_matrix(scale=None, shear=None, angles=None, translate=None,
             R = quaternion.matrix(angles)
         else:
             R = euler.matrix(angles[0], angles[1], angles[2], axes)
-            R = new2old(R)
         M = np.dot(M, R)
     if shear is not None:
         Z = np.identity(4)
@@ -441,4 +427,3 @@ def is_same_transform(matrix0, matrix1):
     matrix1 = np.array(matrix1, dtype=np.float64, copy=True)
     matrix1 /= matrix1[3, 3]
     return np.allclose(matrix0, matrix1)
-
diff --git a/navipy/maths/quaternion.py b/navipy/maths/quaternion.py
index 80ed1bb..8082a41 100644
--- a/navipy/maths/quaternion.py
+++ b/navipy/maths/quaternion.py
@@ -1,15 +1,29 @@
 """
 
+
 """
 import numpy as np
 import navipy.maths.constants as constants
 from navipy.maths.tools import vector_norm
 # import navipy.maths.new_euler as new_euler
 
-def qat(a,n):
-    tmp = np.zeros((1,4))
-    tmp[0,0]= np.cos((1/2)*a)
-    tmp[0,1:4] = np.dot(n,np.sin((1/2)*a))
+
+def qat(a, n):
+    """ axis-angle quaternion function
+    Returns a unit quaternion
+
+    :param a: angle in degrees
+    :param n: unit canonical vector for a specific axis
+    :returns: a vector
+    :rtype: (np.ndarray)
+    ..ref: James Diebel
+           "Representing Attitude: Euler Angles, Unit Quaternions, and Rotation
+           Vectors."
+           (2006): p. 17 (6.12)
+    """
+    tmp = np.zeros((4))
+    tmp[0] = np.cos((1/2)*a)
+    tmp[1:4] = np.dot(n, np.sin((1/2)*a))
     return tmp
 
 
@@ -17,17 +31,28 @@ def from_euler(ai, aj, ak, axes='sxyz'):
     """Return quaternion from Euler angles and axis sequence.
     ai, aj, ak : Euler's roll, pitch and yaw angles
     axes : One of 24 axis sequences as string or encoded tuple
+    :param ai: angle in degrees to be rotated about the first axis
+    :param aj: angle in degrees to be rotated about the second axis
+    :param ak: angle in degrees to be rotated about the third axis
+    :param axes: string that encodes the order of the axes and
+                 whether rotational or stationary axes should be used
+    :returns: a vector
+    :rtype: (np.ndarray)
+    ..ref: James Diebel
+           "Representing Attitude: Euler Angles, Unit Quaternions, and Rotation
+           Vectors."
+           (2006): p. 18 (6.14)
     """
     # shouldnt there be a difference between rotating
     # and stationary frame?
     vects = [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
     r, i, j, k = constants._AXES2TUPLE_new[axes]
-    q1 = qat(ai,vects[i])
-    q2 = qat(aj,vects[j])
-    q3 = qat(ak,vects[k])
-    print("first", np.dot(np.transpose(q1),q2))
-    print("sec", np.dot(np.dot(np.transpose(q1),q2),np.transpose(q3)))
-    qijk = np.dot(np.dot(np.transpose(q1), q2), np.transpose(q3))
+    q1 = qat(ai, vects[i])
+    q2 = qat(aj, vects[j])
+    q3 = qat(ak, vects[k])
+    quaternion = multiply(np.array(q1), np.array(q2))
+    qijk = multiply(quaternion, np.array(q3))
+    # qijk = np.dot(np.dot(np.transpose(q1), q2), np.transpose(q3))
     return qijk
 
     """
@@ -86,40 +111,53 @@ def about_axis(angle, axis):
     q = np.array([0.0, axis[0], axis[1], axis[2]])
     qlen = vector_norm(q)
     if qlen > constants._EPS:
-        q *= np.sin(angle / 2.0) / qlen
-    q[0] = np.cos(angle / 2.0)
+        q *= np.sin(angle/2.0)/qlen
+    q[0] = np.cos(angle/2.0)
     return q
 
 
-def matrix(quaternion, axes = 'rxyz'):
+def matrix(quaternion, axes='sxyz'):
     """Return homogeneous rotation matrix from quaternion.
+    :param quaternion : vector with at least 3 entrences (unit quaternion)
+    :param axes: string that encodes the order of the axes and
+                 whether rotational or stationary axes should be used
+    :returns: a matrix
+    :rtype: (np.ndarray)
+    ..ref: James Diebel
+           "Representing Attitude: Euler Angles, Unit Quaternions, and Rotation
+           Vectors."
+           (2006): p. 15 (6.4)
     """
     q = np.array(quaternion, dtype=np.float64, copy=True)
-    #n = np.dot(q, q)
-    #if n < constants._EPS:
+    # n = np.dot(q, q)
+    # if n < constants._EPS:
     #    return np.identity(4)
-    #q *= np.sqrt(2.0 / n)
-    #print("within q2", q)
-    #q = np.outer(q, q)
-    #print("within q3", q)
+    # q *= np.sqrt(2.0 / n)
+    # print("within q2", q)
+    # q = np.outer(q, q)
+    # print("within q3", q)
     # return np.array([
-    #    [1.0 - q[2, 2] - q[3, 3],	 q[1, 2] - q[3, 0],	 q[1, 3] + q[2, 0], 0.0],
-    #    [	q[1, 2] + q[3, 0], 1.0 - q[1, 1] - q[3, 3],	 q[2, 3] - q[1, 0], 0.0],
-    #    [	q[1, 3] - q[2, 0],	 q[2, 3] + q[1, 0], 1.0 - q[1, 1] - q[2, 2], 0.0],
+    #    [1.0 - q[2, 2] - q[3, 3],	 q[1, 2] - q[3, 0],
+    #     q[1, 3] + q[2, 0], 0.0],
+    #    [	q[1, 2] + q[3, 0], 1.0 - q[1, 1] - q[3, 3],
+    #     q[2, 3] - q[1, 0], 0.0],
+    #    [	q[1, 3] - q[2, 0],	 q[2, 3] + q[1, 0],
+    #     1.0 - q[1, 1] - q[2, 2], 0.0],
     #    [0.0,			 0.0,			 0.0, 1.0]])
-    r , _, _, _ = constants._AXES2TUPLE_new[axes] 
-    #q=q[0]
-    print("within q4", q)
+    r, _, _, _ = constants._AXES2TUPLE_new[axes]
+    # q=q[0]
     if np.sqrt(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]+q[3]*q[3]) != 1:
         q = q/np.sqrt(q[0]*q[0]+q[1]*q[1]+q[2]*q[2]+q[3]*q[3])
-    mat = np.array([[1, 0, 0, 0],
-          [0, q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3],
-                     2*q[1]* q[2] + 2*q[0]*q[3], 2*q[1]*q[3] - 2*q[0] *q[2]],
-          [0, 2*q[1]* q[2] - 2*q[0]*q[3],
+    mat = np.array([[q[0]*q[0] + q[1]*q[1] - q[2]*q[2] - q[3]*q[3],
+                     2*q[1] * q[2] + 2*q[0]*q[3],
+                     2*q[1]*q[3] - 2*q[0] * q[2], 0],
+                    [2 * q[1] * q[2] - 2 * q[0] * q[3],
                      q[0]*q[0] - q[1]*q[1] + q[2]*q[2] - q[3]*q[3],
-                     2*q[2]* q[3] + 2*q[0]*q[1]],
-          [0, 2*q[1]*q[3] + 2*q[0]*q[2], 2*q[2]* q[3] - 2*q[0]* q[1],
-                     q[0]*q[0] - q[1]*q[1] - q[2]*q[2] + q[3]*q[3]]])
+                     2 * q[2] * q[3] + 2 * q[0] * q[1], 0],
+                    [2 * q[1] * q[3] + 2 * q[0] * q[2],
+                     2 * q[2] * q[3] - 2 * q[0] * q[1],
+                     q[0] * q[0] - q[1] * q[1] - q[2] * q[2] + q[3]*q[3], 0],
+                    [0, 0, 0, 1]])
     if not r:
         return mat
     else:
@@ -130,30 +168,51 @@ def from_matrix(matrix, isprecise=False):
     """Return quaternion from rotation matrix.
     If isprecise is True, the input matrix is assumed to be a precise rotation
     matrix and a faster algorithm is used.
+    :param matrix: a rotation matrix
+    :param axes: string that encodes the order of the axes and
+                 whether rotational or stationary axes should be used
+    :returns: a vector
+    :rtype: (np.ndarray)
+    ..ref: James Diebel
+           "Representing Attitude: Euler Angles, Unit Quaternions, and Rotation
+           Vectors."
+           (2006): p. 15 (6.5)
     """
     r = matrix
 
-    if r[2, 2] > -r[3, 3] or r[1, 1] > -r[2, 2] or r[1, 1] > -r[3, 3]:
-        return (1/2)*np.array([np.sqrt(1 + r[1, 1] + r[2, 2] + r[3,3]),
-           (r[2, 3] - r[3, 2] )/np.sqrt(1 + r[1, 1] + r[2, 2] + r[3,3]),
-           (r[3, 1] - r[1, 3] )/np.sqrt(1 + r[1, 1] + r[2, 2] + r[3, 3]),
-           (r[1, 2] - r[2, 1] )/np.sqrt(1 + r[1, 1] + r[2, 2] + r[3, 3])])
-
-    if r[2, 2] < -r[3, 3] or r[1, 1] > r[2, 2] or r[1, 1] > r[3, 3]:
-        return (1/2)*np.array([(r[2, 3] - r[3, 2] )/np.sqrt(1 + r[1, 1] - r[2, 2] - r[3,3]),
-           np.sqrt(1 + r[1, 1] - r[2, 2] - r[3,3]),
-           (r[1, 2] + r[2, 1] )/np.sqrt(1 + r[1, 1] - r[2, 2] - r[3, 3]),
-           (r[3, 1] + r[1, 3] )/np.sqrt(1 + r[1, 1] - r[2, 2] - r[3, 3])])
-    if r[2, 2] > r[3, 3] or r[1, 1] < r[2, 2] or r[1, 1] < -r[3, 3]:
-        return (1/2)*np.array([(r[3, 1] - r[1, 3] )/np.sqrt(1 - r[1, 1] + r[2, 2] - r[3, 3]),
-           (r[1, 2] + r[2, 1] )/np.sqrt(1 - r[1, 1] + r[2, 2] - r[3, 3]),
-           np.sqrt(1 - r[1, 1] + r[2, 2] - r[3,3] ),
-           (r[2, 3] + r[3, 2] )/np.sqrt(1 - r[1, 1] + r[2, 2] - r[3,3] )])
-    if r[2, 2] < r[3, 3] or r[1, 1] < -r[2, 2] or r[1, 1] < r[3, 3]:
-        return (1/2)*np.array([(r[1, 2] - r[2, 1] )/np.sqrt(1 - r[1, 1] - r[2, 2] + r[3, 3]),
-           (r[3, 1] + r[1, 3] )/np.sqrt(1 - r[1, 1] - r[2, 2] + r[3, 3]),
-           (r[2, 3] + r[3, 2] )/np.sqrt(1 - r[1, 1] - r[2, 2] + r[3,3]),
-           np.sqrt(1 - r[1, 1] - r[2, 2] + r[3,3])])
+    if r[1, 1] > -r[2, 2] or r[0, 0] > -r[1, 1] or r[0, 0] > -r[2, 2]:
+        return (1/2)*np.array([np.sqrt(1 + r[0, 0] + r[1, 1] + r[2, 2]),
+                               (r[1, 2] - r[2, 1]) /
+                               np.sqrt(1 + r[0, 0] + r[1, 1] + r[2, 2]),
+                               (r[2, 0] - r[0, 2]) /
+                               np.sqrt(1 + r[0, 0] + r[1, 1] + r[2, 2]),
+                               (r[0, 1] - r[1, 0]) /
+                               np.sqrt(1 + r[0, 0] + r[1, 1] + r[2, 2])])
+
+    if r[1, 1] < -r[2, 2] or r[0, 0] > r[1, 1] or r[0, 0] > r[2, 2]:
+        return (1/2)*np.array([(r[1, 2] - r[2, 1]) /
+                               np.sqrt(1 + r[0, 0] - r[1, 1] - r[2, 2]),
+                               np.sqrt(1 + r[0, 0] - r[1, 1] - r[2, 2]),
+                               (r[0, 1] + r[1, 0]) /
+                               np.sqrt(1 + r[0, 0] - r[1, 1] - r[2, 2]),
+                               (r[2, 0] + r[0, 2]) /
+                               np.sqrt(1 + r[0, 0] - r[1, 1] - r[2, 2])])
+    if r[1, 1] > r[2, 2] or r[0, 0] < r[1, 1] or r[0, 0] < -r[2, 2]:
+        return (1/2)*np.array([(r[2, 0] - r[0, 2]) /
+                               np.sqrt(1 - r[0, 0] + r[1, 1] - r[2, 2]),
+                               (r[0, 1] + r[1, 0]) /
+                               np.sqrt(1 - r[0, 0] + r[1, 1] - r[2, 2]),
+                               np.sqrt(1 - r[0, 0] + r[1, 1] - r[2, 2]),
+                               (r[1, 2] + r[2, 1]) /
+                               np.sqrt(1 - r[0, 0] + r[1, 1] - r[2, 2])])
+    if r[1, 1] < r[2, 2] or r[0, 0] < -r[1, 1] or r[0, 0] < r[2, 2]:
+        return (1/2)*np.array([(r[0, 1] - r[1, 0]) /
+                               np.sqrt(1 - r[0, 0] - r[1, 1] + r[2, 2]),
+                               (r[2, 0] + r[0, 2]) /
+                               np.sqrt(1 - r[0, 0] - r[1, 1] + r[2, 2]),
+                               (r[1, 2] + r[2, 1]) /
+                               np.sqrt(1 - r[0, 0] - r[1, 1] + r[2, 2]),
+                               np.sqrt(1 - r[0, 0] - r[1, 1] + r[2, 2])])
 
     """
     M = np.array(matrix, dtype=np.float64, copy=False)[:4, :4]
diff --git a/navipy/processing/mcode.py b/navipy/processing/mcode.py
index b90e02f..dcb20d4 100644
--- a/navipy/processing/mcode.py
+++ b/navipy/processing/mcode.py
@@ -10,7 +10,7 @@ from navipy.maths.coordinates\
     import spherical_to_cartesian, cartesian_to_spherical_vectors
 import numpy as np
 import pandas as pd
-from navipy.maths.new_euler import angular_velocity
+from navipy.maths.euler import angular_velocity
 
 
 def spherical_to_Vector(sp):
@@ -431,7 +431,7 @@ def optic_flow(scene, viewing_directions, velocity, distance_channel=3):
         for j, e in enumerate(elevation):
             spline = OFSubroutineSpToVector([a, e])
             M = compose_matrix(angles=[yaw, pitch, roll], translate=None,
-                               perspective=None, axes='rzyx')[:3, :3]
+                               perspective=None, axes=convention)[:3, :3]
             spline = np.dot(M, spline)
             p = distance[j, i]
             if(p == 0):
diff --git a/navipy/processing/test_OpticFlow.py b/navipy/processing/test_OpticFlow.py
index 0dd5c32..81afce8 100644
--- a/navipy/processing/test_OpticFlow.py
+++ b/navipy/processing/test_OpticFlow.py
@@ -291,17 +291,15 @@ class TestCase(unittest.TestCase):
                                          self.viewing_directions,
                                          self.velocity)
 
-        testrof = [[0.00000000e+00, 0.00000000e+00, 0.00000000e+00],
-                   [-8.07793567e-28, 0.00000000e+00, 6.77626358e-21],
-                   [-1.61558713e-27, 0.00000000e+00, 0.00000000e+00]]
-
+        testrof = [[2.88404299e-34, 8.22008280e-20, 1.64376617e-19],
+                   [2.34252646e-05, 4.64310635e-05, 6.94159777e-05],
+                   [9.36297157e-05, 1.38760866e-04, 1.83822856e-04]]
         testhof = [[0.07692013, 0.07690842, 0.07687327],
-                   [0.07692013, 0.07690842, 0.07687327],
-                   [0.07692013, 0.07690842,  0.07687327]]
-
-        testvof = [[2.46536979e-10, 1.34244164e-03, 2.68447412e-03],
-                   [2.46499430e-10, 1.34223718e-03, 2.68406526e-03],
-                   [2.46386795e-10, 1.34162387e-03, 2.68283881e-03]]
+                   [0.0768967, 0.07686158, 0.07680307],
+                   [0.07682644, 0.07676796, 0.07668619]]
+        testvof = [[2.46536984e-10, 1.34244164e-03, 2.68447412e-03],
+                   [1.34203275e-03, 2.68345936e-03, 4.02366094e-03],
+                   [2.68120450e-03, 4.02043495e-03, 5.35762781e-03]]
 
         assert np.all(np.isclose(rof, testrof))
         assert np.all(np.isclose(hof, testhof))
@@ -524,15 +522,15 @@ class TestCase(unittest.TestCase):
         rof, hof, vof = mcode.optic_flow(self.scene, self.viewing_directions,
                                          self.velocity)
 
-        testrof = [[-0.1, -0.1, -0.1],
-                   [-0.09998477, -0.09998477, -0.09998477],
-                   [-0.09993908, -0.09993908, -0.09993908]]
+        testrof = [[0.1, 0.1, 0.1],
+                   [0.09998477, 0.09998477, 0.09998477],
+                   [0.09993908, 0.09993908, 0.09993908]]
         testhof = [[1.02726882e-18, 5.59367784e-12, 1.11856508e-11],
                    [1.02726882e-18, 5.59367784e-12, 1.11856508e-11],
                    [1.02726882e-18, 5.59367784e-12, 1.11856508e-11]]
-        testvof = [[-3.20510339e-10, -3.20461524e-10, -3.20315093e-10],
-                   [1.74524032e-03, 1.74524032e-03, 1.74524032e-03],
-                   [3.48994935e-03, 3.48994935e-03, 3.48994935e-03]]
+        testvof = [[-3.20510352e-10, -3.20461536e-10, -3.20315105e-10],
+                   [-1.74524096e-03, -1.74524096e-03, -1.74524096e-03],
+                   [-3.48994999e-03, -3.48994999e-03, -3.48994999e-03]]
 
         assert np.all(np.isclose(rof, testrof))
         assert np.all(np.isclose(hof, testhof))
@@ -702,15 +700,15 @@ class TestCase(unittest.TestCase):
         self.velocity[self.convention]['dalpha_2'] = droll
         rof, hof, vof = mcode.optic_flow(self.scene, self.viewing_directions,
                                          self.velocity)
-        testrof = [[0.00888717, 0.00888717, 0.00888717],
-                   [0.00715797, 0.00715491, 0.00715237],
-                   [0.00542659, 0.00542046, 0.00541539]]
-        testhof = [[-0.01092516, -0.00919565, -0.00746334],
-                   [-0.01092516, -0.00919565, -0.00746334],
-                   [-0.01092516, -0.00919565, -0.00746334]]
-        testvof = [[-0.09900333, -0.09917892, -0.0993243],
-                   [-0.09914335, -0.09931892, -0.09946428],
-                   [-0.09925318, -0.09942866, -0.09957395]]
+        testrof = [[0.01088051, 0.01088051, 0.01088051],
+                   [0.0126067, 0.01260916, 0.01261109],
+                   [0.01432905, 0.01433397, 0.01433784]]
+        testhof = [[0.00894177, 0.00721257, 0.00548116],
+                   [0.00894177, 0.00721257, 0.00548116],
+                   [0.00894177, 0.00721257, 0.00548116]]
+        testvof = [[0.09900333, 0.09914431, 0.09925508],
+                   [0.09879836, 0.09893931, 0.09905007],
+                   [0.09856329, 0.09870419, 0.09881489]]
 
         assert np.all(np.isclose(rof, testrof))
         assert np.all(np.isclose(hof, testhof))
-- 
GitLab