diff options
author | Brecht Van Lommel <brechtvanlommel@pandora.be> | 2009-11-10 01:42:41 +0300 |
---|---|---|
committer | Brecht Van Lommel <brechtvanlommel@pandora.be> | 2009-11-10 01:42:41 +0300 |
commit | 60ea7456137a019c2ddad75f14b3b6a0892f7b56 (patch) | |
tree | 7bcd9059496194f5af8d1c671998996a15e38f15 /source/blender/blenlib/intern/math_geom.c | |
parent | 5935ef004935b27fc5795349aed32f87cf637049 (diff) |
Math Lib Reorganization
* New header and source files.
* Still need a few tweaks before switching code to use them.
Diffstat (limited to 'source/blender/blenlib/intern/math_geom.c')
-rw-r--r-- | source/blender/blenlib/intern/math_geom.c | 1598 |
1 files changed, 1598 insertions, 0 deletions
diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c new file mode 100644 index 00000000000..5c7890ffee4 --- /dev/null +++ b/source/blender/blenlib/intern/math_geom.c @@ -0,0 +1,1598 @@ +/** + * $Id$ + * + * ***** BEGIN GPL LICENSE BLOCK ***** + * + * This program is free software; you can redistribute it and/or + * modify it under the terms of the GNU General Public License + * as published by the Free Software Foundation; either version 2 + * of the License, or (at your option) any later version. + * + * This program is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software Foundation, + * Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. + * + * The Original Code is Copyright (C) 2001-2002 by NaN Holding BV. + * All rights reserved. + + * The Original Code is: some of this file. + * + * ***** END GPL LICENSE BLOCK ***** + * */ + +#include <float.h> +#include <stdio.h> +#include <stdlib.h> +#include <string.h> + +#include "BLI_math.h" +#include "BLI_memarena.h" + +/********************************** Polygons *********************************/ + +void cent_tri_v3(float *cent, float *v1, float *v2, float *v3) +{ + cent[0]= 0.33333f*(v1[0]+v2[0]+v3[0]); + cent[1]= 0.33333f*(v1[1]+v2[1]+v3[1]); + cent[2]= 0.33333f*(v1[2]+v2[2]+v3[2]); +} + +void cent_quad_v3(float *cent, float *v1, float *v2, float *v3, float *v4) +{ + cent[0]= 0.25f*(v1[0]+v2[0]+v3[0]+v4[0]); + cent[1]= 0.25f*(v1[1]+v2[1]+v3[1]+v4[1]); + cent[2]= 0.25f*(v1[2]+v2[2]+v3[2]+v4[2]); +} + +float normal_tri_v3(float *n, float *v1, float *v2, float *v3) +{ + float n1[3],n2[3]; + + n1[0]= v1[0]-v2[0]; + n2[0]= v2[0]-v3[0]; + n1[1]= v1[1]-v2[1]; + n2[1]= v2[1]-v3[1]; + n1[2]= v1[2]-v2[2]; + n2[2]= v2[2]-v3[2]; + n[0]= n1[1]*n2[2]-n1[2]*n2[1]; + n[1]= n1[2]*n2[0]-n1[0]*n2[2]; + n[2]= n1[0]*n2[1]-n1[1]*n2[0]; + return normalize_v3(n); +} + +float normal_quad_v3(float *n, float *v1, float *v2, float *v3, float *v4) +{ + /* real cross! */ + float n1[3],n2[3]; + + n1[0]= v1[0]-v3[0]; + n1[1]= v1[1]-v3[1]; + n1[2]= v1[2]-v3[2]; + + n2[0]= v2[0]-v4[0]; + n2[1]= v2[1]-v4[1]; + n2[2]= v2[2]-v4[2]; + + n[0]= n1[1]*n2[2]-n1[2]*n2[1]; + n[1]= n1[2]*n2[0]-n1[0]*n2[2]; + n[2]= n1[0]*n2[1]-n1[1]*n2[0]; + + return normalize_v3(n); +} + +float area_tri_v2(float *v1, float *v2, float *v3) +{ + return (float)(0.5*fabs((v1[0]-v2[0])*(v2[1]-v3[1]) + (v1[1]-v2[1])*(v3[0]-v2[0]))); +} + + +float area_quad_v3(float *v1, float *v2, float *v3, float *v4) /* only convex Quadrilaterals */ +{ + float len, vec1[3], vec2[3], n[3]; + + sub_v3_v3v3(vec1, v2, v1); + sub_v3_v3v3(vec2, v4, v1); + cross_v3_v3v3(n, vec1, vec2); + len= normalize_v3(n); + + sub_v3_v3v3(vec1, v4, v3); + sub_v3_v3v3(vec2, v2, v3); + cross_v3_v3v3(n, vec1, vec2); + len+= normalize_v3(n); + + return (len/2.0f); +} + +float area_tri_v3(float *v1, float *v2, float *v3) /* Triangles */ +{ + float len, vec1[3], vec2[3], n[3]; + + sub_v3_v3v3(vec1, v3, v2); + sub_v3_v3v3(vec2, v1, v2); + cross_v3_v3v3(n, vec1, vec2); + len= normalize_v3(n); + + return (len/2.0f); +} + +#define MAX2(x,y) ((x)>(y) ? (x) : (y)) +#define MAX3(x,y,z) MAX2(MAX2((x),(y)) , (z)) + + +float area_poly_v3(int nr, float *verts, float *normal) +{ + float x, y, z, area, max; + float *cur, *prev; + int a, px=0, py=1; + + /* first: find dominant axis: 0==X, 1==Y, 2==Z */ + x= (float)fabs(normal[0]); + y= (float)fabs(normal[1]); + z= (float)fabs(normal[2]); + max = MAX3(x, y, z); + if(max==y) py=2; + else if(max==x) { + px=1; + py= 2; + } + + /* The Trapezium Area Rule */ + prev= verts+3*(nr-1); + cur= verts; + area= 0; + for(a=0; a<nr; a++) { + area+= (cur[px]-prev[px])*(cur[py]+prev[py]); + prev= cur; + cur+=3; + } + + return (float)fabs(0.5*area/max); +} + +/********************************* Distance **********************************/ + +/* distance v1 to line v2-v3 */ +/* using Hesse formula, NO LINE PIECE! */ +float dist_to_line_v2(float *v1, float *v2, float *v3) +{ + float a[2],deler; + + a[0]= v2[1]-v3[1]; + a[1]= v3[0]-v2[0]; + deler= (float)sqrt(a[0]*a[0]+a[1]*a[1]); + if(deler== 0.0f) return 0; + + return (float)(fabs((v1[0]-v2[0])*a[0]+(v1[1]-v2[1])*a[1])/deler); + +} + +/* distance v1 to line-piece v2-v3 */ +float dist_to_line_segment_v2(float *v1, float *v2, float *v3) +{ + float labda, rc[2], pt[2], len; + + rc[0]= v3[0]-v2[0]; + rc[1]= v3[1]-v2[1]; + len= rc[0]*rc[0]+ rc[1]*rc[1]; + if(len==0.0) { + rc[0]= v1[0]-v2[0]; + rc[1]= v1[1]-v2[1]; + return (float)(sqrt(rc[0]*rc[0]+ rc[1]*rc[1])); + } + + labda= (rc[0]*(v1[0]-v2[0]) + rc[1]*(v1[1]-v2[1]))/len; + if(labda<=0.0) { + pt[0]= v2[0]; + pt[1]= v2[1]; + } + else if(labda>=1.0) { + pt[0]= v3[0]; + pt[1]= v3[1]; + } + else { + pt[0]= labda*rc[0]+v2[0]; + pt[1]= labda*rc[1]+v2[1]; + } + + rc[0]= pt[0]-v1[0]; + rc[1]= pt[1]-v1[1]; + return (float)sqrt(rc[0]*rc[0]+ rc[1]*rc[1]); +} + +/* point closest to v1 on line v2-v3 in 3D */ +void closest_to_line_segment_v3(float *closest, float v1[3], float v2[3], float v3[3]) +{ + float lambda, cp[3]; + + lambda= closest_to_line_v3(cp,v1, v2, v3); + + if(lambda <= 0.0f) + copy_v3_v3(closest, v2); + else if(lambda >= 1.0f) + copy_v3_v3(closest, v3); + else + copy_v3_v3(closest, cp); +} + +/* distance v1 to line-piece v2-v3 in 3D */ +float dist_to_line_segment_v3(float *v1, float *v2, float *v3) +{ + float closest[3]; + + closest_to_line_segment_v3(closest, v1, v2, v3); + + return len_v3v3(closest, v1); +} + +/******************************* Intersection ********************************/ + +/* intersect Line-Line, shorts */ +short isect_line_line_v2_short(short *v1, short *v2, short *v3, 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; + + 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; + } + return 0; +} + +/* intersect Line-Line, floats */ +short isect_line_line_v2(float *v1, float *v2, float *v3, 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; + + 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; + } + return 0; +} + +/* +-1: colliniar + 1: intersection + +*/ +static short IsectLLPt2Df(float x0,float y0,float x1,float y1, + float x2,float y2,float x3,float y3, float *xi,float *yi) + +{ + /* + this function computes the intersection of the sent lines + and returns the intersection point, note that the function assumes + the lines intersect. the function can handle vertical as well + as horizontal lines. note the function isn't very clever, it simply + applies the math, but we don't need speed since this is a + pre-processing step + */ + float c1,c2, // constants of linear equations + det_inv, // the inverse of the determinant of the coefficient + m1,m2; // the slopes of each line + /* + compute slopes, note the cludge for infinity, however, this will + be close enough + */ + if (fabs(x1-x0) > 0.000001) + m1 = (y1-y0) / (x1-x0); + else + return -1; /*m1 = (float) 1e+10;*/ // close enough to infinity + + if (fabs(x3-x2) > 0.000001) + m2 = (y3-y2) / (x3-x2); + else + return -1; /*m2 = (float) 1e+10;*/ // close enough to infinity + + if (fabs(m1-m2) < 0.000001) + return -1; /* paralelle lines */ + +// compute constants + + c1 = (y0-m1*x0); + c2 = (y2-m2*x2); + +// compute the inverse of the determinate + + det_inv = 1.0f / (-m1 + m2); + +// use Kramers rule to compute xi and yi + + *xi= ((-c2 + c1) *det_inv); + *yi= ((m2*c1 - m1*c2) *det_inv); + + 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 */ +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) { + return 1; + } + } + } else { + if (! (SIDE_OF_LINE(v2,v3,pt)>=0.0)) { + if (! (SIDE_OF_LINE(v3,v1,pt)>=0.0)) { + return -1; + } + } + } + + return 0; +} +/* 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) { + 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)) { + return -1; + } + } + } + } + + return 0; +} + +/* moved from effect.c + test if the line starting at p1 ending at p2 intersects the triangle v0..v2 + return non zero if it does +*/ +int isect_line_tri_v3(float p1[3], float p2[3], float v0[3], float v1[3], float v2[3], float *lambda, float *uv) +{ + + float p[3], s[3], d[3], e1[3], e2[3], q[3]; + float a, f, u, v; + + sub_v3_v3v3(e1, v1, v0); + sub_v3_v3v3(e2, v2, v0); + sub_v3_v3v3(d, p2, p1); + + cross_v3_v3v3(p, d, e2); + a = dot_v3v3(e1, p); + if ((a > -0.000001) && (a < 0.000001)) return 0; + f = 1.0f/a; + + sub_v3_v3v3(s, p1, v0); + + cross_v3_v3v3(q, s, e1); + *lambda = f * dot_v3v3(e2, q); + if ((*lambda < 0.0)||(*lambda > 1.0)) return 0; + + u = f * dot_v3v3(s, p); + if ((u < 0.0)||(u > 1.0)) return 0; + + v = f * dot_v3v3(d, q); + if ((v < 0.0)||((u + v) > 1.0)) return 0; + + if(uv) { + uv[0]= u; + uv[1]= v; + } + + return 1; +} + +/* moved from effect.c + test if the ray starting at p1 going in d direction intersects the triangle v0..v2 + return non zero if it does +*/ +int isect_ray_tri_v3(float p1[3], float d[3], float v0[3], float v1[3], float v2[3], float *lambda, float *uv) +{ + float p[3], s[3], e1[3], e2[3], q[3]; + float a, f, u, v; + + sub_v3_v3v3(e1, v1, v0); + sub_v3_v3v3(e2, v2, v0); + + cross_v3_v3v3(p, d, e2); + a = dot_v3v3(e1, p); + if ((a > -0.000001) && (a < 0.000001)) return 0; + f = 1.0f/a; + + sub_v3_v3v3(s, p1, v0); + + cross_v3_v3v3(q, s, e1); + *lambda = f * dot_v3v3(e2, q); + if ((*lambda < 0.0)) return 0; + + u = f * dot_v3v3(s, p); + if ((u < 0.0)||(u > 1.0)) return 0; + + v = f * dot_v3v3(d, q); + if ((v < 0.0)||((u + v) > 1.0)) return 0; + + if(uv) { + uv[0]= u; + uv[1]= v; + } + + return 1; +} + +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) +{ + float p[3], s[3], e1[3], e2[3], q[3]; + float a, f, u, v; + float du = 0, dv = 0; + + sub_v3_v3v3(e1, v1, v0); + sub_v3_v3v3(e2, v2, v0); + + cross_v3_v3v3(p, d, e2); + a = dot_v3v3(e1, p); + if ((a > -0.000001) && (a < 0.000001)) return 0; + f = 1.0f/a; + + sub_v3_v3v3(s, p1, v0); + + cross_v3_v3v3(q, s, e1); + *lambda = f * dot_v3v3(e2, q); + if ((*lambda < 0.0)) return 0; + + u = f * dot_v3v3(s, p); + v = f * dot_v3v3(d, q); + + if (u < 0) du = u; + if (u > 1) du = u - 1; + if (v < 0) dv = v; + if (v > 1) dv = v - 1; + if (u > 0 && v > 0 && u + v > 1) + { + float t = u + v - 1; + du = u - t/2; + dv = v - t/2; + } + + mul_v3_fl(e1, du); + mul_v3_fl(e2, dv); + + if (dot_v3v3(e1, e1) + dot_v3v3(e2, e2) > threshold * threshold) + { + return 0; + } + + if(uv) { + uv[0]= u; + uv[1]= v; + } + + return 1; +} + + +/* Adapted from the paper by Kasper Fauerby */ +/* "Improved Collision detection and Response" */ +static int getLowestRoot(float a, float b, float c, float maxR, float* root) +{ + // Check if a solution exists + float determinant = b*b - 4.0f*a*c; + + // If determinant is negative it means no solutions. + if (determinant >= 0.0f) + { + // calculate the two roots: (if determinant == 0 then + // x1==x2 but let’s disregard that slight optimization) + float sqrtD = (float)sqrt(determinant); + float r1 = (-b - sqrtD) / (2.0f*a); + float r2 = (-b + sqrtD) / (2.0f*a); + + // Sort so x1 <= x2 + if (r1 > r2) + SWAP(float, r1, r2); + + // Get lowest root: + if (r1 > 0.0f && r1 < maxR) + { + *root = r1; + return 1; + } + + // It is possible that we want x2 - this can happen + // if x1 < 0 + if (r2 > 0.0f && r2 < maxR) + { + *root = r2; + return 1; + } + } + // No (valid) solutions + return 0; +} + +int isect_sweeping_sphere_tri_v3(float p1[3], float p2[3], float radius, float v0[3], float v1[3], float v2[3], float *lambda, float *ipoint) +{ + 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 newLambda; + int found_by_sweep=0; + + sub_v3_v3v3(e1,v1,v0); + sub_v3_v3v3(e2,v2,v0); + sub_v3_v3v3(vel,p2,p1); + +/*---test plane of tri---*/ + cross_v3_v3v3(nor,e1,e2); + normalize_v3(nor); + + /* flip normal */ + if(dot_v3v3(nor,vel)>0.0f) negate_v3(nor); + + a=dot_v3v3(p1,nor)-dot_v3v3(v0,nor); + nordotv=dot_v3v3(nor,vel); + + if (fabs(nordotv) < 0.000001) + { + if(fabs(a)>=radius) + { + return 0; + } + } + else + { + float t0=(-a+radius)/nordotv; + float t1=(-a-radius)/nordotv; + + if(t0>t1) + SWAP(float, t0, t1); + + if(t0>1.0f || t1<0.0f) return 0; + + /* clamp to [0,1] */ + CLAMP(t0, 0.0f, 1.0f); + CLAMP(t1, 0.0f, 1.0f); + + /*---test inside of tri---*/ + /* plane intersection point */ + + point[0] = p1[0] + vel[0]*t0 - nor[0]*radius; + point[1] = p1[1] + vel[1]*t0 - nor[1]*radius; + point[2] = p1[2] + vel[2]*t0 - nor[2]*radius; + + + /* is the point in the tri? */ + a=dot_v3v3(e1,e1); + b=dot_v3v3(e1,e2); + c=dot_v3v3(e2,e2); + + sub_v3_v3v3(temp,point,v0); + d=dot_v3v3(temp,e1); + e=dot_v3v3(temp,e2); + + x=d*c-e*b; + y=e*a-d*b; + z=x+y-(a*c-b*b); + + + if(z <= 0.0f && (x >= 0.0f && y >= 0.0f)) + { + //(((unsigned int)z)& ~(((unsigned int)x)|((unsigned int)y))) & 0x80000000){ + *lambda=t0; + copy_v3_v3(ipoint,point); + return 1; + } + } + + + *lambda=1.0f; + +/*---test points---*/ + a=vel2=dot_v3v3(vel,vel); + + /*v0*/ + sub_v3_v3v3(temp,p1,v0); + b=2.0f*dot_v3v3(vel,temp); + c=dot_v3v3(temp,temp)-radius2; + + if(getLowestRoot(a, b, c, *lambda, lambda)) + { + copy_v3_v3(ipoint,v0); + found_by_sweep=1; + } + + /*v1*/ + sub_v3_v3v3(temp,p1,v1); + b=2.0f*dot_v3v3(vel,temp); + c=dot_v3v3(temp,temp)-radius2; + + if(getLowestRoot(a, b, c, *lambda, lambda)) + { + copy_v3_v3(ipoint,v1); + found_by_sweep=1; + } + + /*v2*/ + sub_v3_v3v3(temp,p1,v2); + b=2.0f*dot_v3v3(vel,temp); + c=dot_v3v3(temp,temp)-radius2; + + if(getLowestRoot(a, b, c, *lambda, lambda)) + { + copy_v3_v3(ipoint,v2); + found_by_sweep=1; + } + +/*---test edges---*/ + sub_v3_v3v3(e3,v2,v1); //wasnt yet calculated + + + /*e1*/ + sub_v3_v3v3(bv,v0,p1); + + elen2 = dot_v3v3(e1,e1); + edotv = dot_v3v3(e1,vel); + edotbv = dot_v3v3(e1,bv); + + a=elen2*(-dot_v3v3(vel,vel))+edotv*edotv; + b=2.0f*(elen2*dot_v3v3(vel,bv)-edotv*edotbv); + c=elen2*(radius2-dot_v3v3(bv,bv))+edotbv*edotbv; + + if(getLowestRoot(a, b, c, *lambda, &newLambda)) + { + e=(edotv*newLambda-edotbv)/elen2; + + if(e >= 0.0f && e <= 1.0f) + { + *lambda = newLambda; + copy_v3_v3(ipoint,e1); + mul_v3_fl(ipoint,e); + add_v3_v3v3(ipoint,ipoint,v0); + found_by_sweep=1; + } + } + + /*e2*/ + /*bv is same*/ + elen2 = dot_v3v3(e2,e2); + edotv = dot_v3v3(e2,vel); + edotbv = dot_v3v3(e2,bv); + + a=elen2*(-dot_v3v3(vel,vel))+edotv*edotv; + b=2.0f*(elen2*dot_v3v3(vel,bv)-edotv*edotbv); + c=elen2*(radius2-dot_v3v3(bv,bv))+edotbv*edotbv; + + if(getLowestRoot(a, b, c, *lambda, &newLambda)) + { + e=(edotv*newLambda-edotbv)/elen2; + + if(e >= 0.0f && e <= 1.0f) + { + *lambda = newLambda; + copy_v3_v3(ipoint,e2); + mul_v3_fl(ipoint,e); + add_v3_v3v3(ipoint,ipoint,v0); + found_by_sweep=1; + } + } + + /*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,v1,p1); + elen2 = dot_v3v3(e3,e3); + edotv = dot_v3v3(e3,vel); + edotbv = dot_v3v3(e3,bv); + + a=elen2*(-dot_v3v3(vel,vel))+edotv*edotv; + b=2.0f*(elen2*dot_v3v3(vel,bv)-edotv*edotbv); + c=elen2*(radius2-dot_v3v3(bv,bv))+edotbv*edotbv; + + if(getLowestRoot(a, b, c, *lambda, &newLambda)) + { + e=(edotv*newLambda-edotbv)/elen2; + + if(e >= 0.0f && e <= 1.0f) + { + *lambda = newLambda; + copy_v3_v3(ipoint,e3); + mul_v3_fl(ipoint,e); + add_v3_v3v3(ipoint,ipoint,v1); + found_by_sweep=1; + } + } + + + return found_by_sweep; +} +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) +{ + float p[3], e1[3], e2[3]; + float u, v, f; + int a0=axis, a1=(axis+1)%3, a2=(axis+2)%3; + + //return isect_line_tri_v3(p1,p2,v0,v1,v2,lambda); + + ///* first a simple bounding box test */ + //if(MIN3(v0[a1],v1[a1],v2[a1]) > p1[a1]) return 0; + //if(MIN3(v0[a2],v1[a2],v2[a2]) > p1[a2]) return 0; + //if(MAX3(v0[a1],v1[a1],v2[a1]) < p1[a1]) return 0; + //if(MAX3(v0[a2],v1[a2],v2[a2]) < p1[a2]) return 0; + + ///* then a full intersection test */ + + sub_v3_v3v3(e1,v1,v0); + sub_v3_v3v3(e2,v2,v0); + sub_v3_v3v3(p,v0,p1); + + f= (e2[a1]*e1[a2]-e2[a2]*e1[a1]); + if ((f > -0.000001) && (f < 0.000001)) return 0; + + v= (p[a2]*e1[a1]-p[a1]*e1[a2])/f; + if ((v < 0.0)||(v > 1.0)) return 0; + + f= e1[a1]; + if((f > -0.000001) && (f < 0.000001)){ + f= e1[a2]; + if((f > -0.000001) && (f < 0.000001)) return 0; + u= (-p[a2]-v*e2[a2])/f; + } + else + u= (-p[a1]-v*e2[a1])/f; + + if ((u < 0.0)||((u + v) > 1.0)) return 0; + + *lambda = (p[a0]+u*e1[a0]+v*e2[a0])/(p2[a0]-p1[a0]); + + if ((*lambda < 0.0)||(*lambda > 1.0)) return 0; + + return 1; +} + +/* Returns the number of point of interests + * 0 - lines are colinear + * 1 - lines are coplanar, i1 is set to intersection + * 2 - i1 and i2 are the nearest points on line 1 (v1, v2) and line 2 (v3, v4) respectively + * */ +int isect_line_line_v3(float v1[3], float v2[3], float v3[3], float v4[3], float i1[3], float i2[3]) +{ + float a[3], b[3], c[3], ab[3], cb[3], dir1[3], dir2[3]; + float d; + + sub_v3_v3v3(c, v3, v1); + sub_v3_v3v3(a, v2, v1); + sub_v3_v3v3(b, v4, v3); + + copy_v3_v3(dir1, a); + normalize_v3(dir1); + copy_v3_v3(dir2, b); + normalize_v3(dir2); + d = dot_v3v3(dir1, dir2); + if (d == 1.0f || d == -1.0f) { + /* colinear */ + return 0; + } + + cross_v3_v3v3(ab, a, b); + d = dot_v3v3(c, ab); + + /* test if the two lines are coplanar */ + if (d > -0.000001f && d < 0.000001f) { + cross_v3_v3v3(cb, c, b); + + mul_v3_fl(a, dot_v3v3(cb, ab) / dot_v3v3(ab, ab)); + add_v3_v3v3(i1, v1, a); + copy_v3_v3(i2, i1); + + return 1; /* one intersection only */ + } + /* if not */ + else { + float n[3], t[3]; + float v3t[3], v4t[3]; + sub_v3_v3v3(t, v1, v3); + + /* offset between both plane where the lines lies */ + cross_v3_v3v3(n, a, b); + project_v3_v3v3(t, t, n); + + /* for the first line, offset the second line until it is coplanar */ + add_v3_v3v3(v3t, v3, t); + add_v3_v3v3(v4t, v4, t); + + sub_v3_v3v3(c, v3t, v1); + sub_v3_v3v3(a, v2, v1); + sub_v3_v3v3(b, v4t, v3t); + + cross_v3_v3v3(ab, a, b); + cross_v3_v3v3(cb, c, b); + + mul_v3_fl(a, dot_v3v3(cb, ab) / dot_v3v3(ab, ab)); + add_v3_v3v3(i1, v1, a); + + /* for the second line, just substract the offset from the first intersection point */ + sub_v3_v3v3(i2, i1, t); + + return 2; /* two nearest points */ + } +} + +/* Intersection point strictly between the two lines + * 0 when no intersection is found + * */ +int isect_line_line_strict_v3(float v1[3], float v2[3], float v3[3], float v4[3], float vi[3], float *lambda) +{ + 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); + sub_v3_v3v3(b, v4, v3); + + copy_v3_v3(dir1, a); + normalize_v3(dir1); + copy_v3_v3(dir2, b); + normalize_v3(dir2); + d = dot_v3v3(dir1, dir2); + if (d == 1.0f || d == -1.0f || d == 0) { + /* colinear or one vector is zero-length*/ + return 0; + } + + d1 = d; + + cross_v3_v3v3(ab, a, b); + d = dot_v3v3(c, ab); + + /* test if the two lines are coplanar */ + if (d > -0.000001f && d < 0.000001f) { + float f1, f2; + cross_v3_v3v3(cb, c, b); + cross_v3_v3v3(ca, c, a); + + f1 = dot_v3v3(cb, ab) / dot_v3v3(ab, ab); + f2 = dot_v3v3(ca, ab) / dot_v3v3(ab, ab); + + if (f1 >= 0 && f1 <= 1 && + f2 >= 0 && f2 <= 1) + { + mul_v3_fl(a, f1); + add_v3_v3v3(vi, v1, a); + + if (lambda != NULL) + { + *lambda = f1; + } + + return 1; /* intersection found */ + } + else + { + return 0; + } + } + else + { + return 0; + } +} + +int isect_aabb_aabb_v3(float min1[3], float max1[3], float min2[3], float max2[3]) +{ + return (min1[0]<max2[0] && min1[1]<max2[1] && min1[2]<max2[2] && + min2[0]<max1[0] && min2[1]<max1[1] && min2[2]<max1[2]); +} + +/* 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 h[3],u[3],lambda; + sub_v3_v3v3(u, l2, l1); + sub_v3_v3v3(h, p, l1); + lambda =dot_v3v3(u,h)/dot_v3v3(u,u); + cp[0] = l1[0] + u[0] * lambda; + cp[1] = l1[1] + u[1] * lambda; + cp[2] = l1[2] + u[2] * 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]) +{ + float h[3],u[3]; + sub_v3_v3v3(u, l2, l1); + sub_v3_v3v3(h, p, l1); + return(dot_v3v3(u,h)/dot_v3v3(u,u)); +} +#endif + +/* Similar to LineIntersectsTriangleUV, except it operates on a quad and in 2d, assumes point is in quad */ +void isect_point_quad_uv_v2(float v0[2], float v1[2], float v2[2], float v3[2], float pt[2], float *uv) +{ + float x0,y0, x1,y1, wtot, v2d[2], w1, w2; + + /* used for paralelle lines */ + float pt3d[3], l1[3], l2[3], pt_on_line[3]; + + /* compute 2 edges of the quad intersection point */ + if (IsectLLPt2Df(v0[0],v0[1],v1[0],v1[1], v2[0],v2[1],v3[0],v3[1], &x0,&y0) == 1) { + /* the intersection point between the quad-edge intersection and the point in the quad we want the uv's for */ + /* should never be paralle !! */ + /*printf("\tnot paralelle 1\n");*/ + IsectLLPt2Df(pt[0],pt[1],x0,y0, v0[0],v0[1],v3[0],v3[1], &x1,&y1); + + /* Get the weights from the new intersection point, to each edge */ + v2d[0] = x1-v0[0]; + v2d[1] = y1-v0[1]; + w1 = len_v2(v2d); + + v2d[0] = x1-v3[0]; /* some but for the other vert */ + v2d[1] = y1-v3[1]; + w2 = len_v2(v2d); + wtot = w1+w2; + /*w1 = w1/wtot;*/ + /*w2 = w2/wtot;*/ + uv[0] = w1/wtot; + } else { + /* lines are paralelle, lambda_cp_line_ex is 3d grrr */ + /*printf("\tparalelle1\n");*/ + pt3d[0] = pt[0]; + pt3d[1] = pt[1]; + pt3d[2] = l1[2] = l2[2] = 0.0f; + + l1[0] = v0[0]; l1[1] = v0[1]; + l2[0] = v1[0]; l2[1] = v1[1]; + closest_to_line_v3(pt_on_line,pt3d, l1, l2); + v2d[0] = pt[0]-pt_on_line[0]; /* same, for the other vert */ + v2d[1] = pt[1]-pt_on_line[1]; + w1 = len_v2(v2d); + + l1[0] = v2[0]; l1[1] = v2[1]; + l2[0] = v3[0]; l2[1] = v3[1]; + closest_to_line_v3(pt_on_line,pt3d, l1, l2); + v2d[0] = pt[0]-pt_on_line[0]; /* same, for the other vert */ + v2d[1] = pt[1]-pt_on_line[1]; + w2 = len_v2(v2d); + wtot = w1+w2; + uv[0] = w1/wtot; + } + + /* Same as above to calc the uv[1] value, alternate calculation */ + + 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*/ + /* never paralle if above was not */ + /*printf("\tnot paralelle2\n");*/ + IsectLLPt2Df(pt[0],pt[1],x0,y0, v0[0],v0[1],v1[0],v1[1], &x1,&y1);/* was v0,v3 now v0,v1*/ + + v2d[0] = x1-v0[0]; + v2d[1] = y1-v0[1]; + w1 = len_v2(v2d); + + v2d[0] = x1-v1[0]; + v2d[1] = y1-v1[1]; + w2 = len_v2(v2d); + wtot = w1+w2; + uv[1] = w1/wtot; + } else { + /* lines are paralelle, lambda_cp_line_ex is 3d grrr */ + /*printf("\tparalelle2\n");*/ + pt3d[0] = pt[0]; + pt3d[1] = pt[1]; + pt3d[2] = l1[2] = l2[2] = 0.0f; + + + l1[0] = v0[0]; l1[1] = v0[1]; + l2[0] = v3[0]; l2[1] = v3[1]; + closest_to_line_v3(pt_on_line,pt3d, l1, l2); + v2d[0] = pt[0]-pt_on_line[0]; /* some but for the other vert */ + v2d[1] = pt[1]-pt_on_line[1]; + w1 = len_v2(v2d); + + l1[0] = v1[0]; l1[1] = v1[1]; + l2[0] = v2[0]; l2[1] = v2[1]; + closest_to_line_v3(pt_on_line,pt3d, l1, l2); + v2d[0] = pt[0]-pt_on_line[0]; /* some but for the other vert */ + v2d[1] = pt[1]-pt_on_line[1]; + w2 = len_v2(v2d); + wtot = w1+w2; + uv[1] = w1/wtot; + } + /* may need to flip UV's here */ +} + +/* same as above but does tri's and quads, tri's are a bit of a hack */ +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) +{ + if (isquad) { + isect_point_quad_uv_v2(v0, v1, v2, v3, pt, uv); + } + else { + /* not for quads, use for our abuse of LineIntersectsTriangleUV */ + float p1_3d[3], p2_3d[3], v0_3d[3], v1_3d[3], v2_3d[3], lambda; + + p1_3d[0] = p2_3d[0] = uv[0]; + p1_3d[1] = p2_3d[1] = uv[1]; + p1_3d[2] = 1.0f; + p2_3d[2] = -1.0f; + v0_3d[2] = v1_3d[2] = v2_3d[2] = 0.0; + + /* generate a new fuv, (this is possibly a non optimal solution, + * since we only need 2d calculation but use 3d func's) + * + * this method makes an imaginary triangle in 2d space using the UV's from the derived mesh face + * Then find new uv coords using the fuv and this face with LineIntersectsTriangleUV. + * This means the new values will be correct in relation to the derived meshes face. + */ + copy_v2_v2(v0_3d, v0); + copy_v2_v2(v1_3d, v1); + copy_v2_v2(v2_3d, v2); + + /* Doing this in 3D is not nice */ + isect_line_tri_v3(p1_3d, p2_3d, v0_3d, v1_3d, v2_3d, &lambda, uv); + } +} + +#if 0 +int isect_point_tri_v2(float v1[2], float v2[2], float v3[2], float pt[2]) +{ + float inp1, inp2, inp3; + + inp1= (v2[0]-v1[0])*(v1[1]-pt[1]) + (v1[1]-v2[1])*(v1[0]-pt[0]); + inp2= (v3[0]-v2[0])*(v2[1]-pt[1]) + (v2[1]-v3[1])*(v2[0]-pt[0]); + inp3= (v1[0]-v3[0])*(v3[1]-pt[1]) + (v3[1]-v1[1])*(v3[0]-pt[0]); + + if(inp1<=0.0f && inp2<=0.0f && inp3<=0.0f) return 1; + if(inp1>=0.0f && inp2>=0.0f && inp3>=0.0f) return 1; + + return 0; +} +#endif + +#if 0 +int isect_point_tri_v2(float v0[2], float v1[2], float v2[2], float pt[2]) +{ + /* not for quads, use for our abuse of LineIntersectsTriangleUV */ + float p1_3d[3], p2_3d[3], v0_3d[3], v1_3d[3], v2_3d[3]; + /* not used */ + float lambda, uv[3]; + + p1_3d[0] = p2_3d[0] = uv[0]= pt[0]; + p1_3d[1] = p2_3d[1] = uv[1]= uv[2]= pt[1]; + p1_3d[2] = 1.0f; + p2_3d[2] = -1.0f; + v0_3d[2] = v1_3d[2] = v2_3d[2] = 0.0; + + /* generate a new fuv, (this is possibly a non optimal solution, + * since we only need 2d calculation but use 3d func's) + * + * this method makes an imaginary triangle in 2d space using the UV's from the derived mesh face + * Then find new uv coords using the fuv and this face with LineIntersectsTriangleUV. + * This means the new values will be correct in relation to the derived meshes face. + */ + copy_v2_v2(v0_3d, v0); + copy_v2_v2(v1_3d, v1); + copy_v2_v2(v2_3d, v2); + + /* Doing this in 3D is not nice */ + return isect_line_tri_v3(p1_3d, p2_3d, v0_3d, v1_3d, v2_3d, &lambda, uv); +} +#endif + +/* + + x1,y2 + | \ + | \ .(a,b) + | \ + x1,y1-- x2,y1 + +*/ +int isect_point_tri_v2_int(int x1, int y1, int x2, int y2, int a, int b) +{ + float v1[2], v2[2], v3[2], p[2]; + + v1[0]= (float)x1; + v1[1]= (float)y1; + + v2[0]= (float)x1; + v2[1]= (float)y2; + + v3[0]= (float)x2; + v3[1]= (float)y1; + + p[0]= (float)a; + p[1]= (float)b; + + return isect_point_tri_v2(v1, v2, v3, p); +} + +static int point_in_slice(float p[3], float v1[3], float l1[3], float l2[3]) +{ +/* +what is a slice ? +some maths: +a line including l1,l2 and a point not on the line +define a subset of R3 delimeted by planes parallel to the line and orthogonal +to the (point --> line) distance vector,one plane on the line one on the point, +the room inside usually is rather small compared to R3 though still infinte +useful for restricting (speeding up) searches +e.g. all points of triangular prism are within the intersection of 3 'slices' +onother trivial case : cube +but see a 'spat' which is a deformed cube with paired parallel planes needs only 3 slices too +*/ + float h,rp[3],cp[3],q[3]; + + closest_to_line_v3(cp,v1,l1,l2); + sub_v3_v3v3(q,cp,v1); + + sub_v3_v3v3(rp,p,v1); + h=dot_v3v3(q,rp)/dot_v3v3(q,q); + if (h < 0.0f || h > 1.0f) return 0; + return 1; +} + +#if 0 +/*adult sister defining the slice planes by the origin and the normal +NOTE |normal| may not be 1 but defining the thickness of the slice*/ +static int point_in_slice_as(float p[3],float origin[3],float normal[3]) +{ + float h,rp[3]; + sub_v3_v3v3(rp,p,origin); + h=dot_v3v3(normal,rp)/dot_v3v3(normal,normal); + if (h < 0.0f || h > 1.0f) return 0; + return 1; +} + +/*mama (knowing the squared lenght of the normal)*/ +static int point_in_slice_m(float p[3],float origin[3],float normal[3],float lns) +{ + float h,rp[3]; + sub_v3_v3v3(rp,p,origin); + h=dot_v3v3(normal,rp)/lns; + if (h < 0.0f || h > 1.0f) return 0; + return 1; +} +#endif + +int isect_point_tri_prism_v3(float p[3], float v1[3], float v2[3], float v3[3]) +{ + if(!point_in_slice(p,v1,v2,v3)) return 0; + if(!point_in_slice(p,v2,v3,v1)) return 0; + if(!point_in_slice(p,v3,v1,v2)) return 0; + return 1; +} + +/****************************** Interpolation ********************************/ + +static float tri_signed_area(float *v1, float *v2, float *v3, int i, int j) +{ + return 0.5f*((v1[i]-v2[i])*(v2[j]-v3[j]) + (v1[j]-v2[j])*(v3[i]-v2[i])); +} + +static int barycentric_weights(float *v1, float *v2, float *v3, float *co, float *n, float *w) +{ + float xn, yn, zn, a1, a2, a3, asum; + short i, j; + + /* find best projection of face XY, XZ or YZ: barycentric weights of + the 2d projected coords are the same and faster to compute */ + xn= (float)fabs(n[0]); + yn= (float)fabs(n[1]); + zn= (float)fabs(n[2]); + if(zn>=xn && zn>=yn) {i= 0; j= 1;} + else if(yn>=xn && yn>=zn) {i= 0; j= 2;} + else {i= 1; j= 2;} + + a1= tri_signed_area(v2, v3, co, i, j); + a2= tri_signed_area(v3, v1, co, i, j); + a3= tri_signed_area(v1, v2, co, i, j); + + asum= a1 + a2 + a3; + + if (fabs(asum) < FLT_EPSILON) { + /* zero area triangle */ + w[0]= w[1]= w[2]= 1.0f/3.0f; + return 1; + } + + asum= 1.0f/asum; + w[0]= a1*asum; + w[1]= a2*asum; + w[2]= a3*asum; + + return 0; +} + +void interp_weights_face_v3(float *w,float *v1, float *v2, float *v3, float *v4, float *co) +{ + float w2[3]; + + w[0]= w[1]= w[2]= w[3]= 0.0f; + + /* first check for exact match */ + if(equals_v3v3(co, v1)) + w[0]= 1.0f; + else if(equals_v3v3(co, v2)) + w[1]= 1.0f; + else if(equals_v3v3(co, v3)) + w[2]= 1.0f; + else if(v4 && equals_v3v3(co, v4)) + w[3]= 1.0f; + else { + /* otherwise compute barycentric interpolation weights */ + float n1[3], n2[3], n[3]; + int degenerate; + + sub_v3_v3v3(n1, v1, v3); + if (v4) { + sub_v3_v3v3(n2, v2, v4); + } + else { + sub_v3_v3v3(n2, v2, v3); + } + cross_v3_v3v3(n, n1, n2); + + /* OpenGL seems to split this way, so we do too */ + if (v4) { + degenerate= barycentric_weights(v1, v2, v4, co, n, w); + SWAP(float, w[2], w[3]); + + if(degenerate || (w[0] < 0.0f)) { + /* if w[1] is negative, co is on the other side of the v1-v3 edge, + so we interpolate using the other triangle */ + degenerate= barycentric_weights(v2, v3, v4, co, n, w2); + + if(!degenerate) { + w[0]= 0.0f; + w[1]= w2[0]; + w[2]= w2[1]; + w[3]= w2[2]; + } + } + } + else + barycentric_weights(v1, v2, v3, co, n, w); + } +} + +/* Mean value weights - smooth interpolation weights for polygons with + * more than 3 vertices */ +static float mean_value_half_tan(float *v1, float *v2, float *v3) +{ + float d2[3], d3[3], cross[3], area, dot, len; + + sub_v3_v3v3(d2, v2, v1); + sub_v3_v3v3(d3, v3, v1); + cross_v3_v3v3(cross, d2, d3); + + area= len_v3(cross); + dot= dot_v3v3(d2, d3); + len= len_v3(d2)*len_v3(d3); + + if(area == 0.0f) + return 0.0f; + else + return (len - dot)/area; +} + +void interp_weights_poly_v3(float *w,float v[][3], int n, float *co) +{ + float totweight, t1, t2, len, *vmid, *vprev, *vnext; + int i; + + totweight= 0.0f; + + for(i=0; i<n; i++) { + vmid= v[i]; + vprev= (i == 0)? v[n-1]: v[i-1]; + vnext= (i == n-1)? v[0]: v[i+1]; + + t1= mean_value_half_tan(co, vprev, vmid); + t2= mean_value_half_tan(co, vmid, vnext); + + len= len_v3v3(co, vmid); + w[i]= (t1+t2)/len; + totweight += w[i]; + } + + if(totweight != 0.0f) + for(i=0; i<n; i++) + w[i] /= totweight; +} + +/* (x1,v1)(t1=0)------(x2,v2)(t2=1), 0<t<1 --> (x,v)(t) */ +void interp_cubic_v3(float *x, float *v,float *x1, float *v1, float *x2, float *v2, float t) +{ + float a[3],b[3]; + float t2= t*t; + float t3= t2*t; + + /* cubic interpolation */ + a[0]= v1[0] + v2[0] + 2*(x1[0] - x2[0]); + a[1]= v1[1] + v2[1] + 2*(x1[1] - x2[1]); + a[2]= v1[2] + v2[2] + 2*(x1[2] - x2[2]); + + b[0]= -2*v1[0] - v2[0] - 3*(x1[0] - x2[0]); + b[1]= -2*v1[1] - v2[1] - 3*(x1[1] - x2[1]); + b[2]= -2*v1[2] - v2[2] - 3*(x1[2] - x2[2]); + + x[0]= a[0]*t3 + b[0]*t2 + v1[0]*t + x1[0]; + x[1]= a[1]*t3 + b[1]*t2 + v1[1]*t + x1[1]; + x[2]= a[2]*t3 + b[2]*t2 + v1[2]*t + x1[2]; + + v[0]= 3*a[0]*t2 + 2*b[0]*t + v1[0]; + v[1]= 3*a[1]*t2 + 2*b[1]*t + v1[1]; + v[2]= 3*a[2]*t2 + 2*b[2]*t + v1[2]; +} + +/***************************** View & Projection *****************************/ + +void orthographic_m4(float matrix[][4],float left, float right, float bottom, float top, float nearClip, float farClip) +{ + float Xdelta, Ydelta, Zdelta; + + Xdelta = right - left; + Ydelta = top - bottom; + Zdelta = farClip - nearClip; + if (Xdelta == 0.0 || Ydelta == 0.0 || Zdelta == 0.0) { + return; + } + unit_m4(matrix); + matrix[0][0] = 2.0f/Xdelta; + matrix[3][0] = -(right + left)/Xdelta; + matrix[1][1] = 2.0f/Ydelta; + matrix[3][1] = -(top + bottom)/Ydelta; + matrix[2][2] = -2.0f/Zdelta; /* note: negate Z */ + matrix[3][2] = -(farClip + nearClip)/Zdelta; +} + +void perspective_m4(float mat[][4],float left, float right, float bottom, float top, float nearClip, float farClip) +{ + float Xdelta, Ydelta, Zdelta; + + Xdelta = right - left; + Ydelta = top - bottom; + Zdelta = farClip - nearClip; + + if (Xdelta == 0.0 || Ydelta == 0.0 || Zdelta == 0.0) { + return; + } + mat[0][0] = nearClip * 2.0f/Xdelta; + mat[1][1] = nearClip * 2.0f/Ydelta; + mat[2][0] = (right + left)/Xdelta; /* note: negate Z */ + mat[2][1] = (top + bottom)/Ydelta; + mat[2][2] = -(farClip + nearClip)/Zdelta; + mat[2][3] = -1.0f; + mat[3][2] = (-2.0f * nearClip * farClip)/Zdelta; + mat[0][1] = mat[0][2] = mat[0][3] = + mat[1][0] = mat[1][2] = mat[1][3] = + mat[3][0] = mat[3][1] = mat[3][3] = 0.0; + +} + +static void i_multmatrix(float icand[][4], float Vm[][4]) +{ + int row, col; + float temp[4][4]; + + for(row=0 ; row<4 ; row++) + for(col=0 ; col<4 ; col++) + temp[row][col] = icand[row][0] * Vm[0][col] + + icand[row][1] * Vm[1][col] + + icand[row][2] * Vm[2][col] + + icand[row][3] * Vm[3][col]; + copy_m4_m4(Vm, temp); +} + + +void polarview_m4(float Vm[][4],float dist, float azimuth, float incidence, float twist) +{ + + unit_m4(Vm); + + translate_m4(Vm,0.0, 0.0, -dist); + rotate_m4(Vm,'z',-twist); + rotate_m4(Vm,'x',-incidence); + rotate_m4(Vm,'z',-azimuth); +} + +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]; + + unit_m4(mat); + unit_m4(mat1); + + rotate_m4(mat,'z',-twist); + + dx = px - vx; + dy = py - vy; + dz = pz - vz; + hyp = dx * dx + dz * dz; /* hyp squared */ + hyp1 = (float)sqrt(dy*dy + hyp); + hyp = (float)sqrt(hyp); /* the real hyp */ + + if (hyp1 != 0.0) { /* rotate X */ + sine = -dy / hyp1; + cosine = hyp /hyp1; + } else { + sine = 0; + cosine = 1.0f; + } + mat1[1][1] = cosine; + mat1[1][2] = sine; + mat1[2][1] = -sine; + mat1[2][2] = cosine; + + i_multmatrix(mat1, mat); + + mat1[1][1] = mat1[2][2] = 1.0f; /* be careful here to reinit */ + mat1[1][2] = mat1[2][1] = 0.0; /* those modified by the last */ + + /* paragraph */ + if (hyp != 0.0f) { /* rotate Y */ + sine = dx / hyp; + cosine = -dz / hyp; + } else { + sine = 0; + cosine = 1.0f; + } + mat1[0][0] = cosine; + mat1[0][2] = -sine; + mat1[2][0] = sine; + mat1[2][2] = cosine; + + i_multmatrix(mat1, mat); + translate_m4(mat,-vx,-vy,-vz); /* translate viewpoint to origin */ +} + +/********************************** Mapping **********************************/ + +void map_to_tube(float *u, float *v,float x, float y, float z) +{ + float len; + + *v = (z + 1.0f) / 2.0f; + + len= (float)sqrt(x*x+y*y); + if(len > 0.0f) + *u = (float)((1.0 - (atan2(x/len,y/len) / M_PI)) / 2.0); + else + *v = *u = 0.0f; /* to avoid un-initialized variables */ +} + +void map_to_sphere(float *u, float *v,float x, float y, float z) +{ + float len; + + len= (float)sqrt(x*x+y*y+z*z); + if(len > 0.0f) { + if(x==0.0f && y==0.0f) *u= 0.0f; /* othwise domain error */ + else *u = (float)((1.0 - (float)atan2(x,y) / M_PI) / 2.0); + + z/=len; + *v = 1.0f - (float)saacos(z)/(float)M_PI; + } else { + *v = *u = 0.0f; /* to avoid un-initialized variables */ + } +} + +/********************************************************/ + +/* Tangents */ + +/* For normal map tangents we need to detect uv boundaries, and only average + * tangents in case the uvs are connected. Alternative would be to store 1 + * tangent per face rather than 4 per face vertex, but that's not compatible + * with games */ + + +/* from BKE_mesh.h */ +#define STD_UV_CONNECT_LIMIT 0.0001f + +#if 0 +void sum_or_add_vertex_tangent(void *arena, VertexTangent **vtang, float *tang, float *uv) +{ + VertexTangent *vt; + + /* find a tangent with connected uvs */ + for(vt= *vtang; vt; vt=vt->next) { + if(fabs(uv[0]-vt->uv[0]) < STD_UV_CONNECT_LIMIT && fabs(uv[1]-vt->uv[1]) < STD_UV_CONNECT_LIMIT) { + add_v3_v3v3(vt->tang, vt->tang, tang); + return; + } + } + + /* if not found, append a new one */ + vt= BLI_memarena_alloc((MemArena *)arena, sizeof(VertexTangent)); + copy_v3_v3(vt->tang, tang); + vt->uv[0]= uv[0]; + vt->uv[1]= uv[1]; + + if(*vtang) + vt->next= *vtang; + *vtang= vt; +} + +float *find_vertex_tangent(VertexTangent *vtang, float *uv) +{ + VertexTangent *vt; + static float nulltang[3] = {0.0f, 0.0f, 0.0f}; + + for(vt= vtang; vt; vt=vt->next) + if(fabs(uv[0]-vt->uv[0]) < STD_UV_CONNECT_LIMIT && fabs(uv[1]-vt->uv[1]) < STD_UV_CONNECT_LIMIT) + return vt->tang; + + return nulltang; /* shouldn't happen, except for nan or so */ +} + +void tangent_from_uv(float *uv1, float *uv2, float *uv3, float *co1, float *co2, float *co3, float *n, float *tang) +{ + float tangv[3], ct[3], e1[3], e2[3], s1, t1, s2, t2, det; + + s1= uv2[0] - uv1[0]; + s2= uv3[0] - uv1[0]; + t1= uv2[1] - uv1[1]; + t2= uv3[1] - uv1[1]; + det= 1.0f / (s1 * t2 - s2 * t1); + + /* normals in render are inversed... */ + sub_v3_v3v3(e1, co1, co2); + sub_v3_v3v3(e2, co1, co3); + tang[0] = (t2*e1[0] - t1*e2[0])*det; + tang[1] = (t2*e1[1] - t1*e2[1])*det; + tang[2] = (t2*e1[2] - t1*e2[2])*det; + tangv[0] = (s1*e2[0] - s2*e1[0])*det; + tangv[1] = (s1*e2[1] - s2*e1[1])*det; + tangv[2] = (s1*e2[2] - s2*e1[2])*det; + cross_v3_v3v3(ct, tang, tangv); + + /* check flip */ + if ((ct[0]*n[0] + ct[1]*n[1] + ct[2]*n[2]) < 0.0f) + negate_v3(tang); +} +#endif + |