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 index 39fe7bee..b605bac8 100644 --- 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 index a05960e3..595e1092 100644 --- a/whitematteranalysis/fibers.pyx +++ b/whitematteranalysis/fibers.pyx @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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() @@ -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 " Line:", lidx, "/", self.number_of_fibers + print " 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: diff --git a/whitematteranalysis/similarity.pyx b/whitematteranalysis/similarity.pyx index d6a2d934..3d4612f3 100644 --- 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,40 @@ 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 +115,28 @@ 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 + + orientation_for = fiber_array.fiber_orien_for_r * fiber.orien_for_r + \ + fiber_array.fiber_orien_for_a * fiber.orien_for_a + \ + fiber_array.fiber_orien_for_s * fiber.orien_for_s + + orientation_back = fiber_array.fiber_orien_back_r * fiber.orien_back_r + \ + fiber_array.fiber_orien_back_a * fiber.orien_back_a + \ + fiber_array.fiber_orien_back_s * fiber.orien_back_s + + orientation_for[orientation_for > 1] = 1 + orientation_for[orientation_for < -1] = -1 + orientation_for = numpy.arccos(orientation_for) + + orientation_back[orientation_back > 1] = 1 + orientation_back[orientation_back < -1] = -1 + orientation_back = numpy.arccos(orientation_back) + + orientation = orientation_for + orientation = numpy.max(orientation, 1) +# orientation = numpy.sum(orientation, 1) +# npts = float(fiber_array.points_per_fiber) +# orientation = orientation / npts # threshold if requested if threshold: @@ -184,7 +215,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