A moment (M) is a specific quantitative measure of a given distribution. The distribution can be one or multi-dimensional (wikipedia definition).
Now let's try to have an intuitive view of what it is in the case of an image.
Remember that an image is only a 2D distribution of a set of points. Each point being localized by two coordinates $x$ and $y$. So fo each pixel we have the following $f(x, y) = I_{x, y}$. Where $I_{x, y}$ corresponds to the intensity at the coordinate $x$ and $y$.
Before going to intuition let's do some easy math. The formal definition of the image moment (2D and continuous) is:
$$M_{pq} = \int_{-\infty}^{\infty} x^p y^q f(x, y) dx dy$$And the 2D discrete equation is :
$$M_{pq} = \sum_x \sum_y x^p y^q I(x,y)$$So we see here the moment (M) depends on two parameters $p$ and $q$. They are called the order of the moment. Different order will give a different measure of the distribution we are looking at.
For example let's define a simple 1D distribution of length $N$ $f(x) = A$.
Now if we divide the first order moment by the zero order moment we have the mean :
$$\bar{A} = \frac{1}{N} \sum_x x = \frac{M_1}{M_0}$$
In [28]:
# Do some import
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from tifffile import TiffFile
In [48]:
# Load the image
tf = TiffFile("binary_cell.tif")
# Get the numpy array
a = tf.asarray()
# Replace all 255 to 1 so the image is now made of "0" and "1"
a[a == 255] = 1
print(np.unique(a))
_ = plt.imshow(a, interpolation="none", cmap="gray", origin='lower')
The goal here is to use the image moment to fit an ellipse. The ellipse can be defined by four parameters : the center, the orientation, the major axis and the minor axis.
Note it exists many different method to fit an ellipse.
First we are going to detect the centroid (the center of the ellipse) of our object by computing the zero and first order of the image.
Remember the discrete equation :
$$M_{pq} = \sum_x \sum_y x^p y^q I(x,y)$$Since we are working on a binary image (only $0$ and $1$) values we remove $I(x,y)$.
The zero order moment is then :
$$M_{00} = \sum_x \sum_y x^0 y^0$$Let's do it in python.
In [49]:
# The sum of all value of 1
M_00 = np.sum(a)
M_00
Out[49]:
Now let's compute the first order moment for $x$ :
$$M_{10} = \sum_x \sum_y x^1 y^0$$
In [50]:
# Here we get all the coordinates of pixel equal to 1
xx, yy = np.where(a == 1)
M_10 = np.sum(xx)
M_10
Out[50]:
Now let's compute the first order moment for $y$ :
In [51]:
M_01 = np.sum(yy)
M_01
Out[51]:
So the centroid $C$ is given by :
$$C_x = \bar{x} = \frac{M_{10}}{M_{00}}$$$$C_y = \bar{y} = \frac{M_{01}}{M_{00}}$$
In [68]:
C_x = M_10 / M_00
C_y = M_01 / M_00
print("C =", (C_x, C_y))
Let's verify it visually :
In [69]:
plt.scatter(C_x, C_y, color='red', marker='+', s=100)
_ = plt.imshow(a, interpolation="none", cmap="gray", origin='lower')
Well it seems to be exact !
Here it becomes a little bit tricky.
The information about the object orientation can be derived using the second order central moments to construct a covariance matrix (see Wikipedia for more details).
Here is a new concept : central moment ($\mu$) which describes distribution of mean (unlike moment (M) which describes distribution only).
Note that sometime moment is also called raw moment.
The discrete equation of the central moment is the :
$$\mu_{pq} = \sum_{x} \sum_{y} (x - \bar{x})^p(y - \bar{y})^q f(x,y)$$with $\bar{x}=\frac{M_{10}}{M_{00}}$ and $\bar{y}=\frac{M_{01}}{M_{00}}$
Now it becomes difficult to get the intuition of what's going next. From Wikipedia :
The covariance matrix of the image $I(x,y)$ is :
$$\operatorname{cov}[I(x,y)] = \begin{bmatrix} \mu'_{20} & \mu'_{11} \\ \mu'_{11} & \mu'_{02} \end{bmatrix}$$The eigenvectors of the covariance matrix of the image correspond to the major and minor axes of the image intensity, so the orientation can thus be extracted from the angle of the eigenvector associated with the largest eigenvalue. It can be shown that this angle Θ is given by the following formula:
$$\Theta = \frac{1}{2} \arctan \left( \frac{2\mu'_{11}}{\mu'_{20} - \mu'_{02}} \right)$$Where :
$$\mu'_{20} = \mu_{20} / \mu_{00} = M_{20}/M_{00} - \bar{x}^2$$$$\mu'_{02} = \mu_{02} / \mu_{00} = M_{02}/M_{00} - \bar{y}^2$$$$\mu'_{11} = \mu_{11} / \mu_{00} = M_{11}/M_{00} - \bar{x}\bar{y}$$
Let's first compute the second order raw moment ($M_{20}$, $M_{02}$ and $M_{11}$) :
In [70]:
M_20 = np.sum(xx ** 2)
M_02 = np.sum(yy ** 2)
M_11 = np.sum(xx * yy)
Compute $\mu'_{20}$, $\mu'_{02}$ and $\mu'_{11}$ :
In [71]:
mu_20 = M_20 / M_00 - C_x ** 2
mu_02 = M_20 / M_00 - C_y ** 2
mu_11 = M_11 / M_00 - C_x * C_y
Get the orientation $\Theta$ :
$$\Theta = \frac{1}{2} \arctan \left( \frac{2\mu'_{11}}{\mu'_{20} - \mu'_{02}} \right)$$
In [72]:
theta = 1/2 * np.arctan((2 * mu_11) / (mu_20 - mu_02))
# Convert it in degree
angle = np.rad2deg(theta)
print("angle = {}°".format(angle))
Get the eigenvalues with this equation :
$$\Delta = \sqrt{4{\mu'}_{11}^2 + ({\mu'}_{20}-{\mu'}_{02})^2}$$$$\lambda_i = \frac{\mu'_{20} + \mu'_{02}}{2} \pm \frac{\Delta}{2}$$
In [73]:
delta = np.sqrt(4 * mu_11 ** 2 + (mu_20 - mu_02) ** 2)
lambda_1 = ((mu_20 + mu_02) + delta) / 2
lambda_2 = ((mu_20 + mu_02) - delta) / 2
Get the major and minor axes from the eigenvalues :
In [74]:
semi_major_length = np.sqrt(np.abs(lambda_1)) * 2
semi_minor_length = np.sqrt(np.abs(lambda_2)) * 2
major_length = semi_major_length * 2
minor_length = semi_minor_length * 2
print("Major axis = {}".format(major_length))
print("Minor axis = {}".format(minor_length))
Alternatively we can also compute (with Numpy) the orientation and axes length from the eigenvectors of the covariance matrix.
$$\operatorname{cov}[I(x,y)] = \begin{bmatrix} \mu'_{20} & \mu'_{11} \\ \mu'_{11} & \mu'_{02} \end{bmatrix}$$
In [75]:
cov = np.asarray([[mu_20, mu_11],
[mu_11, mu_02]])
eigvalues, eigvectors = np.linalg.eig(cov)
# Get the associated eigenvectors and eigenvalues
eigval_1, eigval_2 = eigvalues
eigvec_1, eigvec_2 = eigvectors[:, 0], eigvectors[:, 1]
Get the orientation from the first eigenvector with $ \Theta = \operatorname{atan2}(y, x) \quad$ and the axes length from the eigenvalues.
In [76]:
theta = np.arctan2(eigvec_1[1], eigvec_1[0])
angle = np.rad2deg(theta)
print("angle = {}°".format(angle))
semi_major_length = np.sqrt(np.abs(eigval_1)) * 2
semi_minor_length = np.sqrt(np.abs(eigval_2)) * 2
major_length = semi_major_length * 2
minor_length = semi_minor_length * 2
print("Major axis = {}".format(major_length))
print("Minor axis = {}".format(minor_length))
We find the same value as above.
Now let's see how our vector looks (which is supposed to define the major orientation of our object).
In [77]:
plt.figure()
# Plot the centroid in red
plt.scatter(C_x, C_y, color='red', marker='+', s=100)
# Plot the first eigenvector
scale = 100
x1, x2 = [C_x, eigvec_1[0] * scale]
y1, y2 = [C_y, eigvec_1[1] * scale]
plt.arrow(x1, y1, x2, y2, color='green', lw=1, head_width=20)
# Show the image
_ = plt.imshow(a, interpolation="none", cmap="gray", origin='lower')
The orientation seems to be correct.
Now we draw the major and minor axes.
In [78]:
fig, ax = plt.subplots()
# Plot the centroid in red
ax.scatter(C_x, C_y, color='red', marker='+', s=100)
# Plot the major axis
x1 = C_x + semi_major_length * np.cos(theta)
y1 = C_y + semi_major_length * np.sin(theta)
x2 = C_x - semi_major_length * np.cos(theta)
y2 = C_y - semi_major_length * np.sin(theta)
ax.plot([x1, x2], [y1, y2], color='green', lw=1)
# Plot the minor axis
x1 = C_x + semi_minor_length * np.cos(theta + np.pi/2)
y1 = C_y + semi_minor_length * np.sin(theta + np.pi/2)
x2 = C_x - semi_minor_length * np.cos(theta + np.pi/2)
y2 = C_y - semi_minor_length * np.sin(theta + np.pi/2)
ax.plot([x1, x2], [y1, y2], color='green', lw=1)
# Plot the ellipse
angles = np.arange(0, 360, 1) * np.pi / 180
x = 0.5 * major_length * np.cos(angles)
y = 0.5 * minor_length * np.sin(angles)
R = np.array([[np.cos(theta), -np.sin(theta)],
[np.sin(theta), np.cos(theta)]])
x, y = np.dot(R, np.array([x, y]))
x += C_x
y += C_y
ax.plot(x, y, lw=1, color='red')
# Show the image
_ = ax.imshow(a, interpolation="none", cmap="gray", origin='lower', aspect='equal')
It seems to work but there is an underestimation of the axes length...
The method is not exactly the same but this understimation does not appear in the ImageJ/Fiji algoritm.
Any modifications are welcome !
Some links which help me :
Another method to detect ellipse :