simulate_points_dem_16_precision

% generate a set of points for dem-precision analysis

 N           = number of points
 sigma_h     = precision of points
 Resolution  = Number of grid cells in each coordinate axis

 points      = Nx4 matrix of points
               [X,Y,Z,sigma]
 BB          = bounding box
 dx          = sampling distance
 dem         = intended dem (not used here)
 bin_im      = binary image of points
 dist_t      = distance transformation of points

 wf 8/2015

0001 %% generate a set of points for dem-precision analysis
0002 %
0003 % N           = number of points
0004 % sigma_h     = precision of points
0005 % Resolution  = Number of grid cells in each coordinate axis
0006 %
0007 % points      = Nx4 matrix of points
0008 %               [X,Y,Z,sigma]
0009 % BB          = bounding box
0010 % dx          = sampling distance
0011 % dem         = intended dem (not used here)
0012 % bin_im      = binary image of points
0013 % dist_t      = distance transformation of points
0014 %
0015 % wf 8/2015
0016 
0017 
0018 function [points,BB,dx,dem,bin_im,dist_t]=simulate_points_dem_16_precision(N,sigma_h,R)
0019 
0020 
0021 %% Example PCV (do not change)
0022 % bounding box
0023 BB=[0,0,R,R];
0024 % grid size
0025 factor=1;
0026 dx = 1/factor;
0027 
0028 p    = (1+rand(N,2)*(R-2));
0029 M=randn(4);
0030 a =[ones(N,1) (2*p(:,1)-R)/R (4*p(:,1).^2-4*p(:,1)*R-R^2)/(2*R^2) ...
0031     (2*p(:,1)-R).*(2*p(:,1).^2-2*p(:,1)*R-R^2)/R^3];
0032 b =[ones(N,1) (2*p(:,2)-R)/R (4*p(:,2).^2-4*p(:,2)*R-R^2)/(2*R^2) ...
0033     (2*p(:,2)-R).*(2*p(:,2).^2-2*p(:,2)*R-R^2)/R^3];
0034 z = 2*diag(a*M*b');
0035 points = [p,z,ones(N,1)*sigma_h];
0036  
0037 k=0;
0038 for i=[1,round(R/2),R-1]
0039     for j=[1,round(R/2),R-1]
0040         k=k+1;
0041         points(k,:)=[i,j,2*[1 (2*i-R)/(2*R) 1/2*(4*i^2-4*i*R-R^2)/R^2 ...
0042             (2*i-R)*(2*i^2-2*i*R-R^2)/R^3]*M*...
0043             [1 (2*j-R)/(2*R) 1/2*(4*j^2-4*j*R-R^2)/R^2 ...
0044             (2*j-R)*(2*j^2-2*j*R-R^2)/R^3]',sigma_h];
0045     end
0046 end
0047 
0048 
0049 
0050 %% plot points and distance transform
0051 
0052 xmin = BB(1);
0053 ymin = BB(2);
0054 xmax = BB(3);
0055 ymax = BB(4);
0056 
0057 Nr = ceil((xmax-xmin)/dx)+1;
0058 Mc = ceil((ymax-ymin)/dx)+1;
0059 
0060 dem = zeros(Nr,Mc);
0061 
0062 bin_im = dem;
0063 for n=1:N
0064     i=ceil(points(n,1)*factor);
0065     j=ceil(points(n,2)*factor);
0066     bin_im(R+1-j,i)=1;
0067 end
0068 
0069 ss = plot_init;
0070 figure('Color','w','Position',[0.1*ss(1),0.2*ss(2),0.7*ss(1),0.6*ss(2)])
0071 subplot(1,2,1)
0072 imshow(bin_im)
0073 axis equal
0074 dist_t=bwdist(bin_im,'euclidean');
0075 subplot(1,2,2)
0076 imagesc(sqrt(dist_t));
0077 title('points and distance transform')
0078 axis equal
0079 axis off
0080 
0081 
0082 
0083 return