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.
Here's one idea, trying hard not to have any bias stemming from the initial orientation of the sphere: