frcs-colour-transfer

1
% rotations = find_all(ndim, NbRotations)
2
%
3
% ndim = 2 or 3 (but the code can be changed)
4
%
5
%
6
% code for generating an optimised sequence of rotations
7
% for the IDT pdf transfer algorithm
8
% although the code is not beautiful, it does the job.
9
%
10
function rotations = generate_rotations(ndim, NbRotations)
11

12
if (ndim == 2)
13
    l = [0 pi/2];
14
elseif (ndim == 3)
15
    l = [0 0 pi/2 0 pi/2 pi/2];
16
else % put here initialisation for higher orders
17
end
18

19
fprintf('rotation ');
20
for i = 1:(NbRotations-1)
21
    fprintf('%d ...', i );
22
    l = [l ; find_next(l, ndim)]; %l(end,:)+ones(1,ndim-1)*pi/2)]
23
    fprintf('\b\b\b', i );
24
end
25

26
M = ndim;
27

28
rotations = cell(1,NbRotations);
29
for i=1:size(l, 1)
30
    for j=1:M
31
        b_prev(j,:) = hyperspherical2cartesianT(l(i,(1:ndim-1) + (j-1)*(ndim-1)));
32
    end
33
    b_prev = grams(b_prev')';
34
    rotations{i} = b_prev;
35
end
36

37

38
end
39
%
40
%
41
function [x] = find_next( list_prev_x, ndim)
42

43
prevx = list_prev_x; % in hyperspherical coordinates
44
nprevx = size(prevx,1);
45
hdim = ndim - 1;
46
M = ndim;
47

48
% convert points to cartesian coordinates
49
c_prevx = zeros(nprevx*M, ndim);
50
c_prevx = [];
51
for i=1:nprevx
52
    for j=1:M
53
        b_prev(j,:) = hyperspherical2cartesianT(prevx(i,(1:hdim) + (j-1)*hdim));
54
    end
55
    b_prev = grams(b_prev')';
56
    c_prevx = [c_prevx; b_prev];
57
end
58

59
c_prevx;
60

61
options = optimset('TolX', 1e-10);
62
options = optimset(options,'Display','off');
63

64
minf = inf;
65
for i=1:10
66
    x0 = rand(1, hdim*M)*pi - pi/2;
67
    x = fminsearch(@myfun, x0, options);
68
    f = myfun(x);
69
    if f < minf
70
        minf = f;
71
        mix = x;
72
    end
73
end
74

75
%%
76
% f - Compute the function value at x
77
    function [f] = myfun(x)
78
        % compute the objective function
79
        c_x = zeros(M, ndim);
80
        for i=1:M
81
            c_x(i, :) = hyperspherical2cartesianT(x((1:hdim) + (i-1)*hdim));
82
        end
83
        c_x = grams(c_x')';
84
        f = 0;
85
        for i=1:M
86
            for p=1:size(c_prevx, 1)
87
                d = (c_prevx(p,:) - c_x(i, :)) * (c_prevx(p,:) - c_x(i, :))';
88
                f = f + 1/(1 + d);
89
                d = (c_prevx(p,:) + c_x(i, :)) * (c_prevx(p,:) + c_x(i, :))';
90
                f = f + 1/(1 + d);
91
            end
92
        end
93
    end
94
%%
95

96
end
97

98

99
%
100
%
101
function c = hyperspherical2cartesianT(x)
102

103
c = zeros(1, length(x)+1);
104
sk = 1;
105
for k=1:length(x)
106
    c(k) = sk*cos(x(k));
107
    sk = sk*sin(x(k));
108
end
109
c(end) = sk;
110

111
end
112

113
% Gram-Schmidt orthogonalization of the columns of A.
114
% The columns of A are assumed to be linearly independent.
115
function [Q, R] = grams(A)
116

117
[m, n] = size(A);
118
Asave = A;
119
for j = 1:n
120
    for k = 1:j-1
121
        mult = (A(:, j)'*A(:, k)) / (A(:, k)'*A(:, k));
122
        A(:, j) = A(:, j) - mult*A(:, k);
123
    end
124
end
125
for j = 1:n
126
    if norm(A(:, j)) < sqrt(eps)
127
        error('Columns of A are linearly dependent.')
128
    end
129
    Q(:, j) = A(:, j) / norm(A(:, j));
130
end
131
R = Q'*Asave;
132
end
133

134