def apply_color_matrix(im, m):
    """Transform a color array with the given 3x3 matrix.

    Parameters
    ==========
    im : array of shape (...,3)
        Can be an image or a 1D array, as long as the last
        dimension is a 3-channels color.

    m : array of shape (3,3)
        Color matrix to apply.

    Returns
    =======
    im : array of shape (...,3)
        Output array, where each input color vector was multiplied by m.
    """
    # Another option is np.einsum('ij, ...j', m, im), but it can be much
    # slower, especially on on float32 types because the matrix multiplication
    # is heavily optimized.
    # So the matmul is generally (much) faster, but we need to take the
    # transpose of m as it gets applied on the right side. Indeed for each
    # column color vector v we wanted $v' = m . v$ . To flip the side we can
    # use $m . v = (v^T . m^T)^T$ . The transposes on the 1d vector are implicit
    # and can be ignored, so we just need to compute $v . m^T$. This is what
    # numpy matmul will do for all the vectors thanks to its broadcasting rules
    # that pick the last 2 dimensions of each array, so it will actually compute
    # matrix multiplications of shape (M,3) x (3,3) with M the penultimate dimension
    # of m. That will write a matrix of shape (M,3) with each row storing the
    # result of $v' = v . M^T$.
    return im @ m.T