问题描述
我正在使用OpenCV-Python使椭圆适合液滴的形状。
然后，我选择一条线，该线代表液滴所停留的表面。
我计算表面和椭圆相交处的切线以获得液滴的接触角。
它在大多数情况下都有效，但是在某些情况下，切线会上下翻转或出现错误。
切线斜率的计算似乎失败了。
有人可以告诉我为什么会这样吗？
在这里您可以看到它的外观（y = 250处的表面）：
这是当我选择y = 47的表面水平时的结果：
我做了一些研究，我需要检测两个maj_ax,min_ax中的哪个平行于x轴，然后椭圆phi旋转，否则斜率计算算法将失败。
我在做什么错了？
这是一个最小的可重现示例：
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)
self.tan_r_m = 0
self.int_r = (0,0)
self.line_r = (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,x_int_l,y_base)
angle_l = pi - atan2(m_t_l,1)
m_t_r = calc_slope_of_ellipse((x0,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))),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))),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,(255,0),2)
img = cv2.ellipse(img,(int(round(x0)),int(round(y0))),(int(round(a)),int(round(b))),int(round(phi*180/pi)),360,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))),(int(round(x_int_l + ((height - y_int)/m_t_l))),int(round(height))),lineType=cv2.LINE_AA)
img = cv2.line(img,(int(round(x_int_r - (y_int/m_t_r))),(int(round(x_int_r + ((height - y_int)/m_t_r))),lineType=cv2.LINE_AA)
img = cv2.ellipse(img,(int(round(x_int_l)),y_int),(20,20),-int(round(angle_l*180/pi)),(int(round(x_int_r)),180,180 + int(round(angle_r*180/pi)),(width,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))
img = cv2.putText(img,'<' + str(round(angle_r*180/pi,(width - 80,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,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,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,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,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)
原始图片：
解决方法
-
我在
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