69
71
float normal_quad_v3(float n[3], const float v1[3], const float v2[3], const float v3[3], const float v4[3])
82
n[0]= n1[1]*n2[2]-n1[2]*n2[1];
83
n[1]= n1[2]*n2[0]-n1[0]*n2[2];
84
n[2]= n1[0]*n2[1]-n1[1]*n2[0];
76
n1[0] = v1[0] - v3[0];
77
n1[1] = v1[1] - v3[1];
78
n1[2] = v1[2] - v3[2];
80
n2[0] = v2[0] - v4[0];
81
n2[1] = v2[1] - v4[1];
82
n2[2] = v2[2] - v4[2];
84
n[0] = n1[1] * n2[2] - n1[2] * n2[1];
85
n[1] = n1[2] * n2[0] - n1[0] * n2[2];
86
n[2] = n1[0] * n2[1] - n1[1] * n2[0];
86
88
return normalize_v3(n);
89
91
float area_tri_v2(const float v1[2], const float v2[2], const float v3[2])
91
return (float)(0.5f*fabs((v1[0]-v2[0])*(v2[1]-v3[1]) + (v1[1]-v2[1])*(v3[0]-v2[0])));
93
return 0.5f * fabsf((v1[0] - v2[0]) * (v2[1] - v3[1]) + (v1[1] - v2[1]) * (v3[0] - v2[0]));
94
96
float area_tri_signed_v2(const float v1[2], const float v2[2], const float v3[2])
96
return (float)(0.5f*((v1[0]-v2[0])*(v2[1]-v3[1]) + (v1[1]-v2[1])*(v3[0]-v2[0])));
98
return 0.5f * ((v1[0] - v2[0]) * (v2[1] - v3[1]) + (v1[1] - v2[1]) * (v3[0] - v2[0]));
99
float area_quad_v3(const float v1[3], const float v2[3], const float v3[3], const float v4[3]) /* only convex Quadrilaterals */
101
/* only convex Quadrilaterals */
102
float area_quad_v3(const float v1[3], const float v2[3], const float v3[3], const float v4[3])
101
104
float len, vec1[3], vec2[3], n[3];
103
106
sub_v3_v3v3(vec1, v2, v1);
104
107
sub_v3_v3v3(vec2, v4, v1);
105
108
cross_v3_v3v3(n, vec1, vec2);
106
len= normalize_v3(n);
109
len = normalize_v3(n);
108
111
sub_v3_v3v3(vec1, v4, v3);
109
112
sub_v3_v3v3(vec2, v2, v3);
110
113
cross_v3_v3v3(n, vec1, vec2);
111
len+= normalize_v3(n);
114
len += normalize_v3(n);
116
float area_tri_v3(const float v1[3], const float v2[3], const float v3[3]) /* Triangles */
120
float area_tri_v3(const float v1[3], const float v2[3], const float v3[3])
118
122
float len, vec1[3], vec2[3], n[3];
120
124
sub_v3_v3v3(vec1, v3, v2);
121
125
sub_v3_v3v3(vec2, v1, v2);
122
126
cross_v3_v3v3(n, vec1, vec2);
123
len= normalize_v3(n);
127
len = normalize_v3(n);
128
float area_poly_v3(int nr, float verts[][3], float *normal)
132
float area_poly_v3(int nr, float verts[][3], const float normal[3])
130
134
float x, y, z, area, max;
131
135
float *cur, *prev;
136
int a, px = 0, py = 1;
134
/* first: find dominant axis: 0==X, 1==Y, 2==Z */
135
x= (float)fabs(normal[0]);
136
y= (float)fabs(normal[1]);
137
z= (float)fabs(normal[2]);
138
/* first: find dominant axis: 0==X, 1==Y, 2==Z
139
* don't use 'axis_dominant_v3()' because we need max axis too */
140
x = fabsf(normal[0]);
141
y = fabsf(normal[1]);
142
z = fabsf(normal[2]);
138
143
max = MAX3(x, y, z);
144
if (max == y) py = 2;
145
150
/* The Trapezium Area Rule */
149
for(a=0; a<nr; a++) {
150
area+= (cur[px]-prev[px])*(cur[py]+prev[py]);
151
prev = verts[nr - 1];
154
for (a = 0; a < nr; a++) {
155
area += (cur[px] - prev[px]) * (cur[py] + prev[py]);
155
return (float)fabs(0.5*area/max);
160
return fabsf(0.5f * area / max);
158
163
/********************************* Distance **********************************/
160
/* distance v1 to line v2-v3 */
161
/* using Hesse formula, NO LINE PIECE! */
162
float dist_to_line_v2(float *v1, float *v2, float *v3)
165
/* distance v1 to line v2-v3
166
* using Hesse formula, NO LINE PIECE! */
167
float dist_to_line_v2(const float v1[2], const float v2[2], const float v3[2])
168
deler= (float)sqrt(a[0]*a[0]+a[1]*a[1]);
169
if(deler== 0.0f) return 0;
171
return (float)(fabs((v1[0]-v2[0])*a[0]+(v1[1]-v2[1])*a[1])/deler);
171
a[0] = v2[1] - v3[1];
172
a[1] = v3[0] - v2[0];
173
deler = (float)sqrt(a[0] * a[0] + a[1] * a[1]);
174
if (deler == 0.0f) return 0;
176
return fabsf((v1[0] - v2[0]) * a[0] + (v1[1] - v2[1]) * a[1]) / deler;
175
180
/* distance v1 to line-piece v2-v3 */
176
float dist_to_line_segment_v2(float *v1, float *v2, float *v3)
181
float dist_to_line_segment_v2(const float v1[2], const float v2[2], const float v3[2])
178
183
float labda, rc[2], pt[2], len;
182
len= rc[0]*rc[0]+ rc[1]*rc[1];
186
return (float)(sqrt(rc[0]*rc[0]+ rc[1]*rc[1]));
189
labda= (rc[0]*(v1[0]-v2[0]) + rc[1]*(v1[1]-v2[1]))/len;
194
else if(labda>=1.0) {
185
rc[0] = v3[0] - v2[0];
186
rc[1] = v3[1] - v2[1];
187
len = rc[0] * rc[0] + rc[1] * rc[1];
189
rc[0] = v1[0] - v2[0];
190
rc[1] = v1[1] - v2[1];
191
return (float)(sqrt(rc[0] * rc[0] + rc[1] * rc[1]));
194
labda = (rc[0] * (v1[0] - v2[0]) + rc[1] * (v1[1] - v2[1])) / len;
199
else if (labda >= 1.0f) {
199
pt[0]= labda*rc[0]+v2[0];
200
pt[1]= labda*rc[1]+v2[1];
204
pt[0] = labda * rc[0] + v2[0];
205
pt[1] = labda * rc[1] + v2[1];
205
return (float)sqrt(rc[0]*rc[0]+ rc[1]*rc[1]);
208
rc[0] = pt[0] - v1[0];
209
rc[1] = pt[1] - v1[1];
210
return sqrtf(rc[0] * rc[0] + rc[1] * rc[1]);
213
/* point closest to v1 on line v2-v3 in 2D */
214
void closest_to_line_segment_v2(float close_r[2], const float p[2], const float l1[2], const float l2[2])
218
lambda = closest_to_line_v2(cp, p, l1, l2);
221
copy_v2_v2(close_r, l1);
222
else if (lambda >= 1.0f)
223
copy_v2_v2(close_r, l2);
225
copy_v2_v2(close_r, cp);
208
228
/* point closest to v1 on line v2-v3 in 3D */
209
void closest_to_line_segment_v3(float *closest, float v1[3], float v2[3], float v3[3])
229
void closest_to_line_segment_v3(float close_r[3], const float v1[3], const float v2[3], const float v3[3])
211
231
float lambda, cp[3];
213
lambda= closest_to_line_v3(cp,v1, v2, v3);
233
lambda = closest_to_line_v3(cp, v1, v2, v3);
216
copy_v3_v3(closest, v2);
217
else if(lambda >= 1.0f)
218
copy_v3_v3(closest, v3);
236
copy_v3_v3(close_r, v2);
237
else if (lambda >= 1.0f)
238
copy_v3_v3(close_r, v3);
220
copy_v3_v3(closest, cp);
240
copy_v3_v3(close_r, cp);
243
/* find the closest point on a plane to another point and store it in close_r
244
* close_r: return coordinate
245
* plane_co: a point on the plane
246
* plane_no_unit: the plane's normal, and d is the last number in the plane equation 0 = ax + by + cz + d
247
* pt: the point that you want the nearest of
250
void closest_to_plane_v3(float close_r[3], const float plane_co[3], const float plane_no_unit[3], const float pt[3])
255
sub_v3_v3v3(temp, pt, plane_co);
256
dotprod = dot_v3v3(temp, plane_no_unit);
258
close_r[0] = pt[0] - (plane_no_unit[0] * dotprod);
259
close_r[1] = pt[1] - (plane_no_unit[1] * dotprod);
260
close_r[2] = pt[2] - (plane_no_unit[2] * dotprod);
263
/* signed distance from the point to the plane in 3D */
264
float dist_to_plane_normalized_v3(const float p[3], const float plane_co[3], const float plane_no_unit[3])
266
float plane_co_other[3];
268
add_v3_v3v3(plane_co_other, plane_co, plane_no_unit);
270
return line_point_factor_v3(p, plane_co, plane_co_other);
273
float dist_to_plane_v3(const float p[3], const float plane_co[3], const float plane_no[3])
275
float plane_no_unit[3];
276
float plane_co_other[3];
278
normalize_v3_v3(plane_no_unit, plane_no);
279
add_v3_v3v3(plane_co_other, plane_co, plane_no_unit);
281
return line_point_factor_v3(p, plane_co, plane_co_other);
223
284
/* distance v1 to line-piece v2-v3 in 3D */
224
float dist_to_line_segment_v3(float *v1, float *v2, float *v3)
285
float dist_to_line_segment_v3(const float v1[3], const float v2[3], const float v3[3])
226
287
float closest[3];
233
294
/******************************* Intersection ********************************/
235
296
/* intersect Line-Line, shorts */
236
int isect_line_line_v2_short(short *v1, short *v2, short *v3, short *v4)
297
int isect_line_line_v2_int(const int v1[2], const int v2[2], const int v3[2], const int v4[2])
240
0: no intersection of segments
241
1: exact intersection of segments
242
2: cross-intersection of segments
244
299
float div, labda, mu;
246
div= (float)((v2[0]-v1[0])*(v4[1]-v3[1])-(v2[1]-v1[1])*(v4[0]-v3[0]));
247
if(div==0.0f) return -1;
249
labda= ((float)(v1[1]-v3[1])*(v4[0]-v3[0])-(v1[0]-v3[0])*(v4[1]-v3[1]))/div;
251
mu= ((float)(v1[1]-v3[1])*(v2[0]-v1[0])-(v1[0]-v3[0])*(v2[1]-v1[1]))/div;
253
if(labda>=0.0f && labda<=1.0f && mu>=0.0f && mu<=1.0f) {
254
if(labda==0.0f || labda==1.0f || mu==0.0f || mu==1.0f) return 1;
301
div = (float)((v2[0] - v1[0]) * (v4[1] - v3[1]) - (v2[1] - v1[1]) * (v4[0] - v3[0]));
302
if (div == 0.0f) return ISECT_LINE_LINE_COLINEAR;
304
labda = ((float)(v1[1] - v3[1]) * (v4[0] - v3[0]) - (v1[0] - v3[0]) * (v4[1] - v3[1])) / div;
306
mu = ((float)(v1[1] - v3[1]) * (v2[0] - v1[0]) - (v1[0] - v3[0]) * (v2[1] - v1[1])) / div;
308
if (labda >= 0.0f && labda <= 1.0f && mu >= 0.0f && mu <= 1.0f) {
309
if (labda == 0.0f || labda == 1.0f || mu == 0.0f || mu == 1.0f) return ISECT_LINE_LINE_EXACT;
310
return ISECT_LINE_LINE_CROSS;
312
return ISECT_LINE_LINE_NONE;
260
315
/* intersect Line-Line, floats */
261
int isect_line_line_v2(float *v1, float *v2, float *v3, float *v4)
316
int isect_line_line_v2(const float v1[2], const float v2[2], const float v3[2], const float v4[2])
265
0: no intersection of segments
266
1: exact intersection of segments
267
2: cross-intersection of segments
269
318
float div, labda, mu;
271
div= (v2[0]-v1[0])*(v4[1]-v3[1])-(v2[1]-v1[1])*(v4[0]-v3[0]);
272
if(div==0.0) return -1;
274
labda= ((float)(v1[1]-v3[1])*(v4[0]-v3[0])-(v1[0]-v3[0])*(v4[1]-v3[1]))/div;
276
mu= ((float)(v1[1]-v3[1])*(v2[0]-v1[0])-(v1[0]-v3[0])*(v2[1]-v1[1]))/div;
278
if(labda>=0.0 && labda<=1.0 && mu>=0.0 && mu<=1.0) {
279
if(labda==0.0 || labda==1.0 || mu==0.0 || mu==1.0) return 1;
320
div = (v2[0] - v1[0]) * (v4[1] - v3[1]) - (v2[1] - v1[1]) * (v4[0] - v3[0]);
321
if (div == 0.0f) return ISECT_LINE_LINE_COLINEAR;
323
labda = ((float)(v1[1] - v3[1]) * (v4[0] - v3[0]) - (v1[0] - v3[0]) * (v4[1] - v3[1])) / div;
325
mu = ((float)(v1[1] - v3[1]) * (v2[0] - v1[0]) - (v1[0] - v3[0]) * (v2[1] - v1[1])) / div;
327
if (labda >= 0.0f && labda <= 1.0f && mu >= 0.0f && mu <= 1.0f) {
328
if (labda == 0.0f || labda == 1.0f || mu == 0.0f || mu == 1.0f) return ISECT_LINE_LINE_EXACT;
329
return ISECT_LINE_LINE_CROSS;
331
return ISECT_LINE_LINE_NONE;
334
/* get intersection point of two 2D segments and return intersection type:
338
int isect_seg_seg_v2_point(const float v1[2], const float v2[2], const float v3[2], const float v4[2], float vi[2])
340
float a1, a2, b1, b2, c1, c2, d;
342
const float eps = 0.000001f;
352
d = a1 * b2 - a2 * b1;
355
if (a1 * c2 - a2 * c1 == 0.0f && b1 * c2 - b2 * c1 == 0.0f) { /* equal lines */
356
float a[2], b[2], c[2];
359
if (len_v2v2(v1, v2) == 0.0f) {
360
if (len_v2v2(v3, v4) > eps) {
361
/* use non-point segment as basis */
362
SWAP(const float *, v1, v3);
363
SWAP(const float *, v2, v4);
365
else { /* both of segments are points */
366
if (equals_v2v2(v1, v3)) { /* points are equal */
371
/* two different points */
376
sub_v2_v2v2(a, v3, v1);
377
sub_v2_v2v2(b, v2, v1);
378
sub_v2_v2v2(c, v2, v1);
379
u = dot_v2v2(a, b) / dot_v2v2(c, c);
381
sub_v2_v2v2(a, v4, v1);
382
u2 = dot_v2v2(a, b) / dot_v2v2(c, c);
384
if (u > u2) SWAP(float, u, u2);
386
if (u > 1.0f + eps || u2 < -eps) return -1; /* non-ovlerlapping segments */
387
else if (maxf(0.0f, u) == minf(1.0f, u2)) { /* one common point: can return result */
388
interp_v2_v2v2(vi, v1, v2, maxf(0, u));
393
/* lines are colliniar */
397
u = (c2 * b1 - b2 * c1) / d;
398
v = (c1 * a2 - a1 * c2) / d;
400
if (u >= -eps && u <= 1.0f + eps && v >= -eps && v <= 1.0f + eps) { /* intersection */
401
interp_v2_v2v2(vi, v1, v2, u);
405
/* out of segment intersection */
409
int isect_line_sphere_v3(const float l1[3], const float l2[3],
410
const float sp[3], const float r,
411
float r_p1[3], float r_p2[3])
413
/* l1: coordinates (point of line)
414
* l2: coordinates (point of line)
415
* sp, r: coordinates and radius (sphere)
416
* r_p1, r_p2: return intersection coordinates
420
/* adapted for use in blender by Campbell Barton - 2011
422
* atelier iebele abel - 2001
424
* http://www.iebele.nl
426
* sphere_line_intersection function adapted from:
427
* http://astronomy.swin.edu.au/pbourke/geometry/sphereline
428
* Paul Bourke pbourke@swin.edu.au
431
const float ldir[3] = {
437
const float a = dot_v3v3(ldir, ldir);
439
const float b = 2.0f *
440
(ldir[0] * (l1[0] - sp[0]) +
441
ldir[1] * (l1[1] - sp[1]) +
442
ldir[2] * (l1[2] - sp[2]));
447
(2.0f * dot_v3v3(sp, l1)) -
450
const float i = b * b - 4.0f * a * c;
455
/* no intersections */
458
else if (i == 0.0f) {
459
/* one intersection */
460
mu = -b / (2.0f * a);
461
madd_v3_v3v3fl(r_p1, l1, ldir, mu);
465
const float i_sqrt = sqrt(i); /* avoid calc twice */
467
/* first intersection */
468
mu = (-b + i_sqrt) / (2.0f * a);
469
madd_v3_v3v3fl(r_p1, l1, ldir, mu);
471
/* second intersection */
472
mu = (-b - i_sqrt) / (2.0f * a);
473
madd_v3_v3v3fl(r_p2, l1, ldir, mu);
477
/* math domain error - nan */
482
/* keep in sync with isect_line_sphere_v3 */
483
int isect_line_sphere_v2(const float l1[2], const float l2[2],
484
const float sp[2], const float r,
485
float r_p1[2], float r_p2[2])
487
const float ldir[2] = {l2[0] - l1[0],
490
const float a = dot_v2v2(ldir, ldir);
492
const float b = 2.0f *
493
(ldir[0] * (l1[0] - sp[0]) +
494
ldir[1] * (l1[1] - sp[1]));
499
(2.0f * dot_v2v2(sp, l1)) -
502
const float i = b * b - 4.0f * a * c;
507
/* no intersections */
510
else if (i == 0.0f) {
511
/* one intersection */
512
mu = -b / (2.0f * a);
513
madd_v2_v2v2fl(r_p1, l1, ldir, mu);
517
const float i_sqrt = sqrt(i); /* avoid calc twice */
519
/* first intersection */
520
mu = (-b + i_sqrt) / (2.0f * a);
521
madd_v2_v2v2fl(r_p1, l1, ldir, mu);
523
/* second intersection */
524
mu = (-b - i_sqrt) / (2.0f * a);
525
madd_v2_v2v2fl(r_p2, l1, ldir, mu);
529
/* math domain error - nan */
290
static short IsectLLPt2Df(float x0,float y0,float x1,float y1,
291
float x2,float y2,float x3,float y3, float *xi,float *yi)
538
static short IsectLLPt2Df(const float x0, const float y0, const float x1, const float y1,
539
const float x2, const float y2, const float x3, const float y3, float *xi, float *yi)
295
this function computes the intersection of the sent lines
296
and returns the intersection point, note that the function assumes
297
the lines intersect. the function can handle vertical as well
298
as horizontal lines. note the function isn't very clever, it simply
299
applies the math, but we don't need speed since this is a
302
float c1,c2, // constants of linear equations
303
det_inv, // the inverse of the determinant of the coefficient
304
m1,m2; // the slopes of each line
543
* this function computes the intersection of the sent lines
544
* and returns the intersection point, note that the function assumes
545
* the lines intersect. the function can handle vertical as well
546
* as horizontal lines. note the function isn't very clever, it simply
547
* applies the math, but we don't need speed since this is a
548
* pre-processing step
550
float c1, c2; /* constants of linear equations */
551
float det_inv; /* the inverse of the determinant of the coefficient */
552
float m1, m2; /* the slopes of each line */
306
compute slopes, note the cludge for infinity, however, this will
309
if (fabs(x1-x0) > 0.000001)
310
m1 = (y1-y0) / (x1-x0);
312
return -1; /*m1 = (float) 1e+10;*/ // close enough to infinity
314
if (fabs(x3-x2) > 0.000001)
315
m2 = (y3-y2) / (x3-x2);
317
return -1; /*m2 = (float) 1e+10;*/ // close enough to infinity
319
if (fabs(m1-m2) < 0.000001)
320
return -1; /* parallel lines */
327
// compute the inverse of the determinate
554
* compute slopes, note the cludge for infinity, however, this will
557
if (fabs(x1 - x0) > 0.000001)
558
m1 = (y1 - y0) / (x1 - x0);
560
return -1; /*m1 = (float)1e+10;*/ // close enough to infinity
562
if (fabs(x3 - x2) > 0.000001)
563
m2 = (y3 - y2) / (x3 - x2);
565
return -1; /*m2 = (float)1e+10;*/ // close enough to infinity
567
if (fabs(m1 - m2) < 0.000001)
568
return -1; /* parallel lines */
575
// compute the inverse of the determinate
329
577
det_inv = 1.0f / (-m1 + m2);
331
// use Kramers rule to compute xi and yi
333
*xi= ((-c2 + c1) *det_inv);
334
*yi= ((m2*c1 - m1*c2) *det_inv);
337
} // end Intersect_Lines
339
#define SIDE_OF_LINE(pa,pb,pp) ((pa[0]-pp[0])*(pb[1]-pp[1]))-((pb[0]-pp[0])*(pa[1]-pp[1]))
579
// use Kramers rule to compute xi and yi
581
*xi = ((-c2 + c1) * det_inv);
582
*yi = ((m2 * c1 - m1 * c2) * det_inv);
340
587
/* point in tri */
341
// XXX was called IsectPT2Df
342
int isect_point_tri_v2(float pt[2], float v1[2], float v2[2], float v3[2])
589
int isect_point_tri_v2(const float pt[2], const float v1[2], const float v2[2], const float v3[2])
344
if (SIDE_OF_LINE(v1,v2,pt)>=0.0) {
345
if (SIDE_OF_LINE(v2,v3,pt)>=0.0) {
346
if (SIDE_OF_LINE(v3,v1,pt)>=0.0) {
591
if (line_point_side_v2(v1, v2, pt) >= 0.0f) {
592
if (line_point_side_v2(v2, v3, pt) >= 0.0f) {
593
if (line_point_side_v2(v3, v1, pt) >= 0.0f) {
351
if (! (SIDE_OF_LINE(v2,v3,pt)>=0.0)) {
352
if (! (SIDE_OF_LINE(v3,v1,pt)>=0.0)) {
599
if (!(line_point_side_v2(v2, v3, pt) >= 0.0f)) {
600
if (!(line_point_side_v2(v3, v1, pt) >= 0.0f)) {
360
609
/* point in quad - only convex quads */
361
int isect_point_quad_v2(float pt[2], float v1[2], float v2[2], float v3[2], float v4[2])
610
int isect_point_quad_v2(const float pt[2], const float v1[2], const float v2[2], const float v3[2], const float v4[2])
363
if (SIDE_OF_LINE(v1,v2,pt)>=0.0) {
364
if (SIDE_OF_LINE(v2,v3,pt)>=0.0) {
365
if (SIDE_OF_LINE(v3,v4,pt)>=0.0) {
366
if (SIDE_OF_LINE(v4,v1,pt)>=0.0) {
612
if (line_point_side_v2(v1, v2, pt) >= 0.0f) {
613
if (line_point_side_v2(v2, v3, pt) >= 0.0f) {
614
if (line_point_side_v2(v3, v4, pt) >= 0.0f) {
615
if (line_point_side_v2(v4, v1, pt) >= 0.0f) {
372
if (! (SIDE_OF_LINE(v2,v3,pt)>=0.0)) {
373
if (! (SIDE_OF_LINE(v3,v4,pt)>=0.0)) {
374
if (! (SIDE_OF_LINE(v4,v1,pt)>=0.0)) {
622
if (!(line_point_side_v2(v2, v3, pt) >= 0.0f)) {
623
if (!(line_point_side_v2(v3, v4, pt) >= 0.0f)) {
624
if (!(line_point_side_v2(v4, v1, pt) >= 0.0f)) {
384
634
/* moved from effect.c
385
test if the line starting at p1 ending at p2 intersects the triangle v0..v2
386
return non zero if it does
388
int isect_line_tri_v3(float p1[3], float p2[3], float v0[3], float v1[3], float v2[3], float *lambda, float *uv)
635
* test if the line starting at p1 ending at p2 intersects the triangle v0..v2
636
* return non zero if it does
638
int isect_line_tri_v3(const float p1[3], const float p2[3],
639
const float v0[3], const float v1[3], const float v2[3],
640
float *r_lambda, float r_uv[2])
391
643
float p[3], s[3], d[3], e1[3], e2[3], q[3];
392
644
float a, f, u, v;
394
646
sub_v3_v3v3(e1, v1, v0);
395
647
sub_v3_v3v3(e2, v2, v0);
396
648
sub_v3_v3v3(d, p2, p1);
398
650
cross_v3_v3v3(p, d, e2);
399
651
a = dot_v3v3(e1, p);
400
if ((a > -0.000001) && (a < 0.000001)) return 0;
652
if ((a > -0.000001f) && (a < 0.000001f)) return 0;
403
655
sub_v3_v3v3(s, p1, v0);
405
657
u = f * dot_v3v3(s, p);
406
if ((u < 0.0)||(u > 1.0)) return 0;
658
if ((u < 0.0f) || (u > 1.0f)) return 0;
408
660
cross_v3_v3v3(q, s, e1);
410
v = f * dot_v3v3(d, q);
411
if ((v < 0.0)||((u + v) > 1.0)) return 0;
413
*lambda = f * dot_v3v3(e2, q);
414
if ((*lambda < 0.0)||(*lambda > 1.0)) return 0;
662
v = f * dot_v3v3(d, q);
663
if ((v < 0.0f) || ((u + v) > 1.0f)) return 0;
665
*r_lambda = f * dot_v3v3(e2, q);
666
if ((*r_lambda < 0.0f) || (*r_lambda > 1.0f)) return 0;
424
676
/* moved from effect.c
425
test if the ray starting at p1 going in d direction intersects the triangle v0..v2
426
return non zero if it does
428
int isect_ray_tri_v3(float p1[3], float d[3], float v0[3], float v1[3], float v2[3], float *lambda, float *uv)
430
float p[3], s[3], e1[3], e2[3], q[3];
433
sub_v3_v3v3(e1, v1, v0);
434
sub_v3_v3v3(e2, v2, v0);
436
cross_v3_v3v3(p, d, e2);
438
/* note: these values were 0.000001 in 2.4x but for projection snapping on
439
* a human head (1BU==1m), subsurf level 2, this gave many errors - campbell */
440
if ((a > -0.00000001) && (a < 0.00000001)) return 0;
443
sub_v3_v3v3(s, p1, v0);
445
u = f * dot_v3v3(s, p);
446
if ((u < 0.0)||(u > 1.0)) return 0;
448
cross_v3_v3v3(q, s, e1);
450
v = f * dot_v3v3(d, q);
451
if ((v < 0.0)||((u + v) > 1.0)) return 0;
453
*lambda = f * dot_v3v3(e2, q);
454
if ((*lambda < 0.0)) return 0;
464
int isect_ray_tri_epsilon_v3(float p1[3], float d[3], float v0[3], float v1[3], float v2[3], float *lambda, float *uv, float epsilon)
466
float p[3], s[3], e1[3], e2[3], q[3];
469
sub_v3_v3v3(e1, v1, v0);
470
sub_v3_v3v3(e2, v2, v0);
472
cross_v3_v3v3(p, d, e2);
474
if (a == 0.0f) return 0;
477
sub_v3_v3v3(s, p1, v0);
479
u = f * dot_v3v3(s, p);
480
if ((u < -epsilon)||(u > 1.0f+epsilon)) return 0;
482
cross_v3_v3v3(q, s, e1);
484
v = f * dot_v3v3(d, q);
485
if ((v < -epsilon)||((u + v) > 1.0f+epsilon)) return 0;
487
*lambda = f * dot_v3v3(e2, q);
488
if ((*lambda < 0.0f)) return 0;
498
int isect_ray_tri_threshold_v3(float p1[3], float d[3], float v0[3], float v1[3], float v2[3], float *lambda, float *uv, float threshold)
677
* test if the ray starting at p1 going in d direction intersects the triangle v0..v2
678
* return non zero if it does
680
int isect_ray_tri_v3(const float p1[3], const float d[3],
681
const float v0[3], const float v1[3], const float v2[3],
682
float *r_lambda, float r_uv[2])
684
float p[3], s[3], e1[3], e2[3], q[3];
687
sub_v3_v3v3(e1, v1, v0);
688
sub_v3_v3v3(e2, v2, v0);
690
cross_v3_v3v3(p, d, e2);
692
/* note: these values were 0.000001 in 2.4x but for projection snapping on
693
* a human head (1BU==1m), subsurf level 2, this gave many errors - campbell */
694
if ((a > -0.00000001f) && (a < 0.00000001f)) return 0;
697
sub_v3_v3v3(s, p1, v0);
699
u = f * dot_v3v3(s, p);
700
if ((u < 0.0f) || (u > 1.0f)) return 0;
702
cross_v3_v3v3(q, s, e1);
704
v = f * dot_v3v3(d, q);
705
if ((v < 0.0f) || ((u + v) > 1.0f)) return 0;
707
*r_lambda = f * dot_v3v3(e2, q);
708
if ((*r_lambda < 0.0f)) return 0;
718
int isect_ray_plane_v3(const float p1[3], const float d[3],
719
const float v0[3], const float v1[3], const float v2[3],
720
float *r_lambda, const int clip)
722
float p[3], s[3], e1[3], e2[3], q[3];
724
/* float u, v; */ /*UNUSED*/
726
sub_v3_v3v3(e1, v1, v0);
727
sub_v3_v3v3(e2, v2, v0);
729
cross_v3_v3v3(p, d, e2);
731
/* note: these values were 0.000001 in 2.4x but for projection snapping on
732
* a human head (1BU==1m), subsurf level 2, this gave many errors - campbell */
733
if ((a > -0.00000001f) && (a < 0.00000001f)) return 0;
736
sub_v3_v3v3(s, p1, v0);
738
/* u = f * dot_v3v3(s, p); */ /*UNUSED*/
740
cross_v3_v3v3(q, s, e1);
742
/* v = f * dot_v3v3(d, q); */ /*UNUSED*/
744
*r_lambda = f * dot_v3v3(e2, q);
745
if (clip && (*r_lambda < 0.0f)) return 0;
750
int isect_ray_tri_epsilon_v3(const float p1[3], const float d[3],
751
const float v0[3], const float v1[3], const float v2[3],
752
float *r_lambda, float uv[2], const float epsilon)
754
float p[3], s[3], e1[3], e2[3], q[3];
757
sub_v3_v3v3(e1, v1, v0);
758
sub_v3_v3v3(e2, v2, v0);
760
cross_v3_v3v3(p, d, e2);
762
if (a == 0.0f) return 0;
765
sub_v3_v3v3(s, p1, v0);
767
u = f * dot_v3v3(s, p);
768
if ((u < -epsilon) || (u > 1.0f + epsilon)) return 0;
770
cross_v3_v3v3(q, s, e1);
772
v = f * dot_v3v3(d, q);
773
if ((v < -epsilon) || ((u + v) > 1.0f + epsilon)) return 0;
775
*r_lambda = f * dot_v3v3(e2, q);
776
if ((*r_lambda < 0.0f)) return 0;
786
int isect_ray_tri_threshold_v3(const float p1[3], const float d[3],
787
const float v0[3], const float v1[3], const float v2[3],
788
float *r_lambda, float r_uv[2], const float threshold)
500
790
float p[3], s[3], e1[3], e2[3], q[3];
501
791
float a, f, u, v;
502
792
float du = 0, dv = 0;
504
794
sub_v3_v3v3(e1, v1, v0);
505
795
sub_v3_v3v3(e2, v2, v0);
507
797
cross_v3_v3v3(p, d, e2);
508
798
a = dot_v3v3(e1, p);
509
if ((a > -0.000001) && (a < 0.000001)) return 0;
799
if ((a > -0.000001f) && (a < 0.000001f)) return 0;
512
802
sub_v3_v3v3(s, p1, v0);
514
804
cross_v3_v3v3(q, s, e1);
515
*lambda = f * dot_v3v3(e2, q);
516
if ((*lambda < 0.0)) return 0;
805
*r_lambda = f * dot_v3v3(e2, q);
806
if ((*r_lambda < 0.0f)) return 0;
518
808
u = f * dot_v3v3(s, p);
519
809
v = f * dot_v3v3(d, q);
521
811
if (u < 0) du = u;
522
812
if (u > 1) du = u - 1;
523
813
if (v < 0) dv = v;
524
814
if (v > 1) dv = v - 1;
525
if (u > 0 && v > 0 && u + v > 1)
815
if (u > 0 && v > 0 && u + v > 1) {
527
816
float t = u + v - 1;
532
821
mul_v3_fl(e1, du);
533
822
mul_v3_fl(e2, dv);
535
if (dot_v3v3(e1, e1) + dot_v3v3(e2, e2) > threshold * threshold)
824
if (dot_v3v3(e1, e1) + dot_v3v3(e2, e2) > threshold * threshold) {
836
int isect_line_plane_v3(float out[3],
837
const float l1[3], const float l2[3],
838
const float plane_co[3], const float plane_no[3], const short no_flip)
840
float l_vec[3]; /* l1 -> l2 normalized vector */
841
float p_no[3]; /* 'plane_no' normalized */
844
sub_v3_v3v3(l_vec, l2, l1);
847
normalize_v3_v3(p_no, plane_no);
849
dot = dot_v3v3(l_vec, p_no);
854
float l1_plane[3]; /* line point aligned with the plane */
855
float dist; /* 'plane_no' aligned distance to the 'plane_co' */
857
/* for predictable flipping since the plane is only used to
858
* define a direction, ignore its flipping and aligned with 'l_vec' */
864
add_v3_v3v3(l1_plane, l1, p_no);
866
dist = line_point_factor_v3(plane_co, l1, l1_plane);
868
/* treat line like a ray, when 'no_flip' is set */
869
if (no_flip && dist < 0.0f) {
873
mul_v3_fl(l_vec, dist / dot);
875
add_v3_v3v3(out, l1, l_vec);
881
/* note: return normal isn't unit length */
882
void isect_plane_plane_v3(float r_isect_co[3], float r_isect_no[3],
883
const float plane_a_co[3], const float plane_a_no[3],
884
const float plane_b_co[3], const float plane_b_no[3])
886
float plane_a_co_other[3];
887
cross_v3_v3v3(r_isect_no, plane_a_no, plane_b_no); /* direction is simply the cross product */
888
cross_v3_v3v3(plane_a_co_other, plane_a_no, r_isect_no);
889
add_v3_v3(plane_a_co_other, plane_a_co);
890
isect_line_plane_v3(r_isect_co, plane_a_co, plane_a_co_other, plane_b_co, plane_b_no, FALSE);
549
894
/* Adapted from the paper by Kasper Fauerby */
550
896
/* "Improved Collision detection and Response" */
551
static int getLowestRoot(float a, float b, float c, float maxR, float* root)
897
static int getLowestRoot(const float a, const float b, const float c, const float maxR, float *root)
553
899
// Check if a solution exists
554
float determinant = b*b - 4.0f*a*c;
900
float determinant = b * b - 4.0f * a * c;
556
902
// If determinant is negative it means no solutions.
557
if (determinant >= 0.0f)
903
if (determinant >= 0.0f) {
559
904
// calculate the two roots: (if determinant == 0 then
560
// x1==x2 but let’s disregard that slight optimization)
905
// x1==x2 but lets disregard that slight optimization)
561
906
float sqrtD = (float)sqrt(determinant);
562
float r1 = (-b - sqrtD) / (2.0f*a);
563
float r2 = (-b + sqrtD) / (2.0f*a);
907
float r1 = (-b - sqrtD) / (2.0f * a);
908
float r2 = (-b + sqrtD) / (2.0f * a);
565
910
// Sort so x1 <= x2
567
912
SWAP(float, r1, r2);
569
914
// Get lowest root:
570
if (r1 > 0.0f && r1 < maxR)
915
if (r1 > 0.0f && r1 < maxR) {
576
920
// It is possible that we want x2 - this can happen
578
if (r2 > 0.0f && r2 < maxR)
922
if (r2 > 0.0f && r2 < maxR) {
631
973
/*---test inside of tri---*/
632
974
/* plane intersection point */
634
point[0] = p1[0] + vel[0]*t0 - nor[0]*radius;
635
point[1] = p1[1] + vel[1]*t0 - nor[1]*radius;
636
point[2] = p1[2] + vel[2]*t0 - nor[2]*radius;
976
point[0] = p1[0] + vel[0] * t0 - nor[0] * radius;
977
point[1] = p1[1] + vel[1] * t0 - nor[1] * radius;
978
point[2] = p1[2] + vel[2] * t0 - nor[2] * radius;
639
981
/* is the point in the tri? */
644
sub_v3_v3v3(temp,point,v0);
653
if(z <= 0.0f && (x >= 0.0f && y >= 0.0f))
655
//(((unsigned int)z)& ~(((unsigned int)x)|((unsigned int)y))) & 0x80000000){
657
copy_v3_v3(ipoint,point);
982
a = dot_v3v3(e1, e1);
983
b = dot_v3v3(e1, e2);
984
c = dot_v3v3(e2, e2);
986
sub_v3_v3v3(temp, point, v0);
987
d = dot_v3v3(temp, e1);
988
e = dot_v3v3(temp, e2);
992
z = x + y - (a * c - b * b);
995
if (z <= 0.0f && (x >= 0.0f && y >= 0.0f)) {
996
//(((unsigned int)z)& ~(((unsigned int)x)|((unsigned int)y))) & 0x80000000) {
998
copy_v3_v3(ipoint, point);
665
/*---test points---*/
666
a=vel2=dot_v3v3(vel,vel);
1006
/*---test points---*/
1007
a = dot_v3v3(vel, vel);
669
sub_v3_v3v3(temp,p1,v0);
670
b=2.0f*dot_v3v3(vel,temp);
671
c=dot_v3v3(temp,temp)-radius2;
1010
sub_v3_v3v3(temp, p1, v0);
1011
b = 2.0f * dot_v3v3(vel, temp);
1012
c = dot_v3v3(temp, temp) - radius2;
673
if(getLowestRoot(a, b, c, *lambda, lambda))
675
copy_v3_v3(ipoint,v0);
1014
if (getLowestRoot(a, b, c, *r_lambda, r_lambda)) {
1015
copy_v3_v3(ipoint, v0);
680
sub_v3_v3v3(temp,p1,v1);
681
b=2.0f*dot_v3v3(vel,temp);
682
c=dot_v3v3(temp,temp)-radius2;
1020
sub_v3_v3v3(temp, p1, v1);
1021
b = 2.0f * dot_v3v3(vel, temp);
1022
c = dot_v3v3(temp, temp) - radius2;
684
if(getLowestRoot(a, b, c, *lambda, lambda))
686
copy_v3_v3(ipoint,v1);
1024
if (getLowestRoot(a, b, c, *r_lambda, r_lambda)) {
1025
copy_v3_v3(ipoint, v1);
691
sub_v3_v3v3(temp,p1,v2);
692
b=2.0f*dot_v3v3(vel,temp);
693
c=dot_v3v3(temp,temp)-radius2;
1030
sub_v3_v3v3(temp, p1, v2);
1031
b = 2.0f * dot_v3v3(vel, temp);
1032
c = dot_v3v3(temp, temp) - radius2;
695
if(getLowestRoot(a, b, c, *lambda, lambda))
697
copy_v3_v3(ipoint,v2);
1034
if (getLowestRoot(a, b, c, *r_lambda, r_lambda)) {
1035
copy_v3_v3(ipoint, v2);
702
sub_v3_v3v3(e3,v2,v1); //wasnt yet calculated
1039
/*---test edges---*/
1040
sub_v3_v3v3(e3, v2, v1); //wasnt yet calculated
706
sub_v3_v3v3(bv,v0,p1);
708
elen2 = dot_v3v3(e1,e1);
709
edotv = dot_v3v3(e1,vel);
710
edotbv = dot_v3v3(e1,bv);
712
a=elen2*(-dot_v3v3(vel,vel))+edotv*edotv;
713
b=2.0f*(elen2*dot_v3v3(vel,bv)-edotv*edotbv);
714
c=elen2*(radius2-dot_v3v3(bv,bv))+edotbv*edotbv;
716
if(getLowestRoot(a, b, c, *lambda, &newLambda))
718
e=(edotv*newLambda-edotbv)/elen2;
720
if(e >= 0.0f && e <= 1.0f)
723
copy_v3_v3(ipoint,e1);
1044
sub_v3_v3v3(bv, v0, p1);
1046
elen2 = dot_v3v3(e1, e1);
1047
edotv = dot_v3v3(e1, vel);
1048
edotbv = dot_v3v3(e1, bv);
1050
a = elen2 * (-dot_v3v3(vel, vel)) + edotv * edotv;
1051
b = 2.0f * (elen2 * dot_v3v3(vel, bv) - edotv * edotbv);
1052
c = elen2 * (radius2 - dot_v3v3(bv, bv)) + edotbv * edotbv;
1054
if (getLowestRoot(a, b, c, *r_lambda, &newLambda)) {
1055
e = (edotv * newLambda - edotbv) / elen2;
1057
if (e >= 0.0f && e <= 1.0f) {
1058
*r_lambda = newLambda;
1059
copy_v3_v3(ipoint, e1);
1060
mul_v3_fl(ipoint, e);
725
1061
add_v3_v3(ipoint, v0);
732
elen2 = dot_v3v3(e2,e2);
733
edotv = dot_v3v3(e2,vel);
734
edotbv = dot_v3v3(e2,bv);
736
a=elen2*(-dot_v3v3(vel,vel))+edotv*edotv;
737
b=2.0f*(elen2*dot_v3v3(vel,bv)-edotv*edotbv);
738
c=elen2*(radius2-dot_v3v3(bv,bv))+edotbv*edotbv;
740
if(getLowestRoot(a, b, c, *lambda, &newLambda))
742
e=(edotv*newLambda-edotbv)/elen2;
744
if(e >= 0.0f && e <= 1.0f)
747
copy_v3_v3(ipoint,e2);
1068
elen2 = dot_v3v3(e2, e2);
1069
edotv = dot_v3v3(e2, vel);
1070
edotbv = dot_v3v3(e2, bv);
1072
a = elen2 * (-dot_v3v3(vel, vel)) + edotv * edotv;
1073
b = 2.0f * (elen2 * dot_v3v3(vel, bv) - edotv * edotbv);
1074
c = elen2 * (radius2 - dot_v3v3(bv, bv)) + edotbv * edotbv;
1076
if (getLowestRoot(a, b, c, *r_lambda, &newLambda)) {
1077
e = (edotv * newLambda - edotbv) / elen2;
1079
if (e >= 0.0f && e <= 1.0f) {
1080
*r_lambda = newLambda;
1081
copy_v3_v3(ipoint, e2);
1082
mul_v3_fl(ipoint, e);
749
1083
add_v3_v3(ipoint, v0);
755
sub_v3_v3v3(bv,v0,p1);
756
elen2 = dot_v3v3(e1,e1);
757
edotv = dot_v3v3(e1,vel);
758
edotbv = dot_v3v3(e1,bv);
760
sub_v3_v3v3(bv,v1,p1);
761
elen2 = dot_v3v3(e3,e3);
762
edotv = dot_v3v3(e3,vel);
763
edotbv = dot_v3v3(e3,bv);
765
a=elen2*(-dot_v3v3(vel,vel))+edotv*edotv;
766
b=2.0f*(elen2*dot_v3v3(vel,bv)-edotv*edotbv);
767
c=elen2*(radius2-dot_v3v3(bv,bv))+edotbv*edotbv;
769
if(getLowestRoot(a, b, c, *lambda, &newLambda))
771
e=(edotv*newLambda-edotbv)/elen2;
773
if(e >= 0.0f && e <= 1.0f)
776
copy_v3_v3(ipoint,e3);
1089
/* sub_v3_v3v3(bv,v0,p1); */ /* UNUSED */
1090
/* elen2 = dot_v3v3(e1,e1); */ /* UNUSED */
1091
/* edotv = dot_v3v3(e1,vel); */ /* UNUSED */
1092
/* edotbv = dot_v3v3(e1,bv); */ /* UNUSED */
1094
sub_v3_v3v3(bv, v1, p1);
1095
elen2 = dot_v3v3(e3, e3);
1096
edotv = dot_v3v3(e3, vel);
1097
edotbv = dot_v3v3(e3, bv);
1099
a = elen2 * (-dot_v3v3(vel, vel)) + edotv * edotv;
1100
b = 2.0f * (elen2 * dot_v3v3(vel, bv) - edotv * edotbv);
1101
c = elen2 * (radius2 - dot_v3v3(bv, bv)) + edotbv * edotbv;
1103
if (getLowestRoot(a, b, c, *r_lambda, &newLambda)) {
1104
e = (edotv * newLambda - edotbv) / elen2;
1106
if (e >= 0.0f && e <= 1.0f) {
1107
*r_lambda = newLambda;
1108
copy_v3_v3(ipoint, e3);
1109
mul_v3_fl(ipoint, e);
778
1110
add_v3_v3(ipoint, v1);
784
1116
return found_by_sweep;
786
int isect_axial_line_tri_v3(int axis, float p1[3], float p2[3], float v0[3], float v1[3], float v2[3], float *lambda)
1119
int isect_axial_line_tri_v3(const int axis, const float p1[3], const float p2[3],
1120
const float v0[3], const float v1[3], const float v2[3], float *r_lambda)
788
1122
float p[3], e1[3], e2[3];
790
int a0=axis, a1=(axis+1)%3, a2=(axis+2)%3;
1124
int a0 = axis, a1 = (axis + 1) % 3, a2 = (axis + 2) % 3;
792
1126
//return isect_line_tri_v3(p1,p2,v0,v1,v2,lambda);
794
1128
///* first a simple bounding box test */
795
//if(MIN3(v0[a1],v1[a1],v2[a1]) > p1[a1]) return 0;
796
//if(MIN3(v0[a2],v1[a2],v2[a2]) > p1[a2]) return 0;
797
//if(MAX3(v0[a1],v1[a1],v2[a1]) < p1[a1]) return 0;
798
//if(MAX3(v0[a2],v1[a2],v2[a2]) < p1[a2]) return 0;
1129
//if (MIN3(v0[a1],v1[a1],v2[a1]) > p1[a1]) return 0;
1130
//if (MIN3(v0[a2],v1[a2],v2[a2]) > p1[a2]) return 0;
1131
//if (MAX3(v0[a1],v1[a1],v2[a1]) < p1[a1]) return 0;
1132
//if (MAX3(v0[a2],v1[a2],v2[a2]) < p1[a2]) return 0;
800
1134
///* then a full intersection test */
802
sub_v3_v3v3(e1,v1,v0);
803
sub_v3_v3v3(e2,v2,v0);
804
sub_v3_v3v3(p,v0,p1);
806
f= (e2[a1]*e1[a2]-e2[a2]*e1[a1]);
807
if ((f > -0.000001) && (f < 0.000001)) return 0;
809
v= (p[a2]*e1[a1]-p[a1]*e1[a2])/f;
810
if ((v < 0.0)||(v > 1.0)) return 0;
813
if((f > -0.000001) && (f < 0.000001)){
815
if((f > -0.000001) && (f < 0.000001)) return 0;
816
u= (-p[a2]-v*e2[a2])/f;
1136
sub_v3_v3v3(e1, v1, v0);
1137
sub_v3_v3v3(e2, v2, v0);
1138
sub_v3_v3v3(p, v0, p1);
1140
f = (e2[a1] * e1[a2] - e2[a2] * e1[a1]);
1141
if ((f > -0.000001f) && (f < 0.000001f)) return 0;
1143
v = (p[a2] * e1[a1] - p[a1] * e1[a2]) / f;
1144
if ((v < 0.0f) || (v > 1.0f)) return 0;
1147
if ((f > -0.000001f) && (f < 0.000001f)) {
1149
if ((f > -0.000001f) && (f < 0.000001f)) return 0;
1150
u = (-p[a2] - v * e2[a2]) / f;
819
u= (-p[a1]-v*e2[a1])/f;
821
if ((u < 0.0)||((u + v) > 1.0)) return 0;
823
*lambda = (p[a0]+u*e1[a0]+v*e2[a0])/(p2[a0]-p1[a0]);
825
if ((*lambda < 0.0)||(*lambda > 1.0)) return 0;
1153
u = (-p[a1] - v * e2[a1]) / f;
1155
if ((u < 0.0f) || ((u + v) > 1.0f)) return 0;
1157
*r_lambda = (p[a0] + u * e1[a0] + v * e2[a0]) / (p2[a0] - p1[a0]);
1159
if ((*r_lambda < 0.0f) || (*r_lambda > 1.0f)) return 0;
932
1259
f1 = dot_v3v3(cb, ab) / dot_v3v3(ab, ab);
933
1260
f2 = dot_v3v3(ca, ab) / dot_v3v3(ab, ab);
935
1262
if (f1 >= 0 && f1 <= 1 &&
938
1265
mul_v3_fl(a, f1);
939
1266
add_v3_v3v3(vi, v1, a);
1268
if (r_lambda) *r_lambda = f1;
946
1270
return 1; /* intersection found */
959
int isect_aabb_aabb_v3(float min1[3], float max1[3], float min2[3], float max2[3])
1281
int isect_aabb_aabb_v3(const float min1[3], const float max1[3], const float min2[3], const float max2[3])
961
return (min1[0]<max2[0] && min1[1]<max2[1] && min1[2]<max2[2] &&
962
min2[0]<max1[0] && min2[1]<max1[1] && min2[2]<max1[2]);
1283
return (min1[0] < max2[0] && min1[1] < max2[1] && min1[2] < max2[2] &&
1284
min2[0] < max1[0] && min2[1] < max1[1] && min2[2] < max1[2]);
965
1287
/* find closest point to p on line through l1,l2 and return lambda,
966
1288
* where (0 <= lambda <= 1) when cp is in the line segement l1,l2
968
float closest_to_line_v3(float cp[3],float p[3], float l1[3], float l2[3])
1290
float closest_to_line_v3(float cp[3], const float p[3], const float l1[3], const float l2[3])
970
float h[3],u[3],lambda;
1292
float h[3], u[3], lambda;
971
1293
sub_v3_v3v3(u, l2, l1);
972
1294
sub_v3_v3v3(h, p, l1);
973
lambda =dot_v3v3(u,h)/dot_v3v3(u,u);
1295
lambda = dot_v3v3(u, h) / dot_v3v3(u, u);
974
1296
cp[0] = l1[0] + u[0] * lambda;
975
1297
cp[1] = l1[1] + u[1] * lambda;
976
1298
cp[2] = l1[2] + u[2] * lambda;
1302
float closest_to_line_v2(float cp[2], const float p[2], const float l1[2], const float l2[2])
1304
float h[2], u[2], lambda;
1305
sub_v2_v2v2(u, l2, l1);
1306
sub_v2_v2v2(h, p, l1);
1307
lambda = dot_v2v2(u, h) / dot_v2v2(u, u);
1308
cp[0] = l1[0] + u[0] * lambda;
1309
cp[1] = l1[1] + u[1] * lambda;
981
1313
/* little sister we only need to know lambda */
982
static float lambda_cp_line(float p[3], float l1[3], float l2[3])
1314
float line_point_factor_v3(const float p[3], const float l1[3], const float l2[3])
985
1317
sub_v3_v3v3(u, l2, l1);
986
1318
sub_v3_v3v3(h, p, l1);
987
return(dot_v3v3(u,h)/dot_v3v3(u,u));
1319
return (dot_v3v3(u, h) / dot_v3v3(u, u));
1322
float line_point_factor_v2(const float p[2], const float l1[2], const float l2[2])
1325
sub_v2_v2v2(u, l2, l1);
1326
sub_v2_v2v2(h, p, l1);
1327
return (dot_v2v2(u, h) / dot_v2v2(u, u));
1330
/* ensyre the distance between these points is no greater then 'dist'
1331
* if it is, scale then both into the center */
1332
void limit_dist_v3(float v1[3], float v2[3], const float dist)
1334
const float dist_old = len_v3v3(v1, v2);
1336
if (dist_old > dist) {
1339
float fac = (dist / dist_old) * 0.5f;
1341
copy_v3_v3(v1_old, v1);
1342
copy_v3_v3(v2_old, v2);
1344
interp_v3_v3v3(v1, v1_old, v2_old, 0.5f - fac);
1345
interp_v3_v3v3(v2, v1_old, v2_old, 0.5f + fac);
991
1349
/* Similar to LineIntersectsTriangleUV, except it operates on a quad and in 2d, assumes point is in quad */
992
void isect_point_quad_uv_v2(float v0[2], float v1[2], float v2[2], float v3[2], float pt[2], float *uv)
1350
void isect_point_quad_uv_v2(const float v0[2], const float v1[2], const float v2[2], const float v3[2],
1351
const float pt[2], float r_uv[2])
994
float x0,y0, x1,y1, wtot, v2d[2], w1, w2;
1353
float x0, y0, x1, y1, wtot, v2d[2], w1, w2;
996
1355
/* used for parallel lines */
997
1356
float pt3d[3], l1[3], l2[3], pt_on_line[3];
999
1358
/* compute 2 edges of the quad intersection point */
1000
if (IsectLLPt2Df(v0[0],v0[1],v1[0],v1[1], v2[0],v2[1],v3[0],v3[1], &x0,&y0) == 1) {
1359
if (IsectLLPt2Df(v0[0], v0[1], v1[0], v1[1], v2[0], v2[1], v3[0], v3[1], &x0, &y0) == 1) {
1001
1360
/* the intersection point between the quad-edge intersection and the point in the quad we want the uv's for */
1002
1361
/* should never be paralle !! */
1003
1362
/*printf("\tnot parallel 1\n");*/
1004
IsectLLPt2Df(pt[0],pt[1],x0,y0, v0[0],v0[1],v3[0],v3[1], &x1,&y1);
1363
IsectLLPt2Df(pt[0], pt[1], x0, y0, v0[0], v0[1], v3[0], v3[1], &x1, &y1);
1006
1365
/* Get the weights from the new intersection point, to each edge */
1366
v2d[0] = x1 - v0[0];
1367
v2d[1] = y1 - v0[1];
1009
1368
w1 = len_v2(v2d);
1011
v2d[0] = x1-v3[0]; /* some but for the other vert */
1370
v2d[0] = x1 - v3[0]; /* some but for the other vert */
1371
v2d[1] = y1 - v3[1];
1013
1372
w2 = len_v2(v2d);
1015
1374
/*w1 = w1/wtot;*/
1016
1375
/*w2 = w2/wtot;*/
1376
r_uv[0] = w1 / wtot;
1019
1379
/* lines are parallel, lambda_cp_line_ex is 3d grrr */
1020
1380
/*printf("\tparallel1\n");*/
1021
1381
pt3d[0] = pt[0];
1022
1382
pt3d[1] = pt[1];
1023
1383
pt3d[2] = l1[2] = l2[2] = 0.0f;
1025
l1[0] = v0[0]; l1[1] = v0[1];
1026
l2[0] = v1[0]; l2[1] = v1[1];
1027
closest_to_line_v3(pt_on_line,pt3d, l1, l2);
1028
v2d[0] = pt[0]-pt_on_line[0]; /* same, for the other vert */
1029
v2d[1] = pt[1]-pt_on_line[1];
1389
closest_to_line_v3(pt_on_line, pt3d, l1, l2);
1390
v2d[0] = pt[0] - pt_on_line[0]; /* same, for the other vert */
1391
v2d[1] = pt[1] - pt_on_line[1];
1030
1392
w1 = len_v2(v2d);
1032
l1[0] = v2[0]; l1[1] = v2[1];
1033
l2[0] = v3[0]; l2[1] = v3[1];
1034
closest_to_line_v3(pt_on_line,pt3d, l1, l2);
1035
v2d[0] = pt[0]-pt_on_line[0]; /* same, for the other vert */
1036
v2d[1] = pt[1]-pt_on_line[1];
1398
closest_to_line_v3(pt_on_line, pt3d, l1, l2);
1399
v2d[0] = pt[0] - pt_on_line[0]; /* same, for the other vert */
1400
v2d[1] = pt[1] - pt_on_line[1];
1037
1401
w2 = len_v2(v2d);
1403
r_uv[0] = w1 / wtot;
1042
1406
/* Same as above to calc the uv[1] value, alternate calculation */
1044
if (IsectLLPt2Df(v0[0],v0[1],v3[0],v3[1], v1[0],v1[1],v2[0],v2[1], &x0,&y0) == 1) { /* was v0,v1 v2,v3 now v0,v3 v1,v2*/
1408
if (IsectLLPt2Df(v0[0], v0[1], v3[0], v3[1], v1[0], v1[1], v2[0], v2[1], &x0, &y0) == 1) { /* was v0,v1 v2,v3 now v0,v3 v1,v2*/
1045
1409
/* never paralle if above was not */
1046
1410
/*printf("\tnot parallel2\n");*/
1047
IsectLLPt2Df(pt[0],pt[1],x0,y0, v0[0],v0[1],v1[0],v1[1], &x1,&y1);/* was v0,v3 now v0,v1*/
1411
IsectLLPt2Df(pt[0], pt[1], x0, y0, v0[0], v0[1], v1[0], v1[1], &x1, &y1); /* was v0,v3 now v0,v1*/
1413
v2d[0] = x1 - v0[0];
1414
v2d[1] = y1 - v0[1];
1051
1415
w1 = len_v2(v2d);
1417
v2d[0] = x1 - v1[0];
1418
v2d[1] = y1 - v1[1];
1055
1419
w2 = len_v2(v2d);
1421
r_uv[1] = w1 / wtot;
1059
1424
/* lines are parallel, lambda_cp_line_ex is 3d grrr */
1060
1425
/*printf("\tparallel2\n");*/
1061
1426
pt3d[0] = pt[0];
1062
1427
pt3d[1] = pt[1];
1063
1428
pt3d[2] = l1[2] = l2[2] = 0.0f;
1066
l1[0] = v0[0]; l1[1] = v0[1];
1067
l2[0] = v3[0]; l2[1] = v3[1];
1068
closest_to_line_v3(pt_on_line,pt3d, l1, l2);
1069
v2d[0] = pt[0]-pt_on_line[0]; /* some but for the other vert */
1070
v2d[1] = pt[1]-pt_on_line[1];
1435
closest_to_line_v3(pt_on_line, pt3d, l1, l2);
1436
v2d[0] = pt[0] - pt_on_line[0]; /* some but for the other vert */
1437
v2d[1] = pt[1] - pt_on_line[1];
1071
1438
w1 = len_v2(v2d);
1073
l1[0] = v1[0]; l1[1] = v1[1];
1074
l2[0] = v2[0]; l2[1] = v2[1];
1075
closest_to_line_v3(pt_on_line,pt3d, l1, l2);
1076
v2d[0] = pt[0]-pt_on_line[0]; /* some but for the other vert */
1077
v2d[1] = pt[1]-pt_on_line[1];
1444
closest_to_line_v3(pt_on_line, pt3d, l1, l2);
1445
v2d[0] = pt[0] - pt_on_line[0]; /* some but for the other vert */
1446
v2d[1] = pt[1] - pt_on_line[1];
1078
1447
w2 = len_v2(v2d);
1449
r_uv[1] = w1 / wtot;
1082
1451
/* may need to flip UV's here */
1085
1454
/* same as above but does tri's and quads, tri's are a bit of a hack */
1086
void isect_point_face_uv_v2(int isquad, float v0[2], float v1[2], float v2[2], float v3[2], float pt[2], float *uv)
1455
void isect_point_face_uv_v2(const int isquad,
1456
const float v0[2], const float v1[2], const float v2[2], const float v3[2],
1457
const float pt[2], float r_uv[2])
1089
isect_point_quad_uv_v2(v0, v1, v2, v3, pt, uv);
1460
isect_point_quad_uv_v2(v0, v1, v2, v3, pt, r_uv);
1092
1463
/* not for quads, use for our abuse of LineIntersectsTriangleUV */
1093
1464
float p1_3d[3], p2_3d[3], v0_3d[3], v1_3d[3], v2_3d[3], lambda;
1095
p1_3d[0] = p2_3d[0] = uv[0];
1096
p1_3d[1] = p2_3d[1] = uv[1];
1466
p1_3d[0] = p2_3d[0] = r_uv[0];
1467
p1_3d[1] = p2_3d[1] = r_uv[1];
1097
1468
p1_3d[2] = 1.0f;
1098
1469
p2_3d[2] = -1.0f;
1099
1470
v0_3d[2] = v1_3d[2] = v2_3d[2] = 0.0;
1101
1472
/* generate a new fuv, (this is possibly a non optimal solution,
1102
1473
* since we only need 2d calculation but use 3d func's)
1104
1475
* this method makes an imaginary triangle in 2d space using the UV's from the derived mesh face
1105
1476
* Then find new uv coords using the fuv and this face with LineIntersectsTriangleUV.
1106
* This means the new values will be correct in relation to the derived meshes face.
1477
* This means the new values will be correct in relation to the derived meshes face.
1108
1479
copy_v2_v2(v0_3d, v0);
1109
1480
copy_v2_v2(v1_3d, v1);
1110
1481
copy_v2_v2(v2_3d, v2);
1112
1483
/* Doing this in 3D is not nice */
1113
isect_line_tri_v3(p1_3d, p2_3d, v0_3d, v1_3d, v2_3d, &lambda, uv);
1484
isect_line_tri_v3(p1_3d, p2_3d, v0_3d, v1_3d, v2_3d, &lambda, r_uv);
1117
1488
#if 0 // XXX this version used to be used in isect_point_tri_v2_int() and was called IsPointInTri2D
1118
1490
int isect_point_tri_v2(float pt[2], float v1[2], float v2[2], float v3[2])
1120
1492
float inp1, inp2, inp3;
1122
inp1= (v2[0]-v1[0])*(v1[1]-pt[1]) + (v1[1]-v2[1])*(v1[0]-pt[0]);
1123
inp2= (v3[0]-v2[0])*(v2[1]-pt[1]) + (v2[1]-v3[1])*(v2[0]-pt[0]);
1124
inp3= (v1[0]-v3[0])*(v3[1]-pt[1]) + (v3[1]-v1[1])*(v3[0]-pt[0]);
1126
if(inp1<=0.0f && inp2<=0.0f && inp3<=0.0f) return 1;
1127
if(inp1>=0.0f && inp2>=0.0f && inp3>=0.0f) return 1;
1494
inp1 = (v2[0] - v1[0]) * (v1[1] - pt[1]) + (v1[1] - v2[1]) * (v1[0] - pt[0]);
1495
inp2 = (v3[0] - v2[0]) * (v2[1] - pt[1]) + (v2[1] - v3[1]) * (v2[0] - pt[0]);
1496
inp3 = (v1[0] - v3[0]) * (v3[1] - pt[1]) + (v3[1] - v1[1]) * (v3[0] - pt[0]);
1498
if (inp1 <= 0.0f && inp2 <= 0.0f && inp3 <= 0.0f) return 1;
1499
if (inp1 >= 0.0f && inp2 >= 0.0f && inp3 >= 0.0f) return 1;
1134
1507
int isect_point_tri_v2(float v0[2], float v1[2], float v2[2], float pt[2])
1136
/* not for quads, use for our abuse of LineIntersectsTriangleUV */
1137
float p1_3d[3], p2_3d[3], v0_3d[3], v1_3d[3], v2_3d[3];
1139
float lambda, uv[3];
1141
p1_3d[0] = p2_3d[0] = uv[0]= pt[0];
1142
p1_3d[1] = p2_3d[1] = uv[1]= uv[2]= pt[1];
1145
v0_3d[2] = v1_3d[2] = v2_3d[2] = 0.0;
1147
/* generate a new fuv, (this is possibly a non optimal solution,
1148
* since we only need 2d calculation but use 3d func's)
1150
* this method makes an imaginary triangle in 2d space using the UV's from the derived mesh face
1151
* Then find new uv coords using the fuv and this face with LineIntersectsTriangleUV.
1152
* This means the new values will be correct in relation to the derived meshes face.
1154
copy_v2_v2(v0_3d, v0);
1155
copy_v2_v2(v1_3d, v1);
1156
copy_v2_v2(v2_3d, v2);
1158
/* Doing this in 3D is not nice */
1159
return isect_line_tri_v3(p1_3d, p2_3d, v0_3d, v1_3d, v2_3d, &lambda, uv);
1509
/* not for quads, use for our abuse of LineIntersectsTriangleUV */
1510
float p1_3d[3], p2_3d[3], v0_3d[3], v1_3d[3], v2_3d[3];
1512
float lambda, uv[3];
1514
p1_3d[0] = p2_3d[0] = uv[0] = pt[0];
1515
p1_3d[1] = p2_3d[1] = uv[1] = uv[2] = pt[1];
1518
v0_3d[2] = v1_3d[2] = v2_3d[2] = 0.0;
1520
/* generate a new fuv, (this is possibly a non optimal solution,
1521
* since we only need 2d calculation but use 3d func's)
1523
* this method makes an imaginary triangle in 2d space using the UV's from the derived mesh face
1524
* Then find new uv coords using the fuv and this face with LineIntersectsTriangleUV.
1525
* This means the new values will be correct in relation to the derived meshes face.
1527
copy_v2_v2(v0_3d, v0);
1528
copy_v2_v2(v1_3d, v1);
1529
copy_v2_v2(v2_3d, v2);
1531
/* Doing this in 3D is not nice */
1532
return isect_line_tri_v3(p1_3d, p2_3d, v0_3d, v1_3d, v2_3d, &lambda, uv);
1172
int isect_point_tri_v2_int(int x1, int y1, int x2, int y2, int a, int b)
1543
int isect_point_tri_v2_int(const int x1, const int y1, const int x2, const int y2, const int a, const int b)
1174
1545
float v1[2], v2[2], v3[2], p[2];
1188
1559
return isect_point_tri_v2(p, v1, v2, v3);
1191
static int point_in_slice(float p[3], float v1[3], float l1[3], float l2[3])
1562
static int point_in_slice(const float p[3], const float v1[3], const float l1[3], const float l2[3])
1196
a line including l1,l2 and a point not on the line
1197
define a subset of R3 delimeted by planes parallel to the line and orthogonal
1198
to the (point --> line) distance vector,one plane on the line one on the point,
1199
the room inside usually is rather small compared to R3 though still infinte
1200
useful for restricting (speeding up) searches
1201
e.g. all points of triangular prism are within the intersection of 3 'slices'
1202
onother trivial case : cube
1203
but see a 'spat' which is a deformed cube with paired parallel planes needs only 3 slices too
1205
float h,rp[3],cp[3],q[3];
1207
closest_to_line_v3(cp,v1,l1,l2);
1208
sub_v3_v3v3(q,cp,v1);
1210
sub_v3_v3v3(rp,p,v1);
1211
h=dot_v3v3(q,rp)/dot_v3v3(q,q);
1567
* a line including l1,l2 and a point not on the line
1568
* define a subset of R3 delimited by planes parallel to the line and orthogonal
1569
* to the (point --> line) distance vector,one plane on the line one on the point,
1570
* the room inside usually is rather small compared to R3 though still infinte
1571
* useful for restricting (speeding up) searches
1572
* e.g. all points of triangular prism are within the intersection of 3 'slices'
1573
* onother trivial case : cube
1574
* but see a 'spat' which is a deformed cube with paired parallel planes needs only 3 slices too
1576
float h, rp[3], cp[3], q[3];
1578
closest_to_line_v3(cp, v1, l1, l2);
1579
sub_v3_v3v3(q, cp, v1);
1581
sub_v3_v3v3(rp, p, v1);
1582
h = dot_v3v3(q, rp) / dot_v3v3(q, q);
1212
1583
if (h < 0.0f || h > 1.0f) return 0;
1217
/*adult sister defining the slice planes by the origin and the normal
1218
NOTE |normal| may not be 1 but defining the thickness of the slice*/
1219
static int point_in_slice_as(float p[3],float origin[3],float normal[3])
1589
/* adult sister defining the slice planes by the origin and the normal
1590
* NOTE |normal| may not be 1 but defining the thickness of the slice */
1591
static int point_in_slice_as(float p[3], float origin[3], float normal[3])
1222
sub_v3_v3v3(rp,p,origin);
1223
h=dot_v3v3(normal,rp)/dot_v3v3(normal,normal);
1594
sub_v3_v3v3(rp, p, origin);
1595
h = dot_v3v3(normal, rp) / dot_v3v3(normal, normal);
1224
1596
if (h < 0.0f || h > 1.0f) return 0;
1228
/*mama (knowing the squared lenght of the normal)*/
1229
static int point_in_slice_m(float p[3],float origin[3],float normal[3],float lns)
1600
/*mama (knowing the squared length of the normal)*/
1601
static int point_in_slice_m(float p[3], float origin[3], float normal[3], float lns)
1232
sub_v3_v3v3(rp,p,origin);
1233
h=dot_v3v3(normal,rp)/lns;
1604
sub_v3_v3v3(rp, p, origin);
1605
h = dot_v3v3(normal, rp) / lns;
1234
1606
if (h < 0.0f || h > 1.0f) return 0;
1239
int isect_point_tri_prism_v3(float p[3], float v1[3], float v2[3], float v3[3])
1611
int isect_point_tri_prism_v3(const float p[3], const float v1[3], const float v2[3], const float v3[3])
1241
if(!point_in_slice(p,v1,v2,v3)) return 0;
1242
if(!point_in_slice(p,v2,v3,v1)) return 0;
1243
if(!point_in_slice(p,v3,v1,v2)) return 0;
1613
if (!point_in_slice(p, v1, v2, v3)) return 0;
1614
if (!point_in_slice(p, v2, v3, v1)) return 0;
1615
if (!point_in_slice(p, v3, v1, v2)) return 0;
1247
int clip_line_plane(float p1[3], float p2[3], float plane[4])
1619
int clip_line_plane(float p1[3], float p2[3], const float plane[4])
1249
1621
float dp[3], n[3], div, t, pc[3];
1251
1623
copy_v3_v3(n, plane);
1252
1624
sub_v3_v3v3(dp, p2, p1);
1253
div= dot_v3v3(dp, n);
1625
div = dot_v3v3(dp, n);
1255
if(div == 0.0f) /* parallel */
1627
if (div == 0.0f) /* parallel */
1258
t= -(dot_v3v3(p1, n) + plane[3])/div;
1630
t = -(dot_v3v3(p1, n) + plane[3]) / div;
1261
1633
/* behind plane, completely clipped */
1268
1640
/* intersect plane */
1270
1642
madd_v3_v3v3fl(pc, p1, dp, t);
1271
1643
copy_v3_v3(p1, pc);
1668
void plot_line_v2v2i(const int p1[2], const int p2[2], int (*callback)(int, int, void *), void *userData)
1678
// if x1 == x2 or y1 == y2, then it does not matter what we set here
1679
int delta_x = (x2 > x1 ? (ix = 1, x2 - x1) : (ix = -1, x1 - x2)) << 1;
1680
int delta_y = (y2 > y1 ? (iy = 1, y2 - y1) : (iy = -1, y1 - y2)) << 1;
1682
if (callback(x1, y1, userData) == 0) {
1686
if (delta_x >= delta_y) {
1687
// error may go below zero
1688
int error = delta_y - (delta_x >> 1);
1692
if (error || (ix > 0)) {
1703
if (callback(x1, y1, userData) == 0) {
1709
// error may go below zero
1710
int error = delta_x - (delta_y >> 1);
1714
if (error || (iy > 0)) {
1725
if (callback(x1, y1, userData) == 0) {
1296
1732
/****************************** Interpolation ********************************/
1298
static float tri_signed_area(float *v1, float *v2, float *v3, int i, int j)
1300
return 0.5f*((v1[i]-v2[i])*(v2[j]-v3[j]) + (v1[j]-v2[j])*(v3[i]-v2[i]));
1303
static int barycentric_weights(float *v1, float *v2, float *v3, float *co, float *n, float *w)
1305
float xn, yn, zn, a1, a2, a3, asum;
1308
/* find best projection of face XY, XZ or YZ: barycentric weights of
1309
the 2d projected coords are the same and faster to compute */
1310
xn= (float)fabs(n[0]);
1311
yn= (float)fabs(n[1]);
1312
zn= (float)fabs(n[2]);
1313
if(zn>=xn && zn>=yn) {i= 0; j= 1;}
1314
else if(yn>=xn && yn>=zn) {i= 0; j= 2;}
1317
a1= tri_signed_area(v2, v3, co, i, j);
1318
a2= tri_signed_area(v3, v1, co, i, j);
1319
a3= tri_signed_area(v1, v2, co, i, j);
1323
if (fabs(asum) < FLT_EPSILON) {
1734
/* get the 2 dominant axis values, 0==X, 1==Y, 2==Z */
1735
void axis_dominant_v3(int *axis_a, int *axis_b, const float axis[3])
1737
const float xn = fabsf(axis[0]);
1738
const float yn = fabsf(axis[1]);
1739
const float zn = fabsf(axis[2]);
1741
if (zn >= xn && zn >= yn) { *axis_a= 0; *axis_b = 1; }
1742
else if (yn >= xn && yn >= zn) { *axis_a= 0; *axis_b = 2; }
1743
else { *axis_a= 1; *axis_b = 2; }
1746
static float tri_signed_area(const float v1[3], const float v2[3], const float v3[3], const int i, const int j)
1748
return 0.5f * ((v1[i] - v2[i]) * (v2[j] - v3[j]) + (v1[j] - v2[j]) * (v3[i] - v2[i]));
1751
/* return 1 when degenerate */
1752
static int barycentric_weights(const float v1[3], const float v2[3], const float v3[3], const float co[3], const float n[3], float w[3])
1757
axis_dominant_v3(&i, &j, n);
1759
w[0] = tri_signed_area(v2, v3, co, i, j);
1760
w[1] = tri_signed_area(v3, v1, co, i, j);
1761
w[2] = tri_signed_area(v1, v2, co, i, j);
1763
wtot = w[0] + w[1] + w[2];
1765
if (fabsf(wtot) > FLT_EPSILON) {
1766
mul_v3_fl(w, 1.0f / wtot);
1324
1770
/* zero area triangle */
1325
w[0]= w[1]= w[2]= 1.0f/3.0f;
1771
copy_v3_fl(w, 1.0f / 3.0f);
1337
void interp_weights_face_v3(float *w,float *v1, float *v2, float *v3, float *v4, float *co)
1776
void interp_weights_face_v3(float w[4], const float v1[3], const float v2[3], const float v3[3], const float v4[3], const float co[3])
1341
w[0]= w[1]= w[2]= w[3]= 0.0f;
1780
w[0] = w[1] = w[2] = w[3] = 0.0f;
1343
1782
/* first check for exact match */
1344
if(equals_v3v3(co, v1))
1346
else if(equals_v3v3(co, v2))
1348
else if(equals_v3v3(co, v3))
1350
else if(v4 && equals_v3v3(co, v4))
1783
if (equals_v3v3(co, v1))
1785
else if (equals_v3v3(co, v2))
1787
else if (equals_v3v3(co, v3))
1789
else if (v4 && equals_v3v3(co, v4))
1353
1792
/* otherwise compute barycentric interpolation weights */
1354
1793
float n1[3], n2[3], n[3];
1546
1984
sub_v3_v3v3(d3, v3, v1);
1547
1985
cross_v3_v3v3(cross, d2, d3);
1549
area= len_v3(cross);
1550
dot= dot_v3v3(d2, d3);
1551
len= len_v3(d2)*len_v3(d3);
1987
area = len_v3(cross);
1988
dot = dot_v3v3(d2, d3);
1989
len = len_v3(d2) * len_v3(d3);
1556
return (len - dot)/area;
1994
return (len - dot) / area;
1559
void interp_weights_poly_v3(float *w,float v[][3], int n, float *co)
1997
void interp_weights_poly_v3(float *w, float v[][3], const int n, const float co[3])
1561
1999
float totweight, t1, t2, len, *vmid, *vprev, *vnext;
1566
for(i=0; i<n; i++) {
1568
vprev= (i == 0)? v[n-1]: v[i-1];
1569
vnext= (i == n-1)? v[0]: v[i+1];
1571
t1= mean_value_half_tan(co, vprev, vmid);
1572
t2= mean_value_half_tan(co, vmid, vnext);
1574
len= len_v3v3(co, vmid);
2004
for (i = 0; i < n; i++) {
2006
vprev = (i == 0) ? v[n - 1] : v[i - 1];
2007
vnext = (i == n - 1) ? v[0] : v[i + 1];
2009
t1 = mean_value_half_tan(co, vprev, vmid);
2010
t2 = mean_value_half_tan(co, vmid, vnext);
2012
len = len_v3v3(co, vmid);
2013
w[i] = (t1 + t2) / len;
1576
2014
totweight += w[i];
1579
if(totweight != 0.0f)
2017
if (totweight != 0.0f)
2018
for (i = 0; i < n; i++)
1581
2019
w[i] /= totweight;
1584
2022
/* (x1,v1)(t1=0)------(x2,v2)(t2=1), 0<t<1 --> (x,v)(t) */
1585
void interp_cubic_v3(float *x, float *v,float *x1, float *v1, float *x2, float *v2, float t)
2023
void interp_cubic_v3(float x[3], float v[3], const float x1[3], const float v1[3], const float x2[3], const float v2[3], const float t)
1591
2029
/* cubic interpolation */
1592
a[0]= v1[0] + v2[0] + 2*(x1[0] - x2[0]);
1593
a[1]= v1[1] + v2[1] + 2*(x1[1] - x2[1]);
1594
a[2]= v1[2] + v2[2] + 2*(x1[2] - x2[2]);
1596
b[0]= -2*v1[0] - v2[0] - 3*(x1[0] - x2[0]);
1597
b[1]= -2*v1[1] - v2[1] - 3*(x1[1] - x2[1]);
1598
b[2]= -2*v1[2] - v2[2] - 3*(x1[2] - x2[2]);
1600
x[0]= a[0]*t3 + b[0]*t2 + v1[0]*t + x1[0];
1601
x[1]= a[1]*t3 + b[1]*t2 + v1[1]*t + x1[1];
1602
x[2]= a[2]*t3 + b[2]*t2 + v1[2]*t + x1[2];
1604
v[0]= 3*a[0]*t2 + 2*b[0]*t + v1[0];
1605
v[1]= 3*a[1]*t2 + 2*b[1]*t + v1[1];
1606
v[2]= 3*a[2]*t2 + 2*b[2]*t + v1[2];
2030
a[0] = v1[0] + v2[0] + 2 * (x1[0] - x2[0]);
2031
a[1] = v1[1] + v2[1] + 2 * (x1[1] - x2[1]);
2032
a[2] = v1[2] + v2[2] + 2 * (x1[2] - x2[2]);
2034
b[0] = -2 * v1[0] - v2[0] - 3 * (x1[0] - x2[0]);
2035
b[1] = -2 * v1[1] - v2[1] - 3 * (x1[1] - x2[1]);
2036
b[2] = -2 * v1[2] - v2[2] - 3 * (x1[2] - x2[2]);
2038
x[0] = a[0] * t3 + b[0] * t2 + v1[0] * t + x1[0];
2039
x[1] = a[1] * t3 + b[1] * t2 + v1[1] * t + x1[1];
2040
x[2] = a[2] * t3 + b[2] * t2 + v1[2] * t + x1[2];
2042
v[0] = 3 * a[0] * t2 + 2 * b[0] * t + v1[0];
2043
v[1] = 3 * a[1] * t2 + 2 * b[1] * t + v1[1];
2044
v[2] = 3 * a[2] * t2 + 2 * b[2] * t + v1[2];
2047
/* unfortunately internal calculations have to be done at double precision to achieve correct/stable results. */
2049
#define IS_ZERO(x) ((x > (-DBL_EPSILON) && x < DBL_EPSILON) ? 1 : 0)
2051
/* Barycentric reverse */
2052
void resolve_tri_uv(float r_uv[2], const float st[2], const float st0[2], const float st1[2], const float st2[2])
2054
/* find UV such that
2055
* t = u * t0 + v * t1 + (1 - u - v) * t2
2056
* u * (t0 - t2) + v * (t1 - t2) = t - t2 */
2057
const double a = st0[0] - st2[0], b = st1[0] - st2[0];
2058
const double c = st0[1] - st2[1], d = st1[1] - st2[1];
2059
const double det = a * d - c * b;
2061
if (IS_ZERO(det) == 0) { /* det should never be zero since the determinant is the signed ST area of the triangle. */
2062
const double x[] = {st[0] - st2[0], st[1] - st2[1]};
2064
r_uv[0] = (float)((d * x[0] - b * x[1]) / det);
2065
r_uv[1] = (float)(((-c) * x[0] + a * x[1]) / det);
2070
/* bilinear reverse */
2071
void resolve_quad_uv(float r_uv[2], const float st[2], const float st0[2], const float st1[2], const float st2[2], const float st3[2])
2073
const double signed_area = (st0[0] * st1[1] - st0[1] * st1[0]) + (st1[0] * st2[1] - st1[1] * st2[0]) +
2074
(st2[0] * st3[1] - st2[1] * st3[0]) + (st3[0] * st0[1] - st3[1] * st0[0]);
2076
/* X is 2D cross product (determinant)
2077
* A= (p0-p) X (p0-p3)*/
2078
const double a = (st0[0] - st[0]) * (st0[1] - st3[1]) - (st0[1] - st[1]) * (st0[0] - st3[0]);
2080
/* B= ( (p0-p) X (p1-p2) + (p1-p) X (p0-p3) ) / 2 */
2081
const double b = 0.5 * (((st0[0] - st[0]) * (st1[1] - st2[1]) - (st0[1] - st[1]) * (st1[0] - st2[0])) +
2082
((st1[0] - st[0]) * (st0[1] - st3[1]) - (st1[1] - st[1]) * (st0[0] - st3[0])));
2084
/* C = (p1-p) X (p1-p2) */
2085
const double fC = (st1[0] - st[0]) * (st1[1] - st2[1]) - (st1[1] - st[1]) * (st1[0] - st2[0]);
2086
const double denom = a - 2 * b + fC;
2091
if (IS_ZERO(denom) != 0) {
2092
const double fDen = a - fC;
2093
if (IS_ZERO(fDen) == 0)
2094
r_uv[0] = (float)(a / fDen);
2097
const double desc_sq = b * b - a * fC;
2098
const double desc = sqrt(desc_sq < 0.0 ? 0.0 : desc_sq);
2099
const double s = signed_area > 0 ? (-1.0) : 1.0;
2101
r_uv[0] = (float)(((a - b) + s * desc) / denom);
2104
/* find UV such that
2105
* fST = (1-u)(1-v) * ST0 + u * (1-v) * ST1 + u * v * ST2 + (1-u) * v * ST3 */
2107
const double denom_s = (1 - r_uv[0]) * (st0[0] - st3[0]) + r_uv[0] * (st1[0] - st2[0]);
2108
const double denom_t = (1 - r_uv[0]) * (st0[1] - st3[1]) + r_uv[0] * (st1[1] - st2[1]);
2110
double denom = denom_s;
2112
if (fabs(denom_s) < fabs(denom_t)) {
2117
if (IS_ZERO(denom) == 0)
2118
r_uv[1] = (float)(((1.0f - r_uv[0]) * (st0[i] - st[i]) + r_uv[0] * (st1[i] - st[i])) / denom);
1609
2124
/***************************** View & Projection *****************************/
1611
void orthographic_m4(float matrix[][4],float left, float right, float bottom, float top, float nearClip, float farClip)
2126
void orthographic_m4(float matrix[][4], const float left, const float right, const float bottom, const float top, const float nearClip, const float farClip)
1613
2128
float Xdelta, Ydelta, Zdelta;
1615
2130
Xdelta = right - left;
1616
2131
Ydelta = top - bottom;
1617
2132
Zdelta = farClip - nearClip;
1618
if (Xdelta == 0.0 || Ydelta == 0.0 || Zdelta == 0.0) {
2133
if (Xdelta == 0.0f || Ydelta == 0.0f || Zdelta == 0.0f) {
1621
2136
unit_m4(matrix);
1622
matrix[0][0] = 2.0f/Xdelta;
1623
matrix[3][0] = -(right + left)/Xdelta;
1624
matrix[1][1] = 2.0f/Ydelta;
1625
matrix[3][1] = -(top + bottom)/Ydelta;
1626
matrix[2][2] = -2.0f/Zdelta; /* note: negate Z */
1627
matrix[3][2] = -(farClip + nearClip)/Zdelta;
2137
matrix[0][0] = 2.0f / Xdelta;
2138
matrix[3][0] = -(right + left) / Xdelta;
2139
matrix[1][1] = 2.0f / Ydelta;
2140
matrix[3][1] = -(top + bottom) / Ydelta;
2141
matrix[2][2] = -2.0f / Zdelta; /* note: negate Z */
2142
matrix[3][2] = -(farClip + nearClip) / Zdelta;
1630
void perspective_m4(float mat[][4],float left, float right, float bottom, float top, float nearClip, float farClip)
2145
void perspective_m4(float mat[4][4], const float left, const float right, const float bottom, const float top, const float nearClip, const float farClip)
1632
2147
float Xdelta, Ydelta, Zdelta;
1789
2333
/********************************** Mapping **********************************/
1791
void map_to_tube(float *u, float *v,float x, float y, float z)
1795
*v = (z + 1.0f) / 2.0f;
1797
len= (float)sqrt(x*x+y*y);
1799
*u = (float)((1.0 - (atan2(x/len,y/len) / M_PI)) / 2.0);
1801
*v = *u = 0.0f; /* to avoid un-initialized variables */
1804
void map_to_sphere(float *u, float *v,float x, float y, float z)
1808
len= (float)sqrt(x*x+y*y+z*z);
1810
if(x==0.0f && y==0.0f) *u= 0.0f; /* othwise domain error */
1811
else *u = (float)((1.0 - (float)atan2(x,y) / M_PI) / 2.0);
1814
*v = 1.0f - (float)saacos(z)/(float)M_PI;
1816
*v = *u = 0.0f; /* to avoid un-initialized variables */
2335
void map_to_tube(float *r_u, float *r_v, const float x, const float y, const float z)
2339
*r_v = (z + 1.0f) / 2.0f;
2341
len = sqrtf(x * x + y * y);
2343
*r_u = (float)((1.0 - (atan2(x / len, y / len) / M_PI)) / 2.0);
2346
*r_v = *r_u = 0.0f; /* to avoid un-initialized variables */
2350
void map_to_sphere(float *r_u, float *r_v, const float x, const float y, const float z)
2354
len = sqrtf(x * x + y * y + z * z);
2356
if (x == 0.0f && y == 0.0f) *r_u = 0.0f; /* othwise domain error */
2357
else *r_u = (1.0f - atan2f(x, y) / (float)M_PI) / 2.0f;
2359
*r_v = 1.0f - (float)saacos(z / len) / (float)M_PI;
2362
*r_v = *r_u = 0.0f; /* to avoid un-initialized variables */
2366
/********************************* Normals **********************************/
2368
void accumulate_vertex_normals(float n1[3], float n2[3], float n3[3],
2369
float n4[3], const float f_no[3], const float co1[3], const float co2[3],
2370
const float co3[3], const float co4[3])
2373
const int nverts = (n4 != NULL && co4 != NULL) ? 4 : 3;
2375
/* compute normalized edge vectors */
2376
sub_v3_v3v3(vdiffs[0], co2, co1);
2377
sub_v3_v3v3(vdiffs[1], co3, co2);
2380
sub_v3_v3v3(vdiffs[2], co1, co3);
2383
sub_v3_v3v3(vdiffs[2], co4, co3);
2384
sub_v3_v3v3(vdiffs[3], co1, co4);
2385
normalize_v3(vdiffs[3]);
2388
normalize_v3(vdiffs[0]);
2389
normalize_v3(vdiffs[1]);
2390
normalize_v3(vdiffs[2]);
2392
/* accumulate angle weighted face normal */
2394
float *vn[] = {n1, n2, n3, n4};
2395
const float *prev_edge = vdiffs[nverts - 1];
2398
for (i = 0; i < nverts; i++) {
2399
const float *cur_edge = vdiffs[i];
2400
const float fac = saacos(-dot_v3v3(cur_edge, prev_edge));
2403
madd_v3_v3fl(vn[i], f_no, fac);
2404
prev_edge = cur_edge;
2409
/* Add weighted face normal component into normals of the face vertices.
2410
* Caller must pass pre-allocated vdiffs of nverts length. */
2411
void accumulate_vertex_normals_poly(float **vertnos, float polyno[3],
2412
float **vertcos, float vdiffs[][3], int nverts)
2416
/* calculate normalized edge directions for each edge in the poly */
2417
for (i = 0; i < nverts; i++) {
2418
sub_v3_v3v3(vdiffs[i], vertcos[(i + 1) % nverts], vertcos[i]);
2419
normalize_v3(vdiffs[i]);
2422
/* accumulate angle weighted face normal */
2424
const float *prev_edge = vdiffs[nverts - 1];
2427
for (i = 0; i < nverts; i++) {
2428
const float *cur_edge = vdiffs[i];
2430
/* calculate angle between the two poly edges incident on
2432
const float fac = saacos(-dot_v3v3(cur_edge, prev_edge));
2435
madd_v3_v3fl(vertnos[i], polyno, fac);
2436
prev_edge = cur_edge;
1820
2441
/********************************* Tangents **********************************/
1822
2443
/* For normal map tangents we need to detect uv boundaries, and only average
1823
* tangents in case the uvs are connected. Alternative would be to store 1
2444
* tangents in case the uvs are connected. Alternative would be to store 1
1824
2445
* tangent per face rather than 4 per face vertex, but that's not compatible
1825
2446
* with games */
1828
2449
/* from BKE_mesh.h */
1829
#define STD_UV_CONNECT_LIMIT 0.0001f
2450
#define STD_UV_CONNECT_LIMIT 0.0001f
1831
void sum_or_add_vertex_tangent(void *arena, VertexTangent **vtang, float *tang, float *uv)
2452
void sum_or_add_vertex_tangent(void *arena, VertexTangent **vtang, const float tang[3], const float uv[2])
1833
2454
VertexTangent *vt;
1835
2456
/* find a tangent with connected uvs */
1836
for(vt= *vtang; vt; vt=vt->next) {
1837
if(fabs(uv[0]-vt->uv[0]) < STD_UV_CONNECT_LIMIT && fabs(uv[1]-vt->uv[1]) < STD_UV_CONNECT_LIMIT) {
2457
for (vt = *vtang; vt; vt = vt->next) {
2458
if (fabsf(uv[0] - vt->uv[0]) < STD_UV_CONNECT_LIMIT && fabsf(uv[1] - vt->uv[1]) < STD_UV_CONNECT_LIMIT) {
1838
2459
add_v3_v3(vt->tang, tang);
1843
2464
/* if not found, append a new one */
1844
vt= BLI_memarena_alloc((MemArena *)arena, sizeof(VertexTangent));
2465
vt = BLI_memarena_alloc((MemArena *) arena, sizeof(VertexTangent));
1845
2466
copy_v3_v3(vt->tang, tang);
1854
float *find_vertex_tangent(VertexTangent *vtang, float *uv)
2475
float *find_vertex_tangent(VertexTangent *vtang, const float uv[2])
1856
2477
VertexTangent *vt;
1857
2478
static float nulltang[3] = {0.0f, 0.0f, 0.0f};
1859
for(vt= vtang; vt; vt=vt->next)
1860
if(fabs(uv[0]-vt->uv[0]) < STD_UV_CONNECT_LIMIT && fabs(uv[1]-vt->uv[1]) < STD_UV_CONNECT_LIMIT)
2480
for (vt = vtang; vt; vt = vt->next)
2481
if (fabsf(uv[0] - vt->uv[0]) < STD_UV_CONNECT_LIMIT && fabsf(uv[1] - vt->uv[1]) < STD_UV_CONNECT_LIMIT)
1861
2482
return vt->tang;
1863
return nulltang; /* shouldn't happen, except for nan or so */
2484
return nulltang; /* shouldn't happen, except for nan or so */
1866
void tangent_from_uv(float *uv1, float *uv2, float *uv3, float *co1, float *co2, float *co3, float *n, float *tang)
2487
void tangent_from_uv(float uv1[2], float uv2[2], float uv3[3], float co1[3], float co2[3], float co3[3], float n[3], float tang[3])
1868
float s1= uv2[0] - uv1[0];
1869
float s2= uv3[0] - uv1[0];
1870
float t1= uv2[1] - uv1[1];
1871
float t2= uv3[1] - uv1[1];
1872
float det= (s1 * t2 - s2 * t1);
2489
float s1 = uv2[0] - uv1[0];
2490
float s2 = uv3[0] - uv1[0];
2491
float t1 = uv2[1] - uv1[1];
2492
float t2 = uv3[1] - uv1[1];
2493
float det = (s1 * t2 - s2 * t1);
1874
if(det != 0.0f) { /* otherwise 'tang' becomes nan */
2495
if (det != 0.0f) { /* otherwise 'tang' becomes nan */
1875
2496
float tangv[3], ct[3], e1[3], e2[3];
1879
2500
/* normals in render are inversed... */
1880
2501
sub_v3_v3v3(e1, co1, co2);
1881
2502
sub_v3_v3v3(e2, co1, co3);
1882
tang[0] = (t2*e1[0] - t1*e2[0])*det;
1883
tang[1] = (t2*e1[1] - t1*e2[1])*det;
1884
tang[2] = (t2*e1[2] - t1*e2[2])*det;
1885
tangv[0] = (s1*e2[0] - s2*e1[0])*det;
1886
tangv[1] = (s1*e2[1] - s2*e1[1])*det;
1887
tangv[2] = (s1*e2[2] - s2*e1[2])*det;
2503
tang[0] = (t2 * e1[0] - t1 * e2[0]) * det;
2504
tang[1] = (t2 * e1[1] - t1 * e2[1]) * det;
2505
tang[2] = (t2 * e1[2] - t1 * e2[2]) * det;
2506
tangv[0] = (s1 * e2[0] - s2 * e1[0]) * det;
2507
tangv[1] = (s1 * e2[1] - s2 * e1[1]) * det;
2508
tangv[2] = (s1 * e2[2] - s2 * e1[2]) * det;
1888
2509
cross_v3_v3v3(ct, tang, tangv);
1890
2511
/* check flip */
1891
if ((ct[0]*n[0] + ct[1]*n[1] + ct[2]*n[2]) < 0.0f)
2512
if (dot_v3v3(ct, n) < 0.0f) {
1892
2513
negate_v3(tang);
1895
tang[0]= tang[1]= tang[2]= 0.0;
2517
tang[0] = tang[1] = tang[2] = 0.0;
added
removed
source/blender/blenlib/intern/math_geom.c
