Home > 16-Surfacereconstruction > Functions > test_smooth_dem_robust_bilinear.m

test_smooth_dem_robust_bilinear

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 = 0;
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 end
0103 
0104 Np = size(points,1);
0105 
0106 
0107 %% interpolate
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 % figure
0117 %        sc = max(ds(:))-min(ds(:));
0118 %        imshow((ds-min(ds(:)))/sc*0.95);
0119 %
0120 % if out_C > 0
0121 %     figure
0122 %     f=max(S(:));
0123 %     mesh(sqrt(S/f*min(Nr,Mc)/3)*10)
0124 %     hold on
0125 %     title('standard deviations')
0126 %     if type_data ~= 5
0127 %         axis equal
0128 %     end
0129 % end
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     %col=(1+col./sqrt(1+col.^2))/2;
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     %mesh(X,Y,ds);
0147     
0148     plot3(points(:,1),points(:,2),points(:,3),'.r','MarkerSize',15)
0149     view=[-20,30]
0150 else
0151     figure %subplot(1,3,3)
0152     %imagesc([dem,ds])
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); %,'EdgeColor',[0,0,0])% ,129*ones(size(S))
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 % figure
0214 % hold on
0215 % mesh_edge_color(ds)
0216 % view([-68,57])
0217 % axis equal
0218 
0219 %plot_edges_dem(Np, Nr, Mc, points, out_in, w, w_f, Tol);
0220 % for p=1:Np
0221 %     if out_in(p)
0222 %         [p,points(p,1),points(p,2),w(p),w_f(p),abs(dz(p))/(3*sigma)];
0223 %     end
0224 % end
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;                    % first index for column curvatures
0260     Irr = Np + (Nr-2)*Mc;        % first index for row curvature
0261     Irc = Np + (Nr-2)*Mc + Nr*(Mc-2);    % first index for torsion
0262     % plot residuals
0263     figure
0264     % heights
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     % curvatures cc
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     % curvatures rr
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     % torsion rc
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

Generated on Sat 01-Oct-2016 21:05:04 by m2html © 2005