Let $AB$, $CD$ be the principal axes of the ellipse, intersecting at center $O$, and let $P$ be any point on the ellipse, $H$ and $K$ its projections on the axes (see figure below). We have:
$$
{OH^2\over OB^2}+{OK^2\over OD^2}=1.
$$
Let now $A'$, $B'$, ... be the perpendicular projections of points $A$, $B$, ... on a given plane. As perpendicular projections preserve the ratios of segments on a line, we have:
$$
{O'H'^2\over O'B'^2}+{O'K'^2\over O'D'^2}=1.
$$
But this equation means that point $P'$ belongs to the ellipse having $A'B'$ and $C'D'$ as conjugate diameters. You can find the principal axes of that ellipse, for instance, following the method explained here.
The following is based on a common proof method for Green's theorem.
I assume you already have a coordinate system or a labelling system for the white points in your diagram. When I say "label", I mean this set of coordinates or this label.
We are going to make a bag of directed edges between pairs of white points that lie on a side of a square. Directed edges are pairs, where the pair $(a,b)$ means the edge from the vertex labelled "$a$" to the vertex labelled "$b$". For each square, add its four edges to the bag:
- (northeast label, southeast label),
- (southeast label, southwest label),
- (southwest label, northwest label),
- (northwest label, northeast label).
This ensures that each boundary edge is included exactly once in the bag and each interior edge is included twice, once as $(a,b)$ and once as $(b,a)$. So we need to scan through the bag deleting interior edges.
For each edge, $e = (a,b)$, determine whether $(b,a)$ is also in the bag. If so delete both of them, otherwise, proceed to the next edge. (Notice that we fetch the edges from the bag in whatever order an iterator of the bag hands us the edges -- we have made no attempt to string the perimeter together into a path yet.)
Now, use min and max functions to find the right-most occupied column of squares (coordinate of right-most blue circle) and of those, the upper-most member of the column. Take the northeast corner of that square and call it $a$. (There is no square more north or more east than this square, so this vertex is on the perimeter.)
Search the bag for the edge (and there is only one) having $a$ as its first member. That edge is $(a,b)$, for some $b$, and is the first edge in the perimeter. Now search for the edge with first member $b$, obtaining the edge $(b,c)$, the second edge of the perimeter. Continue stepping along the alternating path of edge, vertex, edge, vertex, ..., until the search turns up an edge with $a$ as its second member, $(z,a)$. This is the last edge of the perimeter and closes up the perimeter.
Best Answer
Revised: From the clarifications it appears that the conical approximation of data points is only a rough one. However the apex is known, which leads to a revised suggestion of projecting points onto a sphere centered at the apex, which would lead to a rough circle of points on the sphere (these points would be exactly on a "small" circle of the sphere if the original data were actually on a cone).
Given the points on a sphere, one might proceed by fitting a small circle (great circles on a sphere being a special case) in a total least squares sense, i.e. minimizing the sum of squares of spherical distances from points to the fitted circle. This is a nonlinear problem, but it involves only three parameters (the radius and two coordinates for center of the circle). See Fujiki and Akaho, 2007 for exposition of both the 2D and higher dimensional analogs, and its references to the literature.
Your interest seems primarily to be identification of the axis of rotation, i.e. the pole that corresponds to the spherical circle. See here for some software tailored to this problem.
I think we should be careful about an approach that assumes points having a roughly uniform distribution over (around) the surface of a cone. If they were, it would make spotting a "V" shape in a slice of the data more likely, but perhaps there are more robust approaches (less dependent on an even sampling of points).
[Something else to keep in mind is that mathematically a cone proceeds in both directions of the axis through its apex. Therefore an ideal slice of the data through a plane containing the apex might reveal an "X" rather than a "V" shape.]
One possibility would be to assess for a candidate "apex" point $A$ of a cone, what axis gives the most tightly clustered slopes of lines passing through the sample points and the "apex". The axis is generated by a unit vector $U$ from $A$, and the slopes through sample point $P$ can be associated with a dot product of $U$ with $P - A$ divided by $||P - A||$.
In fact we might weight these slopes so that samples close to $A$ are less important than points far away (since the points close to $A$ are close to the cone regardless of what axis is chosen), and this suggests weighting by $||P - A||$, thereby cancelling the division used above.
This replaces the search for a "good slice" of the data with a search for a good approximation of the cone's apex.
The problem is of some importance in the field of computer vision and robotics. In this book chapter by Chernov and Ma, Least Squares Fitting of Quadratic Curves and Surfaces (2011), some improvements are sought using a "fast projection" method of Dave Eberly. The authors treat three non-degenerate quadratic surfaces in 3D, ellipsoids, hyperbolic paraboloids, and hyperboloids of one sheet, and mention without details that elliptic paraboloids and hyperboloids of two sheets can be treated similarly. A cone, however, is a degenerate quadric surface, intermediate between the hyperboloids of one and two sheets. I'd expect that the procedures described by Chernov and Ma would adapt readily to this specialization.