1 | %L2depths Compute depths of PRMM from basis L. |
---|
2 | % |
---|
3 | % No known depths are exploited. These can be different because of noise |
---|
4 | % etc. |
---|
5 | % |
---|
6 | % Parameters: |
---|
7 | % opt.verbose(1) .. whether display info |
---|
8 | |
---|
9 | function [Mdepths, lambda] = L2depths(L, M, Idepths, opt) |
---|
10 | |
---|
11 | if nargin < 4, opt = []; end |
---|
12 | if ~isfield(opt,'verbose') |
---|
13 | opt.verbose = 1; end |
---|
14 | |
---|
15 | if opt.verbose, fprintf('Computing depths...'); tic; end |
---|
16 | |
---|
17 | Mdepths = M; |
---|
18 | |
---|
19 | [m n] = size(M); m = m/3; |
---|
20 | lambda(m, n) = 0; % memory allocation |
---|
21 | |
---|
22 | for j = 1:n |
---|
23 | full = find(~isnan(M(1:3:end,j))); |
---|
24 | mis_rows = intersect(find(Idepths(:,j)==0),full); |
---|
25 | if length(mis_rows) > 0 |
---|
26 | submatrix = spread_depths_col(M(k2i(full),j),Idepths(full,j)); |
---|
27 | |
---|
28 | % We want submatrix to be in the space L -> |
---|
29 | % we search for combination of columns of the base of L i.e. |
---|
30 | % L(k2i(full),:)*res(1:4)-submatrix*[1 res(5:length(res))] = 0 |
---|
31 | right = submatrix(:,1); |
---|
32 | A = [ L(k2i(full),:) -(submatrix(:,2:size(submatrix,2))) ]; |
---|
33 | if rank(A) < size(A, 2) % depths cannot be computed => kill the data |
---|
34 | kill = full(~Idepths(full,j)); |
---|
35 | Mdepths(k2i(kill),j) = NaN; lambda(kill,j) = NaN; |
---|
36 | else |
---|
37 | res = A \ right; |
---|
38 | |
---|
39 | %test: er should be near to zero |
---|
40 | %er=L(k2i(full),:)*res(1:4)-submatrix*[1 res(5:length(res))']' |
---|
41 | |
---|
42 | % depth corresponding to right is/are set to 1 |
---|
43 | i = full(find(right(1:3:end))); lambda(i,j) = 1; |
---|
44 | Mdepths(k2i(i),j) = M(k2i(i),j); |
---|
45 | |
---|
46 | for ii = 1:size(submatrix,2)-1 |
---|
47 | i = full(find(submatrix(1:3:end,1+ii))); lambda(i,j) = res(4+ii); |
---|
48 | Mdepths(k2i(i),j) = M(k2i(i),j)*lambda(i,j); |
---|
49 | end |
---|
50 | end |
---|
51 | end |
---|
52 | end |
---|
53 | |
---|
54 | if opt.verbose, disp(['(' num2str(toc) ' sec)']); end |
---|