Skip to content

Huan orientation #78

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

Open
luoh15 wants to merge 2 commits into
base: master
Choose a base branch
from
Open

Huan orientation #78

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
c6fcf4b
add orientation information
Sep 21, 2018
74e8563
Fix path to python
Sep 21, 2018
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
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
164 changes: 163 additions & 1 deletion whitematteranalysis/fibers.pyx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,16 @@ class Fiber:
self.r = None
self.a = None
self.s = None

self.orien_for_r = None
self.orien_for_a = None
self.orien_for_s = None

self.orien_back_r = None
self.orien_back_a = None
self.orien_back_s = None


self.points_per_fiber = None
self.hemisphere_percent_threshold = 0.95

Expand All @@ -31,7 +41,15 @@ class Fiber:
fiber.r = self.r[::-1]
fiber.a = self.a[::-1]
fiber.s = self.s[::-1]


fiber.orien_for_r = self.orien_back_r[::-1]
fiber.orien_for_a = self.orien_back_a[::-1]
fiber.orien_for_s = self.orien_back_s[::-1]

fiber.orien_back_r = self.orien_for_r[::-1]
fiber.orien_back_a = self.orien_for_a[::-1]
fiber.orien_back_s = self.orien_for_s[::-1]

fiber.points_per_fiber = self.points_per_fiber

return fiber
Expand All @@ -45,6 +63,14 @@ class Fiber:
fiber.r = - self.r
fiber.a = self.a
fiber.s = self.s

fiber.orien_for_r = -self.orien_for_r
fiber.orien_for_a = self.orien_for_a
fiber.orien_for_s = self.orien_for_s

fiber.orien_back_r = -self.orien_back_r
fiber.orien_back_a = self.orien_back_a
fiber.orien_back_s = self.orien_back_s

fiber.points_per_fiber = self.points_per_fiber

Expand Down Expand Up @@ -118,6 +144,16 @@ class FiberArray:
self.fiber_array_r = None
self.fiber_array_a = None
self.fiber_array_s = None

# fiber forward orientation
self.fiber_orien_for_r = None
self.fiber_orien_for_a = None
self.fiber_orien_for_s = None

# fiber backward orientation
self.fiber_orien_back_r = None
self.fiber_orien_back_a = None
self.fiber_orien_back_s = None

# output arrays indicating hemisphere/callosal (L,C,R= -1, 0, 1)
self.fiber_hemisphere = None
Expand Down Expand Up @@ -190,6 +226,14 @@ class FiberArray:
fiber.r = self.fiber_array_r[fiber_index, :]
fiber.a = self.fiber_array_a[fiber_index, :]
fiber.s = self.fiber_array_s[fiber_index, :]

fiber.orien_for_r = self.fiber_orien_for_r[fiber_index, :]
fiber.orien_for_a = self.fiber_orien_for_a[fiber_index, :]
fiber.orien_for_s = self.fiber_orien_for_s[fiber_index, :]

fiber.orien_back_r = self.fiber_orien_back_r[fiber_index, :]
fiber.orien_back_a = self.fiber_orien_back_a[fiber_index, :]
fiber.orien_back_s = self.fiber_orien_back_s[fiber_index, :]

fiber.points_per_fiber = self.points_per_fiber

Expand Down Expand Up @@ -346,6 +390,18 @@ class FiberArray:
self.points_per_fiber))
self.fiber_array_s = numpy.zeros((self.number_of_fibers,
self.points_per_fiber))
self.fiber_orien_for_r = numpy.zeros((self.number_of_fibers,
self.points_per_fiber))
self.fiber_orien_for_a = numpy.zeros((self.number_of_fibers,
self.points_per_fiber))
self.fiber_orien_for_s = numpy.zeros((self.number_of_fibers,
self.points_per_fiber))
self.fiber_orien_back_r = numpy.zeros((self.number_of_fibers,
self.points_per_fiber))
self.fiber_orien_back_a = numpy.zeros((self.number_of_fibers,
self.points_per_fiber))
self.fiber_orien_back_s = numpy.zeros((self.number_of_fibers,
self.points_per_fiber))

# loop over lines
input_vtk_polydata.GetLines().InitTraversal()
Expand All @@ -371,11 +427,117 @@ class FiberArray:
ptidx = line_ptids.GetId(int(round(line_index)))

point = inpoints.GetPoint(ptidx)

# if int(round(line_index)) == line_length - 1:
# point_near = inpoints.GetPoint(ptidx - 1)
# l = numpy.sqrt(numpy.square(point[0] - point_near[0]) + \
# numpy.square(point[1] - point_near[1]) + \
# numpy.square(point[2] - point_near[2]))
#
# self.fiber_orien_for_r[lidx, pidx] = (point[0] - point_near[0]) / l
# self.fiber_orien_for_a[lidx, pidx] = (point[1] - point_near[1]) / l
# self.fiber_orien_for_s[lidx, pidx] = (point[2] - point_near[2]) / l
#
# self.fiber_orien_back_r[lidx, pidx] = (point[0] - point_near[0]) / l
# self.fiber_orien_back_a[lidx, pidx] = (point[1] - point_near[1]) / l
# self.fiber_orien_back_s[lidx, pidx] = (point[2] - point_near[2]) / l
#
# elif int(round(line_index)) == 0:
# point_near = inpoints.GetPoint(ptidx + 1)
# l = numpy.sqrt(numpy.square(point_near[0] - point[0]) + \
# numpy.square(point_near[1] - point[1]) + \
# numpy.square(point_near[2] - point[2]))
#
# self.fiber_orien_for_r[lidx, pidx] = (point[0] - point_near[0]) / l
# self.fiber_orien_for_a[lidx, pidx] = (point[1] - point_near[1]) / l
# self.fiber_orien_for_s[lidx, pidx] = (point[2] - point_near[2]) / l
#
# self.fiber_orien_back_r[lidx, pidx] = (point[0] - point_near[0]) / l
# self.fiber_orien_back_a[lidx, pidx] = (point[1] - point_near[1]) / l
# self.fiber_orien_back_s[lidx, pidx] = (point[2] - point_near[2]) / l
# else:
# point_next = inpoints.GetPoint(ptidx + 1)
# point_pre = inpoints.GetPoint(ptidx - 1)
# l1 = numpy.sqrt(numpy.square(point[0] - point_next[0]) + \
# numpy.square(point[1] - point_next[1]) + \
# numpy.square(point[2] - point_next[2]))
#
# l2 = numpy.sqrt(numpy.square(point[0] - point_pre[0]) + \
# numpy.square(point[1] - point_pre[1]) + \
# numpy.square(point[2] - point_pre[2]))
#
# self.fiber_orien_for_r[lidx, pidx] = (point[0] - point_next[0]) / l1
# self.fiber_orien_for_a[lidx, pidx] = (point[1] - point_next[1]) / l1
# self.fiber_orien_for_s[lidx, pidx] = (point[2] - point_next[2]) / l1
#
# self.fiber_orien_back_r[lidx, pidx] = (point[0] - point_pre[0]) / l2
# self.fiber_orien_back_a[lidx, pidx] = (point[1] - point_pre[1]) / l2
# self.fiber_orien_back_s[lidx, pidx] = (point[2] - point_pre[2]) / l2

self.fiber_array_r[lidx, pidx] = point[0]
self.fiber_array_a[lidx, pidx] = point[1]
self.fiber_array_s[lidx, pidx] = point[2]
pidx = pidx + 1

input_vtk_polydata.GetLines().InitTraversal()
line_ptids = vtk.vtkIdList()
inpoints = input_vtk_polydata.GetPoints()

for lidx in range(0, self.number_of_fibers):

input_vtk_polydata.GetLines().GetNextCell(line_ptids)
line_length = line_ptids.GetNumberOfIds()

if self.verbose:
if lidx % 100 == 0:
print "<fibers.py> Line:", lidx, "/", self.number_of_fibers
print "<fibers.py> number of points:", line_length

# loop over the indices that we want and get those points
for pidx in range(0, self.points_per_fiber):

# do nearest neighbor interpolation: round index
if pidx == self.points_per_fiber - 1:
l = numpy.sqrt(numpy.square(self.fiber_array_r[lidx, pidx] - self.fiber_array_r[lidx, pidx - 1]) + \
numpy.square(self.fiber_array_a[lidx, pidx] - self.fiber_array_a[lidx, pidx - 1]) + \
numpy.square(self.fiber_array_s[lidx, pidx] - self.fiber_array_s[lidx, pidx - 1]))

self.fiber_orien_for_r[lidx, pidx] = (self.fiber_array_r[lidx, pidx] - self.fiber_array_r[lidx, pidx - 1]) / l
self.fiber_orien_for_a[lidx, pidx] = (self.fiber_array_a[lidx, pidx] - self.fiber_array_a[lidx, pidx - 1]) / l
self.fiber_orien_for_s[lidx, pidx] = (self.fiber_array_s[lidx, pidx] - self.fiber_array_s[lidx, pidx - 1]) / l

self.fiber_orien_back_r[lidx, pidx] = (self.fiber_array_r[lidx, pidx] - self.fiber_array_r[lidx, pidx - 1]) / l
self.fiber_orien_back_a[lidx, pidx] = (self.fiber_array_a[lidx, pidx] - self.fiber_array_a[lidx, pidx - 1]) / l
self.fiber_orien_back_s[lidx, pidx] = (self.fiber_array_s[lidx, pidx] - self.fiber_array_s[lidx, pidx - 1]) / l

elif pidx == 0:
l = numpy.sqrt(numpy.square(self.fiber_array_r[lidx, pidx] - self.fiber_array_r[lidx, pidx + 1]) + \
numpy.square(self.fiber_array_a[lidx, pidx] - self.fiber_array_a[lidx, pidx + 1]) + \
numpy.square(self.fiber_array_s[lidx, pidx] - self.fiber_array_s[lidx, pidx + 1]))

self.fiber_orien_for_r[lidx, pidx] = (self.fiber_array_r[lidx, pidx + 1] - self.fiber_array_r[lidx, pidx]) / l
self.fiber_orien_for_a[lidx, pidx] = (self.fiber_array_a[lidx, pidx + 1] - self.fiber_array_a[lidx, pidx]) / l
self.fiber_orien_for_s[lidx, pidx] = (self.fiber_array_s[lidx, pidx + 1] - self.fiber_array_s[lidx, pidx]) / l

self.fiber_orien_back_r[lidx, pidx] = (self.fiber_array_r[lidx, pidx + 1] - self.fiber_array_r[lidx, pidx]) / l
self.fiber_orien_back_a[lidx, pidx] = (self.fiber_array_a[lidx, pidx + 1] - self.fiber_array_a[lidx, pidx]) / l
self.fiber_orien_back_s[lidx, pidx] = (self.fiber_array_s[lidx, pidx + 1] - self.fiber_array_s[lidx, pidx]) / l
else:
l1 = numpy.sqrt(numpy.square(self.fiber_array_r[lidx, pidx] - self.fiber_array_r[lidx, pidx + 1]) + \
numpy.square(self.fiber_array_a[lidx, pidx] - self.fiber_array_a[lidx, pidx + 1]) + \
numpy.square(self.fiber_array_s[lidx, pidx] - self.fiber_array_s[lidx, pidx + 1]))

l2 = numpy.sqrt(numpy.square(self.fiber_array_r[lidx, pidx] - self.fiber_array_r[lidx, pidx - 1]) + \
numpy.square(self.fiber_array_a[lidx, pidx] - self.fiber_array_a[lidx, pidx - 1]) + \
numpy.square(self.fiber_array_s[lidx, pidx] - self.fiber_array_s[lidx, pidx - 1]))

self.fiber_orien_for_r[lidx, pidx] = (self.fiber_array_r[lidx, pidx + 1] - self.fiber_array_r[lidx, pidx]) / l1
self.fiber_orien_for_a[lidx, pidx] = (self.fiber_array_a[lidx, pidx + 1] - self.fiber_array_a[lidx, pidx]) / l1
self.fiber_orien_for_s[lidx, pidx] = (self.fiber_array_s[lidx, pidx + 1] - self.fiber_array_s[lidx, pidx]) / l1

self.fiber_orien_back_r[lidx, pidx] = (self.fiber_array_r[lidx, pidx] - self.fiber_array_r[lidx, pidx - 1]) / l2
self.fiber_orien_back_a[lidx, pidx] = (self.fiber_array_a[lidx, pidx] - self.fiber_array_a[lidx, pidx - 1]) / l2
self.fiber_orien_back_s[lidx, pidx] = (self.fiber_array_s[lidx, pidx] - self.fiber_array_s[lidx, pidx - 1]) / l2

# initialize hemisphere info
if self.hemispheres:
Expand Down
Loading