<jupyter_start><jupyter_text>What is the image moment (or how to fit and ellipse)A moment (M) is a specific quantitative measure of a given distribution. The distribution can be one or multi-dimensional ([wikipedia definition](https://en.wikipedia.org/wiki/Moment_(mathematics)).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$.- the zero order ($p=0$) moment of $A$ is **the number of element** : $M_0 = \sum_x x^0 = len(A)$ (for pythonist)- the first order ($p=1$) moment of a $A$ if the **sum** of all the element : $M_1 = \sum_x x^1 = np.sum(A)$ (for pythonist)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}$$  The image moment Now let's try to apply this 1D example to an image.<jupyter_code># Do some import
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
from tifffile import TiffFile
# 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')<jupyter_output>[0 1]<jupyter_text>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. Find the centroidFirst 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.<jupyter_code># The sum of all value of 1
M_00 = np.sum(a)
M_00<jupyter_output><empty_output><jupyter_text>Now let's compute the first order moment for $x$ :$$M_{10} = \sum_x \sum_y x^1 y^0$$<jupyter_code># Here we get all the coordinates of pixel equal to 1
xx, yy = np.where(a == 1)

M_10 = np.sum(xx)
M_10<jupyter_output><empty_output><jupyter_text>Now let's compute the first order moment for $y$ :<jupyter_code>M_01 = np.sum(yy)
M_01<jupyter_output><empty_output><jupyter_text>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}}$$<jupyter_code>C_x = M_10 / M_00
C_y = M_01 / M_00
print("C =", (C_x, C_y))<jupyter_output>C = (327.43648547085388, 337.30288897161603)<jupyter_text>Let's verify it visually :<jupyter_code>plt.scatter(C_x, C_y, color='red', marker='+', s=100)
_ = plt.imshow(a, interpolation="none", cmap="gray", origin='lower')<jupyter_output><empty_output><jupyter_text>Well it seems to be exact ! Find the major and minor axisHere 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](https://en.wikipedia.org/wiki/Image_moment) 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](https://en.wikipedia.org/wiki/Image_moment) :  > 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}$) :<jupyter_code>M_20 = np.sum(xx ** 2)
M_02 = np.sum(yy ** 2)
M_11 = np.sum(xx * yy)<jupyter_output><empty_output><jupyter_text>Compute $\mu'_{20}$, $\mu'_{02}$ and $\mu'_{11}$ :<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Get the **orientation** $\Theta$ : $$\Theta = \frac{1}{2} \arctan \left( \frac{2\mu'_{11}}{\mu'_{20} - \mu'_{02}} \right)$$<jupyter_code>theta = 1/2 * np.arctan((2 * mu_11) / (mu_20 - mu_02))

# Convert it in degree
angle = np.rad2deg(theta)
print("angle = {}°".format(angle))<jupyter_output>angle = -36.72080335329064°<jupyter_text>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}$$<jupyter_code>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<jupyter_output><empty_output><jupyter_text>Get the major and minor axes from the eigenvalues :<jupyter_code>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))<jupyter_output>Major axis = 590.4638690697747
Minor axis = 139.87004447009642<jupyter_text>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}$$<jupyter_code>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]<jupyter_output><empty_output><jupyter_text>Get the orientation from the **first eigenvector** with $ \Theta = \operatorname{atan2}(y, x) \quad$ and the **axes length** from the **eigenvalues**.<jupyter_code>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))<jupyter_output>angle = -36.72080335329064°
Major axis = 590.4638690697747
Minor axis = 139.87004447009653<jupyter_text>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).<jupyter_code>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')<jupyter_output><empty_output><jupyter_text>The **orientation** seems to be correct. Now we draw the major and minor axes.<jupyter_code>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')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.16438497169325786, 'fraction_numerical': 0.04673185795162121, 'mean_word_length': 4.906382978723404, 'pattern_counts': {'":': 0, '<': 56, '<?xml version=': 0, '>': 64, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '510620', 'license': 'MIT', 'n_tokens_mistral': 3944, 'n_tokens_neox': 3590, 'n_words': 1094, 'path': 'Analysis/Fit_Ellipse/notebook.ipynb', 'repo_name': 'hadim/Public-Notebooks', 'size': '9715'}