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

smooth_dem_robust_bilinear

PURPOSE ^

% smooth_dem_robust_bilinear

SYNOPSIS ^

function [dem_smoothed,S,Sigma,Np,Nr,Mc,v,A,weights,weights_f,W] =smooth_dem_robust_bilinear(points,BB,delta_x,sigma_k,out_C,type_robust,type_data,out_in,print_type,plot_type)

DESCRIPTION ^

% smooth_dem_robust_bilinear

 reconstruct surface with 
   energy = sum_m rho(l_m-a_m)/sigma_n) + 
             1/sigma_d sum_ij d_ii^2+2d_ij^2+djj^2
 l_m = observed heights, interpolated bilinearly
 a_m unknown grid heights, d_ij second derivatives
 
 function [dem_smoothed,S,Sigma,Np,Nr,Mc,v,A,weights,W] = ...
    smooth_dem_robust_bilinear...
    (points,BB,delta_x,sigma_k,out_C,out_print,type_robust,type_data)

 points  = Nx4-matrix of (x,y,h,sigma) values, 
 BB      = [xmin,ymin,xmax,ymax] boundin,g box of square grid
 delta_x = grid size, squared 
           xmax, ymax will be adapted, not decreased
           points not incide BB will be discarded
 sigms_k = std. dev. of curvature of surface, referring to delta_x
 out_C   = 0 no covariance matrix as output, (default if U>1600)
           1 covariance matrix as output
 out_print = 0 no output
             1 little output
             2 much output
 type_robust = 0,1,2,3 (00,01,10,11) for points and dem

 dem_smoothe = Nr x Mc matrix of smoothed heigths
 S           = corresponding standard deviations, if required
 Sigma       = full covariance matrix. if required
 Np          = number of points
 Nr, Mc      = number of grid points in x- and y-direction
 v           = residuals (Np + (Nr-2)*Mc + Nr*(Mc-2) + (Nr-1)*(Mc-1)) for
               points
               second derivatives cc
               second derivatives rr
               torsion
 A           = design matrix (sparse) 
 weights     = weight-factors after robust estimation (0 or 1)
 W           = weight of observations including both, uncertainty and robust factor 

 grid is matrix Nr x Mc and stored column wise

 wf 7/2014

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 %% smooth_dem_robust_bilinear
0002 %
0003 % reconstruct surface with
0004 %   energy = sum_m rho(l_m-a_m)/sigma_n) +
0005 %             1/sigma_d sum_ij d_ii^2+2d_ij^2+djj^2
0006 % l_m = observed heights, interpolated bilinearly
0007 % a_m unknown grid heights, d_ij second derivatives
0008 %
0009 % function [dem_smoothed,S,Sigma,Np,Nr,Mc,v,A,weights,W] = ...
0010 %    smooth_dem_robust_bilinear...
0011 %    (points,BB,delta_x,sigma_k,out_C,out_print,type_robust,type_data)
0012 %
0013 % points  = Nx4-matrix of (x,y,h,sigma) values,
0014 % BB      = [xmin,ymin,xmax,ymax] boundin,g box of square grid
0015 % delta_x = grid size, squared
0016 %           xmax, ymax will be adapted, not decreased
0017 %           points not incide BB will be discarded
0018 % sigms_k = std. dev. of curvature of surface, referring to delta_x
0019 % out_C   = 0 no covariance matrix as output, (default if U>1600)
0020 %           1 covariance matrix as output
0021 % out_print = 0 no output
0022 %             1 little output
0023 %             2 much output
0024 % type_robust = 0,1,2,3 (00,01,10,11) for points and dem
0025 %
0026 % dem_smoothe = Nr x Mc matrix of smoothed heigths
0027 % S           = corresponding standard deviations, if required
0028 % Sigma       = full covariance matrix. if required
0029 % Np          = number of points
0030 % Nr, Mc      = number of grid points in x- and y-direction
0031 % v           = residuals (Np + (Nr-2)*Mc + Nr*(Mc-2) + (Nr-1)*(Mc-1)) for
0032 %               points
0033 %               second derivatives cc
0034 %               second derivatives rr
0035 %               torsion
0036 % A           = design matrix (sparse)
0037 % weights     = weight-factors after robust estimation (0 or 1)
0038 % W           = weight of observations including both, uncertainty and robust factor
0039 %
0040 % grid is matrix Nr x Mc and stored column wise
0041 %
0042 % wf 7/2014
0043 %
0044  
0045 function [dem_smoothed,S,Sigma,Np,Nr,Mc,v,A,weights,weights_f,W] = ...
0046     smooth_dem_robust_bilinear...
0047     (points,BB,delta_x,sigma_k,out_C,type_robust,type_data,...
0048     out_in,print_type,plot_type)
0049     
0050 
0051 %% number of points, preliminary --> final (which are in BB)
0052 [Npp,tmp] = size(points);
0053 
0054 xmin = BB(1);
0055 ymin = BB(2);
0056 xmax = BB(3);
0057 ymax = BB(4);
0058 Nr =ceil((xmax-xmin)/delta_x)+1;
0059 Mc =ceil((ymax-ymin)/delta_x)+1;
0060 
0061 
0062 Np = 0;
0063 for p=1:Npp
0064     if points(p,1) >= xmin && ...
0065             points(p,1) < xmax && ...
0066             points(p,2) >= ymin && ...
0067             points(p,2) < ymax
0068         Np=Np+1;
0069         poi(Np,:) = [(points(p,1)-xmin)/delta_x,(points(p,2)-ymin)/delta_x,points(p,3:4)];
0070         
0071     end
0072 end
0073             
0074 %% number of unknowns and observations incl. prior
0075 U = Nr*Mc;
0076 N = Np + (Nr-2)*Mc + Nr*(Mc-2) + (Nr-1)*(Mc-1);
0077             if print_type > 0
0078                 number_of_unknowns = U %#ok<NOPRT,NASGU>
0079             end
0080 
0081 %% do not give covariance matrix for U>10000.
0082             if U > 40401
0083                 out_C = 0;
0084             end
0085             if print_type > 0 
0086                 tic
0087             end
0088 %% estimate iteratively
0089 
0090 iter_switch=3;
0091 if type_robust == 0
0092     N_iter = 1;
0093 else
0094     N_iter=6;
0095 end
0096 
0097 xa = zeros(U,1);               % initial unknown parameters
0098 weights   = ones(N,1);         % weight-factors for robust estimation
0099 weights_f = weights;
0100 sigma_k0  = sigma_k;
0101 gain=1.0;
0102 
0103 for iter = 1:N_iter
0104     L1 = 0;
0105     if type_robust > 0 && iter <= iter_switch
0106         L1=1;
0107     else
0108         L1=0;
0109     end
0110     sigma_k=sigma_k0*gain^(N_iter-iter);     % adapt prior sigma
0111     
0112     [A,weights,dx,v,W,Nmatrix]=estimate_dem_robust...
0113         (N,U,Np,Nr,Mc,poi,xa,weights,sigma_k,iter,type_robust,L1,...
0114         out_in,print_type,plot_type);
0115     
0116     xa = xa + dx;
0117      
0118     if print_type > 0
0119         iter_dx = [iter,norm(dx)] %#ok<NOPRT,NASGU>
0120     end
0121     %%
0122                 if plot_type > 0
0123                     if iter == 1
0124                         if type_robust > 0
0125                         figure
0126                         hold on
0127                         title('surf')
0128                         X=([4*Nr:4*Nr+Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0129                         Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0130 
0131                             if type_data ~=5 
0132                                 colormap(gray)
0133                                 surf(X,Y,reshape(xa,Nr,Mc));
0134                                 alpha(0.3);
0135                                 %mesh(ds);
0136                                 view([-29,65])
0137                                 plot3(points(:,1)+4*Nr*delta_x,points(:,2),points(:,3),'.r','MarkerSize',3);
0138 
0139                             else
0140                                 subplot(1,3,2)
0141                                 %imagesc(reshape(xa,Nr,Mc));
0142                                 imshow(reshape(xa,Nr,Mc));
0143                                 colormap(gray)
0144                                 view(0,90)
0145                                 title('1.iteration')
0146                             end
0147                             title(strcat('iteration=',num2str(iter),' fitted dem'))
0148 
0149                             axis equal
0150                         end
0151                     else
0152                         figure
0153                         subplot(2,2,1)
0154                         hold on
0155                         X=([0:Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0156                         Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0157                         colormap(gray)
0158                         mesh(X,Y,reshape(xa,Nr,Mc),'EdgeColor',[0,0,0])%  %mesh(ds,'EdgeColor',[0,0,0])
0159                         %axis equal
0160                         view([-29,65])
0161                         zlim([0,60])
0162                         grid off
0163                         %alpha(0.3);
0164                         %mesh(ds);
0165                         %plot3(points(:,1),points(:,2),points(:,3),'.r','MarkerSize',5);
0166 
0167                         title(strcat('iteration=',num2str(iter),' fitted dem'))
0168 
0169                         if type_data ~= 5
0170                             axis equal
0171                         end
0172                         Icc = Np;                    % first index for column curvatures
0173                         Irr = Np + (Nr-2)*Mc;        % first index for row curvature
0174                         Irc = Np + (Nr-2)*Mc + Nr*(Mc-2);    % first index for torsion
0175                         subplot(2,2,2)
0176                         imagesc(reshape(weights(Icc+1:Irr),Nr-2,Mc));
0177                         subplot(2,2,3)
0178                         if type_robust ~= 2
0179                             imagesc(reshape(weights(Irr+1:Irc),Nr,Mc-2));
0180                         else
0181                             rr=floor(sqrt(Np));
0182                             imagesc(reshape(1-weights(1:rr^2),rr,rr))
0183                         end
0184                         subplot(2,2,4)
0185                         if type_robust ~= 2
0186                             imagesc(reshape(weights(Irc+1:end),Nr-1,Mc-1));
0187                         else
0188                             rr=floor(sqrt(Np));
0189                             imagesc(reshape(out_in(1:rr^2),rr,rr))
0190                         end
0191                     end
0192                 end
0193 
0194 
0195 end
0196 if type_robust > 0
0197     %final estimate
0198     weights_f = 0.9999*ceil(weights-0.5)+0.0001;
0199     if print_type > 1
0200         weights_final = weights_f(1:20)' %#ok<NOPRT,NASGU>
0201     end
0202     xa=zeros(length(xa),1);
0203     
0204     [A,weights_f,dx,v,W,Nmatrix]=estimate_dem_robust...
0205         (N,U,Np,Nr,Mc,poi,xa,weights_f,sigma_k,1,0,0,...
0206         out_in,print_type,plot_type);
0207     
0208     xa = xa + dx;
0209                 if print_type > 0
0210                     iter_dx = [iter+1,norm(dx)] %#ok<NOPRT,NASGU>
0211                 end
0212 
0213                 if plot_type > 0
0214                 figure
0215                     subplot(2,2,1)
0216                     hold on
0217                     X=([0:Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0218                     Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0219                     colormap(gray)
0220                     mesh(X,Y,reshape(xa,Nr,Mc),'EdgeColor',[0,0,0])%  %mesh(ds,'EdgeColor',[0,0,0])
0221                     %axis equal
0222                     view([-29,65])
0223                     zlim([0,60])
0224                     grid off
0225                     %alpha(0.3);
0226                     %mesh(ds);
0227                     %plot3(points(:,1),points(:,2),points(:,3),'.r','MarkerSize',5);
0228 
0229                     title(strcat('last iteration, fitted dem'))
0230 
0231                     if type_data ~= 5
0232                         axis equal
0233                     end
0234                     Icc = Np;                    % first index for column curvatures
0235                     Irr = Np + (Nr-2)*Mc;        % first index for row curvature
0236                     Irc = Np + (Nr-2)*Mc + Nr*(Mc-2);    % first index for torsion
0237                     subplot(2,2,2)
0238                     imagesc(reshape(weights_f(Icc+1:Irr),Nr-2,Mc));
0239                     subplot(2,2,3)
0240                     if type_robust ~= 2
0241                         imagesc(reshape(weights(Irr+1:Irc),Nr,Mc-2));
0242                     else
0243                        rr=floor(sqrt(Np));
0244                        imagesc(reshape(1-weights(1:rr^2),rr,rr))
0245                     end
0246                     subplot(2,2,4)       
0247                     if type_robust ~= 2
0248                        imagesc(reshape(weights(Irc+1:end),Nr-1,Mc-1));
0249                     else
0250                        rr=floor(sqrt(Np));
0251                        imagesc(reshape(out_in(1:rr^2),rr,rr))
0252                     end
0253                 end
0254 end
0255 
0256 dem_smoothed = reshape(xa,Nr,Mc);
0257 
0258 %% determine covariance matrix for quality analysis
0259 Sigma=0;
0260 S = 0;
0261 
0262 if out_C > 0    
0263                   
0264     Sigma = sparseinv(Nmatrix);       % use sparse inverse
0265                 
0266     S     = full(reshape(diag(Sigma),Nr,Mc));
0267 end
0268 
0269 
0270 
0271 return

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