Home > .. > 16-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 % if rob_p_d > 0
0048 %     gain=(sigma_n/sigma_e0)^(1/(Niter-1));
0049 % else
0050 %     gain = 1;
0051 % end
0052 gain=1;
0053 
0054 %% robust estimation
0055 for iter = 1:Niter+1
0056     if plot_type> 0
0057         figure
0058     end
0059     
0060     sigma_e = sigma_e0*gain^(Niter-iter);
0061     
0062     if print_type > 0
0063         iter__sigmae_sigma_n_ratio = [iter,sigma_e,sigma_n,sigma_e/sigma_n] %#ok<NOPRT,NASGU>
0064         weight_before_iteration = weights(1:20)' %#ok<NOPRT,NASGU>
0065     end
0066     
0067     b=zeros(M+N-2,1);
0068     % set A and b for observed points
0069     for m = 1:M
0070         w = 1/sigma_n*sqrt(weights(m));
0071         A(m,ts(m)) = w;
0072         dl      = ys(m)-xa(ts(m));
0073         b(m)    = dl*w;
0074     end
0075 %     if iter == Niter+1
0076 %         rhs_w=[b(1:M)';weights(1:M)'];
0077 %     end
0078     
0079     kk = M;
0080     % set A and b for regularizing observations
0081     for n = 2:N-1
0082         kk = kk+1;
0083         w = 1/sigma_e*sqrt(weights(kk));
0084         A(kk,n-1:n+1)= [1 -2 1]*w;
0085         dl           = -(xa(n+1)-2*xa(n)+xa(n-1));
0086         b(kk)        = dl*w;
0087     end
0088     
0089     % normal equation system
0090     Nmatrix = A'*A;
0091     nvector = A'*b;
0092     if N <= 1000
0093         Cov = inv(Nmatrix);
0094     end
0095     
0096     % solution of correction
0097     dx = Nmatrix\nvector;
0098     % corrected estimate
0099     xa = xa+dx;
0100 
0101     if iter == Niter+1 && print_type > 0
0102         xa_dx = [xa';dx'] %#ok<NOPRT,NASGU>
0103     end
0104     v = A*dx-b;                             % estimated standardized residuals
0105     
0106     eso   = norm(v)/sqrt(M-2);              % estimated sigma0
0107 
0108     vori  = (v*sigma_e./sqrt(weights));
0109     % vs    = conv2(abs(vori),[1,2,1]'/4,'same');
0110     vs    = vori;
0111     
0112     eso_p = median(abs(vs(1:M)))*1.48;       % estimated RMS of point residuals
0113     
0114     eso_k = median(abs(vs(M+1:end)))*1.48;  % estimated RMS curvature residuals
0115     
0116     if print_type > 0
0117         estimated_s0_s0p_s0k=[eso,eso_p,eso_k] %#ok<NOPRT,NASGU>
0118     end
0119 
0120     %kp = k*eso_p;
0121     kp = k;
0122     kk = k*eso_k;
0123 
0124 %% adapt weights
0125     if iter < Niter
0126         if type_out == 0                            % symmetric
0127             if print_type > 0
0128                 display('rob_p_d')
0129             end
0130             if iter < 4                             % first iterations
0131                 if rob_p_d == 1 || rob_p_d == 3     % dem
0132                     weights(M+1:M+N-2) = ...
0133                         min(1,1./(abs(vs(M+1:M+N-2))/kk+0.0001));
0134                     if print_type > 1
0135                         iter_kk_max_v_min_w = ...
0136                             [iter,kk,max(abs(vs(M+1:M+N-2))),...
0137                             min(weights(M+1:M+N-2))]; %#ok<NASGU>
0138                         display('iter_kk_max_v_min_w')
0139                         display('modify dem weights')
0140                         residuals_new_weights_of_dem=...
0141                             [vs(M+1:M+N-2)';weights(M+1:M+N-2)']; %#ok<NASGU>
0142                         display('residuals_new_weights_of_dem')
0143                     end
0144                     if plot_type > 0
0145                         subplot(3,1,1)
0146                         plot(2:N-1,xa(2:N-1),'-b')
0147                         title(num2str(sigma_e/sigma_n));
0148                         subplot(3,1,2);
0149                         hold on
0150                         plot(2:N-1,vs(M+1:M+N-2),'-b')
0151                         plot(2:N-1,vori(M+1:M+N-2),'-r')
0152                         plot(2:N-1,ones(N-2,1)*kk,'-b')
0153                         plot(2:N-1,-ones(N-2,1)*kk,'-b')
0154                         title(num2str(kk))
0155                         subplot(3,1,3);
0156                         plot(2:N-1,weights(M+1:M+N-2),'-b')
0157                         title(num2str(iter))
0158                     end
0159                 end
0160                 if rob_p_d == 2 || rob_p_d == 3     % points
0161                     weights(1:M) = min(1,1./(abs(vs(1:M))/kp+0.0001));
0162                     if print_type > 1
0163                         residuals_new_weights_of_points_L1 = ...
0164                             [(1:20);v(1:20)';vs(1:20)';weights(1:20)'] %#ok<NOPRT,NASGU>
0165                     end
0166                 end               
0167             else                                    % iteration 4 following
0168                 if rob_p_d == 1 || rob_p_d == 3     % dem
0169                     weights(M+1:M+N-2) = exp(-abs(vs(M+1:M+N-2)).^2./kk^2/2)+0.0001;
0170                     if print_type > 1
0171                         iter_kk_max_v_min_w = ...
0172                             [iter,kk,max(abs(vs(M+1:M+N-2))),min(weights(M+1:M+N-2))] %#ok<NOPRT,NASGU>
0173                         display('modify dem weights')
0174                         residuals_new_weights_of_dem = ...
0175                             [vs(M+1:M+N-2)';weights(M+1:M+N-2)'] %#ok<NOPRT,NASGU>
0176                     end
0177                     %
0178                     if plot_type > 0
0179                         subplot(3,1,1)
0180                         plot(2:N-1,xa(2:N-1),'-b')
0181                         title(num2str(sigma_e/sigma_n));
0182                         subplot(3,1,2);
0183                         hold on
0184                         plot(2:N-1,vs(M+1:M+N-2),'-b')
0185                         plot(2:N-1,vori(M+1:M+N-2),'-r')
0186                         plot(2:N-1,ones(N-2,1)*kk,'-b')
0187                         plot(2:N-1,-ones(N-2,1)*kk,'-b')
0188                         title(num2str(kk))
0189                         subplot(3,1,3);
0190                         plot(2:N-1,weights(M+1:M+N-2),'-b')
0191                         title(num2str(iter))
0192                     end
0193                 end
0194                 if rob_p_d == 2 || rob_p_d == 3         %points
0195                     weights(1:M) = exp(-abs(vs(1:M)).^2/kp^2/2)+0.0001;
0196                     if print_type > 1
0197                         residuals_new_weights_of_points_exp = ...
0198                             [(1:20);v(1:20)';vs(1:20)'/kp;weights(1:20)'] %#ok<NOPRT,NASGU>
0199                     end
0200                 end                
0201             end
0202         else  %  aymmetric
0203             k = 1;
0204             g = g_factor*k;
0205             w = w_factor*k;
0206             if type_robust==0   % asymmetric L1
0207                 weights(1:M) = max(min(1,1./abs(v(1:M))/k+0.0001),((sign(v(1:M))/k)+1)/2);
0208             else                % Kraus/Pfeifer
0209                 for m=1:M
0210                     if v(m) < g-w
0211                         weights(m) = 0;
0212                     elseif v(m) > g
0213                         weights(m)=1;
0214                     else
0215                         weights(m)=1/(1+(v(m)-g)^4);
0216                     end
0217                 end
0218             end
0219         end
0220     else
0221 
0222         weights = 0.9999*ceil(weights-0.1)+0.0001;
0223         % initiate last iteration (Niter+1) with approximate values 0
0224         if iter == Niter
0225             xa = zeros(N,1);
0226         end
0227     end
0228 end
0229 
0230 xest = xa;
0231 
0232 
0233

Generated on Sat 17-Sep-2016 20:58:55 by m2html © 2005