generate_points_on_plane

% generate_points_on_plane(r,A);

 N     = number of points
 A     = plane (assumed to be not vertical, A3\=0))
 sigma = standard deviation (in all coordinates)
 r     = [xmin,xmax,ymin,ymax]
 RL    = 0 left plane
       = 1 right plane
 rPCV  = 0 generate points
       = 1 take coordiantes from book

 X      = list of noisy coordinates (structs)

 Wolfgang Förstner
 wfoerstn@uni-bonn.de

0001 %% generate_points_on_plane(r,A);
0002 %
0003 % N     = number of points
0004 % A     = plane (assumed to be not vertical, A3\=0))
0005 % sigma = standard deviation (in all coordinates)
0006 % r     = [xmin,xmax,ymin,ymax]
0007 % RL    = 0 left plane
0008 %       = 1 right plane
0009 % rPCV  = 0 generate points
0010 %       = 1 take coordiantes from book
0011 %
0012 % X      = list of noisy coordinates (structs)
0013 %
0014 % Wolfgang Förstner
0015 % wfoerstn@uni-bonn.de
0016 
0017 function X = generate_points_on_plane(N,A,sigma,r,RL,rPCV)
0018 
0019 global plot_option
0020 
0021 % Margin
0022 m = 0.15;
0023 % AX+BY+CZ+D=0;
0024 Xlist = (1-2*m)*rand(N,1)*(r(2)-r(1))+r(1)+m*(r(2)-r(1));
0025 Ylist = (1-2*m)*rand(N,1)*(r(4)-r(3))+r(3)+m*(r(4)-r(3));
0026 Zlist = -(A(1)*Xlist + A(2)*Ylist + A(4)) / A(3);
0027 
0028 X_true = [Xlist,Ylist,Zlist];
0029 
0030 %% redundancy numbers for non-orthogonal fit
0031 % 1-1/N-xq/sumxq-yq/sumyq
0032 
0033 meanX = mean(Xlist);
0034 meanY = mean(Ylist);
0035 Xq = Xlist-meanX;
0036 Yq = Ylist-meanY;
0037 C = (N-1)*cov([Xq,Yq]);
0038 D = [Xq,Yq]*inv(C)*[Xq,Yq]';                                               %#ok<MINV>
0039 ri = ones(N,1)*(1-1/N)- diag(D);
0040 deltai = 4./sqrt(ri);
0041 sum(ri);
0042 disp(['Redundancy numbers  : ',num2str(ri')])
0043 disp(['Sensitivity factors : ',num2str(deltai')])
0044 
0045 %% noise
0046 if rPCV == 0
0047     Xe = X_true + rand_gauss([0,0,0]',sigma^2*eye(3),N)';
0048 else    
0049     XLRn=[-0.8749 3.8996 2.3502 0.5370 2.1825 2.6343;...
0050         -1.6218 3.2287 2.0395 0.8666 1.3679 2.2913;...
0051         -1.6059 1.6108 1.8873 0.9658 3.6694 2.4097;...
0052         -1.4275 0.0609 2.0249 1.0134 2.3504 2.2827;...
0053         -1.2100 0.8542 2.4765 -0.0517 2.4755 2.9952;...
0054         0 0 0                 1.6559 -0.0646 2.0355];
0055     
0056     if RL==0
0057         Xe = XLRn(1:5,1:3);
0058     else
0059         Xe = XLRn(1:6,4:6);
0060     end
0061 end
0062 
0063 X.h   = zeros(N,4);
0064 X.Crr = zeros(N,3,3);
0065 for n=1:N
0066     Xc = sugr_Point_3D(Xe(n,:)',sigma^2*eye(3));
0067     X.h(n,:)     = Xc.h';
0068     X.Crr(n,:,:) = Xc.Crr;
0069 end
0070 
0071 if plot_option > 0
0072    z = -(A(1)*r(1:2) + A(2)*r(3:4) + A(4)) / A(3);
0073    
0074    for n=1:N
0075        plot3(Xe(n,1),Xe(n,2),Xe(n,3),'.k','MarkerSize',15);
0076    end
0077    plot3([r(1),r(2),r(2),r(1),r(1)],[r(3),r(3),r(4),r(4),r(3)],[z(1),z(2),z(2),z(1),z(1)],'--k');
0078    xlabel('X')
0079    ylabel('Y')
0080    zlabel('Z')
0081    axis equal
0082 end