def test_correct_mesh_orientation():
    sphere_small = ellipsoid(1, 1, 1, levelset=True)

    # Mesh with incorrectly oriented faces which was previously returned from
    # `marching_cubes`, before it guaranteed correct mesh orientation
    verts = np.array([[1., 2., 2.],
                      [2., 2., 1.],
                      [2., 1., 2.],
                      [2., 2., 3.],
                      [2., 3., 2.],
                      [3., 2., 2.]])

    faces = np.array([[0, 1, 2],
                      [2, 0, 3],
                      [1, 0, 4],
                      [4, 0, 3],
                      [1, 2, 5],
                      [2, 3, 5],
                      [1, 4, 5],
                      [5, 4, 3]])

    # Correct mesh orientation - descent
    corrected_faces1 = correct_mesh_orientation(sphere_small, verts, faces,
                                                gradient_direction='descent')
    corrected_faces2 = correct_mesh_orientation(sphere_small, verts, faces,
                                                gradient_direction='ascent')

    # Ensure ascent is opposite of descent for all faces
    assert_array_equal(corrected_faces1, corrected_faces2[:, ::-1])

    # Ensure correct faces have been reversed: 1, 4, and 5
    idx = [1, 4, 5]
    expected = faces.copy()
    expected[idx] = expected[idx, ::-1]
    assert_array_equal(expected, corrected_faces1)

def test_marching_cubes_isotropic():
    ellipsoid_isotropic = ellipsoid(6, 10, 16, levelset=True)
    _, surf = ellipsoid_stats(6, 10, 16)
    verts, faces = marching_cubes(ellipsoid_isotropic, 0.)
    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():
    spacing = (1.0, 10 / 6.0, 16 / 6.0)
    ellipsoid_anisotropic = ellipsoid(6, 10, 16, spacing=spacing, levelset=True)
    _, surf = ellipsoid_stats(6, 10, 16)
    verts, faces = marching_cubes(ellipsoid_anisotropic, 0.0, spacing=spacing)
    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

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 main(select=3, **kwargs):
    """Script main function.

    select: int
        1: Medical data
        2: Blocky data, different every time
        3: Two donuts
        4: Ellipsoid

    """
    import visvis as vv  # noqa: delay import visvis and GUI libraries

    # Create test volume
    if select == 1:
        vol = vv.volread('stent')
        isovalue = kwargs.pop('level', 800)
    elif select == 2:
        vol = vv.aVolume(20, 128)
        isovalue = kwargs.pop('level', 0.2)
    elif select == 3:
        with timer('computing donuts'):
            vol = donuts()
        isovalue = kwargs.pop('level', 0.0)
        # Uncommenting the line below will yield different results for
        # classic MC
        # vol *= -1
    elif select == 4:
        vol = ellipsoid(4, 3, 2, levelset=True)
        isovalue = kwargs.pop('level', 0.0)
    else:
        raise ValueError('invalid selection')

    # Get surface meshes
    with timer('finding surface lewiner'):
        vertices1, faces1 = marching_cubes_lewiner(vol, isovalue, **kwargs)[:2]

    with timer('finding surface classic'):
        vertices2, faces2 = marching_cubes_classic(vol, isovalue, **kwargs)

    # Show
    vv.figure(1)
    vv.clf()
    a1 = vv.subplot(121)
    vv.title('Lewiner')
    m1 = vv.mesh(np.fliplr(vertices1), faces1)
    a2 = vv.subplot(122)
    vv.title('Classic')
    m2 = vv.mesh(np.fliplr(vertices2), faces2)
    a1.camera = a2.camera

    # visvis uses right-hand rule, gradient_direction param uses left-hand rule
    m1.cullFaces = m2.cullFaces = 'front'  # None, front or back

    vv.use().Run()

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_classic(sphere_small, 0)[:2]
    vertices2, faces2 = marching_cubes_lewiner(sphere_small, 0, allow_degenerate=False)[:2]
    vertices3, faces3 = marching_cubes_lewiner(sphere_small, 0, allow_degenerate=False, use_classic=True)[: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 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_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_correct_mesh_orientation():
    sphere_small = ellipsoid(1, 1, 1, levelset=True)
    verts, faces = marching_cubes(sphere_small, 0.)

    # Correct mesh orientation - descent
    corrected_faces1 = correct_mesh_orientation(sphere_small, verts, faces,
                                                gradient_direction='descent')
    corrected_faces2 = correct_mesh_orientation(sphere_small, verts, faces,
                                                gradient_direction='ascent')

    # Ensure ascent is opposite of descent for all faces
    np.testing.assert_array_equal(corrected_faces1, corrected_faces2[:, ::-1])

    # Ensure correct faces have been reversed: 1, 4, and 5
    idx = [1, 4, 5]
    expected = faces.copy()
    expected[idx] = expected[idx, ::-1]
    np.testing.assert_array_equal(expected, corrected_faces1)

    def test_mc_skimage_orig_example(self):
        import numpy as np
        import matplotlib.pyplot as plt
        from mpl_toolkits.mplot3d.art3d import Poly3DCollection

        from skimage import measure
        from skimage.draw import ellipsoid


        # Generate a level set about zero of two identical ellipsoids in 3D
        ellip_base = ellipsoid(6, 10, 16, levelset=True)
        ellip_double = np.concatenate((ellip_base[:-1, ...],
                                       ellip_base[2:, ...]), axis=0)

        # Use marching cubes to obtain the surface mesh of these ellipsoids
        # outs = measure.marching_cubes(ellip_double, 0)
        # verts, faces, normals, values = measure.marching_cubes(ellip_double, 0)
        verts, faces = measure.marching_cubes(ellip_double, 0)

        # Display resulting triangular mesh using Matplotlib. This can also be done
        # with mayavi (see skimage.measure.marching_cubes docstring).
        fig = plt.figure(figsize=(10, 10))
        ax = fig.add_subplot(111, projection='3d')

        # Fancy indexing: `verts[faces]` to generate a collection of triangles
        mesh = Poly3DCollection(verts[faces])
        mesh.set_edgecolor('k')
        ax.add_collection3d(mesh)

        ax.set_xlabel("x-axis: a = 6 per ellipsoid")
        ax.set_ylabel("y-axis: b = 10")
        ax.set_zlabel("z-axis: c = 16")

        ax.set_xlim(0, 24)  # a = 6 (times two for 2nd ellipsoid)
        ax.set_ylim(0, 20)  # b = 10
        ax.set_zlim(0, 32)  # c = 16

        plt.tight_layout()
        plt.show()

def test_marching_cubes_anisotropic():
    spacing = (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_classic(ellipsoid_anisotropic, 0.,
                                          spacing=spacing)
    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_lewiner(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 spacing together with allow_degenerate=False
    marching_cubes_lewiner(ellipsoid_anisotropic, 0, spacing=spacing,
                           allow_degenerate=False)

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_ellipsoid_sign_parameters1():
    ellipsoid(-1, 2, 2)

def test_ellipsoid_sign_parameters3():
    ellipsoid(-3, -2, 2)

def test_ellipsoid_sign_parameters2():
    ellipsoid(0, 2, 2)

    for iz in range(vol.shape[0]):
        for iy in range(vol.shape[1]):
            for ix in range(vol.shape[2]):
                z, y, x = float(iz)*a+b, float(iy)*a+b, float(ix)*a+b
                vol[iz, iy, ix] = ( (
                    (8*x)**2 + (8*y-2)**2 + (8*z)**2 + 16 - 1.85*1.85 ) * ( (8*x)**2 +
                    (8*y-2)**2 + (8*z)**2 + 16 - 1.85*1.85 ) - 64 * ( (8*x)**2 + (8*y-2)**2 )
                    ) * ( ( (8*x)**2 + ((8*y-2)+4)*((8*y-2)+4) + (8*z)**2 + 16 - 1.85*1.85 )
                    * ( (8*x)**2 + ((8*y-2)+4)*((8*y-2)+4) + (8*z)**2 + 16 - 1.85*1.85 ) -
                    64 * ( ((8*y-2)+4)*((8*y-2)+4) + (8*z)**2
                    ) ) + 1025
    # Uncommenting the line below will yield different results for classic MC
    #vol = -vol

elif SELECT == 4:
    vol = ellipsoid(4, 3, 2, levelset=True)
    isovalue = 0

# Get surface meshes
t0 = time.time()
vertices1, faces1, _ = marching_cubes_lewiner(vol, isovalue, gradient_direction=gradient_dir, use_classic=False)
print('finding surface lewiner took %1.0f ms' % (1000*(time.time()-t0)) )

t0 = time.time()
vertices2, faces2, _ = marching_cubes_classic(vol, isovalue, gradient_direction=gradient_dir)
print('finding surface classic took %1.0f ms' % (1000*(time.time()-t0)) )

# Show
vv.figure(1); vv.clf()
a1 = vv.subplot(121); m1 = vv.mesh(np.fliplr(vertices1), faces1)
a2 = vv.subplot(122); m2 = vv.mesh(np.fliplr(vertices2), faces2)

a mesh for regions of bone or bone-like density.

This implementation also works correctly on anisotropic datasets, where the
voxel spacing is not equal for every spatial dimension, through use of the
`spacing` kwarg.

"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.art3d import Poly3DCollection

from skimage import measure
from skimage.draw import ellipsoid

# Generate a level set about zero of two identical ellipsoids in 3D
ellip_base = ellipsoid(6, 10, 16, levelset=True)
ellip_double = np.concatenate((ellip_base[:-1, ...],
                               ellip_base[2:, ...]), axis=0)

# Use marching cubes to obtain the surface mesh of these ellipsoids
verts, faces, normals, values = measure.marching_cubes(ellip_double, 0)

# Display resulting triangular mesh using Matplotlib. This can also be done
# with mayavi or visvis (see skimage.measure.marching_cubes docstring).
fig = plt.figure(figsize=(10, 12))
ax = fig.add_subplot(111, projection='3d')

# Fancy indexing: `verts[faces]` to generate a collection of triangles
mesh = Poly3DCollection(verts[faces])
ax.add_collection3d(mesh)

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

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

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