1 /*
2     Copyright (C) 1998  Dennis Roddeman
3     email: dennis.roddeman@feat.nl
4 
5     This program is free software; you can redistribute it and/or modify
6     it under the terms of the GNU General Public License as published by
7     the Free Software Foundation; either version 2 of the License, or
8     (at your option) any later version.
9 
10     This program is distributed in the hope that it will be useful,
11     but WITHOUT ANY WARRANTY; without even the implied warranty of
12     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
13     GNU General Public License for more details.
14 
15 
16     You should have received a copy of the GNU General Public License
17     along with this program; if not, write to the Free Software Foundation
18     59 Temple Place, Suite 330, Boston, MA, 02111-1307, USA
19 */
20 
21 #include "tochnog.h"
22 
intersect_line_with_line(double line_a0[],double line_a1[],double line_b0[],double line_b1[],double & iso_line_a,double & iso_line_b)23 long int intersect_line_with_line( double line_a0[], double line_a1[],
24   double line_b0[], double line_b1[], double &iso_line_a,
25   double &iso_line_b )
26 
27 {
28   double rdum=0., vec[MDIM], mat[MDIM*MDIM],
29     inv_mat[MDIM*MDIM], rh[MDIM], iso[MDIM];
30 
31   assert( ndim==2 );
32 
33   array_subtract( line_a1, line_a0, vec, ndim );
34   mat[0*ndim+0] = vec[0];
35   mat[1*ndim+0] = vec[1];
36   array_subtract( line_b1, line_b0, vec, ndim );
37   mat[0*ndim+1] = -vec[0];
38   mat[1*ndim+1] = -vec[1];
39   if ( matrix_inverse( mat, inv_mat, rdum, ndim ) ) {
40     array_subtract( line_b0, line_a0, rh, ndim );
41     matrix_ab( inv_mat, rh, iso, ndim, ndim, 1 );
42     iso_line_a = iso[0];
43     iso_line_b = iso[1];
44     return 1;
45   }
46   else
47     return 0;
48 }
49 
intersect_line_with_point(double line0[],double line1[],double point[],double & iso_line)50 long int intersect_line_with_point( double line0[], double line1[],
51   double point[], double &iso_line )
52 
53 {
54   assert( ndim==1 );
55 
56   iso_line = (point[0]-line0[0]) / (line1[0]-line0[0]);
57   return 1;
58 }
59 
intersect_line_with_triangle(double line0[],double line1[],double triangle0[],double triangle1[],double triangle2[],double & iso_line,double iso_triangle[])60 long int intersect_line_with_triangle( double line0[], double line1[],
61   double triangle0[], double triangle1[], double triangle2[],
62   double &iso_line, double iso_triangle[] )
63 
64 {
65   double rdum=0., vec[MDIM], mat[MDIM*MDIM],
66     inv_mat[MDIM*MDIM], rh[MDIM], iso[MDIM];
67 
68   assert( ndim==3 );
69 
70   array_subtract( line1, line0, vec, ndim );
71   mat[0*ndim+0] = vec[0];
72   mat[1*ndim+0] = vec[1];
73   mat[2*ndim+0] = vec[2];
74   array_subtract( triangle0, triangle2, vec, ndim );
75   mat[0*ndim+1] = -vec[0];
76   mat[1*ndim+1] = -vec[1];
77   mat[2*ndim+1] = -vec[2];
78   array_subtract( triangle1, triangle2, vec, ndim );
79   mat[0*ndim+2] = -vec[0];
80   mat[1*ndim+2] = -vec[1];
81   mat[2*ndim+2] = -vec[2];
82   if ( matrix_inverse( mat, inv_mat, rdum, ndim ) ) {
83     array_subtract( triangle2, line0, rh, ndim );
84     matrix_ab( inv_mat, rh, iso, ndim, ndim, 1 );
85     iso_line = iso[0];
86     iso_triangle[0] = iso[1];
87     iso_triangle[1] = iso[2];
88     return 1;
89   }
90   else
91     return 0;
92 }
93