diff --git a/ntrfc/turbo/pointcloud_methods.py b/ntrfc/turbo/pointcloud_methods.py
index be56ca0ed18e469b84c9bce5f2d7837511185bd4..86fd377fe6111c4eb738e7018873f12222ee8866 100644
--- a/ntrfc/turbo/pointcloud_methods.py
+++ b/ntrfc/turbo/pointcloud_methods.py
@@ -4,11 +4,11 @@ import pyvista as pv
from ntrfc.geometry.line import refine_spline, polyline_from_points
-def midline_from_sides(ps_poly, ss_poly):
+def midline_from_sides(ps_poly, ss_poly,res =100):
x_ps, y_ps = ps_poly.points[::, 0], ps_poly.points[::, 1]
x_ss, y_ss = ss_poly.points[::, 0], ss_poly.points[::, 1]
z = ps_poly.points[0][2]
- midsres = 100
+ midsres = res
if x_ps[0] > x_ps[-1]:
ax, ay = refine_spline(x_ps[::-1], y_ps[::-1], midsres)
else:
@@ -23,12 +23,16 @@ def midline_from_sides(ps_poly, ss_poly):
return midsPoly
-def extractSidePolys(ind_2, sortedPoly):
+def extractSidePolys(ind_1, ind_2, sortedPoly):
# xs, ys = list(sortedPoly.points[::, 0]), list(sortedPoly.points[::, 1])
indices = np.arange(0, sortedPoly.number_of_points)
- side_one_idx = indices[0:ind_2 + 1]
- side_two_idx = np.concatenate([[0], indices[ind_2:][::-1]])
+ if ind_2 > ind_1:
+ side_one_idx = indices[ind_1:ind_2 + 1]
+ side_two_idx = np.concatenate((indices[:ind_1 + 1][::-1], indices[ind_2:][::-1]))
+ elif ind_1 > ind_2:
+ side_one_idx = indices[ind_2:ind_1 + 1]
+ side_two_idx = np.concatenate((indices[:ind_2 + 1][::-1], indices[ind_1:][::-1]))
side_one = pv.PolyData()
for i in side_one_idx:
diff --git a/ntrfc/turbo/profile_tele_extraction.py b/ntrfc/turbo/profile_tele_extraction.py
index 78672afd9019c39bd92743d2afb39a6ea031aedf..89192b0bdccbfbd536067230bdbb200d6af3f7c4 100644
--- a/ntrfc/turbo/profile_tele_extraction.py
+++ b/ntrfc/turbo/profile_tele_extraction.py
@@ -8,9 +8,10 @@ from scipy.spatial import Voronoi
from scipy.spatial.distance import cdist
from shapely.geometry import Polygon, LineString
from skimage.morphology import skeletonize
-
+from ntrfc.geometry.line import lines_from_points
+from ntrfc.turbo.pointcloud_methods import extractSidePolys, midline_from_sides
from ntrfc.math.vectorcalc import findNearest, vecDir
-
+from tqdm import tqdm
def detect_inliers_tukey(data):
"""
@@ -74,6 +75,10 @@ def clean_sites(voronoi_sites, skeletonize_sites):
inliers_turkey = detect_inliers_tukey(min_distances)
inliers_mad = detect_inliers_mad(min_distances)
inliers_zscore = detect_inliers_zscore(min_distances)
+ # Compute the intersection with the polygon of all rays
+ intersections_positive = [LineString([centers[i], ray_positive[i]]).intersection(shapelypoly) for i in range(len(centers))]
+ intersections_negative = [LineString([centers[i], ray_negative[i]]).intersection(shapelypoly) for i in range(len(centers))]
+
inlier_indices = np.where(inliers_turkey * inliers_zscore * inliers_mad)[0]
@@ -82,29 +87,90 @@ def clean_sites(voronoi_sites, skeletonize_sites):
def extract_vk_hk(sortedPoly: pv.PolyData) -> (int, int):
- voronoires = 16000
- skeletonres = 6000
+ voronoires = 64000
points_orig = sortedPoly.points
+ npts = len(points_orig)
points_2d_closed_refined_voronoi = pointcloud_to_profile(points_orig, voronoires)
- points_2d_closed_refined_skeletonize = pointcloud_to_profile(points_orig, skeletonres)
- voronoi_sites = voronoi_skeleton_sites(points_2d_closed_refined_voronoi)
- skeletonize_sites = skeletonize_skeleton_sites(points_2d_closed_refined_skeletonize)
+ #create sheply polygon from refined pointset
+ shapelypoly = Polygon(points_2d_closed_refined_voronoi)
+ raylength = 1000
+ # todo; estimate le and te regions, only check combinations of these
+
+ min_x = np.min(points_orig[:, 0])
+ max_x = np.max(points_orig[:, 0])
+
+ le_est_xmax = min_x + 0.01 * (max_x - min_x)
+ te_est_xmin = max_x - 0.01 * (max_x - min_x)
+
+ le_indices = np.where(points_orig[:, 0] < le_est_xmax)[0]
+ te_indices = np.where(points_orig[:, 0] > te_est_xmin)[0]
+
+
+ combinations = [(i, j) for i in le_indices for j in te_indices]
+
+
+
+ rms_errors = []
+ #max_errors = []
+
+ midsp = []
+ for i, j in tqdm(combinations):
+ ind_1_start = i
+ ind_2_start = j
+
+ ps_pv, ss_pv = extractSidePolys(ind_1_start,ind_2_start, sortedPoly)
+
+ mid = lines_from_points(midline_from_sides(ps_pv, ss_pv,res=200).points)
+
+ mid_normals = mid.extrude([0,0,1]).compute_normals().slice("z").compute_cell_sizes()
+
+ normals = mid_normals["Normals"]
+ lengths = mid_normals["Length"]
+ centers = mid_normals.cell_centers().points
+
+ ray_positive = centers + raylength * normals
+ ray_negative = centers - raylength * normals
+
+ # Compute the intersection with the polygon of all rays
+ intersections_positive = [LineString([centers[i], ray_positive[i]]).intersection(shapelypoly) for i in range(len(centers))]
+ intersections_negative = [LineString([centers[i], ray_negative[i]]).intersection(shapelypoly) for i in range(len(centers))]
+
+ # compute length of each intersection
+ lengths_positive = np.array([line.length for line in intersections_positive])
+ lengths_negative = np.array([line.length for line in intersections_negative])
+
+ differences = np.abs(lengths_positive - lengths_negative) /((lengths_positive + lengths_negative)/2)
+
+ rms_error = np.sqrt(np.mean(differences**2))
+
+ rms_errors.append(rms_error)
+ # max_errors.append(np.max(differences))
+
+ # mid_normals["diff"]=differences
+ # midsp.append(mid_normals)
- valid_midline = clean_sites(voronoi_sites, skeletonize_sites)
- sort_indices = np.argsort(valid_midline[:, 0])
- valid_midline_sorted = valid_midline[sort_indices]
+ # verbose = False
+ # if verbose:
+ # p = pv.Plotter()
+ # mid_normals = mid_normals.translate([0, 0, -mid_normals.bounds[5]])
+ # p.add_mesh(mid_normals, scalars="diff")
+ # p.add_mesh(ps_pv, color="red")
+ # p.add_mesh(ss_pv, color="green")
+ # for ineg in intersections_negative:
+ # p.add_mesh(pv.Line(*np.array(ineg.coords) - ineg.coords[0][2]), color="blue")
+ # for ipos in intersections_positive:
+ # p.add_mesh(pv.Line(*np.array(ipos.coords) - ipos.coords[0][2]), color="green")
+ # p.add_mesh(pv.PolyData(np.zeros(3)), color="black")
+ # p.view_xy()
+ # p.show_axes()
+ # p.show()
- smoothing = 0
- u_new = np.arange(0, 1, 1 / 1024)
- (tck, u), fp, ier, msg = splprep((valid_midline_sorted[::, 0], valid_midline_sorted[::, 1]), u=None, per=0, k=3,
- s=smoothing,
- full_output=True) # s = optional parameter (default used here)
- x_new, y_new = splev(u_new, tck, der=0)
- le_ind, te_ind, skeletonline_complete = skeletonline_completion(2, points_orig, points_2d_closed_refined_voronoi,
- np.stack([x_new[1:-1], y_new[1:-1]]).T)
+ inds = combinations[np.argmin(rms_errors)]
+ le_ind = inds[np.argmin(points_orig[inds][:, 0])]
+ te_ind = inds[np.argmax(points_orig[inds][:, 0])]
return le_ind, te_ind