randomAffine3d is not properly random

Question

Carson Purnell 2024-11-7，19:53

2
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2164715-randomaffine3d-is-not-properly-random

编辑： Bruno Luong 2024-11-8，19:10

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

2 个评论
显示无隐藏无

Cris LaPierre 2024-11-7，20:10

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 2024-11-7，21:24

在 MATLAB Online 中打开

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 2024-11-7，20:30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2164715-randomaffine3d-is-not-properly-random#answer_1542565

在 MATLAB Online 中打开

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

Bruno Luong 2024-11-8，7:31

编辑：Bruno Luong 2024-11-8，8:38

在 MATLAB Online 中打开

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 2024-11-8，13:06

编辑：Bruno Luong 2024-11-8，13:07

在 MATLAB Online 中打开

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 2024-11-7，20:56

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2164715-randomaffine3d-is-not-properly-random#answer_1542580

编辑：Matt J 2024-11-7，23:24

在 MATLAB Online 中打开

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

4 个评论
显示 2更早的评论隐藏 2更早的评论

Bruno Luong 2024-11-8，12:54

编辑：Bruno Luong 2024-11-8，13:03

在 MATLAB Online 中打开

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 2024-11-8，16:49

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 2024-11-8，19:06

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2164715-randomaffine3d-is-not-properly-random#answer_1543083

编辑：Bruno Luong 2024-11-8，19:10

在 MATLAB Online 中打开

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?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

randomAffine3d is not properly random

2 个评论
显示无隐藏无

回答（3 个）

5 个评论
显示 3更早的评论隐藏 3更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

randomAffine3d is not properly random

2 个评论 显示 无隐藏 无

回答（3 个）

5 个评论 显示 3更早的评论隐藏 3更早的评论

4 个评论 显示 2更早的评论隐藏 2更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

2 个评论
显示无隐藏无

5 个评论
显示 3更早的评论隐藏 3更早的评论

4 个评论
显示 2更早的评论隐藏 2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论