diff options
author | Brecht Van Lommel <brechtvanlommel@pandora.be> | 2008-11-13 00:16:53 +0300 |
---|---|---|
committer | Brecht Van Lommel <brechtvanlommel@pandora.be> | 2008-11-13 00:16:53 +0300 |
commit | bdfe7d89e2f1292644577972c716931b4ce3c6c3 (patch) | |
tree | d00eb50b749cb001e2b08272c91791e66740b05d /source/blender/blenlib/intern/arithb.c | |
parent | 78a1c27c4a6abe0ed31ca93ad21910f3df04da56 (diff) | |
parent | 7e4db234cee71ead34ee81a12e27da4bd548eb4b (diff) |
Merge of trunk into blender 2.5:
svn merge https://svn.blender.org/svnroot/bf-blender/trunk/blender -r12987:17416
Issues:
* GHOST/X11 had conflicting changes. Some code was added in 2.5, which was
later added in trunk also, but reverted partially, specifically revision
16683. I have left out this reversion in the 2.5 branch since I think it is
needed there.
http://projects.blender.org/plugins/scmsvn/viewcvs.php?view=rev&root=bf-blender&revision=16683
* Scons had various conflicting changes, I decided to go with trunk version
for everything except priorities and some library renaming.
* In creator.c, there were various fixes and fixes for fixes related to the -w
-W and -p options. In 2.5 -w and -W is not coded yet, and -p is done
differently. Since this is changed so much, and I don't think those fixes
would be needed in 2.5, I've left them out.
* Also in creator.c: there was code for a python bugfix where the screen was not
initialized when running with -P. The code that initializes the screen there
I had to disable, that can't work in 2.5 anymore but left it commented as a
reminder.
Further I had to disable some new function calls. using src/ and python/, as
was done already in this branch, disabled function calls:
* bpath.c: error reporting
* BME_conversions.c: editmesh conversion functions.
* SHD_dynamic: disabled almost completely, there is no python/.
* KX_PythonInit.cpp and Ketsji/ build files: Mathutils is not there, disabled.
* text.c: clipboard copy call.
* object.c: OB_SUPPORT_MATERIAL.
* DerivedMesh.c and subsurf_ccg, stipple_quarttone.
Still to be done:
* Go over files and functions that were moved to a different location but could
still use changes that were done in trunk.
Diffstat (limited to 'source/blender/blenlib/intern/arithb.c')
-rw-r--r-- | source/blender/blenlib/intern/arithb.c | 761 |
1 files changed, 574 insertions, 187 deletions
diff --git a/source/blender/blenlib/intern/arithb.c b/source/blender/blenlib/intern/arithb.c index 735ff2e65f7..0f2226fc872 100644 --- a/source/blender/blenlib/intern/arithb.c +++ b/source/blender/blenlib/intern/arithb.c @@ -54,11 +54,13 @@ #include <stdio.h> #include "BLI_arithb.h" +#include "BLI_memarena.h" /* A few small defines. Keep'em local! */ #define SMALL_NUMBER 1.e-8 #define ABS(x) ((x) < 0 ? -(x) : (x)) #define SWAP(type, a, b) { type sw_ap; sw_ap=(a); (a)=(b); (b)=sw_ap; } +#define CLAMP(a, b, c) if((a)<(b)) (a)=(b); else if((a)>(c)) (a)=(c) #if defined(WIN32) || defined(__APPLE__) @@ -757,6 +759,28 @@ void Mat4MulSerie(float answ[][4], float m1[][4], } } +void Mat3BlendMat3(float out[][3], float dst[][3], float src[][3], float srcweight) +{ + float squat[4], dquat[4], fquat[4]; + float ssize[3], dsize[3], fsize[4]; + float rmat[3][3], smat[3][3]; + + Mat3ToQuat(dst, dquat); + Mat3ToSize(dst, dsize); + + Mat3ToQuat(src, squat); + Mat3ToSize(src, ssize); + + /* do blending */ + QuatInterpol(fquat, dquat, squat, srcweight); + VecLerpf(fsize, dsize, ssize, srcweight); + + /* compose new matrix */ + QuatToMat3(fquat, rmat); + SizeToMat3(fsize, smat); + Mat3MulMat3(out, rmat, smat); +} + void Mat4BlendMat4(float out[][4], float dst[][4], float src[][4], float srcweight) { float squat[4], dquat[4], fquat[4]; @@ -1002,6 +1026,19 @@ int FloatCompare( float *v1, float *v2, float limit) return 0; } +int FloatCompare4( float *v1, float *v2, float limit) +{ + + if( fabs(v1[0]-v2[0])<limit ) { + if( fabs(v1[1]-v2[1])<limit ) { + if( fabs(v1[2]-v2[2])<limit ) { + if( fabs(v1[3]-v2[3])<limit ) return 1; + } + } + } + return 0; +} + float FloatLerpf( float target, float origin, float fac) { return (fac*target) + (1.0f-fac)*origin; @@ -1334,9 +1371,36 @@ void NormalQuat(float *q) } } -float *vectoquat( float *vec, short axis, short upflag) +void RotationBetweenVectorsToQuat(float *q, float v1[3], float v2[3]) +{ + float axis[3]; + float angle; + + Crossf(axis, v1, v2); + + angle = NormalizedVecAngle2(v1, v2); + + AxisAngleToQuat(q, axis, angle); +} + +void AxisAngleToQuat(float *q, float *axis, float angle) +{ + float nor[3]; + float si; + + VecCopyf(nor, axis); + Normalize(nor); + + angle /= 2; + si = (float)sin(angle); + q[0] = (float)cos(angle); + q[1] = nor[0] * si; + q[2] = nor[1] * si; + q[3] = nor[2] * si; +} + +void vectoquat(float *vec, short axis, short upflag, float *q) { - static float q1[4]; float q2[4], nor[3], *fp, mat[3][3], angle, si, co, x2, y2, z2, len1; /* first rotate to axis */ @@ -1348,11 +1412,11 @@ float *vectoquat( float *vec, short axis, short upflag) x2= -vec[0] ; y2= -vec[1] ; z2= -vec[2]; } - q1[0]=1.0; - q1[1]=q1[2]=q1[3]= 0.0; + q[0]=1.0; + q[1]=q[2]=q[3]= 0.0; len1= (float)sqrt(x2*x2+y2*y2+z2*z2); - if(len1 == 0.0) return(q1); + if(len1 == 0.0) return; /* nasty! I need a good routine for this... * problem is a rotation of an Y axis to the negative Y-axis for example. @@ -1363,9 +1427,8 @@ float *vectoquat( float *vec, short axis, short upflag) nor[1]= -z2; nor[2]= y2; - if( fabs(y2)+fabs(z2)<0.0001 ) { + if(fabs(y2)+fabs(z2)<0.0001) nor[1]= 1.0; - } co= x2; } @@ -1374,9 +1437,8 @@ float *vectoquat( float *vec, short axis, short upflag) nor[1]= 0.0; nor[2]= -x2; - if( fabs(x2)+fabs(z2)<0.0001 ) { + if(fabs(x2)+fabs(z2)<0.0001) nor[2]= 1.0; - } co= y2; } @@ -1385,9 +1447,8 @@ float *vectoquat( float *vec, short axis, short upflag) nor[1]= x2; nor[2]= 0.0; - if( fabs(x2)+fabs(y2)<0.0001 ) { + if(fabs(x2)+fabs(y2)<0.0001) nor[0]= 1.0; - } co= z2; } @@ -1397,13 +1458,13 @@ float *vectoquat( float *vec, short axis, short upflag) angle= 0.5f*saacos(co); si= (float)sin(angle); - q1[0]= (float)cos(angle); - q1[1]= nor[0]*si; - q1[2]= nor[1]*si; - q1[3]= nor[2]*si; + q[0]= (float)cos(angle); + q[1]= nor[0]*si; + q[2]= nor[1]*si; + q[3]= nor[2]*si; if(axis!=upflag) { - QuatToMat3(q1, mat); + QuatToMat3(q, mat); fp= mat[2]; if(axis==0) { @@ -1426,10 +1487,8 @@ float *vectoquat( float *vec, short axis, short upflag) q2[2]= y2*si; q2[3]= z2*si; - QuatMul(q1,q2,q1); + QuatMul(q,q2,q); } - - return(q1); } void VecUpMat3old( float *vec, float mat[][3], short axis) @@ -1729,9 +1788,13 @@ void DQuatToMat4(DualQuat *dq, float mat[][4]) void DQuatAddWeighted(DualQuat *dqsum, DualQuat *dq, float weight) { + int flipped= 0; + /* make sure we interpolate quats in the right direction */ - if (QuatDot(dq->quat, dqsum->quat) < 0) - weight = -weight; + if (QuatDot(dq->quat, dqsum->quat) < 0) { + flipped= 1; + weight= -weight; + } /* interpolate rotation and translation */ dqsum->quat[0] += weight*dq->quat[0]; @@ -1748,6 +1811,9 @@ void DQuatAddWeighted(DualQuat *dqsum, DualQuat *dq, float weight) if (dq->scale_weight) { float wmat[4][4]; + if(flipped) /* we don't want negative weights for scaling */ + weight= -weight; + Mat4CpyMat4(wmat, dq->scale); Mat4MulFloat((float*)wmat, weight); Mat4AddMat4(dqsum->scale, dqsum->scale, wmat); @@ -2094,6 +2160,14 @@ void VecLerpf(float *target, float *a, float *b, float t) target[2]= s*a[2] + t*b[2]; } +void Vec2Lerpf(float *target, float *a, float *b, float t) +{ + float s = 1.0f-t; + + target[0]= s*a[0] + t*b[0]; + target[1]= s*a[1] + t*b[1]; +} + void VecMidf(float *v, float *v1, float *v2) { v[0]= 0.5f*(v1[0]+ v2[0]); @@ -2111,8 +2185,9 @@ void VecMulf(float *v1, float f) void VecOrthoBasisf(float *v, float *v1, float *v2) { - if (v[0] == 0.0f && v[1] == 0.0f) - { + float f = sqrt(v[0]*v[0] + v[1]*v[1]); + + if (f < 1e-35f) { // degenerate case v1[0] = 0.0f; v1[1] = 1.0f; v1[2] = 0.0f; if (v[2] > 0.0f) { @@ -2122,9 +2197,8 @@ void VecOrthoBasisf(float *v, float *v1, float *v2) v2[0] = -1.0f; v2[1] = v2[2] = 0.0f; } } - else - { - float f = 1.0f/sqrt(v[0]*v[0] + v[1]*v[1]); + else { + f = 1.0f/f; v1[0] = v[1]*f; v1[1] = -v[0]*f; v1[2] = 0.0f; @@ -2157,6 +2231,11 @@ int VecEqual(float *v1, float *v2) return ((v1[0]==v2[0]) && (v1[1]==v2[1]) && (v1[2]==v2[2])); } +int VecIsNull(float *v) +{ + return (v[0] == 0 && v[1] == 0 && v[2] == 0); +} + void CalcNormShort( short *v1, short *v2, short *v3, float *n) /* is also cross product */ { float n1[3],n2[3]; @@ -2256,6 +2335,20 @@ double Sqrt3d(double d) else return exp(log(d)/3); } +void NormalShortToFloat(float *out, short *in) +{ + out[0] = in[0] / 32767.0; + out[1] = in[1] / 32767.0; + out[2] = in[2] / 32767.0; +} + +void NormalFloatToShort(short *out, float *in) +{ + out[0] = (short)(in[0] * 32767.0); + out[1] = (short)(in[1] * 32767.0); + out[2] = (short)(in[2] * 32767.0); +} + /* distance v1 to line v2-v3 */ /* using Hesse formula, NO LINE PIECE! */ float DistVL2Dfl( float *v1, float *v2, float *v3) { @@ -2427,7 +2520,7 @@ short IsectLL2Df(float *v1, float *v2, float *v3, float *v4) 1: intersection */ -short IsectLLPt2Df(float x0,float y0,float x1,float y1, +static short IsectLLPt2Df(float x0,float y0,float x1,float y1, float x2,float y2,float x3,float y3, float *xi,float *yi) { @@ -2477,37 +2570,50 @@ short IsectLLPt2Df(float x0,float y0,float x1,float y1, } // 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])) -#define ISECT_EPSILON 1e-6 - /* point in tri */ int IsectPT2Df(float pt[2], float v1[2], float v2[2], float v3[2]) { - if ((SIDE_OF_LINE(v1,v2,pt)>=-ISECT_EPSILON) && - (SIDE_OF_LINE(v2,v3,pt)>=-ISECT_EPSILON) && - (SIDE_OF_LINE(v3,v1,pt)>=-ISECT_EPSILON)) - return 1; - else { - return 0; + 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 IsectPQ2Df(float pt[2], float v1[2], float v2[2], float v3[2], float v4[2]) { - if ((SIDE_OF_LINE(v1,v2,pt)>=-ISECT_EPSILON) && - (SIDE_OF_LINE(v2,v3,pt)>=-ISECT_EPSILON) && - (SIDE_OF_LINE(v3,v4,pt)>=-ISECT_EPSILON) && - (SIDE_OF_LINE(v4,v1,pt)>=-ISECT_EPSILON)) - return 1; - else - return 0; + 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; } - - /* copied from Geometry.c - todo - move to arithb.c or some other generic place we can reuse */ -#define SIDE_OF_LINE(pa,pb,pp) ((pa[0]-pp[0])*(pb[1]-pp[1]))-((pb[0]-pp[0])*(pa[1]-pp[1])) -#define POINT_IN_TRI(p0,p1,p2,p3) ((SIDE_OF_LINE(p1,p2,p0)>=0) && (SIDE_OF_LINE(p2,p3,p0)>=0) && (SIDE_OF_LINE(p3,p1,p0)>=0)) - /** * * @param min @@ -2871,6 +2977,22 @@ float VecAngle3(float *v1, float *v2, float *v3) return NormalizedVecAngle2(vec1, vec2) * 180.0/M_PI; } +float VecAngle3_2D(float *v1, float *v2, float *v3) +{ + float vec1[2], vec2[2]; + + vec1[0] = v2[0]-v1[0]; + vec1[1] = v2[1]-v1[1]; + + vec2[0] = v2[0]-v3[0]; + vec2[1] = v2[1]-v3[1]; + + Normalize2(vec1); + Normalize2(vec2); + + return NormalizedVecAngle2_2D(vec1, vec2) * 180.0/M_PI; +} + /* Return the shortest angle in degrees between the 2 vectors */ float VecAngle2(float *v1, float *v2) { @@ -2900,6 +3022,21 @@ float NormalizedVecAngle2(float *v1, float *v2) return 2.0f*saasin(VecLenf(v2, v1)/2.0); } +float NormalizedVecAngle2_2D(float *v1, float *v2) +{ + /* this is the same as acos(Inpf(v1, v2)), but more accurate */ + if (Inp2f(v1, v2) < 0.0f) { + float vec[2]; + + vec[0]= -v2[0]; + vec[1]= -v2[1]; + + return (float)M_PI - 2.0f*saasin(Vec2Lenf(vec, v1)/2.0f); + } + else + return 2.0f*saasin(Vec2Lenf(v2, v1)/2.0); +} + void euler_rot(float *beul, float ang, char axis) { float eul[3], mat1[3][3], mat2[3][3], totmat[3][3]; @@ -3351,6 +3488,80 @@ void rgb_to_hsv(float r, float g, float b, float *lh, float *ls, float *lv) *lv = v; } +/*http://brucelindbloom.com/index.html?Eqn_RGB_XYZ_Matrix.html */ + +void xyz_to_rgb(float xc, float yc, float zc, float *r, float *g, float *b, int colorspace) +{ + switch (colorspace) { + case BLI_CS_SMPTE: + *r = (3.50570 * xc) + (-1.73964 * yc) + (-0.544011 * zc); + *g = (-1.06906 * xc) + (1.97781 * yc) + (0.0351720 * zc); + *b = (0.0563117 * xc) + (-0.196994 * yc) + (1.05005 * zc); + break; + case BLI_CS_REC709: + *r = (3.240476 * xc) + (-1.537150 * yc) + (-0.498535 * zc); + *g = (-0.969256 * xc) + (1.875992 * yc) + (0.041556 * zc); + *b = (0.055648 * xc) + (-0.204043 * yc) + (1.057311 * zc); + break; + case BLI_CS_CIE: + *r = (2.28783848734076f * xc) + (-0.833367677835217f * yc) + (-0.454470795871421f * zc); + *g = (-0.511651380743862f * xc) + (1.42275837632178f * yc) + (0.0888930017552939f * zc); + *b = (0.00572040983140966f * xc) + (-0.0159068485104036f * yc) + (1.0101864083734f * zc); + break; + } +} + +/*If the requested RGB shade contains a negative weight for + one of the primaries, it lies outside the colour gamut + accessible from the given triple of primaries. Desaturate + it by adding white, equal quantities of R, G, and B, enough + to make RGB all positive. The function returns 1 if the + components were modified, zero otherwise.*/ +int constrain_rgb(float *r, float *g, float *b) +{ + float w; + + /* Amount of white needed is w = - min(0, *r, *g, *b) */ + + w = (0 < *r) ? 0 : *r; + w = (w < *g) ? w : *g; + w = (w < *b) ? w : *b; + w = -w; + + /* Add just enough white to make r, g, b all positive. */ + + if (w > 0) { + *r += w; *g += w; *b += w; + return 1; /* Colour modified to fit RGB gamut */ + } + + return 0; /* Colour within RGB gamut */ +} + +/*Transform linear RGB values to nonlinear RGB values. Rec. + 709 is ITU-R Recommendation BT. 709 (1990) ``Basic + Parameter Values for the HDTV Standard for the Studio and + for International Programme Exchange'', formerly CCIR Rec. + 709.*/ +static void gamma_correct(float *c) +{ + /* Rec. 709 gamma correction. */ + float cc = 0.018; + + if (*c < cc) { + *c *= ((1.099 * pow(cc, 0.45)) - 0.099) / cc; + } else { + *c = (1.099 * pow(*c, 0.45)) - 0.099; + } +} + +void gamma_correct_rgb(float *r, float *g, float *b) +{ + gamma_correct(r); + gamma_correct(g); + gamma_correct(b); +} + /* we define a 'cpack' here as a (3 byte color code) number that can be expressed like 0xFFAA66 or so. for that reason it is sensitive for endianness... with this function it works correctly @@ -3412,6 +3623,8 @@ void tubemap(float x, float y, float z, float *u, float *v) len= sqrt(x*x+y*y); if(len>0) { *u = (1.0 - (atan2(x/len,y/len) / M_PI)) / 2.0; + } else { + *v = *u = 0.0f; /* to avoid un-initialized variables */ } } @@ -3429,11 +3642,15 @@ void spheremap(float x, float y, float z, float *u, float *v) z/=len; *v = 1.0- saacos(z)/M_PI; + } else { + *v = *u = 0.0f; /* to avoid un-initialized variables */ } } /* ------------------------------------------------------------------------- */ +/* proposed api by ton and zr, not used yet */ +#if 0 /* ***************** m1 = m2 ***************** */ void cpy_m3_m3(float m1[][3], float m2[][3]) { @@ -3457,7 +3674,6 @@ void ident_m4(float m[][4]) m[3][0]= m[3][1]= m[3][2]= 0.0; } - /* ***************** m1 = m2 (pre) * m3 (post) ***************** */ void mul_m3_m3m3(float m1[][3], float m2[][3], float m3[][3]) { @@ -3595,6 +3811,8 @@ void mul_v3_v3m4(float *v1, float *v2, float mat[][4]) } +#endif + /* 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 @@ -3634,14 +3852,89 @@ int LineIntersectsTriangle(float p1[3], float p2[3], float v0[3], float v1[3], f 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 RayIntersectsTriangle(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; + + VecSubf(e1, v1, v0); + VecSubf(e2, v2, v0); + + Crossf(p, d, e2); + a = Inpf(e1, p); + if ((a > -0.000001) && (a < 0.000001)) return 0; + f = 1.0f/a; + + VecSubf(s, p1, v0); + + Crossf(q, s, e1); + *lambda = f * Inpf(e2, q); + if ((*lambda < 0.0)) return 0; + + u = f * Inpf(s, p); + if ((u < 0.0)||(u > 1.0)) return 0; + + v = f * Inpf(d, q); + if ((v < 0.0)||((u + v) > 1.0)) 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 = 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 SweepingSphereIntersectsTriangleUV(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, t, t0, t1, radius2=radius*radius; + float a, b, c, d, e, x, y, z, radius2=radius*radius; float elen2,edotv,edotbv,nordotv,vel2; - int embedded_in_plane=0, found_by_sweep=0; + float newLambda; + int found_by_sweep=0; VecSubf(e1,v1,v0); VecSubf(e2,v2,v0); @@ -3650,44 +3943,41 @@ int SweepingSphereIntersectsTriangleUV(float p1[3], float p2[3], float radius, f /*---test plane of tri---*/ Crossf(nor,e1,e2); Normalize(nor); + /* flip normal */ if(Inpf(nor,vel)>0.0f) VecMulf(nor,-1.0f); a=Inpf(p1,nor)-Inpf(v0,nor); - nordotv=Inpf(nor,vel); - if ((nordotv > -0.000001) && (nordotv < 0.000001)) { - if(fabs(a)>=1.0f) + if (fabs(nordotv) < 0.000001) + { + if(fabs(a)>=radius) + { return 0; - else{ - embedded_in_plane=1; - t0=0.0f; - t1=1.0f; } } - else{ - t0=(radius-a)/nordotv; - t1=(-radius-a)/nordotv; - /* make t0<t1 */ - if(t0>t1){b=t1; t1=t0; t0=b;} + 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] */ - t0=(t0<0.0f)?0.0f:((t0>1.0f)?1.0:t0); - t1=(t1<0.0f)?0.0f:((t1>1.0f)?1.0:t1); - } + CLAMP(t0, 0.0f, 1.0f); + CLAMP(t1, 0.0f, 1.0f); -/*---test inside of tri---*/ - if(embedded_in_plane==0){ + /*---test inside of tri---*/ /* plane intersection point */ - VecCopyf(point,vel); - VecMulf(point,t0); - VecAddf(point,point,p1); - VecCopyf(temp,nor); - VecMulf(temp,radius); - VecSubf(point,point,temp); + + 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=Inpf(e1,e1); @@ -3702,14 +3992,19 @@ int SweepingSphereIntersectsTriangleUV(float p1[3], float p2[3], float radius, f y=e*a-d*b; z=x+y-(a*c-b*b); - if(( ((unsigned int)z)& ~(((unsigned int)x)|((unsigned int)y)) ) & 0x80000000){ + + if( z <= 0.0f && (x >= 0.0f && y >= 0.0f)) + { + //( ((unsigned int)z)& ~(((unsigned int)x)|((unsigned int)y)) ) & 0x80000000){ *lambda=t0; VecCopyf(ipoint,point); return 1; } } + *lambda=1.0f; + /*---test points---*/ a=vel2=Inpf(vel,vel); @@ -3717,73 +4012,42 @@ int SweepingSphereIntersectsTriangleUV(float p1[3], float p2[3], float radius, f VecSubf(temp,p1,v0); b=2.0f*Inpf(vel,temp); c=Inpf(temp,temp)-radius2; - d=b*b-4*a*c; - - if(d>=0.0f){ - if(d==0.0f) - t=-b/2*a; - else{ - z=sqrt(d); - x=(-b-z)*0.5/a; - y=(-b+z)*0.5/a; - t=x<y?x:y; - } - if(t>0.0 && t < *lambda){ - *lambda=t; - VecCopyf(ipoint,v0); - found_by_sweep=1; - } + if(getLowestRoot(a, b, c, *lambda, lambda)) + { + VecCopyf(ipoint,v0); + found_by_sweep=1; } /*v1*/ VecSubf(temp,p1,v1); b=2.0f*Inpf(vel,temp); c=Inpf(temp,temp)-radius2; - d=b*b-4*a*c; - - if(d>=0.0f){ - if(d==0.0f) - t=-b/2*a; - else{ - z=sqrt(d); - x=(-b-z)*0.5/a; - y=(-b+z)*0.5/a; - t=x<y?x:y; - } - if(t>0.0 && t < *lambda){ - *lambda=t; - VecCopyf(ipoint,v1); - found_by_sweep=1; - } + if(getLowestRoot(a, b, c, *lambda, lambda)) + { + VecCopyf(ipoint,v1); + found_by_sweep=1; } + /*v2*/ VecSubf(temp,p1,v2); b=2.0f*Inpf(vel,temp); c=Inpf(temp,temp)-radius2; - d=b*b-4*a*c; - - if(d>=0.0f){ - if(d==0.0f) - t=-b/2*a; - else{ - z=sqrt(d); - x=(-b-z)*0.5/a; - y=(-b+z)*0.5/a; - t=x<y?x:y; - } - if(t>0.0 && t < *lambda){ - *lambda=t; - VecCopyf(ipoint,v2); - found_by_sweep=1; - } + if(getLowestRoot(a, b, c, *lambda, lambda)) + { + VecCopyf(ipoint,v2); + found_by_sweep=1; } /*---test edges---*/ + VecSubf(e3,v2,v1); //wasnt yet calculated + + /*e1*/ VecSubf(bv,v0,p1); + elen2 = Inpf(e1,e1); edotv = Inpf(e1,vel); edotbv = Inpf(e1,bv); @@ -3791,27 +4055,18 @@ int SweepingSphereIntersectsTriangleUV(float p1[3], float p2[3], float radius, f a=elen2*(-Inpf(vel,vel))+edotv*edotv; b=2.0f*(elen2*Inpf(vel,bv)-edotv*edotbv); c=elen2*(radius2-Inpf(bv,bv))+edotbv*edotbv; - d=b*b-4*a*c; - if(d>=0.0f){ - if(d==0.0f) - t=-b/2*a; - else{ - z=sqrt(d); - x=(-b-z)*0.5/a; - y=(-b+z)*0.5/a; - t=x<y?x:y; - } - - e=(edotv*t-edotbv)/elen2; - if((e>=0.0f) && (e<=1.0f)){ - if(t>0.0 && t < *lambda){ - *lambda=t; - VecCopyf(ipoint,e1); - VecMulf(ipoint,e); - VecAddf(ipoint,ipoint,v0); - found_by_sweep=1; - } + if(getLowestRoot(a, b, c, *lambda, &newLambda)) + { + e=(edotv*newLambda-edotbv)/elen2; + + if(e >= 0.0f && e <= 1.0f) + { + *lambda = newLambda; + VecCopyf(ipoint,e1); + VecMulf(ipoint,e); + VecAddf(ipoint,ipoint,v0); + found_by_sweep=1; } } @@ -3824,32 +4079,27 @@ int SweepingSphereIntersectsTriangleUV(float p1[3], float p2[3], float radius, f a=elen2*(-Inpf(vel,vel))+edotv*edotv; b=2.0f*(elen2*Inpf(vel,bv)-edotv*edotbv); c=elen2*(radius2-Inpf(bv,bv))+edotbv*edotbv; - d=b*b-4*a*c; - if(d>=0.0f){ - if(d==0.0f) - t=-b/2*a; - else{ - z=sqrt(d); - x=(-b-z)*0.5/a; - y=(-b+z)*0.5/a; - t=x<y?x:y; - } - e=(edotv*t-edotbv)/elen2; - - if((e>=0.0f) && (e<=1.0f)){ - if(t>0.0 && t < *lambda){ - *lambda=t; - VecCopyf(ipoint,e2); - VecMulf(ipoint,e); - VecAddf(ipoint,ipoint,v0); - found_by_sweep=1; - } + if(getLowestRoot(a, b, c, *lambda, &newLambda)) + { + e=(edotv*newLambda-edotbv)/elen2; + + if(e >= 0.0f && e <= 1.0f) + { + *lambda = newLambda; + VecCopyf(ipoint,e2); + VecMulf(ipoint,e); + VecAddf(ipoint,ipoint,v0); + found_by_sweep=1; } } /*e3*/ - VecSubf(e3,v2,v1); + VecSubf(bv,v0,p1); + elen2 = Inpf(e1,e1); + edotv = Inpf(e1,vel); + edotbv = Inpf(e1,bv); + VecSubf(bv,v1,p1); elen2 = Inpf(e3,e3); edotv = Inpf(e3,vel); @@ -3858,30 +4108,22 @@ int SweepingSphereIntersectsTriangleUV(float p1[3], float p2[3], float radius, f a=elen2*(-Inpf(vel,vel))+edotv*edotv; b=2.0f*(elen2*Inpf(vel,bv)-edotv*edotbv); c=elen2*(radius2-Inpf(bv,bv))+edotbv*edotbv; - d=b*b-4*a*c; - if(d>=0.0f){ - if(d==0.0f) - t=-b/2*a; - else{ - z=sqrt(d); - x=(-b-z)*0.5/a; - y=(-b+z)*0.5/a; - t=x<y?x:y; - } - e=(edotv*t-edotbv)/elen2; - - if((e>=0.0f) && (e<=1.0f)){ - if(t>0.0 && t < *lambda){ - *lambda=t; - VecCopyf(ipoint,e3); - VecMulf(ipoint,e); - VecAddf(ipoint,ipoint,v1); - found_by_sweep=1; - } + if(getLowestRoot(a, b, c, *lambda, &newLambda)) + { + e=(edotv*newLambda-edotbv)/elen2; + + if(e >= 0.0f && e <= 1.0f) + { + *lambda = newLambda; + VecCopyf(ipoint,e3); + VecMulf(ipoint,e); + VecAddf(ipoint,ipoint,v1); + found_by_sweep=1; } } + return found_by_sweep; } int AxialLineIntersectsTriangle(int axis, float p1[3], float p2[3], float v0[3], float v1[3], float v2[3], float *lambda) @@ -3928,6 +4170,74 @@ int AxialLineIntersectsTriangle(int axis, float p1[3], float p2[3], float v0[3], 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 LineIntersectLine(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; + + VecSubf(c, v3, v1); + VecSubf(a, v2, v1); + VecSubf(b, v4, v3); + + VecCopyf(dir1, a); + Normalize(dir1); + VecCopyf(dir2, b); + Normalize(dir2); + d = Inpf(dir1, dir2); + if (d == 1.0f || d == -1.0f) { + /* colinear */ + return 0; + } + + Crossf(ab, a, b); + d = Inpf(c, ab); + + /* test if the two lines are coplanar */ + if (d > -0.000001f && d < 0.000001f) { + Crossf(cb, c, b); + + VecMulf(a, Inpf(cb, ab) / Inpf(ab, ab)); + VecAddf(i1, v1, a); + VecCopyf(i2, i1); + + return 1; /* one intersection only */ + } + /* if not */ + else { + float n[3], t[3]; + float v3t[3], v4t[3]; + VecSubf(t, v1, v3); + + /* offset between both plane where the lines lies */ + Crossf(n, a, b); + Projf(t, t, n); + + /* for the first line, offset the second line until it is coplanar */ + VecAddf(v3t, v3, t); + VecAddf(v4t, v4, t); + + VecSubf(c, v3t, v1); + VecSubf(a, v2, v1); + VecSubf(b, v4t, v3); + + Crossf(ab, a, b); + Crossf(cb, c, b); + + VecMulf(a, Inpf(cb, ab) / Inpf(ab, ab)); + VecAddf(i1, v1, a); + + /* for the second line, just substract the offset from the first intersection point */ + VecSubf(i2, i1, t); + + return 2; /* two nearest points */ + } +} + int AabbIntersectAabb(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] && @@ -3950,7 +4260,7 @@ float lambda_cp_line_ex(float p[3], float l1[3], float l2[3], float cp[3]) } /* little sister we only need to know lambda */ -float lambda_cp_line(float p[3], float l1[3], float l2[3]) +static float lambda_cp_line(float p[3], float l1[3], float l2[3]) { float h[3],u[3]; VecSubf(u, l2, l1); @@ -4166,7 +4476,7 @@ void VecfCubicInterpol(float *x1, float *v1, float *x2, float *v2, float t, floa v[2]= 3*a[2]*t2 + 2*b[2]*t + v1[2]; } -int point_in_slice(float p[3], float v1[3], float l1[3], float l2[3]) +static int point_in_slice(float p[3], float v1[3], float l1[3], float l2[3]) { /* what is a slice ? @@ -4193,7 +4503,7 @@ but see a 'spat' which is a deformed cube with paired parallel planes needs only /*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*/ -int point_in_slice_as(float p[3],float origin[3],float normal[3]) +static int point_in_slice_as(float p[3],float origin[3],float normal[3]) { float h,rp[3]; VecSubf(rp,p,origin); @@ -4203,7 +4513,7 @@ int point_in_slice_as(float p[3],float origin[3],float normal[3]) } /*mama (knowing the squared lenght of the normal)*/ -int point_in_slice_m(float p[3],float origin[3],float normal[3],float lns) +static int point_in_slice_m(float p[3],float origin[3],float normal[3],float lns) { float h,rp[3]; VecSubf(rp,p,origin); @@ -4293,3 +4603,80 @@ void LocQuatSizeToMat4(float mat[][4], float loc[3], float quat[4], float size[3 mat[3][1] = loc[1]; mat[3][2] = loc[2]; } + +/* 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 + +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) { + VecAddf(vt->tang, vt->tang, tang); + return; + } + } + + /* if not found, append a new one */ + vt= BLI_memarena_alloc((MemArena *)arena, sizeof(VertexTangent)); + VecCopyf(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... */ + VecSubf(e1, co1, co2); + VecSubf(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; + Crossf(ct, tang, tangv); + + /* check flip */ + if ((ct[0]*n[0] + ct[1]*n[1] + ct[2]*n[2]) < 0.0f) + VecMulf(tang, -1.0f); +} + +/* used for zoom values*/ +float power_of_2(float val) { + return pow(2, ceil(log(val) / log(2))); +} |