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