-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathincg_tri.c
160 lines (136 loc) · 4.51 KB
/
incg_tri.c
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
#include <stdio.h>
#include <stdlib.h>
#include <unistd.h>
#include <math.h>
#ifdef __cplusplus
extern "C" {
#endif
#include "incg_tri.h"
#include "incg_utils.h"
//
// Function to test whether a point falls within a triangle
// Assumes that all points lie (approximately) on the same plane.
//
int incg_Tri_PointInside(
const double p1[3],
const double p2[3],
const double p3[3],
const double xp[3] )
{
double dx1[3], dx2[3], proj;
#ifdef _DEBUG2_
FILE *fp=fopen("data.dat","w");
test_write_basis(fp);
fprintf(fp,"variables = x y z u v w \n");
fprintf(fp,"zone T=\"triangle\", N=3,E=1,F=FEPOINT,ET=TRIANGLE\n");
fprintf(fp,"%lf %lf %lf 0 0 0 \n", p1[0],p1[1],p1[2] );
fprintf(fp,"%lf %lf %lf 0 0 0 \n", p2[0],p2[1],p2[2] );
fprintf(fp,"%lf %lf %lf 0 0 0 \n", p3[0],p3[1],p3[2] );
fprintf(fp,"1 2 3 \n");
fprintf(fp,"variables = x y z u v w \n");
fprintf(fp,"zone T=\"point\", N=1,E=1,F=FEPOINT,ET=TRIANGLE\n");
fprintf(fp,"%lf %lf %lf 0 0 0 \n", xp[0],xp[1],xp[2] );
fprintf(fp,"1 1 1 \n");
#endif
// side of edge 1
dx1[0] = p2[0] - p1[0];
dx1[1] = p2[1] - p1[1];
dx1[2] = p2[2] - p1[2];
dx2[0] = p1[0] - p3[0];
dx2[1] = p1[1] - p3[1];
dx2[2] = p1[2] - p3[2];
incg_Vec_Normalize3( dx1 );
proj = incg_Vec_DotProduct( dx1, dx2 );
//printf("Projection %lf \n", proj );//HACK
dx2[0] = -(dx2[0] - proj*dx1[0]);
dx2[1] = -(dx2[1] - proj*dx1[1]);
dx2[2] = -(dx2[2] - proj*dx1[2]);
//printf("Vector perpend. to edge: %lf %lf %lf \n", dx2[0], dx2[1], dx2[2] );
#ifdef _DEBUG2_
fprintf(fp,"variables = x y z u v w \n");
fprintf(fp,"zone T=\"vector\", N=2,E=1,F=FEPOINT,ET=TRIANGLE\n");
fprintf(fp,"%lf %lf %lf 0 0 0 \n", 0.0, 0.0, 0.0 );
fprintf(fp,"%lf %lf %lf 0 0 0 \n", dx2[0],dx2[1],dx2[2] );
fprintf(fp,"1 2 1 \n");
#endif
proj = (xp[0] - p1[0])*dx2[0] +
(xp[1] - p1[1])*dx2[1] +
(xp[2] - p1[2])*dx2[2];
//printf("Projection two %lf \n", proj );//HACK
if( proj < 0.0 ) return(0);
// side of edge 2
dx1[0] = p3[0] - p2[0];
dx1[1] = p3[1] - p2[1];
dx1[2] = p3[2] - p2[2];
dx2[0] = p2[0] - p1[0];
dx2[1] = p2[1] - p1[1];
dx2[2] = p2[2] - p1[2];
incg_Vec_Normalize3( dx1 );
proj = incg_Vec_DotProduct( dx1, dx2 );
//printf("Projection %lf \n", proj );//HACK
dx2[0] = -(dx2[0] - proj*dx1[0]);
dx2[1] = -(dx2[1] - proj*dx1[1]);
dx2[2] = -(dx2[2] - proj*dx1[2]);
//printf("Vector perpend. to edge: %lf %lf %lf \n", dx2[0], dx2[1], dx2[2] );
#ifdef _DEBUG2_
fprintf(fp,"variables = x y z u v w \n");
fprintf(fp,"zone T=\"vector\", N=2,E=1,F=FEPOINT,ET=TRIANGLE\n");
fprintf(fp,"%lf %lf %lf 0 0 0 \n", 0.0, 0.0, 0.0 );
fprintf(fp,"%lf %lf %lf 0 0 0 \n", dx2[0],dx2[1],dx2[2] );
fprintf(fp,"1 2 1 \n");
#endif
proj = (xp[0] - p2[0])*dx2[0] +
(xp[1] - p2[1])*dx2[1] +
(xp[2] - p2[2])*dx2[2];
//printf("Projection two %lf \n", proj );//HACK
if( proj < 0.0 ) return(0);
// side of edge 2
dx1[0] = p1[0] - p3[0];
dx1[1] = p1[1] - p3[1];
dx1[2] = p1[2] - p3[2];
dx2[0] = p3[0] - p2[0];
dx2[1] = p3[1] - p2[1];
dx2[2] = p3[2] - p2[2];
incg_Vec_Normalize3( dx1 );
proj = incg_Vec_DotProduct( dx1, dx2 );
//printf("Projection %lf \n", proj );//HACK
dx2[0] = -(dx2[0] - proj*dx1[0]);
dx2[1] = -(dx2[1] - proj*dx1[1]);
dx2[2] = -(dx2[2] - proj*dx1[2]);
//printf("Vector perpend. to edge: %lf %lf %lf \n", dx2[0], dx2[1], dx2[2] );
#ifdef _DEBUG2_
fprintf(fp,"variables = x y z u v w \n");
fprintf(fp,"zone T=\"vector\", N=2,E=1,F=FEPOINT,ET=TRIANGLE\n");
fprintf(fp,"%lf %lf %lf 0 0 0 \n", 0.0, 0.0, 0.0 );
fprintf(fp,"%lf %lf %lf 0 0 0 \n", dx2[0],dx2[1],dx2[2] );
fprintf(fp,"1 2 1 \n");
#endif
proj = (xp[0] - p3[0])*dx2[0] +
(xp[1] - p3[1])*dx2[1] +
(xp[2] - p3[2])*dx2[2];
//printf("Projection two %lf \n", proj );//HACK
if( proj < 0.0 ) return(0);
return(1);
}
//
// Function to pick a random point inside the triangle with uniform probability
// given two numbers randomly drawn from a uniform distribution.
//
void incg_Tri_MakeRandomPoint(
const double p1[3],
const double p2[3],
const double p3[3],
double xp[3] )
{
double r0 = 1.0/((double) RAND_MAX); // should be static or a const
double r1 = ((double) rand())*r0;
double r2 = ((double) rand())*r0;
r1 = sqrt(r1);
double t1 = 1.0 - r1, t2 = r1*(1.0 - r2), t3 = r1*r2;
xp[0] = t1 * p1[0] + t2 * p2[0] + t3 * p3[0];
xp[1] = t1 * p1[1] + t2 * p2[1] + t3 * p3[1];
xp[2] = t1 * p1[0] + t2 * p2[2] + t3 * p3[2];
}
#ifdef __cplusplus
}
#endif