0001
0002
0003 clear all
0004 close all
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025 init_rand = 6;
0026 type_data = 16;
0027 type_robust = 0
0028
0029
0030
0031
0032
0033 out_C = 0;
0034 out_print = 0;
0035 print_type = 0;
0036 plot_type = 0;
0037 Tol =0.15;
0038
0039 if init_rand==0
0040 init_rand=round(rand( sum(100*clock)))
0041 r=rand(3,1);
0042 end;
0043 rand('state',init_rand);
0044 randn('state',init_rand);
0045
0046
0047
0048 switch type_data
0049 case 0
0050 [points,BB,delta_x,sigma_k,sigma_h,dem]=simulate_points_dem_0
0051 out_in =ones(size(points),1);
0052 case 1
0053 [points,BB,delta_x,sigma_k,dem]=simulate_points_dem_1;
0054 case 2
0055 [points,BB,delta_x,sigma_k,dem]=simulate_points_dem_2;
0056 case 3
0057 [points,BB,delta_x,sigma_k,dem]=simulate_points_dem_3;
0058 case 4
0059 [points,BB,delta_x,sigma_k,dem]=simulate_points_dem_4;
0060 case 5
0061 im_name = 'IMG_2092_4.JPG';
0062 im_name = 'IMG_2726-2.JPG';
0063 im_name = 'IMG_8004-4.JPG';
0064
0065
0066
0067 im_name = 'IMG_8942-4.JPG';
0068 [points,BB,delta_x,sigma_k,sigma_h,dem]=simulate_points_dem_5(im_name);
0069 case 6
0070 [points,BB,delta_x,sigma_k,sigma_s,sigma_h,out_in,dem]=simulate_points_dem_6(70,10,10,1,0.05);
0071 [points,BB,delta_x,sigma_k,sigma_s,sigma_h,out_in,dem]=simulate_points_dem_6(70,10,10,0.5,0.05);
0072 case 7
0073 [points,BB,delta_x,sigma_k,dem]=simulate_points_dem_7(30,15,15,1);
0074 dem0=dem;
0075 case 8
0076 [points,BB,delta_x,sigma_k,dem]=simulate_points_dem_8(10,15,15,1);
0077 dem0=dem;
0078 case 9
0079 sigma=0.001;
0080 [points,BB,delta_x,sigma_k,dem,out_in,dz]=simulate_points_dem_9(2000,0.15,sigma);
0081 dem0=dem;
0082 case 10
0083 [points,BB,delta_x,sigma_k,XYZo]=...
0084 simulate_points_dem_10('fa2_aurelo_result_pyra0_ausgeschnitten-3.ply',40);
0085
0086 case 11
0087 [points,BB,delta_x,sigma_k,sigma_s,dem]=simulate_points_dem_0_flat
0088 case 12
0089 [points,BB,delta_x,sigma_k,sigma_s,dem]=simulate_points_dem_cone
0090 case 13
0091 [points,BB,delta_x,sigma_k,sigma_s,dem]=simulate_points_dem_steps(60,4,5)
0092 case 14
0093 sigma=0.001;
0094 [points,BB,delta_x,sigma_k,dem,out_in,dem]=simulate_points_dem_14(2000,0.15,sigma);
0095 dem0=dem;
0096
0097 case 15
0098 d=8;
0099 n=4;
0100 [points,BB,delta_x,sigma_k,sigma_s,out_in,dem]=simulate_points_dem_15_flat(d,n)
0101 out_C = 1;
0102 case 16
0103 [points,BB,delta_x,sigma_k,sigma_s,out_in,dem]=simulate_points_dem_square
0104 end
0105
0106 Np = size(points,1);
0107
0108
0109
0110
0111 starttime = cputime
0112 [ds,S,Sigma,Np,Nr,Mc,ver,A,w,w_f,W] =...
0113 smooth_dem_robust_bilinear...
0114 (points,BB,delta_x,sigma_k,out_C,type_robust,type_data,...
0115 out_in,print_type,plot_type)
0116 complete_time_for_solution=cputime-starttime
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131
0132
0133 figure
0134 hold on
0135 X=([0:Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0136 Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0137
0138 if type_data ~= 5
0139 colormap(gray)
0140 shade = [-1 0; 0 1]/sqrt(2);
0141 shade = -[ 0 -1; 1 0]/sqrt(2);
0142 col=conv2(ds,shade,'same');
0143
0144 colo=min(1, max(0,tanh(100*(col-0)/(max(col(:))-min(col(:))))));
0145 surf(X,Y,ds,colo,'FaceLighting','gouraud','EdgeColor','none')
0146 alpha(0.3);
0147 shading interp
0148
0149
0150 plot3(points(:,1),points(:,2),points(:,3),'.r','MarkerSize',15)
0151 view=[-20,30]
0152 else
0153 figure
0154
0155 imshow([dem,ds])
0156 colormap(gray)
0157 end
0158
0159 title('fitted dem with given points')
0160
0161 if type_data ~= 5
0162 axis equal
0163 end
0164
0165 figure
0166 if type_data == 0
0167 subplot(1,2,1)
0168
0169 hold on
0170 surf(X,Y,ds);
0171 colormap(gray)
0172
0173 title('fitted dem - min gradient')
0174 alpha(0.3)
0175 col=zeros(Np,2);
0176 for n=1:Np
0177 plot3(points(n,1),points(n,2),points(n,3),'.k','MarkerSize',15)
0178 end
0179 view([-29,65])
0180 xlim([BB(1),BB(3)])
0181 ylim([BB(2),BB(4)])
0182 zlim([0,10])
0183 axis equal
0184
0185 subplot(1,2,2)
0186 hold on
0187 surf(X,Y,sqrt(S)*sigma_h,'EdgeColor',[0,0,0])
0188 colormap(cool);
0189 title('fitted dem - sigmas')
0190 col=zeros(Np,2);
0191 alpha(0.3)
0192 for n=1:Np
0193 z=interpolate_bilinear(S,points(n,1),points(n,2),delta_x,BB,sigma_h);
0194 plot3(points(n,1),points(n,2),z,'.k','MarkerSize',15)
0195 end
0196 plot3(0,0,0,'.k')
0197 view([-29,65])
0198 xlim([BB(1),BB(3)])
0199 ylim([BB(2),BB(4)])
0200 zlim([-1,7])
0201 axis equal
0202
0203 end
0204
0205
0206
0207 figure
0208
0209 mesh(ds,'EdgeColor',[0,0,0])
0210 axis equal
0211 view([-33,63])
0212 zlim([-30,70])
0213 grid off
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227
0228
0229 figure
0230 hold on
0231 mesh(S);
0232 title('Heigh \sigma_z')
0233 view([-33,63])
0234 sigma_centre=S(1+(n-1)*d,1+(n-1)*d)
0235 return
0236
0237 if type_data ~= 5
0238 axis equal
0239 end
0240 figure
0241 imagesc(ds-dem)
0242 colormap(gray)
0243
0244 figure
0245 mesh(ds-dem)
0246 colormap(gray)
0247
0248 full(inv(W)*A);
0249
0250 figure
0251 mesh(ds)
0252 colormap(gray);
0253 axis equal
0254
0255 number_of_outliers = sum(out_in)
0256
0257
0258 return
0259
0260 if type_data == 1
0261 Icc = Np;
0262 Irr = Np + (Nr-2)*Mc;
0263 Irc = Np + (Nr-2)*Mc + Nr*(Mc-2);
0264
0265 figure
0266
0267 subplot(1,4,1)
0268 res=zeros(Nr,Mc);
0269 vg = reshape(ver(1:Np),Nr-2,Mc-2);
0270 for i=1:Nr-2
0271 for j=1:Mc-2
0272 res(i+1,j+1)=vg(i,j);
0273 end
0274 end
0275 mesh(-res(:)');
0276 title('residuals - heights')
0277
0278
0279 subplot(1,4,2)
0280 res=zeros(Nr,Mc);
0281 vg = reshape(ver(Icc+1:Irr),Nr-2,Mc);
0282 for i=1:Nr-2
0283 for j=1:Mc
0284 res(i+1,j)=vg(i,j);
0285 end
0286 end
0287 mesh(-res'*sqrt(sigma_k));
0288 title('residuals - curvatures cc')
0289
0290
0291
0292 subplot(1,4,3)
0293 res=zeros(Nr,Mc);
0294 vg = reshape(ver(Irr+1:Irc),Nr,Mc-2);
0295 for i=1:Nr
0296 for j=1:Mc-2
0297 res(i,j+1)=vg(i,j);
0298 end
0299 end
0300 mesh(-res'*sqrt(sigma_k));
0301 title('residuals - curvatures rr')
0302
0303
0304
0305 subplot(1,4,4)
0306 res=zeros(Nr,Mc);
0307 vg = reshape(ver(Irc+1:end),Nr-1,Mc-1);
0308 for i=1:Nr-1
0309 for j=1:Mc-1
0310 res(i,j)=vg(i,j);
0311 end
0312 end
0313 mesh(-res'*sqrt(sigma_k));
0314 title('residuals - torsion rc')
0315 end