Подбор и измерение функций изображения с помощью cv2 ⇐ Python

[Anonymous]

Я провожу некоторые измерения на капельнице, находящейся на ограниченной поверхности, на avi-видео после обработки. У меня есть два видео: одно для эллиптической шапки, а другое для нижней части. Основная цель — измерить длину плоской базовой линии там, где встречаются крышка и часть под ней, касательные линии эллипса, где базовая линия пересекается с установленным эллипсом, и нарисовать эти линии, чтобы визуализировать качество измерений. Я написал код на Python, но не понимаю, почему две точки пересечения каким-то образом отклоняются от того места, где они на самом деле пересекаются, и касательные линии тоже отклоняются. Вот мой код и два видео. Спасибо!
Примечание: координата может быть обратной в случае, как в opencv2(x,y) и numpy(row, столбец)
ссылка на Dropbox для двух видео:
https://docsend.com/view/s/a33pt6nscx8t7z4e
текущий результат, не очень точный:
введите здесь описание изображения
Код:

Код: Выделить всё

import cv2
import numpy as np
import math
import pandas as pd
from skimage.measure import EllipseModel, ransac
import os
from scipy.optimize import fsolve

def ellipse_line_intersection(xc, yc, a, b, theta, vx, vy, x0, y0):
# Rotate line into ellipse's local coordinates
cos_theta = np.cos(-theta)
sin_theta = np.sin(-theta)

# Transform point and direction
x0_local = cos_theta * (x0 - xc) + sin_theta * (y0 - yc)
y0_local = -sin_theta * (x0 - xc) + cos_theta * (y0 - yc)
vx_local = cos_theta * vx + sin_theta * vy
vy_local = -sin_theta * vx + cos_theta * vy

# Quadratic coefficients
A = (vx_local / a)**2 + (vy_local / b)**2
B = 2 * (x0_local * vx_local / a**2 + y0_local * vy_local / b**2)
C = (x0_local / a)**2 + (y0_local / b)**2 - 1

# Solve quadratic equation
discriminant = B**2 - 4 * A * C
if discriminant <  0:
return []  # No intersection

t1 = (-B + np.sqrt(discriminant)) / (2 * A)
t2 = (-B - np.sqrt(discriminant)) / (2 * A)

# Intersection points in local coordinates
x1_local = x0_local + vx_local * t1
y1_local = y0_local + vy_local * t1
x2_local = x0_local + vx_local * t2
y2_local = y0_local + vy_local * t2

# Transform back to global coordinates
x1_global = xc + cos_theta * x1_local - sin_theta * y1_local
y1_global = yc + sin_theta * x1_local + cos_theta * y1_local
x2_global = xc + cos_theta * x2_local - sin_theta * y2_local
y2_global = yc + sin_theta * x2_local + cos_theta * y2_local

return [(x1_global, y1_global), (x2_global, y2_global)]

def calculate_tangent_slope(x, y, xc, yc, a, b, theta):  # x, y: intersection points

cos_theta = np.cos(-theta)
sin_theta = np.sin(-theta)
x_local = cos_theta * (x - xc) + sin_theta * (y - yc)
y_local = -sin_theta * (x - xc) + cos_theta * (y - yc)

# Gradient in local coordinates
dx = 2 * x_local / a**2
dy = 2 * y_local / b**2

# Transform gradient back to global coordinates
dx_global = cos_theta * dx - sin_theta * dy
dy_global = sin_theta * dx + cos_theta * dy

# Slope of the tangent line (perpendicular to the gradient)
if dx_global == 0:
return float('inf')  # Vertical line
return -dy_global / dx_global

# def ellipse_line_intersection(xc, yc, a, b, theta, vx, vy, x0, y0):  # fsolve to find two intersection points
#     # Parametric equations of the ellipse
#     def parametric_ellipse(t):
#         x = xc + a * np.cos(t) * np.cos(theta) + b * np.sin(t) * np.sin(theta)
#         y = yc - a * np.cos(t) * np.sin(theta) + b * np.sin(t) * np.cos(theta)
#         return x, y

#     # Function to solve for t (intersection condition)
#     def intersection_equation(t):
#         x, y = parametric_ellipse(t)
#         return (x - x0) * vy - (y - y0) * vx

#     # Solve for t values (two solutions)
#     t_values = fsolve(intersection_equation, [0, np.pi])  # Initial guesses pi
#     intersection_points = [parametric_ellipse(t) for t in t_values]
#     return intersection_points

# Load the video
video_path = "C:\\Users\\qiang.yu\\Desktop\\droplet2avi.avi"
cap = cv2.VideoCapture(video_path)
cap_beneath = cv2.VideoCapture("C:\\Users\\qiang.yu\\Desktop\\NW2avi.avi")
output_folder = "C:\\Users\\qiang.yu\\Desktop\\processed_frames7"
excel_path = "C:\\Users\\qiang.yu\\Desktop\\processed_frames7\\ellipse_baseline_data.xlsx"
os.makedirs(output_folder, exist_ok=True)

results = []
frame_count = 0
while cap.isOpened():
ret, frame = cap.read()
if ret is False or frame_count == 30:  # test for 30 frames only
break
frame_count += 1
kernel1 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (4, 4))

binary = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
_, binary = cv2.threshold(binary, 100, 255, cv2.THRESH_BINARY)
cleaned_binary1 = cv2.morphologyEx(binary, cv2.MORPH_OPEN, kernel1)  # cleaned binary is the smoothed ellipse cap
cleaned_binary = cv2.morphologyEx(cleaned_binary1, cv2.MORPH_CLOSE, kernel1)
cleaned_binary = cv2.morphologyEx(cleaned_binary, cv2.MORPH_DILATE, kernel1)
contours, _ = cv2.findContours(cleaned_binary, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)

ret, frame_beneath = cap_beneath.read()
binary_beneath = cv2.cvtColor(frame_beneath, cv2.COLOR_BGR2GRAY)
_, binary_beneath = cv2.threshold(binary_beneath, 100, 255, cv2.THRESH_BINARY)
# binary_beneath1 = cv2.morphologyEx(binary_beneath, cv2.MORPH_OPEN, kernel1)
cleaned_binary_beneath = cv2.morphologyEx(binary_beneath, cv2.MORPH_CLOSE, kernel1)
contours_beneath, _ = cv2.findContours(cleaned_binary_beneath, cv2.RETR_TREE, cv2.CHAIN_APPROX_SIMPLE)

if contours:
contours = sorted(contours, key=cv2.contourArea, reverse=True)
largest_contour = contours[0]
points = np.squeeze(largest_contour)

ellipse_model = EllipseModel()
if points.shape[0] >= 5:   # At least 5 points needed for ellipse fitting
model_robust, inliers = ransac(
points, EllipseModel, min_samples=5, residual_threshold=4, max_trials=3000
)

if model_robust:
# Extract ellipse parameters
xc, yc, a, b, theta = model_robust.params
center = (int(xc), int(yc))
axes = (int(a), int(b))
angle = np.degrees(theta)
cv2.ellipse(frame, center, axes, angle, 0, 360, (0, 255, 0), 1) # full (0,360deg) thickness 1 fitted ellipse in green 0,255,0
a_x1 = xc + a * np.cos(theta)
a_y1 = yc + a * np.sin(theta)
a_x2 = xc - a * np.cos(theta)
a_y2 = yc - a * np.sin(theta)

# Calculate b axis endpoints
b_x1 = xc + b * np.sin(theta)
b_y1 = yc - b * np.cos(theta)
b_x2 = xc - b * np.sin(theta)
b_y2 = yc + b * np.cos(theta)

if a > b: # draw major axis only
# Draw the a axis (in red),
cv2.line(frame, (int(a_x1), int(a_y1)), (int(a_x2), int(a_y2)), (0, 0, 255), 2)
else: # Draw the b axis (in blue)
cv2.line(frame, (int(b_x1), int(b_y1)), (int(b_x2), int(b_y2)), (255, 0, 0), 2)
# print(f"a: {a}, b: {b}, theta: {theta}")

kernel2 = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (8, 8))

cleaned_binary_dilated = cv2.morphologyEx(cleaned_binary, cv2.MORPH_DILATE, kernel2) #Dilated region
cleaned_binary_beneath_dilated = cv2.morphologyEx(cleaned_binary_beneath, cv2.MORPH_DILATE, kernel2)

intersection_like1 = cv2.bitwise_and(cleaned_binary_dilated, cleaned_binary_beneath_dilated)

red_mask = np.zeros_like(frame)  # Initialize an all-black frame
red_mask[cleaned_binary_beneath_dilated == 255] = (0, 0, 255)  #red color in NW
frame = cv2.addWeighted(frame, 1, red_mask, 1, 0)  # Add red mask with 0 transparency

blue_mask = np.zeros_like(frame)  # Initialize an all-black frame
blue_mask[intersection_like1 == 255] = (255, 0, 0)  # blue color in intersection region
frame = cv2.addWeighted(frame, 1, blue_mask, 1, 0)  # Add blue mask with 0 transparency

points_intersection = np.column_stack(np.where(intersection_like1 == 255))
points_intersection = points_intersection[:, [1, 0]]  # Switch to (x, y)
if points_intersection.shape[0] >  0:
# print(points_intersection.shape[0])
[vx, vy, x0, y0] = cv2.fitLine(points_intersection, cv2.DIST_L2, 0, 0.01, 0.01)
vx, vy, x0, y0 = vx.item(), vy.item(), x0.item(), y0.item()

left_point = (int(x0 - 1000), int(y0 - 1000 * vy/vx))
right_point = (int(x0 + 1000), int(y0 + 1000 * vy/vx))
# Draw the fitted line (baseline) on the frame in green
cv2.line(frame, left_point, right_point, (0, 255, 0), 2)

intersections = ellipse_line_intersection(xc, yc, a, b, theta, vx, vy, x0, y0)
x1, y1, x2, y2 = intersections[0][0], intersections[0][1], intersections[1][0], intersections[1][1]
slope1 = calculate_tangent_slope(x1, y1, xc, yc, a, b, theta)
if slope1 == float('inf'):  # Vertical tangent line
start_point = (int(x1), int(y1 - 1000))
end_point = (int(x1), int(y1 + 1000))
else:
start_point = (int(x1 - 1000), int(y1 - 1000 * slope1))
end_point = (int(x1 + 1000), int(y1 + 1000 * slope1))
cv2.line(frame, start_point, end_point, (255, 255, 0), 2)  # Tangent in cyan

slope2 = calculate_tangent_slope(x2, y2, xc, yc, a, b, theta)
if slope1 == float('inf'):  # Vertical tangent line
start_point = (int(x2), int(y2 - 1000))
end_point = (int(x2), int(y2 + 1000))
else:
start_point = (int(x2 - 1000), int(y2 - 1000 * slope2))
end_point = (int(x2 + 1000), int(y2 + 1000 * slope2))
# Draw the tangent line
cv2.line(frame, start_point, end_point, (255, 255, 0), 2)  # Tangent in cyan
cv2.circle(frame, (int(x1), int(y1)), 5, (0, 255, 255), -1)  # Draw intersection points in yellow
cv2.circle(frame, (int(x2), int(y2)), 5, (0, 255, 255), -1)  # Draw intersection points in yellow

contact_angle1 = math.degrees(math.atan(abs(slope1)) - math.atan(abs(vy/vx)))
contact_angle2 = math.degrees(math.atan(abs(slope2)) - math.atan(abs(vy/vx)))
contact_angle = (contact_angle1 + contact_angle2)/2
print("Contact angle:", contact_angle)

cv2.circle(frame, (int((intersections[0][0]+intersections[1][0])/2), int((intersections[0][1]+intersections[1][1])/2)), 5, (0, 0, 255), -1)  # center point in red
baseline_length = (math.sqrt((intersections[0][0]-intersections[1][0])**2+(intersections[0][1]-intersections[1][1])**2))

cv2.imshow("Red Highlighted Difference", frame)

# Save results
results.append({
"Frame": frame_count,
"Major Axis": a,
"Minor Axis": b,
"Baseline Length": baseline_length,
"Contact angle": 180-contact_angle
})

if cv2.waitKey(1) & 0xFF == ord('q'):
break
output_image_path = os.path.join(output_folder, f"frame_{frame_count:04d}.png")
cv2.imwrite(output_image_path, frame)

# Save measured data to Excel
df = pd.DataFrame(results)
df.to_excel(excel_path, index=False)
print(f"Results saved to {excel_path}")

# Release resources
cap.release()
cv2.destroyAllWindows()

измерьте длину плоской базовой линии там, где встречаются крышка и часть под ней, касательные линии эллипса, где базовая линия пересекается с установленным эллипсом, и нарисуйте эти линии, чтобы визуализировать качество измерений
п>

Подробнее здесь: https://stackoverflow.com/questions/792 ... ent-by-cv2