Anyway here is the rotation code as well.
% rotations = find_all(ndim, NbRotations)
%
% ndim = 2 or 3 (but the code can be changed)
%
%
% code for generating an optimised sequence of rotations
% for the IDT pdf transfer algorithm
% although the code is not beautiful, it does the job.
%
function rotations = find_all(ndim, NbRotations)
if (ndim == 2)
l = [0 pi/2];
elseif (ndim == 3)
l = [0 0 pi/2 0 pi/2 pi/2];
else % put here initialisation stuff for higher orders
end
fprintf('rotation ');
for i = 1
NbRotations-1)fprintf('%d ...', i );
l = [l ; find_next(l, ndim)]; %l(end,
+ones(1,ndim-1)*pi/2)]fprintf('\b\b\b', i );
end
M = ndim;
rotations = cell(1,NbRotations);
for i=1:size(l, 1)
for j=1:M
b_prev(j,
= hyperspherical2cartesianT(l(i,(1:ndim-1) + (j-1)*(ndim-1)));end
b_prev = grams(b_prev')';
rotations{i} = b_prev;
end
end
%
%
function [x] = find_next( list_prev_x, ndim)
prevx = list_prev_x; % in hyperspherical coordinates
nprevx = size(prevx,1);
hdim = ndim - 1;
M = ndim;
% convert points to cartesian coordinates
c_prevx = zeros(nprevx*M, ndim);
c_prevx = [];
for i=1:nprevx
for j=1:M
b_prev(j,
= hyperspherical2cartesianT(prevx(i,(1:hdim) + (j-1)*hdim));end
b_prev = grams(b_prev')';
c_prevx = [c_prevx; b_prev];
end
c_prevx;
options = optimset('TolX', 1e-10);
options = optimset(options,'Display','off');
minf = inf;
for i=1:10
x0 = rand(1, hdim*M)*pi - pi/2;
x = fminsearch(@myfun, x0, options);
f = myfun(x);
if f < minf
minf = f;
mix = x;
end
end
%%
% f - Compute the function value at x
function [f] = myfun(x)
% compute the objective function
c_x = zeros(M, ndim);
for i=1:M
c_x(i,
= hyperspherical2cartesianT(x((1:hdim) + (i-1)*hdim));end
c_x = grams(c_x')';
f = 0;
for i=1:M
for p=1:size(c_prevx, 1)
d = (c_prevx(p,
- c_x(i, ) * (c_prevx(p, - c_x(i, )';f = f + 1/(1 + d);
d = (c_prevx(p,
+ c_x(i, ) * (c_prevx(p, + c_x(i, )';f = f + 1/(1 + d);
end
end
end
%%
end
%
%
function c = hyperspherical2cartesianT(x)
c = zeros(1, length(x)+1);
sk = 1;
for k=1:length(x)
c(k) = sk*cos(x(k));
sk = sk*sin(x(k));
end
c(end) = sk;
end
% Gram-Schmidt orthogonalization of the columns of A.
% The columns of A are assumed to be linearly independent.
function [Q, R] = grams(A)
[m, n] = size(A);
Asave = A;
for j = 1:n
for k = 1:j-1
mult = (A(:, j)'*A(:, k)) / (A(:, k)'*A(:, k));
A(:, j) = A(:, j) - mult*A(:, k);
end
end
for j = 1:n
if norm(A(:, j)) < sqrt(eps)
error('Columns of A are linearly dependent.')
end
Q(:, j) = A(:, j) / norm(A(:, j));
end
R = Q'*Asave;
end