Can 3D rotation be inferred from its effect on directional observations? I am thinking about a camera-based gyro that can infer the 3D orientation of an observation platform, by locking onto two (or more) known directional beacons.
Assume there is a reference coordinate axis system (with an origin and defined axis-directions for x-y and z.)  Now, assume a mobile platform (mobile in 3D such as a spaceship) having its own local x-y-z axes.  A full-spherical camera system on the platform can sense the direction of any observed directional beacon (a visible point like a very-distant galaxy, so far removed that it's apparent direction does not change observed over the entire nav-space.  Let there be a lookup table of reference directions for these beacons (directions in reference coordinates).  You can think of each beacon having a unique 3D direction vector d =  [ x y z ], or if you prefer, equivalent azimuth-elevation angle-pair [ phi theta ].
Can you solve this problem?  Using as few as two recognized beacons, you acquire (measure) the directions of these beacons in your local axis system:
d1 ,  d2   (measured directions with platform's 3D orientation superimposed)
You lookup from the reference galaxy chart the known reference directions for these two beacons:
d1_ref ,  d2_ref
Can you suggest a math formula or algorithm that crunches these 4 pieces of data into the 3D orientation of your observational platform (in relation to reference axes)?  I.e., can you turn your spherical hi-res camera into a gyroscope through software that combines your local directional observations to their known directions in reference coordinates?
An acceptable output would be a 3x3 rotational matrix, equivalent quaternion, or equivalent (roll, pitch, yaw) angles.
 A: Many problems of this type can be solved by remembering that the columns of a transformation matrix are the images of the basis vectors. Using this fact, you can easily construct a pair of matrices that map between the reference and camera coordinate systems, respectively, and some common intermediate coordinate system. You then compose one with the inverse of the other to get the direct transformation.  
In this case, you can use the coordinate system defined by the two measured reference objects as the common intermediate. For the third coordinate axis, take the cross product of the two vectors, which is guaranteed to be linearly independent of them (assuming that the two beacon vectors are linearly independent, i.e., not scalar multiples of each other). If $\mathbf d_1$ and $\mathbf d_2$ are the coordinates of the beacons in the reference frame and $\mathbf d_1'$, $\mathbf d_2'$ the coordinates measured by the camera, then the matrix that maps from the intermediate coordinate system to the reference coordinate system is simply: $$M_{\text{ref}}=\begin{bmatrix}\mathbf d_1 & \mathbf d_2 & \mathbf d_1\times\mathbf d_2 \end{bmatrix}$$ and from the intermediate coordinate system to the camera’s: $$M_{\text{cam}}=\begin{bmatrix} \mathbf d_1' & \mathbf d_2' & \mathbf d_1'\times\mathbf d_2' \end{bmatrix}.$$ We want to know what the the camera’s coordinate axes are in the reference frame, i.e., the images of the unit basis vectors under the camera-to-reference-frame mapping. Per the observation made at the start of this answer, these are the columns of $M_{\text{ref}}M_{\text{cam}}^{-1}$. This matrix maps from one orthonormal coordinate system to another, and, as is easily verified, $\det(M_{\text{ref}}M_{\text{cam}}^{-1})=1$, so it is indeed a rotation. Note that you don’t need unit vectors for this, nor does the order in which you take them matter, as long as it’s consistent between the two matrices.  
This is a theoretical result that assumes perfect data, but will work well with high-resolution measurements. In practice, to compensate for noise and other inaccuracies in the measurements, you’d more likely instead use observed positions of many more reference objects and compute a best-fit rotation from them. 
