计算图像的形心主惯性轴

import numpy as np
import math
import numba as nb
from matplotlib import pyplot as plt

# 要求是二值图，前景为1，背景为0
@nb.jit(nb.types.UniTuple(nb.float32, 2)(nb.uint8[:, :]), nopython=True)
def gravity_center(array):
    """计算重心，array的元素都是1或者0,计算1部分的重心"""
    rows, cols = array.shape
    area = 0
    sum_x = 0
    sum_y = 0
    for x in range(rows):
        for y in range(cols):
            if array[x, y] != 0:
                area += 1
                sum_x += x
                sum_y += y
    return sum_x / area, sum_y / area  # x 对应行 ， y 对于列！

# 取图像的竖直方向为x方向

@nb.jit(nb.float32(nb.uint8[:, :], nb.float32), nopython=True)
def calc_Ix(array, axis=0.0):  # moments of inertia
    # 计算 对于 y= axis （某列）的惯性矩
    rows, cols = array.shape
    sum = 0
    for x in range(rows):
        for y in range(cols):
            if array[x, y] != 0:
                sum += (y - axis) ** 2
    return sum


@nb.jit(nb.float32(nb.uint8[:, :], nb.float32), nopython=True)
def calc_Iy(array, axis=0.0):  # moments of inertia
    # 计算 对于 x= axis （某行）的惯性矩
    rows, cols = array.shape
    sum = 0
    for x in range(rows):
        for y in range(cols):
            if array[x, y] != 0:
                sum += (x - axis) ** 2
    return sum


@nb.jit(nb.float32(nb.uint8[:, :], nb.float32, nb.float32), nopython=True)
def calc_Ixy(array, axis_x=0.0, axis_y=0.0):  # product of inertia
    # 计算 对于 某行，某列的惯性积
    rows, cols = array.shape
    sum = 0
    for x in range(rows):
        for y in range(cols):
            if array[x, y] != 0:
                sum += (x - axis_x) * (y - axis_y)
    return sum

@nb.jit(nb.types.UniTuple(nb.float32, 3)(nb.float32, nb.float32, nb.float32), nopython=True)
def calc_I1_I2_alpha(Ix, Iy, Ixy):
    temp1 = 0.5 * (Ix + Iy)
    temp2 = 0.5 * math.sqrt((Ix - Iy) ** 2 + 4.0 * Ixy ** 2)
    I1 = temp1 + temp2
    I2 = temp1 - temp2
    if Ix == Iy:
        alpha = math.pi / 2.0
    else:
        alpha = 0.5 * math.atan(-2.0 * Ixy / (Ix - Iy))  # 值域 （-pi/2，+pi/2）
    return I1, I2, alpha

if __name__ == "__main__":
    img2 = np.ones((1000,1200),dtype=np.uint8)

    cx,cy = gravity_center(img2)
    Ix = calc_Ix(img2, axis=cy)
    Iy = calc_Iy(img2, axis=cx)
    Ixy = calc_Ixy(img2, axis_x=cx, axis_y=cy)
    I1, I2, alpha = calc_I1_I2_alpha(Ix, Iy, Ixy)
    print(f"Ix = {Ix}, Iy = {Iy}, Ixy = {Ixy}")
    print(f"I1 = {I1}, I2 = {I2}, alpha = {alpha}, angle ={alpha/math.pi *180} degree")

    plt.imshow(img2, cmap ="gray")
    plt.show()

Ix = 143999893504.0, Iy = 99999899648.0, Ixy = 0.0
I1 = 143999893504.0, I2 = 99999899648.0, alpha = -0.0, angle =-0.0 degree

if __name__ == "__main__":
    from PIL import Image
    img = Image.open("cat3.png")
    img = np.array(img)
    img = img[:,:,0]
    img2 = np.zeros(img.shape, dtype = np.uint8)  # reversed
    img2[np.where(img<150)] = 1

    # img2 = np.ones((1000,1200),dtype=np.uint8)
    cx,cy = gravity_center(img2)
    Ix = calc_Ix(img2, axis=cy)
    Iy = calc_Iy(img2, axis=cx)
    Ixy = calc_Ixy(img2, axis_x=cx, axis_y=cy)
    I1, I2, alpha = calc_I1_I2_alpha(Ix, Iy, Ixy)
    print(f"Ix = {Ix}, Iy = {Iy}, Ixy = {Ixy}")
    print(f"I1 = {I1}, I2 = {I2}, alpha = {alpha}, angle ={alpha/math.pi *180} degree")

    plt.imshow(img2, cmap ="gray")
    plt.show()

Ix = 21522734.0, Iy = 21877528.0, Ixy = 5528713.5
I1 = 27231690.0, I2 = 16168573.0, alpha = 0.7693604230880737, angle =44.081105167346 degree