Python affines.apply_affine函数代码示例

def transform_streamlines(streamlines, mat, in_place=False):
    """ Apply affine transformation to streamlines

    Parameters
    ----------
    streamlines : Streamlines
        Streamlines object
    mat : array, (4, 4)
        transformation matrix
    in_place : bool
        If True then change data in place.
        Be careful changes input streamlines.

    Returns
    -------
    new_streamlines : Streamlines
        Sequence transformed 2D ndarrays of shape[-1]==3
    """
    # using new Streamlines API
    if isinstance(streamlines, Streamlines):
        if in_place:
            streamlines._data = apply_affine(mat, streamlines._data)
            return streamlines
        new_streamlines = streamlines.copy()
        new_streamlines._data = apply_affine(mat, new_streamlines._data)
        return new_streamlines
    # supporting old data structure of streamlines
    return [apply_affine(mat, s) for s in streamlines]

    def slice_volume(self, projection: str=SAGITTAL, ras: Union[numpy.ndarray, list]=ORIGIN)\
            -> (numpy.ndarray, numpy.ndarray, numpy.ndarray):
        """
        This determines slice colors and axes coordinates for the slice.
        :param projection: one of sagittal, axial or coronal
        :param ras: 3D point where to do the slicing
        :return: X, Y, 2D data matrix
        """

        affine_inverse = numpy.linalg.inv(self.affine_matrix)
        ijk_ras = numpy.round(apply_affine(affine_inverse, ras)).astype('i')

        slice_index_1, slice_index_2 = X_Y_INDEX[projection]

        slice_data = numpy.zeros(
            (self.dimensions[slice_index_1], self.dimensions[slice_index_2]))
        x_axis_coords = numpy.zeros_like(slice_data)
        y_axis_coords = numpy.zeros_like(slice_data)

        for i in range(self.dimensions[slice_index_1]):
            for j in range(self.dimensions[slice_index_2]):
                ijk_ras[slice_index_1] = i
                ijk_ras[slice_index_2] = j

                ras_coordinates = apply_affine(self.affine_matrix, ijk_ras)
                x_axis_coords[i, j] = ras_coordinates[slice_index_1]
                y_axis_coords[i, j] = ras_coordinates[slice_index_2]

                color = self.data[ijk_ras[0], ijk_ras[1], ijk_ras[2]]
                if isinstance(color, (list, numpy.ndarray)):
                    color = color[0]
                slice_data[i][j] = color

        return x_axis_coords, y_axis_coords, slice_data

    def _register_neighb_to_model(self, model_bundle, neighb_streamlines,
                                  metric=None, x0=None, bounds=None,
                                  select_model=400, select_target=600,
                                  method='L-BFGS-B',
                                  nb_pts=20, num_threads=None):

        if self.verbose:
            print('# Local SLR of neighb_streamlines to model')
            t = time()

        if metric is None or metric == 'symmetric':
            metric = BundleMinDistanceMetric(num_threads=num_threads)
        if metric == 'asymmetric':
            metric = BundleMinDistanceAsymmetricMetric()
        if metric == 'diagonal':
            metric = BundleSumDistanceMatrixMetric()

        if x0 is None:
            x0 = 'similarity'

        if bounds is None:
            bounds = [(-30, 30), (-30, 30), (-30, 30),
                      (-45, 45), (-45, 45), (-45, 45), (0.8, 1.2)]

        # TODO this can be speeded up by using directly the centroids
        static = select_random_set_of_streamlines(model_bundle,
                                                  select_model, rng=self.rng)
        moving = select_random_set_of_streamlines(neighb_streamlines,
                                                  select_target, rng=self.rng)

        static = set_number_of_points(static, nb_pts)
        moving = set_number_of_points(moving, nb_pts)

        slr = StreamlineLinearRegistration(metric=metric, x0=x0,
                                           bounds=bounds,
                                           method=method)
        slm = slr.optimize(static, moving)

        transf_streamlines = neighb_streamlines.copy()
        transf_streamlines._data = apply_affine(
            slm.matrix, transf_streamlines._data)

        transf_matrix = slm.matrix
        slr_bmd = slm.fopt
        slr_iterations = slm.iterations

        if self.verbose:
            print(' Square-root of BMD is %.3f' % (np.sqrt(slr_bmd),))
            if slr_iterations is not None:
                print(' Number of iterations %d' % (slr_iterations,))
            print(' Matrix size {}'.format(slm.matrix.shape))
            original = np.get_printoptions()
            np.set_printoptions(3, suppress=True)
            print(transf_matrix)
            print(slm.xopt)
            np.set_printoptions(**original)

            print(' Duration %0.3f sec. \n' % (time() - t,))

        return transf_streamlines, slr_bmd

def _get_cut_slices(stat_map_img, cut_coords=None, threshold=None):
    """ For internal use.
        Find slice numbers for the cut.
        Based on find_xyz_cut_coords
    """
    # Select coordinates for the cut
    if cut_coords is None:
        cut_coords = find_xyz_cut_coords(
            stat_map_img, activation_threshold=threshold)

    # Convert cut coordinates into cut slices
    try:
        cut_slices = apply_affine(np.linalg.inv(stat_map_img.affine),
                                  cut_coords)
    except ValueError:
        raise ValueError(
            "The input given for display_mode='ortho' needs to be "
            "a list of 3d world coordinates in (x, y, z). "
            "You provided cut_coords={0}".format(cut_coords))
    except IndexError:
        raise ValueError(
            "The input given for display_mode='ortho' needs to be "
            "a list of 3d world coordinates in (x, y, z). "
            "You provided single cut, cut_coords={0}".format(cut_coords))

    return cut_slices

        def display_extent(self, x1, x2, y1, y2, z1, z2):

            tmp_mask = np.zeros(grid_shape, dtype=np.bool)
            tmp_mask[x1:x2 + 1, y1:y2 + 1, z1:z2 + 1] = True
            tmp_mask = np.bitwise_and(tmp_mask, mask)

            ijk = np.ascontiguousarray(np.array(np.nonzero(tmp_mask)).T)
            if len(ijk) == 0:
                self.SetMapper(None)
                return
            if affine is not None:
                ijk_trans = np.ascontiguousarray(apply_affine(affine, ijk))
            list_dirs = []
            for index, center in enumerate(ijk):
                # center = tuple(center)
                if affine is None:
                    xyz = center[:, None]
                else:
                    xyz = ijk_trans[index][:, None]
                xyz = xyz.T
                for i in range(peaks_dirs[tuple(center)].shape[-2]):

                    if peaks_values is not None:
                        pv = peaks_values[tuple(center)][i]
                    else:
                        pv = 1.
                    symm = np.vstack((-peaks_dirs[tuple(center)][i] * pv + xyz,
                                      peaks_dirs[tuple(center)][i] * pv + xyz))
                    list_dirs.append(symm)

            self.mapper = line(list_dirs, colors=colors,
                               opacity=opacity, linewidth=linewidth,
                               lod=lod, lod_points=lod_points,
                               lod_points_size=lod_points_size).GetMapper()
            self.SetMapper(self.mapper)

def get_MRI_values(path2mri_file,function_space,mesh, offset=None):
		"""
		==============	 ==============================================================
		Argument                Explanation
		==============	 ==============================================================
		path2_mri_file	The path to the MRI file with extension mgz. typically orig.mgz
		function_space	The function space of the mesh.		
		mesh 	 	The mesh of the brain	 
		function	Instead of a return value, update the function inside. 	
		offset 		The translation of the origin of the mesh , need to be np.array
		"""
		import nibabel
		from nibabel.affines import apply_affine
		import numpy.linalg as npl
		import numpy as np

		img= nibabel.load(path2mri_file) 
		inv_aff = npl.inv ( img.get_header().get_vox2ras_tkr() )
		data = img.get_data()
		if offset==None :
			offset=np.array([0,0,0])

		xyz = function_space.dofmap().tabulate_all_coordinates(mesh).reshape((function_space.dim(),-1)) - offset

		i,j,k = apply_affine(inv_aff,xyz).T
	
		# 4.7 = 5 and not 4
		i= map(round,i) 
		j= map(round,j) 
		k= map(round,k) 
	
		return np.array(map(data.item,i,j,k),dtype=float)

def get_data():

    #dname = '/home/eleftherios/Data/Test_data_Jasmeen/Elef_Test_RecoBundles/'
    dname = '/home/eleftherios/Data/Elef_Test_RecoBundles/'
    fname = dname + 'tracts.trk'
    fname_npz = dname + 'tracts.npz'

    streamlines = nib.streamlines.compact_list.load_compact_list(fname_npz)

    streamlines = streamlines[::10].copy()
    streamlines._data -= np.mean(streamlines._data, axis=0)

    # Rotate brain to see a sagital view.
    R1 = np.eye(4)
    R1[:3, :3] = rodrigues_axis_rotation((0, 1, 0), theta=90)
    R2 = np.eye(4)
    R2[:3, :3] = rodrigues_axis_rotation((0, 0, 1), theta=90)
    R = np.dot(R2, R1)
    streamlines._data = apply_affine(R, streamlines._data)

#    renderer = window.Renderer()
#    bundle_actor = actor.line(streamlines)
#    renderer.add(bundle_actor)
#    window.show(renderer)

    return streamlines

 def main(self,X,affine):
   inv_affine = numpy.linalg.inv(affine)
   shift05 = numpy.array([0.5,0.5,0.5])
   points = []
   for i,line in enumerate(X):
     points.append((apply_affine(inv_affine,line)+shift05).astype(int).tolist())
   return FancyDict(points=points)

 def set_space_pos(self, new_space_coord):
     """Set current cursor position based on RAS coordinates."""
     new_coord = apply_affine(npl.inv(self._affine), new_space_coord)
     new_coord = np.floor(new_coord)
     new_coord = [int(item) for item in new_coord]
     if self.is_valid_coord(new_coord):
         self.set_cross_pos(new_coord)

def extract_cube(vol_target,vol_tdata,affine_ref,initial_position,epi_coord,cube_size):
        # Calculate transformation from ref to target vol
        # calculate the affine transf matrix from EPI to anat
        epi_vox2anat_vox = npl.inv(vol_target.get_affine()).dot(affine_ref)
        # x,y,z
        target_anat = apply_affine(epi_vox2anat_vox,epi_coord)
        # get the cube
        side_size = int(np.floor(cube_size/2.))

        # target spacing
        side_size_targ = (apply_affine(epi_vox2anat_vox,np.array(epi_coord)+side_size) - target_anat)
        side_size_targ = side_size_targ[0]

        x_s,y_s,z_s = vol_target.get_data().shape
        x_interval = np.arange(target_anat[0]-side_size_targ,target_anat[0]+side_size_targ+1,dtype=int)
        y_interval = np.arange(target_anat[1]-side_size_targ,target_anat[1]+side_size_targ+1,dtype=int)
        z_interval = np.arange(target_anat[2]-side_size_targ,target_anat[2]+side_size_targ+1,dtype=int)

        # normalize the data between 0 and 1
        #norm_vol_target = vol_target.get_data()
        norm_vol_target = vol_tdata#avg_frame_norm(vol_target.get_data())
        #norm_vol_target = rescale(norm_vol_target)
        # check if we need padding
        if (x_interval>=0).all() & (x_interval<(x_s)).all() & (y_interval>=0).all() & (y_interval<(y_s)).all() & (z_interval>=0).all() & (z_interval<(z_s)).all():
            small_cube = norm_vol_target[x_interval,...][:,y_interval,...][:,:,z_interval,...]
        else:
            padded_target = np.lib.pad(norm_vol_target, (side_size,side_size), 'constant', constant_values=(0))
            x_interval = np.arange(target_anat[0]-side_size_targ,target_anat[0]+side_size_targ+1,dtype=int)
            y_interval = np.arange(target_anat[1]-side_size_targ,target_anat[1]+side_size_targ+1,dtype=int)
            z_interval = np.arange(target_anat[2]-side_size_targ,target_anat[2]+side_size_targ+1,dtype=int)
            small_cube = padded_target[x_interval+side_size,...][:,y_interval+side_size,...][:,:,z_interval+side_size,...]

        # pad the vols
        #padded_target = np.lib.pad(vol_target.get_data(), (side_size,side_size), 'constant', constant_values=(0))
        #small_cube = padded_target[x_interval+side_size,...][:,y_interval+side_size,...][:,:,z_interval+side_size,...]

        #x_interval = np.arange(target_anat[0]-side_size,target_anat[0]+side_size,dtype=int)
        #y_interval = np.arange(target_anat[1]-side_size,target_anat[1]+side_size,dtype=int)
        #z_interval = np.arange(target_anat[2]-side_size,target_anat[2]+side_size,dtype=int)

        #padded_pos_vol = np.lib.pad(initial_position, (side_size,side_size), 'constant', constant_values=(0))[...,side_size:-side_size]
        #init_pos = padded_pos_vol[x_interval+side_size,...][:,y_interval+side_size,...][:,:,z_interval+side_size,...]
        #init_pos = initial_position[target_anat[0],target_anat[1],target_anat[2],...]
        init_pos = initial_position[epi_coord[0],epi_coord[1],epi_coord[2],...]
        return small_cube, init_pos

def test_creation():
    # This is the simplest possible example, where there is a thing we are
    # optimizing, and an optional pre and post transform
    # Reset the aff2 object
    aff2_obj = Affine(AFF2.copy())
    ct = ChainTransform(aff2_obj)
    # Check apply gives expected result
    assert_array_equal(ct.apply(POINTS),
                       apply_affine(AFF2, POINTS))
    # Check that result is changed by setting params
    assert_array_equal(ct.param, aff2_obj.param)
    ct.param = np.zeros((12,))
    assert_array_almost_equal(ct.apply(POINTS), POINTS)
    # Does changing params in chain object change components passed in?
    assert_array_almost_equal(aff2_obj.param, np.zeros((12,)))
    # Reset the aff2 object
    aff2_obj = Affine(AFF2.copy())
    # Check apply gives the expected results
    ct = ChainTransform(aff2_obj, pre=AFF1_OBJ)
    assert_array_almost_equal(AFF1_OBJ.as_affine(), AFF1)
    assert_array_almost_equal(aff2_obj.as_affine(), AFF2)
    tmp = np.dot(AFF2, AFF1)
    assert_array_almost_equal(ct.apply(POINTS),
                       apply_affine(np.dot(AFF2, AFF1), POINTS))
    # Check that result is changed by setting params
    assert_array_almost_equal(ct.param, aff2_obj.param)
    ct.param = np.zeros((12,))
    assert_array_almost_equal(ct.apply(POINTS), apply_affine(AFF1, POINTS))
    # Does changing params in chain object change components passed in?
    assert_array_almost_equal(aff2_obj.param, np.zeros((12,)))
    # Reset the aff2 object
    aff2_obj = Affine(AFF2.copy())
    ct = ChainTransform(aff2_obj, pre=AFF1_OBJ, post=AFF3_OBJ)
    assert_array_almost_equal(ct.apply(POINTS),
                       apply_affine(np.dot(AFF3, np.dot(AFF2, AFF1)), POINTS))
    # Check that result is changed by setting params
    assert_array_equal(ct.param, aff2_obj.param)
    ct.param = np.zeros((12,))
    assert_array_almost_equal(ct.apply(POINTS),
                              apply_affine(np.dot(AFF3, AFF1), POINTS))
    # Does changing params in chain object change components passed in?
    assert_array_equal(aff2_obj.param, np.zeros((12,)))

def transform_streamlines(streamlines, mat):
    """ Apply affine transformation to streamlines

    Parameters
    ----------
    streamlines : list
        List of 2D ndarrays of shape[-1]==3

    Returns
    -------
    new_streamlines : list
        List of the transformed 2D ndarrays of shape[-1]==3
    """
    return [apply_affine(mat, s) for s in streamlines]

def cuts_from_bbox(mask_nii, cuts=3):
    """Finds equi-spaced cuts for presenting images"""
    from nibabel.affines import apply_affine

    mask_data = mask_nii.get_data() > 0.0

    # First, project the number of masked voxels on each axes
    ijk_counts = [
        mask_data.sum(2).sum(1),  # project sagittal planes to transverse (i) axis
        mask_data.sum(2).sum(0),  # project coronal planes to to longitudinal (j) axis
        mask_data.sum(1).sum(0),  # project axial planes to vertical (k) axis
    ]

    # If all voxels are masked in a slice (say that happens at k=10),
    # then the value for ijk_counts for the projection to k (ie. ijk_counts[2])
    # at that element of the orthogonal axes (ijk_counts[2][10]) is
    # the total number of voxels in that slice (ie. Ni x Nj).
    # Here we define some thresholds to consider the plane as "masked"
    # The thresholds vary because of the shape of the brain
    # I have manually found that for the axial view requiring 30%
    # of the slice elements to be masked drops almost empty boxes
    # in the mosaic of axial planes (and also addresses #281)
    ijk_th = [
        int((mask_data.shape[1] * mask_data.shape[2]) * 0.2),   # sagittal
        int((mask_data.shape[0] * mask_data.shape[2]) * 0.0),   # coronal
        int((mask_data.shape[0] * mask_data.shape[1]) * 0.3),   # axial
    ]

    vox_coords = []
    for ax, (c, th) in enumerate(zip(ijk_counts, ijk_th)):
        B = np.argwhere(c > th)
        if B.size:
            smin, smax = B.min(), B.max()

        # Avoid too narrow selections of cuts (very small masks)
        if not B.size or (th > 0 and (smin + cuts + 1) >= smax):
            B = np.argwhere(c > 0)

        # Resort to full plane if mask is seemingly empty
        smin, smax = B.min(), B.max() if B.size else (0, mask_data.shape[ax])
        inc = (smax - smin) / (cuts + 1)
        vox_coords.append([smin + (i + 1) * inc for i in range(cuts)])

    ras_coords = []
    for cross in np.array(vox_coords).T:
        ras_coords.append(apply_affine(
            mask_nii.affine, cross).tolist())
    ras_cuts = [list(coords) for coords in np.transpose(ras_coords)]
    return {k: v for k, v in zip(['x', 'y', 'z'], ras_cuts)}

def get_motion_deriv(aff_transforms,tmp_deriv_init):
    world_motion=[]
    point_motion_deriv = []
    tmp_deriv_old = np.copy(tmp_deriv_init)
    for ii in range(aff_transforms.shape[2]):
        world_motion.append(apply_affine(aff_transforms[...,ii],tmp_deriv_old))
    
    # compute the delta
    world_motion = np.array(world_motion)
    deriv_values = world_motion[1:]-world_motion[0:-1]
    
    if len(deriv_values.shape)>3:
        deriv_values = np.squeeze(np.swapaxes(deriv_values[...,np.newaxis],0,5))
    else:
        deriv_values = np.squeeze(np.swapaxes(deriv_values[...,np.newaxis],0,2))
    
    return deriv_values

    def apply_affine(self, affine, lazy=False):
        """ Applies an affine transformation on the points of each streamline.

        If `lazy` is not specified, this is performed *in-place*.

        Parameters
        ----------
        affine : ndarray of shape (4, 4)
            Transformation that will be applied to every streamline.
        lazy : {False, True}, optional
            If True, streamlines are *not* transformed in-place and a
            :class:`LazyTractogram` object is returned. Otherwise, streamlines
            are modified in-place.

        Returns
        -------
        tractogram : :class:`Tractogram` or :class:`LazyTractogram` object
            Tractogram where the streamlines have been transformed according
            to the given affine transformation. If the `lazy` option is true,
            it returns a :class:`LazyTractogram` object, otherwise it returns a
            reference to this :class:`Tractogram` object with updated
            streamlines.
        """
        if lazy:
            lazy_tractogram = LazyTractogram.from_tractogram(self)
            return lazy_tractogram.apply_affine(affine)

        if len(self.streamlines) == 0:
            return self

        if np.all(affine == np.eye(4)):
            return self  # No transformation.

        BUFFER_SIZE = 10000000  # About 128 Mb since pts shape is 3.
        for start in range(0, len(self.streamlines.data), BUFFER_SIZE):
            end = start + BUFFER_SIZE
            pts = self.streamlines._data[start:end]
            self.streamlines.data[start:end] = apply_affine(affine, pts)

        if self.affine_to_rasmm is not None:
            # Update the affine that brings back the streamlines to RASmm.
            self.affine_to_rasmm = np.dot(self.affine_to_rasmm,
                                          np.linalg.inv(affine))

        return self