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
0040
0041
0042
0043
0044
0045 function [dem_smoothed,S,Sigma,Np,Nr,Mc,v,A,weights,weights_f,W] = ...
0046 smooth_dem_robust_bilinear...
0047 (points,BB,delta_x,sigma_k,out_C,type_robust,type_data,...
0048 out_in,print_type,plot_type)
0049
0050
0051
0052 [Npp,tmp] = size(points);
0053
0054 xmin = BB(1);
0055 ymin = BB(2);
0056 xmax = BB(3);
0057 ymax = BB(4);
0058 Nr =ceil((xmax-xmin)/delta_x)+1;
0059 Mc =ceil((ymax-ymin)/delta_x)+1;
0060
0061
0062 Np = 0;
0063 for p=1:Npp
0064 if points(p,1) >= xmin && ...
0065 points(p,1) < xmax && ...
0066 points(p,2) >= ymin && ...
0067 points(p,2) < ymax
0068 Np=Np+1;
0069 poi(Np,:) = [(points(p,1)-xmin)/delta_x,(points(p,2)-ymin)/delta_x,points(p,3:4)];
0070
0071 end
0072 end
0073
0074
0075 U = Nr*Mc;
0076 N = Np + (Nr-2)*Mc + Nr*(Mc-2) + (Nr-1)*(Mc-1);
0077 if print_type > 0
0078 number_of_unknowns = U
0079 end
0080
0081
0082 if U > 40401
0083 out_C = 0;
0084 end
0085 if print_type > 0
0086 tic
0087 end
0088
0089
0090 iter_switch=3;
0091 if type_robust == 0
0092 N_iter = 1;
0093 else
0094 N_iter=6;
0095 end
0096
0097 xa = zeros(U,1);
0098 weights = ones(N,1);
0099 weights_f = weights;
0100 sigma_k0 = sigma_k;
0101 gain=1.0;
0102
0103 for iter = 1:N_iter
0104 L1 = 0;
0105 if type_robust > 0 && iter <= iter_switch
0106 L1=1;
0107 else
0108 L1=0;
0109 end
0110 sigma_k=sigma_k0*gain^(N_iter-iter);
0111
0112 [A,weights,dx,v,W,Nmatrix]=estimate_dem_robust...
0113 (N,U,Np,Nr,Mc,poi,xa,weights,sigma_k,iter,type_robust,L1,...
0114 out_in,print_type,plot_type);
0115
0116 xa = xa + dx;
0117
0118 if print_type > 0
0119 iter_dx = [iter,norm(dx)]
0120 end
0121
0122 if plot_type > 0
0123 if iter == 1
0124 if type_robust > 0
0125 figure
0126 hold on
0127 title('surf')
0128 X=([4*Nr:4*Nr+Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0129 Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0130
0131 if type_data ~=5
0132 colormap(gray)
0133 surf(X,Y,reshape(xa,Nr,Mc));
0134 alpha(0.3);
0135
0136 view([-29,65])
0137 plot3(points(:,1)+4*Nr*delta_x,points(:,2),points(:,3),'.r','MarkerSize',3);
0138
0139 else
0140 subplot(1,3,2)
0141
0142 imshow(reshape(xa,Nr,Mc));
0143 colormap(gray)
0144 view(0,90)
0145 title('1.iteration')
0146 end
0147 title(strcat('iteration=',num2str(iter),' fitted dem'))
0148
0149 axis equal
0150 end
0151 else
0152 figure
0153 subplot(2,2,1)
0154 hold on
0155 X=([0:Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0156 Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0157 colormap(gray)
0158 mesh(X,Y,reshape(xa,Nr,Mc),'EdgeColor',[0,0,0])
0159
0160 view([-29,65])
0161 zlim([0,60])
0162 grid off
0163
0164
0165
0166
0167 title(strcat('iteration=',num2str(iter),' fitted dem'))
0168
0169 if type_data ~= 5
0170 axis equal
0171 end
0172 Icc = Np;
0173 Irr = Np + (Nr-2)*Mc;
0174 Irc = Np + (Nr-2)*Mc + Nr*(Mc-2);
0175 subplot(2,2,2)
0176 imagesc(reshape(weights(Icc+1:Irr),Nr-2,Mc));
0177 subplot(2,2,3)
0178 if type_robust ~= 2
0179 imagesc(reshape(weights(Irr+1:Irc),Nr,Mc-2));
0180 else
0181 rr=floor(sqrt(Np));
0182 imagesc(reshape(1-weights(1:rr^2),rr,rr))
0183 end
0184 subplot(2,2,4)
0185 if type_robust ~= 2
0186 imagesc(reshape(weights(Irc+1:end),Nr-1,Mc-1));
0187 else
0188 rr=floor(sqrt(Np));
0189 imagesc(reshape(out_in(1:rr^2),rr,rr))
0190 end
0191 end
0192 end
0193
0194
0195 end
0196 if type_robust > 0
0197
0198 weights_f = 0.9999*ceil(weights-0.5)+0.0001;
0199 if print_type > 1
0200 weights_final = weights_f(1:20)'
0201 end
0202 xa=zeros(length(xa),1);
0203
0204 [A,weights_f,dx,v,W,Nmatrix]=estimate_dem_robust...
0205 (N,U,Np,Nr,Mc,poi,xa,weights_f,sigma_k,1,0,0,...
0206 out_in,print_type,plot_type);
0207
0208 xa = xa + dx;
0209 if print_type > 0
0210 iter_dx = [iter+1,norm(dx)]
0211 end
0212
0213 if plot_type > 0
0214 figure
0215 subplot(2,2,1)
0216 hold on
0217 X=([0:Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0218 Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0219 colormap(gray)
0220 mesh(X,Y,reshape(xa,Nr,Mc),'EdgeColor',[0,0,0])
0221
0222 view([-29,65])
0223 zlim([0,60])
0224 grid off
0225
0226
0227
0228
0229 title(strcat('last iteration, fitted dem'))
0230
0231 if type_data ~= 5
0232 axis equal
0233 end
0234 Icc = Np;
0235 Irr = Np + (Nr-2)*Mc;
0236 Irc = Np + (Nr-2)*Mc + Nr*(Mc-2);
0237 subplot(2,2,2)
0238 imagesc(reshape(weights_f(Icc+1:Irr),Nr-2,Mc));
0239 subplot(2,2,3)
0240 if type_robust ~= 2
0241 imagesc(reshape(weights(Irr+1:Irc),Nr,Mc-2));
0242 else
0243 rr=floor(sqrt(Np));
0244 imagesc(reshape(1-weights(1:rr^2),rr,rr))
0245 end
0246 subplot(2,2,4)
0247 if type_robust ~= 2
0248 imagesc(reshape(weights(Irc+1:end),Nr-1,Mc-1));
0249 else
0250 rr=floor(sqrt(Np));
0251 imagesc(reshape(out_in(1:rr^2),rr,rr))
0252 end
0253 end
0254 end
0255
0256 dem_smoothed = reshape(xa,Nr,Mc);
0257
0258
0259 Sigma=0;
0260 S = 0;
0261
0262 if out_C > 0
0263
0264 Sigma = sparseinv(Nmatrix);
0265
0266 S = full(reshape(diag(Sigma),Nr,Mc));
0267 end
0268
0269
0270
0271 return