Home > 16-Surface-Reconstruction > Profilereconstruction > Functions > estimate_profile_robust.m

estimate_profile_robust

PURPOSE ^

% estimate_profile_robust

SYNOPSIS ^

function [xest,A,v,weights,Cov] = estimate_profile_robust(N,ts,ys,sigma_e,sigma_n,Niter,type_out,type_r,print_type,plot_type)

DESCRIPTION ^

% estimate_profile_robust

 xest = estimate_profile_robust(N,ts,ys,sigma_e,sigma_n,
                    Niter,type_out,type_r,print_type,plot_type)

 N        = scalar, length of profile
 ts,ys    = Mx1 vectors, observed profile
 sigma_e  = std of curvature
 sigma_n  = std of observations
 Niter    = number of iterations
          = 1 -> no robust estimation
 type_out = type of outlier (0=symm, 1=asymm)
 type_r   (1) > 0 = L1, 1 = Kraus 
          (2) = g for Kraus
          (3) = w for Kraus
          (4) = 1 only dem
                2 only points
                3 both
 plot_type    = 0,1,2 (no, low, high)
 print_type   = 0,1,2 (no, low, high)

 Wolfgang Förstner 2014-10-05
 last changes: Susanne Wenzel 09/16
 wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %% estimate_profile_robust
0002 %
0003 % xest = estimate_profile_robust(N,ts,ys,sigma_e,sigma_n,
0004 %                    Niter,type_out,type_r,print_type,plot_type)
0005 %
0006 % N        = scalar, length of profile
0007 % ts,ys    = Mx1 vectors, observed profile
0008 % sigma_e  = std of curvature
0009 % sigma_n  = std of observations
0010 % Niter    = number of iterations
0011 %          = 1 -> no robust estimation
0012 % type_out = type of outlier (0=symm, 1=asymm)
0013 % type_r   (1) > 0 = L1, 1 = Kraus
0014 %          (2) = g for Kraus
0015 %          (3) = w for Kraus
0016 %          (4) = 1 only dem
0017 %                2 only points
0018 %                3 both
0019 % plot_type    = 0,1,2 (no, low, high)
0020 % print_type   = 0,1,2 (no, low, high)
0021 %
0022 % Wolfgang Förstner 2014-10-05
0023 % last changes: Susanne Wenzel 09/16
0024 % wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de
0025 
0026 function [xest,A,v,weights,Cov] = estimate_profile_robust...
0027     (N,ts,ys,sigma_e,sigma_n,Niter,type_out,type_r,...
0028     print_type,plot_type)
0029 
0030 % robust - explore sensistivity of result wrt k! Choose k in [1:3]
0031 k = 2.0; 
0032 
0033 type_robust = type_r(1);
0034 g_factor    = type_r(2);
0035 w_factor    = type_r(3);
0036 rob_p_d     = type_r(4);
0037 
0038 % Number of observations
0039 M = length(ts);
0040 A = zeros(M+N-2,N);
0041 
0042 % initialize
0043 xa = zeros(N,1);
0044 weights = ones(M+N-2,1);
0045 Cov = 0;
0046 sigma_e0 = sigma_e;
0047 gain = 1;
0048 
0049 %% robust estimation
0050 for iter = 1:Niter+1
0051     if plot_type> 0
0052         figure
0053     end
0054     
0055     sigma_e = sigma_e0*gain^(Niter-iter);
0056     
0057     if print_type > 0
0058         iter__sigmae_sigma_n_ratio = [iter,sigma_e,sigma_n,sigma_e/sigma_n] %#ok<NOPRT,NASGU>
0059         weight_before_iteration = weights(1:20)' %#ok<NOPRT,NASGU>
0060     end
0061     
0062     b=zeros(M+N-2,1);
0063     % set A and b for observed points
0064     for m = 1:M
0065         w = 1/sigma_n*sqrt(weights(m));
0066         A(m,ts(m)) = w;
0067         dl      = ys(m)-xa(ts(m));
0068         b(m)    = dl*w;
0069     end
0070     
0071     kk = M;
0072     % set A and b for regularizing observations
0073     for n = 2:N-1
0074         kk = kk+1;
0075         w = 1/sigma_e*sqrt(weights(kk));
0076         A(kk,n-1:n+1)= [1 -2 1]*w;
0077         dl           = -(xa(n+1)-2*xa(n)+xa(n-1));
0078         b(kk)        = dl*w;
0079     end
0080     
0081     % normal equation system
0082     Nmatrix = A'*A;
0083     nvector = A'*b;
0084     if N <= 1000
0085         Cov = inv(Nmatrix);
0086     end
0087     
0088     % solution of correction
0089     dx = Nmatrix\nvector;
0090     % corrected estimate
0091     xa = xa+dx;
0092 
0093     if iter == Niter+1 && print_type > 0
0094         xa_dx = [xa';dx'] %#ok<NOPRT,NASGU>
0095     end
0096     v = A*dx-b;                             % estimated standardized residuals
0097     
0098     eso   = norm(v)/sqrt(M-2);              % estimated sigma0
0099 
0100     vori  = (v*sigma_e./sqrt(weights));
0101     vs    = vori;
0102     
0103     eso_p = median(abs(vs(1:M)))*1.48;       % estimated RMS of point residuals
0104     
0105     eso_k = median(abs(vs(M+1:end)))*1.48;  % estimated RMS curvature residuals
0106     
0107     if print_type > 0
0108         estimated_s0_s0p_s0k=[eso,eso_p,eso_k] %#ok<NOPRT,NASGU>
0109     end
0110 
0111     kp = k;
0112     kk = k*eso_k;
0113 
0114 %% adapt weights
0115     if iter < Niter
0116         if type_out == 0                            % symmetric
0117             if print_type > 0
0118                 display('rob_p_d')
0119             end
0120             if iter < 4                             % first iterations
0121                 if rob_p_d == 1 || rob_p_d == 3     % dem
0122                     weights(M+1:M+N-2) = ...
0123                         min(1,1./(abs(vs(M+1:M+N-2))/kk+0.0001));
0124                     if print_type > 1
0125                         iter_kk_max_v_min_w = ...
0126                             [iter,kk,max(abs(vs(M+1:M+N-2))),...
0127                             min(weights(M+1:M+N-2))]; %#ok<NASGU>
0128                         display('iter_kk_max_v_min_w')
0129                         display('modify dem weights')
0130                         residuals_new_weights_of_dem=...
0131                             [vs(M+1:M+N-2)';weights(M+1:M+N-2)']; %#ok<NASGU>
0132                         display('residuals_new_weights_of_dem')
0133                     end
0134                     if plot_type > 0
0135                         subplot(3,1,1)
0136                         plot(2:N-1,xa(2:N-1),'-b')
0137                         title(num2str(sigma_e/sigma_n));
0138                         subplot(3,1,2);
0139                         hold on
0140                         plot(2:N-1,vs(M+1:M+N-2),'-b')
0141                         plot(2:N-1,vori(M+1:M+N-2),'-r')
0142                         plot(2:N-1,ones(N-2,1)*kk,'-b')
0143                         plot(2:N-1,-ones(N-2,1)*kk,'-b')
0144                         title(num2str(kk))
0145                         subplot(3,1,3);
0146                         plot(2:N-1,weights(M+1:M+N-2),'-b')
0147                         title(num2str(iter))
0148                     end
0149                 end
0150                 if rob_p_d == 2 || rob_p_d == 3     % points
0151                     weights(1:M) = min(1,1./(abs(vs(1:M))/kp+0.0001));
0152                     if print_type > 1
0153                         residuals_new_weights_of_points_L1 = ...
0154                             [(1:20);v(1:20)';vs(1:20)';weights(1:20)'] %#ok<NOPRT,NASGU>
0155                     end
0156                 end               
0157             else                                    % iteration 4 following
0158                 if rob_p_d == 1 || rob_p_d == 3     % dem
0159                     weights(M+1:M+N-2) = exp(-abs(vs(M+1:M+N-2)).^2./kk^2/2)+0.0001;
0160                     if print_type > 1
0161                         iter_kk_max_v_min_w = ...
0162                             [iter,kk,max(abs(vs(M+1:M+N-2))),min(weights(M+1:M+N-2))] %#ok<NOPRT,NASGU>
0163                         display('modify dem weights')
0164                         residuals_new_weights_of_dem = ...
0165                             [vs(M+1:M+N-2)';weights(M+1:M+N-2)'] %#ok<NOPRT,NASGU>
0166                     end
0167                     %
0168                     if plot_type > 0
0169                         subplot(3,1,1)
0170                         plot(2:N-1,xa(2:N-1),'-b')
0171                         title(num2str(sigma_e/sigma_n));
0172                         subplot(3,1,2);
0173                         hold on
0174                         plot(2:N-1,vs(M+1:M+N-2),'-b')
0175                         plot(2:N-1,vori(M+1:M+N-2),'-r')
0176                         plot(2:N-1,ones(N-2,1)*kk,'-b')
0177                         plot(2:N-1,-ones(N-2,1)*kk,'-b')
0178                         title(num2str(kk))
0179                         subplot(3,1,3);
0180                         plot(2:N-1,weights(M+1:M+N-2),'-b')
0181                         title(num2str(iter))
0182                     end
0183                 end
0184                 if rob_p_d == 2 || rob_p_d == 3         %points
0185                     weights(1:M) = exp(-abs(vs(1:M)).^2/kp^2/2)+0.0001;
0186                     if print_type > 1
0187                         residuals_new_weights_of_points_exp = ...
0188                             [(1:20);v(1:20)';vs(1:20)'/kp;weights(1:20)'] %#ok<NOPRT,NASGU>
0189                     end
0190                 end                
0191             end
0192         else  %  aymmetric
0193             k = 1;
0194             g = g_factor*k;
0195             w = w_factor*k;
0196             if type_robust==0   % asymmetric L1
0197                 weights(1:M) = max(min(1,1./abs(v(1:M))/k+0.0001),((sign(v(1:M))/k)+1)/2);
0198             else                % Kraus/Pfeifer
0199                 for m=1:M
0200                     if v(m) < g-w
0201                         weights(m) = 0;
0202                     elseif v(m) > g
0203                         weights(m)=1;
0204                     else
0205                         weights(m)=1/(1+(v(m)-g)^4);
0206                     end
0207                 end
0208             end
0209         end
0210     else
0211 
0212         weights = 0.9999*ceil(weights-0.1)+0.0001;
0213         % initiate last iteration (Niter+1) with approximate values 0
0214         if iter == Niter
0215             xa = zeros(N,1);
0216         end
0217     end
0218 end
0219 
0220 xest = xa;
0221 
0222 
0223

Generated on Mon 08-Jan-2018 17:21:49 by m2html © 2005