randomAffine3d is not properly random

the transformations are highly biased, and don't cover the unit sphere with any angle coverage. Image shows 4 compass points randomly rotated +-45. Code replicates the image, as well as a 'correctly' random version.

is this a bug, or was it a bug now fixed (i am using 2019b)

ori = [0,0,1;0,0,-1;1,0,0;-1,0,0];

%% randomaffine3d is NOT usefully random (highly biased, missing regions)
o = [];
for i=1:10000
tform = randomAffine3d('Rotation',[-45 45]);
mat = tform.T(1:3,1:3); r = ori*mat; % alternative to prove TPF isn't broken
r = transformPointsForward(tform,ori);
o(end+1:end+size(ori,1),:) = r;
end
plot3(o(:,1),o(:,2),o(:,3),'.'); axis equal

%% manually random rotation matrix - correctly random
b = [];
for i=1:10000
    ax = randn(1,3); ax = ax./norm(ax); %true random unit vectors
    mat = rotmat(ax,rand*pi/4); % random rotation matrix (but not isotropic)
    r = ori*mat;
    b(end+1:end+size(ori,1),:) = r;
end
figure(); plot3(b(:,1),b(:,2),b(:,3),'.'); axis equal

%% internal funct
function t = rotmat(ax,rad)
ax = ax/norm(ax);
x = ax(1); y = ax(2); z = ax(3);
c = cos(rad); s = sin(rad);

t1 = c + x^2*(1-c);
t2 = x*y*(1-c) - z*s;
t3 = x*z*(1-c) + y*s;
t4 = y*x*(1-c) + z*s;
t5 = c + y^2*(1-c);
t6 = y*z*(1-c)-x*s;
t7 = z*x*(1-c)-y*s;
t8 = z*y*(1-c)+x*s;
t9 = c+z^2*(1-c);

t = [t1 t2 t3
    t4 t5 t6
    t7 t8 t9];
end

[Cris LaPierre]

I ran your code so you can see the results in R2024b.

Note that you are asking the community. If you want to report a concern directly to MathWorks, please contact support: https://www.mathworks.com/support/contact_us.html

[Bruno Luong]

Alternative correct code

ori = [0,0,1;
    0,0,-1;
    1,0,0;
    -1,0,0];

%% manually random rotation matrix - correctly random
b = [];
for i=1:10000
    mat = rotmat2(randn(1,3),rand*pi/4); % random rotation matrix (but not isotropic)
    r = ori*mat;
    b(end+1:end+size(ori,1),:) = r;
end
figure(); plot3(b(:,1),b(:,2),b(:,3),'.'); axis equal

%% internal funct
function t = rotmat2(ax,rad)
M = makehgtform("axisrotate",ax,rad);
t = M(1:3,1:3);
end

Answer 1

[Matt J]

If you want an isotropic distribution across the entire sphere, this should do it.

ori = [eye(3);-eye(3)];

N=1000;
o=cell(N,1);
for i=1:N
    [Q,R]=qr(rand(3));
    Q=Q*det(Q);
    r = ori*Q;
   o{i} = r;
end
o=vertcat(o{:});

plot3(o(:,1),o(:,2),o(:,3),'.'); axis equal

[Bruno Luong]

No this is NOT uniform, hard to see until you plot some histogram

ori = [eye(3);-eye(3)];
N=1000000;
o=cell(N,1);
for i=1:N
    [Q,R]=qr(rand(3));
    Q=Q*det(Q);
    r = ori*Q;
   o{i} = r;
end
o=vertcat(o{:});
az = atan2(o(:,2),o(:,1));
histogram(az)

To be uniform one should use RANDN and not RAND

ori = [eye(3);-eye(3)];
N=1000000;
o=cell(N,1);
for i=1:N
    [Q,R]=qr(randn(3));
    Q=Q*det(Q);
    r = ori*Q;
    o{i} = r;
end
o=vertcat(o{:});
az = atan2(o(:,2),o(:,1));
histogram(az)

[Bruno Luong]

Or simply

N=1000000;
o = randn(6*N,3);
o = o ./ vecnorm(o,2,2);
az = atan2(o(:,2),o(:,1));
histogram(az)

[Matt J]

I think I understand why the original version isn't right, but I'm not sure why the az histogram confirms that the counter-proposals are. The results are uniform in az, but not el.

ori = [eye(3);-eye(3)];
N=1e4;
o=cell(N,1);
for i=1:N
    [Q,R]=qr(randn(3));
    Q=Q*det(Q);
    r = ori*Q;
    o{i} = r;
end
o=vertcat(o{:});

[az,el] = cart2sph(o(:,1), o(:,2), o(:,3));

subplot(1,2,1) 
histogram(az)  ; axis square; xlabel AZ
subplot(1,2,2)
histogram(el) ; axis square; xlabel EL

[Bruno Luong]

When you slide the sphere with delta el the area is not uniform as with az.

The delta area is cos(el) so what you see is histogram(el) ~ cos(el) when the points are uniformly distributed on the sphere.

N = 1e6;
o = randn(6*N,3);
o = o ./ vecnorm(o,2,2);

figure
[az,el] = cart2sph(o(:,1), o(:,2), o(:,3));
h = histogram(el,'Normalization','percentage');
mx = max(h.Values);
hold on

EZ = linspace(-pi/2,pi/2);
plot(EZ, mx*cos(EZ),'r','LineWidth',2)

As criteria you could also check uniformity of

figure
histogram(atan2(o(:,3),o(:,1))), 

or any projection of o on two arbitrary orthogonal unit vectors.

[Bruno Luong]

Generate unfiform points on S^2 using spherical coordinates

N = 1e6;
az = 2*pi*rand(1,N);
el = asin(2*rand(1,N)-1);
r = ones(1,N);
[x,y,z] = sph2cart(az,el,r);
figure
plot3(x,y,z,'.')
axis equal

figure
histogram(atan2(z,x))

Answer 2

[Matt J]

EDITED: If you don't like what the default randomizer is doing, randomAffine3d() also let's you define your own, e.g.,

ori = [eye(3);-eye(3)];

N=1e4;
o=cell(N,1);
for i=1:N
    tform = randomAffine3d('Rotation',@selectRotation);
    mat = tform.T(1:3,1:3);
    r = transformPointsForward(tform,ori);
   o{i} = r;
end
o=vertcat(o{:});

[az,el] = cart2sph(o(:,1), o(:,2), o(:,3));

subplot(1,2,1) 
histogram(az)  ; axis square; xlabel AZ
subplot(1,2,2)
histogram(el) ; axis square; xlabel EL

function [rotationAxis,theta] = selectRotation()
 

  [Q,~]=qr(randn(3));
  Q=Q*det(Q);
  [V,d]=eig(Q,'vector');

   [~,i]=max(real(d));

   rotationAxis=real(V(:,i)');

   B=[null(rotationAxis),rotationAxis'];
   B(:,1)=B(:,1)*det(B);
   c=B'*Q*B(:,1);

   theta=atan2d(c(2),c(1));

end

[Bruno Luong]

Using rand on el is again incorrect (will not generated uniform points on sphere)

[Matt J]

I guess it's right now, seeing as it gives the same histograms as the other answer.

[Bruno Luong]

One must convert angle to degree as for interface with randomAffine3d

To convert rotation matrix (Q) to rotation angle/ vector better method is using quaternion or Caykey method:

https://www.mathworks.com/matlabcentral/answers/1994988-difference-between-2-orientations

But in this thread one can generate directly rotation vector uniform on S^2 and rotation angle, uniform as user prescribed as with randomAffine3d('Rotation',[angledeg_lo angle_hi])

ori = [0,0,1;
    0,0,-1;
    1,0,0;
    -1,0,0];

N=1e4;
o=[];
for i=1:N
    tform = randomAffine3d('Rotation',@selectRotation2);
    mat = tform.T(1:3,1:3);
    r = transformPointsForward(tform,ori);
   o(end+1:end+size(ori,1),:) = r;
end
figure(); plot3(o(:,1),o(:,2),o(:,3),'.'); axis equal

%% internal funct
function [rotationAxis,theta] = selectRotation2()
rotationAxis = randn(1,3);
rotationAxis = rotationAxis / norm(rotationAxis);
theta = 45 * (2*rand()-1);
end

[Bruno Luong]

There is a not well specification the the doc of randomAffine3d

In the "Rotation" section one can read

"A function handle of the form

[rotationAxis,theta] = selectRotation

The function selectRotation must accept no input arguments. The function must return two output arguments: rotationAxis, a 3-element vector defining the axis of rotation, and theta, a rotation angle in degrees."

Actually it accepts only a 3-element row vector, a column vector will crash the function.

Answer 3

[Bruno Luong]

This is not an answer but I puspose introduce a bug that generates the rotation axis in the positivve part oof S^2 and it reproduces the same figure as reported in the question:

ori = [0,0,1;
    0,0,-1;
    1,0,0;
    -1,0,0];

N=1e4;
o=[];
for i=1:N
    tform = randomAffine3d('Rotation',@selectRotation2);
    mat = tform.T(1:3,1:3);
    r = transformPointsForward(tform,ori);
   o(end+1:end+size(ori,1),:) = r;
end
figure(); plot3(o(:,1),o(:,2),o(:,3),'.'); axis equal

%% internal funct
function [rotationAxis,theta] = selectRotation2()
az = 2*pi*rand();
el = asin(2*rand()-1);
r = 1;
[x,y,z] = sph2cart(az,el,r);
rotationAxis = [x,y,z]; 
rotationAxis = abs(rotationAxis); % Bug introduced with ABS
theta = 45 * (2*rand()-1);
end

For those who has the toolbox, Is randomAffine3d code is inn an lfiile and visible?