""" Computes the mam distance between two streamlines. """
# For simplicity, features will be the vector between endpoints of a streamline. super(mam, self).__init__(feature=VectorOfEndpointsFeature()) def are_compatible(self, shape1, shape2): """ Checks if two features are vectors of same dimension. Basically this method exists so we don't have to do this check inside the `dist` method (speedup). """ return shape1 == shape2 and shape1 == 1 def dist(self, v1, v2): track1 = np.ascontiguousarray(v1, dtype=np.float32) t1_len = track1.shape track2 = np.ascontiguousarray(v2, dtype=np.float32) t2_len = track2.shape # preallocate buffer array for track distance calculations #distances_buffer = np.zeros((t1_len ,), dtype=np.float32) min_t2t1 = np.zeros((t2_len ,), dtype=np.float32) min_t1t2 = np.zeros((t1_len ,), dtype=np.float32) for t2_pi in range(0,t2_len): min_t2t1[t2_pi] = np.inf for t1_pi in range(0,t1_len): min_t1t2[t1_pi] = np.inf # pointer to current point in track 1 t1_pt = track1 t2_pt = track2 # calculate min squared distance between each point in the two # lines. Squared distance to delay doing the sqrt until after this # speed-critical loop for t1_pi in range(0,t1_len): # pointer to current point in track 2 for t2_pi in range(0,t2_len): d0 = t1_pt[t1_pi] - t2_pt[t2_pi] d1 = t1_pt[t1_pi] - t2_pt[t2_pi] delta2 = d0*d0 + d1*d1 #+ d2*d2 if delta2 < min_t1t2[t1_pi]: min_t1t2[t1_pi]=delta2 for t2_pi in range(0,t2_len): # pointer to current point in track 2 for t1_pi in range(0,t1_len): d0 = t1_pt[t1_pi] - t2_pt[t2_pi] d1 = t1_pt[t1_pi] - t2_pt[t2_pi] delta2 = d0*d0 + d1*d1 #+ d2*d2 if delta2 < min_t2t1[t2_pi]: min_t2t1[t2_pi]=delta2 # sqrt to get Euclidean distance from squared distance for t1_pi in range(0,t1_len): min_t1t2[t1_pi]=math.sqrt(min_t1t2[t1_pi]) for t2_pi in range(0,t2_len): min_t2t1[t2_pi]=math.sqrt(min_t2t1[t2_pi]) mean_t2t1 = 0 mean_t1t2 = 0 for t1_pi in range(0, t1_len): mean_t1t2+=min_t1t2[t1_pi] mean_t1t2=mean_t1t2 / t1_len for t2_pi in range(0, t2_len): mean_t2t1+=min_t2t1[t2_pi] mean_t2t1=mean_t2t1 / t2_len return np.min((mean_t2t1,mean_t1t2))
metric = mam()
qb2 = QuickBundles(threshold=0.15, metric=metric)
clus = qb2.cluster(streamlines)
load_tractogramusing the b0 image as a reference and plotting them with spherical ROIs using
stream_actor, the two are in very different spaces. It looks like the origin is shifted and the spatial scaling is different. Is there some transformation happening under the hood in
load_tractogramsomewhere? If so, I'm guessing applying the same transformation to my coordinates will solve the issue.
VectorOfEndpointsFeature()and then I will remove
shape1 == 1on your function
are_compatible. You do not need vector in your case for MAM. Let me know if it works, otherwise create an issue on DIPY, it will be easier to help you
to_xxxxmethods after loading, but without success. They shift the tracks, but not into the space my ROI coordinates are in.
I would be grateful for your thoughts on what could be wrong when feeding labels in utils.connectivity.matrix: M, grouping = utils.connectivitymatrix(streamlines, labels.astype(int), affine_labels,
File "/Users/alex/py3/lib/python3.7/site-packages/dipy/tracking/utils.py", line 157, in connectivity_matrix
raise ValueError("label_volume must be a 3d integer array with labelvolume must be a 3d integer array withnon-negative label values
labels.shape= (200, 400, 200) when labels.max() = 332 and labels.min() = 0. Thanks much!
labelsin your function arguments. To make sure that is correct, be explicit on your function:
utils.connectivitymatrix(streamlines, label_volume=labels.astype(int), affine=affine_labels, return_mapping=False, mapping_as_streamlines=False)
Does anybody know where the source code for this denoise library can be found: https://dipy.org/documentation/1.0.0./reference/dipy.denoise/. Particularly the nl means filters
The source code is found here: https://github.com/nipy/dipy/blob/master/dipy/denoise/
I have a look, but I want to display my TRK file and t1 data, I don't know what to change?Can you help me @arokem
DIPY examples/tutorials use data that is fetched directly from a few sites that host the data of interest. If you want to adapt the examples to use your own data, then, instead of using the
read_ * commands, you will need to read your local files.
In order to read your t1, you will need to use NiBabel as in
import nibabel as nib t1_img = nib.load(your_t1_filename) t1 = t1_img.get_data() fig = plt.figure() plt.imshow(t1, cmap="gray")
Similarly, to show your TRK file, you can do something similar to this
from dipy.io.streamline import load_tractogram from fury import window, actor tractogram = load_tractogram(your_trk_filename, your_t1_filename) streamlines = tractogram.streamlines renderer = window.Renderer() stream_actor = actor.line(streamlines) renderer.add(stream_actor) window.record(renderer)
If you want to display your streamlines and your T1 at the same time, you will need to add both actors to the same scene, adding to the above something like
t1_actor = actor.slicer(t1) renderer.add(t1_actor)