fig_16_21_quality_surface_reconstruction

% Fig. 16.21 theoretical quality  of surface reconstruction

 Precision as a function of position
 Sensitivity wrt outliers in given heights

 wf 6/2015

 type_data: type of generated data
 0     a set of given points

0001 %% Fig. 16.21 theoretical quality  of surface reconstruction
0002 %
0003 % Precision as a function of position
0004 % Sensitivity wrt outliers in given heights
0005 %
0006 % wf 6/2015
0007 %
0008 % type_data: type of generated data
0009 % 0     a set of given points
0010 
0011 
0012 close all
0013 % clear all
0014 clearvars
0015 
0016 addpath(genpath('../../General-Functions/'));
0017 addpath('Functions')
0018 addpath('Data')
0019 
0020 ss = plot_init;
0021 
0022 
0023 
0024 disp('---------------------------------------------------------')
0025 disp('----- theoretical quality of surface reconstruction -----')
0026 disp('---------------------------------------------------------')
0027 
0028 %% set parameters
0029 init_rand    =   2     % Example: with intit_rand = 2
0030 N            =  25     % Example: with N = 25
0031 sigma_h      =   1     % Example: with 1
0032 sigma_k0     = 250     % Example: with 250
0033 Resolution   =  50     % Example: with 50
0034 sigma_k      = sigma_k0/Resolution   % adapted for resolution
0035 
0036 out_sensitivity = 1  % = 0 if number of observations too large
0037 
0038 type_robust = 0;
0039 %             0 not robust
0040 %             1 only dem
0041 %             2 only points
0042 %             3 both
0043 
0044 out_C = 1;
0045 out_print = 0;
0046 Tol =0.15
0047 
0048 print_type = 1;
0049 plot_type  = 0;
0050 
0051 init_rand_seed(init_rand);
0052 
0053 
0054 %% generate dem point cloud
0055 
0056 [points,BB,delta_x,dem,bi,dt]=simulate_points_dem_16_precision(N,sigma_h,Resolution);
0057 out_in=ones(size(points,1),1);
0058 Np = size(points,1);
0059 
0060 
0061 %% interpolate
0062 
0063 starttime = cputime;
0064 
0065 [ds,S,Sigma,Np,Nr,Mc,ver,A,w,w_f,W] =...
0066     smooth_dem_robust_bilinear...
0067     (points,BB,delta_x,sigma_k,out_C,type_robust,out_in,...
0068     print_type,plot_type);
0069 
0070 disp('---------------------------------------------------------')
0071 complete_time_for_solution=cputime-starttime
0072 
0073 %% sensitivity factors
0074 rk=zeros(Np,1);
0075 mu=zeros(Np,1);
0076 for k=1:Np
0077     ak = A(k,:)';                          % k-th row of Jacobian
0078     sigma_k = ak'*Sigma*ak;                % standard deviation of k-th fitted height
0079     rk(k) = (1-sigma_k)/sigma_h^2;         % redundancy number
0080     mu(k) = sqrt(sigma_k/(1-sigma_k));     % sensitivity factor
0081 end
0082 
0083 
0084 %%
0085 X=([0:Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0086 Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0087 %% plot fitted dem
0088 
0089 figure('Color','w','Position',[20,0.2*ss(2),0.6*ss(1),0.6*ss(2)])
0090 hold on
0091 %subplot(2,2,1)
0092 
0093 hold on
0094 surf(X,Y,ds); %,'EdgeColor',[0,0,0])% ,129*ones(size(S))
0095 colormap(gray)
0096 
0097 title('fitted dem: z')
0098 alpha(0.3)
0099 col=zeros(Np,2);
0100 for n=1:Np
0101     plot3(points(n,1),points(n,2),points(n,3),'.k','MarkerSize',15)
0102 end
0103 view([-29,65])
0104 xlim([BB(1),BB(3)])
0105 ylim([BB(2),BB(4)])
0106 zlim([0,10])
0107 axis equal
0108 
0109 %% plot theoretical quality
0110 figure('Color','w','Position',[30,0.15*ss(2),0.6*ss(1),0.6*ss(2)])
0111 hold on
0112 Smax=max(sqrt(S(:)));
0113 surf(X,Y,sqrt(S)/Smax*10,'EdgeColor',[0,0,0])
0114 colormap(cool);
0115 title('standard deviations $\sigma_z$ - sensitifyty factors $\mu$')
0116 col=zeros(Np,2);
0117 alpha(0.3)
0118 z=zeros(Np,1);
0119 for n=1:Np
0120     z(n)=interpolate_bilinear(S,points(n,1),points(n,2),delta_x,BB,sigma_h);
0121     plot3(points(n,1),points(n,2),sqrt(z(n)),'.k','MarkerSize',15)
0122 end
0123 if out_sensitivity ~= 0
0124     for n=1:Np
0125         plot3([points(n,1),points(n,1)],[points(n,2),points(n,2)],[0,mu(n)],'-r','LineWidth',2)
0126     end
0127 end
0128 plot3(0,0,0,'.k')
0129 view([-29,65])
0130 xlim([BB(1),BB(3)])
0131 ylim([BB(2),BB(4)])
0132 zlim([-1,7])
0133 axis equal
0134 
0135 x_y_sigma_std_fitted_points = [points(:,[1,2,4]),sqrt(z)];
0136 
0137 Std_z=sqrt(S);
0138 
0139 %% plot sensitivity factors
0140 
0141 if out_sensitivity ~= 0
0142     figure('Color','w','Position',[40,0.20*ss(2),0.5*ss(1),0.5*ss(2)])
0143     hold on
0144     R=Resolution;
0145     plot([0,R+5,R+5,0,0],[0,0,R+5,R+5,0],'-k')
0146     for k=1:Np
0147         i=points(k,1)+1;
0148         j=points(k,2)+1;
0149         plot(i,j,'.b','MarkerSize',10)
0150         text(i+0.5,j+0.5,num2str(round(mu(k)*10.001)/10,'%6.1f'))
0151     end
0152     disp('---------------------------------------------------------')
0153     title('sensitivity factors')
0154     axis equal
0155     axis off
0156 
0157     disp('             sensitivity                ')
0158     disp('         i         j        rk        mu')
0159     [points(:,1:2),rk,mu]
0160 end
0161 average_redundancy_observations=sum(rk)/Np
0162 
0163 
0164 return