Home > 16-Surfacereconstruction > fig_16_21_precision_theoretical.m

fig_16_21_precision_theoretical

PURPOSE ^

% Fig. 15. 20b theoretical quality

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

% Fig. 15. 20b theoretical quality

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %% Fig. 15. 20b theoretical quality
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 
0013 close all
0014 
0015 
0016 display('----- theoretical quality of surface reconstruction -----')
0017 init_rand = 2    % Example: with intit_rand = 2
0018 N=25             % Example: with N = 25
0019 sigma_h=1        % Example: with 1
0020 sigma_k=1        % Example: with 5
0021 Resolution=50    % Example: with 50
0022 
0023 addpath(genpath('../General_Functions/'));
0024 addpath('Functions')
0025 
0026 type_data = 0;
0027 type_robust = 0;
0028 %             0 not robust
0029 %             1 only dem
0030 %             2 only points
0031 %             3 both
0032 
0033 out_C = 1;
0034 out_print = 0;
0035 Tol =0.15
0036 
0037 print_type = 0;
0038 plot_type  = 0;
0039 
0040 if init_rand==0
0041     init_rand=round(rand( sum(100*clock)))
0042     r=rand(3,1);
0043 end;
0044 rand('state',init_rand);
0045 randn('state',init_rand);
0046 
0047 
0048 %% generate dem point cloud
0049 
0050 [points,BB,delta_x,dem,bi,dt]=simulate_points_dem_16_precision(N,sigma_h,Resolution);
0051 out_in=zeros(size(points,1),1);
0052 Np = size(points,1);
0053 
0054 
0055 %% interpolate
0056 
0057 starttime = cputime;
0058 
0059 [ds,S,Sigma,Np,Nr,Mc,ver,A,w,w_f,W] =...
0060     smooth_dem_robust_bilinear...
0061     (points,BB,delta_x,sigma_k,out_C,type_robust,type_data,out_in,...
0062     print_type,plot_type);
0063 
0064 complete_time_for_solution=cputime-starttime
0065 
0066 %% sensitivity factors
0067 
0068 rk=zeros(Np,1);
0069 mu=zeros(Np,1);
0070 for k=1:Np
0071     i=points(k,1)+1;
0072     j=points(k,2)+1;
0073     vk = A(k,:)*Sigma*A(k,:)';
0074     rk(k) = (sigma_h^2-vk)/sigma_h^2;
0075     mu(k)=sqrt(vk/(1-vk));
0076 end
0077 
0078 %%
0079 X=([0:Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0080 Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0081 %%
0082 
0083 figure
0084 hold on
0085 %subplot(2,2,1)
0086 
0087 hold on
0088 surf(X,Y,ds); %,'EdgeColor',[0,0,0])% ,129*ones(size(S))
0089 colormap(gray)
0090 
0091 title('fitted dem: z')
0092 alpha(0.3)
0093 col=zeros(Np,2);
0094 for n=1:Np
0095     plot3(points(n,1),points(n,2),points(n,3),'.k','MarkerSize',15)
0096 end
0097 view([-29,65])
0098 xlim([BB(1),BB(3)])
0099 ylim([BB(2),BB(4)])
0100 zlim([0,10])
0101 axis equal
0102 
0103 %subplot(2,2,2)
0104 figure
0105 hold on
0106 Smax=max(sqrt(S(:)));
0107 surf(X,Y,sqrt(S)/Smax*10,'EdgeColor',[0,0,0])
0108 colormap(cool);
0109 title('standard deviations \sigma_z - sensitifyty factors \mu')
0110 col=zeros(Np,2);
0111 alpha(0.3)
0112 z=zeros(Np,1);
0113 for n=1:Np
0114     z(n)=interpolate_bilinear(S,points(n,1),points(n,2),delta_x,BB,sigma_h);
0115     plot3(points(n,1),points(n,2),sqrt(z(n)),'.k','MarkerSize',15) 
0116     %plot3([points(n,1),points(n,1)],[points(n,2),points(n,2)],[0,mu(n)+1],'-w','LineWidth',4)
0117     plot3([points(n,1),points(n,1)],[points(n,2),points(n,2)],[0,mu(n)],'-r','LineWidth',2) 
0118 end
0119 plot3(0,0,0,'.k')
0120 view([-29,65])
0121 xlim([BB(1),BB(3)])
0122 ylim([BB(2),BB(4)])
0123 zlim([-1,7])
0124 axis equal
0125 
0126 x_y_sigma_std_fitted_points = [points(:,[1,2,4]),sqrt(z)];
0127 
0128 Std_z=sqrt(S);
0129 
0130 %% plot sensitivity factors
0131 figure
0132 hold on
0133 R=Resolution;
0134 plot([0,R+5,R+5,0,0],[0,0,R+5,R+5,0],'-k')
0135 for k=1:Np
0136     i=points(k,1)+1;
0137     j=points(k,2)+1;
0138     plot(i,j,'.b','MarkerSize',10)
0139     text(i+0.5,j+0.5,num2str(round(mu(k)*10.001)/10,'%6.1f'))
0140 end
0141 title('sensitivity factors')
0142 axis equal
0143 axis off
0144 
0145 display('             sensitivity                ')
0146 display('         i         j        rk        mu')
0147 [points(:,1:2),rk,mu]
0148 average_redundancy_observations=sum(rk)/Np
0149 
0150 
0151      
0152 return

Generated on Sat 01-Oct-2016 19:56:24 by m2html © 2005