diff --git a/fury/transform.py b/fury/transform.py new file mode 100644 index 000000000..d8818493c --- /dev/null +++ b/fury/transform.py @@ -0,0 +1,400 @@ +import math + +import numpy as np +from scipy.spatial.transform import Rotation as Rot # type: ignore + +from fury.decorators import warn_on_args_to_kwargs + +# axis sequences for Euler angles +_NEXT_AXIS = [1, 2, 0, 1] + +# map axes strings to/from tuples of inner axis, parity, repetition, frame +_AXES2TUPLE = { + "sxyz": (0, 0, 0, 0), + "sxyx": (0, 0, 1, 0), + "sxzy": (0, 1, 0, 0), + "sxzx": (0, 1, 1, 0), + "syzx": (1, 0, 0, 0), + "syzy": (1, 0, 1, 0), + "syxz": (1, 1, 0, 0), + "syxy": (1, 1, 1, 0), + "szxy": (2, 0, 0, 0), + "szxz": (2, 0, 1, 0), + "szyx": (2, 1, 0, 0), + "szyz": (2, 1, 1, 0), + "rzyx": (0, 0, 0, 1), + "rxyx": (0, 0, 1, 1), + "ryzx": (0, 1, 0, 1), + "rxzx": (0, 1, 1, 1), + "rxzy": (1, 0, 0, 1), + "ryzy": (1, 0, 1, 1), + "rzxy": (1, 1, 0, 1), + "ryxy": (1, 1, 1, 1), + "ryxz": (2, 0, 0, 1), + "rzxz": (2, 0, 1, 1), + "rxyz": (2, 1, 0, 1), + "rzyz": (2, 1, 1, 1), +} + +_TUPLE2AXES = {v: k for k, v in _AXES2TUPLE.items()} + + +@warn_on_args_to_kwargs() +def euler_matrix(ai, aj, ak, *, axes="sxyz"): + """Return homogeneous rotation matrix from Euler angles and axis sequence. + + Code modified from the work of Christoph Gohlke link provided here + http://www.lfd.uci.edu/~gohlke/code/transformations.py.html + + Parameters + ---------- + ai, aj, ak : Euler's roll, pitch and yaw angles + axes : One of 24 axis sequences as string or encoded tuple + + Returns + ------- + matrix : ndarray (4, 4) + + Code modified from the work of Christoph Gohlke link provided here + http://www.lfd.uci.edu/~gohlke/code/transformations.py.html + + Examples + -------- + >>> import numpy + >>> R = euler_matrix(1, 2, 3, 'syxz') + >>> numpy.allclose(numpy.sum(R[0]), -1.34786452) + True + >>> R = euler_matrix(1, 2, 3, (0, 1, 0, 1)) + >>> numpy.allclose(numpy.sum(R[0]), -0.383436184) + True + >>> ai, aj, ak = (4.0*math.pi) * (numpy.random.random(3) - 0.5) + >>> for axes in _AXES2TUPLE.keys(): + ... _ = euler_matrix(ai, aj, ak, axes) + >>> for axes in _TUPLE2AXES.keys(): + ... _ = euler_matrix(ai, aj, ak, axes) + + """ + try: + firstaxis, parity, repetition, frame = _AXES2TUPLE[axes] + except (AttributeError, KeyError): + firstaxis, parity, repetition, frame = axes + + i = firstaxis + j = _NEXT_AXIS[i + parity] + k = _NEXT_AXIS[i - parity + 1] + + if frame: + ai, ak = ak, ai + if parity: + ai, aj, ak = -ai, -aj, -ak + + si, sj, sk = math.sin(ai), math.sin(aj), math.sin(ak) + ci, cj, ck = math.cos(ai), math.cos(aj), math.cos(ak) + cc, cs = ci * ck, ci * sk + sc, ss = si * ck, si * sk + + M = np.identity(4) + if repetition: + M[i, i] = cj + M[i, j] = sj * si + M[i, k] = sj * ci + M[j, i] = sj * sk + M[j, j] = -cj * ss + cc + M[j, k] = -cj * cs - sc + M[k, i] = -sj * ck + M[k, j] = cj * sc + cs + M[k, k] = cj * cc - ss + else: + M[i, i] = cj * ck + M[i, j] = sj * sc - cs + M[i, k] = sj * cc + ss + M[j, i] = cj * sk + M[j, j] = sj * ss + cc + M[j, k] = sj * cs - sc + M[k, i] = -sj + M[k, j] = cj * si + M[k, k] = cj * ci + return M + + +def sphere2cart(r, theta, phi): + """Spherical to Cartesian coordinates. + + This is the standard physics convention where `theta` is the + inclination (polar) angle, and `phi` is the azimuth angle. + + Imagine a sphere with center (0,0,0). Orient it with the z axis + running south-north, the y axis running west-east and the x axis + from posterior to anterior. `theta` (the inclination angle) is the + angle to rotate from the z-axis (the zenith) around the y-axis, + towards the x axis. Thus the rotation is counter-clockwise from the + point of view of positive y. `phi` (azimuth) gives the angle of + rotation around the z-axis towards the y axis. The rotation is + counter-clockwise from the point of view of positive z. + + Equivalently, given a point P on the sphere, with coordinates x, y, + z, `theta` is the angle between P and the z-axis, and `phi` is + the angle between the projection of P onto the XY plane, and the X + axis. + + Geographical nomenclature designates theta as 'co-latitude', and phi + as 'longitude' + + Parameters + ---------- + r : array_like + radius + theta : array_like + inclination or polar angle + phi : array_like + azimuth angle + + Returns + ------- + x : array + x coordinate(s) in Cartesian space + y : array + y coordinate(s) in Cartesian space + z : array + z coordinate + + Notes + ----- + See these pages: + + * http://en.wikipedia.org/wiki/Spherical_coordinate_system + * http://mathworld.wolfram.com/SphericalCoordinates.html + + for excellent discussion of the many different conventions + possible. Here we use the physics conventions, used in the + wikipedia page. + + Derivations of the formulae are simple. Consider a vector x, y, z of + length r (norm of x, y, z). The inclination angle (theta) can be + found from: cos(theta) == z / r -> z == r * cos(theta). This gives + the hypotenuse of the projection onto the XY plane, which we will + call Q. Q == r*sin(theta). Now x / Q == cos(phi) -> x == r * + sin(theta) * cos(phi) and so on. + + We have deliberately named this function ``sphere2cart`` rather than + ``sph2cart`` to distinguish it from the Matlab function of that + name, because the Matlab function uses an unusual convention for the + angles that we did not want to replicate. The Matlab function is + trivial to implement with the formulae given in the Matlab help. + + """ + sin_theta = np.sin(theta) + x = r * np.cos(phi) * sin_theta + y = r * np.sin(phi) * sin_theta + z = r * np.cos(theta) + x, y, z = np.broadcast_arrays(x, y, z) + return x, y, z + + +def cart2sphere(x, y, z): + r"""Return angles for Cartesian 3D coordinates `x`, `y`, and `z`. + + See doc for ``sphere2cart`` for angle conventions and derivation + of the formulae. + + $0\le\theta\mathrm{(theta)}\le\pi$ and $-\pi\le\phi\mathrm{(phi)}\le\pi$ + + Parameters + ---------- + x : array_like + x coordinate in Cartesian space + y : array_like + y coordinate in Cartesian space + z : array_like + z coordinate + + Returns + ------- + r : array + radius + theta : array + inclination (polar) angle + phi : array + azimuth angle + + """ + r = np.sqrt(x * x + y * y + z * z) + theta = np.arccos(np.divide(z, r, where=r > 0)) + theta = np.where(r > 0, theta, 0.0) + phi = np.arctan2(y, x) + r, theta, phi = np.broadcast_arrays(r, theta, phi) + return r, theta, phi + + +def translate(translation): + """Return transformation matrix for translation array. + + Parameters + ---------- + translation : ndarray + translation in x, y and z directions. + + Returns + ------- + translation : ndarray (4, 4) + Numpy array of shape 4,4 containing translation parameter in the last + column of the matrix. + + Examples + -------- + >>> import numpy as np; import fury + >>> tran = np.array([0.3, 0.2, 0.25]) + >>> transform = fury.transform.translate(tran) + >>> transform + array([[1. , 0. , 0. , 0.3 ], + [0. , 1. , 0. , 0.2 ], + [0. , 0. , 1. , 0.25], + [0. , 0. , 0. , 1. ]]) + + """ + iden = np.identity(4) + translation = np.append(translation, 0).reshape(-1, 1) + + t = np.array( + [[0, 0, 0, 1], [0, 0, 0, 1], [0, 0, 0, 1], [0, 0, 0, 1]], + np.float32, + ) + translation = np.multiply(t, translation) + translation = np.add(iden, translation) + + return translation + + +def rotate(quat): + """Return transformation matrix for rotation quaternion. + + Parameters + ---------- + quat : ndarray (4, ) + rotation quaternion. + + Returns + ------- + rotation_mat : ndarray (4, 4) + Transformation matrix of shape (4, 4) to rotate a vector. + + Examples + -------- + >>> import numpy as np; import fury + >>> quat = np.array([0.259, 0.0, 0.0, 0.966]) + >>> rotation = fury.transform.rotate(quat) + >>> rotation + array([[ 1. , 0. , 0. , 0. ], + [ 0. , 0.86586979, -0.50026944, 0. ], + [ 0. , 0.50026944, 0.86586979, 0. ], + [ 0. , 0. , 0. , 1. ]]) + + """ + iden = np.identity(3) + rotation_mat = Rot.from_quat(quat).as_matrix() + + iden = np.append(iden, [[0, 0, 0]]).reshape(-1, 3) + + rotation_mat = np.dot(iden, rotation_mat) + iden = np.array([[0, 0, 0, 1]]).reshape(-1, 1) + + rotation_mat = np.concatenate((rotation_mat, iden), axis=1) + return rotation_mat + + +def scale(scales): + """Return transformation matrix for scales array. + + Parameters + ---------- + scales : ndarray + scales in x, y and z directions. + + Returns + ------- + scale_mat : ndarray (4, 4) + Numpy array of shape 4,4 containing elements of scale matrix along + the diagonal. + + Examples + -------- + >>> import numpy as np; import fury + >>> scales = np.array([2.0, 1.0, 0.5]) + >>> transform = fury.transform.scale(scales) + >>> transform + array([[2. , 0. , 0. , 0. ], + [0. , 1. , 0. , 0. ], + [0. , 0. , 0.5, 0. ], + [0. , 0. , 0. , 1. ]]) + + """ + scale_mat = np.identity(4) + scales = np.append(scales, [1]) + + for i in range(len(scales)): + scale_mat[i][i] = scales[i] + + return scale_mat + + +def apply_transformation(vertices, transformation): + """Multiplying transformation matrix with vertices + + Parameters + ---------- + vertices : ndarray (n, 3) + vertices of the mesh + transformation : ndarray (4, 4) + transformation matrix + + Returns + ------- + vertices : ndarray (n, 3) + transformed vertices of the mesh + + """ + shape = vertices.shape + temp = np.full((shape[0], 1), 1) + vertices = np.concatenate((vertices, temp), axis=1) + + vertices = np.dot(transformation, vertices.T) + vertices = vertices.T + vertices = vertices[:, : shape[1]] + + return vertices + + +def transform_from_matrix(matrix): + """Returns translation, rotation and scale arrays from transformation + matrix. + + Parameters + ---------- + matrix : ndarray (4, 4) + the transformation matrix of shape 4*4 + + Returns + ------- + translate : ndarray (3, ) + translation component from the transformation matrix + rotate : ndarray (4, ) + rotation component from the transformation matrix + scale : ndarray (3, ) + scale component from the transformation matrix. + + """ + translate = matrix[:, -1:].reshape((-1,))[:-1] + + temp = matrix[:, :3][:3] + sx = np.linalg.norm(temp[:, :1]) + sy = np.linalg.norm(temp[:, 1:-1]) + sz = np.linalg.norm(temp[:, -1:]) + scale = np.array([sx, sy, sz]) + + rot_matrix = temp / scale[None, :] + rotation = Rot.from_matrix(rot_matrix) + rot_vec = rotation.as_rotvec() + angle = np.linalg.norm(rot_vec) + rotation = [np.rad2deg(angle), *rot_vec] + + return translate, rotation, scale diff --git a/fury/utils.py b/fury/utils.py new file mode 100644 index 000000000..87b6b189c --- /dev/null +++ b/fury/utils.py @@ -0,0 +1,305 @@ +import numpy as np +from scipy.ndimage import map_coordinates + + +def map_coordinates_3d_4d(input_array, indices): + """Evaluate input_array at the given indices using trilinear interpolation. + + Parameters + ---------- + input_array : ndarray, + 3D or 4D array + indices : ndarray + + Returns + ------- + output : ndarray + 1D or 2D array + + """ + if input_array.ndim <= 2 or input_array.ndim >= 5: + raise ValueError("Input array can only be 3d or 4d") + + if input_array.ndim == 3: + return map_coordinates(input_array, indices.T, order=1) + + if input_array.ndim == 4: + values_4d = [] + for i in range(input_array.shape[-1]): + values_tmp = map_coordinates( + input_array[..., i], indices.T, order=1) + values_4d.append(values_tmp) + return np.ascontiguousarray(np.array(values_4d).T) + + +def apply_affine(aff, pts): + """Apply affine matrix `aff` to points `pts`. + + Returns result of application of `aff` to the *right* of `pts`. The + coordinate dimension of `pts` should be the last. + For the 3D case, `aff` will be shape (4,4) and `pts` will have final axis + length 3 - maybe it will just be N by 3. The return value is the + transformed points, in this case:: + res = np.dot(aff[:3,:3], pts.T) + aff[:3,3:4] + transformed_pts = res.T + This routine is more general than 3D, in that `aff` can have any shape + (N,N), and `pts` can have any shape, as long as the last dimension is for + the coordinates, and is therefore length N-1. + + Parameters + ---------- + aff : (N, N) array-like + + Homogeneous affine, for 3D points, will be 4 by 4. Contrary to first + appearance, the affine will be applied on the left of `pts`. + pts : (..., N-1) array-like + Points, where the last dimension contains the coordinates of each + point. For 3D, the last dimension will be length 3. + + Returns + ------- + transformed_pts : (..., N-1) array + transformed points + + Notes + ----- + Copied from nibabel to remove dependency. + + Examples + -------- + >>> aff = np.array([[0,2,0,10],[3,0,0,11],[0,0,4,12],[0,0,0,1]]) + >>> pts = np.array([[1,2,3],[2,3,4],[4,5,6],[6,7,8]]) + >>> apply_affine(aff, pts) #doctest: +ELLIPSIS + array([[14, 14, 24], + [16, 17, 28], + [20, 23, 36], + [24, 29, 44]]...) + >>> # Just to show that in the simple 3D case, it is equivalent to: + >>> (np.dot(aff[:3,:3], pts.T) + aff[:3,3:4]).T #doctest: +ELLIPSIS + array([[14, 14, 24], + [16, 17, 28], + [20, 23, 36], + [24, 29, 44]]...) + >>> # But `pts` can be a more complicated shape: + >>> pts = pts.reshape((2,2,3)) + >>> apply_affine(aff, pts) #doctest: +ELLIPSIS + array([[[14, 14, 24], + [16, 17, 28]], + + [[20, 23, 36], + [24, 29, 44]]]...) + + """ + aff = np.asarray(aff) + pts = np.asarray(pts) + shape = pts.shape + pts = pts.reshape((-1, shape[-1])) + # rzs == rotations, zooms, shears + rzs = aff[:-1, :-1] + trans = aff[:-1, -1] + res = np.dot(pts, rzs.T) + trans[None, :] + return res.reshape(shape) + + +def asbytes(s): + if isinstance(s, bytes): + return s + return s.encode("latin1") + + +def get_grid_cells_position(shapes, *, aspect_ratio=16 / 9.0, dim=None): + """Construct a XY-grid based on the cells content shape. + + This function generates the coordinates of every grid cell. The width and + height of every cell correspond to the largest width and the largest height + respectively. The grid dimensions will automatically be adjusted to respect + the given aspect ratio unless they are explicitly specified. + + The grid follows a row-major order with the top left corner being at + coordinates (0,0,0) and the bottom right corner being at coordinates + (nb_cols*cell_width, -nb_rows*cell_height, 0). Note that the X increases + while the Y decreases. + + Parameters + ---------- + shapes : list of tuple of int + The shape (width, height) of every cell content. + aspect_ratio : float (optional) + Aspect ratio of the grid (width/height). Default: 16:9. + dim : tuple of int (optional) + Dimension (nb_rows, nb_cols) of the grid, if provided. + + Returns + ------- + ndarray + 3D coordinates of every grid cell. + + """ + cell_shape = np.r_[np.max(shapes, axis=0), 0] + cell_aspect_ratio = cell_shape[0] / cell_shape[1] + + count = len(shapes) + if dim is None: + # Compute the number of rows and columns. + n_cols = np.ceil(np.sqrt(count * aspect_ratio / cell_aspect_ratio)) + n_rows = np.ceil(count / n_cols) + else: + n_rows, n_cols = dim + + if n_cols * n_rows < count: + msg = "Size is too small, it cannot contain at least {} elements." + raise ValueError(msg.format(count)) + + # Use indexing="xy" so the cells are in row-major (C-order). Also, + # the Y coordinates are negative so the cells are order from top to bottom. + X, Y, Z = np.meshgrid(np.arange(n_cols), - + np.arange(n_rows), [0], indexing="xy") + return cell_shape * np.array([X.flatten(), Y.flatten(), Z.flatten()]).T + + +def normalize_v3(arr): + """Normalize a numpy array of 3 component vectors shape=(N, 3). + + Parameters + ---------- + array : ndarray + Shape (N, 3) + + Returns + ------- + norm_array + + """ + lens = np.sqrt(arr[:, 0] ** 2 + arr[:, 1] ** 2 + arr[:, 2] ** 2) + arr[:, 0] /= lens + arr[:, 1] /= lens + arr[:, 2] /= lens + return arr + + +def normals_from_v_f(vertices, faces): + """Calculate normals from vertices and faces. + + Parameters + ---------- + verices : ndarray + faces : ndarray + + Returns + ------- + normals : ndarray + Shape same as vertices + + """ + norm = np.zeros(vertices.shape, dtype=vertices.dtype) + tris = vertices[faces] + n = np.cross(tris[::, 1] - tris[::, 0], tris[::, 2] - tris[::, 0]) + normalize_v3(n) + norm[faces[:, 0]] += n + norm[faces[:, 1]] += n + norm[faces[:, 2]] += n + normalize_v3(norm) + return norm + + +def tangents_from_direction_of_anisotropy(normals, direction): + """Calculate tangents from normals and a 3D vector representing the + direction of anisotropy. + + Parameters + ---------- + normals : normals, represented as 2D ndarrays (Nx3) (one per vertex) + direction : tuple (3,) or array (3,) + + Returns + ------- + output : array (N, 3) + Tangents, represented as 2D ndarrays (Nx3). + + """ + tangents = np.cross(normals, direction) + binormals = normalize_v3(np.cross(normals, tangents)) + return normalize_v3(np.cross(normals, binormals)) + + +def triangle_order(vertices, faces): + """Determine the winding order of a given set of vertices and a triangle. + + Parameters + ---------- + vertices : ndarray + array of vertices making up a shape + faces : ndarray + array of triangles + + Returns + ------- + order : int + If the order is counter clockwise (CCW), returns True. + Otherwise, returns False. + + """ + v1 = vertices[faces[0]] + v2 = vertices[faces[1]] + v3 = vertices[faces[2]] + + # https://stackoverflow.com/questions/40454789/computing-face-normals-and-winding + m_orient = np.ones((4, 4)) + m_orient[0, :3] = v1 + m_orient[1, :3] = v2 + m_orient[2, :3] = v3 + m_orient[3, :3] = 0 + + val = np.linalg.det(m_orient) + + return bool(val > 0) + + +def change_vertices_order(triangle): + """Change the vertices order of a given triangle. + + Parameters + ---------- + triangle : ndarray, shape(1, 3) + array of 3 vertices making up a triangle + + Returns + ------- + new_triangle : ndarray, shape(1, 3) + new array of vertices making up a triangle in the opposite winding + order of the given parameter + + """ + return np.array([triangle[2], triangle[1], triangle[0]]) + + +def fix_winding_order(vertices, triangles, *, clockwise=False): + """Return corrected triangles. + + Given an ordering of the triangle's three vertices, a triangle can appear + to have a clockwise winding or counter-clockwise winding. + Clockwise means that the three vertices, in order, rotate clockwise around + the triangle's center. + + Parameters + ---------- + vertices : ndarray + array of vertices corresponding to a shape + triangles : ndarray + array of triangles corresponding to a shape + clockwise : bool + triangle order type: clockwise (default) or counter-clockwise. + + Returns + ------- + corrected_triangles : ndarray + The corrected order of the vert parameter + + """ + corrected_triangles = triangles.copy() + desired_order = clockwise + for nb, face in enumerate(triangles): + current_order = triangle_order(vertices, face) + if desired_order != current_order: + corrected_triangles[nb] = change_vertices_order(face) + return corrected_triangles