def make_transform_matrix(center, y=None, z=None):
    """Make 4x4 homogenous transformation matrix.

    Given Nx4 array A where A[:, 4] = 1., the transform matrix M should be
    used with dot(M, A.T).T. Order of operations is 1. translation, 2. align
    `y` x `z` plane to yz-plane 3. align `y` to y-axis.

    Parameters
    ----------
    center : 1x3 array of float
        Coordinate that should be centered after transformation.
    y : None or 1x3 array of float
        Vector that should lie on the y-axis after transformation
    z : None or 1x3 array of float
        Vector that after transformation should lie on yz-plane in direction
        of z-axis.

    Returns
    -------
    4x4 array of float
        4x4 homogenous transformation matrix.
    """
    translate = np.identity(4, dtype=np.float64)
    translate[:3, 3] = -np.asarray(center, dtype=np.float64)
    if y is not None:
        y = np.atleast_2d(y)
        if z is None:
            rotate = np.identity(4, dtype=np.float64)
            rotate[:3, :3] = make_rotation_matrix(y, Y_AXIS)
        else:
            z = np.atleast_2d(z)
            rotate_norm = np.identity(4, dtype=np.float64)
            x_unit = as_unit(np.cross(y, z))
            rotate_norm[:3, :3] = make_rotation_matrix(x_unit, X_AXIS)
            new_y = np.dot(rotate_norm[:3, :3], y.flatten())
            rotate_y = np.identity(4, dtype=np.float64)
            rotate_y[:3, :3] = make_rotation_matrix(new_y.flatten(), Y_AXIS)
            rotate = np.dot(rotate_y, rotate_norm)
        transform = np.dot(rotate, translate)
    else:
        transform = translate
    return transform