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