<p>I think we can compute an inertia-like term of the cloud, the eigenvectors will tell us the two (major and minor) axes in 2D, and the eigenvalues will tell us the extents along the two axes.
This is probably the &quot;PCA idea&quot; mentioned by lhf in the comments.</p>
<p>I tried this idea out in Matlab. The code should work fine in Octave too.</p>
<p>First seed some points in an ellipse with scaling/rotation/translation. In your case you should use the real inputs. I only use this for testing.</p>
<pre><code>function [X,Y] = sample(a, b, s, q, tx, ty, n)
    X = (rand(n, 1) - 0.5) * s;
    Y = (rand(n, 1) - 0.5) * s;
    X2 = X .* X / (a * a);
    Y2 = Y .* Y / (b * b);
    i = (X2 + Y2) &lt; 1;
    X0 = X(i);
    Y0 = Y(i);
    X = cos(q) * X0 - sin(q) * Y0 + tx;
    Y = sin(q) * X0 + cos(q) * Y0 + ty;
end
</code></pre>
<p>On my side, I did this to generate the cloud: (also take the total number of points)</p>
<pre><code>[X,Y]=sample(2,1,4,pi/6, 1.9, 4.7,2000);
N=size(X)
N=N(1)
</code></pre>
<p>Then compute the center of the points and normalize the inputs:</p>
<pre><code>C=[sum(X)/N, sum(Y)/N];
nX=X-C(1);
nY=Y-C(2);
</code></pre>
<p>I will compute the inertia-like term using <code>nX</code> and <code>nY</code>.
Given a vector <span class="math-container">$\mathbf{r}$</span>, each inertia term is <span class="math-container">$\mathbf{I} \mathbf{r}^T \mathbf{r} - \mathbf{r} \mathbf{r}^T$</span>. Here <span class="math-container">$\mathbf{I}$</span> is the identity matrix.
In Matlab/Octave:</p>
<pre><code>function M = ptinertia2(X, Y)
    rtr = (X' * X + Y' * Y) * eye(2);
    P = [X, Y];
    rrt = P' * P;
    M = rtr - rrt;
end
</code></pre>
<p>The inertia and its eigen pairs are:</p>
<pre><code>M=ptinertia2(nX,nY);
[V,D]=eig(M)
f=sqrt(D(2,2)/D(1,1))
</code></pre>
<p>(<code>f</code> is the scaling factor on the major axis).</p>
<p>The eigenvalues are in the diagonal matrix <code>D</code>, and the eigenvectors are in the orthonormal matrix <code>V</code>.
To visualize the inputs and the identified axes:</p>
<pre><code>scatter(nX,nY,'.')
hold on
line([0, V(1,2)],[0,V(2,2)],'Color','red')
line([0, V(1,1)*f],[0,V(2,1)*f],'Color','green')
</code></pre>
<p><a href="https://i.sstatic.net/k2H7K.jpg" rel="nofollow noreferrer"><img src="https://i.sstatic.net/k2H7K.jpg" alt="sample output of the code" /></a></p>

CC BY-SA 4.0
2024-04-15T01:32:07.667Z
2024-04-15T01:32:07.667Z
TheBusyTypist
<p>I have a set of particles in a gravitationally bound system. The shape of the particles change with time, I have used a convex hull algorithm to &quot;enclose&quot; the particles of the cloud. I want to find the major and the minor axes of the cloud.</p>
<pre><code>hulls = [hull1, hull2, hull3, hull4, hull5]

for hull in hulls:
    for simplex in hull.simplices:
        plt.plot(hull.points[simplex, 0], hull.points[simplex, 1], '-', color='black')
    
    points = hull.points[hull.vertices] # shape: (#vertices, #coordinates)
    max_distance = 0
    major_axis = None
    for i in range(len(points)):
        for j in range(i + 1, len(points)):
            distance = np.linalg.norm(points[i] - points[j])
            if distance &gt; max_distance:
                max_distance = distance
                major_axis = (points[i], points[j])

    plt.plot([major_axis[0][0], major_axis[1][0]], [major_axis[0][1], major_axis[1][1]], color='black')
</code></pre>
<p><a href="https://i.sstatic.net/ndE9m.png" rel="nofollow noreferrer"><img src="https://i.sstatic.net/ndE9m.png" alt="enter image description here" /></a></p>
<p>The plot shows the cloud changing with time. Each cloud is for a different time.</p>
<p>I have been able to get the major axes from this. But not the minor axes. How can I do it? Basically just need the distance of the perpendicular through the middle of the major axis. This is a aspect ratio problem. I need the length and width of the cloud. If there is a better way to do this, please let me know.</p>

CC BY-SA 4.0
2024-03-28T14:26:06.710Z
PCA will give a good approximation in linear time.
lhf
CC BY-SA 4.0
2024-03-28T20:35:03.853Z
https://docs.opencv.org/4.x/de/d62/tutorial_bounding_rotated_ellipses.html
shawn_halayka
<computational-geometry>
How to find the major and minor axes of a convex hull?
Penguin1068
# How to find the major and minor axes of a convex hull?

(Asked by: username_0 on 2024-03-28)

I have a set of particles in a gravitationally bound system. The shape of the particles change with time, I have used a convex hull algorithm to "enclose" the particles of the cloud. I want to find the major and the minor axes of the cloud.

```
hulls = [hull1, hull2, hull3, hull4, hull5]

for hull in hulls:
    for simplex in hull.simplices:
        plt.plot(hull.points[simplex, 0], hull.points[simplex, 1], '-', color='black')
    
    points = hull.points[hull.vertices] # shape: (#vertices, #coordinates)
    max_distance = 0
    major_axis = None
    for i in range(len(points)):
        for j in range(i + 1, len(points)):
            distance = np.linalg.norm(points[i] - points[j])
            if distance > max_distance:
                max_distance = distance
                major_axis = (points[i], points[j])

    plt.plot([major_axis[0][0], major_axis[1][0]], [major_axis[0][1], major_axis[1][1]], color='black')

```

[![enter image description here](https://i.sstatic.net/ndE9m.png)](https://i.sstatic.net/ndE9m.png)

The plot shows the cloud changing with time. Each cloud is for a different time.

I have been able to get the major axes from this. But not the minor axes. How can I do it? Basically just need the distance of the perpendicular through the middle of the major axis. This is a aspect ratio problem. I need the length and width of the cloud. If there is a better way to do this, please let me know.

**Tags:** computational-geometry

> **Comments:**
> **username_1**: PCA will give a good approximation in linear time.  
> **username_2**: https://docs.opencv.org/4.x/de/d62/tutorial_bounding_rotated_ellipses.html  

---

## 1 answers

### Answer by username_3 on 2024-04-15. Score: 3 | Accepted answer

I think we can compute an inertia-like term of the cloud, the eigenvectors will tell us the two (major and minor) axes in 2D, and the eigenvalues will tell us the extents along the two axes.
This is probably the "PCA idea" mentioned by lhf in the comments.

I tried this idea out in Matlab. The code should work fine in Octave too.

First seed some points in an ellipse with scaling/rotation/translation. In your case you should use the real inputs. I only use this for testing.

```
function [X,Y] = sample(a, b, s, q, tx, ty, n)
    X = (rand(n, 1) - 0.5) * s;
    Y = (rand(n, 1) - 0.5) * s;
    X2 = X .* X / (a * a);
    Y2 = Y .* Y / (b * b);
    i = (X2 + Y2) < 1;
    X0 = X(i);
    Y0 = Y(i);
    X = cos(q) * X0 - sin(q) * Y0 + tx;
    Y = sin(q) * X0 + cos(q) * Y0 + ty;
end

```

On my side, I did this to generate the cloud: (also take the total number of points)

```
[X,Y]=sample(2,1,4,pi/6, 1.9, 4.7,2000);
N=size(X)
N=N(1)

```

Then compute the center of the points and normalize the inputs:

```
C=[sum(X)/N, sum(Y)/N];
nX=X-C(1);
nY=Y-C(2);

```

I will compute the inertia-like term using `nX` and `nY`.
Given a vector $\mathbf{r}$, each inertia term is $\mathbf{I} \mathbf{r}^T \mathbf{r} - \mathbf{r} \mathbf{r}^T$. Here $\mathbf{I}$ is the identity matrix.
In Matlab/Octave:

```
function M = ptinertia2(X, Y)
    rtr = (X' * X + Y' * Y) * eye(2);
    P = [X, Y];
    rrt = P' * P;
    M = rtr - rrt;
end

```

The inertia and its eigen pairs are:

```
M=ptinertia2(nX,nY);
[V,D]=eig(M)
f=sqrt(D(2,2)/D(1,1))

```

(`f` is the scaling factor on the major axis).

The eigenvalues are in the diagonal matrix `D`, and the eigenvectors are in the orthonormal matrix `V`.
To visualize the inputs and the identified axes:

```
scatter(nX,nY,'.')
hold on
line([0, V(1,2)],[0,V(2,2)],'Color','red')
line([0, V(1,1)*f],[0,V(2,1)*f],'Color','green')

```

[![sample output of the code](https://i.sstatic.net/k2H7K.jpg)](https://i.sstatic.net/k2H7K.jpg)


---

