[37] | 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)); |
---|