simulate_points_dem_10

% read ply-file

 [points,BB,delta_x,sigma_k]=simulate_points_dem_10(file_name);

 file_name     = name of ply-file
                 given in '..\Data'
                 fa2_aurelo_result_pyra0_ausgeschnitten-1.ply
                        (used in Fig. 16.23)
                 fa2_aurelo_result_pyra0_ausgeschnitten-4.ply
                        (larger data set)
 Nbin          = resolution

 points        = N x 4 vector of points and sigma
 BB            = bounding box
 delta_x       = grid distance
 sigma_k       = sigma of curvatures

 Wolfgang Förstner 07/14
 last changes: Susanne Wenzel 09/16
 wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de

0001 %% read ply-file
0002 %
0003 % [points,BB,delta_x,sigma_k]=simulate_points_dem_10(file_name);
0004 %
0005 % file_name     = name of ply-file
0006 %                 given in '..\Data'
0007 %                 fa2_aurelo_result_pyra0_ausgeschnitten-1.ply
0008 %                        (used in Fig. 16.23)
0009 %                 fa2_aurelo_result_pyra0_ausgeschnitten-4.ply
0010 %                        (larger data set)
0011 % Nbin          = resolution
0012 %
0013 % points        = N x 4 vector of points and sigma
0014 % BB            = bounding box
0015 % delta_x       = grid distance
0016 % sigma_k       = sigma of curvatures
0017 %
0018 % Wolfgang Förstner 07/14
0019 % last changes: Susanne Wenzel 09/16
0020 % wfoerstn@uni-bonn.de, wenzel@igg.uni-bonn.de
0021 
0022 function [points,BB,delta_x,sigma_k,XYZ_ori]=simulate_points_dem_10(file_name,Nbin)
0023 
0024 %% plot settings
0025 ss = plot_init;
0026 
0027 %% read file
0028 pts = ply_read(file_name);
0029 
0030 %% extract coordinates
0031 XYZ_ori = [pts.vertex.x,pts.vertex.y,pts.vertex.z];
0032 N = size(XYZ_ori,1);
0033 
0034 %% possibly rotate points to see them from the side
0035 % % determine principal axes
0036 % C = cov(XYZ_ori);
0037 % [R,D] = eigs(C);
0038 % R = R*diag(diag(sign(R)));
0039 % % normalize
0040 % XYZ = XYZ_ori*R;
0041 XYZ = XYZ_ori;
0042 
0043 %%
0044 xmin = min(XYZ(:,1));
0045 xmax = max(XYZ(:,1));
0046 ymin = min(XYZ(:,2));
0047 ymax = max(XYZ(:,2));
0048 frame_size = 1.1;
0049 bord = (frame_size-1)/2;
0050 Range_x = (xmax-xmin);
0051 Range_y = (ymax-ymin);
0052 delta_x = min(Range_x,Range_y)/(Nbin-1);
0053 Nr = ceil(Range_x*frame_size/delta_x);
0054 Mc = ceil(Range_y*frame_size/delta_x);
0055 BB = [xmin-bord*Range_x,               ymin-bord*Range_y,...
0056       xmin-bord*Range_x+(Nr-1)*delta_x,ymin-bord*Range_y+(Mc-1)*delta_x];
0057 
0058 % sigma_n
0059 factor_n = 0.1;
0060 sigma_n = factor_n*max([Range_x,Range_y])/sqrt(N);
0061 sigma_k = 3.0 * sigma_n;
0062 
0063 % prepare points
0064 points = [XYZ(:,1),XYZ(:,2),XYZ(:,3),sigma_n*ones(N,1)];
0065 
0066 %%
0067 figure('name','Given points','color','w',...
0068     'Position',[0.1*ss(1),0.3*ss(2),0.3*ss(1),0.4*ss(2)]);
0069 plot3(points(:,1),points(:,2),points(:,3),'.k')
0070 axis equal;axis off;hold on
0071 xlabel('x');ylabel('y');zlabel('z');
0072 view([0,90])
0073 title('Given points')
0074 plot3([BB(1),BB(1),BB(3),BB(3),BB(1)],[BB(2),BB(4),BB(4),BB(2),BB(2)],[0,0,0,0,0],'-k');
0075 
0076 
0077 
0078 
0079 return