Shapely检测不存在的交点问题排查求助

• 使用Shapely检测蓝色原始精确轮廓与绿色法线直线的交点时，偶尔出现不在任一直线上的异常点，严重干扰表面粗糙度计算。
• 已排查坐标系对齐问题（通过绘图和坐标打印验证），正常检测时坐标值合理，但异常值问题仍存在。
• 元素说明：蓝色=原始精确轮廓，红色=高斯滤波后的基线轮廓，绿色=基于基线点生成的法线直线，品红色=检测出的交点，黄色=基线对应位置。
• 核心疑问：该异常是Shapely的问题，还是代码存在未排查的错误？

import numpy as np
import cv2
import math 
from scipy.ndimage import gaussian_filter1d
import shapely
import matplotlib.pyplot as plt


def euclidDist(x1, y1, x2, y2):
    return math.sqrt((x2-x1)**2 + (y2-y1)**2)


def createNormalLine(x,y, dx, dy):
    xlin = np.linspace(-100,100,50)
    ylin = np.linspace(-100,100,50)
    
    return xlin*-dy  + x, ylin*dx+ y


def isInImage(r, x, y):
    xIn, yIn = False, False 
    
    if( x >= r[0] and x <= r[0]+r[2]):
        xIn = True
    
    if y >= r[1] and y <= r[1]+r[3]:
        yIn = True
    
    return xIn and yIn


#pixel to length ratio
scale = 2.27

# Process image
img = cv2.imread("/Users/v.jayaweera/Pictures/FindingEdgesCutContour/OneFile/14.3a.png")
imgGray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
constant= cv2.copyMakeBorder(imgGray,1,1,1,1,cv2.BORDER_CONSTANT,value=255)
_, imgBin = cv2.threshold(constant, 70, 255, 0)
imgBin = cv2.bitwise_not(imgBin)

cont, hier = cv2.findContours(imgBin, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE) 


#crop image
cv2.namedWindow('Select ROI', cv2.WINDOW_NORMAL)
cv2.resizeWindow('Select ROI', int(img.shape[1]*0.25), int(img.shape[0]*0.25))
cropCoord = cv2.selectROI('Select ROI', img, showCrosshair=True)
cv2.destroyWindow('Select ROI')

#Get contours from image 
large_contours = []
areaThres = 50000
xscipy = []
yscipy = []
realCont = []
sig = 50
rad = 15

#Filter and smooth contours
for k in cont:
    x = k[:,0,0]*scale
    y = k[:,0,1]*scale

    
    if(cv2.contourArea(k) > areaThres):
        plt.plot(x,y,'b.')
        plt.plot(gaussian_filter1d(x, sig),gaussian_filter1d(y, sig), 'r.')
        plt.show()
        large_contours.append(k)
        
        #if contours is enclosed
        if(euclidDist(x[0], y[0],x[-1],y[-1]) <= rad):
            xscipy.extend(gaussian_filter1d(x, sig, mode='wrap'))
            yscipy.extend(gaussian_filter1d(y, sig, mode='wrap'))
        else:
            xscipy.extend(gaussian_filter1d(x, sig))
            yscipy.extend(gaussian_filter1d(y, sig))
         
        rCont = np.squeeze(k, axis=1)
        realCont.extend(rCont)        


realCont = np.array(realCont)*scale

#remove edges of smooth contour 
smoothContX = []
smoothContY = []
border = 20*scale

for i in range(len(xscipy)):
    if(xscipy[i] > border and yscipy[i] > border and xscipy[i] < (img.shape[1]-border)
    and yscipy[i] < (img.shape[0]-border) and isInImage(cropCoord, xscipy[i], yscipy[i])):
        smoothContX.append(xscipy[i])
        smoothContY.append(yscipy[i])


#Calculate surface roughness
dx = np.diff(smoothContX)
dy = np.diff(smoothContY)
polyGon = shapely.geometry.LineString(realCont)
distanceE = []

for i in range(len(dx)):
    xs, ys = createNormalLine(smoothContX[i], smoothContY[i], dx[i], dy[i])

    stack = np.stack((xs,ys), axis=-1)
    line = shapely.geometry.LineString(stack)
    
    #Shapely intersection
    interPoints = polyGon.intersection(line)
    mx, my = None,None
        
    if(interPoints.geom_type == "MultiPoint"):
        interPoints = interPoints.geoms
        mx = interPoints[0].x
        my = interPoints[0].y
    elif (interPoints.geom_type == "Point"):
        mx = interPoints.x
        my = interPoints.y
        
    if(mx is not None and polyGon.intersects(line)):
        euD = euclidDist(smoothContX[i], smoothContY[i], mx, my)
        plt.plot(*polyGon.xy, 'b.')
        plt.plot(smoothContX, smoothContY, 'r.')
        plt.plot(xs, ys,'g-')
        plt.plot(mx, my, 'm.')
        plt.plot(smoothContX[i], smoothContY[i], 'y.')
        plt.show()
        input('Press Enter to continue...')
        print(smoothContX[i], smoothContY[i], mx, my, ' EUD ', euD)
        distanceE.append(euD)
        
    

print(np.average(distanceE))

• 原始图片：
• 选择的ROI：

核心错误（非Shapely问题） #

1. 多轮廓合并生成虚假线段
   代码中将所有面积达标的轮廓点合并到realCont，再转换为单个LineString。这会强行连接多个独立轮廓的首尾，生成原始图像中不存在的虚假线段。法线直线极易与这些虚假线段相交，得到完全不符合预期的异常交点。

2. 平滑点集跨轮廓，法线方向无效
   平滑后的xscipy/yscipy同样合并了多轮廓的点，np.diff计算的dx/dy会包含不同轮廓之间的差分，得到的切线方向完全错误，基于此生成的法线自然也不合法，进一步加剧异常交点问题。

修复建议 #

1. 单个轮廓单独处理
   遍历每个大轮廓，分别对单个轮廓做平滑、生成基线、计算法线，再求该法线与当前轮廓的交点，禁止合并多轮廓数据。修改后的核心逻辑示例：

   # 替换原有的轮廓处理和粗糙度计算部分
   for k in cont:
       x = k[:,0,0]*scale
       y = k[:,0,1]*scale
       
       if cv2.contourArea(k) <= areaThres:
           continue
       
       # 处理单个轮廓的平滑
       if euclidDist(x[0], y[0], x[-1], y[-1]) <= rad:
           smooth_x = gaussian_filter1d(x, sig, mode='wrap')
           smooth_y = gaussian_filter1d(y, sig, mode='wrap')
       else:
           smooth_x = gaussian_filter1d(x, sig)
           smooth_y = gaussian_filter1d(y, sig)
       
       # 过滤平滑后点的边界（仅保留ROI内的点）
       smoothContX = []
       smoothContY = []
       for sx, sy in zip(smooth_x, smooth_y):
           if (sx > border and sy > border and sx < (img.shape[1]-border)
               and sy < (img.shape[0]-border) and isInImage(cropCoord, sx, sy)):
               smoothContX.append(sx)
               smoothContY.append(sy)
       
       # 生成当前轮廓的LineString
       real_cont = np.squeeze(k, axis=1)*scale
       polyGon = shapely.geometry.LineString(real_cont)
       
       # 计算当前轮廓的粗糙度
       dx = np.diff(smoothContX)
       dy = np.diff(smoothContY)
       distanceE = []
       
       for i in range(len(dx)):
           xs, ys = createNormalLine(smoothContX[i], smoothContY[i], dx[i], dy[i])
           line = shapely.geometry.LineString(np.stack((xs, ys), axis=-1))
           interPoints = polyGon.intersection(line)
           
           mx, my = None, None
           if interPoints.geom_type == "MultiPoint":
               # 选择距离基线点最近的交点
               min_dist = float('inf')
               for p in interPoints.geoms:
                   dist = euclidDist(smoothContX[i], smoothContY[i], p.x, p.y)
                   if dist < min_dist:
                       min_dist = dist
                       mx, my = p.x, p.y
           elif interPoints.geom_type == "Point":
               mx, my = interPoints.x, interPoints.y
           
           if mx is not None:
               # 验证交点是否在原始轮廓的线段上（过滤浮点误差导致的虚假点）
               point = shapely.geometry.Point(mx, my)
               if polyGon.distance(point) < 1e-6:
                   euD = euclidDist(smoothContX[i], smoothContY[i], mx, my)
                   distanceE.append(euD)
       
       if distanceE:
           print(f"当前轮廓平均粗糙度: {np.average(distanceE)}")
   

2. 优化法线生成（可选）
   原createNormalLine逻辑正确，但可将线段延伸至图像边界，避免因线段过短漏检交点：

   def createNormalLine(x, y, dx, dy, img_width, img_height, scale):
       # 计算法线单位向量
       norm = math.hypot(dx, dy)
       if norm < 1e-6:
           return np.array([x]), np.array([y])
       nx = -dy / norm
       ny = dx / norm
       
       # 计算线段两端点，延伸到图像边界
       t_vals = [
           (0 - x)/nx, (0 - y)/ny,
           (img_width*scale - x)/nx, (img_height*scale - y)/ny
       ]
       t1 = max(t_vals)
       t2 = min(t_vals)
       
       xs = np.array([x + t2*nx, x + t1*nx])
       ys = np.array([y + t2*ny, y + t1*ny])
       return xs, ys
   

3. 添加交点有效性验证
   利用Shapely的distance方法，验证交点到原始轮廓的距离是否在浮点精度容差内（如1e-6），过滤因浮点误差产生的虚假交点。

内容的提问来源于stack exchange，提问作者user22034139