fitting ellipse in images with poor contrast

Question

I am working on image processing in Python, on the topic of underwater photogrammetry. My goal is to fit an ellipse to fidicual markers, and retrieve its a) center, b) axis and c) orientation.

My markers are

1. radial,
2. white on black background, and some have a
3. binary code:
   

A ML-model delivers a small image snippets for each marker in each image, containting only the center of the marker.

So far, I've implemented these approaches:

1. Using openCV: a) Thresholding, which results in a binary image (cv2.threshold) b) Find Contours (cv2.findContours) c) fit Ellipse (v2.fitEllipse)
2. Using Scikit: a) Detect edge (using Canny) b) Apply hough transform
3. Star operator (work in progress) a) Estimate ellipse center b) Send 360rays in all directions c) Build an array, comprising coordinates of the largest gradient on each ray d) Calculate best-fit ellipse using least-square method e) Use the new center to repeat process (possibly several iterations required)

I perform these methods for each color-channel seperately. So far, the results between channels differ within several pixels for the ellipse center.

1. Do you have any suggestions on what pre-processing methods I should use, prior detecting/fitting the ellipse?
2. Any thoughts on which of the above methods will lead to the most accurate results?

Answer 1

This is amazing! Thank you. I just started to read about moments (e.g. https://www.pythonpool.com/opencv-moments/) and inertia.

However, there is a challange applying your code to this example: As you can see, the image was poorly cropped, and the inertia of the image is more in the image center than in the center of the expected ellipse.

My first attempt to fix this is to binarize the image first:

import cv2 as cv2
T = int(cv2.mean(image)[0])
ret,image = cv2.threshold(image,T,255,0)

Is that a reasonable approach? I fear, that the binarization will have an unwanted impact on the moments of inertia. Thank you for claryfying.

Answer 2

This code finds the center of mass of the image, and the main axis of symmetry by calculating the moments of inertia.

I tried many libraries that calculate moments of inertia of images, but they give strange results (like 4x4 matrix for what should be a 2x2 matrix of inertia. Also, ndimage.measurements.center_of_mass() appears to return (Cy,Cx) (row, column)

So, I resorted to manually calculating the moments of inertia

import matplotlib.pyplot as plt
import numpy as np
from PIL import Image as Pim
from io import BytesIO
import requests
photoURL = "https://i.stack.imgur.com/EcLYk.png"

response = requests.get(photoURL)

image = np.array(Pim.open(BytesIO(response.content)).convert('L'))  # Convert to greyscale

plt.imshow(image)

if True:  # calculate eigen vectors = main axis of inertia
    # xCoord, yCoord are the column and row numbers in image
    xCoord, yCoord = np.meshgrid(np.arange(image.shape[1]), np.arange(image.shape[0]))
    # mass M is the total sum of Image
    M = np.sum(image)
    # Cx, Cy are the coordinates of the center of mass
    #Cx = sum(xCoord * image) / sum(image)
    Cx = np.einsum('ij,ij', xCoord, image)/M
    Cy = np.einsum('ij,ij', yCoord, image)/M
    # Ixx is the second order moment of image respect to the horizontal axis passing through the center of mass
    # Ixx=sum(Image*y^2)
    Ixx = np.einsum('ij,ij,ij', yCoord-Cy, yCoord-Cy, image)
    # Iyy is the second order moment of image respect to the vertical axis passing through the center of mass
    # Iyy=sum(Image*x^2)
    Iyy = np.einsum('ij,ij,ij', xCoord-Cx, xCoord-Cx, image)
    # Ixy is the second order moment of image respect to both axis passing through the center of mass
    # Ixy=sum(Image*x*y)
    Ixy = np.einsum('ij,ij,ij', xCoord-Cx, yCoord-Cy, image)

    inertiaMatrix = np.array([[Ixx, Ixy],
                              [Ixy, Iyy]])
    eigValues, eigVectors = np.linalg.eig(inertiaMatrix)

# Plot center of mass
plt.scatter(Cx, Cy, c='r')
# Plot eigenvectors from center to direction of eigenvectors
plt.quiver(Cx, Cy, eigVectors[0, 0], eigVectors[1, 0], color='r', scale=2)
plt.quiver(Cx, Cy, eigVectors[0, 1], eigVectors[1, 1], color='r', scale=2)

plt.show()

nothing = 0