find angle between major axis of ellipse and x-axis of coordinate (help me implement method from paper)

Question

So I am trying to implement a method from this paper. I am stuck at the part where I have to find the angle between the major axis of the lesion’s best-fit ellipse and the x-axis of the coordinate system.

Here is the sample image:

Here is what I got so far:

Is it possible to find that angle? And after the angle has been found, I have to flip the RoI along x-axis by the angle.

UPDATE ----------

Google drive link to Roi Image: RoI image

Implementing method step by step based on the paper.

First, I should recenter the RoI to the center of the image coordinate. In the paper, they centered the RoI using its centroid. I manage to do it based on this code I found in this answer. The result is fine if my RoI is small and not touching the image border. But if I have large image the result is really bad. So I ended up centering the RoI using boundingRect. Here is the result of centering:

Code for centering RoI:

import math

import cv2
import numpy as np
import matplotlib.pyplot as plt

# read image
cont_img = cv2.imread(r"C:\Users\Pandu\Desktop\IMD064_lesion.bmp", 0)
cont_rgb = cv2.cvtColor(cont_img, cv2.COLOR_GRAY2RGB)

# fit ellipse and find ellipse properties
hh, ww = cont_img.shape
contours, hierarchy = cv2.findContours(cont_img, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
ellipse = cv2.fitEllipse(contours[0])
(xc, yc), (d1, d2), angle = ellipse

# centering by centroid
half_width = int(ww/2)
half_height = int(hh/2)
offset_x = (half_width-xc)
offset_y = (half_height-yc)

T = np.float32([[1, 0, offset_x], [0, 1, offset_y]])
centered_by_centroid = cv2.warpAffine(cont_img.copy(), T, (ww, hh))
plt.imshow(centered_by_centroid, cmap=plt.cm.gray)

# centering by boundingRect
# This centered RoI is (L)
x, y, w, h = cv2.boundingRect(contours[0])
startx = (ww - w)//2
starty = (hh - h)//2
centered_by_boundingRect = np.zeros_like(cont_img)
centered_by_boundingRect[starty:starty+h, startx:startx+w] = cont_img[y:y+h, x:x+w]
plt.imshow(centered_by_boundingRect, cmap=plt.cm.gray)

Second, after centering the RoI, I should find the orientation angel and rotate the RoI based on that angel and then flip . Using code from this answer. (is this the correct way to rotate the RoI?):

# find ellipse properties of centered RoI
contours, hierarchy = cv2.findContours(centered_by_boundingRect, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)
ellipse = cv2.fitEllipse(contours[0])
(xc, yc), (d1, d2), angle = ellipse
roi_centroid = (xc, yc)

rot_angle = 90 - angle
if rot_angle < 0:
    rot_angle += 180

# This rotated RoI is (Lx)
M = cv2.getRotationMatrix2D(roi_centroid, -rot_angle, 1.0)
rot_im = cv2.warpAffine(centered_by_boundingRect, M, (ww, hh))
plt.imshow(rot_im, cmap=plt.cm.gray)

# (Ly)
# by passing 0 to flip() should flip image around x-axis, but I get the same result as the paper
res_flip_y = cv2.flip(rot_im.copy(), 0)
plt.imshow(res_flip_y , cmap=plt.cm.gray)

# (L) (xor) (Lx)
res_x_xor = cv2.bitwise_xor(centered_by_boundingRect, rot_im)
plt.imshow(res_x_xor, cmap=plt.cm.gray)

# (L) (xor) (Ly)
res_y_xor = cv2.bitwise_xor(centered_by_boundingRect, res_flip_x)
plt.imshow(res_y_xor, cmap=plt.cm.gray)

I still can't get the same result as the paper, the rotating operation also produce bad result on large RoI. Help...

UPDATE ---------- 20/03/2021

Small RoI: fine result on rotation and looks similar with the paper, but still not getting the same end result on the L (xor) Lx or L (xor) Ly

Large RoI: bad result on rotation as the RoI get out of border/image

Answer 1

The angle you're looking for is returned from fitEllipse. It's just rotated a bit according to a different reference frame. You can get your counter-clockwise rotation angle by doing 90 - angle. As for rotating the roi you can either use minAreaRect to get a minimum-fit rectangle directly, or you can fit a bounding box to the contour and rotate each point individually.

The green rectangle is the minAreaRect(), the red rectangle is the boundingRect() after it's been rotated.

import cv2
import numpy as np
import math

# rotate point
def rotate2D(point, deg):
    rads = math.radians(deg);
    x, y = point;
    rcos = math.cos(rads);
    rsin = math.sin(rads);
    rx = x * rcos - y * rsin;
    ry = x * rsin + y * rcos;
    rx = round(rx);
    ry = round(ry);
    point[0] = rx;
    point[1] = ry;

# translate point
def translate2D(src, target, sign):
    tx, ty = target;
    src[0] += tx * sign;
    src[1] += ty * sign;

# read image
cont_img = cv2.imread("blob.png", 0)
cont_rgb = cv2.cvtColor(cont_img, cv2.COLOR_GRAY2RGB)

# find contour
_, contours, hierarchy = cv2.findContours(cont_img, cv2.RETR_LIST, cv2.CHAIN_APPROX_SIMPLE)

# fit ellipse and get ellipse properties
ellipse = cv2.fitEllipse(contours[0])
(xc, yc), (d1, d2), angle = ellipse

# -------- NEW STUFF IN HERE --------------
# calculate counter-clockwise angle relative to x-axis
rot_angle = 90 - angle;
if rot_angle < 0:
    rot_angle += 180;
print(rot_angle);

# if you want a rotated ROI I would recommend using minAreaRect rather than rotating a different rectangle
# fit a minrect to the image # this is taken directly from OpenCV's tutorials
rect = cv2.minAreaRect(contours[0]);
box = cv2.boxPoints(rect);
box = np.int0(box);
cv2.drawContours(cont_rgb, [box], 0, (0,255,0), 2);

# but if you really want to use a different rectangle and rotate it, here's how to do it
# create rectangle
x,y,w,h = cv2.boundingRect(contours[0]);
rect = [];
rect.append([x,y]);
rect.append([x+w,y]);
rect.append([x+w,y+h]);
rect.append([x,y+h]);

# rotate it
rotated_rect = [];
center = [x + w/2, y + h/2];
for point in rect:
    # for each point, center -> rotate -> uncenter
    translate2D(point, center, -1);
    rotate2D(point, 90 - rot_angle); # "90 - angle" is because rotation goes clockwise
    translate2D(point, center, 1);
    rotated_rect.append([point]);
rotated_rect = np.array(rotated_rect);
cv2.drawContours(cont_rgb, [rotated_rect.astype(int)], -1, (0,0,255), 2);
# ------------- END OF NEW STUFF -----------------


# draw fitted ellipse and centroid
target_ellipse = cv2.ellipse(cont_rgb.copy(), ellipse, (37, 99, 235), 10)
centroid = cv2.circle(target_ellipse.copy(), (int(xc), int(yc)), 20, (250, 204, 21), -1)

# draw major axis
rmajor = max(d1, d2)/2

if angle > 90:
    angle = angle - 90
else:
    angle = angle + 90

xtop_major = xc + math.cos(math.radians(angle))*rmajor
ytop_major = yc + math.sin(math.radians(angle))*rmajor
xbot_major = xc + math.cos(math.radians(angle+180))*rmajor
ybot_major = yc + math.sin(math.radians(angle+180))*rmajor
top_major = (int(xtop_major), int(ytop_major))
bot_major = (int(xbot_major), int(ybot_major))

target_major_axis = cv2.line(centroid.copy(),
                             top_major, bot_major,
                             (0, 255, 255), 5)

## image center coordinate
hh, ww = target_major_axis.shape[:2];
x_center_start = (0, int(hh/2))
x_center_end = (int(ww), int(hh/2))
y_center_start = (int(ww/2), 0)
y_center_end = (int(ww/2), int(hh))
img_x_middle_coor = cv2.line(target_major_axis.copy(), x_center_start, x_center_end, (219, 39, 119), 10)
img_y_middle_coor = cv2.line(img_x_middle_coor.copy(), y_center_start,
                             y_center_end, (190, 242, 100), 10)

# show
cv2.imshow("image", img_y_middle_coor);
cv2.waitKey(0);

For the future: check that your code runs before pasting it on here. Aside from the missing "import" lines it was also missing this line:

hh, ww = target_major_axis.shape[:2]

If the sample code you paste has errors, then everyone who wants to help will have to waste some time bug-stomping before they can begin working on a solution.