高级示例：方向和纵横比

在这个例子中，我们展示了如何使用一个给定标签对象的惯量矩阵来找到它的方向。

图像0

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
%matplotlib nbagg

path = "balls.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 ['Image_Processing_Tutorial_3.ipynb', 'Orientation.ipynb', '.ipynb_checkpoints', 'Image_Processing_Tutorial_4.ipynb', 'bugs.jpg', 'Image_Processing_Tutorial_2.ipynb', 'balls.jpg', 'Image_Processing_Tutorial_1.ipynb']

为了解释惯性和展弦比的概念，我们使用这幅放大的手绘图：

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

<IPython.core.display.Javascript object>

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 = .3, label = "Red")
plt.hist(G.flatten(), bins = np.arange(256), histtype = "stepfilled", color = "g", alpha = .3, label = "Green")
plt.hist(B.flatten(), bins = np.arange(256), histtype = "stepfilled", color = "b", alpha = .3, label = "Blue")
plt.grid()
plt.legend()
plt.show()

<IPython.core.display.Javascript object>

阈值水平明显：

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

<IPython.core.display.Javascript object>

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


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

<IPython.core.display.Javascript object>

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

物体的惯性矩阵

下面表示的对象向某个方向拉伸。让我们看看如何使用它的惯性矩阵来确定它的拉伸方向和拉伸量。

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

<IPython.core.display.Javascript object>

二维物体的惯性矩阵定义如下：

\[我=\]

这个矩阵是对称的，因此，它可以在一个角度旋转的新框架中对角化。 \(\theta\) 在飞机上。该框架由两个归一化特征向量组成。 \((\vec e_1, \vec e_2)\) 矩阵的。在这个框架中，矩阵有两个特征值 \((I_1, I_2)\) 订购以便 \(I_1 \geq I_2\) .然后： * \(\theta = (\vec x, \vec e_1)\) 而且， * 纵横比 \(a = \sqrt{{I_1 / I_2}}\) .

角度 \(\theta\) 给出了物体的延伸方向和 \(a\) 显示它被拉长了多少。例如，如果 \(a == 1\) ，对象不会被拉长，而如果 \(a=10\) 在惯性视点上，方向1比方向2长10倍。

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)**.5

data[["area","theta", "aspect_ratio"]]

	area	theta	aspect_ratio
1	27036.0	125.441405	3.367242
2	17453.0	99.074781	3.837729
4	4204.0	131.093224	7.616207
6	3039.0	2.379166	1.955832
5	2736.0	82.039928	1.359502
8	2291.0	39.208508	1.710753
0	1438.0	25.883299	1.024153
7	1433.0	117.968405	1.265832
9	905.0	15.912720	1.151798
3	872.0	149.487814	1.112203

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()

<IPython.core.display.Javascript object>