From ac9e65df6179f6a40d541f0b1a45b81d7226a259 Mon Sep 17 00:00:00 2001 From: Huan Luo Date: Fri, 21 Sep 2018 13:47:29 -0400 Subject: [PATCH] Another way to add orientation information --- bin/wm_cluster_atlas.py | 4 +- whitematteranalysis/cluster.py | 59 ++++++++++++++++++++---------- whitematteranalysis/fibers.pyx | 0 whitematteranalysis/similarity.pyx | 48 ++++++++++++++++++++---- 4 files changed, 82 insertions(+), 29 deletions(-) mode change 100644 => 100755 whitematteranalysis/cluster.py mode change 100644 => 100755 whitematteranalysis/fibers.pyx mode change 100644 => 100755 whitematteranalysis/similarity.pyx diff --git a/bin/wm_cluster_atlas.py b/bin/wm_cluster_atlas.py index 3546e0db..8e70432c 100755 --- a/bin/wm_cluster_atlas.py +++ b/bin/wm_cluster_atlas.py @@ -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 diff --git a/whitematteranalysis/cluster.py b/whitematteranalysis/cluster.py old mode 100644 new mode 100755 index 39fe7bee..b605bac8 --- a/whitematteranalysis/cluster.py +++ b/whitematteranalysis/cluster.py @@ -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', @@ -696,7 +706,7 @@ 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': @@ -704,7 +714,7 @@ def _rectangular_similarity_matrix(input_polydata_n, input_polydata_m, threshold 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 @@ -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', @@ -769,7 +788,7 @@ 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': @@ -777,7 +796,7 @@ def _pairwise_similarity_matrix(input_polydata, threshold, sigma, 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)) @@ -785,7 +804,7 @@ def _pairwise_similarity_matrix(input_polydata, threshold, sigma, # 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 " ERROR: On-diagonal elements are not all 1.0." print" Minimum on-diagonal value:", numpy.min(numpy.diag(similarity_matrix)) diff --git a/whitematteranalysis/fibers.pyx b/whitematteranalysis/fibers.pyx old mode 100644 new mode 100755 diff --git a/whitematteranalysis/similarity.pyx b/whitematteranalysis/similarity.pyx old mode 100644 new mode 100755 index d6a2d934..2a2805e0 --- a/whitematteranalysis/similarity.pyx +++ b/whitematteranalysis/similarity.pyx @@ -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 @@ -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. @@ -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: @@ -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