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