0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013 close all
0014
0015
0016 display('----- theoretical quality of surface reconstruction -----')
0017 init_rand = 2
0018 N=25
0019 sigma_h=1
0020 sigma_k=1
0021 Resolution=50
0022
0023 addpath(genpath('../General_Functions/'));
0024 addpath('Functions')
0025
0026 type_data = 0;
0027 type_robust = 0;
0028
0029
0030
0031
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
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
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
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
0086
0087 hold on
0088 surf(X,Y,ds);
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
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
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
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