test_smooth_dem_robust_bilinear_square

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