Python sphere.HemiSphere类代码示例

def test_mirror():
    verts = [[0, 0, 1],
             [0, 1, 0],
             [1, 0, 0],
             [-1, -1, -1]]
    verts = np.array(verts, 'float')
    verts = verts / np.sqrt((verts * verts).sum(-1)[:, None])
    faces = [[0, 1, 3],
             [0, 2, 3],
             [1, 2, 3]]

    h = HemiSphere(xyz=verts, faces=faces)
    s = h.mirror()

    nt.assert_equal(len(s.vertices), 8)
    nt.assert_equal(len(s.faces), 6)
    verts = s.vertices

    def _angle(a, b):
        return np.arccos(np.dot(a, b))

    for triangle in s.faces:
        a, b, c = triangle
        nt.assert_(_angle(verts[a], verts[b]) <= np.pi/2)
        nt.assert_(_angle(verts[a], verts[c]) <= np.pi/2)
        nt.assert_(_angle(verts[b], verts[c]) <= np.pi/2)

def test_hemisphere_subdivide():

    def flip(vertices):
        x, y, z = vertices.T
        f = (z < 0) | ((z == 0) & (y < 0)) | ((z == 0) & (y == 0) & (x < 0))
        return 1 - 2*f[:, None]

    decimals = 6
    # Test HemiSphere.subdivide
    # Create a hemisphere by dividing a hemi-icosahedron
    hemi1 = HemiSphere.from_sphere(unit_icosahedron).subdivide(4)
    vertices1 = np.round(hemi1.vertices, decimals)
    vertices1 *= flip(vertices1)
    order = np.lexsort(vertices1.T)
    vertices1 = vertices1[order]

    # Create a hemisphere from a subdivided sphere
    sphere = unit_icosahedron.subdivide(4)
    hemi2 = HemiSphere.from_sphere(sphere)
    vertices2 = np.round(hemi2.vertices, decimals)
    vertices2 *= flip(vertices2)
    order = np.lexsort(vertices2.T)
    vertices2 = vertices2[order]

    # The two hemispheres should have the same vertices up to their order
    nt.assert_array_equal(vertices1, vertices2)

    # Create a hemisphere from vertices
    hemi3 = HemiSphere(xyz=hemi1.vertices)
    nt.assert_array_equal(hemi1.faces, hemi3.faces)
    nt.assert_array_equal(hemi1.edges, hemi3.edges)

def test_bdg_initial_direction():
    """This test the number of inital direction."
    """

    hsph_updated = HemiSphere.from_sphere(unit_icosahedron).subdivide(2)
    vertices = hsph_updated.vertices
    bvecs = vertices
    bvals = np.ones(len(vertices)) * 1000
    bvecs = np.insert(bvecs, 0, np.array([0, 0, 0]), axis=0)
    bvals = np.insert(bvals, 0, 0)
    gtab = gradient_table(bvals, bvecs)

    # test that we get one direction when we have a single tensor
    voxel = single_tensor(gtab).reshape([1, 1, 1, -1])
    dti_model = dti.TensorModel(gtab)
    boot_dg = BootDirectionGetter.from_data(voxel, dti_model, 30, sh_order=6)
    initial_direction = boot_dg.initial_direction(np.zeros(3))
    npt.assert_equal(len(initial_direction), 1)
    npt.assert_allclose(initial_direction[0], [1, 0, 0], atol=0.1)

    # test that we get multiple directions when we have a multi-tensor
    mevals = np.array([[1.5, 0.4, 0.4], [1.5, 0.4, 0.4]]) * 1e-3
    fracs = [60, 40]
    voxel, primary_evecs = multi_tensor(gtab, mevals, fractions=fracs,
                                        snr=None)
    voxel = voxel.reshape([1, 1, 1, -1])
    response = (np.array([0.0015, 0.0004, 0.0004]), 1)
    csd_model = ConstrainedSphericalDeconvModel(gtab, response=response,
                                                sh_order=4)
    boot_dg = BootDirectionGetter.from_data(voxel, csd_model, 30)
    initial_direction = boot_dg.initial_direction(np.zeros(3))

    npt.assert_equal(len(initial_direction), 2)
    npt.assert_allclose(initial_direction, primary_evecs, atol=0.1)

def test_peak_direction_tracker():
    """This tests that the Peaks And Metrics Direction Getter plays nice
    LocalTracking and produces reasonable streamlines in a simple example.
    """
    sphere = HemiSphere.from_sphere(unit_octahedron)

    # A simple image with three possible configurations, a vertical tract,
    # a horizontal tract and a crossing
    peaks_values_lookup = np.array([[0., 0.],
                                    [1., 0.],
                                    [1., 0.],
                                    [0.5, 0.5]])
    peaks_indices_lookup = np.array([[-1, -1],
                                     [0, -1],
                                     [1, -1],
                                     [0,  1]])
    # PeaksAndMetricsDirectionGetter needs at 3 slices on each axis to work
    simple_image = np.zeros([5, 6, 3], dtype=int)
    simple_image[:, :, 1] = np.array([[0, 1, 0, 1, 0, 0],
                                      [0, 1, 0, 1, 0, 0],
                                      [0, 3, 2, 2, 2, 0],
                                      [0, 1, 0, 0, 0, 0],
                                      [0, 1, 0, 0, 0, 0],
                                      ])

    dg = PeaksAndMetrics()
    dg.sphere = sphere
    dg.peak_values = peaks_values_lookup[simple_image]
    dg.peak_indices = peaks_indices_lookup[simple_image]
    dg.ang_thr = 90

    mask = (simple_image >= 0).astype(float)
    tc = ThresholdTissueClassifier(mask, 0.5)
    seeds = [np.array([1., 1., 1.]),
             np.array([2., 4., 1.]),
             np.array([1., 3., 1.]),
             np.array([4., 4., 1.])]

    streamlines = LocalTracking(dg, tc, seeds, np.eye(4), 1.)

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

    for i, sl in enumerate(streamlines):
        npt.assert_(np.allclose(sl, expected[i]))

def test_sphere_scaling_csdmodel():
    """Check that mirroring regulization sphere does not change the result of
    csddeconv model"""
    _, fbvals, fbvecs = get_data('small_64D')

    bvals = np.load(fbvals)
    bvecs = np.load(fbvecs)

    gtab = gradient_table(bvals, bvecs)
    mevals = np.array(([0.0015, 0.0003, 0.0003],
                       [0.0015, 0.0003, 0.0003]))

    angles = [(0, 0), (60, 0)]

    S, sticks = multi_tensor(gtab, mevals, 100., angles=angles,
                             fractions=[50, 50], snr=None)

    sphere = get_sphere('symmetric362')
    hemi = HemiSphere.from_sphere(sphere)

    response = (np.array([0.0015, 0.0003, 0.0003]), 100)
    model_full = ConstrainedSphericalDeconvModel(gtab, response,
                                                reg_sphere=sphere)
    model_hemi = ConstrainedSphericalDeconvModel(gtab, response,
                                                reg_sphere=hemi)
    csd_fit_full = model_full.fit(S)
    csd_fit_hemi = model_hemi.fit(S)

    assert_array_almost_equal(csd_fit_full.shm_coeff, csd_fit_hemi.shm_coeff)

def test_MaximumDeterministicTracker():
    """This tests that the Maximum Deterministic Direction Getter plays nice
    LocalTracking and produces reasonable streamlines in a simple example.
    """
    sphere = HemiSphere.from_sphere(unit_octahedron)

    # A simple image with three possible configurations, a vertical tract,
    # a horizontal tract and a crossing
    pmf_lookup = np.array([[0., 0., 1.],
                           [1., 0., 0.],
                           [0., 1., 0.],
                           [.4, .6, 0.]])
    simple_image = np.array([[0, 1, 0, 0, 0, 0],
                             [0, 1, 0, 0, 0, 0],
                             [0, 3, 2, 2, 2, 0],
                             [0, 1, 0, 0, 0, 0],
                             [0, 1, 0, 0, 0, 0],
                             ])

    simple_image = simple_image[..., None]
    pmf = pmf_lookup[simple_image]

    seeds = [np.array([1., 1., 0.])] * 30

    mask = (simple_image > 0).astype(float)
    tc = ThresholdTissueClassifier(mask, .5)

    dg = DeterministicMaximumDirectionGetter.from_pmf(pmf, 90, sphere)
    streamlines = LocalTracking(dg, tc, seeds, np.eye(4), 1.)

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

    def allclose(x, y):
        return x.shape == y.shape and np.allclose(x, y)

    for sl in streamlines:
        dir = ( -sphere.vertices[0] ).copy()
        if not allclose(sl, expected[0]):
            raise AssertionError()

    # The first path is not possible if 90 degree turns are excluded
    dg = DeterministicMaximumDirectionGetter.from_pmf(pmf, 80, sphere)
    streamlines = LocalTracking(dg, tc, seeds, np.eye(4), 1.)

    for sl in streamlines:
        npt.assert_(np.allclose(sl, expected[1]))

def test_pmf_from_sh():
    sphere = HemiSphere.from_sphere(unit_octahedron)
    pmfgen = SHCoeffPmfGen(np.ones([2, 2, 2, 28]), sphere, None)

    # Test that the pmf is greater than 0 for a valid point
    pmf = pmfgen.get_pmf(np.array([0, 0, 0], dtype='float'))
    npt.assert_equal(np.sum(pmf) > 0, True)

    # Test that the pmf is 0 for invalid Points
    npt.assert_array_equal(pmfgen.get_pmf(np.array([-1, 0, 0], dtype='float')),
                           np.zeros(len(sphere.vertices)))
    npt.assert_array_equal(pmfgen.get_pmf(np.array([0, 0, 10], dtype='float')),
                           np.zeros(len(sphere.vertices)))

def test_closest_peak_tracker():
    """This tests that the Closest Peak Direction Getter plays nice
    LocalTracking and produces reasonable streamlines in a simple example.
    """
    sphere = HemiSphere.from_sphere(unit_octahedron)

    # A simple image with three possible configurations, a vertical tract,
    # a horizontal tract and a crossing
    pmf_lookup = np.array([[0., 0., 1.],
                           [1., 0., 0.],
                           [0., 1., 0.],
                           [.5, .5, 0.]])
    simple_image = np.array([[0, 1, 0, 0, 0, 0],
                             [0, 1, 0, 0, 0, 0],
                             [2, 3, 2, 2, 2, 0],
                             [0, 1, 0, 0, 0, 0],
                             [0, 1, 0, 0, 0, 0],
                             ])

    simple_image = simple_image[..., None]
    pmf = pmf_lookup[simple_image]

    seeds = [np.array([1., 1., 0.]), np.array([2., 4., 0.])]

    mask = (simple_image > 0).astype(float)
    tc = BinaryTissueClassifier(mask)

    dg = ClosestPeakDirectionGetter.from_pmf(pmf, 90, sphere,
                                             pmf_threshold=0.1)
    streamlines = Streamlines(LocalTracking(dg, tc, seeds, np.eye(4), 1.))

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

    def allclose(x, y):
        return x.shape == y.shape and np.allclose(x, y)

    if not allclose(streamlines[0], expected[0]):
        raise AssertionError()
    if not allclose(streamlines[1], expected[1]):
        raise AssertionError()

def test_pmf_from_array():
    sphere = HemiSphere.from_sphere(unit_octahedron)
    pmfgen = SimplePmfGen(np.ones([2, 2, 2, len(sphere.vertices)]))

    # Test that the pmf is greater than 0 for a valid point
    pmf = pmfgen.get_pmf(np.array([0, 0, 0], dtype='float'))
    npt.assert_equal(np.sum(pmf) > 0, True)

    # Test that the pmf is 0 for invalid Points
    npt.assert_array_equal(pmfgen.get_pmf(np.array([-1, 0, 0], dtype='float')),
                           np.zeros(len(sphere.vertices)))
    npt.assert_array_equal(pmfgen.get_pmf(np.array([0, 0, 10], dtype='float')),
                           np.zeros(len(sphere.vertices)))

    npt.assert_raises(
        ValueError,
        lambda: SimplePmfGen(np.ones([2, 2, 2, len(sphere.vertices)])*-1))

def test_bdg_residual():
    """This tests the bootstrapping residual.
    """

    hsph_updated = HemiSphere.from_sphere(unit_icosahedron).subdivide(2)
    vertices = hsph_updated.vertices
    bvecs = vertices
    bvals = np.ones(len(vertices)) * 1000
    bvecs = np.insert(bvecs, 0, np.array([0, 0, 0]), axis=0)
    bvals = np.insert(bvals, 0, 0)
    gtab = gradient_table(bvals, bvecs)
    r, theta, phi = cart2sphere(*vertices.T)
    B, m, n = shm.real_sym_sh_basis(6, theta, phi)
    shm_coeff = np.random.random(B.shape[1])

    # sphere_func is sampled of the spherical function for each point of
    # the sphere
    sphere_func = np.dot(shm_coeff, B.T)

    voxel = np.concatenate((np.zeros(1), sphere_func))
    data = np.tile(voxel, (3, 3, 3, 1))

    csd_model = ConstrainedSphericalDeconvModel(gtab, None, sh_order=6)
    boot_pmf_gen = BootPmfGen(data, model=csd_model, sphere=hsph_updated,
                              sh_order=6)

    # Two boot samples should be the same
    odf1 = boot_pmf_gen.get_pmf(np.array([1.5, 1.5, 1.5]))
    odf2 = boot_pmf_gen.get_pmf(np.array([1.5, 1.5, 1.5]))
    npt.assert_array_almost_equal(odf1, odf2)

    # A boot sample with less sh coeffs should have residuals, thus the two
    # should be different
    boot_pmf_gen2 = BootPmfGen(data, model=csd_model, sphere=hsph_updated,
                               sh_order=4)
    odf1 = boot_pmf_gen2.get_pmf(np.array([1.5, 1.5, 1.5]))
    odf2 = boot_pmf_gen2.get_pmf(np.array([1.5, 1.5, 1.5]))
    npt.assert_(np.any(odf1 != odf2))

    # test with a gtab with two shells and assert you get an error
    bvals[-1] = 2000
    gtab = gradient_table(bvals, bvecs)
    csd_model = ConstrainedSphericalDeconvModel(gtab, None, sh_order=6)
    npt.assert_raises(ValueError, BootPmfGen, data, csd_model, hsph_updated, 6)

def test_boot_pmf():
    """This tests the local model used for the bootstrapping.
    """
    hsph_updated = HemiSphere.from_sphere(unit_octahedron)
    vertices = hsph_updated.vertices
    bvecs = vertices
    bvals = np.ones(len(vertices)) * 1000
    bvecs = np.insert(bvecs, 0, np.array([0, 0, 0]), axis=0)
    bvals = np.insert(bvals, 0, 0)
    gtab = gradient_table(bvals, bvecs)
    voxel = single_tensor(gtab)
    data = np.tile(voxel, (3, 3, 3, 1))
    point = np.array([1., 1., 1.])
    tensor_model = TensorModel(gtab)

    boot_pmf_gen = BootPmfGen(data, model=tensor_model, sphere=hsph_updated)
    no_boot_pmf = boot_pmf_gen.get_pmf_no_boot(point)

    model_pmf = tensor_model.fit(voxel).odf(hsph_updated)

    npt.assert_equal(len(hsph_updated.vertices), no_boot_pmf.shape[0])
    npt.assert_array_almost_equal(no_boot_pmf, model_pmf)

    # test model sherical harminic order different than bootstrap order
    with warnings.catch_warnings(record=True) as w:
        warnings.simplefilter("always", category=UserWarning)
        csd_model = ConstrainedSphericalDeconvModel(gtab, None, sh_order=6)
        assert_greater(len([lw for lw in w if issubclass(lw.category,
                                                         UserWarning)]), 0)

    boot_pmf_gen_sh4 = BootPmfGen(data, model=csd_model, sphere=hsph_updated,
                                  sh_order=4)
    pmf_sh4 = boot_pmf_gen_sh4.get_pmf(point)
    npt.assert_equal(len(hsph_updated.vertices), pmf_sh4.shape[0])
    npt.assert_(np.sum(pmf_sh4.shape) > 0)

    boot_pmf_gen_sh8 = BootPmfGen(data, model=csd_model, sphere=hsph_updated,
                                  sh_order=8)
    pmf_sh8 = boot_pmf_gen_sh8.get_pmf(point)
    npt.assert_equal(len(hsph_updated.vertices), pmf_sh8.shape[0])
    npt.assert_(np.sum(pmf_sh8.shape) > 0)

def create_unit_hemisphere(recursion_level=2):
    """Creates a unit sphere by subdividing a unit octahedron, returns half
    the sphere.

    Parameters
    -------------
    recursion_level : int
        Level of subdivision, recursion_level=1 will return an octahedron,
        anything bigger will return a more subdivided sphere. The sphere will
        have $(4^recursion_level+2)/2$ vertices.

    Returns
    ---------
    HemiSphere :
        Half of a unit sphere.

    See Also
    ----------
    create_unit_sphere, Sphere, HemiSphere
    """
    sphere = create_unit_sphere(recursion_level)
    return HemiSphere.from_sphere(sphere)

def _create_mt_sim(mevals, angles, fractions, S0, SNR, half_sphere=False):

    _, fbvals, fbvecs = get_data('small_64D')

    bvals = np.load(fbvals)
    bvecs = np.load(fbvecs)

    gtab = gradient_table(bvals, bvecs)

    S, sticks = multi_tensor(gtab, mevals, S0, angles=angles,
                             fractions=fractions, snr=SNR)

    sphere = get_sphere('symmetric724').subdivide(2)

    if half_sphere:

        sphere = HemiSphere.from_sphere(sphere)

    odf_gt = multi_tensor_odf(sphere.vertices, mevals,
                              angles=angles, fractions=fractions)

    return odf_gt, sticks, sphere

def test_bdg_get_direction():
    """This tests the direction found by the bootstrap direction getter.
    """

    sphere = HemiSphere.from_sphere(unit_icosahedron.subdivide())
    two_neighbors = sphere.edges[0]
    direction1 = sphere.vertices[two_neighbors[0]]
    direction2 = sphere.vertices[two_neighbors[1]]
    angle = np.rad2deg(direction1.dot(direction2))
    point = np.zeros(3)
    prev_direction = direction2.copy()

    pmf = np.zeros([1, 1, 1, len(sphere.vertices)])
    pmf[:, :, :, two_neighbors[0]] = 1
    pmf_gen = SimplePmfGen(pmf)

    # test case in which no valid direction is found with default maxdir
    boot_dg = BootDirectionGetter(pmf_gen, angle / 2., sphere=sphere)
    npt.assert_equal(boot_dg.get_direction(point, prev_direction), 1)
    npt.assert_equal(direction2, prev_direction)

    # test case in which no valid direction is found with new max attempts
    boot_dg = BootDirectionGetter(pmf_gen, angle / 2., sphere=sphere,
                                  max_attempts=3)
    npt.assert_equal(boot_dg.get_direction(point, prev_direction), 1)
    npt.assert_equal(direction2, prev_direction)

    # test case in which a valid direction is found
    boot_dg = BootDirectionGetter(pmf_gen, angle * 2., sphere=sphere,
                                  max_attempts=1)
    npt.assert_equal(boot_dg.get_direction(point, prev_direction), 0)
    npt.assert_equal(direction1, prev_direction)

    # test invalid max_attempts parameters
    npt.assert_raises(
        ValueError,
        lambda: BootDirectionGetter(pmf_gen, angle * 2., sphere=sphere,
                                    max_attempts=0))

def test_ProbabilisticOdfWeightedTracker():
    sphere = HemiSphere.from_sphere(unit_octahedron)

    # A simple image with three possible configurations, a vertical tract,
    # a horizontal tract and a crossing
    odf_list = [np.array([0., 0., 0.]),
                np.array([1., 0., 0.]),
                np.array([0., 1., 0.]),
                np.array([1., 1., 0.]),
                ]
    simple_image = np.array([[0, 1, 0, 0, 0, 0],
                             [0, 1, 0, 0, 0, 0],
                             [0, 3, 2, 2, 2, 0],
                             [0, 1, 0, 0, 0, 0],
                             [0, 1, 0, 0, 0, 0],
                             ])
    # Make the image 4d
    simple_image = simple_image[..., None, None]

    # Simple model and fit for this image
    class MyModel():
        def fit(self, data):
            return MyFit(data)

    class MyFit(object):
        def __init__(self, n):
            self.n = n

        def odf(self, sphere):
            return odf_list[self.n]

    seeds = [np.array([1.5, 1.5, .5])] * 30
    model = MyModel()
    mask = np.ones([5, 6, 1], dtype="bool")
    stepper = FixedSizeStepper(1.)
    interpolator = NearestNeighborInterpolator(simple_image, (1, 1, 1))

    # These are the only two possible paths though the simple_image
    pwt = ProbabilisticOdfWeightedTracker(model, interpolator, mask,
                                          stepper, 90, seeds, sphere)
    expected = [np.array([[0.5, 1.5, 0.5],
                          [1.5, 1.5, 0.5],
                          [2.5, 1.5, 0.5],
                          [2.5, 2.5, 0.5],
                          [2.5, 3.5, 0.5],
                          [2.5, 4.5, 0.5],
                          [2.5, 5.5, 0.5]]),
                np.array([[0.5, 1.5, 0.5],
                          [1.5, 1.5, 0.5],
                          [2.5, 1.5, 0.5],
                          [3.5, 1.5, 0.5],
                          [4.5, 1.5, 0.5]])
                ]

    def allclose(x, y):
        return x.shape == y.shape and np.allclose(x, y)

    path = [False, False]
    for streamline in pwt:
        if allclose(streamline, expected[0]):
            path[0] = True
        elif allclose(streamline, expected[1]):
            path[1] = True
        else:
            raise AssertionError()
    assert_(all(path))

    # The first path is not possible if 90 degree turns are excluded
    pwt = ProbabilisticOdfWeightedTracker(model, interpolator, mask,
                                          stepper, 80, seeds, sphere)
    for streamline in pwt:
        assert_(np.allclose(streamline, expected[1]))