Is this fundamentally the same or similar problem as the rectangular stretch mentioned here? But with a "jacobian" stretch as opposed to a linear stretch?
Could you please read the short section https://en.wikipedia.org/wiki/Orthogonal_matrix#Randomizatio... ?
It is very short, first it defines what is meant with a random orthogonal matrix, then it describes the problem with an incorrect approach, ... and then it describes a recursive procedure to generate a random orthogonal matrix of dimension D using D * (D+1) / 2 uniformly distributed numbers:
>Stewart (1980) replaced this with a more efficient idea that Diaconis & Shahshahani (1987) later generalized as the "subgroup algorithm" (in which form it works just as well for permutations and rotations). To generate an (n + 1) × (n + 1) orthogonal matrix, take an n × n one and a uniformly distributed unit vector of dimension n + 1. Construct a Householder reflection from the vector, then apply it to the smaller matrix (embedded in the larger size with a 1 at the bottom right corner).
I predict this will give a better uniform distribution on that sphere, if you share your code for generation and visualization I might try myself :)
Also, I applaud you for the extremely clear exposition of the problem, thought process and analogies presented. The one section I found harder to grasp was the one on dithering, as I don't fully understand the process from input picture to output pixel locations and intensities...