Ctrl+K

Orientation an aspect ratio

Orientation an aspect ratio#

In this example, we show how to use inertia matrix of a given labeled object to find its orientation.

Required files

Before using this notebook, download the image blobs.jpg

%matplotlib widget
from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm
from scipy import ndimage
import pandas as pd
import os

path = "blobs.jpg"
files = os.listdir("./")
if path in files:
    print("Ok, the file is in {0}".format(files))
else:
    print("The file is not in {0} , retry !".format(files))

Ok, the file is in ['.ipynb_checkpoints', 'image_processing_basics.ipynb', 'blobs.jpg', 'blob_orientation.ipynb', 'rabbit.jpg', 'support_material.md']

In order to explain the concepts of inertia and aspect ratio, we use this magnificient hand drawing:

im = Image.open(path)
Nc, Nl = im.size
im = im.resize((Nc // 4, Nl // 4))
fig, ax = plt.subplots()
# ax.axis("off")
plt.imshow(im)
plt.colorbar()
plt.show()

R, G, B = im.split()
R = np.array(R)
G = np.array(G)
B = np.array(B)

plt.figure()
plt.hist(
    R.flatten(),
    bins=np.arange(256),
    histtype="stepfilled",
    color="r",
    alpha=0.3,
    label="Red",
)
plt.hist(
    G.flatten(),
    bins=np.arange(256),
    histtype="stepfilled",
    color="g",
    alpha=0.3,
    label="Green",
)
plt.hist(
    B.flatten(),
    bins=np.arange(256),
    histtype="stepfilled",
    color="b",
    alpha=0.3,
    label="Blue",
)
plt.grid()
plt.legend()
plt.show()

Thresholding level is obvious:

Bt = np.where(B < 50, 1, 0)
plt.figure()
plt.imshow(Bt, cmap=cm.gray)
plt.colorbar()
plt.show()

Btc = ndimage.binary_closing(Bt, structure=np.ones((5, 5)))


Bl, number = ndimage.label(Btc)
plt.figure()
plt.imshow(np.where(Bl != 0, Bl, np.nan), cmap=cm.jet)
plt.show()
number

10

obj = ndimage.find_objects(Bl)
len(obj)

10

Inertia matrix of an object#

The object represented bellow is stretched in a direction. Let’s see how we can use its inertia matrix to determine in which direction and how much it is stretched.

plt.figure()
plt.imshow(np.array(im)[obj[1]])
plt.show()

The inertia matrix of a 2D object can be defined as follows:

\[\begin{split} I = \begin{bmatrix} I_{xx} & -I_{xy} \\ -I_{xy} & I_{yy} \\ \end{bmatrix} \end{split}\]

This matrix is symmetric and, as a consequence, it can be diagonalized in a new frame rotated by an angle \(\theta\) in the plane. This frame is composed of the two normalized eigenvectors \((\vec e_1, \vec e_2)\) of the matrix. In this frame, the matrix has two eigenvalues \((I_1, I_2)\) ordered so that \(I_1 \geq I_2\). Then:

  • \(\theta = (\vec x, \vec e_1)\) and,

  • The aspect ratio \(a = \sqrt{I_1 / I_2}\).

The angle \(\theta\) gives the direction of the elongation of the object and \(a\) shows how much it is elongated. For example, if \(a == 1\), the object is not elongated whereas if \(a=10\) it is 10 times longer in direction 1 than in direction 2 in an inertial point of view.

data = pd.DataFrame(
    columns=["area", "xg", "yg", "Ixx", "Iyy", "Ixy", "I1", "I2", "theta"]
)
for i in range(len(obj)):
    x, y = np.where(Bl == i + 1)
    xg, yg = x.mean(), y.mean()
    x = x - xg
    y = y - yg
    A = len(x)
    Ixx = (y**2).sum()
    Iyy = (x**2).sum()
    Ixy = (x * y).sum()
    I = np.array([[Ixx, -Ixy], [-Ixy, Iyy]])
    eigvals, eigvecs = np.linalg.eig(I)
    eigvals = abs(eigvals)
    loc = np.argsort(eigvals)[::-1]
    d = eigvecs[loc[0]]
    d *= np.sign(d[0])
    theta = np.degrees(np.arccos(d[1]))
    eigvals = eigvals[loc]
    data.loc[i] = [A, xg, yg, Ixx, Iyy, Ixy, eigvals[0], eigvals[1], theta]
data.sort_values("area", inplace=True, ascending=False)
data["aspect_ratio"] = (data.I1 / data.I2) ** 0.5

data[["area", "theta", "aspect_ratio"]]
area theta aspect_ratio
1 2970.0 116.026308 3.969835
2 1891.0 96.371877 5.339632
4 442.0 121.962612 8.767714
6 327.0 3.551442 1.442622
5 296.0 87.196715 1.879190
8 246.0 57.994151 1.817390
0 153.0 91.013709 1.439440
7 150.0 100.404572 1.656024
9 94.0 82.012160 1.342536
3 90.0 97.955535 1.375880 
fig = plt.figure()
counter = 1
for i in data.index.values:
    ax = fig.add_subplot(3, 4, counter)
    z = Image.fromarray(np.array(im)[obj[i]])
    z = z.rotate(-data.loc[i, "theta"] + 90, expand=True)
    z = np.array(z)
    plt.imshow(z)
    plt.title(str(i))
    ax.axis("off")
    counter += 1
    # plt.grid()
plt.show()
Contents