Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Another way to add orientation information #80

Open
wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Another way to add orientation information #80

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions bin/wm_cluster_atlas.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -509,8 +509,8 @@
cluster_distances = wma.cluster._pairwise_distance_matrix(pd_c, 0.0, number_of_jobs=1, bilateral=bilateral, distance_method=distance_method, sigmasq = cluster_local_sigma * cluster_local_sigma)
cluster_similarity = cluster_distances
else:
cluster_distances = wma.cluster._pairwise_distance_matrix(pd_c, 0.0, number_of_jobs=1, bilateral=bilateral, distance_method=distance_method)
cluster_similarity = wma.similarity.distance_to_similarity(cluster_distances, cluster_local_sigma * cluster_local_sigma)
cluster_distances, cluster_orientations = wma.cluster._pairwise_distance_matrix(pd_c, 0.0, number_of_jobs=1, bilateral=bilateral, distance_method=distance_method)
cluster_similarity = wma.similarity.distance_to_similarity(cluster_distances, cluster_orientations, cluster_local_sigma * cluster_local_sigma, sigmasq2 = 10)

#p(f1) = sum over all f2 of p(f1|f2) * p(f2)
# by using sample we estimate expected value of the above
Expand Down
59 changes: 39 additions & 20 deletions whitematteranalysis/cluster.py
100644 → 100755
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -665,21 +665,31 @@ def _rectangular_distance_matrix(input_polydata_n, input_polydata_m, threshold,
landmarks_n = numpy.zeros((fiber_array_n.number_of_fibers,3))

# pairwise distance matrix
all_fibers_n = range(0, fiber_array_n.number_of_fibers)

distances = Parallel(n_jobs=number_of_jobs,
verbose=0)(
delayed(similarity.fiber_distance)(
# all_fibers_n = range(0, fiber_array_n.number_of_fibers)
#
# distances = Parallel(n_jobs=number_of_jobs,
# verbose=0)(
# delayed(similarity.fiber_distance)(
# fiber_array_n.get_fiber(lidx),
# fiber_array_m,
# threshold, distance_method=distance_method,
# fiber_landmarks=landmarks_n[lidx,:],
# landmarks=landmarks_m, bilateral=bilateral)
# for lidx in all_fibers_n)
distances = numpy.zeros([fiber_array_n.number_of_fibers,fiber_array_m.number_of_fibers])
orientations = numpy.zeros([fiber_array_n.number_of_fibers,fiber_array_m.number_of_fibers])
for lidx in xrange(0,fiber_array_n.number_of_fibers):
distances[lidx,:], orientations[lidx,:] = similarity.fiber_distance(
fiber_array_n.get_fiber(lidx),
fiber_array_m,
threshold, distance_method=distance_method,
fiber_landmarks=landmarks_n[lidx,:],
landmarks=landmarks_m, bilateral=bilateral)
for lidx in all_fibers_n)

distances = numpy.array(distances).T
orientations = numpy.array(orientations).T

return distances
return distances, orientations

def _rectangular_similarity_matrix(input_polydata_n, input_polydata_m, threshold, sigma,
number_of_jobs=3, landmarks_n=None, landmarks_m=None, distance_method='Hausdorff',
Expand All @@ -696,15 +706,15 @@ def _rectangular_similarity_matrix(input_polydata_n, input_polydata_m, threshold

"""

distances = _rectangular_distance_matrix(input_polydata_n, input_polydata_m, threshold,
distances, orientations = _rectangular_distance_matrix(input_polydata_n, input_polydata_m, threshold,
number_of_jobs, landmarks_n, landmarks_m, distance_method, bilateral=bilateral)

if distance_method == 'StrictSimilarity':
similarity_matrix = distances
else:
# similarity matrix
sigmasq = sigma * sigma
similarity_matrix = similarity.distance_to_similarity(distances, sigmasq)
similarity_matrix = similarity.distance_to_similarity(distances, orientations, sigmasq, sigmasq2 = 10)

return similarity_matrix

Expand All @@ -731,28 +741,37 @@ def _pairwise_distance_matrix(input_polydata, threshold,
fiber_array.convert_from_polydata(input_polydata, points_per_fiber=15)

# pairwise distance matrix
all_fibers = range(0, fiber_array.number_of_fibers)
# all_fibers = range(0, fiber_array.number_of_fibers)

if landmarks is None:
landmarks2 = numpy.zeros((fiber_array.number_of_fibers,3))
else:
landmarks2 = landmarks

distances = Parallel(n_jobs=number_of_jobs,
verbose=0)(
delayed(similarity.fiber_distance)(
# distances = Parallel(n_jobs=number_of_jobs,
# verbose=0)(
# delayed(similarity.fiber_distance)(
# fiber_array.get_fiber(lidx),
# fiber_array,
# threshold, distance_method=distance_method,
# fiber_landmarks=landmarks2[lidx,:],
# landmarks=landmarks, bilateral=bilateral, sigmasq=sigmasq)
# for lidx in all_fibers)
#
# distances = numpy.array(distances)
distances = numpy.zeros([fiber_array.number_of_fibers,fiber_array.number_of_fibers])
orientations = numpy.zeros([fiber_array.number_of_fibers,fiber_array.number_of_fibers])
for lidx in xrange(0,fiber_array.number_of_fibers):
distances[lidx,:], orientations[lidx,:] = similarity.fiber_distance(
fiber_array.get_fiber(lidx),
fiber_array,
threshold, distance_method=distance_method,
fiber_landmarks=landmarks2[lidx,:],
landmarks=landmarks, bilateral=bilateral, sigmasq=sigmasq)
for lidx in all_fibers)

distances = numpy.array(distances)

# remove outliers if desired????

return distances
return distances, orientations

def _pairwise_similarity_matrix(input_polydata, threshold, sigma,
number_of_jobs=3, landmarks=None, distance_method='Hausdorff',
Expand All @@ -769,23 +788,23 @@ def _pairwise_similarity_matrix(input_polydata, threshold, sigma,

"""

distances = _pairwise_distance_matrix(input_polydata, threshold,
distances, orientations = _pairwise_distance_matrix(input_polydata, threshold,
number_of_jobs, landmarks, distance_method, bilateral=bilateral)

if distance_method == 'StrictSimilarity':
similarity_matrix = distances
else:
# similarity matrix
sigmasq = sigma * sigma
similarity_matrix = similarity.distance_to_similarity(distances, sigmasq)
similarity_matrix = similarity.distance_to_similarity(distances, orientations, sigmasq, sigmasq2 = 10)

# sanity check that on-diagonal elements are all 1
#print "This should be 1.0: ", numpy.min(numpy.diag(similarity_matrix))
#print numpy.min(numpy.diag(similarity_matrix)) == 1.0
# test
if __debug__:
# this tests that on-diagonal elements are all 1
test = numpy.min(numpy.diag(similarity_matrix)) == 1.0
test = numpy.min(numpy.diag(similarity_matrix)) >= 0.9
if not test:
print "<cluster.py> ERROR: On-diagonal elements are not all 1.0."
print" Minimum on-diagonal value:", numpy.min(numpy.diag(similarity_matrix))
Expand Down
Empty file modified whitematteranalysis/fibers.pyx
100644 → 100755
View file
Open in desktop
Empty file.
48 changes: 41 additions & 7 deletions whitematteranalysis/similarity.pyx
100644 → 100755
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,10 +6,10 @@ import sys

sys.setrecursionlimit(1000000)

def distance_to_similarity(distance, sigmasq=100):
def distance_to_similarity(distance, orientation, sigmasq1=100, sigmasq2=100):

# compute the similarities using Gaussian kernel
similarities = numpy.exp(-distance / (sigmasq))
similarities = numpy.exp(-(distance / (sigmasq1) + numpy.square(orientation)/ (sigmasq2)))

return similarities

Expand All @@ -27,31 +27,41 @@ def fiber_distance(fiber, fiber_array, threshold=0, distance_method='MeanSquared
fiber_equiv = fiber.get_equivalent_fiber()

# compute pairwise fiber distances along fibers
distance_1 = _fiber_distance_internal_use(fiber, fiber_array, threshold, distance_method, fiber_landmarks, landmarks, sigmasq)
distance_2 = _fiber_distance_internal_use(fiber_equiv, fiber_array, threshold, distance_method, fiber_landmarks, landmarks, sigmasq)
distance_1, orientation_1 = _fiber_distance_internal_use(fiber, fiber_array, threshold, distance_method, fiber_landmarks, landmarks, sigmasq)
distance_2, orientation_2 = _fiber_distance_internal_use(fiber_equiv, fiber_array, threshold, distance_method, fiber_landmarks, landmarks, sigmasq)

# choose the lowest distance, corresponding to the optimal fiber
# representation (either forward or reverse order)
orientation = numpy.zeros(orientation_1.shape)
if distance_method == 'StrictSimilarity':
# for use in laterality
# this is the product of all similarity values along the fiber
distance = numpy.maximum(distance_1, distance_2)
orientation[(distance_1 - distance_2) >= 0] = orientation_1[(distance_1 - distance_2) >= 0]
orientation[(distance_1 - distance_2) < 0] = orientation_2[(distance_1 - distance_2) < 0]
else:
distance = numpy.minimum(distance_1, distance_2)
orientation[(distance_1 - distance_2) >= 0] = orientation_2[(distance_1 - distance_2) >= 0]
orientation[(distance_1 - distance_2) < 0] = orientation_1[(distance_1 - distance_2) < 0]

if bilateral:
fiber_reflect = fiber.get_reflected_fiber()
# call this function again with the reflected fiber. Do NOT reflect again (bilateral=False) to avoid infinite loop.
distance_reflect = fiber_distance(fiber_reflect, fiber_array, threshold, distance_method, fiber_landmarks, landmarks, sigmasq, bilateral=False)
distance_reflect, orientation_reflect = fiber_distance(fiber_reflect, fiber_array, threshold, distance_method, fiber_landmarks, landmarks, sigmasq, bilateral=False)
# choose the best distance, corresponding to the optimal fiber
# representation (either reflected or not)
if distance_method == 'StrictSimilarity':
# this is the product of all similarity values along the fiber

orientation[(distance - distance_reflect) >= 0] = orientation[(distance - distance_reflect) >= 0]
orientation[(distance - distance_reflect) < 0] = orientation_reflect[(distance - distance_reflect) < 0]
distance = numpy.maximum(distance, distance_reflect)
else:
orientation[(distance - distance_reflect) >= 0] = orientation_reflect[(distance - distance_reflect) >= 0]
orientation[(distance - distance_reflect) < 0] = orientation[(distance - distance_reflect) < 0]
distance = numpy.minimum(distance, distance_reflect)

return distance
return distance, orientation

def fiber_distance_oriented(fiber, fiber_array, threshold=0, distance_method='MeanSquared', fiber_landmarks=None, landmarks=None, sigmasq=6400, bilateral=False):
"""Find pairwise fiber distance from fiber to all fibers in fiber_array, without testing both fiber orders.
Expand Down Expand Up @@ -106,6 +116,30 @@ def _fiber_distance_internal_use(fiber, fiber_array, threshold=0, distance_metho
# sum dx dx dz at each point on the fiber and sqrt for threshold
#distance = numpy.sqrt(dx + dy + dz)
distance = dx + dy + dz

array_orien_r = numpy.delete(fiber_array.fiber_array_r, -1, axis = 1) - numpy.delete(fiber_array.fiber_array_r, 0, axis = 1)
array_orien_a = numpy.delete(fiber_array.fiber_array_a, -1, axis = 1) - numpy.delete(fiber_array.fiber_array_a, 0, axis = 1)
array_orien_s = numpy.delete(fiber_array.fiber_array_s, -1, axis = 1) - numpy.delete(fiber_array.fiber_array_s, 0, axis = 1)
l1 = numpy.sqrt(numpy.square(array_orien_r) + numpy.square(array_orien_a) + numpy.square(array_orien_s))
array_orien_r = array_orien_r / l1
array_orien_a = array_orien_a / l1
array_orien_s = array_orien_s / l1

fiber_orien_r = numpy.delete(fiber.r, -1) - numpy.delete(fiber.r, 0)
fiber_orien_a = numpy.delete(fiber.a, -1) - numpy.delete(fiber.a, 0)
fiber_orien_s = numpy.delete(fiber.s, -1) - numpy.delete(fiber.s, 0)
l2 = numpy.sqrt(numpy.square(fiber_orien_r) + numpy.square(fiber_orien_a) + numpy.square(fiber_orien_s))
fiber_orien_r = fiber_orien_r / l2
fiber_orien_a = fiber_orien_a / l2
fiber_orien_s = fiber_orien_s / l2

array_length = numpy.sum(l1, 1)
fiber_length = numpy.sum(l2)
length = array_length - fiber_length

orientation = array_orien_r * fiber_orien_r + array_orien_a * fiber_orien_a + array_orien_s * fiber_orien_s
orientation = - orientation + 1
orientation = numpy.max(orientation, 1)

# threshold if requested
if threshold:
Expand Down Expand Up @@ -184,7 +218,7 @@ def _fiber_distance_internal_use(fiber, fiber_array, threshold=0, distance_metho
raise Exception("unknown distance method")


return distance
return distance, orientation

def rectangular_frechet_distances(input_vtk_polydata_m,input_vtk_polydata_n):
# Computes distance matrix (nxm) for all n+m fibers in input
Expand Down