smooth_dem_robust_bilinear_flat

% smooth_dem_robust_bilinear, sparse data on grid, flatness

 
 function [dem_smoothed,S,Sigma,Np,Nr,Mc,v,A,weights,W] = ...
    smooth_dem_robust_bilinear_flat...
    (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_s = std. dev. of slope 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)  for
               points
               first derivatives c
               first derivatives r
 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 10/2014

0001 %% smooth_dem_robust_bilinear, sparse data on grid, flatness
0002 %
0003 %
0004 % function [dem_smoothed,S,Sigma,Np,Nr,Mc,v,A,weights,W] = ...
0005 %    smooth_dem_robust_bilinear_flat...
0006 %    (points,BB,delta_x,sigma_k,out_C,out_print,type_robust,type_data)
0007 %
0008 % points  = Nx4-matrix of (x,y,h,sigma) values,
0009 % BB      = [xmin,ymin,xmax,ymax] boundin,g box of square grid
0010 % delta_x = grid size, squared
0011 %           xmax, ymax will be adapted, not decreased
0012 %           points not incide BB will be discarded
0013 % sigms_s = std. dev. of slope of surface, referring to delta_x
0014 % out_C   = 0 no covariance matrix as output, (default if U>1600)
0015 %           1 covariance matrix as output
0016 % out_print = 0 no output
0017 %             1 little output
0018 %             2 much output
0019 % type_robust = 0,1,2,3 (00,01,10,11) for points and dem
0020 %
0021 % dem_smoothe = Nr x Mc matrix of smoothed heigths
0022 % S           = corresponding standard deviations, if required
0023 % Sigma       = full covariance matrix. if required
0024 % Np          = number of points
0025 % Nr, Mc      = number of grid points in x- and y-direction
0026 % v           = residuals (Np + (Nr-2)*Mc + Nr*(Mc-2)  for
0027 %               points
0028 %               first derivatives c
0029 %               first derivatives r
0030 % A           = design matrix (sparse)
0031 % weights     = weight-factors after robust estimation (0 or 1)
0032 % W           = weight of observations including both, uncertainty and robust factor
0033 %
0034 % grid is matrix Nr x Mc and stored column wise
0035 %
0036 % wf 10/2014
0037 %
0038 
0039 function [dem_smoothed,S,Sigma,Np,Nr,Mc,v,A,weights,W] = ...
0040     smooth_dem_robust_bilinear_flat...
0041     (points,BB,delta_x,sigma_s,out_C,type_robust,...
0042     print_type,plot_type)
0043 
0044 % number of points, preliminary
0045 [Npp,dim] = size(points);
0046 
0047 if print_type > 0
0048     number_of_points_given = Npp
0049 end
0050 
0051 
0052 xmin = BB(1);
0053 ymin = BB(2);
0054 xmax = BB(3);
0055 ymax = BB(4);
0056 Nr =ceil((xmax-xmin)/delta_x)+1;
0057 Mc =ceil((ymax-ymin)/delta_x)+1;
0058 
0059 
0060 Np = 0;
0061 for p=1:Npp
0062     if points(p,1) >= xmin && ...
0063             points(p,1) < xmax && ...
0064             points(p,2) >= ymin && ...
0065             points(p,2) < ymax
0066         Np=Np+1;
0067         poi(Np,:) = [(points(p,1)-xmin)/delta_x,(points(p,2)-ymin)/delta_x,points(p,3:4)];
0068     end
0069 end
0070 
0071 if print_type > 1
0072     number_of_points_in_BB = Np
0073 end
0074 
0075 
0076 % number of unknowns and observations incl. prior
0077 
0078 U = Nr*Mc;
0079 N = Np + (Nr-1)*Mc + Nr*(Mc-1);
0080 if print_type > 0
0081     number_of_unknowns = U
0082 end
0083 
0084 % do not give covariance matrix for U>1600.
0085 if U > 1600
0086     out_C = 0;
0087 end
0088 if print_type > 0
0089     tic
0090 end
0091 %% estimate iteratively
0092 
0093 iter_switch=3;
0094 if type_robust == 0
0095     N_iter = 1;
0096     first = 1;
0097 else
0098     N_iter=6;
0099     first=1;
0100 end
0101 xa = zeros(U,1);
0102 weights  = ones(N,1);
0103 sigma_s0=sigma_s;
0104 for iter = 1:N_iter
0105     if print_type > 0
0106         disp('Call estimate_dem_robust_flat')
0107         iter_typerobust=[iter,type_robust]
0108     end
0109     L1=0;
0110     if type_robust >0 && iter <= iter_switch
0111         L1=1;
0112         sigma_s=sigma_s0*4;
0113     else
0114         L1=0;
0115         sigma_s=sigma_s0;
0116     end
0117     sigma_s=sigma_s0*1.4^(N_iter-iter);
0118     [A,weights,dx,v,W]=estimate_dem_robust_flat...
0119         (N,U,Np,Nr,Mc,poi,xa,weights,sigma_s,iter,type_robust,L1,...
0120         print_type,plot_type);
0121     first=0;
0122     xa = xa + dx;
0123     
0124 
0125     if print_type > 0
0126         iter_dx = [iter,norm(dx)]
0127     end
0128     %%
0129     if  plot_type > 0
0130         if iter == 1
0131             figure(20)
0132             hold on
0133             X=([4*Nr:4*Nr+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 && type_robust > 0
0137                 colormap(gray)
0138                 surf(X,Y,reshape(xa,Nr,Mc));
0139                 alpha(0.3);
0140                 %mesh(ds);
0141 
0142                 plot3(points(:,1)+4*Nr*delta_x,points(:,2),points(:,3),'.r','MarkerSize',20);
0143 
0144             else
0145                 subplot(1,3,2)
0146                 %imagesc(reshape(xa,Nr,Mc));
0147                 imshow(reshape(xa,Nr,Mc));
0148                 colormap(gray)
0149                 view(0,90)
0150                 title('1.iteration')
0151             end
0152             title(strcat('iteration=',num2str(iter),' fitted dem'))
0153 
0154             axis equal
0155         else
0156             figure
0157             hold on
0158             X=([0:Nr-1]'*ones(1,Mc))*delta_x+BB(1);
0159             Y=(ones(Nr,1)*[0:Mc-1])*delta_x+BB(2);
0160             colormap(gray)
0161             surf(X,Y,reshape(xa,Nr,Mc));
0162             alpha(0.3);
0163             %mesh(ds);
0164             plot3(points(:,1),points(:,2),points(:,3),'.r','MarkerSize',20);
0165 
0166             title(strcat('iteration=',num2str(iter),' fitted dem'))
0167 
0168             if type_data ~= 5
0169                 axis equal
0170             end
0171         end
0172     end
0173 
0174 end
0175 % if type_robust > 0
0176 %     %final estimate
0177 %     weights = 0.9999*ceil(weights-0.5)+0.0001;
0178 %     weights_final = weights(1:20)'
0179 %     xa=zeros(length(xa),1);
0180 %     [A,weights,dx,v]=estimate_dem_robust_flat(N,U,Np,Nr,Mc,poi,xa,weights,sigma_s,1,out_print,0);
0181 %     xa = xa + dx;
0182 %     iter_dx = [iter+1,norm(dx)]
0183 % end
0184 
0185 dem_smoothed = reshape(xa,Nr,Mc);
0186 
0187 %% determine covariance matrix
0188 Sigma=0;
0189 S = 0;
0190 if print_type > 0
0191     tic
0192 end
0193 if out_C > 0
0194     Sigma = full(inv(A'*A));
0195     S     = full(reshape(diag(Sigma),Nr,Mc));
0196 end
0197 if print_type > 0
0198     time_for_inverse=toc
0199 end
0200 
0201 
0202 return