fig_16_18_step_break

% Fig. 15.18 page 758
 Robust estimation for tackling steps and break points

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

0001 %% Fig. 15.18 page 758
0002 % Robust estimation for tackling steps and break points
0003 %
0004 % Wolfgang Förstner 2014-08-06
0005 % last changes: Susanne Wenzel 09/16
0006 % wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de
0007 
0008 
0009 
0010 close all
0011 clearvars
0012  
0013 addpath(genpath('../../General-Functions/'));
0014 addpath('Functions')
0015  
0016 
0017 %% plot settings
0018 ss = plot_init;
0019 
0020 %% set parameters
0021 
0022 % params for generating the data
0023 
0024 init_rand   = 6;                % just changes the noise, may be changed for other example
0025 Nd          = 20;               % third of length of profile
0026 N           = 3*Nd;             % total number of points
0027 dens        = 1.0;              % percentage of observations
0028 sigma_e     = 0.15;             % process noise
0029 sigma_n     = 0.9;              % observation noise
0030 Poutlier    = 0.00;             % P(outlier) - fraction of outliers
0031 Max_outlier = 25;               % max(outlier)
0032 factor = 1.0;                   % vertical scaling of true signal
0033                                 % (for changing signal to noise ratio)
0034 type_outlier = 0;               % type of outlier (0=symm, 1=asymm)
0035 
0036 % params for estimation
0037 type_robust  = [0,...           %  [0 = L1, 1 = Kraus,g_factor,w_factor]
0038                 2,...           %  g
0039                 2,...           %  w
0040                 1];             %  [0: not robust,
0041                                 %   1: only points,
0042                                 %   2: only dem]
0043                         
0044 print_type = 0;
0045 plot_type  = 0;
0046                         
0047 if type_outlier ==0
0048     Niter = 1;
0049     if type_robust(4) > 0
0050         Niter = 6;
0051     end
0052     factor_sigma = 1;
0053 else
0054     Niter = 3;
0055     if type_robust(1) == 0
0056         factor_sigma= 8;
0057     else
0058         factor_sigma= 1;
0059     end
0060 end
0061 
0062 
0063 %% initialize random number generation by fixed seed
0064 init_rand_seed(init_rand);
0065 
0066 %% generate profile
0067 [x,y,out_in,select,xs,ys] = ...
0068     generate_step_corner(Nd,sigma_e,sigma_n,Poutlier,Max_outlier,type_outlier,dens,factor);
0069 
0070 %% reconstruct profile non robust
0071 type_robust(4) = 0;   % reconstruction non robust
0072 [xest_nonrob,~,~,weigths_nonrob] = estimate_profile_robust...
0073     (N,select,ys,sigma_e/factor_sigma,sigma_n,Niter,type_outlier,...
0074     type_robust,print_type,plot_type);
0075 
0076 %% reconstruct profile robust
0077 type_robust(4) = 1;     % reconstruction robust, only points
0078 [xest_robust,~,v,weights] = estimate_profile_robust...
0079     (N,select,ys,sigma_e/factor_sigma,sigma_n,Niter,type_outlier,...
0080     type_robust,print_type,plot_type);
0081 
0082 %% plot
0083 figure('name','Fig. 16.18: Robust estimation for tackling steps and break points',...
0084     'color','w', 'Position',[0.1*ss(1),0.2*ss(2),0.5*ss(1),0.7*ss(2)]);
0085 
0086 subplot(2,1,1);hold on
0087 plot(1:N,x,'--r','LineWidth',2)
0088 plot(1:N,xest_nonrob,'-k','LineWidth',2)
0089 plot(1:N,ones(N,1),'-k')
0090 plot(select,ys,'.b','MarkerSize',15)
0091 plot(2:N-1,weigths_nonrob(2:N-1)+1,'.k','Markersize',15)
0092 ylim([0,N/3+8]); xlim([0,N+1]); axis equal
0093 xlabel('$x$'); ylabel('$z$');
0094 title('Fig. 16.18a: Without robust estimation');
0095 
0096 
0097 subplot(2,1,2);hold on
0098 plot(1:N,x,'--r','LineWidth',2)
0099 plot(1:N,xest_robust,'-k','LineWidth',2)
0100 plot(1:N,ones(N,1),'-k')
0101 plot(select,ys,'.b','MarkerSize',15)
0102 plot(2:N-1,weights(end-(N-3):end)+1,'.k','Markersize',15)
0103 xlim([0,N+1]); ylim([0,N/3+8]); axis equal
0104 xlabel('$x$'); ylabel('$z$');
0105 title({'Fig. 15.18b: With robust estimation.'...
0106     ['Iterations $1$ to $3$ use the $L_{12}$-norm,'...
0107     ' iterations $4$ to $6$ use the exponential weight function']});
0108 
0109 
0110 
0111 
0112 
0113