"Normalizing" Points on a Sphere I have a set of points on a unit sphere representing different orientations:

Now I need to apply rotation(s) such that the points will lay on the horizon as tightly as possible:

The ideal representation for rotation would be YPR (Yaw-Pitch-Roll), but any axis-angle representation would suffice.
Here is where I got up to now:


*

*Find "centroid" orientation for the given orientations (How to find centroid of points laying on sphere?) The result would be some vector $u=(x,y,z)^{\mathbb{T}}$, $||u||=1$

*Compute normalizing yaw and pitch angles (maybe using atan2 ?) so that the rotated vector $u'=(0,0,1)$

*Compute roll angle (around $u'$) so that the points lay on the horizon line (the horizon is intersection of the sphere and a plane $y=0$). Maybe a least squares approach? Or compute average of normalizing angles for the points?
The problem arises in panoramic image stitching where I need to remove "wavy" effect from the projected mosaic (every image is represented by rotation around common center):

Note that the full 360° mosaic is a special case. I need to "straighten" partial mosaic as well, which are like the bunch of points depicted above.
 A: Here's one idea, trying hard not to have any bias stemming from the initial orientation of the sphere:


*

*Represent your points in euclidean 3-d coordinates with the unit sphere centered at $(0,0,0)$.

*Take cross product of all pairs of points (or a random sample of pairs if you have many points and $O(n^2)$ is too much work).

*Select the numerically largest cross product and adjust the sign of all others such that each makes an acute angle with it. 

*Add the sign-adjusted cross products. The direction of the resulting vector is the one you want to move to the north pole by your rotation.

A: A different idea, supposing that in your application the points already almost line up with the equator, such that they can be normalized by a small rotation around some horizontal axis:
Assume we have the points given in "geographical" coordinates with latitude $\phi$ and longitude $\lambda$. For points close to the equator, a small rotation about a horizontal axis will act like
$$ (\phi,\lambda) \mapsto (\phi - \alpha\cos(\lambda+\psi), \lambda) $$
where $\alpha$ measures the amount of rotation and $\psi$ is a phase term determining the orientation of the rotation axis.
Now if $\alpha$ and $\psi$ are both unknown we can represent the north-south adjustment instead as $a\cos \lambda + b\sin\lambda$ for some (also unknown) coefficients $a$ and $b$. However in this form we can estimate the best $a$ and $b$ by a standard least-squares fit to the $\phi$ coordinates of our known points.
Once $a$ and $b$ are known, it's a matter of standard trigonometry to derive $\alpha$ and $\psi$ from them. We can then do the actual rotation of the sphere by rotating by angle $\psi$ around the $z$ axis and then by $\alpha$ around the $x$ axis.
