test_smooth_dem_robust_bilinear

% test dem smoothing

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