diff options
Diffstat (limited to 'source/blender/blenlib/intern/math_geom.c')
| -rw-r--r-- | source/blender/blenlib/intern/math_geom.c | 455 |
1 files changed, 446 insertions, 9 deletions
diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c index 27f36752f80..af408e6bcc9 100644 --- a/source/blender/blenlib/intern/math_geom.c +++ b/source/blender/blenlib/intern/math_geom.c @@ -26,9 +26,12 @@ * */ +#include "MEM_guardedalloc.h" + #include "BLI_math.h" #include "BLI_memarena.h" -#include "MEM_guardedalloc.h" + +#include "BKE_utildefines.h" /********************************** Polygons *********************************/ @@ -121,10 +124,6 @@ float area_tri_v3(const float v1[3], const float v2[3], const float v3[3]) /* T 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[][3], float *normal) { float x, y, z, area, max; @@ -1763,6 +1762,27 @@ int box_clip_bounds_m4(float boundbox[2][3], float bounds[4], float winmat[4][4] return flag; } +void box_minmax_bounds_m4(float min[3], float max[3], float boundbox[2][3], float mat[4][4]) +{ + float mn[3], mx[3], vec[3]; + int a; + + copy_v3_v3(mn, min); + copy_v3_v3(mx, max); + + for(a=0; a<8; a++) { + vec[0]= (a & 1)? boundbox[0][0]: boundbox[1][0]; + vec[1]= (a & 2)? boundbox[0][1]: boundbox[1][1]; + vec[2]= (a & 4)? boundbox[0][2]: boundbox[1][2]; + + mul_m4_v3(mat, vec); + DO_MINMAX(vec, mn, mx); + } + + copy_v3_v3(min, mn); + copy_v3_v3(max, mx); +} + /********************************** Mapping **********************************/ void map_to_tube(float *u, float *v,float x, float y, float z) @@ -1794,9 +1814,7 @@ void map_to_sphere(float *u, float *v,float x, float y, float z) } } -/********************************************************/ - -/* Tangents */ +/********************************* 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 @@ -1875,7 +1893,7 @@ void tangent_from_uv(float *uv1, float *uv2, float *uv3, float *co1, float *co2, } } -/********************************************************/ +/****************************** Vector Clouds ********************************/ /* vector clouds */ /* void vcloud_estimate_transform(int list_size, float (*pos)[3], float *weight,float (*rpos)[3], float *rweight, @@ -2034,3 +2052,422 @@ void vcloud_estimate_transform(int list_size, float (*pos)[3], float *weight,flo } } +/******************************* Form Factor *********************************/ + +static void vec_add_dir(float r[3], float v1[3], float v2[3], float fac) +{ + r[0]= v1[0] + fac*(v2[0] - v1[0]); + r[1]= v1[1] + fac*(v2[1] - v1[1]); + r[2]= v1[2] + fac*(v2[2] - v1[2]); +} + +static int ff_visible_quad(float p[3], float n[3], float v0[3], float v1[3], float v2[3], float q0[3], float q1[3], float q2[3], float q3[3]) +{ + static const float epsilon = 1e-6f; + float c, sd[3]; + + c= dot_v3v3(n, p); + + /* signed distances from the vertices to the plane. */ + sd[0]= dot_v3v3(n, v0) - c; + sd[1]= dot_v3v3(n, v1) - c; + sd[2]= dot_v3v3(n, v2) - c; + + if(fabsf(sd[0]) < epsilon) sd[0] = 0.0f; + if(fabsf(sd[1]) < epsilon) sd[1] = 0.0f; + if(fabsf(sd[2]) < epsilon) sd[2] = 0.0f; + + if(sd[0] > 0) { + if(sd[1] > 0) { + if(sd[2] > 0) { + // +++ + copy_v3_v3(q0, v0); + copy_v3_v3(q1, v1); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + else if(sd[2] < 0) { + // ++- + copy_v3_v3(q0, v0); + copy_v3_v3(q1, v1); + vec_add_dir(q2, v1, v2, (sd[1]/(sd[1]-sd[2]))); + vec_add_dir(q3, v0, v2, (sd[0]/(sd[0]-sd[2]))); + } + else { + // ++0 + copy_v3_v3(q0, v0); + copy_v3_v3(q1, v1); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + } + else if(sd[1] < 0) { + if(sd[2] > 0) { + // +-+ + copy_v3_v3(q0, v0); + vec_add_dir(q1, v0, v1, (sd[0]/(sd[0]-sd[1]))); + vec_add_dir(q2, v1, v2, (sd[1]/(sd[1]-sd[2]))); + copy_v3_v3(q3, v2); + } + else if(sd[2] < 0) { + // +-- + copy_v3_v3(q0, v0); + vec_add_dir(q1, v0, v1, (sd[0]/(sd[0]-sd[1]))); + vec_add_dir(q2, v0, v2, (sd[0]/(sd[0]-sd[2]))); + copy_v3_v3(q3, q2); + } + else { + // +-0 + copy_v3_v3(q0, v0); + vec_add_dir(q1, v0, v1, (sd[0]/(sd[0]-sd[1]))); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + } + else { + if(sd[2] > 0) { + // +0+ + copy_v3_v3(q0, v0); + copy_v3_v3(q1, v1); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + else if(sd[2] < 0) { + // +0- + copy_v3_v3(q0, v0); + copy_v3_v3(q1, v1); + vec_add_dir(q2, v0, v2, (sd[0]/(sd[0]-sd[2]))); + copy_v3_v3(q3, q2); + } + else { + // +00 + copy_v3_v3(q0, v0); + copy_v3_v3(q1, v1); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + } + } + else if(sd[0] < 0) { + if(sd[1] > 0) { + if(sd[2] > 0) { + // -++ + vec_add_dir(q0, v0, v1, (sd[0]/(sd[0]-sd[1]))); + copy_v3_v3(q1, v1); + copy_v3_v3(q2, v2); + vec_add_dir(q3, v0, v2, (sd[0]/(sd[0]-sd[2]))); + } + else if(sd[2] < 0) { + // -+- + vec_add_dir(q0, v0, v1, (sd[0]/(sd[0]-sd[1]))); + copy_v3_v3(q1, v1); + vec_add_dir(q2, v1, v2, (sd[1]/(sd[1]-sd[2]))); + copy_v3_v3(q3, q2); + } + else { + // -+0 + vec_add_dir(q0, v0, v1, (sd[0]/(sd[0]-sd[1]))); + copy_v3_v3(q1, v1); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + } + else if(sd[1] < 0) { + if(sd[2] > 0) { + // --+ + vec_add_dir(q0, v0, v2, (sd[0]/(sd[0]-sd[2]))); + vec_add_dir(q1, v1, v2, (sd[1]/(sd[1]-sd[2]))); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + else if(sd[2] < 0) { + // --- + return 0; + } + else { + // --0 + return 0; + } + } + else { + if(sd[2] > 0) { + // -0+ + vec_add_dir(q0, v0, v2, (sd[0]/(sd[0]-sd[2]))); + copy_v3_v3(q1, v1); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + else if(sd[2] < 0) { + // -0- + return 0; + } + else { + // -00 + return 0; + } + } + } + else { + if(sd[1] > 0) { + if(sd[2] > 0) { + // 0++ + copy_v3_v3(q0, v0); + copy_v3_v3(q1, v1); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + else if(sd[2] < 0) { + // 0+- + copy_v3_v3(q0, v0); + copy_v3_v3(q1, v1); + vec_add_dir(q2, v1, v2, (sd[1]/(sd[1]-sd[2]))); + copy_v3_v3(q3, q2); + } + else { + // 0+0 + copy_v3_v3(q0, v0); + copy_v3_v3(q1, v1); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + } + else if(sd[1] < 0) { + if(sd[2] > 0) { + // 0-+ + copy_v3_v3(q0, v0); + vec_add_dir(q1, v1, v2, (sd[1]/(sd[1]-sd[2]))); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + else if(sd[2] < 0) { + // 0-- + return 0; + } + else { + // 0-0 + return 0; + } + } + else { + if(sd[2] > 0) { + // 00+ + copy_v3_v3(q0, v0); + copy_v3_v3(q1, v1); + copy_v3_v3(q2, v2); + copy_v3_v3(q3, q2); + } + else if(sd[2] < 0) { + // 00- + return 0; + } + else { + // 000 + return 0; + } + } + } + + return 1; +} + +/* altivec optimization, this works, but is unused */ + +#if 0 +#include <Accelerate/Accelerate.h> + +typedef union { + vFloat v; + float f[4]; +} vFloatResult; + +static vFloat vec_splat_float(float val) +{ + return (vFloat){val, val, val, val}; +} + +static float ff_quad_form_factor(float *p, float *n, float *q0, float *q1, float *q2, float *q3) +{ + vFloat vcos, rlen, vrx, vry, vrz, vsrx, vsry, vsrz, gx, gy, gz, vangle; + vUInt8 rotate = (vUInt8){4,5,6,7,8,9,10,11,12,13,14,15,0,1,2,3}; + vFloatResult vresult; + float result; + + /* compute r* */ + vrx = (vFloat){q0[0], q1[0], q2[0], q3[0]} - vec_splat_float(p[0]); + vry = (vFloat){q0[1], q1[1], q2[1], q3[1]} - vec_splat_float(p[1]); + vrz = (vFloat){q0[2], q1[2], q2[2], q3[2]} - vec_splat_float(p[2]); + + /* normalize r* */ + rlen = vec_rsqrte(vrx*vrx + vry*vry + vrz*vrz + vec_splat_float(1e-16f)); + vrx = vrx*rlen; + vry = vry*rlen; + vrz = vrz*rlen; + + /* rotate r* for cross and dot */ + vsrx= vec_perm(vrx, vrx, rotate); + vsry= vec_perm(vry, vry, rotate); + vsrz= vec_perm(vrz, vrz, rotate); + + /* cross product */ + gx = vsry*vrz - vsrz*vry; + gy = vsrz*vrx - vsrx*vrz; + gz = vsrx*vry - vsry*vrx; + + /* normalize */ + rlen = vec_rsqrte(gx*gx + gy*gy + gz*gz + vec_splat_float(1e-16f)); + gx = gx*rlen; + gy = gy*rlen; + gz = gz*rlen; + + /* angle */ + vcos = vrx*vsrx + vry*vsry + vrz*vsrz; + vcos= vec_max(vec_min(vcos, vec_splat_float(1.0f)), vec_splat_float(-1.0f)); + vangle= vacosf(vcos); + + /* dot */ + vresult.v = (vec_splat_float(n[0])*gx + + vec_splat_float(n[1])*gy + + vec_splat_float(n[2])*gz)*vangle; + + result= (vresult.f[0] + vresult.f[1] + vresult.f[2] + vresult.f[3])*(0.5f/(float)M_PI); + result= MAX2(result, 0.0f); + + return result; +} + +#endif + +/* SSE optimization, acos code doesn't work */ + +#if 0 + +#include <xmmintrin.h> + +static __m128 sse_approx_acos(__m128 x) +{ + /* needs a better approximation than taylor expansion of acos, since that + * gives big erros for near 1.0 values, sqrt(2*x)*acos(1-x) should work + * better, see http://www.tom.womack.net/projects/sse-fast-arctrig.html */ + + return _mm_set_ps1(1.0f); +} + +static float ff_quad_form_factor(float *p, float *n, float *q0, float *q1, float *q2, float *q3) +{ + float r0[3], r1[3], r2[3], r3[3], g0[3], g1[3], g2[3], g3[3]; + float a1, a2, a3, a4, dot1, dot2, dot3, dot4, result; + float fresult[4] __attribute__((aligned(16))); + __m128 qx, qy, qz, rx, ry, rz, rlen, srx, sry, srz, gx, gy, gz, glen, rcos, angle, aresult; + + /* compute r */ + qx = _mm_set_ps(q3[0], q2[0], q1[0], q0[0]); + qy = _mm_set_ps(q3[1], q2[1], q1[1], q0[1]); + qz = _mm_set_ps(q3[2], q2[2], q1[2], q0[2]); + + rx = qx - _mm_set_ps1(p[0]); + ry = qy - _mm_set_ps1(p[1]); + rz = qz - _mm_set_ps1(p[2]); + + /* normalize r */ + rlen = _mm_rsqrt_ps(rx*rx + ry*ry + rz*rz + _mm_set_ps1(1e-16f)); + rx = rx*rlen; + ry = ry*rlen; + rz = rz*rlen; + + /* cross product */ + srx = _mm_shuffle_ps(rx, rx, _MM_SHUFFLE(0,3,2,1)); + sry = _mm_shuffle_ps(ry, ry, _MM_SHUFFLE(0,3,2,1)); + srz = _mm_shuffle_ps(rz, rz, _MM_SHUFFLE(0,3,2,1)); + + gx = sry*rz - srz*ry; + gy = srz*rx - srx*rz; + gz = srx*ry - sry*rx; + + /* normalize g */ + glen = _mm_rsqrt_ps(gx*gx + gy*gy + gz*gz + _mm_set_ps1(1e-16f)); + gx = gx*glen; + gy = gy*glen; + gz = gz*glen; + + /* compute angle */ + rcos = rx*srx + ry*sry + rz*srz; + rcos= _mm_max_ps(_mm_min_ps(rcos, _mm_set_ps1(1.0f)), _mm_set_ps1(-1.0f)); + + angle = sse_approx_cos(rcos); + aresult = (_mm_set_ps1(n[0])*gx + _mm_set_ps1(n[1])*gy + _mm_set_ps1(n[2])*gz)*angle; + + /* sum together */ + result= (fresult[0] + fresult[1] + fresult[2] + fresult[3])*(0.5f/(float)M_PI); + result= MAX2(result, 0.0f); + + return result; +} + +#endif + +static void ff_normalize(float n[3]) +{ + float d; + + d= dot_v3v3(n, n); + + if(d > 1.0e-35F) { + d= 1.0f/sqrtf(d); + + n[0] *= d; + n[1] *= d; + n[2] *= d; + } +} + +static float ff_quad_form_factor(float *p, float *n, float *q0, float *q1, float *q2, float *q3) +{ + float r0[3], r1[3], r2[3], r3[3], g0[3], g1[3], g2[3], g3[3]; + float a1, a2, a3, a4, dot1, dot2, dot3, dot4, result; + + sub_v3_v3v3(r0, q0, p); + sub_v3_v3v3(r1, q1, p); + sub_v3_v3v3(r2, q2, p); + sub_v3_v3v3(r3, q3, p); + + ff_normalize(r0); + ff_normalize(r1); + ff_normalize(r2); + ff_normalize(r3); + + cross_v3_v3v3(g0, r1, r0); ff_normalize(g0); + cross_v3_v3v3(g1, r2, r1); ff_normalize(g1); + cross_v3_v3v3(g2, r3, r2); ff_normalize(g2); + cross_v3_v3v3(g3, r0, r3); ff_normalize(g3); + + a1= saacosf(dot_v3v3(r0, r1)); + a2= saacosf(dot_v3v3(r1, r2)); + a3= saacosf(dot_v3v3(r2, r3)); + a4= saacosf(dot_v3v3(r3, r0)); + + dot1= dot_v3v3(n, g0); + dot2= dot_v3v3(n, g1); + dot3= dot_v3v3(n, g2); + dot4= dot_v3v3(n, g3); + + result= (a1*dot1 + a2*dot2 + a3*dot3 + a4*dot4)*0.5f/(float)M_PI; + result= MAX2(result, 0.0f); + + return result; +} + +float form_factor_hemi_poly(float p[3], float n[3], float v1[3], float v2[3], float v3[3], float v4[3]) +{ + /* computes how much hemisphere defined by point and normal is + covered by a quad or triangle, cosine weighted */ + float q0[3], q1[3], q2[3], q3[3], contrib= 0.0f; + + if(ff_visible_quad(p, n, v1, v2, v3, q0, q1, q2, q3)) + contrib += ff_quad_form_factor(p, n, q0, q1, q2, q3); + + if(v4 && ff_visible_quad(p, n, v1, v3, v4, q0, q1, q2, q3)) + contrib += ff_quad_form_factor(p, n, q0, q1, q2, q3); + + return contrib; +} + |