0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039 function [dem_smoothed,S,Sigma,Np,Nr,Mc,v,A,weights,W] = ...
0040 smooth_dem_robust_bilinear_flat...
0041 (points,BB,delta_x,sigma_s,out_C,type_robust,type_data,...
0042 out_in,print_type,plot_type);
0043
0044
0045 [Npp,dim] = size(points);
0046
0047 if print_type > 0
0048 number_of_points_given = Npp
0049 end
0050
0051
0052 xmin = BB(1);
0053 ymin = BB(2);
0054 xmax = BB(3);
0055 ymax = BB(4);
0056 Nr =ceil((xmax-xmin)/delta_x)+1;
0057 Mc =ceil((ymax-ymin)/delta_x)+1;
0058
0059
0060 Np = 0;
0061 for p=1:Npp
0062 if points(p,1) >= xmin && ...
0063 points(p,1) < xmax && ...
0064 points(p,2) >= ymin && ...
0065 points(p,2) < ymax
0066 Np=Np+1;
0067 poi(Np,:) = [(points(p,1)-xmin)/delta_x,(points(p,2)-ymin)/delta_x,points(p,3:4)];
0068 end
0069 end
0070
0071 if print_type > 1
0072 number_of_points_in_BB = Np
0073 end
0074
0075
0076
0077
0078 U = Nr*Mc;
0079 N = Np + (Nr-1)*Mc + Nr*(Mc-1);
0080 if print_type > 0
0081 number_of_unknowns = U
0082 end
0083
0084
0085 if U > 1600
0086 out_C = 0;
0087 end
0088 if print_type > 0
0089 tic
0090 end
0091
0092
0093 iter_switch=3;
0094 if type_robust == 0
0095 N_iter = 1;
0096 first = 1;
0097 else
0098 N_iter=6;
0099 first=1;
0100 end
0101 xa = zeros(U,1);
0102 weights = ones(N,1);
0103 sigma_s0=sigma_s;
0104 for iter = 1:N_iter
0105 if print_type > 0
0106 display('Aufruf von estimate_dem_robust_flat')
0107 start_ex=[iter,type_robust]
0108 end
0109 L1=0;
0110 if type_robust >0 && iter <= iter_switch
0111 L1=1;
0112 sigma_s=sigma_s0*4;
0113 else
0114 L1=0;
0115 sigma_s=sigma_s0;
0116 end
0117 sigma_s=sigma_s0*1.4^(N_iter-iter);
0118 [A,weights,dx,v,W]=estimate_dem_robust_flat...
0119 (N,U,Np,Nr,Mc,poi,xa,weights,sigma_s,iter,type_robust,L1,...
0120 out_in,print_type,plot_type);
0121 first=0;
0122 xa = xa + dx;
0123
0124 if print_type > 0
0125 iter_dx = [iter,norm(dx)]
0126 end
0127
0128 if plot_type > 0
0129 if iter == 1
0130 figure(20)
0131 hold on
0132 X=([4*Nr:4*Nr+Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0133 Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0134
0135 if type_data ~=5 && type_robust > 0
0136 colormap(gray)
0137 surf(X,Y,reshape(xa,Nr,Mc));
0138 alpha(0.3);
0139
0140
0141 plot3(points(:,1)+4*Nr*delta_x,points(:,2),points(:,3),'.r','MarkerSize',20);
0142
0143 else
0144 subplot(1,3,2)
0145
0146 imshow(reshape(xa,Nr,Mc));
0147 colormap(gray)
0148 view(0,90)
0149 title('1.iteration')
0150 end
0151 title(strcat('iteration=',num2str(iter),' fitted dem'))
0152
0153 axis equal
0154 else
0155 figure
0156 hold on
0157 X=([0:Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0158 Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0159 colormap(gray)
0160 surf(X,Y,reshape(xa,Nr,Mc));
0161 alpha(0.3);
0162
0163 plot3(points(:,1),points(:,2),points(:,3),'.r','MarkerSize',20);
0164
0165 title(strcat('iteration=',num2str(iter),' fitted dem'))
0166
0167 if type_data ~= 5
0168 axis equal
0169 end
0170 end
0171 end
0172
0173 end
0174
0175
0176
0177
0178
0179
0180
0181
0182
0183
0184 dem_smoothed = reshape(xa,Nr,Mc);
0185
0186
0187 Sigma=0;
0188 S = 0;
0189 if print_type > 0
0190 tic
0191 end
0192 if out_C > 0
0193 Sigma = full(inv(A'*A));
0194 S = full(reshape(diag(Sigma),Nr,Mc));
0195 end
0196 if print_type > 0
0197 time_for_inverse=toc
0198 end
0199
0200
0201 return