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 = 0;
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 end
0103
0104 Np = size(points,1);
0105
0106
0107
0108
0109 starttime = cputime
0110 [ds,S,Sigma,Np,Nr,Mc,ver,A,w,w_f,W] =...
0111 smooth_dem_robust_bilinear...
0112 (points,BB,delta_x,sigma_k,out_C,type_robust,type_data,...
0113 out_in,print_type,plot_type)
0114 complete_time_for_solution=cputime-starttime
0115
0116
0117
0118
0119
0120
0121
0122
0123
0124
0125
0126
0127
0128
0129
0130
0131 figure
0132 hold on
0133 X=([0:Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0134 Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0135
0136 if type_data ~= 5
0137 colormap(gray)
0138 shade = [-1 0; 0 1]/sqrt(2);
0139 shade = -[ 0 -1; 1 0]/sqrt(2);
0140 col=conv2(ds,shade,'same');
0141
0142 colo=min(1, max(0,tanh(100*(col-0)/(max(col(:))-min(col(:))))));
0143 surf(X,Y,ds,colo,'FaceLighting','gouraud','EdgeColor','none')
0144 alpha(0.3);
0145 shading interp
0146
0147
0148 plot3(points(:,1),points(:,2),points(:,3),'.r','MarkerSize',15)
0149 view=[-20,30]
0150 else
0151 figure
0152
0153 imshow([dem,ds])
0154 colormap(gray)
0155 end
0156
0157 title('fitted dem with given points')
0158
0159 if type_data ~= 5
0160 axis equal
0161 end
0162
0163 figure
0164 if type_data == 0
0165 subplot(1,2,1)
0166
0167 hold on
0168 surf(X,Y,ds);
0169 colormap(gray)
0170
0171 title('fitted dem - min gradient')
0172 alpha(0.3)
0173 col=zeros(Np,2);
0174 for n=1:Np
0175 plot3(points(n,1),points(n,2),points(n,3),'.k','MarkerSize',15)
0176 end
0177 view([-29,65])
0178 xlim([BB(1),BB(3)])
0179 ylim([BB(2),BB(4)])
0180 zlim([0,10])
0181 axis equal
0182
0183 subplot(1,2,2)
0184 hold on
0185 surf(X,Y,sqrt(S)*sigma_h,'EdgeColor',[0,0,0])
0186 colormap(cool);
0187 title('fitted dem - sigmas')
0188 col=zeros(Np,2);
0189 alpha(0.3)
0190 for n=1:Np
0191 z=interpolate_bilinear(S,points(n,1),points(n,2),delta_x,BB,sigma_h);
0192 plot3(points(n,1),points(n,2),z,'.k','MarkerSize',15)
0193 end
0194 plot3(0,0,0,'.k')
0195 view([-29,65])
0196 xlim([BB(1),BB(3)])
0197 ylim([BB(2),BB(4)])
0198 zlim([-1,7])
0199 axis equal
0200
0201 end
0202
0203
0204
0205 figure
0206
0207 mesh(ds,'EdgeColor',[0,0,0])
0208 axis equal
0209 VIEW([-33,63])
0210 zlim([-30,70])
0211 grid off
0212
0213
0214
0215
0216
0217
0218
0219
0220
0221
0222
0223
0224
0225
0226
0227 figure
0228 hold on
0229 mesh(S);
0230 title('Heigh \sigma_z')
0231 view([-33,63])
0232 sigma_centre=S(1+(n-1)*d,1+(n-1)*d)
0233 return
0234
0235 if type_data ~= 5
0236 axis equal
0237 end
0238 figure
0239 imagesc(ds-dem)
0240 colormap(gray)
0241
0242 figure
0243 mesh(ds-dem)
0244 colormap(gray)
0245
0246 full(inv(W)*A);
0247
0248 figure
0249 mesh(ds)
0250 colormap(gray);
0251 axis equal
0252
0253 number_of_outliers = sum(out_in)
0254
0255
0256 return
0257
0258 if type_data == 1
0259 Icc = Np;
0260 Irr = Np + (Nr-2)*Mc;
0261 Irc = Np + (Nr-2)*Mc + Nr*(Mc-2);
0262
0263 figure
0264
0265 subplot(1,4,1)
0266 res=zeros(Nr,Mc);
0267 vg = reshape(ver(1:Np),Nr-2,Mc-2);
0268 for i=1:Nr-2
0269 for j=1:Mc-2
0270 res(i+1,j+1)=vg(i,j);
0271 end
0272 end
0273 mesh(-res(:)');
0274 title('residuals - heights')
0275
0276
0277 subplot(1,4,2)
0278 res=zeros(Nr,Mc);
0279 vg = reshape(ver(Icc+1:Irr),Nr-2,Mc);
0280 for i=1:Nr-2
0281 for j=1:Mc
0282 res(i+1,j)=vg(i,j);
0283 end
0284 end
0285 mesh(-res'*sqrt(sigma_k));
0286 title('residuals - curvatures cc')
0287
0288
0289
0290 subplot(1,4,3)
0291 res=zeros(Nr,Mc);
0292 vg = reshape(ver(Irr+1:Irc),Nr,Mc-2);
0293 for i=1:Nr
0294 for j=1:Mc-2
0295 res(i,j+1)=vg(i,j);
0296 end
0297 end
0298 mesh(-res'*sqrt(sigma_k));
0299 title('residuals - curvatures rr')
0300
0301
0302
0303 subplot(1,4,4)
0304 res=zeros(Nr,Mc);
0305 vg = reshape(ver(Irc+1:end),Nr-1,Mc-1);
0306 for i=1:Nr-1
0307 for j=1:Mc-1
0308 res(i,j)=vg(i,j);
0309 end
0310 end
0311 mesh(-res'*sqrt(sigma_k));
0312 title('residuals - torsion rc')
0313 end