diff options
Diffstat (limited to 'source/blender/blenlib/intern/math_geom.c')
-rw-r--r-- | source/blender/blenlib/intern/math_geom.c | 263 |
1 files changed, 212 insertions, 51 deletions
diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c index 68b1feea632..d7a71f8567c 100644 --- a/source/blender/blenlib/intern/math_geom.c +++ b/source/blender/blenlib/intern/math_geom.c @@ -30,8 +30,9 @@ #include "BLI_math.h" #include "BLI_memarena.h" +#include "BLI_utildefines.h" + -#include "BKE_utildefines.h" /********************************** Polygons *********************************/ @@ -233,53 +234,114 @@ float dist_to_line_segment_v3(float *v1, float *v2, float *v3) /******************************* Intersection ********************************/ /* intersect Line-Line, shorts */ -int isect_line_line_v2_short(short *v1, short *v2, short *v3, short *v4) +int isect_line_line_v2_short(const short *v1, const short *v2, const short *v3, const short *v4) { - /* return: - -1: colliniar - 0: no intersection of segments - 1: exact intersection of segments - 2: cross-intersection of segments - */ float div, labda, mu; div= (float)((v2[0]-v1[0])*(v4[1]-v3[1])-(v2[1]-v1[1])*(v4[0]-v3[0])); - if(div==0.0f) return -1; + if(div==0.0f) return ISECT_LINE_LINE_COLINEAR; labda= ((float)(v1[1]-v3[1])*(v4[0]-v3[0])-(v1[0]-v3[0])*(v4[1]-v3[1]))/div; mu= ((float)(v1[1]-v3[1])*(v2[0]-v1[0])-(v1[0]-v3[0])*(v2[1]-v1[1]))/div; if(labda>=0.0f && labda<=1.0f && mu>=0.0f && mu<=1.0f) { - if(labda==0.0f || labda==1.0f || mu==0.0f || mu==1.0f) return 1; - return 2; + if(labda==0.0f || labda==1.0f || mu==0.0f || mu==1.0f) return ISECT_LINE_LINE_EXACT; + return ISECT_LINE_LINE_CROSS; } - return 0; + return ISECT_LINE_LINE_NONE; } /* intersect Line-Line, floats */ -int isect_line_line_v2(float *v1, float *v2, float *v3, float *v4) +int isect_line_line_v2(const float *v1, const float *v2, const float *v3, const float *v4) { - /* return: - -1: colliniar -0: no intersection of segments -1: exact intersection of segments -2: cross-intersection of segments - */ float div, labda, mu; div= (v2[0]-v1[0])*(v4[1]-v3[1])-(v2[1]-v1[1])*(v4[0]-v3[0]); - if(div==0.0) return -1; + if(div==0.0) return ISECT_LINE_LINE_COLINEAR; labda= ((float)(v1[1]-v3[1])*(v4[0]-v3[0])-(v1[0]-v3[0])*(v4[1]-v3[1]))/div; mu= ((float)(v1[1]-v3[1])*(v2[0]-v1[0])-(v1[0]-v3[0])*(v2[1]-v1[1]))/div; if(labda>=0.0 && labda<=1.0 && mu>=0.0 && mu<=1.0) { - if(labda==0.0 || labda==1.0 || mu==0.0 || mu==1.0) return 1; - return 2; + if(labda==0.0 || labda==1.0 || mu==0.0 || mu==1.0) return ISECT_LINE_LINE_EXACT; + return ISECT_LINE_LINE_CROSS; } - return 0; + return ISECT_LINE_LINE_NONE; +} + +/* get intersection point of two 2D segments and return intersection type: + -1: colliniar + 1: intersection */ +int isect_seg_seg_v2_point(const float *v1, const float *v2, const float *v3, const float *v4, float vi[2]) +{ + float a1, a2, b1, b2, c1, c2, d; + float u, v; + const float eps= 0.000001f; + + a1= v2[0]-v1[0]; + b1= v4[0]-v3[0]; + c1= v1[0]-v4[0]; + + a2= v2[1]-v1[1]; + b2= v4[1]-v3[1]; + c2= v1[1]-v4[1]; + + d= a1*b2-a2*b1; + + if(d==0) { + if(a1*c2-a2*c1==0.0f && b1*c2-b2*c1==0.0f) { /* equal lines */ + float a[2], b[2], c[2]; + float u2; + + if(len_v2v2(v1, v2)==0.0f) { + if(len_v2v2(v3, v4)>eps) { + /* use non-point segment as basis */ + SWAP(const float *, v1, v3); + SWAP(const float *, v2, v4); + } else { /* both of segments are points */ + if(equals_v2v2(v1, v3)) { /* points are equal */ + copy_v2_v2(vi, v1); + return 1; + } + + /* two different points */ + return -1; + } + } + + sub_v2_v2v2(a, v3, v1); + sub_v2_v2v2(b, v2, v1); + sub_v2_v2v2(c, v2, v1); + u= dot_v2v2(a, b) / dot_v2v2(c, c); + + sub_v2_v2v2(a, v4, v1); + u2= dot_v2v2(a, b) / dot_v2v2(c, c); + + if(u>u2) SWAP(float, u, u2); + + if(u>1.0f+eps || u2<-eps) return -1; /* non-ovlerlapping segments */ + else if(maxf(0.0f, u) == minf(1.0f, u2)){ /* one common point: can return result */ + interp_v2_v2v2(vi, v1, v2, maxf(0, u)); + return 1; + } + } + + /* lines are colliniar */ + return -1; + } + + u= (c2*b1-b2*c1)/d; + v= (c1*a2-a1*c2)/d; + + if(u>=-eps && u<=1.0f+eps && v>=-eps && v<=1.0f+eps) { /* intersection */ + interp_v2_v2v2(vi, v1, v2, u); + return 1; + } + + /* out of segment intersection */ + return -1; } /* @@ -336,20 +398,19 @@ static short IsectLLPt2Df(float x0,float y0,float x1,float y1, return 1; } // end Intersect_Lines -#define SIDE_OF_LINE(pa,pb,pp) ((pa[0]-pp[0])*(pb[1]-pp[1]))-((pb[0]-pp[0])*(pa[1]-pp[1])) /* point in tri */ -// XXX was called IsectPT2Df + int isect_point_tri_v2(float pt[2], float v1[2], float v2[2], float v3[2]) { - if (SIDE_OF_LINE(v1,v2,pt)>=0.0) { - if (SIDE_OF_LINE(v2,v3,pt)>=0.0) { - if (SIDE_OF_LINE(v3,v1,pt)>=0.0) { + if (line_point_side_v2(v1,v2,pt)>=0.0) { + if (line_point_side_v2(v2,v3,pt)>=0.0) { + if (line_point_side_v2(v3,v1,pt)>=0.0) { return 1; } } } else { - if (! (SIDE_OF_LINE(v2,v3,pt)>=0.0)) { - if (! (SIDE_OF_LINE(v3,v1,pt)>=0.0)) { + if (! (line_point_side_v2(v2,v3,pt)>=0.0)) { + if (! (line_point_side_v2(v3,v1,pt)>=0.0)) { return -1; } } @@ -360,18 +421,18 @@ int isect_point_tri_v2(float pt[2], float v1[2], float v2[2], float v3[2]) /* point in quad - only convex quads */ int isect_point_quad_v2(float pt[2], float v1[2], float v2[2], float v3[2], float v4[2]) { - if (SIDE_OF_LINE(v1,v2,pt)>=0.0) { - if (SIDE_OF_LINE(v2,v3,pt)>=0.0) { - if (SIDE_OF_LINE(v3,v4,pt)>=0.0) { - if (SIDE_OF_LINE(v4,v1,pt)>=0.0) { + if (line_point_side_v2(v1,v2,pt)>=0.0) { + if (line_point_side_v2(v2,v3,pt)>=0.0) { + if (line_point_side_v2(v3,v4,pt)>=0.0) { + if (line_point_side_v2(v4,v1,pt)>=0.0) { return 1; } } } } else { - if (! (SIDE_OF_LINE(v2,v3,pt)>=0.0)) { - if (! (SIDE_OF_LINE(v3,v4,pt)>=0.0)) { - if (! (SIDE_OF_LINE(v4,v1,pt)>=0.0)) { + if (! (line_point_side_v2(v2,v3,pt)>=0.0)) { + if (! (line_point_side_v2(v3,v4,pt)>=0.0)) { + if (! (line_point_side_v2(v4,v1,pt)>=0.0)) { return -1; } } @@ -557,7 +618,7 @@ static int getLowestRoot(float a, float b, float c, float maxR, float* root) if (determinant >= 0.0f) { // calculate the two roots: (if determinant == 0 then - // x1==x2 but let’s disregard that slight optimization) + // x1==x2 but lets disregard that slight optimization) float sqrtD = (float)sqrt(determinant); float r1 = (-b - sqrtD) / (2.0f*a); float r2 = (-b + sqrtD) / (2.0f*a); @@ -589,7 +650,7 @@ int isect_sweeping_sphere_tri_v3(float p1[3], float p2[3], float radius, float v { float e1[3], e2[3], e3[3], point[3], vel[3], /*dist[3],*/ nor[3], temp[3], bv[3]; float a, b, c, d, e, x, y, z, radius2=radius*radius; - float elen2,edotv,edotbv,nordotv,vel2; + float elen2,edotv,edotbv,nordotv; float newLambda; int found_by_sweep=0; @@ -663,7 +724,7 @@ int isect_sweeping_sphere_tri_v3(float p1[3], float p2[3], float radius, float v *lambda=1.0f; /*---test points---*/ - a=vel2=dot_v3v3(vel,vel); + a=dot_v3v3(vel,vel); /*v0*/ sub_v3_v3v3(temp,p1,v0); @@ -752,10 +813,10 @@ int isect_sweeping_sphere_tri_v3(float p1[3], float p2[3], float radius, float v } /*e3*/ - sub_v3_v3v3(bv,v0,p1); - elen2 = dot_v3v3(e1,e1); - edotv = dot_v3v3(e1,vel); - edotbv = dot_v3v3(e1,bv); + /* sub_v3_v3v3(bv,v0,p1); */ /* UNUSED */ + /* elen2 = dot_v3v3(e1,e1); */ /* UNUSED */ + /* edotv = dot_v3v3(e1,vel); */ /* UNUSED */ + /* edotbv = dot_v3v3(e1,bv); */ /* UNUSED */ sub_v3_v3v3(bv,v1,p1); elen2 = dot_v3v3(e3,e3); @@ -900,7 +961,6 @@ int isect_line_line_strict_v3(float v1[3], float v2[3], float v3[3], float v4[3] { float a[3], b[3], c[3], ab[3], cb[3], ca[3], dir1[3], dir2[3]; float d; - float d1; sub_v3_v3v3(c, v3, v1); sub_v3_v3v3(a, v2, v1); @@ -913,8 +973,6 @@ int isect_line_line_strict_v3(float v1[3], float v2[3], float v3[3], float v4[3] /* colinear or one vector is zero-length*/ return 0; } - - d1 = d; cross_v3_v3v3(ab, a, b); d = dot_v3v3(c, ab); @@ -961,7 +1019,7 @@ int isect_aabb_aabb_v3(float min1[3], float max1[3], float min2[3], float max2[3 /* find closest point to p on line through l1,l2 and return lambda, * where (0 <= lambda <= 1) when cp is in the line segement l1,l2 */ -float closest_to_line_v3(float cp[3],float p[3], float l1[3], float l2[3]) +float closest_to_line_v3(float cp[3], const float p[3], const float l1[3], const float l2[3]) { float h[3],u[3],lambda; sub_v3_v3v3(u, l2, l1); @@ -973,6 +1031,17 @@ float closest_to_line_v3(float cp[3],float p[3], float l1[3], float l2[3]) return lambda; } +float closest_to_line_v2(float cp[2],const float p[2], const float l1[2], const float l2[2]) +{ + float h[2],u[2],lambda; + sub_v2_v2v2(u, l2, l1); + sub_v2_v2v2(h, p, l1); + lambda =dot_v2v2(u,h)/dot_v2v2(u,u); + cp[0] = l1[0] + u[0] * lambda; + cp[1] = l1[1] + u[1] * lambda; + return lambda; +} + #if 0 /* little sister we only need to know lambda */ static float lambda_cp_line(float p[3], float l1[3], float l2[3]) @@ -1289,6 +1358,71 @@ int clip_line_plane(float p1[3], float p2[3], float plane[4]) } } + +void plot_line_v2v2i(int p1[2], int p2[2], int (*callback)(int, int, void *), void *userData) +{ + int x1= p1[0]; + int y1= p1[1]; + int x2= p2[0]; + int y2= p2[1]; + + signed char ix; + signed char iy; + + // if x1 == x2 or y1 == y2, then it does not matter what we set here + int delta_x = (x2 > x1?(ix = 1, x2 - x1):(ix = -1, x1 - x2)) << 1; + int delta_y = (y2 > y1?(iy = 1, y2 - y1):(iy = -1, y1 - y2)) << 1; + + if(callback(x1, y1, userData) == 0) { + return; + } + + if (delta_x >= delta_y) { + // error may go below zero + int error = delta_y - (delta_x >> 1); + + while (x1 != x2) { + if (error >= 0) { + if (error || (ix > 0)) { + y1 += iy; + error -= delta_x; + } + // else do nothing + } + // else do nothing + + x1 += ix; + error += delta_y; + + if(callback(x1, y1, userData) == 0) { + return ; + } + } + } + else { + // error may go below zero + int error = delta_x - (delta_y >> 1); + + while (y1 != y2) { + if (error >= 0) { + if (error || (iy > 0)) { + x1 += ix; + error -= delta_y; + } + // else do nothing + } + // else do nothing + + y1 += iy; + error += delta_x; + + if(callback(x1, y1, userData) == 0) { + return; + } + } + } +} + /****************************** Interpolation ********************************/ static float tri_signed_area(float *v1, float *v2, float *v3, int i, int j) @@ -1647,6 +1781,35 @@ void perspective_m4(float mat[][4],float left, float right, float bottom, float } +/* translate a matrix created by orthographic_m4 or perspective_m4 in XY coords (used to jitter the view) */ +void window_translate_m4(float winmat[][4], float perspmat[][4], float x, float y) +{ + if(winmat[2][3] == -1.0f) { + /* in the case of a win-matrix, this means perspective always */ + float v1[3]; + float v2[3]; + float len1, len2; + + v1[0]= perspmat[0][0]; + v1[1]= perspmat[1][0]; + v1[2]= perspmat[2][0]; + + v2[0]= perspmat[0][1]; + v2[1]= perspmat[1][1]; + v2[2]= perspmat[2][1]; + + len1= (1.0f / len_v3(v1)); + len2= (1.0f / len_v3(v2)); + + winmat[2][0] += len1 * winmat[0][0] * x; + winmat[2][1] += len2 * winmat[1][1] * y; + } + else { + winmat[3][0] += x; + winmat[3][1] += y; + } +} + static void i_multmatrix(float icand[][4], float Vm[][4]) { int row, col; @@ -1676,10 +1839,9 @@ void polarview_m4(float Vm[][4],float dist, float azimuth, float incidence, floa void lookat_m4(float mat[][4],float vx, float vy, float vz, float px, float py, float pz, float twist) { float sine, cosine, hyp, hyp1, dx, dy, dz; - float mat1[4][4]; + float mat1[4][4]= MAT4_UNITY; unit_m4(mat); - unit_m4(mat1); rotate_m4(mat, 'Z', -twist); @@ -2469,4 +2631,3 @@ float form_factor_hemi_poly(float p[3], float n[3], float v1[3], float v2[3], fl return contrib; } - |