难以计算一组旋转和移动的椭圆的斜率，有时是倒置的，有时是完全错误的

Question

我正在使用 OpenCV-Python 使椭圆适合水滴的形状。然后我选择一条线，代表液滴停留的表面。我计算曲面和椭圆相交处的切线，得到液滴的接触角。
它在大多数情况下都有效，但在某些情况下，切线会上下颠倒或完全错误。看来切线的斜率计算失败了。
谁能告诉我为什么会这样？

在这里你可以看到它应该是什么样子（表面在 y=250）：

这是当我选择 y=47 的表面水平时的结果：

我做了一些研究，我需要在椭圆旋转 phi 之前检测两个 maj_ax, min_ax 中的哪一个与 x-Axis 平行，否则斜率计算算法会失败.

我做错了什么？

这是一个最小的可重现示例：

from math import cos, sin, pi, sqrt, tan, atan2, radians
import cv2

class Droplet():
    def __init__(self):
        self.is_valid = False
        self.angle_l = 0
        self.angle_r = 0
        self.center = (0,0)
        self.maj = 0
        self.min = 0
        self.phi = 0.0
        self.tilt_deg = 0
        self.foc_pt1 = (0,0)
        self.foc_pt2 = (0,0)
        self.tan_l_m = 0
        self.int_l = (0,0)
        self.line_l = (0,0,0,0)
        self.tan_r_m = 0
        self.int_r = (0,0)
        self.line_r = (0,0,0,0)
        self.base_diam = 0
        
def evaluate_droplet(img, y_base) -> Droplet:
    drplt = Droplet()
    crop_img = img[:y_base,:]
    shape = img.shape
    height = shape[0]
    width = shape[1]
    # values only for 8bit images!
    bw_edges = cv2.Canny(crop_img, 76, 179)

    contours, hierarchy = cv2.findContours(bw_edges, cv2.RETR_TREE, cv2.CHAIN_APPROX_NONE)

    if len(contours) == 0:
        raise ValueError('No contours found!')

    edge = max(contours, key=cv2.contourArea)
    (x0,y0),(maj_ax,min_ax),phi_deg = cv2.fitEllipse(edge)

    phi = radians(phi_deg) # to radians
    a = maj_ax/2
    b = min_ax/2

    intersection = calc_intersection_line_ellipse((x0,y0,a,b,phi),(0,y_base))
    if intersection is None:
        raise ValueError('No intersections found')
    # select left and right intersection points
    x_int_l = min(intersection)
    x_int_r = max(intersection)

    foc_len = sqrt(abs(a**2 - b**2))

    # calc slope and angle of tangent
    m_t_l = calc_slope_of_ellipse((x0,y0,a,b,phi), x_int_l, y_base)
    angle_l = pi - atan2(m_t_l,1)
 
    m_t_r = calc_slope_of_ellipse((x0,y0,a,b,phi), x_int_r, y_base)
    angle_r = atan2(m_t_r,1) + pi

    drplt.angle_l = angle_l
    drplt.angle_r = angle_r
    drplt.maj = maj_ax
    drplt.min = min_ax
    drplt.center = (x0, y0)
    drplt.phi = phi
    drplt.tilt_deg = phi_deg
    drplt.tan_l_m = m_t_l
    drplt.tan_r_m = m_t_r
    drplt.line_l = (int(round(x_int_l - (int(round(y_base))/m_t_l))), 0, int(round(x_int_l + ((height - int(round(y_base)))/m_t_l))), int(round(height)))
    drplt.line_r = (int(round(x_int_r - (int(round(y_base))/m_t_r))), 0, int(round(x_int_r + ((height - int(round(y_base)))/m_t_r))), int(round(height)))
    drplt.int_l = (x_int_l, y_base)
    drplt.int_r = (x_int_r, y_base)
    drplt.foc_pt1 = (x0 + foc_len*cos(phi), y0 + foc_len*sin(phi))
    drplt.foc_pt2 = (x0 - foc_len*cos(phi), y0 - foc_len*sin(phi))
    drplt.base_diam = x_int_r - x_int_l
    drplt.is_valid = True
    
    # draw ellipse and lines
    img = cv2.drawContours(img,contours,-1,(100,100,255),2)
    img = cv2.drawContours(img,edge,-1,(255,0,0),2)
    img = cv2.ellipse(img, (int(round(x0)),int(round(y0))), (int(round(a)),int(round(b))), int(round(phi*180/pi)), 0, 360, (255,0,255), thickness=1, lineType=cv2.LINE_AA)
    y_int = int(round(y_base))
    img = cv2.line(img, (int(round(x_int_l - (y_int/m_t_l))), 0), (int(round(x_int_l + ((height - y_int)/m_t_l))), int(round(height))), (255,0,255), thickness=1, lineType=cv2.LINE_AA)
    img = cv2.line(img, (int(round(x_int_r - (y_int/m_t_r))), 0), (int(round(x_int_r + ((height - y_int)/m_t_r))), int(round(height))), (255,0,255), thickness=1, lineType=cv2.LINE_AA)
    img = cv2.ellipse(img, (int(round(x_int_l)),y_int), (20,20), 0, 0, -int(round(angle_l*180/pi)), (255,0,255), thickness=1, lineType=cv2.LINE_AA)
    img = cv2.ellipse(img, (int(round(x_int_r)),y_int), (20,20), 0, 180, 180 + int(round(angle_r*180/pi)), (255,0,255), thickness=1, lineType=cv2.LINE_AA)
    img = cv2.line(img, (0,y_int), (width, y_int), (255,0,0), thickness=2, lineType=cv2.LINE_AA)
    img = cv2.putText(img, '<' + str(round(angle_l*180/pi,1)), (5,y_int-5), cv2.FONT_HERSHEY_COMPLEX, .5, (0,0,0))
    img = cv2.putText(img, '<' + str(round(angle_r*180/pi,1)), (width - 80,y_int-5), cv2.FONT_HERSHEY_COMPLEX, .5, (0,0,0))
    cv2.imshow('Test',img)
    cv2.waitKey(0)
    return drplt

def calc_intersection_line_ellipse(ellipse_pars, line_pars):
    """
    calculates intersection(s) of an ellipse with a line  
    :param ellipse_pars: tuple of (x0,y0,a,b,phi): x0,y0 center of ellipse; a,b sem-axis of ellipse; phi tilt rel to x axis  
    :param line_pars: tuple of (m,t): m is the slope and t is intercept of the intersecting line  
    :returns: x-coordinate(s) of intesection as list or float or none if none found
    """
    ## -->> http://quickcalcbasic.com/ellipse%20line%20intersection.pdf
    (x0, y0, h, v, phi) = ellipse_pars
    (m, t) = line_pars
    y = t - y0
    try:
        a = v**2 * cos(phi)**2 + h**2 * sin(phi)**2
        b = 2*y*cos(phi)*sin(phi) * (v**2 - h**2)
        c = y**2 * (v**2 * sin(phi)**2 + h**2 * cos(phi)**2) - (h**2 * v**2)
        det = b**2 - 4*a*c
        if det > 0:
            x1 = int(round((-b - sqrt(det))/(2*a) + x0))
            x2 = int(round((-b + sqrt(det))/(2*a) + x0))
            return x1,x2
        elif det == 0:
            x = int(round(-b / (2*a)))
            return x
        else:
            return None
    except Exception as ex:
        raise ex
 
def calc_slope_of_ellipse(ellipse_pars, x, y):
    """
    calculates the slope of the tangent at point x,y, the point needs to be on the ellipse!
    :param ellipse_params: tuple of (x0,y0,a,b,phi): x0,y0 center of ellipse; a,b sem-axis of ellipse; phi tilt rel to x axis  
    :param x: x-coord where the slope will be calculated  
    :returns: the slope of the tangent 
    """
    (x0, y0, a, b, phi) = ellipse_pars
    # transform to non-rotated ellipse
    x_rot = (x - x0)*cos(phi) + (y - y0)*sin(phi)
    y_rot = (x - x0)*sin(phi) + (y - y0)*cos(phi)
    m_rot = -(b**2 * x_rot)/(a**2 * y_rot) # slope of tangent to unrotated ellipse
    #rotate tangent line back to angle of the rotated ellipse
    m_tan = tan(atan2(m_rot,1) + phi)

    return m_tan

if __name__ == "__main__":
    im = cv2.imread('untitled1.png')
    # any value below 250 is just the droplet without the substrate
    drp = evaluate_droplet(im, 250)

原图：

Answer 1

我写错了calc_slope_ellipse:
x_rot = (x - x0)*cos(phi) + (y - y0)*sin(phi)
应该是
x_rot = (x - x0)*cos(phi) - (y - y0)*sin(phi)
这修复了 y=47 处斜率的错误符号。
我换了atan2:

m_rot = -(b**2 * x_rot)/(a**2 * y_rot) # slope of tangent to unrotated ellipse
#rotate tangent line back to angle of the rotated ellipse
m_tan = tan(atan2(m_rot,1) + phi)

和

tan_a = x_rot/a**2
tan_b = y_rot/b**2
#rotate tangent line back to angle of the rotated ellipse
tan_a_r = tan_a*cos(phi) + tan_b*sin(phi)
tan_b_r = tan_b*cos(phi) - tan_a*sin(phi)
m_tan = - (tan_a_r / tan_b_r)

这修复了某些情况下的奇怪行为 (y=62)。

完整 fcn:

def calc_slope_of_ellipse(ellipse_pars, x, y):
    (x0, y0, a, b, phi) = ellipse_pars

    # transform to non-rotated ellipse centered to origin
    x_rot = (x - x0)*cos(phi) - (y - y0)*sin(phi)
    y_rot = (x - x0)*sin(phi) + (y - y0)*cos(phi)
    # Ax + By = C
    tan_a = x_rot/a**2
    tan_b = y_rot/b**2
    #rotate tangent line back to angle of the rotated ellipse
    tan_a_r = tan_a*cos(phi) + tan_b*sin(phi)
    tan_b_r = tan_b*cos(phi) - tan_a*sin(phi)
    m_tan = - (tan_a_r / tan_b_r)

    return m_tan

难以计算一组旋转和移动的椭圆的斜率，有时是倒置的，有时是完全错误的

Difficulty calculating slope for a set of rotated and shifted ellipses, sometimes inverted, sometimes completely wrong

python

math

trigonometry

curve-fitting

opencv-python