[FIXED] An algorithm to sort top and bottom slices of curved surfaces

Issue

I try to do:

1. Cut STL file https://www.dropbox.com/s/pex20yqfgmxgt0w/wing_fish.stl?dl=0 at Z-coordinate using PyVsita https://docs.pyvista.org/ )
2. Extract point's coordinates X, Y at given section Z
3. Sort points to Upper and Down groups for further manipulation

Here is my code:

import pyvista as pv
import matplotlib.pylab as plt
import numpy as np
import math

mesh = pv.read('wing_fish.stl')
z_slice = [0, 0, 1]     # normal to cut at

single_slice = mesh.slice(normal=z_slice, origin=[0, 0, 200]) #  slicing
a = single_slice.points                    #  choose only points

# p = pv.Plotter()          #show section
# p.add_mesh(single_slice)
# p.show()

a = a[a[:,0].astype(float).argsort()]       #  sort all points by Х coord
#  X min of all points
x0 = a[0][0]
#  Y min of all points
y0 = a[0][1] 

#  X tail 1 of 2 
xn = a[-1][0]
#  Y tail 1 of 2 
yn = a[-1][1]

#  X tail 2 of 2
xn2 = a[-2][0]
#  Y tail 2 of 2
yn2 = a[-2][1]

def line_y(x, x0, y0, xn, yn):
    # return y coord at arbitary x coord of x0, y0  xn, yn LINE
    return ((x - x0)*(yn-y0))/(xn-x0)+y0

def line_c(x0, y0, xn, yn):
    # return x, y middle points of LINE
    xc = (x0+xn)/2
    yc = (y0+yn)/2
    return xc, yc

def chord(P1, P2):
    return math.sqrt((P2[1] - P1[1])**2 + (P2[0] - P1[0])**2)

xc_end, yc_end = line_c(xn, yn, xn2, yn2)   # return midle at trailing edge

midLine = np.array([[x0,y0],[xc_end,yc_end]],dtype='float32')

c_temp_x_d = []
c_temp_y_d = []

c_temp_x_u = []
c_temp_y_u = []

isUp = None
isDown = None

for i in a:
    if i[1] == line_y(i[0], x0=x0, y0=y0, xn=xc_end, yn=yc_end):
        continue

    elif i[1] < line_y(i[0], x0=x0, y0=y0, xn=xc_end, yn=yc_end):
       
        c_temp_y_d.append(i[1])
        c_temp_x_d.append(i[0])
        isDown = True
    else:
        c_temp_y_u.append(i[1])
        c_temp_x_u.append(i[0])
        isUp = True
    
    if len(c_temp_y_d) != 0 and len(c_temp_y_u) != 0:
        print(c_temp_y_d[-1])

plt.plot(c_temp_x_d, c_temp_y_d, label='suppose to be down points')
plt.plot(c_temp_x_u, c_temp_y_u, label='suppose to be upper points')
plt.plot(midLine[:,0], midLine[:,1], label='Chord')
plt.scatter(a[:,0],a[:,1], label='raw points')

plt.legend();plt.grid();plt.show()

What I have:

What I want:

I would highly appreciate for any help and advises! Thanks in advance!

Solution

You are discarding precious connectivity information that is already there in your STL mesh and in your slice!

I couldn't think of a more idiomatic solution within PyVista, but at worst you can take the cell (line) information from the slice and start walking your shape (that is topologically equivalent to a circle) from its left side to its right, and vice versa. Here's one way:

import numpy as np
import matplotlib.pyplot as plt
import pyvista as pv

mesh = pv.read('../wing_fish.stl')
z_slice = [0, 0, 1]     # normal to cut at

single_slice = mesh.slice(normal=z_slice, origin=[0, 0, 200]) #  slicing

# find points with smallest and largest x coordinate
points = single_slice.points
left_ind = points[:, 0].argmin()
right_ind = points[:, 0].argmax()

# sanity check for what we're about to do:
# 1. all cells are lines
assert single_slice.n_cells == single_slice.n_points
assert (single_slice.lines[::3] == 2).all()

# 2. all points appear exactly once as segment start and end
lines = single_slice.lines.reshape(-1, 3)  # each row: [2, i_from, i_to]
assert len(set(lines[:, 1])) == lines.shape[0]

# create an auxiliary dict with from -> to index mappings
conn = dict(lines[:, 1:])

# and a function that walks this connectivity graph
def walk_connectivity(connectivity, start, end):
    this_ind = start
    path_inds = [this_ind]
    while True:
        next_ind = connectivity[this_ind]
        path_inds.append(next_ind)
        this_ind = next_ind
        if this_ind == end:
            # we're done
            return path_inds

# start walking at point left_ind, walk until right_ind
first_side_inds = walk_connectivity(conn, left_ind, right_ind)

# now walk forward for the other half curve
second_side_inds = walk_connectivity(conn, right_ind, left_ind)

# get the point coordinates for plotting
first_side_points = points[first_side_inds, :-1]
second_side_points = points[second_side_inds, :-1]

# plot the two sides
fig, ax = plt.subplots()
ax.plot(*first_side_points.T)
ax.plot(*second_side_points.T)
plt.show()

In order to avoid using an O(n^2) algorithm, I defined an auxiliary dict that maps line segment start indices to end indices. In order for this to work we need some sanity checks, namely that the cells are all simple line segments, and that each segment has the same orientation (i.e. each start point is unique, and each end point is unique). Once we have this it's easy to start from the left edge of your wing profile and walk each line segment until we find the right edge.

The nature of this approach implies that we can't know a priori whether the path from left to right goes on the upper or the lower path. This needs experimentation on your part; name the two paths in whatever way you see fit.

And of course there's always room for fine tuning. For instance, the above implementation creates two paths that both start and end with the left and right-side boundary points of the mesh. If you want the top and bottom curves to share no points, you'll have to adjust the algorithm accordingly. And if the end point is not found on the path then the current implementation will give you an infinite loop with a list growing beyond all available memory. Consider adding some checks in the implementation to avoid this.

Anyway, this is what we get from the above:

Answered By - Andras Deak