1 | %balance_triplets Balance PRMM by column-wise and "triplet-of-rows"-wise |
---|
2 | % scalar multiplications. |
---|
3 | % |
---|
4 | % After balancing, overall weight of M will be m*n where 3m*n is size of M |
---|
5 | % (i.e. 3 coordinates of each image point will give together 1 in average). |
---|
6 | % |
---|
7 | % Parameters: |
---|
8 | % opt.info_separately .. 1(default) .. display info on separate row |
---|
9 | % 0 .. display info in brackets |
---|
10 | % .verbose(1) .. whether display info |
---|
11 | |
---|
12 | function B = balance_triplets(M, opt); |
---|
13 | |
---|
14 | if nargin < 2, opt = []; end |
---|
15 | if ~isfield(opt,'info_separately') |
---|
16 | opt.info_separately = 1; end |
---|
17 | if ~isfield(opt,'verbose') |
---|
18 | opt.verbose = 1; end |
---|
19 | |
---|
20 | if opt.verbose |
---|
21 | if opt.info_separately, fprintf(1,'Balancing PRMM...'); tic |
---|
22 | else, fprintf(1,'(balancing PRMM...'); tic; end; end |
---|
23 | |
---|
24 | m=size(M,1)/3; %number of cameras |
---|
25 | n=size(M,2); %number of points |
---|
26 | |
---|
27 | B=M; change=inf; diff_rows = inf; iteration = 0; |
---|
28 | while (change > 0.01 | diff_rows > 1 | diff_cols > 1) & iteration <= 20 |
---|
29 | Bold=B; |
---|
30 | |
---|
31 | % 1. rescale each column l so that sum w_r_l^2=1 i.e. column of unit |
---|
32 | % length. However, due to the missing data, the length must be smaller |
---|
33 | % in (linear) dependance on amount of missing data. |
---|
34 | diff_cols = -inf; |
---|
35 | for l=1:n |
---|
36 | rows = find(~isnan(M(1:3:end,l))); |
---|
37 | if length(rows) > 0 |
---|
38 | rowsb = k2i(rows); |
---|
39 | s = sum(B(rowsb,l) .^ 2); |
---|
40 | supposed_weight = length(rows); % the less data, the slighter impact |
---|
41 | diff_cols = max(abs(s-supposed_weight), diff_cols); |
---|
42 | B(rowsb,l) = B(rowsb,l) ./ sqrt(s/supposed_weight); |
---|
43 | end |
---|
44 | end |
---|
45 | |
---|
46 | % 2. rescale each triplet of rows so that it's sum w_i_l^2=1 i.e. unit |
---|
47 | % area. However, due to the missing data, the area must be smaller in |
---|
48 | % (linear) dependance on amount of missing data. |
---|
49 | diff_rows = -inf; |
---|
50 | for k=1:m |
---|
51 | cols = find(~isnan(M(3*k,:))); |
---|
52 | if length(cols) > 0 |
---|
53 | s = sum(sum( B(3*k-2:3*k,cols) .^ 2 )); |
---|
54 | supposed_weight = length(cols); % the less data, the slighter impact |
---|
55 | diff_rows = max(abs(s-supposed_weight), diff_rows); |
---|
56 | B(3*k-2:3*k,cols) = B(3*k-2:3*k,cols) ./ sqrt(s/supposed_weight); |
---|
57 | end |
---|
58 | end |
---|
59 | |
---|
60 | % repeat steps 1 and 2 if W changed significantly |
---|
61 | % Note: It is not ensured that sums (s) would not change significantly |
---|
62 | % in the (hypothetical) next step. The reason is that rescaling |
---|
63 | % No. 1 rescales n columns to overall weight n whereas rescaling |
---|
64 | % No. 2 rescales m triplets of rows to overall weight m. |
---|
65 | change = 0; |
---|
66 | for k=1:m |
---|
67 | cols = find(~isnan(M(3*k,:))); |
---|
68 | change = change + sum(sum( (B (3*k-2:3*k,cols)-... |
---|
69 | Bold(3*k-2:3*k,cols)) .^ 2 )); |
---|
70 | end |
---|
71 | %disp(sprintf('change=%f, diff_cols=%f, diff_rows=%f', ... |
---|
72 | % change, diff_cols, diff_rows)); |
---|
73 | iteration = iteration +1; |
---|
74 | if opt.verbose, fprintf(1,'.'); end |
---|
75 | end |
---|
76 | |
---|
77 | |
---|
78 | if opt.verbose |
---|
79 | if opt.info_separately, disp(['(' num2str(toc) ' sec)']); |
---|
80 | else, fprintf(1,[num2str(toc) ' sec)']); end; end |
---|
81 | |
---|
82 | %disp(sprintf('change: %f, diff_rows: %f, diff_cols: %f, iterations: %d', ... |
---|
83 | % change, diff_rows, diff_cols, iteration)); |
---|