demo_singularity_line_through_four_lines

% Test singularity of line through four lines
 Given four lines L_j,j=1,2,3,4, determine lines M_i passing these

 See PCV Section 7.2.1.3 Three and More Entities

 0 or 2 solutions

 Cases:
   random
   three parallel, one other line
   three intersecting lines, one other
   three coplanar intersecting lines, one other
   two intersecting lines, two others

 Output per case
   3D lines
   nullspace
   coefficients (should not be zero) see PCV (7.53)
   solutions (may be identical)
   two lines, if possible

 Wolfgang Förstner 10/2016
 wfoerstn@uni-bonn.de

 last changes: Susanne Wenzel 12/17
 wenzel@igg.uni-bonn.de

0001 %% Test singularity of line through four lines
0002 % Given four lines L_j,j=1,2,3,4, determine lines M_i passing these
0003 %
0004 % See PCV Section 7.2.1.3 Three and More Entities
0005 %
0006 % 0 or 2 solutions
0007 %
0008 % Cases:
0009 %   random
0010 %   three parallel, one other line
0011 %   three intersecting lines, one other
0012 %   three coplanar intersecting lines, one other
0013 %   two intersecting lines, two others
0014 %
0015 % Output per case
0016 %   3D lines
0017 %   nullspace
0018 %   coefficients (should not be zero) see PCV (7.53)
0019 %   solutions (may be identical)
0020 %   two lines, if possible
0021 %
0022 % Wolfgang Förstner 10/2016
0023 % wfoerstn@uni-bonn.de
0024 %
0025 % last changes: Susanne Wenzel 12/17
0026 % wenzel@igg.uni-bonn.de
0027 
0028 % clear all
0029 clearvars
0030 clc
0031 close all
0032 
0033 addpath(genpath('../General-Functions'))
0034 
0035 fprintf('\n------ DEMO Line Through Four Lines ------\n')
0036 
0037 fprintf('\nIn order to shorten the output, we display transposed vectors\n')
0038 
0039 
0040 % four random lines
0041 disp('------------- four random lines --------------')
0042 L1 = calc_Pi(rand(4,1))*randn(4,1);
0043 L2 = calc_Pi(rand(4,1))*randn(4,1);
0044 L3 = calc_Pi(rand(4,1))*randn(4,1);
0045 L4 = calc_Pi(rand(4,1))*randn(4,1);
0046 M = [L1,L2,L3,L4]                                                           %#ok<NOPTS>
0047 disp(['rank(M) = ', num2str(rank(M))])
0048 
0049 % determine polynomial for intersecting line
0050 fprintf('\ndetermine polynomial for intersecting line\n')
0051 NM = null(M'*calc_Dual);
0052 null_space_M = NM                                                           %#ok<NOPTS>
0053 
0054 N1 = NM(:,1);
0055 N2 = NM(:,2);
0056 a = (N1-N2)'*calc_Dual*(N1-N2);
0057 b = 2*(N1-N2)'*calc_Dual*N2;
0058 c = N2'*calc_Dual*N2;
0059 coefficients = [a,b,c]                                                       %#ok<NOPTS>
0060 rr = roots([a,b,c]);
0061 disp(['solutions: ',num2str(rr')])
0062 
0063 M1 = rr(1)*N1+(1-rr(1))*N2;
0064 M2 = rr(2)*N1+(1-rr(2))*N2;
0065 two_lines = [M1';M2'];
0066 fprintf('\ntwo lines \n')
0067 disp(['M_1 = [',num2str(M1'),']'])
0068 disp(['M_2 = [',num2str(M2'),']'])
0069 incidences = two_lines*calc_Dual*M                                          %#ok<NASGU,NOPTS>
0070 
0071 % three parallel
0072 fprintf('\n------------- three random but parallel lines, one other --------------\n')
0073 % randomly generate two directions
0074 D123 = [rand(3,1);0];
0075 D4 = [rand(3,1);0];
0076 % generate 3 3D lines having the same normal direction
0077 L1 = calc_Pi(rand(4,1))*D123;
0078 L2 = calc_Pi(rand(4,1))*D123;
0079 L3 = calc_Pi(rand(4,1))*D123;
0080 % generate another 3D line with different normal direction
0081 L4 = calc_Pi(rand(4,1))*D4;
0082 L1 = L1/L1(1);
0083 L2 = L2/L2(1);
0084 L3 = L3/L3(1);
0085 L4 = L4/L4(1);
0086 M = [L1,L2,L3,L4]                                                           %#ok<NOPTS>
0087 disp(['rank(M) = ', num2str(rank(M))])
0088 
0089 NM = null(M'*calc_Dual);
0090 null_space = NM                                                             %#ok<NASGU,NOPTS>
0091 N1 = NM(:,1);
0092 N2 = NM(:,2);
0093 a = (N1-N2)'*calc_Dual*(N1-N2);
0094 b = 2*(N1-N2)'*calc_Dual*N2;
0095 c = N2'*calc_Dual*N2;
0096 coefficients_polynomial_no_solution = [a,b,c]                               %#ok<NASGU,NOPTS>
0097 solutions = roots([a,b,c])';
0098 disp(['solutions: ',num2str(solutions)])
0099 
0100 % three through one point
0101 fprintf('\n------------- three random but intersecting lines, one other --------------\n')
0102 % define 2 3D points
0103 X1 = rand(4,1);
0104 X2 = rand(4,1);
0105 
0106 % generate random 3d lines passing through the same point X1
0107 L1 = calc_Pi(rand(4,1))*X1;
0108 L2 = calc_Pi(rand(4,1))*X1;
0109 L3 = calc_Pi(rand(4,1))*X1;
0110 % generate another random 3d lines passing through X2
0111 L4 = calc_Pi(rand(4,1))*X2;
0112 M = [L1,L2,L3,L4]                                                           %#ok<NOPTS>
0113 disp(['rank(M) = ', num2str(rank(M))])
0114 
0115 NM = null(M'*calc_Dual);
0116 nullspace = NM                                                              %#ok<NOPTS,NASGU>
0117 N1 = NM(:,1);
0118 N2 = NM(:,2);
0119 a = (N1-N2)'*calc_Dual*(N1-N2);
0120 b = 2*(N1-N2)'*calc_Dual*N2;
0121 c = N2'*calc_Dual*N2;
0122 coefficients_polynomial_no_solution = [a,b,c]                               %#ok<NOPTS,NASGU>
0123 solutions = roots([a,b,c])';
0124 disp(['solutions: ',num2str(solutions)])
0125 
0126 %% three coplanar three through one point
0127 fprintf('\n------------- three coplanar intersecting lines, one other --------------\n')
0128 % define three points on the plane
0129 X = [1,0,0,1]';
0130 Y = [0,1,0,1]';
0131 Z = [0,0,1,1]';
0132 % define another arbitrary point
0133 T = [randn(1,3),1]';
0134 % check whether they occasionally lie on a plane
0135 disp(['Determinant should not be zero: ', num2str(det([X,Y,Z,T]))])
0136 
0137 % 3D lines trough points on plane (X,Y,Z)
0138 L1 = calc_Pi(X)*Y;
0139 L2 = calc_Pi(Y)*Z;
0140 L3 = calc_Pi(Z)*X;
0141 % 3D line not on the same plane
0142 L4 = calc_Pi(X)*T;
0143 
0144 M = [L1,L2,L3,L4]                                                           %#ok<NOPTS>
0145 disp(['rank(M) = ', num2str(rank(M))])
0146 
0147 NM = null(M'*calc_Dual);
0148 nullspace = NM                                                              %#ok<NOPTS>
0149 N1 = NM(:,1);
0150 N2 = NM(:,2);
0151 as = (N1-N2)'*calc_Dual*(N1-N2);
0152 bs = 2*(N1-N2)'*calc_Dual*N2;
0153 cs = N2'*calc_Dual*N2;
0154 coefficients_polynomial_no_solution = [as,bs,cs]                            %#ok<NOPTS>
0155 solutions = roots([as,bs,cs])';
0156 disp(['solutions: ',num2str(solutions)])
0157 
0158 
0159 %% two intersecting lines
0160 fprintf('\n------------- two intersecting lines, two others --------------\n')
0161 % define 4 3D points as points on target lines
0162 X = [-1,0,2,1]';
0163 Y = [1,0,3,1]';
0164 Z = [0,1,2,1]';
0165 T = [1,-1,3,1]';
0166 % Line in direction of X-axis
0167 L1 = [1,0,0,0,0,0]';
0168 % Line in direction of Y-axis
0169 L2 = [0,1,0,0,0,0]';
0170 % Line through points X and Y, intersection the X-axis
0171 L3 = calc_Pi(X)*Y;
0172 % Line through points Z and T, intersecting the Y-axis
0173 L4 = calc_Pi(Z)*T;
0174 
0175 M = [L1,L2,L3,L4]                                                           %#ok<NOPTS>
0176 disp(['rank(M) = ', num2str(rank(M))])
0177 
0178 NM = null(M'*calc_Dual);
0179 null_space = NM                                                             %#ok<NOPTS>
0180 N1 = NM(:,1);
0181 N2 = NM(:,2);
0182 a = (N1-N2)'*calc_Dual*(N1-N2);
0183 b = 2*(N1-N2)'*calc_Dual*N2;
0184 c = N2'*calc_Dual*N2;
0185 coefficients_polynomial = [a,b,c]                                           %#ok<NOPTS>
0186 solutions = roots([a,b,c])';
0187 disp(['solutions: ',num2str(solutions)])
0188 
0189 M1 = rr(1)*N1+(1-rr(1))*N2;
0190 M2 = rr(2)*N1+(1-rr(2))*N2;
0191 two_lines = [M1';M2'];
0192 fprintf('\ntwo lines \n')
0193 disp(['M_1 = [',num2str(M1'),']'])
0194 disp(['M_2 = [',num2str(M2'),']'])
0195 incidences = two_lines*calc_Dual*M                                          %#ok<NOPTS>