Home > 16-Surface-Reconstruction > Surfacereconstruction > Functions > test_smooth_dem_robust_bilinear_square.m

test_smooth_dem_robust_bilinear_square

PURPOSE ^

% test dem smoothing

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

% test dem smoothing

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %% test dem smoothing
0002 
0003 clear all
0004 close all
0005 % type_data: type of generated data
0006 % 0     a set of given points
0007 % 1     a roof parallel to the grid
0008 % 2     a roof skew to the grid
0009 % 3     a roof skew to the grid + a flat roof, points can be thinned out
0010 % 4     same as 3, but with points not sitting on a grid
0011 % 5     read an image
0012 % 6     Uwe's test image
0013 % 7     cake
0014 % 8     1/4 cake
0015 % 9     constant height with outliers
0016 % 10    read point cloud
0017 % 11    flat
0018 % 12    cone
0019 % 13    cake
0020 % 14    flat with outliers
0021 % 15    flat 4x4 points
0022 %
0023 % wf 7/2014
0024 
0025 init_rand = 6;
0026 type_data = 16;
0027 type_robust = 0
0028 %             0 not robust
0029 %             1 only dem
0030 %             2 only points
0031 %             3 both
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 %% generate dem point cloud
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         % im_name = 'IMG_2425-2.JPG';  % Figur Paris
0065         %im_name = 'IMG_5977-4.JPG'; % Zebra
0066         %im_name = 'IMG_7452-4.JPG'; % Katze
0067         im_name = 'IMG_8942-4.JPG'; % Fassade Prag
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 %% interpolate
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 % figure
0119 %        sc = max(ds(:))-min(ds(:));
0120 %        imshow((ds-min(ds(:)))/sc*0.95);
0121 %
0122 % if out_C > 0
0123 %     figure
0124 %     f=max(S(:));
0125 %     mesh(sqrt(S/f*min(Nr,Mc)/3)*10)
0126 %     hold on
0127 %     title('standard deviations')
0128 %     if type_data ~= 5
0129 %         axis equal
0130 %     end
0131 % end
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     %col=(1+col./sqrt(1+col.^2))/2;
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     %mesh(X,Y,ds);
0149     
0150     plot3(points(:,1),points(:,2),points(:,3),'.r','MarkerSize',15)
0151     view=[-20,30]
0152 else
0153     figure %subplot(1,3,3)
0154     %imagesc([dem,ds])
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); %,'EdgeColor',[0,0,0])% ,129*ones(size(S))
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 % figure
0216 % hold on
0217 % mesh_edge_color(ds)
0218 % view([-68,57])
0219 % axis equal
0220 
0221 %plot_edges_dem(Np, Nr, Mc, points, out_in, w, w_f, Tol);
0222 % for p=1:Np
0223 %     if out_in(p)
0224 %         [p,points(p,1),points(p,2),w(p),w_f(p),abs(dz(p))/(3*sigma)];
0225 %     end
0226 % end
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;                    % first index for column curvatures
0262     Irr = Np + (Nr-2)*Mc;        % first index for row curvature
0263     Irc = Np + (Nr-2)*Mc + Nr*(Mc-2);    % first index for torsion
0264     % plot residuals
0265     figure
0266     % heights
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     % curvatures cc
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     % curvatures rr
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     % torsion rc
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

Generated on Mon 08-Jan-2018 17:21:49 by m2html © 2005