diff options
Diffstat (limited to 'source/blender/blenlib')
19 files changed, 975 insertions, 1094 deletions
diff --git a/source/blender/blenlib/BLI_arithb.h b/source/blender/blenlib/BLI_arithb.h index 1502c4870be..71604758b80 100644 --- a/source/blender/blenlib/BLI_arithb.h +++ b/source/blender/blenlib/BLI_arithb.h @@ -174,7 +174,39 @@ void CalcNormShort(short *v1, short *v2, short *v3, float *n); float power_of_2(float val); /** - * @section Euler conversion routines + * @section Euler conversion routines (With Custom Order) + */ + +/* Defines for rotation orders + * WARNING: must match the ePchan_RotMode in DNA_action_types.h + * order matters - types are saved to file! + */ +typedef enum eEulerRotationOrders { + EULER_ORDER_DEFAULT = 1, /* Blender 'default' (classic) is basically XYZ */ + EULER_ORDER_XYZ = 1, /* Blender 'default' (classic) - must be as 1 to sync with PoseChannel rotmode */ + EULER_ORDER_XZY, + EULER_ORDER_YXZ, + EULER_ORDER_YZX, + EULER_ORDER_ZXY, + EULER_ORDER_ZYX, + /* NOTE: there are about 6 more entries when including duplicated entries too */ +} eEulerRotationOrders; + +void EulOToQuat(float eul[3], short order, float quat[4]); +void QuatToEulO(float quat[4], float eul[3], short order); + +void EulOToMat3(float eul[3], short order, float Mat[3][3]); +void EulOToMat4(float eul[3], short order, float Mat[4][4]); + +void Mat3ToEulO(float Mat[3][3], float eul[3], short order); +void Mat4ToEulO(float Mat[4][4], float eul[3], short order); + +void Mat3ToCompatibleEulO(float mat[3][3], float eul[3], float oldrot[3], short order); + +void eulerO_rot(float beul[3], float ang, char axis, short order); + +/** + * @section Euler conversion routines (Blender XYZ) */ void EulToMat3(float *eul, float mat[][3]); @@ -185,11 +217,14 @@ void Mat4ToEul(float tmat[][4],float *eul); void EulToQuat(float *eul, float *quat); -void compatible_eul(float *eul, float *oldrot); - void Mat3ToCompatibleEul(float mat[][3], float *eul, float *oldrot); + +void compatible_eul(float *eul, float *oldrot); +void euler_rot(float *beul, float ang, char axis); + + /** * @section Quaternion arithmetic routines */ @@ -216,6 +251,8 @@ void printquat(char *str, float q[4]); void QuatInterpol(float *result, float *quat1, float *quat2, float t); void QuatAdd(float *result, float *quat1, float *quat2, float t); +void QuatToMat3(float *q, float m[][3]); +void QuatToMat4(float *q, float m[][4]); /** * @section matrix multiplication and copying routines @@ -325,6 +362,7 @@ void printvec4f(char *str, float v[4]); void VecAddf(float *v, float *v1, float *v2); void VecSubf(float *v, float *v1, float *v2); +void VecMulVecf(float *v, float *v1, float *v2); void VecLerpf(float *target, float *a, float *b, float t); void VecMidf(float *v, float *v1, float *v2); @@ -339,17 +377,18 @@ void Vec2Copyf(float *v1, float *v2); void Vec2Lerpf(float *target, float *a, float *b, float t); void AxisAngleToQuat(float *q, float *axis, float angle); +void QuatToAxisAngle(float *q, float *axis, float *angle); void RotationBetweenVectorsToQuat(float *q, float v1[3], float v2[3]); void vectoquat(float *vec, short axis, short upflag, float *q); +void VecReflect(float *out, float *v1, float *v2); +void VecBisect3(float *v, float *v1, float *v2, float *v3); float VecAngle2(float *v1, float *v2); float VecAngle3(float *v1, float *v2, float *v3); float NormalizedVecAngle2(float *v1, float *v2); float VecAngle3_2D(float *v1, float *v2, float *v3); float NormalizedVecAngle2_2D(float *v1, float *v2); - -void euler_rot(float *beul, float ang, char axis); void NormalShortToFloat(float *out, short *in); void NormalFloatToShort(short *out, float *in); @@ -421,9 +460,6 @@ void VecStar(float mat[][3],float *vec); short EenheidsMat(float mat[][3]); -void QuatToMat3(float *q, float m[][3]); -void QuatToMat4(float *q, float m[][4]); - void Mat3ToQuat_is_ok(float wmat[][3], float *q); void i_ortho(float left, float right, float bottom, float top, float nearClip, float farClip, float matrix[][4]); @@ -434,8 +470,6 @@ void i_rotate(float angle, char axis, float mat[][4]); - - void MinMax3(float *min, float *max, float *vec); void SizeToMat3(float *size, float mat[][3]); void SizeToMat4(float *size, float mat[][4]); @@ -455,8 +489,9 @@ void Mat4ToSize(float mat[][4], float *size); void triatoquat(float *v1, float *v2, float *v3, float *quat); -void LocEulSizeToMat4(float mat[][4], float loc[3], float eul[3], float size[3]); -void LocQuatSizeToMat4(float mat[][4], float loc[3], float quat[4], float size[3]); +void LocEulSizeToMat4(float mat[4][4], float loc[3], float eul[3], float size[3]); +void LocEulOSizeToMat4(float mat[4][4], float loc[3], float eul[3], float size[3], short rotOrder); +void LocQuatSizeToMat4(float mat[4][4], float loc[3], float quat[4], float size[3]); void tubemap(float x, float y, float z, float *u, float *v); void spheremap(float x, float y, float z, float *u, float *v); diff --git a/source/blender/blenlib/BLI_fileops.h b/source/blender/blenlib/BLI_fileops.h index d2dd40b21a5..0c1cdbc4d3a 100644 --- a/source/blender/blenlib/BLI_fileops.h +++ b/source/blender/blenlib/BLI_fileops.h @@ -36,8 +36,6 @@ #ifndef BLI_FILEOPS_H #define BLI_FILEOPS_H - - void BLI_recurdir_fileops(char *dirname); int BLI_link(char *file, char *to); int BLI_is_writable(char *filename); diff --git a/source/blender/blenlib/BLI_kdopbvh.h b/source/blender/blenlib/BLI_kdopbvh.h index 50462d531ef..fe6bc576fbd 100644 --- a/source/blender/blenlib/BLI_kdopbvh.h +++ b/source/blender/blenlib/BLI_kdopbvh.h @@ -71,6 +71,8 @@ typedef void (*BVHTree_NearestPointCallback) (void *userdata, int index, const f /* callback must update hit in case it finds a nearest successful hit */ typedef void (*BVHTree_RayCastCallback) (void *userdata, int index, const BVHTreeRay *ray, BVHTreeRayHit *hit); +/* callback to range search query */ +typedef void (*BVHTree_RangeQuery) (void *userdata, int index, float squared_dist); BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis); void BLI_bvhtree_free(BVHTree *tree); @@ -95,5 +97,9 @@ int BLI_bvhtree_ray_cast(BVHTree *tree, const float *co, const float *dir, float float BLI_bvhtree_bb_raycast(float *bv, float *light_start, float *light_end, float *pos); +/* range query */ +int BLI_bvhtree_range_query(BVHTree *tree, const float *co, float radius, BVHTree_RangeQuery callback, void *userdata); + + #endif // BLI_KDOPBVH_H diff --git a/source/blender/blenlib/MTC_vectorops.h b/source/blender/blenlib/BLI_voxel.h index 4fec93b37b9..934bc820259 100644 --- a/source/blender/blenlib/MTC_vectorops.h +++ b/source/blender/blenlib/BLI_voxel.h @@ -1,7 +1,4 @@ -/* - * vectorops.h - * - * $Id$ +/** * * ***** BEGIN GPL LICENSE BLOCK ***** * @@ -24,35 +21,21 @@ * * The Original Code is: all of this file. * - * Contributor(s): none yet. + * Contributor(s): Matt Ebb, Raul Fernandez Hernandez (Farsthary). * * ***** END GPL LICENSE BLOCK ***** */ -#ifndef VECTOROPS_H -#define VECTOROPS_H - -/* ------------------------------------------------------------------------- */ - -void MTC_diff3Int(int v1[3], int v2[3], int v3[3]); -void MTC_cross3Int(int v1[3], int v2[3], int v3[3]); -int MTC_dot3Int(int v1[3], int v2[3]); - -void MTC_diff3Float(float v1[3], float v2[3], float v3[3]); -void MTC_cross3Float(float v1[3], float v2[3], float v3[3]); -float MTC_dot3Float(float v1[3], float v2[3]); -void MTC_cp3Float(float v1[3], float v2[3]); -/** - * Copy vector with a minus sign (so a = -b) - */ -void MTC_cp3FloatInv(float v1[3], float v2[3]); - -void MTC_swapInt(int *i1, int *i2); +#ifndef BLI_VOXEL_H +#define BLI_VOXEL_H -void MTC_diff3DFF(double v1[3], float v2[3], float v3[3]); -void MTC_cross3Double(double v1[3], double v2[3], double v3[3]); -float MTC_normalize3DF(float n[3]); +/* find the index number of a voxel, given x/y/z integer coords and resolution vector */ +#define V_I(x, y, z, res) ( (z)*(res)[1]*(res)[0] + (y)*(res)[0] + (x) ) -/* ------------------------------------------------------------------------- */ -#endif /* VECTOROPS_H */ +/* all input coordinates must be in bounding box 0.0 - 1.0 */ +float voxel_sample_nearest(float *data, int *res, float *co); +float voxel_sample_trilinear(float *data, int *res, float *co); +float voxel_sample_triquadratic(float *data, int *res, float *co); +float voxel_sample_tricubic(float *data, int *res, float *co, int bspline); +#endif /* BLI_VOXEL_H */ diff --git a/source/blender/blenlib/BLI_winstuff.h b/source/blender/blenlib/BLI_winstuff.h index 3e1b73e51be..b46ebebacd5 100644 --- a/source/blender/blenlib/BLI_winstuff.h +++ b/source/blender/blenlib/BLI_winstuff.h @@ -65,22 +65,7 @@ extern "C" { #endif -#ifndef M_PI -#define M_PI 3.14159265358979323846 -#endif -#ifndef M_PI_2 -#define M_PI_2 1.57079632679489661923 -#endif -#ifndef M_SQRT2 -#define M_SQRT2 1.41421356237309504880 -#endif -#ifndef M_SQRT1_2 -#define M_SQRT1_2 0.70710678118654752440 -#endif -#ifndef M_1_PI -#define M_1_PI 0.318309886183790671538 -#endif - +#define _USE_MATH_DEFINES #define MAXPATHLEN MAX_PATH #ifndef S_ISREG @@ -90,6 +75,18 @@ extern "C" { #define S_ISDIR(x) ((x&S_IFMT) == S_IFDIR) #endif +/* defines for using ISO C++ conformant names */ +#define open _open +#define close _close +#define write _write +#define read _read +#define getcwd _getcwd +#define chdir _chdir +#define strdup _strdup +#define lseek _lseek +#define getpid _getpid +#define snprintf _snprintf + #ifndef FREE_WINDOWS typedef unsigned int mode_t; #endif diff --git a/source/blender/blenlib/CMakeLists.txt b/source/blender/blenlib/CMakeLists.txt index a0bf2367b98..4ed9eb4b007 100644 --- a/source/blender/blenlib/CMakeLists.txt +++ b/source/blender/blenlib/CMakeLists.txt @@ -27,20 +27,20 @@ FILE(GLOB SRC intern/*.c) SET(INC - . ../makesdna ../blenkernel ../../../intern/guardedalloc ../include - ${FREETYPE_INCLUDE_DIRS} - ${ZLIB_INC} + . ../makesdna ../blenkernel ../../../intern/guardedalloc ../include + ${FREETYPE_INCLUDE_DIRS} + ${ZLIB_INC} ) IF(CMAKE_SYSTEM_NAME MATCHES "Linux") SET(INC - ${INC} - ${BINRELOC_INC} + ${INC} + ${BINRELOC_INC} ) ENDIF(CMAKE_SYSTEM_NAME MATCHES "Linux") IF(WIN32) - SET(INC ${INC} ${PTHREADS_INC}) + SET(INC ${INC} ${PTHREADS_INC}) ENDIF(WIN32) BLENDERLIB(bf_blenlib "${SRC}" "${INC}") diff --git a/source/blender/blenlib/MTC_matrixops.h b/source/blender/blenlib/MTC_matrixops.h deleted file mode 100644 index 2bc644be564..00000000000 --- a/source/blender/blenlib/MTC_matrixops.h +++ /dev/null @@ -1,162 +0,0 @@ -/* - * matrixops.h - * - * $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: all of this file. - * - * Contributor(s): none yet. - * - * ***** END GPL LICENSE BLOCK ***** - */ - -#ifndef MATRIXOPS_H -#define MATRIXOPS_H - -#ifdef __cplusplus -extern "C" { -#endif - -/* ------------------------------------------------------------------------- */ -/* need rewriting: */ -/** - * copy the left upp3 3 by 3 of m2 to m1 - */ -void MTC_Mat3CpyMat4(float m1[][3], float m2[][4]); - -/* ------------------------------------------------------------------------- */ -/* operations based on 4 by 4 matrices */ - -/** - * Copy m1 to m2 - */ -void MTC_Mat4CpyMat4(float m1[][4], float m2[][4]); - -/** - * Multiply all matrices after the first, leave the result in the - * first argument - */ -void MTC_Mat4MulSerie(float answ[][4], - float m1[][4], float m2[][4], float m3[][4], - float m4[][4], float m5[][4], float m6[][4], - float m7[][4], float m8[][4]); - -/** - * m1 = m2 matprod m3 - */ -void MTC_Mat4MulMat4(float m1[][4], float m2[][4], float m3[][4]); - -/** - * Do vec^t prod mat, result in vec. Ignore vec[3] (vec is a - * float[3]) - */ -void MTC_Mat4MulVecfl(float mat[][4], float *vec); - -/** - * Invert mat, result in inverse. Always returns 1 - */ -int MTC_Mat4Invert(float inverse[][4], float mat[][4]); - -/** - * Make the set of mat orthonormal (mat should already be orthogonal)? - * (doesn't appear to normalize properly?) - */ -void MTC_Mat4Ortho(float mat[][4]); - -/** - * vec = mat prod vec, result in vec, ignore fourth component entirely - * (4th component is _not_ accessed!!! vec is 3d) - */ -void MTC_Mat4Mul3Vecfl(float mat[][4], float *vec); - -/** - * vec = mat prod vec, result in vec - */ -void MTC_Mat4MulVec4fl(float mat[][4], float *vec); - -/** - * Set <m> to the 4-D unity matrix - */ -void MTC_Mat4One(float m[][4]); - -/** - * Swap matrices m1 and m2 - */ -void MTC_Mat4SwapMat4(float m1[][4], float m2[][4]); - -/** - * Copy m2 to the top-left 3x3 of m1, don't touch the remaining elements. - */ -void MTC_Mat4CpyMat3nc(float m1[][4], float m2[][3]); - -/** - * m1 = m2 * m3, but only the top-left 3x3 - */ -void MTC_Mat4MulMat33(float m1[][3], float m2[][4], float m3[][3]); - -/* ------------------------------------------------------------------------- */ -/* Operations based on 3 by 3 matrices */ -/** - * Do vec^t prod mat, result in vec.(vex is 3d) - */ -void MTC_Mat3MulVecfl(float mat[][3], float *vec); - -/** - * Copy m1 to m2 - */ -void MTC_Mat3CpyMat3(float m1[][3], float m2[][3]); - -/** - * m1 = m2 prod m3 - */ -void MTC_Mat3MulMat3(float m1[][3], float m3[][3], float m2[][3]); - -/** - * vec = vec prod mat - */ -void MTC_Mat3MulVecd(float mat[][3], double *vec); - -/** - * Guess: invert matrix - * result goes to m1 - */ -void MTC_Mat3Inv(float m1[][3], float m2[][3]); - -/** - * Sort of a determinant matrix? Doesn't seem very adjoint to me... - * result goes to m1 - */ -void MTC_Mat3Adj(float m1[][3], float m[][3]); - -/** - * Set <m> to the 3D unity matrix - */ -void MTC_Mat3One(float m[][3]); - -/* ------------------------------------------------------------------------- */ - -#ifdef __cplusplus -} -#endif - -#endif /* MATRIXOPS_H */ - diff --git a/source/blender/blenlib/SConscript b/source/blender/blenlib/SConscript index 3d7d6b63e64..fc586de5085 100644 --- a/source/blender/blenlib/SConscript +++ b/source/blender/blenlib/SConscript @@ -16,4 +16,4 @@ if env['OURPLATFORM'] == 'linux2': if env['OURPLATFORM'] in ('win32-vc', 'win32-mingw', 'linuxcross', 'win64-vc'): incs += ' ' + env['BF_PTHREADS_INC'] -env.BlenderLib ( 'bf_blenlib', sources, Split(incs), Split(defs), libtype=['core'], priority = [180], compileflags =cflags ) +env.BlenderLib ( 'bf_blenlib', sources, Split(incs), Split(defs), libtype=['core','player'], priority = [363,170], compileflags =cflags ) diff --git a/source/blender/blenlib/intern/BLI_kdopbvh.c b/source/blender/blenlib/intern/BLI_kdopbvh.c index 07e81b291f5..61d9cce1a58 100644 --- a/source/blender/blenlib/intern/BLI_kdopbvh.c +++ b/source/blender/blenlib/intern/BLI_kdopbvh.c @@ -72,10 +72,10 @@ struct BVHTree char start_axis, stop_axis; // KDOP_AXES array indices according to axis }; -typedef struct BVHOverlapData -{ - BVHTree *tree1, *tree2; - BVHTreeOverlap *overlap; +typedef struct BVHOverlapData +{ + BVHTree *tree1, *tree2; + BVHTreeOverlap *overlap; int i, max_overlap; /* i is number of overlaps */ int start_axis, stop_axis; } BVHOverlapData; @@ -109,7 +109,7 @@ typedef struct BVHRayCastData //////////////////////////////////////////////////////////////////////// // Bounding Volume Hierarchy Definition -// +// // Notes: From OBB until 26-DOP --> all bounding volumes possible, just choose type below // Notes: You have to choose the type at compile time ITM // Notes: You can choose the tree type --> binary, quad, octree, choose below @@ -188,10 +188,10 @@ int ADJUST_MEMORY(void *local_memblock, void **memblock, int new_size, int *max_ ////////////////////////////////////////////////////////////////////////////////////////////////////// -// Introsort +// Introsort // with permission deriven from the following Java code: // http://ralphunden.net/content/tutorials/a-guide-to-introsort/ -// and he derived it from the SUN STL +// and he derived it from the SUN STL ////////////////////////////////////////////////////////////////////////////////////////////////////// static int size_threshold = 16; /* @@ -362,7 +362,7 @@ static void create_kdop_hull(BVHTree *tree, BVHNode *node, float *co, int numpoi float newminmax; float *bv = node->bv; int i, k; - + // don't init boudings for the moving case if(!moving) { @@ -372,7 +372,7 @@ static void create_kdop_hull(BVHTree *tree, BVHNode *node, float *co, int numpoi bv[2*i + 1] = -FLT_MAX; } } - + for(k = 0; k < numpoints; k++) { // for all Axes. @@ -394,7 +394,7 @@ static void refit_kdop_hull(BVHTree *tree, BVHNode *node, int start, int end) int i, j; float *bv = node->bv; - + for (i = tree->start_axis; i < tree->stop_axis; i++) { bv[2*i] = FLT_MAX; @@ -406,10 +406,10 @@ static void refit_kdop_hull(BVHTree *tree, BVHNode *node, int start, int end) // for all Axes. for (i = tree->start_axis; i < tree->stop_axis; i++) { - newmin = tree->nodes[j]->bv[(2 * i)]; + newmin = tree->nodes[j]->bv[(2 * i)]; if ((newmin < bv[(2 * i)])) bv[(2 * i)] = newmin; - + newmax = tree->nodes[j]->bv[(2 * i) + 1]; if ((newmax > bv[(2 * i) + 1])) bv[(2 * i) + 1] = newmax; @@ -427,14 +427,14 @@ static char get_largest_axis(float *bv) middle_point[0] = (bv[1]) - (bv[0]); // x axis middle_point[1] = (bv[3]) - (bv[2]); // y axis middle_point[2] = (bv[5]) - (bv[4]); // z axis - if (middle_point[0] > middle_point[1]) + if (middle_point[0] > middle_point[1]) { if (middle_point[0] > middle_point[2]) return 1; // max x axis else return 5; // max z axis } - else + else { if (middle_point[1] > middle_point[2]) return 3; // max y axis @@ -448,24 +448,24 @@ static char get_largest_axis(float *bv) static void node_join(BVHTree *tree, BVHNode *node) { int i, j; - + for (i = tree->start_axis; i < tree->stop_axis; i++) { node->bv[2*i] = FLT_MAX; node->bv[2*i + 1] = -FLT_MAX; } - + for (i = 0; i < tree->tree_type; i++) { - if (node->children[i]) + if (node->children[i]) { for (j = tree->start_axis; j < tree->stop_axis; j++) { - // update minimum - if (node->children[i]->bv[(2 * j)] < node->bv[(2 * j)]) + // update minimum + if (node->children[i]->bv[(2 * j)] < node->bv[(2 * j)]) node->bv[(2 * j)] = node->children[i]->bv[(2 * j)]; - - // update maximum + + // update maximum if (node->children[i]->bv[(2 * j) + 1] > node->bv[(2 * j) + 1]) node->bv[(2 * j) + 1] = node->children[i]->bv[(2 * j) + 1]; } @@ -518,7 +518,7 @@ static void bvhtree_info(BVHTree *tree) static void verify_tree(BVHTree *tree) { int i, j, check = 0; - + // check the pointer list for(i = 0; i < tree->totleaf; i++) { @@ -538,7 +538,7 @@ static void verify_tree(BVHTree *tree) check = 0; } } - + // check the leaf list for(i = 0; i < tree->totleaf; i++) { @@ -558,7 +558,7 @@ static void verify_tree(BVHTree *tree) check = 0; } } - + printf("branches: %d, leafs: %d, total: %d\n", tree->totbranch, tree->totleaf, tree->totbranch + tree->totleaf); } #endif @@ -703,7 +703,7 @@ static void non_recursive_bvh_div_nodes(BVHTree *tree, BVHNode *branches_array, BVHBuildHelper data; int depth; - + // set parent from root node to NULL BVHNode *tmp = branches_array+0; tmp->parent = NULL; @@ -722,7 +722,7 @@ static void non_recursive_bvh_div_nodes(BVHTree *tree, BVHNode *branches_array, } branches_array--; //Implicit trees use 1-based indexs - + build_implicit_tree_helper(tree, &data); //Loop tree levels (log N) loops @@ -806,11 +806,11 @@ BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis) { BVHTree *tree; int numnodes, i; - + // theres not support for trees below binary-trees :P if(tree_type < 2) return NULL; - + if(tree_type > MAX_TREETYPE) return NULL; @@ -820,13 +820,13 @@ BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis) //so that tangent rays can still hit a bounding volume.. //this bug would show up when casting a ray aligned with a kdop-axis and with an edge of 2 faces epsilon = MAX2(FLT_EPSILON, epsilon); - + if(tree) { tree->epsilon = epsilon; - tree->tree_type = tree_type; + tree->tree_type = tree_type; tree->axis = axis; - + if(axis == 26) { tree->start_axis = 0; @@ -863,13 +863,13 @@ BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis) numnodes = maxsize + implicit_needed_branches(tree_type, maxsize) + tree_type; tree->nodes = (BVHNode **)MEM_callocN(sizeof(BVHNode *)*numnodes, "BVHNodes"); - + if(!tree->nodes) { MEM_freeN(tree); return NULL; } - + tree->nodebv = (float*)MEM_callocN(sizeof(float)* axis * numnodes, "BVHNodeBV"); if(!tree->nodebv) { @@ -886,7 +886,7 @@ BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis) } tree->nodearray = (BVHNode *)MEM_callocN(sizeof(BVHNode)* numnodes, "BVHNodeArray"); - + if(!tree->nodearray) { MEM_freeN(tree->nodechild); @@ -902,14 +902,14 @@ BVHTree *BLI_bvhtree_new(int maxsize, float epsilon, char tree_type, char axis) tree->nodearray[i].bv = tree->nodebv + i * axis; tree->nodearray[i].children = tree->nodechild + i * tree_type; } - + } return tree; } void BLI_bvhtree_free(BVHTree *tree) -{ +{ if(tree) { MEM_freeN(tree->nodes); @@ -946,27 +946,27 @@ int BLI_bvhtree_insert(BVHTree *tree, int index, float *co, int numpoints) { int i; BVHNode *node = NULL; - + // insert should only possible as long as tree->totbranch is 0 if(tree->totbranch > 0) return 0; - + if(tree->totleaf+1 >= MEM_allocN_len(tree->nodes)/sizeof(*(tree->nodes))) return 0; - + // TODO check if have enough nodes in array - + node = tree->nodes[tree->totleaf] = &(tree->nodearray[tree->totleaf]); tree->totleaf++; - + create_kdop_hull(tree, node, co, numpoints, 0); node->index= index; - + // inflate the bv with some epsilon for (i = tree->start_axis; i < tree->stop_axis; i++) { - node->bv[(2 * i)] -= tree->epsilon; // minimum - node->bv[(2 * i) + 1] += tree->epsilon; // maximum + node->bv[(2 * i)] -= tree->epsilon; // minimum + node->bv[(2 * i) + 1] += tree->epsilon; // maximum } return 1; @@ -978,23 +978,23 @@ int BLI_bvhtree_update_node(BVHTree *tree, int index, float *co, float *co_movin { int i; BVHNode *node= NULL; - + // check if index exists if(index > tree->totleaf) return 0; - + node = tree->nodearray + index; - + create_kdop_hull(tree, node, co, numpoints, 0); - + if(co_moving) create_kdop_hull(tree, node, co_moving, numpoints, 1); - + // inflate the bv with some epsilon for (i = tree->start_axis; i < tree->stop_axis; i++) { - node->bv[(2 * i)] -= tree->epsilon; // minimum - node->bv[(2 * i) + 1] += tree->epsilon; // maximum + node->bv[(2 * i)] -= tree->epsilon; // minimum + node->bv[(2 * i) + 1] += tree->epsilon; // maximum } return 1; @@ -1030,24 +1030,24 @@ static int tree_overlap(BVHNode *node1, BVHNode *node2, int start_axis, int stop float *bv2 = node2->bv; float *bv1_end = bv1 + (stop_axis<<1); - + bv1 += start_axis<<1; bv2 += start_axis<<1; - + // test all axis if min + max overlap for (; bv1 != bv1_end; bv1+=2, bv2+=2) { - if ((*(bv1) > *(bv2 + 1)) || (*(bv2) > *(bv1 + 1))) + if ((*(bv1) > *(bv2 + 1)) || (*(bv2) > *(bv1 + 1))) return 0; } - + return 1; } static void traverse(BVHOverlapData *data, BVHNode *node1, BVHNode *node2) { int j; - + if(tree_overlap(node1, node2, data->start_axis, data->stop_axis)) { // check if node1 is a leaf @@ -1056,17 +1056,17 @@ static void traverse(BVHOverlapData *data, BVHNode *node1, BVHNode *node2) // check if node2 is a leaf if(!node2->totnode) { - + if(node1 == node2) { return; } - + if(data->i >= data->max_overlap) - { + { // try to make alloc'ed memory bigger data->overlap = realloc(data->overlap, sizeof(BVHTreeOverlap)*data->max_overlap*2); - + if(!data->overlap) { printf("Out of Memory in traverse\n"); @@ -1074,7 +1074,7 @@ static void traverse(BVHOverlapData *data, BVHNode *node1, BVHNode *node2) } data->max_overlap *= 2; } - + // both leafs, insert overlap! data->overlap[data->i].indexA = node1->index; data->overlap[data->i].indexB = node2->index; @@ -1092,7 +1092,7 @@ static void traverse(BVHOverlapData *data, BVHNode *node1, BVHNode *node2) } else { - + for(j = 0; j < data->tree2->tree_type; j++) { if(node1->children[j]) @@ -1108,21 +1108,21 @@ BVHTreeOverlap *BLI_bvhtree_overlap(BVHTree *tree1, BVHTree *tree2, int *result) int j, total = 0; BVHTreeOverlap *overlap = NULL, *to = NULL; BVHOverlapData **data; - + // check for compatibility of both trees (can't compare 14-DOP with 18-DOP) if((tree1->axis != tree2->axis) && (tree1->axis == 14 || tree2->axis == 14) && (tree1->axis == 18 || tree2->axis == 18)) return 0; - + // fast check root nodes for collision before doing big splitting + traversal if(!tree_overlap(tree1->nodes[tree1->totleaf], tree2->nodes[tree2->totleaf], MIN2(tree1->start_axis, tree2->start_axis), MIN2(tree1->stop_axis, tree2->stop_axis))) return 0; data = MEM_callocN(sizeof(BVHOverlapData *)* tree1->tree_type, "BVHOverlapData_star"); - + for(j = 0; j < tree1->tree_type; j++) { data[j] = (BVHOverlapData *)MEM_callocN(sizeof(BVHOverlapData), "BVHOverlapData"); - + // init BVHOverlapData data[j]->overlap = (BVHTreeOverlap *)malloc(sizeof(BVHTreeOverlap)*MAX2(tree1->totleaf, tree2->totleaf)); data[j]->tree1 = tree1; @@ -1138,25 +1138,25 @@ BVHTreeOverlap *BLI_bvhtree_overlap(BVHTree *tree1, BVHTree *tree2, int *result) { traverse(data[j], tree1->nodes[tree1->totleaf]->children[j], tree2->nodes[tree2->totleaf]); } - + for(j = 0; j < tree1->tree_type; j++) total += data[j]->i; - + to = overlap = (BVHTreeOverlap *)MEM_callocN(sizeof(BVHTreeOverlap)*total, "BVHTreeOverlap"); - + for(j = 0; j < tree1->tree_type; j++) { memcpy(to, data[j]->overlap, data[j]->i*sizeof(BVHTreeOverlap)); to+=data[j]->i; } - + for(j = 0; j < tree1->tree_type; j++) { free(data[j]->overlap); MEM_freeN(data[j]); } MEM_freeN(data); - + (*result) = total; return overlap; } @@ -1173,7 +1173,7 @@ static float squared_dist(const float *a, const float *b) } //Determines the nearest point of the given node BV. Returns the squared distance to that point. -static float calc_nearest_point(BVHNearestData *data, BVHNode *node, float *nearest) +static float calc_nearest_point(const float *proj, BVHNode *node, float *nearest) { int i; const float *bv = node->bv; @@ -1181,12 +1181,12 @@ static float calc_nearest_point(BVHNearestData *data, BVHNode *node, float *near //nearest on AABB hull for(i=0; i != 3; i++, bv += 2) { - if(bv[0] > data->proj[i]) + if(bv[0] > proj[i]) nearest[i] = bv[0]; - else if(bv[1] < data->proj[i]) + else if(bv[1] < proj[i]) nearest[i] = bv[1]; else - nearest[i] = data->proj[i]; + nearest[i] = proj[i]; } /* @@ -1208,7 +1208,7 @@ static float calc_nearest_point(BVHNearestData *data, BVHNode *node, float *near } } */ - return squared_dist(data->co, nearest); + return squared_dist(proj, nearest); } @@ -1231,7 +1231,7 @@ static void dfs_find_nearest_dfs(BVHNearestData *data, BVHNode *node) else { data->nearest.index = node->index; - data->nearest.dist = calc_nearest_point(data, node, data->nearest.co); + data->nearest.dist = calc_nearest_point(data->proj, node, data->nearest.co); } } else @@ -1240,12 +1240,12 @@ static void dfs_find_nearest_dfs(BVHNearestData *data, BVHNode *node) int i; float nearest[3]; - if(data->proj[ (int)node->main_axis ] <= node->children[0]->bv[(int)node->main_axis*2+1]) + if(data->proj[ node->main_axis ] <= node->children[0]->bv[node->main_axis*2+1]) { for(i=0; i != node->totnode; i++) { - if( calc_nearest_point(data, node->children[i], nearest) >= data->nearest.dist) continue; + if( calc_nearest_point(data->proj, node->children[i], nearest) >= data->nearest.dist) continue; dfs_find_nearest_dfs(data, node->children[i]); } } @@ -1253,7 +1253,7 @@ static void dfs_find_nearest_dfs(BVHNearestData *data, BVHNode *node) { for(i=node->totnode-1; i >= 0 ; i--) { - if( calc_nearest_point(data, node->children[i], nearest) >= data->nearest.dist) continue; + if( calc_nearest_point(data->proj, node->children[i], nearest) >= data->nearest.dist) continue; dfs_find_nearest_dfs(data, node->children[i]); } } @@ -1263,7 +1263,7 @@ static void dfs_find_nearest_dfs(BVHNearestData *data, BVHNode *node) static void dfs_find_nearest_begin(BVHNearestData *data, BVHNode *node) { float nearest[3], sdist; - sdist = calc_nearest_point(data, node, nearest); + sdist = calc_nearest_point(data->proj, node, nearest); if(sdist >= data->nearest.dist) return; dfs_find_nearest_dfs(data, node); } @@ -1301,7 +1301,7 @@ static void bfs_find_nearest(BVHNearestData *data, BVHNode *node) } current.node = node; - current.dist = calc_nearest_point(data, node, nearest); + current.dist = calc_nearest_point(data->proj, node, nearest); while(current.dist < data->nearest.dist) { @@ -1329,7 +1329,7 @@ static void bfs_find_nearest(BVHNearestData *data, BVHNode *node) } heap[heap_size].node = current.node->children[i]; - heap[heap_size].dist = calc_nearest_point(data, current.node->children[i], nearest); + heap[heap_size].dist = calc_nearest_point(data->proj, current.node->children[i], nearest); if(heap[heap_size].dist >= data->nearest.dist) continue; heap_size++; @@ -1339,7 +1339,7 @@ static void bfs_find_nearest(BVHNearestData *data, BVHNode *node) push_heaps++; } } - + if(heap_size == 0) break; current = heap[0]; @@ -1355,6 +1355,7 @@ static void bfs_find_nearest(BVHNearestData *data, BVHNode *node) } #endif + int BLI_bvhtree_find_nearest(BVHTree *tree, const float *co, BVHTreeNearest *nearest, BVHTree_NearestPointCallback callback, void *userdata) { int i; @@ -1435,7 +1436,7 @@ static float ray_nearest_hit(BVHRayCastData *data, float *bv) if(lu > low) low = lu; if(ll < upper) upper = ll; } - + if(low > upper) return FLT_MAX; } } @@ -1532,28 +1533,115 @@ float BLI_bvhtree_bb_raycast(float *bv, float *light_start, float *light_end, fl float dist = 0.0; data.hit.dist = FLT_MAX; - + // get light direction data.ray.direction[0] = light_end[0] - light_start[0]; data.ray.direction[1] = light_end[1] - light_start[1]; data.ray.direction[2] = light_end[2] - light_start[2]; - + data.ray.radius = 0.0; - + data.ray.origin[0] = light_start[0]; data.ray.origin[1] = light_start[1]; data.ray.origin[2] = light_start[2]; - + Normalize(data.ray.direction); VECCOPY(data.ray_dot_axis, data.ray.direction); - + dist = ray_nearest_hit(&data, bv); - + if(dist > 0.0) { VECADDFAC(pos, light_start, data.ray.direction, dist); } return dist; + +} + +/* + * Range Query - as request by broken :P + * + * Allocs and fills an array with the indexs of node that are on the given spherical range (center, radius) + * Returns the size of the array. + */ +typedef struct RangeQueryData +{ + BVHTree *tree; + const float *center; + float radius; //squared radius + + int hits; + + BVHTree_RangeQuery callback; + void *userdata; + + +} RangeQueryData; + + +static void dfs_range_query(RangeQueryData *data, BVHNode *node) +{ + if(node->totnode == 0) + { + + //Calculate the node min-coords (if the node was a point then this is the point coordinates) + float co[3]; + co[0] = node->bv[0]; + co[1] = node->bv[2]; + co[2] = node->bv[4]; + } + else + { + int i; + for(i=0; i != node->totnode; i++) + { + float nearest[3]; + float dist = calc_nearest_point(data->center, node->children[i], nearest); + if(dist < data->radius) + { + //Its a leaf.. call the callback + if(node->children[i]->totnode == 0) + { + data->hits++; + data->callback( data->userdata, node->children[i]->index, dist ); + } + else + dfs_range_query( data, node->children[i] ); + } + } + } } +int BLI_bvhtree_range_query(BVHTree *tree, const float *co, float radius, BVHTree_RangeQuery callback, void *userdata) +{ + BVHNode * root = tree->nodes[tree->totleaf]; + + RangeQueryData data; + data.tree = tree; + data.center = co; + data.radius = radius*radius; + data.hits = 0; + + data.callback = callback; + data.userdata = userdata; + + if(root != NULL) + { + float nearest[3]; + float dist = calc_nearest_point(data.center, root, nearest); + if(dist < data.radius) + { + //Its a leaf.. call the callback + if(root->totnode == 0) + { + data.hits++; + data.callback( data.userdata, root->index, dist ); + } + else + dfs_range_query( &data, root ); + } + } + + return data.hits; +} diff --git a/source/blender/blenlib/intern/arithb.c b/source/blender/blenlib/intern/arithb.c index a26e333e095..96056ba7783 100644 --- a/source/blender/blenlib/intern/arithb.c +++ b/source/blender/blenlib/intern/arithb.c @@ -1408,22 +1408,6 @@ void RotationBetweenVectorsToQuat(float *q, float v1[3], float v2[3]) 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) { float q2[4], nor[3], *fp, mat[3][3], angle, si, co, x2, y2, z2, len1; @@ -1610,7 +1594,7 @@ void VecUpMat3(float *vec, float mat[][3], short axis) } /* A & M Watt, Advanced animation and rendering techniques, 1992 ACM press */ -void QuatInterpolW(float *, float *, float *, float ); +void QuatInterpolW(float *, float *, float *, float ); // XXX why this? void QuatInterpolW(float *result, float *quat1, float *quat2, float t) { @@ -2176,6 +2160,13 @@ void VecSubf(float *v, float *v1, float *v2) v[2]= v1[2]- v2[2]; } +void VecMulVecf(float *v, float *v1, float *v2) +{ + v[0] = v1[0] * v2[0]; + v[1] = v1[1] * v2[1]; + v[2] = v1[2] * v2[2]; +} + void VecLerpf(float *target, float *a, float *b, float t) { float s = 1.0f-t; @@ -2800,6 +2791,241 @@ void MeanValueWeights(float v[][3], int n, float *co, float *w) /* ************ EULER *************** */ +/* Euler Rotation Order Code: + * was adapted from + ANSI C code from the article + "Euler Angle Conversion" + by Ken Shoemake, shoemake@graphics.cis.upenn.edu + in "Graphics Gems IV", Academic Press, 1994 + * for use in Blender + */ + +/* Type for rotation order info - see wiki for derivation details */ +typedef struct RotOrderInfo { + short i; /* first axis index */ + short j; /* second axis index */ + short k; /* third axis index */ + short parity; /* parity of axis permuation (even=0, odd=1) - 'n' in original code */ +} RotOrderInfo; + +/* Array of info for Rotation Order calculations + * WARNING: must be kept in same order as eEulerRotationOrders + */ +static RotOrderInfo rotOrders[]= { + /* i, j, k, n */ + {0, 1, 2, 0}, // XYZ + {0, 2, 1, 1}, // XZY + {1, 0, 2, 1}, // YXZ + {1, 2, 0, 0}, // YZX + {2, 0, 1, 0}, // ZXY + {2, 1, 0, 1} // ZYZ +}; + +/* Get relevant pointer to rotation order set from the array + * NOTE: since we start at 1 for the values, but arrays index from 0, + * there is -1 factor involved in this process... + */ +#define GET_ROTATIONORDER_INFO(order) (((order)>=1) ? &rotOrders[(order)-1] : &rotOrders[0]) + +/* Construct quaternion from Euler angles (in radians). */ +void EulOToQuat(float e[3], short order, float q[4]) +{ + RotOrderInfo *R= GET_ROTATIONORDER_INFO(order); + short i=R->i, j=R->j, k=R->k; + double ti, tj, th, ci, cj, ch, si, sj, sh, cc, cs, sc, ss; + double a[3]; + + ti = e[i]/2; tj = e[j]/2; th = e[k]/2; + + if (R->parity) e[j] = -e[j]; + + ci = cos(ti); cj = cos(tj); ch = cos(th); + si = sin(ti); sj = sin(tj); sh = sin(th); + + cc = ci*ch; cs = ci*sh; + sc = si*ch; ss = si*sh; + + a[i] = cj*sc - sj*cs; + a[j] = cj*ss + sj*cc; + a[k] = cj*cs - sj*sc; + + q[0] = cj*cc + sj*ss; + q[1] = a[0]; + q[2] = a[1]; + q[3] = a[2]; + + if (R->parity) q[j] = -q[j]; +} + +/* Convert quaternion to Euler angles (in radians). */ +void QuatToEulO(float q[4], float e[3], short order) +{ + float M[3][3]; + + QuatToMat3(q, M); + Mat3ToEulO(M, e, order); +} + +/* Construct 3x3 matrix from Euler angles (in radians). */ +void EulOToMat3(float e[3], short order, float M[3][3]) +{ + RotOrderInfo *R= GET_ROTATIONORDER_INFO(order); + short i=R->i, j=R->j, k=R->k; + double ti, tj, th, ci, cj, ch, si, sj, sh, cc, cs, sc, ss; + + if (R->parity) { + ti = -e[i]; tj = -e[j]; th = -e[k]; + } + else { + ti = e[i]; tj = e[j]; th = e[k]; + } + + ci = cos(ti); cj = cos(tj); ch = cos(th); + si = sin(ti); sj = sin(tj); sh = sin(th); + + cc = ci*ch; cs = ci*sh; + sc = si*ch; ss = si*sh; + + M[i][i] = cj*ch; M[j][i] = sj*sc-cs; M[k][i] = sj*cc+ss; + M[i][j] = cj*sh; M[j][j] = sj*ss+cc; M[k][j] = sj*cs-sc; + M[i][k] = -sj; M[j][k] = cj*si; M[k][k] = cj*ci; +} + +/* Construct 4x4 matrix from Euler angles (in radians). */ +void EulOToMat4(float e[3], short order, float M[4][4]) +{ + float m[3][3]; + + /* for now, we'll just do this the slow way (i.e. copying matrices) */ + Mat3Ortho(m); + EulOToMat3(e, order, m); + Mat4CpyMat3(M, m); +} + +/* Convert 3x3 matrix to Euler angles (in radians). */ +void Mat3ToEulO(float M[3][3], float e[3], short order) +{ + RotOrderInfo *R= GET_ROTATIONORDER_INFO(order); + short i=R->i, j=R->j, k=R->k; + double cy = sqrt(M[i][i]*M[i][i] + M[i][j]*M[i][j]); + + if (cy > 16*FLT_EPSILON) { + e[i] = atan2(M[j][k], M[k][k]); + e[j] = atan2(-M[i][k], cy); + e[k] = atan2(M[i][j], M[i][i]); + } + else { + e[i] = atan2(-M[k][j], M[j][j]); + e[j] = atan2(-M[i][k], cy); + e[k] = 0; + } + + if (R->parity) { + e[0] = -e[0]; + e[1] = -e[1]; + e[2] = -e[2]; + } +} + +/* Convert 4x4 matrix to Euler angles (in radians). */ +void Mat4ToEulO(float M[4][4], float e[3], short order) +{ + float m[3][3]; + + /* for now, we'll just do this the slow way (i.e. copying matrices) */ + Mat3CpyMat4(m, M); + Mat3Ortho(m); + Mat3ToEulO(m, e, order); +} + +/* returns two euler calculation methods, so we can pick the best */ +static void mat3_to_eulo2(float M[3][3], float *e1, float *e2, short order) +{ + RotOrderInfo *R= GET_ROTATIONORDER_INFO(order); + short i=R->i, j=R->j, k=R->k; + float m[3][3]; + double cy; + + /* process the matrix first */ + Mat3CpyMat3(m, M); + Mat3Ortho(m); + + cy= sqrt(m[i][i]*m[i][i] + m[i][j]*m[i][j]); + + if (cy > 16*FLT_EPSILON) { + e1[i] = atan2(m[j][k], m[k][k]); + e1[j] = atan2(-m[i][k], cy); + e1[k] = atan2(m[i][j], m[i][i]); + + e2[i] = atan2(-m[j][k], -m[k][k]); + e2[j] = atan2(-m[i][k], -cy); + e2[k] = atan2(-m[i][j], -m[i][i]); + } + else { + e1[i] = atan2(-m[k][j], m[j][j]); + e1[j] = atan2(-m[i][k], cy); + e1[k] = 0; + + VecCopyf(e2, e1); + } + + if (R->parity) { + e1[0] = -e1[0]; + e1[1] = -e1[1]; + e1[2] = -e1[2]; + + e2[0] = -e2[0]; + e2[1] = -e2[1]; + e2[2] = -e2[2]; + } +} + +/* uses 2 methods to retrieve eulers, and picks the closest */ +void Mat3ToCompatibleEulO(float mat[3][3], float eul[3], float oldrot[3], short order) +{ + float eul1[3], eul2[3]; + float d1, d2; + + mat3_to_eulo2(mat, eul1, eul2, order); + + compatible_eul(eul1, oldrot); + compatible_eul(eul2, oldrot); + + d1= (float)fabs(eul1[0]-oldrot[0]) + (float)fabs(eul1[1]-oldrot[1]) + (float)fabs(eul1[2]-oldrot[2]); + d2= (float)fabs(eul2[0]-oldrot[0]) + (float)fabs(eul2[1]-oldrot[1]) + (float)fabs(eul2[2]-oldrot[2]); + + /* return best, which is just the one with lowest difference */ + if (d1 > d2) + VecCopyf(eul, eul2); + else + VecCopyf(eul, eul1); +} + +/* rotate the given euler by the given angle on the specified axis */ +// NOTE: is this safe to do with different axis orders? +void eulerO_rot(float beul[3], float ang, char axis, short order) +{ + float eul[3], mat1[3][3], mat2[3][3], totmat[3][3]; + + eul[0]= eul[1]= eul[2]= 0.0f; + if (axis=='x') + eul[0]= ang; + else if (axis=='y') + eul[1]= ang; + else + eul[2]= ang; + + EulOToMat3(eul, order, mat1); + EulOToMat3(beul, order, mat2); + + Mat3MulMat3(totmat, mat2, mat1); + + Mat3ToEulO(totmat, beul, order); +} + +/* ************ EULER (old XYZ) *************** */ + +/* XYZ order */ void EulToMat3( float *eul, float mat[][3]) { double ci, cj, ch, si, sj, sh, cc, cs, sc, ss; @@ -2827,6 +3053,7 @@ void EulToMat3( float *eul, float mat[][3]) } +/* XYZ order */ void EulToMat4( float *eul,float mat[][4]) { double ci, cj, ch, si, sj, sh, cc, cs, sc, ss; @@ -2858,6 +3085,7 @@ void EulToMat4( float *eul,float mat[][4]) } /* returns two euler calculation methods, so we can pick the best */ +/* XYZ order */ static void mat3_to_eul2(float tmat[][3], float *eul1, float *eul2) { float cy, quat[4], mat[3][3]; @@ -2888,6 +3116,7 @@ static void mat3_to_eul2(float tmat[][3], float *eul1, float *eul2) } } +/* XYZ order */ void Mat3ToEul(float tmat[][3], float *eul) { float eul1[3], eul2[3]; @@ -2903,6 +3132,7 @@ void Mat3ToEul(float tmat[][3], float *eul) } } +/* XYZ order */ void Mat4ToEul(float tmat[][4], float *eul) { float tempMat[3][3]; @@ -2912,6 +3142,7 @@ void Mat4ToEul(float tmat[][4], float *eul) Mat3ToEul(tempMat, eul); } +/* XYZ order */ void QuatToEul(float *quat, float *eul) { float mat[3][3]; @@ -2920,7 +3151,7 @@ void QuatToEul(float *quat, float *eul) Mat3ToEul(mat, eul); } - +/* XYZ order */ void EulToQuat(float *eul, float *quat) { float ti, tj, th, ci, cj, ch, si, sj, sh, cc, cs, sc, ss; @@ -2936,6 +3167,155 @@ void EulToQuat(float *eul, float *quat) quat[3] = cj*cs - sj*sc; } +/* XYZ order */ +void euler_rot(float *beul, float ang, char axis) +{ + float eul[3], mat1[3][3], mat2[3][3], totmat[3][3]; + + eul[0]= eul[1]= eul[2]= 0.0f; + if(axis=='x') eul[0]= ang; + else if(axis=='y') eul[1]= ang; + else eul[2]= ang; + + EulToMat3(eul, mat1); + EulToMat3(beul, mat2); + + Mat3MulMat3(totmat, mat2, mat1); + + Mat3ToEul(totmat, beul); + +} + +/* exported to transform.c */ +/* order independent! */ +void compatible_eul(float *eul, float *oldrot) +{ + float dx, dy, dz; + + /* correct differences of about 360 degrees first */ + dx= eul[0] - oldrot[0]; + dy= eul[1] - oldrot[1]; + dz= eul[2] - oldrot[2]; + + while(fabs(dx) > 5.1) { + if(dx > 0.0f) eul[0] -= 2.0f*(float)M_PI; else eul[0]+= 2.0f*(float)M_PI; + dx= eul[0] - oldrot[0]; + } + while(fabs(dy) > 5.1) { + if(dy > 0.0f) eul[1] -= 2.0f*(float)M_PI; else eul[1]+= 2.0f*(float)M_PI; + dy= eul[1] - oldrot[1]; + } + while(fabs(dz) > 5.1) { + if(dz > 0.0f) eul[2] -= 2.0f*(float)M_PI; else eul[2]+= 2.0f*(float)M_PI; + dz= eul[2] - oldrot[2]; + } + + /* is 1 of the axis rotations larger than 180 degrees and the other small? NO ELSE IF!! */ + if( fabs(dx) > 3.2 && fabs(dy)<1.6 && fabs(dz)<1.6 ) { + if(dx > 0.0) eul[0] -= 2.0f*(float)M_PI; else eul[0]+= 2.0f*(float)M_PI; + } + if( fabs(dy) > 3.2 && fabs(dz)<1.6 && fabs(dx)<1.6 ) { + if(dy > 0.0) eul[1] -= 2.0f*(float)M_PI; else eul[1]+= 2.0f*(float)M_PI; + } + if( fabs(dz) > 3.2 && fabs(dx)<1.6 && fabs(dy)<1.6 ) { + if(dz > 0.0) eul[2] -= 2.0f*(float)M_PI; else eul[2]+= 2.0f*(float)M_PI; + } + + /* the method below was there from ancient days... but why! probably because the code sucks :) + */ +#if 0 + /* calc again */ + dx= eul[0] - oldrot[0]; + dy= eul[1] - oldrot[1]; + dz= eul[2] - oldrot[2]; + + /* special case, tested for x-z */ + + if( (fabs(dx) > 3.1 && fabs(dz) > 1.5 ) || ( fabs(dx) > 1.5 && fabs(dz) > 3.1 ) ) { + if(dx > 0.0) eul[0] -= M_PI; else eul[0]+= M_PI; + if(eul[1] > 0.0) eul[1]= M_PI - eul[1]; else eul[1]= -M_PI - eul[1]; + if(dz > 0.0) eul[2] -= M_PI; else eul[2]+= M_PI; + + } + else if( (fabs(dx) > 3.1 && fabs(dy) > 1.5 ) || ( fabs(dx) > 1.5 && fabs(dy) > 3.1 ) ) { + if(dx > 0.0) eul[0] -= M_PI; else eul[0]+= M_PI; + if(dy > 0.0) eul[1] -= M_PI; else eul[1]+= M_PI; + if(eul[2] > 0.0) eul[2]= M_PI - eul[2]; else eul[2]= -M_PI - eul[2]; + } + else if( (fabs(dy) > 3.1 && fabs(dz) > 1.5 ) || ( fabs(dy) > 1.5 && fabs(dz) > 3.1 ) ) { + if(eul[0] > 0.0) eul[0]= M_PI - eul[0]; else eul[0]= -M_PI - eul[0]; + if(dy > 0.0) eul[1] -= M_PI; else eul[1]+= M_PI; + if(dz > 0.0) eul[2] -= M_PI; else eul[2]+= M_PI; + } +#endif +} + +/* uses 2 methods to retrieve eulers, and picks the closest */ +/* XYZ order */ +void Mat3ToCompatibleEul(float mat[][3], float *eul, float *oldrot) +{ + float eul1[3], eul2[3]; + float d1, d2; + + mat3_to_eul2(mat, eul1, eul2); + + compatible_eul(eul1, oldrot); + compatible_eul(eul2, oldrot); + + d1= (float)fabs(eul1[0]-oldrot[0]) + (float)fabs(eul1[1]-oldrot[1]) + (float)fabs(eul1[2]-oldrot[2]); + d2= (float)fabs(eul2[0]-oldrot[0]) + (float)fabs(eul2[1]-oldrot[1]) + (float)fabs(eul2[2]-oldrot[2]); + + /* return best, which is just the one with lowest difference */ + if( d1 > d2) { + VecCopyf(eul, eul2); + } + else { + VecCopyf(eul, eul1); + } + +} + +/* ************ AXIS ANGLE *************** */ + +/* Axis angle to Quaternions */ +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; +} + +/* Quaternions to Axis Angle */ +void QuatToAxisAngle(float q[4], float axis[3], float *angle) +{ + float ha, si; + + /* calculate angle/2, and sin(angle/2) */ + ha= (float)acos(q[0]); + si= (float)sin(ha); + + /* from half-angle to angle */ + *angle= ha * 2; + + /* prevent division by zero for axis conversion */ + if (fabs(si) < 0.0005) + si= 1.0f; + + axis[0]= q[1] / si; + axis[1]= q[2] / si; + axis[2]= q[3] / si; +} + +/* axis angle to 3x3 matrix */ void VecRotToMat3(float *vec, float phi, float mat[][3]) { /* rotation of phi radials around vec */ @@ -2962,6 +3342,7 @@ void VecRotToMat3(float *vec, float phi, float mat[][3]) } +/* axis angle to 4x4 matrix */ void VecRotToMat4(float *vec, float phi, float mat[][4]) { float tmat[3][3]; @@ -2971,6 +3352,7 @@ void VecRotToMat4(float *vec, float phi, float mat[][4]) Mat4CpyMat3(mat, tmat); } +/* axis angle to quaternion */ void VecRotToQuat(float *vec, float phi, float *quat) { /* rotation of phi radials around vec */ @@ -2992,6 +3374,41 @@ void VecRotToQuat(float *vec, float phi, float *quat) } } +/* Returns a vector bisecting the angle at v2 formed by v1, v2 and v3 */ +void VecBisect3(float *out, float *v1, float *v2, float *v3) +{ + float d_12[3], d_23[3]; + VecSubf(d_12, v2, v1); + VecSubf(d_23, v3, v2); + Normalize(d_12); + Normalize(d_23); + VecAddf(out, d_12, d_23); + Normalize(out); +} + +/* Returns a reflection vector from a vector and a normal vector +reflect = vec - ((2 * DotVecs(vec, mirror)) * mirror) +*/ +void VecReflect(float *out, float *v1, float *v2) +{ + float vec[3], normal[3]; + float reflect[3] = {0.0f, 0.0f, 0.0f}; + float dot2; + + VecCopyf(vec, v1); + VecCopyf(normal, v2); + + Normalize(normal); + + dot2 = 2 * Inpf(vec, normal); + + reflect[0] = vec[0] - (dot2 * normal[0]); + reflect[1] = vec[1] - (dot2 * normal[1]); + reflect[2] = vec[2] - (dot2 * normal[2]); + + VecCopyf(out, reflect); +} + /* Return the angle in degrees between vecs 1-2 and 2-3 in degrees If v1 is a shoulder, v2 is the elbow and v3 is the hand, this would return the angle at the elbow */ @@ -3067,111 +3484,6 @@ float NormalizedVecAngle2_2D(float *v1, float *v2) return 2.0f*(float)saasin(Vec2Lenf(v2, v1)/2.0f); } -void euler_rot(float *beul, float ang, char axis) -{ - float eul[3], mat1[3][3], mat2[3][3], totmat[3][3]; - - eul[0]= eul[1]= eul[2]= 0.0f; - if(axis=='x') eul[0]= ang; - else if(axis=='y') eul[1]= ang; - else eul[2]= ang; - - EulToMat3(eul, mat1); - EulToMat3(beul, mat2); - - Mat3MulMat3(totmat, mat2, mat1); - - Mat3ToEul(totmat, beul); - -} - -/* exported to transform.c */ -void compatible_eul(float *eul, float *oldrot) -{ - float dx, dy, dz; - - /* correct differences of about 360 degrees first */ - dx= eul[0] - oldrot[0]; - dy= eul[1] - oldrot[1]; - dz= eul[2] - oldrot[2]; - - while(fabs(dx) > 5.1) { - if(dx > 0.0f) eul[0] -= 2.0f*(float)M_PI; else eul[0]+= 2.0f*(float)M_PI; - dx= eul[0] - oldrot[0]; - } - while(fabs(dy) > 5.1) { - if(dy > 0.0f) eul[1] -= 2.0f*(float)M_PI; else eul[1]+= 2.0f*(float)M_PI; - dy= eul[1] - oldrot[1]; - } - while(fabs(dz) > 5.1) { - if(dz > 0.0f) eul[2] -= 2.0f*(float)M_PI; else eul[2]+= 2.0f*(float)M_PI; - dz= eul[2] - oldrot[2]; - } - - /* is 1 of the axis rotations larger than 180 degrees and the other small? NO ELSE IF!! */ - if( fabs(dx) > 3.2 && fabs(dy)<1.6 && fabs(dz)<1.6 ) { - if(dx > 0.0) eul[0] -= 2.0f*(float)M_PI; else eul[0]+= 2.0f*(float)M_PI; - } - if( fabs(dy) > 3.2 && fabs(dz)<1.6 && fabs(dx)<1.6 ) { - if(dy > 0.0) eul[1] -= 2.0f*(float)M_PI; else eul[1]+= 2.0f*(float)M_PI; - } - if( fabs(dz) > 3.2 && fabs(dx)<1.6 && fabs(dy)<1.6 ) { - if(dz > 0.0) eul[2] -= 2.0f*(float)M_PI; else eul[2]+= 2.0f*(float)M_PI; - } - - /* the method below was there from ancient days... but why! probably because the code sucks :) - */ -#if 0 - /* calc again */ - dx= eul[0] - oldrot[0]; - dy= eul[1] - oldrot[1]; - dz= eul[2] - oldrot[2]; - - /* special case, tested for x-z */ - - if( (fabs(dx) > 3.1 && fabs(dz) > 1.5 ) || ( fabs(dx) > 1.5 && fabs(dz) > 3.1 ) ) { - if(dx > 0.0) eul[0] -= M_PI; else eul[0]+= M_PI; - if(eul[1] > 0.0) eul[1]= M_PI - eul[1]; else eul[1]= -M_PI - eul[1]; - if(dz > 0.0) eul[2] -= M_PI; else eul[2]+= M_PI; - - } - else if( (fabs(dx) > 3.1 && fabs(dy) > 1.5 ) || ( fabs(dx) > 1.5 && fabs(dy) > 3.1 ) ) { - if(dx > 0.0) eul[0] -= M_PI; else eul[0]+= M_PI; - if(dy > 0.0) eul[1] -= M_PI; else eul[1]+= M_PI; - if(eul[2] > 0.0) eul[2]= M_PI - eul[2]; else eul[2]= -M_PI - eul[2]; - } - else if( (fabs(dy) > 3.1 && fabs(dz) > 1.5 ) || ( fabs(dy) > 1.5 && fabs(dz) > 3.1 ) ) { - if(eul[0] > 0.0) eul[0]= M_PI - eul[0]; else eul[0]= -M_PI - eul[0]; - if(dy > 0.0) eul[1] -= M_PI; else eul[1]+= M_PI; - if(dz > 0.0) eul[2] -= M_PI; else eul[2]+= M_PI; - } -#endif -} - -/* uses 2 methods to retrieve eulers, and picks the closest */ -void Mat3ToCompatibleEul(float mat[][3], float *eul, float *oldrot) -{ - float eul1[3], eul2[3]; - float d1, d2; - - mat3_to_eul2(mat, eul1, eul2); - - compatible_eul(eul1, oldrot); - compatible_eul(eul2, oldrot); - - d1= (float)fabs(eul1[0]-oldrot[0]) + (float)fabs(eul1[1]-oldrot[1]) + (float)fabs(eul1[2]-oldrot[2]); - d2= (float)fabs(eul2[0]-oldrot[0]) + (float)fabs(eul2[1]-oldrot[1]) + (float)fabs(eul2[2]-oldrot[2]); - - /* return best, which is just the one with lowest difference */ - if( d1 > d2) { - VecCopyf(eul, eul2); - } - else { - VecCopyf(eul, eul1); - } - -} - /* ******************************************** */ void SizeToMat3( float *size, float mat[][3]) @@ -4694,7 +5006,8 @@ float PdistVL3Dfl(float *v1, float *v2, float *v3) /* make a 4x4 matrix out of 3 transform components */ /* matrices are made in the order: scale * rot * loc */ -void LocEulSizeToMat4(float mat[][4], float loc[3], float eul[3], float size[3]) +// TODO: need to have a version that allows for rotation order... +void LocEulSizeToMat4(float mat[4][4], float loc[3], float eul[3], float size[3]) { float rmat[3][3], smat[3][3], tmat[3][3]; @@ -4717,7 +5030,31 @@ void LocEulSizeToMat4(float mat[][4], float loc[3], float eul[3], float size[3]) /* make a 4x4 matrix out of 3 transform components */ /* matrices are made in the order: scale * rot * loc */ -void LocQuatSizeToMat4(float mat[][4], float loc[3], float quat[4], float size[3]) +void LocEulOSizeToMat4(float mat[4][4], float loc[3], float eul[3], float size[3], short rotOrder) +{ + float rmat[3][3], smat[3][3], tmat[3][3]; + + /* initialise new matrix */ + Mat4One(mat); + + /* make rotation + scaling part */ + EulOToMat3(eul, rotOrder, rmat); + SizeToMat3(size, smat); + Mat3MulMat3(tmat, rmat, smat); + + /* copy rot/scale part to output matrix*/ + Mat4CpyMat3(mat, tmat); + + /* copy location to matrix */ + mat[3][0] = loc[0]; + mat[3][1] = loc[1]; + mat[3][2] = loc[2]; +} + + +/* make a 4x4 matrix out of 3 transform components */ +/* matrices are made in the order: scale * rot * loc */ +void LocQuatSizeToMat4(float mat[4][4], float loc[3], float quat[4], float size[3]) { float rmat[3][3], smat[3][3], tmat[3][3]; @@ -4738,6 +5075,8 @@ void LocQuatSizeToMat4(float mat[][4], float loc[3], float quat[4], float size[3 mat[3][2] = loc[2]; } +/********************************************************/ + /* Tangents */ /* For normal map tangents we need to detect uv boundaries, and only average diff --git a/source/blender/blenlib/intern/bpath.c b/source/blender/blenlib/intern/bpath.c index 6c89afe7173..cadf8d2bdd1 100644 --- a/source/blender/blenlib/intern/bpath.c +++ b/source/blender/blenlib/intern/bpath.c @@ -26,6 +26,24 @@ * ***** END GPL LICENSE BLOCK ***** */ +#include <sys/stat.h> +#include <sys/types.h> + +#include <fcntl.h> +#include <stdio.h> +#include <string.h> +#include <stdlib.h> +#include <string.h> + +/* path/file handeling stuff */ +#ifndef WIN32 + #include <dirent.h> + #include <unistd.h> +#else + #include <io.h> + #include "BLI_winstuff.h" +#endif + #include "MEM_guardedalloc.h" #include "DNA_ID.h" /* Library */ @@ -53,23 +71,7 @@ //XXX define below from BSE_sequence.h - otherwise potentially odd behaviour #define SEQ_HAS_PATH(seq) (seq->type==SEQ_MOVIE || seq->type==SEQ_IMAGE) -/* path/file handeling stuff */ -#ifndef WIN32 - #include <dirent.h> - #include <unistd.h> -#else - #include "BLI_winstuff.h" - #include <io.h> -#endif -#include <sys/stat.h> -#include <sys/types.h> - -#include <fcntl.h> -#include <stdio.h> -#include <string.h> -#include <stdlib.h> -#include <string.h> #define FILE_MAX 240 diff --git a/source/blender/blenlib/intern/fileops.c b/source/blender/blenlib/intern/fileops.c index 0228032df01..e7dc9b0eb1f 100644 --- a/source/blender/blenlib/intern/fileops.c +++ b/source/blender/blenlib/intern/fileops.c @@ -31,27 +31,29 @@ #include <stdio.h> #include <stdlib.h> +#include <sys/types.h> +#include <sys/stat.h> +#include <fcntl.h> + +#include <errno.h> + #include "zlib.h" #ifdef WIN32 -#include "BLI_winstuff.h" #include <io.h> +#include "BLI_winstuff.h" #else #include <unistd.h> // for read close #include <sys/param.h> #endif + #include "BLI_blenlib.h" #include "BLI_storage.h" #include "BLI_fileops.h" #include "BLI_callbacks.h" -#include <sys/types.h> -#include <sys/stat.h> -#include <fcntl.h> - #include "BKE_utildefines.h" -#include <errno.h> #include "BLO_sys_types.h" // for intptr_t support diff --git a/source/blender/blenlib/intern/freetypefont.c b/source/blender/blenlib/intern/freetypefont.c index bde4b561f26..985700efda1 100644 --- a/source/blender/blenlib/intern/freetypefont.c +++ b/source/blender/blenlib/intern/freetypefont.c @@ -146,9 +146,10 @@ void freetypechar_to_vchar(FT_Face face, FT_ULong charcode, VFontData *vfd) bezt = (BezTriple*)MEM_callocN((onpoints[j])* sizeof(BezTriple),"objfnt_bezt") ; BLI_addtail(&che->nurbsbase, nu); - nu->type= CU_BEZIER+CU_2D; + nu->type= CU_BEZIER; nu->pntsu = onpoints[j]; nu->resolu= 8; + nu->flag= CU_2D; nu->flagu= CU_CYCLIC; nu->bezt = bezt; diff --git a/source/blender/blenlib/intern/matrixops.c b/source/blender/blenlib/intern/matrixops.c deleted file mode 100644 index 0f9fc65f016..00000000000 --- a/source/blender/blenlib/intern/matrixops.c +++ /dev/null @@ -1,438 +0,0 @@ -/* - * - * Some matrix operations. - * - * Always use - * - vector with x components : float x[3], int x[3], etc - * - * $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: all of this file. - * - * Contributor(s): none yet. - * - * ***** END GPL LICENSE BLOCK ***** - */ - -/* ------------------------------------------------------------------------- */ -#include <string.h> -#include "MTC_matrixops.h" -#include "MTC_vectorops.h" -/* ------------------------------------------------------------------------- */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#if defined(__sun__) || defined( __sun ) || defined (__sparc) || defined (__sparc__) -#include <strings.h> -#endif - -#define ABS(x) ((x) < 0 ? -(x) : (x)) -#define SWAP(type, a, b) { type sw_ap; sw_ap=(a); (a)=(b); (b)=sw_ap; } - -void MTC_Mat4CpyMat4(float m1[][4], float m2[][4]) -{ - memcpy(m1, m2, 4*4*sizeof(float)); -} - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat4MulSerie(float answ[][4], - float m1[][4], float m2[][4], float m3[][4], - float m4[][4], float m5[][4], float m6[][4], - float m7[][4], float m8[][4]) -{ - float temp[4][4]; - - if(m1==0 || m2==0) return; - - MTC_Mat4MulMat4(answ, m2, m1); - if(m3) { - MTC_Mat4MulMat4(temp, m3, answ); - if(m4) { - MTC_Mat4MulMat4(answ, m4, temp); - if(m5) { - MTC_Mat4MulMat4(temp, m5, answ); - if(m6) { - MTC_Mat4MulMat4(answ, m6, temp); - if(m7) { - MTC_Mat4MulMat4(temp, m7, answ); - if(m8) { - MTC_Mat4MulMat4(answ, m8, temp); - } - else MTC_Mat4CpyMat4(answ, temp); - } - } - else MTC_Mat4CpyMat4(answ, temp); - } - } - else MTC_Mat4CpyMat4(answ, temp); - } -} - -/* ------------------------------------------------------------------------- */ -void MTC_Mat4MulMat4(float m1[][4], float m2[][4], float m3[][4]) -{ - /* matrix product: c[j][k] = a[j][i].b[i][k] */ - - m1[0][0] = m2[0][0]*m3[0][0] + m2[0][1]*m3[1][0] + m2[0][2]*m3[2][0] + m2[0][3]*m3[3][0]; - m1[0][1] = m2[0][0]*m3[0][1] + m2[0][1]*m3[1][1] + m2[0][2]*m3[2][1] + m2[0][3]*m3[3][1]; - m1[0][2] = m2[0][0]*m3[0][2] + m2[0][1]*m3[1][2] + m2[0][2]*m3[2][2] + m2[0][3]*m3[3][2]; - m1[0][3] = m2[0][0]*m3[0][3] + m2[0][1]*m3[1][3] + m2[0][2]*m3[2][3] + m2[0][3]*m3[3][3]; - - m1[1][0] = m2[1][0]*m3[0][0] + m2[1][1]*m3[1][0] + m2[1][2]*m3[2][0] + m2[1][3]*m3[3][0]; - m1[1][1] = m2[1][0]*m3[0][1] + m2[1][1]*m3[1][1] + m2[1][2]*m3[2][1] + m2[1][3]*m3[3][1]; - m1[1][2] = m2[1][0]*m3[0][2] + m2[1][1]*m3[1][2] + m2[1][2]*m3[2][2] + m2[1][3]*m3[3][2]; - m1[1][3] = m2[1][0]*m3[0][3] + m2[1][1]*m3[1][3] + m2[1][2]*m3[2][3] + m2[1][3]*m3[3][3]; - - m1[2][0] = m2[2][0]*m3[0][0] + m2[2][1]*m3[1][0] + m2[2][2]*m3[2][0] + m2[2][3]*m3[3][0]; - m1[2][1] = m2[2][0]*m3[0][1] + m2[2][1]*m3[1][1] + m2[2][2]*m3[2][1] + m2[2][3]*m3[3][1]; - m1[2][2] = m2[2][0]*m3[0][2] + m2[2][1]*m3[1][2] + m2[2][2]*m3[2][2] + m2[2][3]*m3[3][2]; - m1[2][3] = m2[2][0]*m3[0][3] + m2[2][1]*m3[1][3] + m2[2][2]*m3[2][3] + m2[2][3]*m3[3][3]; - - m1[3][0] = m2[3][0]*m3[0][0] + m2[3][1]*m3[1][0] + m2[3][2]*m3[2][0] + m2[3][3]*m3[3][0]; - m1[3][1] = m2[3][0]*m3[0][1] + m2[3][1]*m3[1][1] + m2[3][2]*m3[2][1] + m2[3][3]*m3[3][1]; - m1[3][2] = m2[3][0]*m3[0][2] + m2[3][1]*m3[1][2] + m2[3][2]*m3[2][2] + m2[3][3]*m3[3][2]; - m1[3][3] = m2[3][0]*m3[0][3] + m2[3][1]*m3[1][3] + m2[3][2]*m3[2][3] + m2[3][3]*m3[3][3]; - -} -/* ------------------------------------------------------------------------- */ - -void MTC_Mat4MulVecfl(float mat[][4], float *vec) -{ - float x,y; - - x=vec[0]; - y=vec[1]; - vec[0]=x*mat[0][0] + y*mat[1][0] + mat[2][0]*vec[2] + mat[3][0]; - vec[1]=x*mat[0][1] + y*mat[1][1] + mat[2][1]*vec[2] + mat[3][1]; - vec[2]=x*mat[0][2] + y*mat[1][2] + mat[2][2]*vec[2] + mat[3][2]; -} - -/* ------------------------------------------------------------------------- */ -void MTC_Mat3MulVecfl(float mat[][3], float *vec) -{ - float x,y; - - x=vec[0]; - y=vec[1]; - vec[0]= x*mat[0][0] + y*mat[1][0] + mat[2][0]*vec[2]; - vec[1]= x*mat[0][1] + y*mat[1][1] + mat[2][1]*vec[2]; - vec[2]= x*mat[0][2] + y*mat[1][2] + mat[2][2]*vec[2]; -} - -/* ------------------------------------------------------------------------- */ - -int MTC_Mat4Invert(float inverse[][4], float mat[][4]) -{ - int i, j, k; - double temp; - float tempmat[4][4]; - float max; - int maxj; - - /* Set inverse to identity */ - for (i=0; i<4; i++) - for (j=0; j<4; j++) - inverse[i][j] = 0; - for (i=0; i<4; i++) - inverse[i][i] = 1; - - /* Copy original matrix so we don't mess it up */ - for(i = 0; i < 4; i++) - for(j = 0; j <4; j++) - tempmat[i][j] = mat[i][j]; - - for(i = 0; i < 4; i++) { - /* Look for row with max pivot */ - max = ABS(tempmat[i][i]); - maxj = i; - for(j = i + 1; j < 4; j++) { - if(ABS(tempmat[j][i]) > max) { - max = ABS(tempmat[j][i]); - maxj = j; - } - } - /* Swap rows if necessary */ - if (maxj != i) { - for( k = 0; k < 4; k++) { - SWAP(float, tempmat[i][k], tempmat[maxj][k]); - SWAP(float, inverse[i][k], inverse[maxj][k]); - } - } - - temp = tempmat[i][i]; - if (temp == 0) - return 0; /* No non-zero pivot */ - for(k = 0; k < 4; k++) { - tempmat[i][k] /= temp; - inverse[i][k] /= temp; - } - for(j = 0; j < 4; j++) { - if(j != i) { - temp = tempmat[j][i]; - for(k = 0; k < 4; k++) { - tempmat[j][k] -= tempmat[i][k]*temp; - inverse[j][k] -= inverse[i][k]*temp; - } - } - } - } - return 1; -} - -/* ------------------------------------------------------------------------- */ -void MTC_Mat3CpyMat4(float m1[][3], float m2[][4]) -{ - - m1[0][0]= m2[0][0]; - m1[0][1]= m2[0][1]; - m1[0][2]= m2[0][2]; - - m1[1][0]= m2[1][0]; - m1[1][1]= m2[1][1]; - m1[1][2]= m2[1][2]; - - m1[2][0]= m2[2][0]; - m1[2][1]= m2[2][1]; - m1[2][2]= m2[2][2]; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat3CpyMat3(float m1[][3], float m2[][3]) -{ - memcpy(m1, m2, 3*3*sizeof(float)); -} - -/* ------------------------------------------------------------------------- */ -/* void Mat3MulMat3(float m1[][3], float m3[][3], float m2[][3]) */ -void MTC_Mat3MulMat3(float m1[][3], float m3[][3], float m2[][3]) -{ - /* be careful about this rewrite... */ - /* m1[i][j] = m2[i][k]*m3[k][j], args are flipped! */ - m1[0][0]= m2[0][0]*m3[0][0] + m2[0][1]*m3[1][0] + m2[0][2]*m3[2][0]; - m1[0][1]= m2[0][0]*m3[0][1] + m2[0][1]*m3[1][1] + m2[0][2]*m3[2][1]; - m1[0][2]= m2[0][0]*m3[0][2] + m2[0][1]*m3[1][2] + m2[0][2]*m3[2][2]; - - m1[1][0]= m2[1][0]*m3[0][0] + m2[1][1]*m3[1][0] + m2[1][2]*m3[2][0]; - m1[1][1]= m2[1][0]*m3[0][1] + m2[1][1]*m3[1][1] + m2[1][2]*m3[2][1]; - m1[1][2]= m2[1][0]*m3[0][2] + m2[1][1]*m3[1][2] + m2[1][2]*m3[2][2]; - - m1[2][0]= m2[2][0]*m3[0][0] + m2[2][1]*m3[1][0] + m2[2][2]*m3[2][0]; - m1[2][1]= m2[2][0]*m3[0][1] + m2[2][1]*m3[1][1] + m2[2][2]*m3[2][1]; - m1[2][2]= m2[2][0]*m3[0][2] + m2[2][1]*m3[1][2] + m2[2][2]*m3[2][2]; - -/* m1[0]= m2[0]*m3[0] + m2[1]*m3[3] + m2[2]*m3[6]; */ -/* m1[1]= m2[0]*m3[1] + m2[1]*m3[4] + m2[2]*m3[7]; */ -/* m1[2]= m2[0]*m3[2] + m2[1]*m3[5] + m2[2]*m3[8]; */ -/* m1+=3; */ -/* m2+=3; */ -/* m1[0]= m2[0]*m3[0] + m2[1]*m3[3] + m2[2]*m3[6]; */ -/* m1[1]= m2[0]*m3[1] + m2[1]*m3[4] + m2[2]*m3[7]; */ -/* m1[2]= m2[0]*m3[2] + m2[1]*m3[5] + m2[2]*m3[8]; */ -/* m1+=3; */ -/* m2+=3; */ -/* m1[0]= m2[0]*m3[0] + m2[1]*m3[3] + m2[2]*m3[6]; */ -/* m1[1]= m2[0]*m3[1] + m2[1]*m3[4] + m2[2]*m3[7]; */ -/* m1[2]= m2[0]*m3[2] + m2[1]*m3[5] + m2[2]*m3[8]; */ -} /* end of void Mat3MulMat3(float m1[][3], float m3[][3], float m2[][3]) */ - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat4Ortho(float mat[][4]) -{ - float len; - - len= MTC_normalize3DF(mat[0]); - if(len!=0.0) mat[0][3]/= len; - len= MTC_normalize3DF(mat[1]); - if(len!=0.0) mat[1][3]/= len; - len= MTC_normalize3DF(mat[2]); - if(len!=0.0) mat[2][3]/= len; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat4Mul3Vecfl(float mat[][4], float *vec) -{ - float x,y; - /* vec = mat^T dot vec !!! or vec a row, then vec = vec dot mat*/ - - x= vec[0]; - y= vec[1]; - vec[0]= x*mat[0][0] + y*mat[1][0] + mat[2][0]*vec[2]; - vec[1]= x*mat[0][1] + y*mat[1][1] + mat[2][1]*vec[2]; - vec[2]= x*mat[0][2] + y*mat[1][2] + mat[2][2]*vec[2]; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat4One(float m[][4]) -{ - - m[0][0]= m[1][1]= m[2][2]= m[3][3]= 1.0; - m[0][1]= m[0][2]= m[0][3]= 0.0; - m[1][0]= m[1][2]= m[1][3]= 0.0; - m[2][0]= m[2][1]= m[2][3]= 0.0; - m[3][0]= m[3][1]= m[3][2]= 0.0; -} - - -/* ------------------------------------------------------------------------- */ -/* Result is a 3-vector!*/ -void MTC_Mat3MulVecd(float mat[][3], double *vec) -{ - double x,y; - - /* vec = mat^T dot vec !!! or vec a row, then vec = vec dot mat*/ - x=vec[0]; - y=vec[1]; - vec[0]= x * mat[0][0] + y * mat[1][0] + mat[2][0] * vec[2]; - vec[1]= x * mat[0][1] + y * mat[1][1] + mat[2][1] * vec[2]; - vec[2]= x * mat[0][2] + y * mat[1][2] + mat[2][2] * vec[2]; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat3Inv(float m1[][3], float m2[][3]) -{ - short a,b; - float det; - - /* first adjoint */ - MTC_Mat3Adj(m1,m2); - - /* then determinant old mat! */ - det= m2[0][0]* (m2[1][1]*m2[2][2] - m2[1][2]*m2[2][1]) - -m2[1][0]* (m2[0][1]*m2[2][2] - m2[0][2]*m2[2][1]) - +m2[2][0]* (m2[0][1]*m2[1][2] - m2[0][2]*m2[1][1]); - - if(det==0) det=1; - det= 1/det; - for(a=0;a<3;a++) { - for(b=0;b<3;b++) { - m1[a][b]*=det; - } - } -} - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat3Adj(float m1[][3], float m[][3]) -{ - m1[0][0]=m[1][1]*m[2][2]-m[1][2]*m[2][1]; - m1[0][1]= -m[0][1]*m[2][2]+m[0][2]*m[2][1]; - m1[0][2]=m[0][1]*m[1][2]-m[0][2]*m[1][1]; - - m1[1][0]= -m[1][0]*m[2][2]+m[1][2]*m[2][0]; - m1[1][1]=m[0][0]*m[2][2]-m[0][2]*m[2][0]; - m1[1][2]= -m[0][0]*m[1][2]+m[0][2]*m[1][0]; - - m1[2][0]=m[1][0]*m[2][1]-m[1][1]*m[2][0]; - m1[2][1]= -m[0][0]*m[2][1]+m[0][1]*m[2][0]; - m1[2][2]=m[0][0]*m[1][1]-m[0][1]*m[1][0]; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat3One(float m[][3]) -{ - - m[0][0]= m[1][1]= m[2][2]= 1.0; - m[0][1]= m[0][2]= 0.0; - m[1][0]= m[1][2]= 0.0; - m[2][0]= m[2][1]= 0.0; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat4SwapMat4(float m1[][4], float m2[][4]) -{ - float t; - int i, j; - - for(i = 0; i < 4; i++) { - for (j = 0; j < 4; j++) { - t = m1[i][j]; - m1[i][j] = m2[i][j]; - m2[i][j] = t; - } - } -} - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat4MulVec4fl(float mat[][4], float *vec) -{ - float x,y,z; - - x = vec[0]; - y = vec[1]; - z = vec[2]; - vec[0] = x*mat[0][0] + y*mat[1][0] + z*mat[2][0] + mat[3][0]*vec[3]; - vec[1] = x*mat[0][1] + y*mat[1][1] + z*mat[2][1] + mat[3][1]*vec[3]; - vec[2] = x*mat[0][2] + y*mat[1][2] + z*mat[2][2] + mat[3][2]*vec[3]; - vec[3] = x*mat[0][3] + y*mat[1][3] + z*mat[2][3] + mat[3][3]*vec[3]; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat4CpyMat3nc(float m1[][4], float m2[][3]) /* no clear */ -{ - m1[0][0]= m2[0][0]; - m1[0][1]= m2[0][1]; - m1[0][2]= m2[0][2]; - - m1[1][0]= m2[1][0]; - m1[1][1]= m2[1][1]; - m1[1][2]= m2[1][2]; - - m1[2][0]= m2[2][0]; - m1[2][1]= m2[2][1]; - m1[2][2]= m2[2][2]; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_Mat4MulMat33(float m1[][3], float m2[][4], float m3[][3]) -{ - /* m1_i_j = m2_i_k * m3_k_j ? */ - - m1[0][0] = m2[0][0]*m3[0][0] + m2[0][1]*m3[1][0] + m2[0][2]*m3[2][0]; - m1[0][1] = m2[0][0]*m3[0][1] + m2[0][1]*m3[1][1] + m2[0][2]*m3[2][1]; - m1[0][2] = m2[0][0]*m3[0][2] + m2[0][1]*m3[1][2] + m2[0][2]*m3[2][2]; - - m1[1][0] = m2[1][0]*m3[0][0] + m2[1][1]*m3[1][0] + m2[1][2]*m3[2][0]; - m1[1][1] = m2[1][0]*m3[0][1] + m2[1][1]*m3[1][1] + m2[1][2]*m3[2][1]; - m1[1][2] = m2[1][0]*m3[0][2] + m2[1][1]*m3[1][2] + m2[1][2]*m3[2][2]; - - m1[2][0] = m2[2][0]*m3[0][0] + m2[2][1]*m3[1][0] + m2[2][2]*m3[2][0]; - m1[2][1] = m2[2][0]*m3[0][1] + m2[2][1]*m3[1][1] + m2[2][2]*m3[2][1]; - m1[2][2] = m2[2][0]*m3[0][2] + m2[2][1]*m3[1][2] + m2[2][2]*m3[2][2]; - -} - -/* ------------------------------------------------------------------------- */ - -/* eof */ diff --git a/source/blender/blenlib/intern/storage.c b/source/blender/blenlib/intern/storage.c index e6e37c58805..cdc5cec705f 100644 --- a/source/blender/blenlib/intern/storage.c +++ b/source/blender/blenlib/intern/storage.c @@ -33,13 +33,6 @@ #include <stdio.h> #include <stdlib.h> -#ifdef WIN32 -#include "BLI_winstuff.h" -#include <sys/types.h> -#include <io.h> -#include <direct.h> -#endif - #ifndef WIN32 #include <dirent.h> #endif @@ -70,9 +63,6 @@ #include <fcntl.h> -#if !defined(WIN32) -#include <sys/mtio.h> /* tape comando's */ -#endif #include <string.h> /* strcpy etc.. */ #ifndef WIN32 @@ -85,6 +75,14 @@ #include <malloc.h> #endif +#ifdef WIN32 +#include <sys/types.h> +#include <io.h> +#include <direct.h> +#include "BLI_winstuff.h" +#endif + + /* lib includes */ #include "MEM_guardedalloc.h" diff --git a/source/blender/blenlib/intern/threads.c b/source/blender/blenlib/intern/threads.c index 2812f17d58f..ed3e07b7f7e 100644 --- a/source/blender/blenlib/intern/threads.c +++ b/source/blender/blenlib/intern/threads.c @@ -142,10 +142,10 @@ void BLI_init_threads(ListBase *threadbase, void *(*do_thread)(void *), int tot) tslot->do_thread= do_thread; tslot->avail= 1; } + + MEM_set_lock_callback(BLI_lock_malloc_thread, BLI_unlock_malloc_thread); + thread_levels++; } - - MEM_set_lock_callback(BLI_lock_malloc_thread, BLI_unlock_malloc_thread); - thread_levels++; } /* amount of available threads */ @@ -235,18 +235,21 @@ void BLI_end_threads(ListBase *threadbase) { ThreadSlot *tslot; - if (threadbase) { + /* only needed if there's actually some stuff to end + * this way we don't end up decrementing thread_levels on an empty threadbase + * */ + if (threadbase && threadbase->first != NULL) { for(tslot= threadbase->first; tslot; tslot= tslot->next) { if(tslot->avail==0) { pthread_join(tslot->pthread, NULL); } } BLI_freelistN(threadbase); + + thread_levels--; + if(thread_levels==0) + MEM_set_lock_callback(NULL, NULL); } - - thread_levels--; - if(thread_levels==0) - MEM_set_lock_callback(NULL, NULL); } void BLI_lock_thread(int type) diff --git a/source/blender/blenlib/intern/util.c b/source/blender/blenlib/intern/util.c index 3c441a81d6b..c7bb7a54457 100644 --- a/source/blender/blenlib/intern/util.c +++ b/source/blender/blenlib/intern/util.c @@ -54,6 +54,8 @@ #include "BKE_utildefines.h" + + #ifdef HAVE_CONFIG_H #include <config.h> #endif @@ -75,11 +77,6 @@ #include "BLI_winstuff.h" -/* for duplicate_defgroup */ -#if !(defined vsnprintf) -#define vsnprintf _vsnprintf -#endif - #endif diff --git a/source/blender/blenlib/intern/vectorops.c b/source/blender/blenlib/intern/vectorops.c deleted file mode 100644 index 3bff5235cfd..00000000000 --- a/source/blender/blenlib/intern/vectorops.c +++ /dev/null @@ -1,166 +0,0 @@ -/* - * - * Some vector operations. - * - * Always use - * - vector with x components : float x[3], int x[3], etc - * - * $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: all of this file. - * - * Contributor(s): none yet. - * - * ***** END GPL LICENSE BLOCK ***** - */ - -/* ------------------------------------------------------------------------- */ -/* General format: op(a, b, c): a = b op c */ -/* Copying is done cp <from, to> */ -/* ------------------------------------------------------------------------- */ - -#include "MTC_vectorops.h" -#include <math.h> - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -void MTC_diff3Int(int v1[3], int v2[3], int v3[3]) -{ - v1[0] = v2[0] - v3[0]; - v1[1] = v2[1] - v3[1]; - v1[2] = v2[2] - v3[2]; -} - -/* ------------------------------------------------------------------------- */ -void MTC_diff3Float(float v1[3], float v2[3], float v3[3]) -{ - v1[0] = v2[0] - v3[0]; - v1[1] = v2[1] - v3[1]; - v1[2] = v2[2] - v3[2]; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_cross3Int(int v1[3], int v2[3], int v3[3]) -{ - v1[0] = v2[1]*v3[2] - v2[2]*v3[1]; - v1[1] = v2[2]*v3[0] - v2[0]*v3[2]; - v1[2] = v2[0]*v3[1] - v2[1]*v3[0]; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_cross3Float(float v1[3], float v2[3], float v3[3]) -{ - v1[0] = v2[1]*v3[2] - v2[2]*v3[1]; - v1[1] = v2[2]*v3[0] - v2[0]*v3[2]; - v1[2] = v2[0]*v3[1] - v2[1]*v3[0]; -} -/* ------------------------------------------------------------------------- */ - -void MTC_cross3Double(double v1[3], double v2[3], double v3[3]) -{ - v1[0] = v2[1]*v3[2] - v2[2]*v3[1]; - v1[1] = v2[2]*v3[0] - v2[0]*v3[2]; - v1[2] = v2[0]*v3[1] - v2[1]*v3[0]; -} - -/* ------------------------------------------------------------------------- */ - -int MTC_dot3Int(int v1[3], int v2[3]) -{ - return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]); -} - -/* ------------------------------------------------------------------------- */ - -float MTC_dot3Float(float v1[3], float v2[3]) -{ - return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]); -} - -/* ------------------------------------------------------------------------- */ - -void MTC_cp3Float(float v1[3], float v2[3]) -{ - v2[0] = v1[0]; - v2[1] = v1[1]; - v2[2] = v1[2]; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_cp3FloatInv(float v1[3], float v2[3]) -{ - v2[0] = -v1[0]; - v2[1] = -v1[1]; - v2[2] = -v1[2]; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_swapInt(int *i1, int *i2) -{ - int swap; - swap = *i1; - *i1 = *i2; - *i2 = swap; -} - -/* ------------------------------------------------------------------------- */ - -void MTC_diff3DFF(double v1[3], float v2[3], float v3[3]) -{ - v1[0] = v2[0] - v3[0]; - v1[1] = v2[1] - v3[1]; - v1[2] = v2[2] - v3[2]; -} - -/* ------------------------------------------------------------------------- */ -float MTC_normalize3DF(float n[3]) -{ - float d; - - d= n[0]*n[0]+n[1]*n[1]+n[2]*n[2]; - /* FLT_EPSILON is too large! A larger value causes normalize errors in */ - /* a scaled down utah teapot */ - if(d>0.0000000000001) { - - /* d= sqrt(d); This _should_ be sqrt, but internally it's a double*/ - /* anyway. This is safe. */ - d = sqrt(d); - - n[0]/=d; - n[1]/=d; - n[2]/=d; - } else { - n[0]=n[1]=n[2]= 0.0; - d= 0.0; - } - return d; -} - -/* ------------------------------------------------------------------------- */ - -/* eof */ diff --git a/source/blender/blenlib/intern/voxel.c b/source/blender/blenlib/intern/voxel.c new file mode 100644 index 00000000000..7dad854af3a --- /dev/null +++ b/source/blender/blenlib/intern/voxel.c @@ -0,0 +1,198 @@ +/** + * + * ***** 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: all of this file. + * + * Contributor(s): Matt Ebb, Raul Fernandez Hernandez (Farsthary). + * + * ***** END GPL LICENSE BLOCK ***** + */ +#include <math.h> + +#include "BLI_voxel.h" + +#include "BKE_utildefines.h" + + +#if defined( _MSC_VER ) && !defined( __cplusplus ) +# define inline __inline +#endif // defined( _MSC_VER ) && !defined( __cplusplus ) + +static inline float D(float *data, int *res, int x, int y, int z) +{ + CLAMP(x, 0, res[0]-1); + CLAMP(y, 0, res[1]-1); + CLAMP(z, 0, res[2]-1); + return data[ V_I(x, y, z, res) ]; +} + +/* *** nearest neighbour *** */ +/* input coordinates must be in bounding box 0.0 - 1.0 */ +float voxel_sample_nearest(float *data, int *res, float *co) +{ + int xi, yi, zi; + + xi = co[0] * res[0]; + yi = co[1] * res[1]; + zi = co[2] * res[2]; + + return D(data, res, xi, yi, zi); +} + +// returns highest integer <= x as integer (slightly faster than floor()) +inline int FLOORI(float x) +{ + const int r = (int)x; + return ((x >= 0.f) || (float)r == x) ? r : (r - 1); +} + +// clamp function, cannot use the CLAMPIS macro, it sometimes returns unwanted results apparently related to gcc optimization flag -fstrict-overflow which is enabled at -O2 +// this causes the test (x + 2) < 0 with int x == 2147483647 to return false (x being an integer, x + 2 should wrap around to -2147483647 so the test < 0 should return true, which it doesn't) +inline int _clamp(int a, int b, int c) +{ + return (a < b) ? b : ((a > c) ? c : a); +} + +float voxel_sample_trilinear(float *data, int *res, float *co) +{ + if (data) { + + const float xf = co[0] * res[0] - 0.5f; + const float yf = co[1] * res[1] - 0.5f; + const float zf = co[2] * res[2] - 0.5f; + + const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf); + + const int xc[2] = {_clamp(x, 0, res[0] - 1), _clamp(x + 1, 0, res[0] - 1)}; + const int yc[2] = {res[0] * _clamp(y, 0, res[1] - 1), res[0] * _clamp(y + 1, 0, res[1] - 1)}; + const int zc[2] = {res[0] * res[1] * _clamp(z, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 1, 0, res[2] - 1)}; + + const float dx = xf - (float)x; + const float dy = yf - (float)y; + const float dz = zf - (float)z; + + const float u[2] = {1.f - dx, dx}; + const float v[2] = {1.f - dy, dy}; + const float w[2] = {1.f - dz, dz}; + + return w[0] * ( v[0] * ( u[0] * data[xc[0] + yc[0] + zc[0]] + u[1] * data[xc[1] + yc[0] + zc[0]] ) + + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[0]] + u[1] * data[xc[1] + yc[1] + zc[0]] ) ) + + w[1] * ( v[0] * ( u[0] * data[xc[0] + yc[0] + zc[1]] + u[1] * data[xc[1] + yc[0] + zc[1]] ) + + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[1]] + u[1] * data[xc[1] + yc[1] + zc[1]] ) ); + + } + return 0.f; +} + + +float voxel_sample_triquadratic(float *data, int *res, float *co) +{ + if (data) { + + const float xf = co[0] * res[0], yf = co[1] * res[1], zf = co[2] * res[2]; + const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf); + + const int xc[3] = {_clamp(x - 1, 0, res[0] - 1), _clamp(x, 0, res[0] - 1), _clamp(x + 1, 0, res[0] - 1)}; + const int yc[3] = {res[0] * _clamp(y - 1, 0, res[1] - 1), res[0] * _clamp(y, 0, res[1] - 1), res[0] * _clamp(y + 1, 0, res[1] - 1)}; + const int zc[3] = {res[0] * res[1] * _clamp(z - 1, 0, res[2] - 1), res[0] * res[1] * _clamp(z, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 1, 0, res[2] - 1)}; + + const float dx = xf - (float)x, dy = yf - (float)y, dz = zf - (float)z; + const float u[3] = {dx*(0.5f*dx - 1.f) + 0.5f, dx*(1.f - dx) + 0.5f, 0.5f*dx*dx}; + const float v[3] = {dy*(0.5f*dy - 1.f) + 0.5f, dy*(1.f - dy) + 0.5f, 0.5f*dy*dy}; + const float w[3] = {dz*(0.5f*dz - 1.f) + 0.5f, dz*(1.f - dz) + 0.5f, 0.5f*dz*dz}; + + return w[0] * ( v[0] * ( u[0] * data[xc[0] + yc[0] + zc[0]] + u[1] * data[xc[1] + yc[0] + zc[0]] + u[2] * data[xc[2] + yc[0] + zc[0]] ) + + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[0]] + u[1] * data[xc[1] + yc[1] + zc[0]] + u[2] * data[xc[2] + yc[1] + zc[0]] ) + + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[0]] + u[1] * data[xc[1] + yc[2] + zc[0]] + u[2] * data[xc[2] + yc[2] + zc[0]] ) ) + + w[1] * ( v[0] * ( u[0] * data[xc[0] + yc[0] + zc[1]] + u[1] * data[xc[1] + yc[0] + zc[1]] + u[2] * data[xc[2] + yc[0] + zc[1]] ) + + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[1]] + u[1] * data[xc[1] + yc[1] + zc[1]] + u[2] * data[xc[2] + yc[1] + zc[1]] ) + + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[1]] + u[1] * data[xc[1] + yc[2] + zc[1]] + u[2] * data[xc[2] + yc[2] + zc[1]] ) ) + + w[2] * ( v[0] * ( u[0] * data[xc[0] + yc[0] + zc[2]] + u[1] * data[xc[1] + yc[0] + zc[2]] + u[2] * data[xc[2] + yc[0] + zc[2]] ) + + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[2]] + u[1] * data[xc[1] + yc[1] + zc[2]] + u[2] * data[xc[2] + yc[1] + zc[2]] ) + + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[2]] + u[1] * data[xc[1] + yc[2] + zc[2]] + u[2] * data[xc[2] + yc[2] + zc[2]] ) ); + +} + return 0.f; +} + +float voxel_sample_tricubic(float *data, int *res, float *co, int bspline) +{ + if (data) { + + const float xf = co[0] * res[0] - 0.5f, yf = co[1] * res[1] - 0.5f, zf = co[2] * res[2] - 0.5f; + const int x = FLOORI(xf), y = FLOORI(yf), z = FLOORI(zf); + + const int xc[4] = {_clamp(x - 1, 0, res[0] - 1), _clamp(x, 0, res[0] - 1), _clamp(x + 1, 0, res[0] - 1), _clamp(x + 2, 0, res[0] - 1)}; + const int yc[4] = {res[0] * _clamp(y - 1, 0, res[1] - 1), res[0] * _clamp(y, 0, res[1] - 1), res[0] * _clamp(y + 1, 0, res[1] - 1), res[0] * _clamp(y + 2, 0, res[1] - 1)}; + const int zc[4] = {res[0] * res[1] * _clamp(z - 1, 0, res[2] - 1), res[0] * res[1] * _clamp(z, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 1, 0, res[2] - 1), res[0] * res[1] * _clamp(z + 2, 0, res[2] - 1)}; + + const float dx = xf - (float)x, dy = yf - (float)y, dz = zf - (float)z; + + float u[4], v[4], w[4]; + if (bspline) { // B-Spline + u[0] = (((-1.f/6.f)*dx + 0.5f)*dx - 0.5f)*dx + (1.f/6.f); + u[1] = (( 0.5f*dx - 1.f )*dx )*dx + (2.f/3.f); + u[2] = (( -0.5f*dx + 0.5f)*dx + 0.5f)*dx + (1.f/6.f); + u[3] = ( 1.f/6.f)*dx*dx*dx; + v[0] = (((-1.f/6.f)*dy + 0.5f)*dy - 0.5f)*dy + (1.f/6.f); + v[1] = (( 0.5f*dy - 1.f )*dy )*dy + (2.f/3.f); + v[2] = (( -0.5f*dy + 0.5f)*dy + 0.5f)*dy + (1.f/6.f); + v[3] = ( 1.f/6.f)*dy*dy*dy; + w[0] = (((-1.f/6.f)*dz + 0.5f)*dz - 0.5f)*dz + (1.f/6.f); + w[1] = (( 0.5f*dz - 1.f )*dz )*dz + (2.f/3.f); + w[2] = (( -0.5f*dz + 0.5f)*dz + 0.5f)*dz + (1.f/6.f); + w[3] = ( 1.f/6.f)*dz*dz*dz; + } + else { // Catmull-Rom + u[0] = ((-0.5f*dx + 1.0f)*dx - 0.5f)*dx; + u[1] = (( 1.5f*dx - 2.5f)*dx )*dx + 1.0f; + u[2] = ((-1.5f*dx + 2.0f)*dx + 0.5f)*dx; + u[3] = (( 0.5f*dx - 0.5f)*dx )*dx; + v[0] = ((-0.5f*dy + 1.0f)*dy - 0.5f)*dy; + v[1] = (( 1.5f*dy - 2.5f)*dy )*dy + 1.0f; + v[2] = ((-1.5f*dy + 2.0f)*dy + 0.5f)*dy; + v[3] = (( 0.5f*dy - 0.5f)*dy )*dy; + w[0] = ((-0.5f*dz + 1.0f)*dz - 0.5f)*dz; + w[1] = (( 1.5f*dz - 2.5f)*dz )*dz + 1.0f; + w[2] = ((-1.5f*dz + 2.0f)*dz + 0.5f)*dz; + w[3] = (( 0.5f*dz - 0.5f)*dz )*dz; + } + + return w[0] * ( v[0] * ( u[0] * data[xc[0] + yc[0] + zc[0]] + u[1] * data[xc[1] + yc[0] + zc[0]] + u[2] * data[xc[2] + yc[0] + zc[0]] + u[3] * data[xc[3] + yc[0] + zc[0]] ) + + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[0]] + u[1] * data[xc[1] + yc[1] + zc[0]] + u[2] * data[xc[2] + yc[1] + zc[0]] + u[3] * data[xc[3] + yc[1] + zc[0]] ) + + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[0]] + u[1] * data[xc[1] + yc[2] + zc[0]] + u[2] * data[xc[2] + yc[2] + zc[0]] + u[3] * data[xc[3] + yc[2] + zc[0]] ) + + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[0]] + u[1] * data[xc[1] + yc[3] + zc[0]] + u[2] * data[xc[2] + yc[3] + zc[0]] + u[3] * data[xc[3] + yc[3] + zc[0]] ) ) + + w[1] * ( v[0] * ( u[0] * data[xc[0] + yc[0] + zc[1]] + u[1] * data[xc[1] + yc[0] + zc[1]] + u[2] * data[xc[2] + yc[0] + zc[1]] + u[3] * data[xc[3] + yc[0] + zc[1]] ) + + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[1]] + u[1] * data[xc[1] + yc[1] + zc[1]] + u[2] * data[xc[2] + yc[1] + zc[1]] + u[3] * data[xc[3] + yc[1] + zc[1]] ) + + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[1]] + u[1] * data[xc[1] + yc[2] + zc[1]] + u[2] * data[xc[2] + yc[2] + zc[1]] + u[3] * data[xc[3] + yc[2] + zc[1]] ) + + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[1]] + u[1] * data[xc[1] + yc[3] + zc[1]] + u[2] * data[xc[2] + yc[3] + zc[1]] + u[3] * data[xc[3] + yc[3] + zc[1]] ) ) + + w[2] * ( v[0] * ( u[0] * data[xc[0] + yc[0] + zc[2]] + u[1] * data[xc[1] + yc[0] + zc[2]] + u[2] * data[xc[2] + yc[0] + zc[2]] + u[3] * data[xc[3] + yc[0] + zc[2]] ) + + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[2]] + u[1] * data[xc[1] + yc[1] + zc[2]] + u[2] * data[xc[2] + yc[1] + zc[2]] + u[3] * data[xc[3] + yc[1] + zc[2]] ) + + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[2]] + u[1] * data[xc[1] + yc[2] + zc[2]] + u[2] * data[xc[2] + yc[2] + zc[2]] + u[3] * data[xc[3] + yc[2] + zc[2]] ) + + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[2]] + u[1] * data[xc[1] + yc[3] + zc[2]] + u[2] * data[xc[2] + yc[3] + zc[2]] + u[3] * data[xc[3] + yc[3] + zc[2]] ) ) + + w[3] * ( v[0] * ( u[0] * data[xc[0] + yc[0] + zc[3]] + u[1] * data[xc[1] + yc[0] + zc[3]] + u[2] * data[xc[2] + yc[0] + zc[3]] + u[3] * data[xc[3] + yc[0] + zc[3]] ) + + v[1] * ( u[0] * data[xc[0] + yc[1] + zc[3]] + u[1] * data[xc[1] + yc[1] + zc[3]] + u[2] * data[xc[2] + yc[1] + zc[3]] + u[3] * data[xc[3] + yc[1] + zc[3]] ) + + v[2] * ( u[0] * data[xc[0] + yc[2] + zc[3]] + u[1] * data[xc[1] + yc[2] + zc[3]] + u[2] * data[xc[2] + yc[2] + zc[3]] + u[3] * data[xc[3] + yc[2] + zc[3]] ) + + v[3] * ( u[0] * data[xc[0] + yc[3] + zc[3]] + u[1] * data[xc[1] + yc[3] + zc[3]] + u[2] * data[xc[2] + yc[3] + zc[3]] + u[3] * data[xc[3] + yc[3] + zc[3]] ) ); + + } + return 0.f; +} |