skimage.draw.ellipsoid

def test_ellipsoid_sign_parameters1():
    with testing.raises(ValueError):
        ellipsoid(-1, 2, 2)


def test_ellipsoid_sign_parameters2():

def test_ellipsoid_sign_parameters2():
    with testing.raises(ValueError):
        ellipsoid(0, 2, 2)


def test_ellipsoid_sign_parameters3():

def test_ellipsoid_sign_parameters3():
    with testing.raises(ValueError):
        ellipsoid(-3, -2, 2)


def test_ellipsoid_bool():

def im3d():
    r = 10
    pad = 10
    im3 = draw.ellipsoid(r, r, r)
    im3 = np.pad(im3, pad, mode='constant').astype(np.uint8)
    return im3


def test_structure_tensor():

def test_marching_cubes_isotropic():
    ellipsoid_isotropic = ellipsoid(6, 10, 16, levelset=True)
    _, surf = ellipsoid_stats(6, 10, 16)

    # Classic
    verts, faces = marching_cubes(ellipsoid_isotropic, 0., method='_lorensen')
    surf_calc = mesh_surface_area(verts, faces)
    # Test within 1% tolerance for isotropic. Will always underestimate.
    assert surf > surf_calc and surf_calc > surf * 0.99

    # Lewiner
    verts, faces = marching_cubes(ellipsoid_isotropic, 0.)[:2]
    surf_calc = mesh_surface_area(verts, faces)
    # Test within 1% tolerance for isotropic. Will always underestimate.
    assert surf > surf_calc and surf_calc > surf * 0.99


def test_marching_cubes_anisotropic():

def test_both_algs_same_result_ellipse():
    # Performing this test on data that does not have ambiguities

    sphere_small = ellipsoid(1, 1, 1, levelset=True)

    vertices1, faces1 = marching_cubes(sphere_small, 0, method='_lorensen')[:2]
    vertices2, faces2 = marching_cubes(sphere_small, 0,
                                       allow_degenerate=False)[:2]
    vertices3, faces3 = marching_cubes(sphere_small, 0,
                                       allow_degenerate=False,
                                       method='lorensen')[:2]

    # Order is different, best we can do is test equal shape and same
    # vertices present
    assert _same_mesh(vertices1, faces1, vertices2, faces2)
    assert _same_mesh(vertices1, faces1, vertices3, faces3)


def _same_mesh(vertices1, faces1, vertices2, faces2, tol=1e-10):

def test_masked_marching_cubes():

    ellipsoid_scalar = ellipsoid(6, 10, 16, levelset=True)
    mask = np.ones_like(ellipsoid_scalar, dtype=bool)
    mask[:10, :, :] = False
    mask[:, :, 20:] = False
    ver, faces, _, _ = marching_cubes(ellipsoid_scalar, 0, mask=mask)
    area = mesh_surface_area(ver, faces)

    np.testing.assert_allclose(area, 299.56878662109375, rtol=.01)


def test_masked_marching_cubes_empty():

def test_masked_marching_cubes_empty():
    ellipsoid_scalar = ellipsoid(6, 10, 16, levelset=True)
    mask = np.array([])
    with pytest.raises(ValueError):
        _ = marching_cubes(ellipsoid_scalar, 0, mask=mask)


def test_masked_marching_cubes_old_lewiner():

def test_masked_marching_cubes_old_lewiner():
    ellipsoid_scalar = ellipsoid(6, 10, 16, levelset=True)
    mask = np.array([])
    with pytest.raises(NotImplementedError):
        _ = marching_cubes(ellipsoid_scalar, 0, mask=mask, method='_lorensen')


def test_masked_marching_cubes_all_true():

def test_masked_marching_cubes_all_true():
    ellipsoid_scalar = ellipsoid(6, 10, 16, levelset=True)
    mask = np.ones_like(ellipsoid_scalar, dtype=bool)
    ver_m, faces_m, _, _ = marching_cubes(ellipsoid_scalar, 0, mask=mask)
    ver, faces, _, _ = marching_cubes(ellipsoid_scalar, 0, mask=mask)
    np.testing.assert_allclose(ver_m, ver, rtol=.00001)
    np.testing.assert_allclose(faces_m, faces, rtol=.00001)

def test_moments_normalized_3d():
    image = draw.ellipsoid(1, 1, 10)
    mu_image = moments_central(image)
    nu = moments_normalized(mu_image)
    assert nu[0, 0, 2] > nu[0, 2, 0]
    assert_almost_equal(nu[0, 2, 0], nu[2, 0, 0])

    coords = np.where(image)
    mu_coords = moments_coords_central(coords)
    assert_almost_equal(mu_coords, mu_image)


def test_moments_normalized_invalid():

def test_ellipsoid_bool():
    test = ellipsoid(2, 2, 2)[1:-1, 1:-1, 1:-1]
    test_anisotropic = ellipsoid(2, 2, 4, spacing=(1., 1., 2.))
    test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1]

    expected = np.array([[[0, 0, 0, 0, 0],
                          [0, 0, 0, 0, 0],
                          [0, 0, 1, 0, 0],
                          [0, 0, 0, 0, 0],
                          [0, 0, 0, 0, 0]],

                         [[0, 0, 0, 0, 0],
                          [0, 1, 1, 1, 0],
                          [0, 1, 1, 1, 0],
                          [0, 1, 1, 1, 0],
                          [0, 0, 0, 0, 0]],

                         [[0, 0, 1, 0, 0],
                          [0, 1, 1, 1, 0],
                          [1, 1, 1, 1, 1],
                          [0, 1, 1, 1, 0],
                          [0, 0, 1, 0, 0]],

                         [[0, 0, 0, 0, 0],
                          [0, 1, 1, 1, 0],
                          [0, 1, 1, 1, 0],
                          [0, 1, 1, 1, 0],
                          [0, 0, 0, 0, 0]],

                         [[0, 0, 0, 0, 0],
                          [0, 0, 0, 0, 0],
                          [0, 0, 1, 0, 0],
                          [0, 0, 0, 0, 0],
                          [0, 0, 0, 0, 0]]])

    assert_array_equal(test, expected.astype(bool))
    assert_array_equal(test_anisotropic, expected.astype(bool))


def test_ellipsoid_levelset():

def test_ellipsoid_levelset():
    test = ellipsoid(2, 2, 2, levelset=True)[1:-1, 1:-1, 1:-1]
    test_anisotropic = ellipsoid(2, 2, 4, spacing=(1., 1., 2.),
                                 levelset=True)
    test_anisotropic = test_anisotropic[1:-1, 1:-1, 1:-1]

    expected = np.array([[[ 2.  ,  1.25,  1.  ,  1.25,  2.  ],
                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
                          [ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
                          [ 2.  ,  1.25,  1.  ,  1.25,  2.  ]],

                         [[ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
                          [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
                          [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
                          [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25]],

                         [[ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
                          [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
                          [ 0.  , -0.75, -1.  , -0.75,  0.  ],
                          [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
                          [ 1.  ,  0.25,  0.  ,  0.25,  1.  ]],

                         [[ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
                          [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
                          [ 0.25, -0.5 , -0.75, -0.5 ,  0.25],
                          [ 0.5 , -0.25, -0.5 , -0.25,  0.5 ],
                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25]],

                         [[ 2.  ,  1.25,  1.  ,  1.25,  2.  ],
                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
                          [ 1.  ,  0.25,  0.  ,  0.25,  1.  ],
                          [ 1.25,  0.5 ,  0.25,  0.5 ,  1.25],
                          [ 2.  ,  1.25,  1.  ,  1.25,  2.  ]]])

    assert_allclose(test, expected)
    assert_allclose(test_anisotropic, expected)


def test_ellipsoid_stats():

def test_marching_cubes_anisotropic():
    # test spacing as numpy array (and not just tuple)
    spacing = np.array([1., 10 / 6., 16 / 6.])
    ellipsoid_anisotropic = ellipsoid(6, 10, 16, spacing=spacing,
                                      levelset=True)
    _, surf = ellipsoid_stats(6, 10, 16)

    # Classic
    verts, faces = marching_cubes(ellipsoid_anisotropic, 0.,
                                  spacing=spacing, method='_lorensen')
    surf_calc = mesh_surface_area(verts, faces)
    # Test within 1.5% tolerance for anisotropic. Will always underestimate.
    assert surf > surf_calc and surf_calc > surf * 0.985

    # Lewiner
    verts, faces = marching_cubes(ellipsoid_anisotropic, 0.,
                                  spacing=spacing)[:2]
    surf_calc = mesh_surface_area(verts, faces)
    # Test within 1.5% tolerance for anisotropic. Will always underestimate.
    assert surf > surf_calc and surf_calc > surf * 0.985

    # Test marching cube with mask
    with pytest.raises(ValueError):
        verts, faces = marching_cubes(
            ellipsoid_anisotropic, 0., spacing=spacing,
            mask=np.array([]))[:2]

    # Test spacing together with allow_degenerate=False
    marching_cubes(ellipsoid_anisotropic, 0, spacing=spacing,
                   allow_degenerate=False)


def test_invalid_input():

def test_inertia_tensor_3d():
    image = draw.ellipsoid(10, 5, 3)
    T0 = inertia_tensor(image)
    eig0, V0 = np.linalg.eig(T0)
    # principal axis of ellipse = eigenvector of smallest eigenvalue
    v0 = V0[:, np.argmin(eig0)]

    assert np.allclose(v0, [1, 0, 0]) or np.allclose(-v0, [1, 0, 0])

    imrot = ndi.rotate(image.astype(float), 30, axes=(0, 1), order=1)
    Tr = inertia_tensor(imrot)
    eigr, Vr = np.linalg.eig(Tr)
    vr = Vr[:, np.argmin(eigr)]

    # Check that axis has rotated by expected amount
    pi, cos, sin = np.pi, np.cos, np.sin
    R = np.array([[ cos(pi/6), -sin(pi/6), 0],
                  [ sin(pi/6),  cos(pi/6), 0],
                  [         0,          0, 1]])
    expected_vr = R @ v0
    assert (np.allclose(vr, expected_vr, atol=1e-3, rtol=0.01) or
            np.allclose(-vr, expected_vr, atol=1e-3, rtol=0.01))


def test_inertia_tensor_eigvals():

def volume_ellipsoid_spacing(a, b, c, spacing):
    """

    :param a: major semi-axis of an ellipsoid
    :param b: least semi-axis of an ellipsoid
    :param c:  minor semi-axis of an ellipsoid
    :param spacing: spacing of a grid, tuple like e.g. (1,  1, 1)
    :return: volume of an ellipsoid in ml, taking spacing into account
    """
    ellipsoid_array = ellipsoid(a, b, c, spacing)
    ellipsoid_non_zero = ellipsoid_array.nonzero()
    num_voxels = len(list(zip(ellipsoid_non_zero[0],
                              ellipsoid_non_zero[1],
                              ellipsoid_non_zero[2])))
    volume_object_ml = (num_voxels * spacing[0] * spacing[1] * spacing[2]) / 1000

    return volume_object_ml


def get_surface_points(dcm_img):