Projecting ellipses onto circles using homographies

Question

I am relatively new to computer vision, and I am experimenting with extracting a clock/dial face from an image.

I managed to fit an ellipse onto the dial, and I wanted to fix the perspective such that the dial would face the camera directly. Essentially, I wanted to project an ellipse onto a circle.

I found online that homographies are normally used to accomplish this task.

I tried following the OpenCV tutorial on homographies, but I am having problems with the source poorly matching the destination.

Reading similar questions, a perfect projection seems impossible since perspective-projected circles are not really ellipses. However, I am not looking for a mathematically exact projection, but even the simplest cases I have seem to produce bad results (see example 4 in the images below).

Here are four example input images on which I annotated the ellipses (red), unit circles (blue), source points (yellow), and destination points (cyan) (note: the centers are annotated but are not used to compute the homography). Annotated input image

However, after applying the homography computed using OpenCV (I got the same exact results using Scikit), the ellipses are not projected well onto the unit circles. Homography results

I also tried using more than 4 points to compute the homography, but again, the results were identical.

Here is the code I used:

def persp_transform(orig):
    img = orig.copy()

    ellipses = find_ellipses(img)
    ellipse = ellipses[0] # Use the largest one

    center_x = ellipse[0]
    center_y = ellipse[1]
    center = np.array([center_x, center_y])
    width = ellipse[2]
    height = ellipse[3]
    angle = ellipse[4]
    minor_semiaxis = min(width, height)/2
    major_semiaxis = max(width, height)/2

    # Translate the image to center the ellipse
    img_center = np.array([img.shape[0]//2, img.shape[1]//2])
    translation = center - img_center
    M = np.float32([[1, 0, -translation[1]], [0, 1, -translation[0]]])
    img = cv.warpAffine(img, M, (img.shape[1], img.shape[0]))

    # Draw ellipse before projection
    rr, cc = ellipse_perimeter(img_center[0], img_center[1], int(major_semiaxis), int(minor_semiaxis), orientation=angle, shape=img.shape)
    draw_ellipse(img, rr, cc, color=0, thickness=10)
    
    def sin(a): return np.sin(np.deg2rad(a))
    def cos(a): return np.cos(np.deg2rad(a))
    
    # Source points around the ellipse
    num_points = 5
    rotation = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]])
    from_points = np.array([[major_semiaxis*sin(a), minor_semiaxis*cos(a)] for a in np.linspace(0, 360, num_points)]) @ rotation + img_center
    from_points = from_points.astype(np.float32)

    # Destination points around a centered circle
    to_radius = int((min(img.shape[:2]) / 2) * 0.8)
    to_points = np.array([[to_radius*sin(a), to_radius*cos(a)] for a in np.linspace(0, 360, num_points)]) @ rotation + img_center
    to_points = to_points.astype(np.float32)

    # Draw ellipse center and source points before projection
    cv.circle(img, (img_center[1], img_center[0]), 20, (255, 255, 0), -1)
    for fp in from_points:
        cv.circle(img, (int(fp[1]), int(fp[0])), 20, (255, 255, 0), -1)
    
    # Compute homography and project
    M, _ = cv.findHomography(from_points, to_points, cv.RANSAC, 5.0)
    img = cv.warpPerspective(img, M, (img.shape[1], img.shape[0]))

    # Draw target unit circle and destination points after projection
    cv.circle(img, (img_center[1], img_center[0]), to_radius, (0, 0, 255), 20)
    cv.circle(img, (img_center[1], img_center[0]), 20, (0, 255, 255), -1)
    for tp in to_points:
        cv.circle(img, (int(tp[1]), int(tp[0])), 20, (0, 255, 255), -1)

    return img

EDIT 1

Here are the 4 raw example images I used as input: GDrive Link

To fit the ellipses, I use the method proposed in AAMED, which can be found on GH here.

I manually played with the input parameters until I found satisfactory ones since I didn't need to find all the ellipses in the image, only the main ones.

My ellipse finding code looks like this:

def find_ellipses(orig, scale=0.3, theta_fsa=20, length_fsa=5.4, T_val=0.9):
    img = cv.resize(orig.copy(), None, fx=scale, fy=scale)
    gray_image = cv.cvtColor(img, cv.COLOR_BGR2GRAY)

    aamed = AAMED(img.shape[0]+1, img.shape[1]+1)
    aamed.setParameters(np.deg2rad(theta_fsa), length_fsa, T_val)
    ellipses = aamed.run_AAMED(gray_image)
    aamed.release()

    ellipses = sorted(ellipses, key=lambda e: ellipse_area(e), reverse=True)
    ellipses = [np.array(e) / scale for e in ellipses]

    return ellipses

Then I use skimage.draw.ellipse_perimeter(...) to impose the largest ellipse on the image.

Edit 2

I tried to better visualize the homography and narrowed down the issue possibly to cv.warpPerspective(...).

Firstly, I color-coded each point pair to check that I was mapping the right pairs: Color-coded points

Then, I computed the homography and individually projected each source point using cv.perspectiveTransform(...) without warping the image to see where they would end up, and they seemed to be projected properly: Projected points

However, when I then use cv.warpPerspective(...) to project the image, the projection is inconsistent with the result from cv.perspectiveTransform(...): Wrong projection

Answer 1

Here is one attempt in Python/OpenCV. The idea is the get the ellipse. Then get the 4 points on the ellipse that corresponds to the min and max semi-major radii. Then get the 4 corresponding points on a circle or radius equal to the max semi-major radii. Then compute the homography for the two sets of 4 corresponding points. Finally warp the image.

• Read the input
• Read the ellipse mask
• Get the contour of the mask
• Fit and ellipse to the contour
• Compute the 4 points on the ellipse from the ellipse parameters
• Compute the 4 points on the circle using the max semi-major radius
• Get the homography between the two sets of points
• Warp the input image
• Save the results

Input:

Ellipse Mask (fit manually in Photoshop):

import cv2
import numpy as np
import math

# read the input
img = cv2.imread('ex_4.jpg')

# read ellipse mask as grayscale
mask = cv2.imread('ex_4_mask.png', cv2.IMREAD_GRAYSCALE)

# get ellipse contour
contours = cv2.findContours(mask, cv2.RETR_EXTERNAL, cv2.CHAIN_APPROX_SIMPLE)
contours = contours[0] if len(contours) == 2 else contours[1]
big_contour = max(contours, key=cv2.contourArea)

# draw contour on black background
mask2 = np.zeros_like(mask, dtype=np.uint8)
cv2.drawContours(mask2, [big_contour], 0, 255, 1)

# fit ellipse to mask contour
points = np.argwhere(mask2.transpose()>0)
ellipse = cv2.fitEllipse(points)
(cx,cy), (d1,d2), angle = ellipse

# draw ellipse
ellipse = img.copy()
cv2.ellipse(ellipse, (int(cx),int(cy)), (int(d1/2),int(d2/2)), angle, 0, 360, (0,0,255), 4)

# correct angle
if angle > 90:
    angle = angle - 90
else:
    angle = angle + 90
print('center: (', cx,cy, ')', 'diameters: (', d1, d2, ')', 'angle:', angle)

# get 4 points from min and max semi-major radii
r1 = min(d1,d2)/2
r2 = max(d1,d2)/2

x1 = np.intp( cx + math.cos(math.radians(angle))*r2 )
y1 = np.intp( cy + math.sin(math.radians(angle))*r2 )
x2 = np.intp( cx + math.cos(math.radians(angle+180))*r2 )
y2 = np.intp( cy + math.sin(math.radians(angle+180))*r2 )
x3 = np.intp( cx + math.cos(math.radians(angle+90))*r1 )
y3 = np.intp( cy + math.sin(math.radians(angle+90))*r1 )
x4 = np.intp( cx + math.cos(math.radians(angle+270))*r1 )
y4 = np.intp( cy + math.sin(math.radians(angle+270))*r1 )

print('x1,y1:', x1, y1)
print('x2,y2:', x1, y1)
print('x3,y3:', x1, y1)
print('x4,y4:', x1, y1)

# input points from ellipse
input_pts = np.float32([[x1,y1], [x2,y2], [x3,y3], [x4,y4]])

# output points from circle of radius = max semi-major radii
ox1 = np.intp( cx + math.cos(math.radians(angle))*r2 )
oy1 = np.intp( cy + math.sin(math.radians(angle))*r2 )
ox2 = np.intp( cx + math.cos(math.radians(angle+180))*r2 )
oy2 = np.intp( cy + math.sin(math.radians(angle+180))*r2 )
ox3 = np.intp( cx + math.cos(math.radians(angle+90))*r2 )
oy3 = np.intp( cy + math.sin(math.radians(angle+90))*r2 )
ox4 = np.intp( cx + math.cos(math.radians(angle+270))*r2 )
oy4 = np.intp( cy + math.sin(math.radians(angle+270))*r2 )

cx = np.intp(cx)
cy = np.intp(cy)
r2 = np.intp(r2)
ellipse_pts = ellipse.copy()

# draw white circle on copy of ellipse
cv2.circle(ellipse_pts, (cx,cy), r2, (255,255,255), 4)

output_pts = np.float32([[ox1,oy1], [ox2,oy2], [ox3,oy3], [ox4,oy4]])

# draw output points on copy of ellipse image
cv2.circle(ellipse_pts, (ox1,oy1), 16, (0,255,0), 4)
cv2.circle(ellipse_pts, (ox2,oy2), 16, (0,255,0), 4)
cv2.circle(ellipse_pts, (ox3,oy3), 16, (0,255,0), 4)
cv2.circle(ellipse_pts, (ox4,oy4), 16, (0,255,0), 4)


# draw input points on copy of ellipse image
cv2.circle(ellipse_pts, (x1,y1), 12, (255,0,0), -1)
cv2.circle(ellipse_pts, (x2,y2), 12, (255,0,0), -1)
cv2.circle(ellipse_pts, (x3,y3), 12, (255,0,0), -1)
cv2.circle(ellipse_pts, (x4,y4), 12, (255,0,0), -1)

# get homography when only 4 pts
h = cv2.getPerspectiveTransform(input_pts, output_pts)

# warp image
result = cv2.warpPerspective(img, h, (img.shape[1],img.shape[0]))


# save results (compress to fit 2Mb limit on this forum)
cv2.imwrite('ex_4_ellipse.jpg', ellipse, [int(cv2.IMWRITE_JPEG_QUALITY), 85])
cv2.imwrite('ex_4_ellipse_pts.jpg', ellipse_pts, [int(cv2.IMWRITE_JPEG_QUALITY), 85])
cv2.imwrite('ex_4_ellipse2circle.jpg', result, [int(cv2.IMWRITE_JPEG_QUALITY), 85])

# show results
cv2.imshow('ellipse', ellipse)
cv2.imshow('ellipse_pts', ellipse_pts)
cv2.imshow('result', result)
cv2.waitKey(0)

Ellipse on input:

Ellipse with points and circle:

Warped Result:

Answer 2

Thanks to Christoph Rackwitz and fmw42, I tracked down the issue to improperly formatted coordinates. Be it swapped x and y positions or badly ordered points in arrays. I attribute the cause to inconsistencies in the order of the coordinates required by OpenCV and PLT, which can be confusing.

Regardless, not using PLT for drawing on the image and relying completely on OpenCV made the final code much clearer, and the issue is not present in the final rewrite.

Another imprecision I noticed is that the rotation matrix computed explicitly using Numpy gave worse results than the one obtained using cv.getRotationMatrix2D(...), which further improved the final result.

I guess if there needs to be a final lesson, it should be that even when prototyping, writing clear and unambiguous code is important.