diff options
Diffstat (limited to 'source/blender/blenlib')
21 files changed, 1594 insertions, 1785 deletions
diff --git a/source/blender/blenlib/BLI_callbacks.h b/source/blender/blenlib/BLI_callbacks.h index 3db7d2d87f1..f53a4b385b4 100644 --- a/source/blender/blenlib/BLI_callbacks.h +++ b/source/blender/blenlib/BLI_callbacks.h @@ -53,10 +53,6 @@ typedef enum { BLI_CB_EVT_UNDO_POST, BLI_CB_EVT_REDO_PRE, BLI_CB_EVT_REDO_POST, - BLI_CB_EVT_SCENE_UPDATE_PRE, - BLI_CB_EVT_SCENE_UPDATE_POST, - BLI_CB_EVT_GAME_PRE, - BLI_CB_EVT_GAME_POST, BLI_CB_EVT_VERSION_UPDATE, BLI_CB_EVT_TOT } eCbEvent; diff --git a/source/blender/blenlib/BLI_graph.h b/source/blender/blenlib/BLI_graph.h deleted file mode 100644 index 0b316d3c5bb..00000000000 --- a/source/blender/blenlib/BLI_graph.h +++ /dev/null @@ -1,185 +0,0 @@ -/* - * ***** 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - * - * The Original Code is Copyright (C) 2008 Blender Foundation. - * All rights reserved. - * - * Contributor(s): Joshua Leung - * - * ***** END GPL LICENSE BLOCK ***** - */ -#ifndef __BLI_GRAPH_H__ -#define __BLI_GRAPH_H__ - -/** \file BLI_graph.h - * \ingroup bli - */ - -#include "DNA_listBase.h" - -struct BGraph; -struct BNode; -struct BArc; - -struct RadialArc; - -typedef void (*FreeArc)(struct BArc *); -typedef void (*FreeNode)(struct BNode *); -typedef void (*RadialSymmetry)(struct BNode *root_node, struct RadialArc *ring, int total); -typedef void (*AxialSymmetry)(struct BNode *root_node, struct BNode *node1, struct BNode *node2, struct BArc *arc1, struct BArc *arc2); - -/* IF YOU MODIFY THOSE TYPES, YOU NEED TO UPDATE ALL THOSE THAT "INHERIT" FROM THEM - * - * RigGraph, ReebGraph - * - * */ - -typedef struct BGraph { - ListBase arcs; - ListBase nodes; - - float length; - - /* function pointer to deal with custom fonctionnality */ - FreeArc free_arc; - FreeNode free_node; - RadialSymmetry radial_symmetry; - AxialSymmetry axial_symmetry; -} BGraph; - -typedef struct BNode { - void *next, *prev; - float p[3]; - int flag; - - int degree; - struct BArc **arcs; - - int subgraph_index; - - int symmetry_level; - int symmetry_flag; - float symmetry_axis[3]; -} BNode; - -typedef struct BArc { - void *next, *prev; - struct BNode *head, *tail; - int flag; - - float length; - - int symmetry_level; - int symmetry_group; - int symmetry_flag; -} BArc; - -struct BArcIterator; - -void *IT_head(void *iter); -void *IT_tail(void *iter); -void *IT_peek(void *iter, int n); -void *IT_next(void *iter); -void *IT_nextN(void *iter, int n); -void *IT_previous(void *iter); -int IT_stopped(void *iter); - -typedef void * (*HeadFct)(void *iter); -typedef void * (*TailFct)(void *iter); -typedef void * (*PeekFct)(void *iter, int n); -typedef void * (*NextFct)(void *iter); -typedef void * (*NextNFct)(void *iter, int n); -typedef void * (*PreviousFct)(void *iter); -typedef int (*StoppedFct)(void *iter); - -typedef struct BArcIterator { - HeadFct head; - TailFct tail; - PeekFct peek; - NextFct next; - NextNFct nextN; - PreviousFct previous; - StoppedFct stopped; - - float *p, *no; - float size; - - int length; - int index; -} BArcIterator; - -/* Helper structure for radial symmetry */ -typedef struct RadialArc { - struct BArc *arc; - float n[3]; /* normalized vector joining the nodes of the arc */ -} RadialArc; - -BNode *BLI_otherNode(BArc *arc, BNode *node); - -void BLI_freeNode(BGraph *graph, BNode *node); -void BLI_removeNode(BGraph *graph, BNode *node); - -void BLI_removeArc(BGraph *graph, BArc *arc); - -void BLI_flagNodes(BGraph *graph, int flag); -void BLI_flagArcs(BGraph *graph, int flag); - -bool BLI_hasAdjacencyList(BGraph *rg); -void BLI_buildAdjacencyList(BGraph *rg); -void BLI_rebuildAdjacencyListForNode(BGraph *rg, BNode *node); -void BLI_freeAdjacencyList(BGraph *rg); - -int BLI_FlagSubgraphs(BGraph *graph); -void BLI_ReflagSubgraph(BGraph *graph, int old_subgraph, int new_subgraph); - -#define SHAPE_RADIX 10 /* each shape level is encoded this base */ - -int BLI_subtreeShape(BGraph *graph, BNode *node, BArc *rootArc, int include_root); -float BLI_subtreeLength(BNode *node); -void BLI_calcGraphLength(BGraph *graph); - -void BLI_replaceNode(BGraph *graph, BNode *node_src, BNode *node_replaced); -void BLI_replaceNodeInArc(BGraph *graph, BArc *arc, BNode *node_src, BNode *node_replaced); -void BLI_removeDoubleNodes(BGraph *graph, float limit); -BNode *BLI_FindNodeByPosition(BGraph *graph, const float p[3], const float limit); - -BArc *BLI_findConnectedArc(BGraph *graph, BArc *arc, BNode *v); - -bool BLI_isGraphCyclic(BGraph *graph); - -/*------------ Symmetry handling ------------*/ -void BLI_markdownSymmetry(BGraph *graph, BNode *root_node, float limit); - -void BLI_mirrorAlongAxis(float v[3], float center[3], float axis[3]); - -/* BNode symmetry flags */ -#define SYM_TOPOLOGICAL 1 -#define SYM_PHYSICAL 2 - -/* the following two are exclusive */ -#define SYM_AXIAL 4 -#define SYM_RADIAL 8 - -/* BArc symmetry flags - * - * axial symmetry sides */ -#define SYM_SIDE_POSITIVE 1 -#define SYM_SIDE_NEGATIVE 2 -/* Anything higher is the order in radial symmetry */ -#define SYM_SIDE_RADIAL 3 - -#endif /*__BLI_GRAPH_H__*/ diff --git a/source/blender/blenlib/BLI_iterator.h b/source/blender/blenlib/BLI_iterator.h new file mode 100644 index 00000000000..d3d375122a1 --- /dev/null +++ b/source/blender/blenlib/BLI_iterator.h @@ -0,0 +1,62 @@ +/* + * ***** 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + * + * Contributor(s): Dalai Felinto + * + * ***** END GPL LICENSE BLOCK ***** + */ + +#ifndef __BLI_ITERATOR_H__ +#define __BLI_ITERATOR_H__ + +/** \file BLI_iterator.h + * \ingroup bli + */ + +typedef struct BLI_Iterator { + void *current; /* current pointer we iterate over */ + void *data; /* stored data required for this iterator */ + bool skip; + bool valid; +} BLI_Iterator; + +typedef void (*IteratorCb)(BLI_Iterator *iter); +typedef void (*IteratorBeginCb)(BLI_Iterator *iter, void *data_in); + +#define ITER_BEGIN(callback_begin, callback_next, callback_end, _data_in, _type, _instance) \ +{ \ + _type _instance; \ + IteratorCb callback_end_func = callback_end; \ + BLI_Iterator iter_macro; \ + iter_macro.skip = false; \ + iter_macro.valid = true; \ + for (callback_begin(&iter_macro, (_data_in)); \ + iter_macro.valid; \ + callback_next(&iter_macro)) \ + { \ + if (iter_macro.skip) { \ + iter_macro.skip = false; \ + continue; \ + } \ + _instance = (_type ) iter_macro.current; + +#define ITER_END \ + } \ + callback_end_func(&iter_macro); \ +} ((void)0) + +#endif /* __BLI_ITERATOR_H__ */ diff --git a/source/blender/blenlib/BLI_kdopbvh.h b/source/blender/blenlib/BLI_kdopbvh.h index c92f40c67bf..e7cbe05d713 100644 --- a/source/blender/blenlib/BLI_kdopbvh.h +++ b/source/blender/blenlib/BLI_kdopbvh.h @@ -40,6 +40,8 @@ extern "C" { #endif struct BVHTree; +struct DistProjectedAABBPrecalc; + typedef struct BVHTree BVHTree; #define USE_KDOPBVH_WATERTIGHT @@ -101,6 +103,13 @@ typedef bool (*BVHTree_OverlapCallback)(void *userdata, int index_a, int index_b /* callback to range search query */ typedef void (*BVHTree_RangeQuery)(void *userdata, int index, const float co[3], float dist_sq); +/* callback to find nearest projected */ +typedef void (*BVHTree_NearestProjectedCallback)( + void *userdata, int index, + const struct DistProjectedAABBPrecalc *precalc, + const float (*clip_plane)[4], const int clip_plane_len, + BVHTreeNearest *nearest); + /* callbacks to BLI_bvhtree_walk_dfs */ /* return true to traverse into this nodes children, else skip. */ @@ -162,6 +171,12 @@ int BLI_bvhtree_range_query( BVHTree *tree, const float co[3], float radius, BVHTree_RangeQuery callback, void *userdata); +int BLI_bvhtree_find_nearest_projected( + BVHTree *tree, float projmat[4][4], float winsize[2], float mval[2], + float clip_planes[6][4], int clip_num, + BVHTreeNearest *nearest, + BVHTree_NearestProjectedCallback callback, void *userdata); + void BLI_bvhtree_walk_dfs( BVHTree *tree, BVHTree_WalkParentCallback walk_parent_cb, diff --git a/source/blender/blenlib/BLI_link_utils.h b/source/blender/blenlib/BLI_link_utils.h index d469b105f93..5322547ac08 100644 --- a/source/blender/blenlib/BLI_link_utils.h +++ b/source/blender/blenlib/BLI_link_utils.h @@ -35,6 +35,18 @@ list = link; \ } (void)0 +/* Use for append (single linked list, storing the last element). */ +#define BLI_LINKS_APPEND(list, link) { \ + (link)->next = NULL; \ + if ((list)->first) { \ + (list)->last->next = link; \ + } \ + else { \ + (list)->first = link; \ + } \ + (list)->last = link; \ +} (void)0 + #define BLI_LINKS_FREE(list) { \ while (list) { \ void *next = list->next; \ diff --git a/source/blender/blenlib/BLI_math_geom.h b/source/blender/blenlib/BLI_math_geom.h index b3a95d65752..640f3143009 100644 --- a/source/blender/blenlib/BLI_math_geom.h +++ b/source/blender/blenlib/BLI_math_geom.h @@ -121,11 +121,14 @@ float dist_squared_ray_to_seg_v3( const float v0[3], const float v1[3], float r_point[3], float *r_depth); +void aabb_get_near_far_from_plane( + const float plane_no[3], const float bbmin[3], const float bbmax[3], + float bb_near[3], float bb_afar[3]); + struct DistRayAABB_Precalc { float ray_origin[3]; float ray_direction[3]; float ray_inv_dir[3]; - bool sign[3]; }; void dist_squared_ray_to_aabb_v3_precalc( struct DistRayAABB_Precalc *neasrest_precalc, @@ -140,6 +143,24 @@ float dist_squared_ray_to_aabb_v3_simple( const float bb_min[3], const float bb_max[3], float r_point[3], float *r_depth); +struct DistProjectedAABBPrecalc { + float ray_origin[3]; + float ray_direction[3]; + float ray_inv_dir[3]; + float pmat[4][4]; + float mval[2]; +}; +void dist_squared_to_projected_aabb_precalc( + struct DistProjectedAABBPrecalc *precalc, + const float projmat[4][4], const float winsize[2], const float mval[2]); +float dist_squared_to_projected_aabb( + struct DistProjectedAABBPrecalc *data, + const float bbmin[3], const float bbmax[3], + bool r_axis_closest[3]); +float dist_squared_to_projected_aabb_simple( + const float projmat[4][4], const float winsize[2], const float mval[2], + const float bbmin[3], const float bbmax[3]); + float closest_to_line_v2(float r_close[2], const float p[2], const float l1[2], const float l2[2]); float closest_to_line_v3(float r_close[3], const float p[3], const float l1[3], const float l2[3]); void closest_to_line_segment_v2(float r_close[2], const float p[2], const float l1[2], const float l2[2]); @@ -223,6 +244,9 @@ bool isect_ray_plane_v3( float *r_lambda, const bool clip); bool isect_point_planes_v3(float (*planes)[4], int totplane, const float p[3]); +bool isect_point_planes_v3_negated( + const float (*planes)[4], const int totplane, const float p[3]); + bool isect_line_plane_v3( float r_isect_co[3], const float l1[3], const float l2[3], const float plane_co[3], const float plane_no[3]) ATTR_WARN_UNUSED_RESULT; @@ -291,6 +315,11 @@ bool isect_ray_seg_v2( const float v0[2], const float v1[2], float *r_lambda, float *r_u); +bool isect_ray_seg_v3( + const float ray_origin[3], const float ray_direction[3], + const float v0[3], const float v1[3], + float *r_lambda); + /* point in polygon */ bool isect_point_poly_v2(const float pt[2], const float verts[][2], const unsigned int nr, const bool use_holes); bool isect_point_poly_v2_int(const int pt[2], const int verts[][2], const unsigned int nr, const bool use_holes); @@ -326,6 +355,14 @@ bool isect_ray_aabb_v3_simple( float *tmin, float *tmax); /* other */ +#define ISECT_AABB_PLANE_BEHIND_ANY 0 +#define ISECT_AABB_PLANE_CROSS_ANY 1 +#define ISECT_AABB_PLANE_IN_FRONT_ALL 2 + +int isect_aabb_planes_v3( + const float (*planes)[4], const int totplane, + const float bbmin[3], const float bbmax[3]); + bool isect_sweeping_sphere_tri_v3(const float p1[3], const float p2[3], const float radius, const float v0[3], const float v1[3], const float v2[3], float *r_lambda, float ipoint[3]); diff --git a/source/blender/blenlib/BLI_math_matrix.h b/source/blender/blenlib/BLI_math_matrix.h index 9d418749bd1..d4d498590b5 100644 --- a/source/blender/blenlib/BLI_math_matrix.h +++ b/source/blender/blenlib/BLI_math_matrix.h @@ -46,81 +46,112 @@ void unit_m2(float R[2][2]); void unit_m3(float R[3][3]); void unit_m4(float R[4][4]); -void copy_m2_m2(float R[2][2], float A[2][2]); -void copy_m3_m3(float R[3][3], float A[3][3]); -void copy_m4_m4(float R[4][4], float A[4][4]); -void copy_m3_m4(float R[3][3], float A[4][4]); -void copy_m4_m3(float R[4][4], float A[3][3]); +void copy_m2_m2(float R[2][2], const float A[2][2]); +void copy_m3_m3(float R[3][3], const float A[3][3]); +void copy_m4_m4(float R[4][4], const float A[4][4]); +void copy_m3_m4(float R[3][3], const float A[4][4]); +void copy_m4_m3(float R[4][4], const float A[3][3]); /* double->float */ -void copy_m3_m3d(float R[3][3], double A[3][3]); +void copy_m3_m3d(float R[3][3], const double A[3][3]); void swap_m3m3(float A[3][3], float B[3][3]); void swap_m4m4(float A[4][4], float B[4][4]); /******************************** Arithmetic *********************************/ -void add_m3_m3m3(float R[3][3], float A[3][3], float B[3][3]); -void add_m4_m4m4(float R[4][4], float A[4][4], float B[4][4]); +void add_m3_m3m3(float R[3][3], const float A[3][3], const float B[3][3]); +void add_m4_m4m4(float R[4][4], const float A[4][4], const float B[4][4]); -void sub_m3_m3m3(float R[3][3], float A[3][3], float B[3][3]); -void sub_m4_m4m4(float R[4][4], float A[4][4], float B[4][4]); +void sub_m3_m3m3(float R[3][3], const float A[3][3], const float B[3][3]); +void sub_m4_m4m4(float R[4][4], const float A[4][4], const float B[4][4]); -void mul_m3_m3m3(float R[3][3], float A[3][3], float B[3][3]); -void mul_m4_m3m4(float R[4][4], float A[3][3], float B[4][4]); -void mul_m4_m4m3(float R[4][4], float A[4][4], float B[3][3]); -void mul_m4_m4m4(float R[4][4], float A[4][4], float B[4][4]); -void mul_m3_m3m4(float R[3][3], float A[4][4], float B[3][3]); +void mul_m3_m3m3(float R[3][3], const float A[3][3], const float B[3][3]); +void mul_m4_m3m4(float R[4][4], const float A[3][3], const float B[4][4]); +void mul_m4_m4m3(float R[4][4], const float A[4][4], const float B[3][3]); +void mul_m4_m4m4(float R[4][4], const float A[4][4], const float B[4][4]); +void mul_m3_m3m4(float R[3][3], const float A[4][4], const float B[3][3]); + +/* special matrix multiplies + * uniq: R <-- AB, R is neither A nor B + * pre: R <-- AR + * post: R <-- RB + */ +void mul_m3_m3m3_uniq(float R[3][3], const float A[3][3], const float B[3][3]); +void mul_m3_m3_pre(float R[3][3], const float A[3][3]); +void mul_m3_m3_post(float R[3][3], const float B[3][3]); +void mul_m4_m4m4_uniq(float R[4][4], const float A[4][4], const float B[4][4]); +void mul_m4_m4_pre(float R[4][4], const float A[4][4]); +void mul_m4_m4_post(float R[4][4], const float B[4][4]); /* mul_m3_series */ -void _va_mul_m3_series_3(float R[3][3], float M1[3][3], float M2[3][3]) ATTR_NONNULL(); -void _va_mul_m3_series_4(float R[3][3], float M1[3][3], float M2[3][3], float M3[3][3]) ATTR_NONNULL(); -void _va_mul_m3_series_5(float R[3][3], float M1[3][3], float M2[3][3], float M3[3][3], float M4[3][3]) ATTR_NONNULL(); -void _va_mul_m3_series_6(float R[3][3], float M1[3][3], float M2[3][3], float M3[3][3], float M4[3][3], - float M5[3][3]) ATTR_NONNULL(); -void _va_mul_m3_series_7(float R[3][3], float M1[3][3], float M2[3][3], float M3[3][3], float M4[3][3], - float M5[3][3], float M6[3][3]) ATTR_NONNULL(); -void _va_mul_m3_series_8(float R[3][3], float M1[3][3], float M2[3][3], float M3[3][3], float M4[3][3], - float M5[3][3], float M6[3][3], float M7[3][3]) ATTR_NONNULL(); -void _va_mul_m3_series_9(float R[3][3], float M1[3][3], float M2[3][3], float M3[3][3], float M4[3][3], - float M5[3][3], float M6[3][3], float M7[3][3], float M8[3][3]) ATTR_NONNULL(); +void _va_mul_m3_series_3( + float R[3][3], const float M1[3][3], const float M2[3][3]) ATTR_NONNULL(); +void _va_mul_m3_series_4( + float R[3][3], const float M1[3][3], const float M2[3][3], const float M3[3][3]) ATTR_NONNULL(); +void _va_mul_m3_series_5( + float R[3][3], const float M1[3][3], const float M2[3][3], const float M3[3][3], + const float M4[3][3]) ATTR_NONNULL(); +void _va_mul_m3_series_6( + float R[3][3], const float M1[3][3], const float M2[3][3], const float M3[3][3], + const float M4[3][3], const float M5[3][3]) ATTR_NONNULL(); +void _va_mul_m3_series_7( + float R[3][3], const float M1[3][3], const float M2[3][3], const float M3[3][3], + const float M4[3][3], const float M5[3][3], const float M6[3][3]) ATTR_NONNULL(); +void _va_mul_m3_series_8( + float R[3][3], const float M1[3][3], const float M2[3][3], const float M3[3][3], + const float M4[3][3], const float M5[3][3], const float M6[3][3], const float M7[3][3]) ATTR_NONNULL(); +void _va_mul_m3_series_9( + float R[3][3], const float M1[3][3], const float M2[3][3], const float M3[3][3], + const float M4[3][3], const float M5[3][3], const float M6[3][3], const float M7[3][3], + const float M8[3][3]) ATTR_NONNULL(); /* mul_m4_series */ -void _va_mul_m4_series_3(float R[4][4], float M1[4][4], float M2[4][4]) ATTR_NONNULL(); -void _va_mul_m4_series_4(float R[4][4], float M1[4][4], float M2[4][4], float M3[4][4]) ATTR_NONNULL(); -void _va_mul_m4_series_5(float R[4][4], float M1[4][4], float M2[4][4], float M3[4][4], float M4[4][4]) ATTR_NONNULL(); -void _va_mul_m4_series_6(float R[4][4], float M1[4][4], float M2[4][4], float M3[4][4], float M4[4][4], - float M5[4][4]) ATTR_NONNULL(); -void _va_mul_m4_series_7(float R[4][4], float M1[4][4], float M2[4][4], float M3[4][4], float M4[4][4], - float M5[4][4], float M6[4][4]) ATTR_NONNULL(); -void _va_mul_m4_series_8(float R[4][4], float M1[4][4], float M2[4][4], float M3[4][4], float M4[4][4], - float M5[4][4], float M6[4][4], float M7[4][4]) ATTR_NONNULL(); -void _va_mul_m4_series_9(float R[4][4], float M1[4][4], float M2[4][4], float M3[4][4], float M4[4][4], - float M5[4][4], float M6[4][4], float M7[4][4], float M8[4][4]) ATTR_NONNULL(); +void _va_mul_m4_series_3( + float R[4][4], const float M1[4][4], const float M2[4][4]) ATTR_NONNULL(); +void _va_mul_m4_series_4( + float R[4][4], const float M1[4][4], const float M2[4][4], const float M3[4][4]) ATTR_NONNULL(); +void _va_mul_m4_series_5( + float R[4][4], const float M1[4][4], const float M2[4][4], const float M3[4][4], + const float M4[4][4]) ATTR_NONNULL(); +void _va_mul_m4_series_6( + float R[4][4], const float M1[4][4], const float M2[4][4], const float M3[4][4], + const float M4[4][4], const float M5[4][4]) ATTR_NONNULL(); +void _va_mul_m4_series_7( + float R[4][4], const float M1[4][4], const float M2[4][4], const float M3[4][4], + const float M4[4][4], const float M5[4][4], const float M6[4][4]) ATTR_NONNULL(); +void _va_mul_m4_series_8( + float R[4][4], const float M1[4][4], const float M2[4][4], const float M3[4][4], + const float M4[4][4], const float M5[4][4], const float M6[4][4], const float M7[4][4]) ATTR_NONNULL(); +void _va_mul_m4_series_9( + float R[4][4], const float M1[4][4], const float M2[4][4], const float M3[4][4], + const float M4[4][4], const float M5[4][4], const float M6[4][4], const float M7[4][4], + const float M8[4][4]) ATTR_NONNULL(); #define mul_m3_series(...) VA_NARGS_CALL_OVERLOAD(_va_mul_m3_series_, __VA_ARGS__) #define mul_m4_series(...) VA_NARGS_CALL_OVERLOAD(_va_mul_m4_series_, __VA_ARGS__) -void mul_m4_v3(float M[4][4], float r[3]); -void mul_v3_m4v3(float r[3], float M[4][4], const float v[3]); -void mul_v2_m4v3(float r[2], float M[4][4], const float v[3]); -void mul_v2_m2v2(float r[2], float M[2][2], const float v[2]); -void mul_m2v2(float M[2][2], float v[2]); -void mul_mat3_m4_v3(float M[4][4], float r[3]); -void mul_v3_mat3_m4v3(float r[3], float M[4][4], const float v[3]); -void mul_m4_v4(float M[4][4], float r[4]); -void mul_v4_m4v4(float r[4], float M[4][4], const float v[4]); -void mul_project_m4_v3(float M[4][4], float vec[3]); -void mul_v3_project_m4_v3(float r[3], float mat[4][4], const float vec[3]); -void mul_v2_project_m4_v3(float r[2], float M[4][4], const float vec[3]); - -void mul_m3_v2(float m[3][3], float r[2]); -void mul_v2_m3v2(float r[2], float m[3][3], float v[2]); -void mul_m3_v3(float M[3][3], float r[3]); -void mul_v3_m3v3(float r[3], float M[3][3], const float a[3]); -void mul_v2_m3v3(float r[2], float M[3][3], const float a[3]); -void mul_transposed_m3_v3(float M[3][3], float r[3]); -void mul_transposed_mat3_m4_v3(float M[4][4], float r[3]); -void mul_m3_v3_double(float M[3][3], double r[3]); +void mul_m4_v3(const float M[4][4], float r[3]); +void mul_v3_m4v3(float r[3], const float M[4][4], const float v[3]); +void mul_v2_m4v3(float r[2], const float M[4][4], const float v[3]); +void mul_v2_m2v2(float r[2], const float M[2][2], const float v[2]); +void mul_m2v2(const float M[2][2], float v[2]); +void mul_mat3_m4_v3(const float M[4][4], float r[3]); +void mul_v3_mat3_m4v3(float r[3], const float M[4][4], const float v[3]); +void mul_m4_v4(const float M[4][4], float r[4]); +void mul_v4_m4v4(float r[4], const float M[4][4], const float v[4]); +void mul_v4_m4v3(float r[4], const float M[4][4], const float v[3]); /* v has implicit w = 1.0f */ +void mul_project_m4_v3(const float M[4][4], float vec[3]); +void mul_v3_project_m4_v3(float r[3], const float mat[4][4], const float vec[3]); +void mul_v2_project_m4_v3(float r[2], const float M[4][4], const float vec[3]); + +void mul_m3_v2(const float m[3][3], float r[2]); +void mul_v2_m3v2(float r[2], const float m[3][3], const float v[2]); +void mul_m3_v3(const float M[3][3], float r[3]); +void mul_v3_m3v3(float r[3], const float M[3][3], const float a[3]); +void mul_v2_m3v3(float r[2], const float M[3][3], const float a[3]); +void mul_transposed_m3_v3(const float M[3][3], float r[3]); +void mul_transposed_mat3_m4_v3(const float M[4][4], float r[3]); +void mul_m3_v3_double(const float M[3][3], double r[3]); void mul_m3_fl(float R[3][3], float f); void mul_m4_fl(float R[4][4], float f); @@ -131,103 +162,103 @@ void negate_mat3_m4(float R[4][4]); void negate_m4(float R[4][4]); bool invert_m3_ex(float m[3][3], const float epsilon); -bool invert_m3_m3_ex(float m1[3][3], float m2[3][3], const float epsilon); +bool invert_m3_m3_ex(float m1[3][3], const float m2[3][3], const float epsilon); bool invert_m3(float R[3][3]); -bool invert_m3_m3(float R[3][3], float A[3][3]); +bool invert_m3_m3(float R[3][3], const float A[3][3]); bool invert_m4(float R[4][4]); -bool invert_m4_m4(float R[4][4], float A[4][4]); +bool invert_m4_m4(float R[4][4], const float A[4][4]); /* double arithmetic (mixed float/double) */ -void mul_m4_v4d(float M[4][4], double r[4]); -void mul_v4d_m4v4d(double r[4], float M[4][4], double v[4]); +void mul_m4_v4d(const float M[4][4], double r[4]); +void mul_v4d_m4v4d(double r[4], const float M[4][4], const double v[4]); /* double matrix functions (no mixing types) */ -void mul_v3_m3v3_db(double r[3], double M[3][3], const double a[3]); -void mul_m3_v3_db(double M[3][3], double r[3]); +void mul_v3_m3v3_db(double r[3], const double M[3][3], const double a[3]); +void mul_m3_v3_db(const double M[3][3], double r[3]); /****************************** Linear Algebra *******************************/ void transpose_m3(float R[3][3]); -void transpose_m3_m3(float R[3][3], float A[3][3]); -void transpose_m3_m4(float R[3][3], float A[4][4]); +void transpose_m3_m3(float R[3][3], const float A[3][3]); +void transpose_m3_m4(float R[3][3], const float A[4][4]); void transpose_m4(float R[4][4]); -void transpose_m4_m4(float R[4][4], float A[4][4]); +void transpose_m4_m4(float R[4][4], const float A[4][4]); -int compare_m4m4(float mat1[4][4], float mat2[4][4], float limit); +int compare_m4m4(const float mat1[4][4], const float mat2[4][4], float limit); void normalize_m3_ex(float R[3][3], float r_scale[3]) ATTR_NONNULL(); void normalize_m3(float R[3][3]) ATTR_NONNULL(); -void normalize_m3_m3_ex(float R[3][3], float A[3][3], float r_scale[3]) ATTR_NONNULL(); -void normalize_m3_m3(float R[3][3], float A[3][3]) ATTR_NONNULL(); +void normalize_m3_m3_ex(float R[3][3], const float A[3][3], float r_scale[3]) ATTR_NONNULL(); +void normalize_m3_m3(float R[3][3], const float A[3][3]) ATTR_NONNULL(); void normalize_m4_ex(float R[4][4], float r_scale[3]) ATTR_NONNULL(); void normalize_m4(float R[4][4]) ATTR_NONNULL(); -void normalize_m4_m4_ex(float R[4][4], float A[4][4], float r_scale[3]) ATTR_NONNULL(); -void normalize_m4_m4(float R[4][4], float A[4][4]) ATTR_NONNULL(); +void normalize_m4_m4_ex(float R[4][4], const float A[4][4], float r_scale[3]) ATTR_NONNULL(); +void normalize_m4_m4(float R[4][4], const float A[4][4]) ATTR_NONNULL(); void orthogonalize_m3(float R[3][3], int axis); void orthogonalize_m4(float R[4][4], int axis); -bool is_orthogonal_m3(float mat[3][3]); -bool is_orthogonal_m4(float mat[4][4]); -bool is_orthonormal_m3(float mat[3][3]); -bool is_orthonormal_m4(float mat[4][4]); +bool is_orthogonal_m3(const float mat[3][3]); +bool is_orthogonal_m4(const float mat[4][4]); +bool is_orthonormal_m3(const float mat[3][3]); +bool is_orthonormal_m4(const float mat[4][4]); -bool is_uniform_scaled_m3(float mat[3][3]); -bool is_uniform_scaled_m4(float m[4][4]); +bool is_uniform_scaled_m3(const float mat[3][3]); +bool is_uniform_scaled_m4(const float m[4][4]); /* Note: 'adjoint' here means the adjugate (adjunct, "classical adjoint") matrix! * Nowadays 'adjoint' usually refers to the conjugate transpose, * which for real-valued matrices is simply the transpose. */ -void adjoint_m2_m2(float R[2][2], float A[2][2]); -void adjoint_m3_m3(float R[3][3], float A[3][3]); -void adjoint_m4_m4(float R[4][4], float A[4][4]); +void adjoint_m2_m2(float R[2][2], const float A[2][2]); +void adjoint_m3_m3(float R[3][3], const float A[3][3]); +void adjoint_m4_m4(float R[4][4], const float A[4][4]); float determinant_m2(float a, float b, float c, float d); float determinant_m3(float a, float b, float c, float d, float e, float f, float g, float h, float i); -float determinant_m3_array(float m[3][3]); -float determinant_m4(float A[4][4]); +float determinant_m3_array(const float m[3][3]); +float determinant_m4(const float A[4][4]); #define PSEUDOINVERSE_EPSILON 1e-8f void svd_m4(float U[4][4], float s[4], float V[4][4], float A[4][4]); -void pseudoinverse_m4_m4(float Ainv[4][4], float A[4][4], float epsilon); -void pseudoinverse_m3_m3(float Ainv[3][3], float A[3][3], float epsilon); +void pseudoinverse_m4_m4(float Ainv[4][4], const float A[4][4], float epsilon); +void pseudoinverse_m3_m3(float Ainv[3][3], const float A[3][3], float epsilon); -bool has_zero_axis_m4(float matrix[4][4]); +bool has_zero_axis_m4(const float matrix[4][4]); -void invert_m4_m4_safe(float Ainv[4][4], float A[4][4]); +void invert_m4_m4_safe(float Ainv[4][4], const float A[4][4]); /****************************** Transformations ******************************/ void scale_m3_fl(float R[3][3], float scale); void scale_m4_fl(float R[4][4], float scale); -float mat3_to_scale(float M[3][3]); -float mat4_to_scale(float M[4][4]); -float mat4_to_xy_scale(float M[4][4]); +float mat3_to_scale(const float M[3][3]); +float mat4_to_scale(const float M[4][4]); +float mat4_to_xy_scale(const float M[4][4]); void size_to_mat3(float R[3][3], const float size[3]); void size_to_mat4(float R[4][4], const float size[3]); -void mat3_to_size(float r[3], float M[3][3]); -void mat4_to_size(float r[3], float M[4][4]); +void mat3_to_size(float r[3], const float M[3][3]); +void mat4_to_size(float r[3], const float M[4][4]); void translate_m4(float mat[4][4], float tx, float ty, float tz); void rotate_m4(float mat[4][4], const char axis, const float angle); void transform_pivot_set_m4(float mat[4][4], const float pivot[3]); -void mat3_to_rot_size(float rot[3][3], float size[3], float mat3[3][3]); -void mat4_to_loc_rot_size(float loc[3], float rot[3][3], float size[3], float wmat[4][4]); -void mat4_to_loc_quat(float loc[3], float quat[4], float wmat[4][4]); -void mat4_decompose(float loc[3], float quat[4], float size[3], float wmat[4][4]); +void mat3_to_rot_size(float rot[3][3], float size[3], const float mat3[3][3]); +void mat4_to_loc_rot_size(float loc[3], float rot[3][3], float size[3], const float wmat[4][4]); +void mat4_to_loc_quat(float loc[3], float quat[4], const float wmat[4][4]); +void mat4_decompose(float loc[3], float quat[4], float size[3], const float wmat[4][4]); -void mat3_polar_decompose(float mat3[3][3], float r_U[3][3], float r_P[3][3]); +void mat3_polar_decompose(const float mat3[3][3], float r_U[3][3], float r_P[3][3]); void loc_eul_size_to_mat4(float R[4][4], const float loc[3], const float eul[3], const float size[3]); @@ -238,20 +269,20 @@ void loc_quat_size_to_mat4(float R[4][4], void loc_axisangle_size_to_mat4(float R[4][4], const float loc[3], const float axis[4], const float angle, const float size[3]); -void blend_m3_m3m3(float R[3][3], float A[3][3], float B[3][3], const float t); -void blend_m4_m4m4(float R[4][4], float A[4][4], float B[4][4], const float t); +void blend_m3_m3m3(float R[3][3], const float A[3][3], const float B[3][3], const float t); +void blend_m4_m4m4(float R[4][4], const float A[4][4], const float B[4][4], const float t); -void interp_m3_m3m3(float R[3][3], float A[3][3], float B[3][3], const float t); -void interp_m4_m4m4(float R[4][4], float A[4][4], float B[4][4], const float t); +void interp_m3_m3m3(float R[3][3], const float A[3][3], const float B[3][3], const float t); +void interp_m4_m4m4(float R[4][4], const float A[4][4], const float B[4][4], const float t); -bool is_negative_m3(float mat[3][3]); -bool is_negative_m4(float mat[4][4]); +bool is_negative_m3(const float mat[3][3]); +bool is_negative_m4(const float mat[4][4]); -bool is_zero_m3(float mat[3][3]); -bool is_zero_m4(float mat[4][4]); +bool is_zero_m3(const float mat[3][3]); +bool is_zero_m4(const float mat[4][4]); -bool equals_m3m3(float mat1[3][3], float mat2[3][3]); -bool equals_m4m4(float mat1[4][4], float mat2[4][4]); +bool equals_m3m3(const float mat1[3][3], const float mat2[3][3]); +bool equals_m4m4(const float mat1[4][4], const float mat2[4][4]); /* SpaceTransform helper */ typedef struct SpaceTransform { @@ -260,8 +291,8 @@ typedef struct SpaceTransform { } SpaceTransform; -void BLI_space_transform_from_matrices(struct SpaceTransform *data, float local[4][4], float target[4][4]); -void BLI_space_transform_global_from_matrices(struct SpaceTransform *data, float local[4][4], float target[4][4]); +void BLI_space_transform_from_matrices(struct SpaceTransform *data, const float local[4][4], const float target[4][4]); +void BLI_space_transform_global_from_matrices(struct SpaceTransform *data, const float local[4][4], const float target[4][4]); void BLI_space_transform_apply(const struct SpaceTransform *data, float co[3]); void BLI_space_transform_invert(const struct SpaceTransform *data, float co[3]); void BLI_space_transform_apply_normal(const struct SpaceTransform *data, float no[3]); @@ -272,8 +303,8 @@ void BLI_space_transform_invert_normal(const struct SpaceTransform *data, float /*********************************** Other ***********************************/ -void print_m3(const char *str, float M[3][3]); -void print_m4(const char *str, float M[4][4]); +void print_m3(const char *str, const float M[3][3]); +void print_m4(const char *str, const float M[4][4]); #define print_m3_id(M) print_m3(STRINGIFY(M), M) #define print_m4_id(M) print_m4(STRINGIFY(M), M) diff --git a/source/blender/blenlib/BLI_math_vector.h b/source/blender/blenlib/BLI_math_vector.h index 39625346756..ec16a6854c4 100644 --- a/source/blender/blenlib/BLI_math_vector.h +++ b/source/blender/blenlib/BLI_math_vector.h @@ -123,13 +123,13 @@ MINLINE void mul_v4_fl(float r[4], float f); MINLINE void mul_v4_v4fl(float r[3], const float a[3], float f); MINLINE void mul_v2_v2_cw(float r[2], const float mat[2], const float vec[2]); MINLINE void mul_v2_v2_ccw(float r[2], const float mat[2], const float vec[2]); -MINLINE float mul_project_m4_v3_zfac(float mat[4][4], const float co[3]) ATTR_WARN_UNUSED_RESULT; -MINLINE float dot_m3_v3_row_x(float M[3][3], const float a[3]) ATTR_WARN_UNUSED_RESULT; -MINLINE float dot_m3_v3_row_y(float M[3][3], const float a[3]) ATTR_WARN_UNUSED_RESULT; -MINLINE float dot_m3_v3_row_z(float M[3][3], const float a[3]) ATTR_WARN_UNUSED_RESULT; -MINLINE float dot_m4_v3_row_x(float M[4][4], const float a[3]) ATTR_WARN_UNUSED_RESULT; -MINLINE float dot_m4_v3_row_y(float M[4][4], const float a[3]) ATTR_WARN_UNUSED_RESULT; -MINLINE float dot_m4_v3_row_z(float M[4][4], const float a[3]) ATTR_WARN_UNUSED_RESULT; +MINLINE float mul_project_m4_v3_zfac(const float mat[4][4], const float co[3]) ATTR_WARN_UNUSED_RESULT; +MINLINE float dot_m3_v3_row_x(const float M[3][3], const float a[3]) ATTR_WARN_UNUSED_RESULT; +MINLINE float dot_m3_v3_row_y(const float M[3][3], const float a[3]) ATTR_WARN_UNUSED_RESULT; +MINLINE float dot_m3_v3_row_z(const float M[3][3], const float a[3]) ATTR_WARN_UNUSED_RESULT; +MINLINE float dot_m4_v3_row_x(const float M[4][4], const float a[3]) ATTR_WARN_UNUSED_RESULT; +MINLINE float dot_m4_v3_row_y(const float M[4][4], const float a[3]) ATTR_WARN_UNUSED_RESULT; +MINLINE float dot_m4_v3_row_z(const float M[4][4], const float a[3]) ATTR_WARN_UNUSED_RESULT; MINLINE void madd_v2_v2fl(float r[2], const float a[2], float f); MINLINE void madd_v3_v3fl(float r[3], const float a[3], float f); diff --git a/source/blender/blenlib/BLI_memiter.h b/source/blender/blenlib/BLI_memiter.h new file mode 100644 index 00000000000..55eb023313d --- /dev/null +++ b/source/blender/blenlib/BLI_memiter.h @@ -0,0 +1,73 @@ +/* + * ***** 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + * + * ***** END GPL LICENSE BLOCK ***** + */ + +#ifndef __BLI_MEMITER_H__ +#define __BLI_MEMITER_H__ + +/** \file BLI_memiter.h + * \ingroup bli + */ + +#ifdef __cplusplus +extern "C" +{ +#endif + +#include "BLI_sys_types.h" +#include "BLI_compiler_attrs.h" +#include "BLI_compiler_compat.h" + +/* 512kb, good default for small elems. */ +#define BLI_MEMITER_DEFAULT_SIZE (1 << 19) + +struct BLI_memiter; +struct BLI_memiter_chunk; + +typedef struct BLI_memiter BLI_memiter; + +/* warning, ATTR_MALLOC flag on BLI_memiter_alloc causes crash, see: D2756 */ +BLI_memiter *BLI_memiter_create(unsigned int chunk_size) ATTR_MALLOC ATTR_WARN_UNUSED_RESULT; +void *BLI_memiter_alloc(BLI_memiter *mi, unsigned int size) ATTR_RETURNS_NONNULL ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(1); +void BLI_memiter_alloc_from(BLI_memiter *mi, uint elem_size, const void *data_from) ATTR_NONNULL(1, 3); +void *BLI_memiter_calloc(BLI_memiter *mi, unsigned int size) ATTR_RETURNS_NONNULL ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(1); +void BLI_memiter_destroy(BLI_memiter *mi) ATTR_NONNULL(1); +void BLI_memiter_clear(BLI_memiter *mi) ATTR_NONNULL(1); +unsigned int BLI_memiter_count(const BLI_memiter *mi) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(1); + +/* utils */ +void *BLI_memiter_elem_first(BLI_memiter *mi); +void *BLI_memiter_elem_first_size(BLI_memiter *mi, unsigned int *r_size); + +/* private structure */ +typedef struct BLI_memiter_handle { + struct BLI_memiter_elem *elem; + uint elem_left; +} BLI_memiter_handle; + +void BLI_memiter_iter_init(BLI_memiter *mi, BLI_memiter_handle *iter) ATTR_NONNULL(); +bool BLI_memiter_iter_done(const BLI_memiter_handle *iter) ATTR_NONNULL(); +void *BLI_memiter_iter_step(BLI_memiter_handle *iter) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(); +void *BLI_memiter_iter_step_size(BLI_memiter_handle *iter, uint *r_size) ATTR_WARN_UNUSED_RESULT ATTR_NONNULL(); + +#ifdef __cplusplus +} +#endif + +#endif /* __BLI_MEMITER_H__ */ diff --git a/source/blender/blenlib/BLI_rand.h b/source/blender/blenlib/BLI_rand.h index 5bb7ab391b9..612151b7ea2 100644 --- a/source/blender/blenlib/BLI_rand.h +++ b/source/blender/blenlib/BLI_rand.h @@ -64,16 +64,6 @@ void BLI_rng_shuffle_array(struct RNG *rng, void *data, unsigned int elem /** Note that skipping is as slow as generating n numbers! */ void BLI_rng_skip(struct RNG *rng, int n) ATTR_NONNULL(1); -/** Seed for the random number generator, using noise.c hash[] */ -void BLI_srandom(unsigned int seed); - -/** Return a pseudo-random number N where 0<=N<(2^31) */ -int BLI_rand(void) ATTR_WARN_UNUSED_RESULT; - -/** Return a pseudo-random number N where 0.0f<=N<1.0f */ -float BLI_frand(void) ATTR_WARN_UNUSED_RESULT; -void BLI_frand_unit_v3(float v[3]); - /** Return a pseudo-random (hash) float from an integer value */ float BLI_hash_frand(unsigned int seed) ATTR_WARN_UNUSED_RESULT; diff --git a/source/blender/blenlib/BLI_utildefines.h b/source/blender/blenlib/BLI_utildefines.h index 3d4b227ffa7..286e1cc6dd5 100644 --- a/source/blender/blenlib/BLI_utildefines.h +++ b/source/blender/blenlib/BLI_utildefines.h @@ -274,11 +274,15 @@ extern "C" { #define SQUARE(a) ({ \ typeof(a) a_ = (a); \ ((a_) * (a_)); }) +#define CUBE(a) ({ \ + typeof(a) a_ = (a); \ + ((a_) * (a_) * (a_)); }) #else #define ABS(a) ((a) < 0 ? (-(a)) : (a)) #define SQUARE(a) ((a) * (a)) +#define CUBE(a) ((a) * (a) * (a)) #endif @@ -629,6 +633,14 @@ extern bool BLI_memory_is_zero(const void *arr, const size_t arr_size); # define UNLIKELY(x) (x) #endif +/* Expands to an integer constant expression evaluating to a close upper bound + * on the number the number of decimal digits in a value expressible in the + * integer type given by the argument (if it is a type name) or the the integer + * type of the argument (if it is an expression). The meaning of the resulting + * expression is unspecified for other arguments. + * i.e: DECIMAL_DIGITS_BOUND(uchar) is equal to 3. */ +#define DECIMAL_DIGITS_BOUND(t) (241 * sizeof(t) / 100 + 1) + #ifdef __cplusplus } #endif diff --git a/source/blender/blenlib/CMakeLists.txt b/source/blender/blenlib/CMakeLists.txt index e3f5773b1e4..9621a759f3c 100644 --- a/source/blender/blenlib/CMakeLists.txt +++ b/source/blender/blenlib/CMakeLists.txt @@ -52,6 +52,7 @@ set(SRC intern/BLI_linklist.c intern/BLI_linklist_lockfree.c intern/BLI_memarena.c + intern/BLI_memiter.c intern/BLI_mempool.c intern/DLRB_tree.c intern/array_store.c @@ -70,7 +71,6 @@ set(SRC intern/fileops.c intern/fnmatch.c intern/freetypefont.c - intern/graph.c intern/gsqueue.c intern/hash_md5.c intern/hash_mm2a.c @@ -155,13 +155,13 @@ set(SRC BLI_fileops_types.h BLI_fnmatch.h BLI_ghash.h - BLI_graph.h BLI_gsqueue.h BLI_hash.h BLI_hash_md5.h BLI_hash_mm2a.h BLI_hash_mm3.h BLI_heap.h + BLI_iterator.h BLI_jitter_2d.h BLI_kdopbvh.h BLI_kdtree.h @@ -185,6 +185,7 @@ set(SRC BLI_math_statistics.h BLI_math_vector.h BLI_memarena.h + BLI_memiter.h BLI_memory_utils.h BLI_mempool.h BLI_noise.h diff --git a/source/blender/blenlib/intern/BLI_kdopbvh.c b/source/blender/blenlib/intern/BLI_kdopbvh.c index 1676bf5d779..ddfb75fc2ce 100644 --- a/source/blender/blenlib/intern/BLI_kdopbvh.c +++ b/source/blender/blenlib/intern/BLI_kdopbvh.c @@ -167,6 +167,18 @@ typedef struct BVHRayCastData { BVHTreeRayHit hit; } BVHRayCastData; +typedef struct BVHNearestProjectedData { + const BVHTree *tree; + struct DistProjectedAABBPrecalc precalc; + bool closest_axis[3]; + float clip_plane[6][4]; + int clip_plane_len; + BVHTree_NearestProjectedCallback callback; + void *userdata; + BVHTreeNearest nearest; + +} BVHNearestProjectedData; + /** \} */ @@ -501,25 +513,27 @@ static void create_kdop_hull(const BVHTree *tree, BVHNode *node, const float *co } /** - * \note depends on the fact that the BVH's for each face is already build + * \note depends on the fact that the BVH's for each face is already built */ static void refit_kdop_hull(const BVHTree *tree, BVHNode *node, int start, int end) { float newmin, newmax; - float *bv = node->bv; + float *__restrict bv = node->bv; int j; axis_t axis_iter; node_minmax_init(tree, node); for (j = start; j < end; j++) { + float *__restrict node_bv = tree->nodes[j]->bv; + /* for all Axes. */ for (axis_iter = tree->start_axis; axis_iter < tree->stop_axis; axis_iter++) { - newmin = tree->nodes[j]->bv[(2 * axis_iter)]; + newmin = node_bv[(2 * axis_iter)]; if ((newmin < bv[(2 * axis_iter)])) bv[(2 * axis_iter)] = newmin; - newmax = tree->nodes[j]->bv[(2 * axis_iter) + 1]; + newmax = node_bv[(2 * axis_iter) + 1]; if ((newmax > bv[(2 * axis_iter) + 1])) bv[(2 * axis_iter) + 1] = newmax; } @@ -2020,6 +2034,198 @@ int BLI_bvhtree_range_query( /* -------------------------------------------------------------------- */ +/** \name BLI_bvhtree_nearest_projected +* \{ */ + +static void bvhtree_nearest_projected_dfs_recursive( + BVHNearestProjectedData *__restrict data, const BVHNode *node) +{ + if (node->totnode == 0) { + if (data->callback) { + data->callback( + data->userdata, node->index, &data->precalc, + NULL, 0, + &data->nearest); + } + else { + data->nearest.index = node->index; + data->nearest.dist_sq = dist_squared_to_projected_aabb( + &data->precalc, + (float[3]) {node->bv[0], node->bv[2], node->bv[4]}, + (float[3]) {node->bv[1], node->bv[3], node->bv[5]}, + data->closest_axis); + } + } + else { + /* First pick the closest node to recurse into */ + if (data->closest_axis[node->main_axis]) { + for (int i = 0; i != node->totnode; i++) { + const float *bv = node->children[i]->bv; + + if (dist_squared_to_projected_aabb( + &data->precalc, + (float[3]) {bv[0], bv[2], bv[4]}, + (float[3]) {bv[1], bv[3], bv[5]}, + data->closest_axis) <= data->nearest.dist_sq) + { + bvhtree_nearest_projected_dfs_recursive(data, node->children[i]); + } + } + } + else { + for (int i = node->totnode; i--;) { + const float *bv = node->children[i]->bv; + + if (dist_squared_to_projected_aabb( + &data->precalc, + (float[3]) {bv[0], bv[2], bv[4]}, + (float[3]) {bv[1], bv[3], bv[5]}, + data->closest_axis) <= data->nearest.dist_sq) + { + bvhtree_nearest_projected_dfs_recursive(data, node->children[i]); + } + } + } + } +} + +static void bvhtree_nearest_projected_with_clipplane_test_dfs_recursive( + BVHNearestProjectedData *__restrict data, const BVHNode *node) +{ + if (node->totnode == 0) { + if (data->callback) { + data->callback( + data->userdata, node->index, &data->precalc, + data->clip_plane, data->clip_plane_len, + &data->nearest); + } + else { + data->nearest.index = node->index; + data->nearest.dist_sq = dist_squared_to_projected_aabb( + &data->precalc, + (float[3]) {node->bv[0], node->bv[2], node->bv[4]}, + (float[3]) {node->bv[1], node->bv[3], node->bv[5]}, + data->closest_axis); + } + } + else { + /* First pick the closest node to recurse into */ + if (data->closest_axis[node->main_axis]) { + for (int i = 0; i != node->totnode; i++) { + const float *bv = node->children[i]->bv; + const float bb_min[3] = {bv[0], bv[2], bv[4]}; + const float bb_max[3] = {bv[1], bv[3], bv[5]}; + + int isect_type = isect_aabb_planes_v3(data->clip_plane, data->clip_plane_len, bb_min, bb_max); + + if ((isect_type != ISECT_AABB_PLANE_BEHIND_ANY) && dist_squared_to_projected_aabb( + &data->precalc, bb_min, bb_max, + data->closest_axis) <= data->nearest.dist_sq) + { + if (isect_type == ISECT_AABB_PLANE_CROSS_ANY) { + bvhtree_nearest_projected_with_clipplane_test_dfs_recursive(data, node->children[i]); + } + else { + /* ISECT_AABB_PLANE_IN_FRONT_ALL */ + bvhtree_nearest_projected_dfs_recursive(data, node->children[i]); + } + } + } + } + else { + for (int i = node->totnode; i--;) { + const float *bv = node->children[i]->bv; + const float bb_min[3] = {bv[0], bv[2], bv[4]}; + const float bb_max[3] = {bv[1], bv[3], bv[5]}; + + int isect_type = isect_aabb_planes_v3(data->clip_plane, data->clip_plane_len, bb_min, bb_max); + + if (isect_type != ISECT_AABB_PLANE_BEHIND_ANY && dist_squared_to_projected_aabb( + &data->precalc, bb_min, bb_max, + data->closest_axis) <= data->nearest.dist_sq) + { + if (isect_type == ISECT_AABB_PLANE_CROSS_ANY) { + bvhtree_nearest_projected_with_clipplane_test_dfs_recursive(data, node->children[i]); + } + else { + /* ISECT_AABB_PLANE_IN_FRONT_ALL */ + bvhtree_nearest_projected_dfs_recursive(data, node->children[i]); + } + } + } + } + } +} + +int BLI_bvhtree_find_nearest_projected( + BVHTree *tree, float projmat[4][4], float winsize[2], float mval[2], + float clip_plane[6][4], int clip_plane_len, + BVHTreeNearest *nearest, + BVHTree_NearestProjectedCallback callback, void *userdata) +{ + BVHNode *root = tree->nodes[tree->totleaf]; + if (root != NULL) { + BVHNearestProjectedData data; + dist_squared_to_projected_aabb_precalc( + &data.precalc, projmat, winsize, mval); + + data.callback = callback; + data.userdata = userdata; + + if (clip_plane) { + data.clip_plane_len = clip_plane_len; + for (int i = 0; i < data.clip_plane_len; i++) { + copy_v4_v4(data.clip_plane[i], clip_plane[i]); + } + } + else { + data.clip_plane_len = 1; + planes_from_projmat( + projmat, + NULL, NULL, NULL, NULL, + data.clip_plane[0], NULL); + } + + if (nearest) { + memcpy(&data.nearest, nearest, sizeof(*nearest)); + } + else { + data.nearest.index = -1; + data.nearest.dist_sq = FLT_MAX; + } + { + const float bb_min[3] = {root->bv[0], root->bv[2], root->bv[4]}; + const float bb_max[3] = {root->bv[1], root->bv[3], root->bv[5]}; + + int isect_type = isect_aabb_planes_v3(data.clip_plane, data.clip_plane_len, bb_min, bb_max); + + if (isect_type != 0 && dist_squared_to_projected_aabb( + &data.precalc, bb_min, bb_max, + data.closest_axis) <= data.nearest.dist_sq) + { + if (isect_type == 1) { + bvhtree_nearest_projected_with_clipplane_test_dfs_recursive(&data, root); + } + else { + bvhtree_nearest_projected_dfs_recursive(&data, root); + } + } + } + + if (nearest) { + memcpy(nearest, &data.nearest, sizeof(*nearest)); + } + + return data.nearest.index; + } + return -1; +} + +/** \} */ + + +/* -------------------------------------------------------------------- */ + /** \name BLI_bvhtree_walk_dfs * \{ */ diff --git a/source/blender/blenlib/intern/BLI_memiter.c b/source/blender/blenlib/intern/BLI_memiter.c new file mode 100644 index 00000000000..9c5f026f836 --- /dev/null +++ b/source/blender/blenlib/intern/BLI_memiter.c @@ -0,0 +1,357 @@ +/* + * ***** 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + * + * ***** END GPL LICENSE BLOCK ***** + */ + +/** \file blender/blenlib/intern/BLI_memiter.c + * \ingroup bli + * + * Simple, fast memory allocator for allocating many small elements of different sizes + * in fixed size memory chunks, + * although allocations bigger than the chunk size are supported. + * They will reduce the efficiency of this data-structure. + * Elements are pointer aligned. + * + * Supports: + * + * - Allocation of mixed sizes. + * - Iterating over allocations in-order. + * - Clearing for re-use. + * + * Unsupported: + * + * - Freeing individual elements. + * + * \note We could inline iteration stepping, + * but tests show this doesn't give noticeable speedup. + */ + +#include <string.h> +#include <stdlib.h> + +#include "BLI_utildefines.h" + +#include "BLI_memiter.h" /* own include */ + +#include "MEM_guardedalloc.h" + +#include "BLI_strict_flags.h" /* keep last */ + +typedef uintptr_t data_t; +typedef intptr_t offset_t; + +/* Write the chunk terminator on adding each element. + * typically we rely on the 'count' to avoid iterating past the end. */ +// #define USE_TERMINATE_PARANOID + +/* Currently totalloc isnt used. */ + // #define USE_TOTALLOC + +/* pad must be power of two */ +#define PADUP(num, pad) (((num) + ((pad) - 1)) & ~((pad) - 1)) + +typedef struct BLI_memiter_elem { + offset_t size; + data_t data[0]; +} BLI_memiter_elem; + +typedef struct BLI_memiter_chunk { + struct BLI_memiter_chunk *next; + /** + * internal format is: + * ``[next_pointer, size:data, size:data, ..., negative_offset]`` + * + * Where negative offset rewinds to the start. + */ + data_t data[0]; +} BLI_memiter_chunk; + +typedef struct BLI_memiter { + /* A pointer to 'head' is needed so we can iterate in the order allocated. */ + struct BLI_memiter_chunk *head, *tail; + data_t *data_curr; + data_t *data_last; + /* Used unless a large element is requested. + * (which should be very rare!). */ + uint chunk_size_in_bytes_min; + uint count; +#ifdef USE_TOTALLOC + uint totalloc; +#endif +} BLI_memiter; + + +BLI_INLINE uint data_offset_from_size(uint size) +{ + return (PADUP(size, (uint)sizeof(data_t))) / (uint)sizeof(data_t); +} + +static void memiter_set_rewind_offset(BLI_memiter *mi) +{ + BLI_memiter_elem *elem = (BLI_memiter_elem *)mi->data_curr; + elem->size = (offset_t)(((data_t *)mi->tail) - mi->data_curr); + BLI_assert(elem->size < 0); +} + +static void memiter_init(BLI_memiter *mi) +{ + mi->head = NULL; + mi->tail = NULL; + mi->data_curr = NULL; + mi->data_last = NULL; + mi->count = 0; +#ifdef USE_TOTALLOC + mi->totalloc = 0; +#endif +} + +/* -------------------------------------------------------------------- */ + +/** \name Public API's + * \{ */ + +/** + * \param chunk_size_min: Should be a power of two and + * significantly larger than the average element size used. + * + * While allocations of any size are supported, they won't be efficient + * (effectively becoming a single-linked list). + * + * Its intended that many elements can be stored per chunk. + */ +BLI_memiter *BLI_memiter_create(uint chunk_size_min) +{ + BLI_memiter *mi = MEM_mallocN(sizeof(BLI_memiter), "BLI_memiter"); + memiter_init(mi); + + /* Small values are used for tests to check for correctness, + * but otherwise not that useful. */ + const uint slop_space = (sizeof(BLI_memiter_chunk) + MEM_SIZE_OVERHEAD); + if (chunk_size_min >= 1024) { + /* As long as the input is a power of 2, this will give efficient sizes. */ + chunk_size_min -= slop_space; + } + + mi->chunk_size_in_bytes_min = chunk_size_min; + return mi; +} + +void *BLI_memiter_alloc(BLI_memiter *mi, uint elem_size) +{ + const uint data_offset = data_offset_from_size(elem_size); + data_t *data_curr_next = mi->data_curr + (1 + data_offset); + + if (UNLIKELY(mi->data_curr == NULL) || (data_curr_next > mi->data_last)) { + +#ifndef USE_TERMINATE_PARANOID + if (mi->data_curr != NULL) { + memiter_set_rewind_offset(mi); + } +#endif + + uint chunk_size_in_bytes = mi->chunk_size_in_bytes_min; + if (UNLIKELY(chunk_size_in_bytes < elem_size + (uint)sizeof(data_t[2]))) { + chunk_size_in_bytes = elem_size + (uint)sizeof(data_t[2]); + } + uint chunk_size = data_offset_from_size(chunk_size_in_bytes); + BLI_memiter_chunk *chunk = MEM_mallocN( + sizeof(BLI_memiter_chunk) + + (chunk_size * sizeof(data_t)), + "BLI_memiter_chunk"); + + if (mi->head == NULL) { + BLI_assert(mi->tail == NULL); + mi->head = chunk; + } + else { + mi->tail->next = chunk; + } + mi->tail = chunk; + chunk->next = NULL; + + mi->data_curr = chunk->data; + mi->data_last = chunk->data + (chunk_size - 1); + data_curr_next = mi->data_curr + (1 + data_offset); + } + + BLI_assert(data_curr_next <= mi->data_last); + + BLI_memiter_elem *elem = (BLI_memiter_elem *)mi->data_curr; + elem->size = (offset_t)elem_size; + mi->data_curr = data_curr_next; + +#ifdef USE_TERMINATE_PARANOID + memiter_set_rewind_offset(mi); +#endif + + mi->count += 1; + +#ifdef USE_TOTALLOC + mi->totalloc += elem_size; +#endif + + return elem->data; +} + +void *BLI_memiter_calloc(BLI_memiter *mi, uint elem_size) +{ + void *data = BLI_memiter_alloc(mi, elem_size); + memset(data, 0, elem_size); + return data; +} + +void BLI_memiter_alloc_from(BLI_memiter *mi, uint elem_size, const void *data_from) +{ + void *data = BLI_memiter_alloc(mi, elem_size); + memcpy(data, data_from, elem_size); +} + +static void memiter_free_data(BLI_memiter *mi) +{ + BLI_memiter_chunk *chunk = mi->head; + while (chunk) { + BLI_memiter_chunk *chunk_next = chunk->next; + MEM_freeN(chunk); + chunk = chunk_next; + } +} + +void BLI_memiter_destroy(BLI_memiter *mi) +{ + memiter_free_data(mi); + MEM_freeN(mi); +} + +void BLI_memiter_clear(BLI_memiter *mi) +{ + memiter_free_data(mi); + memiter_init(mi); +} + +uint BLI_memiter_count(const BLI_memiter *mi) +{ + return mi->count; +} + +/** \} */ + + +/* -------------------------------------------------------------------- */ + +/** \name Helper API's + * \{ */ + +/* Support direct lookup for first. */ +void *BLI_memiter_elem_first(BLI_memiter *mi) +{ + if (mi->head != NULL) { + BLI_memiter_chunk *chunk = mi->head; + BLI_memiter_elem *elem = (BLI_memiter_elem *)chunk->data; + return elem->data; + } + else { + return NULL; + } +} + +void *BLI_memiter_elem_first_size(BLI_memiter *mi, uint *r_size) +{ + if (mi->head != NULL) { + BLI_memiter_chunk *chunk = mi->head; + BLI_memiter_elem *elem = (BLI_memiter_elem *)chunk->data; + *r_size = (uint)elem->size; + return elem->data; + } + else { + return NULL; + } +} + +/** \} */ + + +/* -------------------------------------------------------------------- */ + +/** \name Iterator API's + * + * \note We could loop over elements until a NULL chunk is found, + * however this means every allocation needs to preemptively run + * #memiter_set_rewind_offset (see #USE_TERMINATE_PARANOID). + * Unless we have a call to finalize allocation (which complicates usage). + * So use a counter instead. + * + * \{ */ + +void BLI_memiter_iter_init(BLI_memiter *mi, BLI_memiter_handle *iter) +{ + iter->elem = mi->head ? (BLI_memiter_elem *)mi->head->data : NULL; + iter->elem_left = mi->count; +} + +bool BLI_memiter_iter_done(const BLI_memiter_handle *iter) +{ + return iter->elem_left != 0; +} + +BLI_INLINE void memiter_chunk_step(BLI_memiter_handle *iter) +{ + BLI_assert(iter->elem->size < 0); + BLI_memiter_chunk *chunk = (BLI_memiter_chunk *)(((data_t *)iter->elem) + iter->elem->size); + chunk = chunk->next; + iter->elem = chunk ? (BLI_memiter_elem *)chunk->data : NULL; + BLI_assert(iter->elem == NULL || iter->elem->size >= 0); +} + +void *BLI_memiter_iter_step_size(BLI_memiter_handle *iter, uint *r_size) +{ + if (iter->elem_left != 0) { + iter->elem_left -= 1; + if (UNLIKELY(iter->elem->size < 0)) { + memiter_chunk_step(iter); + } + BLI_assert(iter->elem->size >= 0); + uint size = (uint)iter->elem->size; + *r_size = size; /* <-- only difference */ + data_t *data = iter->elem->data; + iter->elem = (BLI_memiter_elem *)&data[data_offset_from_size(size)]; + return (void *)data; + } + else { + return NULL; + } +} + +void *BLI_memiter_iter_step(BLI_memiter_handle *iter) +{ + if (iter->elem_left != 0) { + iter->elem_left -= 1; + if (UNLIKELY(iter->elem->size < 0)) { + memiter_chunk_step(iter); + } + BLI_assert(iter->elem->size >= 0); + uint size = (uint)iter->elem->size; + data_t *data = iter->elem->data; + iter->elem = (BLI_memiter_elem *)&data[data_offset_from_size(size)]; + return (void *)data; + } + else { + return NULL; + } +} + +/** \} */ diff --git a/source/blender/blenlib/intern/graph.c b/source/blender/blenlib/intern/graph.c deleted file mode 100644 index e346b8ec003..00000000000 --- a/source/blender/blenlib/intern/graph.c +++ /dev/null @@ -1,1016 +0,0 @@ -/* - * ***** 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., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. - * - * Contributor(s): Martin Poirier - * - * ***** END GPL LICENSE BLOCK ***** - * graph.c: Common graph interface and methods - */ - -/** \file blender/blenlib/intern/graph.c - * \ingroup bli - */ - -#include <float.h> -#include <math.h> - -#include "MEM_guardedalloc.h" - -#include "BLI_utildefines.h" -#include "BLI_listbase.h" -#include "BLI_graph.h" -#include "BLI_math.h" - - -static void testRadialSymmetry(BGraph *graph, BNode *root_node, RadialArc *ring, int total, float axis[3], float limit, int group); - -static void handleAxialSymmetry(BGraph *graph, BNode *root_node, int depth, float axis[3], float limit); -static void testAxialSymmetry(BGraph *graph, BNode *root_node, BNode *node1, BNode *node2, BArc *arc1, BArc *arc2, float axis[3], float limit, int group); -static void flagAxialSymmetry(BNode *root_node, BNode *end_node, BArc *arc, int group); - -void BLI_freeNode(BGraph *graph, BNode *node) -{ - if (node->arcs) { - MEM_freeN(node->arcs); - } - - if (graph->free_node) { - graph->free_node(node); - } -} - -void BLI_removeNode(BGraph *graph, BNode *node) -{ - BLI_freeNode(graph, node); - BLI_freelinkN(&graph->nodes, node); -} - -BNode *BLI_otherNode(BArc *arc, BNode *node) -{ - return (arc->head == node) ? arc->tail : arc->head; -} - -void BLI_removeArc(BGraph *graph, BArc *arc) -{ - if (graph->free_arc) { - graph->free_arc(arc); - } - - BLI_freelinkN(&graph->arcs, arc); -} - -void BLI_flagNodes(BGraph *graph, int flag) -{ - BNode *node; - - for (node = graph->nodes.first; node; node = node->next) { - node->flag = flag; - } -} - -void BLI_flagArcs(BGraph *graph, int flag) -{ - BArc *arc; - - for (arc = graph->arcs.first; arc; arc = arc->next) { - arc->flag = flag; - } -} - -static void addArcToNodeAdjacencyList(BNode *node, BArc *arc) -{ - node->arcs[node->flag] = arc; - node->flag++; -} - -void BLI_buildAdjacencyList(BGraph *graph) -{ - BNode *node; - BArc *arc; - - for (node = graph->nodes.first; node; node = node->next) { - if (node->arcs != NULL) { - MEM_freeN(node->arcs); - } - - node->arcs = MEM_callocN((node->degree) * sizeof(BArc *), "adjacency list"); - - /* temporary use to indicate the first index available in the lists */ - node->flag = 0; - } - - for (arc = graph->arcs.first; arc; arc = arc->next) { - addArcToNodeAdjacencyList(arc->head, arc); - addArcToNodeAdjacencyList(arc->tail, arc); - } - - for (node = graph->nodes.first; node; node = node->next) { - if (node->degree != node->flag) { - printf("error in node [%p]. Added only %i arcs out of %i\n", (void *)node, node->flag, node->degree); - } - } -} - -void BLI_rebuildAdjacencyListForNode(BGraph *graph, BNode *node) -{ - BArc *arc; - - if (node->arcs != NULL) { - MEM_freeN(node->arcs); - } - - node->arcs = MEM_callocN((node->degree) * sizeof(BArc *), "adjacency list"); - - /* temporary use to indicate the first index available in the lists */ - node->flag = 0; - - for (arc = graph->arcs.first; arc; arc = arc->next) { - if (arc->head == node) { - addArcToNodeAdjacencyList(arc->head, arc); - } - else if (arc->tail == node) { - addArcToNodeAdjacencyList(arc->tail, arc); - } - } - - if (node->degree != node->flag) { - printf("error in node [%p]. Added only %i arcs out of %i\n", (void *)node, node->flag, node->degree); - } -} - -void BLI_freeAdjacencyList(BGraph *graph) -{ - BNode *node; - - for (node = graph->nodes.first; node; node = node->next) { - if (node->arcs != NULL) { - MEM_freeN(node->arcs); - node->arcs = NULL; - } - } -} - -bool BLI_hasAdjacencyList(BGraph *graph) -{ - BNode *node; - - for (node = graph->nodes.first; node; node = node->next) { - if (node->arcs == NULL) { - return false; - } - } - - return true; -} - -void BLI_replaceNodeInArc(BGraph *graph, BArc *arc, BNode *node_src, BNode *node_replaced) -{ - if (arc->head == node_replaced) { - arc->head = node_src; - node_src->degree++; - } - - if (arc->tail == node_replaced) { - arc->tail = node_src; - node_src->degree++; - } - - if (arc->head == arc->tail) { - node_src->degree -= 2; - - graph->free_arc(arc); - BLI_freelinkN(&graph->arcs, arc); - } - - if (node_replaced->degree == 0) { - BLI_removeNode(graph, node_replaced); - } -} - -void BLI_replaceNode(BGraph *graph, BNode *node_src, BNode *node_replaced) -{ - BArc *arc, *next_arc; - - for (arc = graph->arcs.first; arc; arc = next_arc) { - next_arc = arc->next; - - if (arc->head == node_replaced) { - arc->head = node_src; - node_replaced->degree--; - node_src->degree++; - } - - if (arc->tail == node_replaced) { - arc->tail = node_src; - node_replaced->degree--; - node_src->degree++; - } - - if (arc->head == arc->tail) { - node_src->degree -= 2; - - graph->free_arc(arc); - BLI_freelinkN(&graph->arcs, arc); - } - } - - if (node_replaced->degree == 0) { - BLI_removeNode(graph, node_replaced); - } -} - -void BLI_removeDoubleNodes(BGraph *graph, float limit) -{ - const float limit_sq = limit * limit; - BNode *node_src, *node_replaced; - - for (node_src = graph->nodes.first; node_src; node_src = node_src->next) { - for (node_replaced = graph->nodes.first; node_replaced; node_replaced = node_replaced->next) { - if (node_replaced != node_src && len_squared_v3v3(node_replaced->p, node_src->p) <= limit_sq) { - BLI_replaceNode(graph, node_src, node_replaced); - } - } - } - -} - -BNode *BLI_FindNodeByPosition(BGraph *graph, const float p[3], const float limit) -{ - const float limit_sq = limit * limit; - BNode *closest_node = NULL, *node; - float min_distance = 0.0f; - - for (node = graph->nodes.first; node; node = node->next) { - float distance = len_squared_v3v3(p, node->p); - if (distance <= limit_sq && (closest_node == NULL || distance < min_distance)) { - closest_node = node; - min_distance = distance; - } - } - - return closest_node; -} -/************************************* SUBGRAPH DETECTION **********************************************/ - -static void flagSubgraph(BNode *node, int subgraph) -{ - if (node->subgraph_index == 0) { - BArc *arc; - int i; - - node->subgraph_index = subgraph; - - for (i = 0; i < node->degree; i++) { - arc = node->arcs[i]; - flagSubgraph(BLI_otherNode(arc, node), subgraph); - } - } -} - -int BLI_FlagSubgraphs(BGraph *graph) -{ - BNode *node; - int subgraph = 0; - - if (BLI_hasAdjacencyList(graph) == 0) { - BLI_buildAdjacencyList(graph); - } - - for (node = graph->nodes.first; node; node = node->next) { - node->subgraph_index = 0; - } - - for (node = graph->nodes.first; node; node = node->next) { - if (node->subgraph_index == 0) { - subgraph++; - flagSubgraph(node, subgraph); - } - } - - return subgraph; -} - -void BLI_ReflagSubgraph(BGraph *graph, int old_subgraph, int new_subgraph) -{ - BNode *node; - - for (node = graph->nodes.first; node; node = node->next) { - if (node->flag == old_subgraph) { - node->flag = new_subgraph; - } - } -} - -/*************************************** CYCLE DETECTION ***********************************************/ - -static bool detectCycle(BNode *node, BArc *src_arc) -{ - bool value = false; - - if (node->flag == 0) { - int i; - - /* mark node as visited */ - node->flag = 1; - - for (i = 0; i < node->degree && value == 0; i++) { - BArc *arc = node->arcs[i]; - - /* don't go back on the source arc */ - if (arc != src_arc) { - value = detectCycle(BLI_otherNode(arc, node), arc); - } - } - } - else { - value = true; - } - - return value; -} - -bool BLI_isGraphCyclic(BGraph *graph) -{ - BNode *node; - bool value = false; - - /* NEED TO CHECK IF ADJACENCY LIST EXIST */ - - /* Mark all nodes as not visited */ - BLI_flagNodes(graph, 0); - - /* detectCycles in subgraphs */ - for (node = graph->nodes.first; node && value == false; node = node->next) { - /* only for nodes in subgraphs that haven't been visited yet */ - if (node->flag == 0) { - value = value || detectCycle(node, NULL); - } - } - - return value; -} - -BArc *BLI_findConnectedArc(BGraph *graph, BArc *arc, BNode *v) -{ - BArc *nextArc; - - for (nextArc = graph->arcs.first; nextArc; nextArc = nextArc->next) { - if (arc != nextArc && (nextArc->head == v || nextArc->tail == v)) { - break; - } - } - - return nextArc; -} - -/*********************************** GRAPH AS TREE FUNCTIONS *******************************************/ - -static int subtreeShape(BNode *node, BArc *rootArc, int include_root) -{ - int depth = 0; - - node->flag = 1; - - if (include_root) { - BNode *newNode = BLI_otherNode(rootArc, node); - return subtreeShape(newNode, rootArc, 0); - } - else { - /* Base case, no arcs leading away */ - if (node->arcs == NULL || *(node->arcs) == NULL) { - return 0; - } - else { - int i; - - for (i = 0; i < node->degree; i++) { - BArc *arc = node->arcs[i]; - BNode *newNode = BLI_otherNode(arc, node); - - /* stop immediate and cyclic backtracking */ - if (arc != rootArc && newNode->flag == 0) { - depth += subtreeShape(newNode, arc, 0); - } - } - } - - return SHAPE_RADIX * depth + 1; - } -} - -int BLI_subtreeShape(BGraph *graph, BNode *node, BArc *rootArc, int include_root) -{ - BLI_flagNodes(graph, 0); - return subtreeShape(node, rootArc, include_root); -} - -float BLI_subtreeLength(BNode *node) -{ - float length = 0; - int i; - - node->flag = 0; /* flag node as visited */ - - for (i = 0; i < node->degree; i++) { - BArc *arc = node->arcs[i]; - BNode *other_node = BLI_otherNode(arc, node); - - if (other_node->flag != 0) { - float subgraph_length = arc->length + BLI_subtreeLength(other_node); - length = MAX2(length, subgraph_length); - } - } - - return length; -} - -void BLI_calcGraphLength(BGraph *graph) -{ - float length = 0; - int nb_subgraphs; - int i; - - nb_subgraphs = BLI_FlagSubgraphs(graph); - - for (i = 1; i <= nb_subgraphs; i++) { - BNode *node; - - for (node = graph->nodes.first; node; node = node->next) { - /* start on an external node of the subgraph */ - if (node->subgraph_index == i && node->degree == 1) { - float subgraph_length = BLI_subtreeLength(node); - length = MAX2(length, subgraph_length); - break; - } - } - } - - graph->length = length; -} - -/********************************* SYMMETRY DETECTION **************************************************/ - -static void markdownSymmetryArc(BGraph *graph, BArc *arc, BNode *node, int level, float limit); - -void BLI_mirrorAlongAxis(float v[3], float center[3], float axis[3]) -{ - float dv[3], pv[3]; - - sub_v3_v3v3(dv, v, center); - project_v3_v3v3(pv, dv, axis); - mul_v3_fl(pv, -2); - add_v3_v3(v, pv); -} - -static void testRadialSymmetry(BGraph *graph, BNode *root_node, RadialArc *ring, int total, float axis[3], float limit, int group) -{ - const float limit_sq = limit * limit; - int symmetric = 1; - int i; - - /* sort ring by angle */ - for (i = 0; i < total - 1; i++) { - float minAngle = FLT_MAX; - int minIndex = -1; - int j; - - for (j = i + 1; j < total; j++) { - float angle = dot_v3v3(ring[i].n, ring[j].n); - - /* map negative values to 1..2 */ - if (angle < 0) { - angle = 1 - angle; - } - - if (angle < minAngle) { - minIndex = j; - minAngle = angle; - } - } - - /* swap if needed */ - if (minIndex != i + 1) { - RadialArc tmp; - tmp = ring[i + 1]; - ring[i + 1] = ring[minIndex]; - ring[minIndex] = tmp; - } - } - - for (i = 0; i < total && symmetric; i++) { - BNode *node1, *node2; - float tangent[3]; - float normal[3]; - float p[3]; - int j = (i + 1) % total; /* next arc in the circular list */ - - add_v3_v3v3(tangent, ring[i].n, ring[j].n); - cross_v3_v3v3(normal, tangent, axis); - - node1 = BLI_otherNode(ring[i].arc, root_node); - node2 = BLI_otherNode(ring[j].arc, root_node); - - copy_v3_v3(p, node2->p); - BLI_mirrorAlongAxis(p, root_node->p, normal); - - /* check if it's within limit before continuing */ - if (len_squared_v3v3(node1->p, p) > limit_sq) { - symmetric = 0; - } - - } - - if (symmetric) { - /* mark node as symmetric physically */ - copy_v3_v3(root_node->symmetry_axis, axis); - root_node->symmetry_flag |= SYM_PHYSICAL; - root_node->symmetry_flag |= SYM_RADIAL; - - /* FLAG SYMMETRY GROUP */ - for (i = 0; i < total; i++) { - ring[i].arc->symmetry_group = group; - ring[i].arc->symmetry_flag = SYM_SIDE_RADIAL + i; - } - - if (graph->radial_symmetry) { - graph->radial_symmetry(root_node, ring, total); - } - } -} - -static void handleRadialSymmetry(BGraph *graph, BNode *root_node, int depth, float axis[3], float limit) -{ - RadialArc *ring = NULL; - RadialArc *unit; - int total = 0; - int group; - int first; - int i; - - /* mark topological symmetry */ - root_node->symmetry_flag |= SYM_TOPOLOGICAL; - - /* total the number of arcs in the symmetry ring */ - for (i = 0; i < root_node->degree; i++) { - BArc *connectedArc = root_node->arcs[i]; - - /* depth is store as a negative in flag. symmetry level is positive */ - if (connectedArc->symmetry_level == -depth) { - total++; - } - } - - ring = MEM_callocN(sizeof(RadialArc) * total, "radial symmetry ring"); - unit = ring; - - /* fill in the ring */ - for (i = 0; i < root_node->degree; i++) { - BArc *connectedArc = root_node->arcs[i]; - - /* depth is store as a negative in flag. symmetry level is positive */ - if (connectedArc->symmetry_level == -depth) { - BNode *otherNode = BLI_otherNode(connectedArc, root_node); - float vec[3]; - - unit->arc = connectedArc; - - /* project the node to node vector on the symmetry plane */ - sub_v3_v3v3(unit->n, otherNode->p, root_node->p); - project_v3_v3v3(vec, unit->n, axis); - sub_v3_v3v3(unit->n, unit->n, vec); - - normalize_v3(unit->n); - - unit++; - } - } - - /* sort ring by arc length - * using a rather bogus insertion sort - * but rings will never get too big to matter - * */ - for (i = 0; i < total; i++) { - int j; - - for (j = i - 1; j >= 0; j--) { - BArc *arc1, *arc2; - - arc1 = ring[j].arc; - arc2 = ring[j + 1].arc; - - if (arc1->length > arc2->length) { - /* swap with smaller */ - RadialArc tmp; - - tmp = ring[j + 1]; - ring[j + 1] = ring[j]; - ring[j] = tmp; - } - else { - break; - } - } - } - - /* Dispatch to specific symmetry tests */ - first = 0; - group = 0; - - for (i = 1; i < total; i++) { - int dispatch = 0; - int last = i - 1; - - if (fabsf(ring[first].arc->length - ring[i].arc->length) > limit) { - dispatch = 1; - } - - /* if not dispatching already and on last arc - * Dispatch using current arc as last - */ - if (dispatch == 0 && i == total - 1) { - last = i; - dispatch = 1; - } - - if (dispatch) { - int sub_total = last - first + 1; - - group += 1; - - if (sub_total == 1) { - group -= 1; /* not really a group so decrement */ - /* NOTHING TO DO */ - } - else if (sub_total == 2) { - BArc *arc1, *arc2; - BNode *node1, *node2; - - arc1 = ring[first].arc; - arc2 = ring[last].arc; - - node1 = BLI_otherNode(arc1, root_node); - node2 = BLI_otherNode(arc2, root_node); - - testAxialSymmetry(graph, root_node, node1, node2, arc1, arc2, axis, limit, group); - } - else if (sub_total != total) /* allocate a new sub ring if needed */ { - RadialArc *sub_ring = MEM_callocN(sizeof(RadialArc) * sub_total, "radial symmetry ring"); - int sub_i; - - /* fill in the sub ring */ - for (sub_i = 0; sub_i < sub_total; sub_i++) { - sub_ring[sub_i] = ring[first + sub_i]; - } - - testRadialSymmetry(graph, root_node, sub_ring, sub_total, axis, limit, group); - - MEM_freeN(sub_ring); - } - else if (sub_total == total) { - testRadialSymmetry(graph, root_node, ring, total, axis, limit, group); - } - - first = i; - } - } - - - MEM_freeN(ring); -} - -static void flagAxialSymmetry(BNode *root_node, BNode *end_node, BArc *arc, int group) -{ - float vec[3]; - - arc->symmetry_group = group; - - sub_v3_v3v3(vec, end_node->p, root_node->p); - - if (dot_v3v3(vec, root_node->symmetry_axis) < 0) { - arc->symmetry_flag |= SYM_SIDE_NEGATIVE; - } - else { - arc->symmetry_flag |= SYM_SIDE_POSITIVE; - } -} - -static void testAxialSymmetry(BGraph *graph, BNode *root_node, BNode *node1, BNode *node2, BArc *arc1, BArc *arc2, float axis[3], float limit, int group) -{ - const float limit_sq = limit * limit; - float nor[3], vec[3], p[3]; - - sub_v3_v3v3(p, node1->p, root_node->p); - cross_v3_v3v3(nor, p, axis); - - sub_v3_v3v3(p, root_node->p, node2->p); - cross_v3_v3v3(vec, p, axis); - add_v3_v3(vec, nor); - - cross_v3_v3v3(nor, vec, axis); - - if (fabsf(nor[0]) > fabsf(nor[1]) && fabsf(nor[0]) > fabsf(nor[2]) && nor[0] < 0) { - negate_v3(nor); - } - else if (fabsf(nor[1]) > fabsf(nor[0]) && fabsf(nor[1]) > fabsf(nor[2]) && nor[1] < 0) { - negate_v3(nor); - } - else if (fabsf(nor[2]) > fabsf(nor[1]) && fabsf(nor[2]) > fabsf(nor[0]) && nor[2] < 0) { - negate_v3(nor); - } - - /* mirror node2 along axis */ - copy_v3_v3(p, node2->p); - BLI_mirrorAlongAxis(p, root_node->p, nor); - - /* check if it's within limit before continuing */ - if (len_squared_v3v3(node1->p, p) <= limit_sq) { - /* mark node as symmetric physically */ - copy_v3_v3(root_node->symmetry_axis, nor); - root_node->symmetry_flag |= SYM_PHYSICAL; - root_node->symmetry_flag |= SYM_AXIAL; - - /* flag side on arcs */ - flagAxialSymmetry(root_node, node1, arc1, group); - flagAxialSymmetry(root_node, node2, arc2, group); - - if (graph->axial_symmetry) { - graph->axial_symmetry(root_node, node1, node2, arc1, arc2); - } - } - else { - /* NOT SYMMETRIC */ - } -} - -static void handleAxialSymmetry(BGraph *graph, BNode *root_node, int depth, float axis[3], float limit) -{ - BArc *arc1 = NULL, *arc2 = NULL; - BNode *node1 = NULL, *node2 = NULL; - int i; - - /* mark topological symmetry */ - root_node->symmetry_flag |= SYM_TOPOLOGICAL; - - for (i = 0; i < root_node->degree; i++) { - BArc *connectedArc = root_node->arcs[i]; - - /* depth is store as a negative in flag. symmetry level is positive */ - if (connectedArc->symmetry_level == -depth) { - if (arc1 == NULL) { - arc1 = connectedArc; - node1 = BLI_otherNode(arc1, root_node); - } - else { - arc2 = connectedArc; - node2 = BLI_otherNode(arc2, root_node); - break; /* Can stop now, the two arcs have been found */ - } - } - } - - /* shouldn't happen, but just to be sure */ - if (node1 == NULL || node2 == NULL) { - return; - } - - testAxialSymmetry(graph, root_node, node1, node2, arc1, arc2, axis, limit, 1); -} - -static void markdownSecondarySymmetry(BGraph *graph, BNode *node, int depth, int level, float limit) -{ - float axis[3] = {0, 0, 0}; - int count = 0; - int i; - - /* count the number of branches in this symmetry group - * and determinate the axis of symmetry - */ - for (i = 0; i < node->degree; i++) { - BArc *connectedArc = node->arcs[i]; - - /* depth is store as a negative in flag. symmetry level is positive */ - if (connectedArc->symmetry_level == -depth) { - count++; - } - /* If arc is on the axis */ - else if (connectedArc->symmetry_level == level) { - add_v3_v3(axis, connectedArc->head->p); - sub_v3_v3v3(axis, axis, connectedArc->tail->p); - } - } - - normalize_v3(axis); - - /* Split between axial and radial symmetry */ - if (count == 2) { - handleAxialSymmetry(graph, node, depth, axis, limit); - } - else { - handleRadialSymmetry(graph, node, depth, axis, limit); - } - - /* markdown secondary symetries */ - for (i = 0; i < node->degree; i++) { - BArc *connectedArc = node->arcs[i]; - - if (connectedArc->symmetry_level == -depth) { - /* markdown symmetry for branches corresponding to the depth */ - markdownSymmetryArc(graph, connectedArc, node, level + 1, limit); - } - } -} - -static void markdownSymmetryArc(BGraph *graph, BArc *arc, BNode *node, int level, float limit) -{ - int i; - - /* if arc is null, we start straight from a node */ - if (arc) { - arc->symmetry_level = level; - - node = BLI_otherNode(arc, node); - } - - for (i = 0; i < node->degree; i++) { - BArc *connectedArc = node->arcs[i]; - - if (connectedArc != arc) { - BNode *connectedNode = BLI_otherNode(connectedArc, node); - - /* symmetry level is positive value, negative values is subtree depth */ - connectedArc->symmetry_level = -BLI_subtreeShape(graph, connectedNode, connectedArc, 0); - } - } - - arc = NULL; - - for (i = 0; i < node->degree; i++) { - int issymmetryAxis = 0; - BArc *connectedArc = node->arcs[i]; - - /* only arcs not already marked as symetric */ - if (connectedArc->symmetry_level < 0) { - int j; - - /* true by default */ - issymmetryAxis = 1; - - for (j = 0; j < node->degree; j++) { - BArc *otherArc = node->arcs[j]; - - /* different arc, same depth */ - if (otherArc != connectedArc && otherArc->symmetry_level == connectedArc->symmetry_level) { - /* not on the symmetry axis */ - issymmetryAxis = 0; - break; - } - } - } - - /* arc could be on the symmetry axis */ - if (issymmetryAxis == 1) { - /* no arc as been marked previously, keep this one */ - if (arc == NULL) { - arc = connectedArc; - } - else if (connectedArc->symmetry_level < arc->symmetry_level) { - /* go with more complex subtree as symmetry arc */ - arc = connectedArc; - } - } - } - - /* go down the arc continuing the symmetry axis */ - if (arc) { - markdownSymmetryArc(graph, arc, node, level, limit); - } - - - /* secondary symmetry */ - for (i = 0; i < node->degree; i++) { - BArc *connectedArc = node->arcs[i]; - - /* only arcs not already marked as symetric and is not the next arc on the symmetry axis */ - if (connectedArc->symmetry_level < 0) { - /* subtree depth is store as a negative value in the symmetry */ - markdownSecondarySymmetry(graph, node, -connectedArc->symmetry_level, level, limit); - } - } -} - -void BLI_markdownSymmetry(BGraph *graph, BNode *root_node, float limit) -{ - BNode *node; - BArc *arc; - - if (root_node == NULL) { - return; - } - - if (BLI_isGraphCyclic(graph)) { - return; - } - - /* mark down all arcs as non-symetric */ - BLI_flagArcs(graph, 0); - - /* mark down all nodes as not on the symmetry axis */ - BLI_flagNodes(graph, 0); - - node = root_node; - - /* sanity check REMOVE ME */ - if (node->degree > 0) { - arc = node->arcs[0]; - - if (node->degree == 1) { - markdownSymmetryArc(graph, arc, node, 1, limit); - } - else { - markdownSymmetryArc(graph, NULL, node, 1, limit); - } - - - - /* mark down non-symetric arcs */ - for (arc = graph->arcs.first; arc; arc = arc->next) { - if (arc->symmetry_level < 0) { - arc->symmetry_level = 0; - } - else { - /* mark down nodes with the lowest level symmetry axis */ - if (arc->head->symmetry_level == 0 || arc->head->symmetry_level > arc->symmetry_level) { - arc->head->symmetry_level = arc->symmetry_level; - } - if (arc->tail->symmetry_level == 0 || arc->tail->symmetry_level > arc->symmetry_level) { - arc->tail->symmetry_level = arc->symmetry_level; - } - } - } - } -} - -void *IT_head(void *arg) -{ - BArcIterator *iter = (BArcIterator *)arg; - return iter->head(iter); -} - -void *IT_tail(void *arg) -{ - BArcIterator *iter = (BArcIterator *)arg; - return iter->tail(iter); -} - -void *IT_peek(void *arg, int n) -{ - BArcIterator *iter = (BArcIterator *)arg; - - if (iter->index + n < 0) { - return iter->head(iter); - } - else if (iter->index + n >= iter->length) { - return iter->tail(iter); - } - else { - return iter->peek(iter, n); - } -} - -void *IT_next(void *arg) -{ - BArcIterator *iter = (BArcIterator *)arg; - return iter->next(iter); -} - -void *IT_nextN(void *arg, int n) -{ - BArcIterator *iter = (BArcIterator *)arg; - return iter->nextN(iter, n); -} - -void *IT_previous(void *arg) -{ - BArcIterator *iter = (BArcIterator *)arg; - return iter->previous(iter); -} - -int IT_stopped(void *arg) -{ - BArcIterator *iter = (BArcIterator *)arg; - return iter->stopped(iter); -} diff --git a/source/blender/blenlib/intern/math_color.c b/source/blender/blenlib/intern/math_color.c index ef9e46fc62d..9399646bedf 100644 --- a/source/blender/blenlib/intern/math_color.c +++ b/source/blender/blenlib/intern/math_color.c @@ -380,8 +380,10 @@ void xyz_to_rgb(float xc, float yc, float zc, float *r, float *g, float *b, int } } -/* we define a 'cpack' here as a (3 byte color code) number that can be expressed like 0xFFAA66 or so. - * for that reason it is sensitive for endianness... with this function it works correctly +/** + * We define a 'cpack' here as a (3 byte color code) number that can be expressed like 0xFFAA66 or so. + * for that reason it is sensitive for endianness... with this function it works correctly. + * \see #imm_cpack */ unsigned int hsv_to_cpack(float h, float s, float v) @@ -675,135 +677,80 @@ void rgb_to_lab(float r, float g, float b, float *ll, float *la, float *lb) xyz_to_lab(x, y, z, ll, la, lb); } -static void xyz_to_lms(float x, float y, float z, float *l, float *m, float *s) -{ - *l = 0.3897f * x + 0.6890f * y - 0.0787f * z; - *m = -0.2298f * x + 1.1834f * y + 0.0464f * z; - *s = z; -} - -static void lms_to_xyz(float l, float m, float s, float *x, float *y, float *z) -{ - *x = 1.9102f * l - 1.1121f * m + 0.2019f * s; - *y = 0.3709f * l + 0.6290f * m + 0.0000f * s; - *z = s; -} - -static void normalize_rgb(float rgb[3]) -{ - const float max = max_fff(rgb[0], rgb[1], rgb[2]); - - if (max > 0.0f) { - mul_v3_fl(rgb, 1.0f / max); - } -} +/* ****************************** blackbody ******************************** */ -/* Color rendering of spectra, adapted from public domain code by John Walker, - * http://www.fourmilab.ch/ +/* Calculate color in range 800..12000 using an approximation + * a/x+bx+c for R and G and ((at + b)t + c)t + d) for B + * Max absolute error for RGB is (0.00095, 0.00077, 0.00057), + * which is enough to get the same 8 bit/channel color. */ -static void spectrum_to_xyz(float temperature, float xyz[3]) -{ - int i; - float lambda, x = 0.0f, y = 0.0f, z = 0.0f, xyz_sum; - - /* CIE colour matching functions xBar, yBar, and zBar for wavelengths from - * 380 through 780 nanometers, every 5 nanometers. - * For a wavelength lambda in this range: - * - * cie_colour_match[(lambda - 380) / 5][0] = xBar - * cie_colour_match[(lambda - 380) / 5][1] = yBar - * cie_colour_match[(lambda - 380) / 5][2] = zBar - */ - - const float cie_colour_match[81][3] = { - {0.0014f, 0.0000f, 0.0065f}, {0.0022f, 0.0001f, 0.0105f}, {0.0042f, 0.0001f, 0.0201f}, - {0.0076f, 0.0002f, 0.0362f}, {0.0143f, 0.0004f, 0.0679f}, {0.0232f, 0.0006f, 0.1102f}, - {0.0435f, 0.0012f, 0.2074f}, {0.0776f, 0.0022f, 0.3713f}, {0.1344f, 0.0040f, 0.6456f}, - {0.2148f, 0.0073f, 1.0391f}, {0.2839f, 0.0116f, 1.3856f}, {0.3285f, 0.0168f, 1.6230f}, - {0.3483f, 0.0230f, 1.7471f}, {0.3481f, 0.0298f, 1.7826f}, {0.3362f, 0.0380f, 1.7721f}, - {0.3187f, 0.0480f, 1.7441f}, {0.2908f, 0.0600f, 1.6692f}, {0.2511f, 0.0739f, 1.5281f}, - {0.1954f, 0.0910f, 1.2876f}, {0.1421f, 0.1126f, 1.0419f}, {0.0956f, 0.1390f, 0.8130f}, - {0.0580f, 0.1693f, 0.6162f}, {0.0320f, 0.2080f, 0.4652f}, {0.0147f, 0.2586f, 0.3533f}, - {0.0049f, 0.3230f, 0.2720f}, {0.0024f, 0.4073f, 0.2123f}, {0.0093f, 0.5030f, 0.1582f}, - {0.0291f, 0.6082f, 0.1117f}, {0.0633f, 0.7100f, 0.0782f}, {0.1096f, 0.7932f, 0.0573f}, - {0.1655f, 0.8620f, 0.0422f}, {0.2257f, 0.9149f, 0.0298f}, {0.2904f, 0.9540f, 0.0203f}, - {0.3597f, 0.9803f, 0.0134f}, {0.4334f, 0.9950f, 0.0087f}, {0.5121f, 1.0000f, 0.0057f}, - {0.5945f, 0.9950f, 0.0039f}, {0.6784f, 0.9786f, 0.0027f}, {0.7621f, 0.9520f, 0.0021f}, - {0.8425f, 0.9154f, 0.0018f}, {0.9163f, 0.8700f, 0.0017f}, {0.9786f, 0.8163f, 0.0014f}, - {1.0263f, 0.7570f, 0.0011f}, {1.0567f, 0.6949f, 0.0010f}, {1.0622f, 0.6310f, 0.0008f}, - {1.0456f, 0.5668f, 0.0006f}, {1.0026f, 0.5030f, 0.0003f}, {0.9384f, 0.4412f, 0.0002f}, - {0.8544f, 0.3810f, 0.0002f}, {0.7514f, 0.3210f, 0.0001f}, {0.6424f, 0.2650f, 0.0000f}, - {0.5419f, 0.2170f, 0.0000f}, {0.4479f, 0.1750f, 0.0000f}, {0.3608f, 0.1382f, 0.0000f}, - {0.2835f, 0.1070f, 0.0000f}, {0.2187f, 0.0816f, 0.0000f}, {0.1649f, 0.0610f, 0.0000f}, - {0.1212f, 0.0446f, 0.0000f}, {0.0874f, 0.0320f, 0.0000f}, {0.0636f, 0.0232f, 0.0000f}, - {0.0468f, 0.0170f, 0.0000f}, {0.0329f, 0.0119f, 0.0000f}, {0.0227f, 0.0082f, 0.0000f}, - {0.0158f, 0.0057f, 0.0000f}, {0.0114f, 0.0041f, 0.0000f}, {0.0081f, 0.0029f, 0.0000f}, - {0.0058f, 0.0021f, 0.0000f}, {0.0041f, 0.0015f, 0.0000f}, {0.0029f, 0.0010f, 0.0000f}, - {0.0020f, 0.0007f, 0.0000f}, {0.0014f, 0.0005f, 0.0000f}, {0.0010f, 0.0004f, 0.0000f}, - {0.0007f, 0.0002f, 0.0000f}, {0.0005f, 0.0002f, 0.0000f}, {0.0003f, 0.0001f, 0.0000f}, - {0.0002f, 0.0001f, 0.0000f}, {0.0002f, 0.0001f, 0.0000f}, {0.0001f, 0.0000f, 0.0000f}, - {0.0001f, 0.0000f, 0.0000f}, {0.0001f, 0.0000f, 0.0000f}, {0.0000f, 0.0000f, 0.0000f} - }; - - for (i = 0, lambda = 380.0f; lambda < 780.1f; i++, lambda += 5.0f) { - /* wavelength in meter */ - const float wlm = lambda * 1e-9f; - const float Me = (3.74183e-16f * powf(wlm, -5.0f)) / (expf(1.4388e-2f / (wlm * temperature)) - 1.0f); - - x += Me * cie_colour_match[i][0]; - y += Me * cie_colour_match[i][1]; - z += Me * cie_colour_match[i][2]; +static const float blackbody_table_r[6][3] = { + { 2.52432244e+03f, -1.06185848e-03f, 3.11067539e+00f }, + { 3.37763626e+03f, -4.34581697e-04f, 1.64843306e+00f }, + { 4.10671449e+03f, -8.61949938e-05f, 6.41423749e-01f }, + { 4.66849800e+03f, 2.85655028e-05f, 1.29075375e-01f }, + { 4.60124770e+03f, 2.89727618e-05f, 1.48001316e-01f }, + { 3.78765709e+03f, 9.36026367e-06f, 3.98995841e-01f }, +}; + +static const float blackbody_table_g[6][3] = { + { -7.50343014e+02f, 3.15679613e-04f, 4.73464526e-01f }, + { -1.00402363e+03f, 1.29189794e-04f, 9.08181524e-01f }, + { -1.22075471e+03f, 2.56245413e-05f, 1.20753416e+00f }, + { -1.42546105e+03f, -4.01730887e-05f, 1.44002695e+00f }, + { -1.18134453e+03f, -2.18913373e-05f, 1.30656109e+00f }, + { -5.00279505e+02f, -4.59745390e-06f, 1.09090465e+00f }, +}; + +static const float blackbody_table_b[6][4] = { + { 0.0f, 0.0f, 0.0f, 0.0f }, + { 0.0f, 0.0f, 0.0f, 0.0f }, + { 0.0f, 0.0f, 0.0f, 0.0f }, + { -2.02524603e-11f, 1.79435860e-07f, -2.60561875e-04f, -1.41761141e-02f }, + { -2.22463426e-13f, -1.55078698e-08f, 3.81675160e-04f, -7.30646033e-01f }, + { 6.72595954e-13f, -2.73059993e-08f, 4.24068546e-04f, -7.52204323e-01f }, +}; + +static void blackbody_temperature_to_rgb(float rgb[3], float t) +{ + if (t >= 12000.0f) { + rgb[0] = 0.826270103f; + rgb[1] = 0.994478524f; + rgb[2] = 1.56626022f; + } + else if (t < 965.0f) { + rgb[0] = 4.70366907f; + rgb[1] = 0.0f; + rgb[2] = 0.0f; + } + else { + int i = (t >= 6365.0f) ? 5 : + (t >= 3315.0f) ? 4 : + (t >= 1902.0f) ? 3 : + (t >= 1449.0f) ? 2 : + (t >= 1167.0f) ? 1 : 0; + + const float *r = blackbody_table_r[i]; + const float *g = blackbody_table_g[i]; + const float *b = blackbody_table_b[i]; + + const float t_inv = 1.0f / t; + rgb[0] = r[0] * t_inv + r[1] * t + r[2]; + rgb[1] = g[0] * t_inv + g[1] * t + g[2]; + rgb[2] = ((b[0] * t + b[1]) * t + b[2]) * t + b[3]; } - - xyz_sum = (x + y + z); - - xyz[0] = x / xyz_sum; - xyz[1] = y / xyz_sum; - xyz[2] = z / xyz_sum; } void blackbody_temperature_to_rgb_table(float *r_table, int width, float min, float max) { - int i, j = 0, dj = 1; - float rgb[3], xyz[3], lms[3], lms_w[3]; - float bb_temp; - - if (min < max) { - SWAP(float, min, max); - j = width - 1; - dj = -1; - } + for (int i = 0; i < width; i++) { + float temperature = min + (max - min) / (float)width * (float)i; - for (i = 0; i < width; i++, j += dj) { - bb_temp = min + (max - min) / (float)width * (float)i; + float rgb[3]; + blackbody_temperature_to_rgb(rgb, temperature); - /* integrate blackbody radiation spectrum to XYZ */ - spectrum_to_xyz(bb_temp, xyz); - - /* normalize highest temperature to white (in LMS system) */ - xyz_to_lms(xyz[0], xyz[1], xyz[2], &lms[0], &lms[1], &lms[2]); - - if (i == 0) { - lms_w[0] = 1.0f / lms[0]; - lms_w[1] = 1.0f / lms[1]; - lms_w[2] = 1.0f / lms[2]; - } - - mul_v3_v3(lms, lms_w); - - lms_to_xyz(lms[0], lms[1], lms[2], &xyz[0], &xyz[1], &xyz[2]); - - /* convert to RGB */ - xyz_to_rgb(xyz[0], xyz[1], xyz[2], &rgb[0], &rgb[1], &rgb[2], BLI_XYZ_CIE); - constrain_rgb(&rgb[0], &rgb[1], &rgb[2]); - normalize_rgb(rgb); - - copy_v3_v3(&r_table[(j << 2)], rgb); - - if (rgb[2] > 0.1f) - r_table[(j << 2) + 3] = rgb[2]; - else - r_table[(j << 2) + 3] = 0.0f; + copy_v3_v3(&r_table[i * 4], rgb); + r_table[i * 4 + 3] = 0.0f; } } diff --git a/source/blender/blenlib/intern/math_geom.c b/source/blender/blenlib/intern/math_geom.c index 61b8f6819cc..9425931133a 100644 --- a/source/blender/blenlib/intern/math_geom.c +++ b/source/blender/blenlib/intern/math_geom.c @@ -570,6 +570,8 @@ float dist_squared_to_ray_v3( *r_depth = dot_v3v3(dvec, ray_direction); return len_squared_v3(dvec) - SQUARE(*r_depth); } + + /** * Find the closest point in a seg to a ray and return the distance squared. * \param r_point: Is the point on segment closest to ray (or to ray_origin if the ray and the segment are parallel). @@ -580,43 +582,68 @@ float dist_squared_ray_to_seg_v3( const float v0[3], const float v1[3], float r_point[3], float *r_depth) { - float a[3], t[3], n[3], lambda; - sub_v3_v3v3(a, v1, v0); - sub_v3_v3v3(t, v0, ray_origin); - cross_v3_v3v3(n, a, ray_direction); - const float nlen = len_squared_v3(n); - - /* if (nlen == 0.0f) the lines are parallel, - * has no nearest point, only distance squared.*/ - if (nlen == 0.0f) { - /* Calculate the distance to the point v0 then */ - copy_v3_v3(r_point, v0); - *r_depth = dot_v3v3(t, ray_direction); - } - else { - float c[3], cray[3]; - sub_v3_v3v3(c, n, t); - cross_v3_v3v3(cray, c, ray_direction); - lambda = dot_v3v3(cray, n) / nlen; - if (lambda <= 0) { + float lambda, depth; + if (isect_ray_seg_v3( + ray_origin, ray_direction, v0, v1, &lambda)) + { + if (lambda <= 0.0f) { copy_v3_v3(r_point, v0); - - *r_depth = dot_v3v3(t, ray_direction); } - else if (lambda >= 1) { + else if (lambda >= 1.0f) { copy_v3_v3(r_point, v1); - - sub_v3_v3v3(t, v1, ray_origin); - *r_depth = dot_v3v3(t, ray_direction); } else { - madd_v3_v3v3fl(r_point, v0, a, lambda); - - sub_v3_v3v3(t, r_point, ray_origin); - *r_depth = dot_v3v3(t, ray_direction); + interp_v3_v3v3(r_point, v0, v1, lambda); } } - return len_squared_v3(t) - SQUARE(*r_depth); + else { + /* has no nearest point, only distance squared. */ + /* Calculate the distance to the point v0 then */ + copy_v3_v3(r_point, v0); + } + + float dvec[3]; + sub_v3_v3v3(dvec, r_point, ray_origin); + depth = dot_v3v3(dvec, ray_direction); + + if (r_depth) { + *r_depth = depth; + } + + return len_squared_v3(dvec) - SQUARE(depth); +} + + +/* Returns the coordinates of the nearest vertex and + * the farthest vertex from a plane (or normal). */ +void aabb_get_near_far_from_plane( + const float plane_no[3], const float bbmin[3], const float bbmax[3], + float bb_near[3], float bb_afar[3]) +{ + if (plane_no[0] < 0.0f) { + bb_near[0] = bbmax[0]; + bb_afar[0] = bbmin[0]; + } + else { + bb_near[0] = bbmin[0]; + bb_afar[0] = bbmax[0]; + } + if (plane_no[1] < 0.0f) { + bb_near[1] = bbmax[1]; + bb_afar[1] = bbmin[1]; + } + else { + bb_near[1] = bbmin[1]; + bb_afar[1] = bbmax[1]; + } + if (plane_no[2] < 0.0f) { + bb_near[2] = bbmax[2]; + bb_afar[2] = bbmin[2]; + } + else { + bb_near[2] = bbmin[2]; + bb_afar[2] = bbmax[2]; + } } /* -------------------------------------------------------------------- */ @@ -634,7 +661,6 @@ void dist_squared_ray_to_aabb_v3_precalc( neasrest_precalc->ray_inv_dir[i] = (neasrest_precalc->ray_direction[i] != 0.0f) ? (1.0f / neasrest_precalc->ray_direction[i]) : FLT_MAX; - neasrest_precalc->sign[i] = (neasrest_precalc->ray_inv_dir[i] < 0.0f); } } @@ -648,30 +674,8 @@ float dist_squared_ray_to_aabb_v3( { // bool r_axis_closest[3]; float local_bvmin[3], local_bvmax[3]; - if (data->sign[0]) { - local_bvmin[0] = bb_max[0]; - local_bvmax[0] = bb_min[0]; - } - else { - local_bvmin[0] = bb_min[0]; - local_bvmax[0] = bb_max[0]; - } - if (data->sign[1]) { - local_bvmin[1] = bb_max[1]; - local_bvmax[1] = bb_min[1]; - } - else { - local_bvmin[1] = bb_min[1]; - local_bvmax[1] = bb_max[1]; - } - if (data->sign[2]) { - local_bvmin[2] = bb_max[2]; - local_bvmax[2] = bb_min[2]; - } - else { - local_bvmin[2] = bb_min[2]; - local_bvmax[2] = bb_max[2]; - } + aabb_get_near_far_from_plane( + data->ray_direction, bb_min, bb_max, local_bvmin, local_bvmax); const float tmin[3] = { (local_bvmin[0] - data->ray_origin[0]) * data->ray_inv_dir[0], @@ -693,38 +697,38 @@ float dist_squared_ray_to_aabb_v3( rtmax = tmax[0]; va[0] = vb[0] = local_bvmax[0]; main_axis = 3; - // r_axis_closest[0] = data->sign[0]; + // r_axis_closest[0] = neasrest_precalc->ray_direction[0] < 0.0f; } else if ((tmax[1] <= tmax[0]) && (tmax[1] <= tmax[2])) { rtmax = tmax[1]; va[1] = vb[1] = local_bvmax[1]; main_axis = 2; - // r_axis_closest[1] = data->sign[1]; + // r_axis_closest[1] = neasrest_precalc->ray_direction[1] < 0.0f; } else { rtmax = tmax[2]; va[2] = vb[2] = local_bvmax[2]; main_axis = 1; - // r_axis_closest[2] = data->sign[2]; + // r_axis_closest[2] = neasrest_precalc->ray_direction[2] < 0.0f; } if ((tmin[0] >= tmin[1]) && (tmin[0] >= tmin[2])) { rtmin = tmin[0]; va[0] = vb[0] = local_bvmin[0]; main_axis -= 3; - // r_axis_closest[0] = !data->sign[0]; + // r_axis_closest[0] = neasrest_precalc->ray_direction[0] >= 0.0f; } else if ((tmin[1] >= tmin[0]) && (tmin[1] >= tmin[2])) { rtmin = tmin[1]; va[1] = vb[1] = local_bvmin[1]; main_axis -= 1; - // r_axis_closest[1] = !data->sign[1]; + // r_axis_closest[1] = neasrest_precalc->ray_direction[1] >= 0.0f; } else { rtmin = tmin[2]; va[2] = vb[2] = local_bvmin[2]; main_axis -= 2; - // r_axis_closest[2] = !data->sign[2]; + // r_axis_closest[2] = neasrest_precalc->ray_direction[2] >= 0.0f; } if (main_axis < 0) { main_axis += 3; @@ -739,14 +743,14 @@ float dist_squared_ray_to_aabb_v3( return 0.0f; } - if (data->sign[main_axis]) { - va[main_axis] = local_bvmax[main_axis]; - vb[main_axis] = local_bvmin[main_axis]; - } - else { + if (data->ray_direction[main_axis] >= 0.0f) { va[main_axis] = local_bvmin[main_axis]; vb[main_axis] = local_bvmax[main_axis]; } + else { + va[main_axis] = local_bvmax[main_axis]; + vb[main_axis] = local_bvmin[main_axis]; + } return dist_squared_ray_to_seg_v3( data->ray_origin, data->ray_direction, va, vb, @@ -765,6 +769,214 @@ float dist_squared_ray_to_aabb_v3_simple( /** \} */ +/* -------------------------------------------------------------------- */ +/** \name dist_squared_to_projected_aabb and helpers +* \{ */ + +/** + * \param projmat: Projection Matrix (usually perspective + * matrix multiplied by object matrix). + */ +void dist_squared_to_projected_aabb_precalc( + struct DistProjectedAABBPrecalc *precalc, + const float projmat[4][4], const float winsize[2], const float mval[2]) +{ + float win_half[2], relative_mval[2], px[4], py[4]; + + mul_v2_v2fl(win_half, winsize, 0.5f); + sub_v2_v2v2(precalc->mval, mval, win_half); + + relative_mval[0] = precalc->mval[0] / win_half[0]; + relative_mval[1] = precalc->mval[1] / win_half[1]; + + copy_m4_m4(precalc->pmat, projmat); + for (int i = 0; i < 4; i++) { + px[i] = precalc->pmat[i][0] - precalc->pmat[i][3] * relative_mval[0]; + py[i] = precalc->pmat[i][1] - precalc->pmat[i][3] * relative_mval[1]; + + precalc->pmat[i][0] *= win_half[0]; + precalc->pmat[i][1] *= win_half[1]; + } +#if 0 + float projmat_trans[4][4]; + transpose_m4_m4(projmat_trans, projmat); + if (!isect_plane_plane_plane_v3( + projmat_trans[0], projmat_trans[1], projmat_trans[3], + precalc->ray_origin)) + { + /* Orthographic projection. */ + isect_plane_plane_v3( + px, py, + precalc->ray_origin, + precalc->ray_direction); + } + else { + /* Perspective projection. */ + cross_v3_v3v3(precalc->ray_direction, py, px); + //normalize_v3(precalc->ray_direction); + } +#else + if (!isect_plane_plane_v3( + px, py, + precalc->ray_origin, + precalc->ray_direction)) + { + /* Matrix with weird coplanar planes. Undetermined origin.*/ + zero_v3(precalc->ray_origin); + precalc->ray_direction[0] = precalc->pmat[0][3]; + precalc->ray_direction[1] = precalc->pmat[1][3]; + precalc->ray_direction[2] = precalc->pmat[2][3]; + } +#endif + + for (int i = 0; i < 3; i++) { + precalc->ray_inv_dir[i] = + (precalc->ray_direction[i] != 0.0f) ? + (1.0f / precalc->ray_direction[i]) : FLT_MAX; + } +} + +/* Returns the distance from a 2d coordinate to a BoundBox (Projected) */ +float dist_squared_to_projected_aabb( + struct DistProjectedAABBPrecalc *data, + const float bbmin[3], const float bbmax[3], + bool r_axis_closest[3]) +{ + float local_bvmin[3], local_bvmax[3]; + aabb_get_near_far_from_plane( + data->ray_direction, bbmin, bbmax, local_bvmin, local_bvmax); + + const float tmin[3] = { + (local_bvmin[0] - data->ray_origin[0]) * data->ray_inv_dir[0], + (local_bvmin[1] - data->ray_origin[1]) * data->ray_inv_dir[1], + (local_bvmin[2] - data->ray_origin[2]) * data->ray_inv_dir[2], + }; + const float tmax[3] = { + (local_bvmax[0] - data->ray_origin[0]) * data->ray_inv_dir[0], + (local_bvmax[1] - data->ray_origin[1]) * data->ray_inv_dir[1], + (local_bvmax[2] - data->ray_origin[2]) * data->ray_inv_dir[2], + }; + /* `va` and `vb` are the coordinates of the AABB edge closest to the ray */ + float va[3], vb[3]; + /* `rtmin` and `rtmax` are the minimum and maximum distances of the ray hits on the AABB */ + float rtmin, rtmax; + int main_axis; + + if ((tmax[0] <= tmax[1]) && (tmax[0] <= tmax[2])) { + rtmax = tmax[0]; + va[0] = vb[0] = local_bvmax[0]; + main_axis = 3; + r_axis_closest[0] = data->ray_direction[0] < 0.0f; + } + else if ((tmax[1] <= tmax[0]) && (tmax[1] <= tmax[2])) { + rtmax = tmax[1]; + va[1] = vb[1] = local_bvmax[1]; + main_axis = 2; + r_axis_closest[1] = data->ray_direction[1] < 0.0f; + } + else { + rtmax = tmax[2]; + va[2] = vb[2] = local_bvmax[2]; + main_axis = 1; + r_axis_closest[2] = data->ray_direction[2] < 0.0f; + } + + if ((tmin[0] >= tmin[1]) && (tmin[0] >= tmin[2])) { + rtmin = tmin[0]; + va[0] = vb[0] = local_bvmin[0]; + main_axis -= 3; + r_axis_closest[0] = data->ray_direction[0] >= 0.0f; + } + else if ((tmin[1] >= tmin[0]) && (tmin[1] >= tmin[2])) { + rtmin = tmin[1]; + va[1] = vb[1] = local_bvmin[1]; + main_axis -= 1; + r_axis_closest[1] = data->ray_direction[1] >= 0.0f; + } + else { + rtmin = tmin[2]; + va[2] = vb[2] = local_bvmin[2]; + main_axis -= 2; + r_axis_closest[2] = data->ray_direction[2] >= 0.0f; + } + if (main_axis < 0) { + main_axis += 3; + } + + /* if rtmin <= rtmax, ray intersect `AABB` */ + if (rtmin <= rtmax) { + return 0; + } + + if (data->ray_direction[main_axis] >= 0.0f) { + va[main_axis] = local_bvmin[main_axis]; + vb[main_axis] = local_bvmax[main_axis]; + } + else { + va[main_axis] = local_bvmax[main_axis]; + vb[main_axis] = local_bvmin[main_axis]; + } + float scale = fabsf(local_bvmax[main_axis] - local_bvmin[main_axis]); + + float va2d[2] = { + (dot_m4_v3_row_x(data->pmat, va) + data->pmat[3][0]), + (dot_m4_v3_row_y(data->pmat, va) + data->pmat[3][1]), + }; + float vb2d[2] = { + (va2d[0] + data->pmat[main_axis][0] * scale), + (va2d[1] + data->pmat[main_axis][1] * scale), + }; + + float w_a = mul_project_m4_v3_zfac(data->pmat, va); + if (w_a != 1.0f) { + /* Perspective Projection. */ + float w_b = w_a + data->pmat[main_axis][3] * scale; + va2d[0] /= w_a; + va2d[1] /= w_a; + vb2d[0] /= w_b; + vb2d[1] /= w_b; + } + + float dvec[2], edge[2], lambda, rdist_sq; + sub_v2_v2v2(dvec, data->mval, va2d); + sub_v2_v2v2(edge, vb2d, va2d); + lambda = dot_v2v2(dvec, edge); + if (lambda != 0.0f) { + lambda /= len_squared_v2(edge); + if (lambda <= 0.0f) { + rdist_sq = len_squared_v2v2(data->mval, va2d); + r_axis_closest[main_axis] = true; + } + else if (lambda >= 1.0f) { + rdist_sq = len_squared_v2v2(data->mval, vb2d); + r_axis_closest[main_axis] = false; + } + else { + madd_v2_v2fl(va2d, edge, lambda); + rdist_sq = len_squared_v2v2(data->mval, va2d); + r_axis_closest[main_axis] = lambda < 0.5f; + } + } + else { + rdist_sq = len_squared_v2v2(data->mval, va2d); + } + + return rdist_sq; +} + +float dist_squared_to_projected_aabb_simple( + const float projmat[4][4], const float winsize[2], const float mval[2], + const float bbmin[3], const float bbmax[3]) +{ + struct DistProjectedAABBPrecalc data; + dist_squared_to_projected_aabb_precalc(&data, projmat, winsize, mval); + + bool dummy[3] = {true, true, true}; + return dist_squared_to_projected_aabb(&data, bbmin, bbmax, dummy); +} +/** \} */ + + /* Adapted from "Real-Time Collision Detection" by Christer Ericson, * published by Morgan Kaufmann Publishers, copyright 2005 Elsevier Inc. * @@ -1769,6 +1981,33 @@ bool isect_ray_seg_v2( return false; } + +bool isect_ray_seg_v3( + const float ray_origin[3], const float ray_direction[3], + const float v0[3], const float v1[3], + float *r_lambda) +{ + float a[3], t[3], n[3]; + sub_v3_v3v3(a, v1, v0); + sub_v3_v3v3(t, v0, ray_origin); + cross_v3_v3v3(n, a, ray_direction); + const float nlen = len_squared_v3(n); + + if (nlen == 0.0f) { + /* the lines are parallel.*/ + return false; + } + + float c[3], cray[3]; + sub_v3_v3v3(c, n, t); + cross_v3_v3v3(cray, c, ray_direction); + + *r_lambda = dot_v3v3(cray, n) / nlen; + + return true; +} + + /** * Check if a point is behind all planes. */ @@ -1786,6 +2025,23 @@ bool isect_point_planes_v3(float (*planes)[4], int totplane, const float p[3]) } /** + * Check if a point is in front all planes. + * Same as isect_point_planes_v3 but with planes facing the opposite direction. + */ +bool isect_point_planes_v3_negated( + const float(*planes)[4], const int totplane, const float p[3]) +{ + for (int i = 0; i < totplane; i++) { + if (plane_point_side_v3(planes[i], p) <= 0.0f) { + return false; + } + } + + return true; +} + + +/** * Intersect line/plane. * * \param r_isect_co The intersection point. @@ -2043,6 +2299,38 @@ static bool getLowestRoot(const float a, const float b, const float c, const flo return false; } + +/** + * Checks status of an AABB in relation to a list of planes. + * + * \returns intersection type: + * - ISECT_AABB_PLANE_BEHIND_ONE (0): AABB is completely behind at least 1 plane; + * - ISECT_AABB_PLANE_CROSS_ANY (1): AABB intersects at least 1 plane; + * - ISECT_AABB_PLANE_IN_FRONT_ALL (2): AABB is completely in front of all planes; + */ +int isect_aabb_planes_v3( + const float (*planes)[4], const int totplane, + const float bbmin[3], const float bbmax[3]) +{ + int ret = ISECT_AABB_PLANE_IN_FRONT_ALL; + + float bb_near[3], bb_far[3]; + for (int i = 0; i < totplane; i++) { + aabb_get_near_far_from_plane(planes[i], bbmin, bbmax, bb_near, bb_far); + + if (plane_point_side_v3(planes[i], bb_far) < 0.0f) { + return ISECT_AABB_PLANE_BEHIND_ANY; + } + else if ((ret != ISECT_AABB_PLANE_CROSS_ANY) && + (plane_point_side_v3(planes[i], bb_near) < 0.0f)) + { + ret = ISECT_AABB_PLANE_CROSS_ANY; + } + } + + return ret; +} + bool isect_sweeping_sphere_tri_v3(const float p1[3], const float p2[3], const float radius, const float v0[3], const float v1[3], const float v2[3], float *r_lambda, float ipoint[3]) diff --git a/source/blender/blenlib/intern/math_matrix.c b/source/blender/blenlib/intern/math_matrix.c index 311d963f64d..48bc1a2928b 100644 --- a/source/blender/blenlib/intern/math_matrix.c +++ b/source/blender/blenlib/intern/math_matrix.c @@ -33,6 +33,8 @@ #include "BLI_strict_flags.h" +#include "eigen_capi.h" + /********************************* Init **************************************/ void zero_m2(float m[2][2]) @@ -74,23 +76,23 @@ void unit_m4(float m[4][4]) m[3][0] = m[3][1] = m[3][2] = 0.0f; } -void copy_m2_m2(float m1[2][2], float m2[2][2]) +void copy_m2_m2(float m1[2][2], const float m2[2][2]) { memcpy(m1, m2, sizeof(float[2][2])); } -void copy_m3_m3(float m1[3][3], float m2[3][3]) +void copy_m3_m3(float m1[3][3], const float m2[3][3]) { /* destination comes first: */ memcpy(m1, m2, sizeof(float[3][3])); } -void copy_m4_m4(float m1[4][4], float m2[4][4]) +void copy_m4_m4(float m1[4][4], const float m2[4][4]) { memcpy(m1, m2, sizeof(float[4][4])); } -void copy_m3_m4(float m1[3][3], float m2[4][4]) +void copy_m3_m4(float m1[3][3], const float m2[4][4]) { m1[0][0] = m2[0][0]; m1[0][1] = m2[0][1]; @@ -105,7 +107,7 @@ void copy_m3_m4(float m1[3][3], float m2[4][4]) m1[2][2] = m2[2][2]; } -void copy_m4_m3(float m1[4][4], float m2[3][3]) /* no clear */ +void copy_m4_m3(float m1[4][4], const float m2[3][3]) /* no clear */ { m1[0][0] = m2[0][0]; m1[0][1] = m2[0][1]; @@ -131,7 +133,7 @@ void copy_m4_m3(float m1[4][4], float m2[3][3]) /* no clear */ } -void copy_m3_m3d(float R[3][3], double A[3][3]) +void copy_m3_m3d(float R[3][3], const double A[3][3]) { /* Keep it stupid simple for better data flow in CPU. */ R[0][0] = (float)A[0][0]; @@ -177,64 +179,127 @@ void swap_m4m4(float m1[4][4], float m2[4][4]) /******************************** Arithmetic *********************************/ -void mul_m4_m4m4(float m1[4][4], float m3_[4][4], float m2_[4][4]) +void mul_m4_m4m4(float R[4][4], const float A[4][4], const float B[4][4]) { - float m2[4][4], m3[4][4]; + if (A == R) + mul_m4_m4_post(R, B); + else if (B == R) + mul_m4_m4_pre(R, A); + else + mul_m4_m4m4_uniq(R, A, B); +} - /* copy so it works when m1 is the same pointer as m2 or m3 */ - copy_m4_m4(m2, m2_); - copy_m4_m4(m3, m3_); +void mul_m4_m4m4_uniq(float R[4][4], const float A[4][4], const float B[4][4]) +{ + BLI_assert(R != A && R != B); - /* matrix product: m1[j][k] = m2[j][i].m3[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]; + /* matrix product: R[j][k] = A[j][i] . B[i][k] */ +#ifdef __SSE2__ + __m128 A0 = _mm_loadu_ps(A[0]); + __m128 A1 = _mm_loadu_ps(A[1]); + __m128 A2 = _mm_loadu_ps(A[2]); + __m128 A3 = _mm_loadu_ps(A[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]; + for (int i = 0; i < 4; i++) { + __m128 B0 = _mm_set1_ps(B[i][0]); + __m128 B1 = _mm_set1_ps(B[i][1]); + __m128 B2 = _mm_set1_ps(B[i][2]); + __m128 B3 = _mm_set1_ps(B[i][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]; + __m128 sum = _mm_add_ps( + _mm_add_ps(_mm_mul_ps(B0, A0), _mm_mul_ps(B1, A1)), + _mm_add_ps(_mm_mul_ps(B2, A2), _mm_mul_ps(B3, A3))); - 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]; + _mm_storeu_ps(R[i], sum); + } +#else + R[0][0] = B[0][0] * A[0][0] + B[0][1] * A[1][0] + B[0][2] * A[2][0] + B[0][3] * A[3][0]; + R[0][1] = B[0][0] * A[0][1] + B[0][1] * A[1][1] + B[0][2] * A[2][1] + B[0][3] * A[3][1]; + R[0][2] = B[0][0] * A[0][2] + B[0][1] * A[1][2] + B[0][2] * A[2][2] + B[0][3] * A[3][2]; + R[0][3] = B[0][0] * A[0][3] + B[0][1] * A[1][3] + B[0][2] * A[2][3] + B[0][3] * A[3][3]; + + R[1][0] = B[1][0] * A[0][0] + B[1][1] * A[1][0] + B[1][2] * A[2][0] + B[1][3] * A[3][0]; + R[1][1] = B[1][0] * A[0][1] + B[1][1] * A[1][1] + B[1][2] * A[2][1] + B[1][3] * A[3][1]; + R[1][2] = B[1][0] * A[0][2] + B[1][1] * A[1][2] + B[1][2] * A[2][2] + B[1][3] * A[3][2]; + R[1][3] = B[1][0] * A[0][3] + B[1][1] * A[1][3] + B[1][2] * A[2][3] + B[1][3] * A[3][3]; + + R[2][0] = B[2][0] * A[0][0] + B[2][1] * A[1][0] + B[2][2] * A[2][0] + B[2][3] * A[3][0]; + R[2][1] = B[2][0] * A[0][1] + B[2][1] * A[1][1] + B[2][2] * A[2][1] + B[2][3] * A[3][1]; + R[2][2] = B[2][0] * A[0][2] + B[2][1] * A[1][2] + B[2][2] * A[2][2] + B[2][3] * A[3][2]; + R[2][3] = B[2][0] * A[0][3] + B[2][1] * A[1][3] + B[2][2] * A[2][3] + B[2][3] * A[3][3]; + + R[3][0] = B[3][0] * A[0][0] + B[3][1] * A[1][0] + B[3][2] * A[2][0] + B[3][3] * A[3][0]; + R[3][1] = B[3][0] * A[0][1] + B[3][1] * A[1][1] + B[3][2] * A[2][1] + B[3][3] * A[3][1]; + R[3][2] = B[3][0] * A[0][2] + B[3][1] * A[1][2] + B[3][2] * A[2][2] + B[3][3] * A[3][2]; + R[3][3] = B[3][0] * A[0][3] + B[3][1] * A[1][3] + B[3][2] * A[2][3] + B[3][3] * A[3][3]; +#endif +} +void mul_m4_m4_pre(float R[4][4], const float A[4][4]) +{ + BLI_assert(A != R); + float B[4][4]; + copy_m4_m4(B, R); + mul_m4_m4m4_uniq(R, A, B); } -void mul_m3_m3m3(float m1[3][3], float m3_[3][3], float m2_[3][3]) +void mul_m4_m4_post(float R[4][4], const float B[4][4]) { - float m2[3][3], m3[3][3]; + BLI_assert(B != R); + float A[4][4]; + copy_m4_m4(A, R); + mul_m4_m4m4_uniq(R, A, B); +} - /* copy so it works when m1 is the same pointer as m2 or m3 */ - copy_m3_m3(m2, m2_); - copy_m3_m3(m3, m3_); +void mul_m3_m3m3(float R[3][3], const float A[3][3], const float B[3][3]) +{ + if (A == R) + mul_m3_m3_post(R, B); + else if (B == R) + mul_m3_m3_pre(R, A); + else + mul_m3_m3m3_uniq(R, A, B); +} - /* 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]; +void mul_m3_m3_pre(float R[3][3], const float A[3][3]) +{ + BLI_assert(A != R); + float B[3][3]; + copy_m3_m3(B, R); + mul_m3_m3m3_uniq(R, A, B); +} - 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]; +void mul_m3_m3_post(float R[3][3], const float B[3][3]) +{ + BLI_assert(B != R); + float A[3][3]; + copy_m3_m3(A, R); + mul_m3_m3m3_uniq(R, A, B); +} - 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]; +void mul_m3_m3m3_uniq(float R[3][3], const float A[3][3], const float B[3][3]) +{ + BLI_assert(R != A && R != B); + + R[0][0] = B[0][0] * A[0][0] + B[0][1] * A[1][0] + B[0][2] * A[2][0]; + R[0][1] = B[0][0] * A[0][1] + B[0][1] * A[1][1] + B[0][2] * A[2][1]; + R[0][2] = B[0][0] * A[0][2] + B[0][1] * A[1][2] + B[0][2] * A[2][2]; + + R[1][0] = B[1][0] * A[0][0] + B[1][1] * A[1][0] + B[1][2] * A[2][0]; + R[1][1] = B[1][0] * A[0][1] + B[1][1] * A[1][1] + B[1][2] * A[2][1]; + R[1][2] = B[1][0] * A[0][2] + B[1][1] * A[1][2] + B[1][2] * A[2][2]; + + R[2][0] = B[2][0] * A[0][0] + B[2][1] * A[1][0] + B[2][2] * A[2][0]; + R[2][1] = B[2][0] * A[0][1] + B[2][1] * A[1][1] + B[2][2] * A[2][1]; + R[2][2] = B[2][0] * A[0][2] + B[2][1] * A[1][2] + B[2][2] * A[2][2]; } -void mul_m4_m4m3(float m1[4][4], float m3_[4][4], float m2_[3][3]) +void mul_m4_m4m3(float m1[4][4], const float m3_[4][4], const float m2_[3][3]) { float m2[3][3], m3[4][4]; /* copy so it works when m1 is the same pointer as m2 or m3 */ + /* TODO: avoid copying when matrices are different */ copy_m3_m3(m2, m2_); copy_m4_m4(m3, m3_); @@ -250,11 +315,12 @@ void mul_m4_m4m3(float m1[4][4], float m3_[4][4], float m2_[3][3]) } /* m1 = m2 * m3, ignore the elements on the 4th row/column of m3 */ -void mul_m3_m3m4(float m1[3][3], float m3_[4][4], float m2_[3][3]) +void mul_m3_m3m4(float m1[3][3], const float m3_[4][4], const float m2_[3][3]) { float m2[3][3], m3[4][4]; /* copy so it works when m1 is the same pointer as m2 or m3 */ + /* TODO: avoid copying when matrices are different */ copy_m3_m3(m2, m2_); copy_m4_m4(m3, m3_); @@ -272,11 +338,12 @@ void mul_m3_m3m4(float m1[3][3], float m3_[4][4], float m2_[3][3]) m1[2][2] = m2[2][0] * m3[0][2] + m2[2][1] * m3[1][2] + m2[2][2] * m3[2][2]; } -void mul_m4_m3m4(float m1[4][4], float m3_[3][3], float m2_[4][4]) +void mul_m4_m3m4(float m1[4][4], const float m3_[3][3], const float m2_[4][4]) { float m2[4][4], m3[3][3]; /* copy so it works when m1 is the same pointer as m2 or m3 */ + /* TODO: avoid copying when matrices are different */ copy_m4_m4(m2, m2_); copy_m3_m3(m3, m3_); @@ -296,20 +363,20 @@ void mul_m4_m3m4(float m1[4][4], float m3_[3][3], float m2_[4][4]) * \{ */ void _va_mul_m3_series_3( float r[3][3], - float m1[3][3], float m2[3][3]) + const float m1[3][3], const float m2[3][3]) { mul_m3_m3m3(r, m1, m2); } void _va_mul_m3_series_4( float r[3][3], - float m1[3][3], float m2[3][3], float m3[3][3]) + const float m1[3][3], const float m2[3][3], const float m3[3][3]) { mul_m3_m3m3(r, m1, m2); mul_m3_m3m3(r, r, m3); } void _va_mul_m3_series_5( float r[3][3], - float m1[3][3], float m2[3][3], float m3[3][3], float m4[3][3]) + const float m1[3][3], const float m2[3][3], const float m3[3][3], const float m4[3][3]) { mul_m3_m3m3(r, m1, m2); mul_m3_m3m3(r, r, m3); @@ -317,8 +384,8 @@ void _va_mul_m3_series_5( } void _va_mul_m3_series_6( float r[3][3], - float m1[3][3], float m2[3][3], float m3[3][3], float m4[3][3], - float m5[3][3]) + const float m1[3][3], const float m2[3][3], const float m3[3][3], const float m4[3][3], + const float m5[3][3]) { mul_m3_m3m3(r, m1, m2); mul_m3_m3m3(r, r, m3); @@ -327,8 +394,8 @@ void _va_mul_m3_series_6( } void _va_mul_m3_series_7( float r[3][3], - float m1[3][3], float m2[3][3], float m3[3][3], float m4[3][3], - float m5[3][3], float m6[3][3]) + const float m1[3][3], const float m2[3][3], const float m3[3][3], const float m4[3][3], + const float m5[3][3], const float m6[3][3]) { mul_m3_m3m3(r, m1, m2); mul_m3_m3m3(r, r, m3); @@ -338,8 +405,8 @@ void _va_mul_m3_series_7( } void _va_mul_m3_series_8( float r[3][3], - float m1[3][3], float m2[3][3], float m3[3][3], float m4[3][3], - float m5[3][3], float m6[3][3], float m7[3][3]) + const float m1[3][3], const float m2[3][3], const float m3[3][3], const float m4[3][3], + const float m5[3][3], const float m6[3][3], const float m7[3][3]) { mul_m3_m3m3(r, m1, m2); mul_m3_m3m3(r, r, m3); @@ -350,8 +417,8 @@ void _va_mul_m3_series_8( } void _va_mul_m3_series_9( float r[3][3], - float m1[3][3], float m2[3][3], float m3[3][3], float m4[3][3], - float m5[3][3], float m6[3][3], float m7[3][3], float m8[3][3]) + const float m1[3][3], const float m2[3][3], const float m3[3][3], const float m4[3][3], + const float m5[3][3], const float m6[3][3], const float m7[3][3], const float m8[3][3]) { mul_m3_m3m3(r, m1, m2); mul_m3_m3m3(r, r, m3); @@ -367,20 +434,20 @@ void _va_mul_m3_series_9( * \{ */ void _va_mul_m4_series_3( float r[4][4], - float m1[4][4], float m2[4][4]) + const float m1[4][4], const float m2[4][4]) { mul_m4_m4m4(r, m1, m2); } void _va_mul_m4_series_4( float r[4][4], - float m1[4][4], float m2[4][4], float m3[4][4]) + const float m1[4][4], const float m2[4][4], const float m3[4][4]) { mul_m4_m4m4(r, m1, m2); mul_m4_m4m4(r, r, m3); } void _va_mul_m4_series_5( float r[4][4], - float m1[4][4], float m2[4][4], float m3[4][4], float m4[4][4]) + const float m1[4][4], const float m2[4][4], const float m3[4][4], const float m4[4][4]) { mul_m4_m4m4(r, m1, m2); mul_m4_m4m4(r, r, m3); @@ -388,8 +455,8 @@ void _va_mul_m4_series_5( } void _va_mul_m4_series_6( float r[4][4], - float m1[4][4], float m2[4][4], float m3[4][4], float m4[4][4], - float m5[4][4]) + const float m1[4][4], const float m2[4][4], const float m3[4][4], const float m4[4][4], + const float m5[4][4]) { mul_m4_m4m4(r, m1, m2); mul_m4_m4m4(r, r, m3); @@ -398,8 +465,8 @@ void _va_mul_m4_series_6( } void _va_mul_m4_series_7( float r[4][4], - float m1[4][4], float m2[4][4], float m3[4][4], float m4[4][4], - float m5[4][4], float m6[4][4]) + const float m1[4][4], const float m2[4][4], const float m3[4][4], const float m4[4][4], + const float m5[4][4], const float m6[4][4]) { mul_m4_m4m4(r, m1, m2); mul_m4_m4m4(r, r, m3); @@ -409,8 +476,8 @@ void _va_mul_m4_series_7( } void _va_mul_m4_series_8( float r[4][4], - float m1[4][4], float m2[4][4], float m3[4][4], float m4[4][4], - float m5[4][4], float m6[4][4], float m7[4][4]) + const float m1[4][4], const float m2[4][4], const float m3[4][4], const float m4[4][4], + const float m5[4][4], const float m6[4][4], const float m7[4][4]) { mul_m4_m4m4(r, m1, m2); mul_m4_m4m4(r, r, m3); @@ -421,8 +488,8 @@ void _va_mul_m4_series_8( } void _va_mul_m4_series_9( float r[4][4], - float m1[4][4], float m2[4][4], float m3[4][4], float m4[4][4], - float m5[4][4], float m6[4][4], float m7[4][4], float m8[4][4]) + const float m1[4][4], const float m2[4][4], const float m3[4][4], const float m4[4][4], + const float m5[4][4], const float m6[4][4], const float m7[4][4], const float m8[4][4]) { mul_m4_m4m4(r, m1, m2); mul_m4_m4m4(r, r, m3); @@ -434,7 +501,7 @@ void _va_mul_m4_series_9( } /** \} */ -void mul_v2_m3v2(float r[2], float m[3][3], float v[2]) +void mul_v2_m3v2(float r[2], const float m[3][3], const float v[2]) { float temp[3], warped[3]; @@ -447,12 +514,12 @@ void mul_v2_m3v2(float r[2], float m[3][3], float v[2]) r[1] = warped[1] / warped[2]; } -void mul_m3_v2(float m[3][3], float r[2]) +void mul_m3_v2(const float m[3][3], float r[2]) { mul_v2_m3v2(r, m, r); } -void mul_m4_v3(float mat[4][4], float vec[3]) +void mul_m4_v3(const float mat[4][4], float vec[3]) { const float x = vec[0]; const float y = vec[1]; @@ -462,7 +529,7 @@ void mul_m4_v3(float mat[4][4], float vec[3]) vec[2] = x * mat[0][2] + y * mat[1][2] + mat[2][2] * vec[2] + mat[3][2]; } -void mul_v3_m4v3(float r[3], float mat[4][4], const float vec[3]) +void mul_v3_m4v3(float r[3], const float mat[4][4], const float vec[3]) { const float x = vec[0]; const float y = vec[1]; @@ -472,7 +539,7 @@ void mul_v3_m4v3(float r[3], float mat[4][4], const float vec[3]) r[2] = x * mat[0][2] + y * mat[1][2] + mat[2][2] * vec[2] + mat[3][2]; } -void mul_v2_m4v3(float r[2], float mat[4][4], const float vec[3]) +void mul_v2_m4v3(float r[2], const float mat[4][4], const float vec[3]) { const float x = vec[0]; @@ -480,7 +547,7 @@ void mul_v2_m4v3(float r[2], float mat[4][4], const float vec[3]) r[1] = x * mat[0][1] + vec[1] * mat[1][1] + mat[2][1] * vec[2] + mat[3][1]; } -void mul_v2_m2v2(float r[2], float mat[2][2], const float vec[2]) +void mul_v2_m2v2(float r[2], const float mat[2][2], const float vec[2]) { const float x = vec[0]; @@ -488,13 +555,13 @@ void mul_v2_m2v2(float r[2], float mat[2][2], const float vec[2]) r[1] = mat[0][1] * x + mat[1][1] * vec[1]; } -void mul_m2v2(float mat[2][2], float vec[2]) +void mul_m2v2(const float mat[2][2], float vec[2]) { mul_v2_m2v2(vec, mat, vec); } /* same as mul_m4_v3() but doesnt apply translation component */ -void mul_mat3_m4_v3(float mat[4][4], float vec[3]) +void mul_mat3_m4_v3(const float mat[4][4], float vec[3]) { const float x = vec[0]; const float y = vec[1]; @@ -504,7 +571,7 @@ void mul_mat3_m4_v3(float mat[4][4], float vec[3]) vec[2] = x * mat[0][2] + y * mat[1][2] + mat[2][2] * vec[2]; } -void mul_v3_mat3_m4v3(float r[3], float mat[4][4], const float vec[3]) +void mul_v3_mat3_m4v3(float r[3], const float mat[4][4], const float vec[3]) { const float x = vec[0]; const float y = vec[1]; @@ -514,7 +581,7 @@ void mul_v3_mat3_m4v3(float r[3], float mat[4][4], const float vec[3]) r[2] = x * mat[0][2] + y * mat[1][2] + mat[2][2] * vec[2]; } -void mul_project_m4_v3(float mat[4][4], float vec[3]) +void mul_project_m4_v3(const float mat[4][4], float vec[3]) { /* absolute value to not flip the frustum upside down behind the camera */ const float w = fabsf(mul_project_m4_v3_zfac(mat, vec)); @@ -525,7 +592,7 @@ void mul_project_m4_v3(float mat[4][4], float vec[3]) vec[2] /= w; } -void mul_v3_project_m4_v3(float r[3], float mat[4][4], const float vec[3]) +void mul_v3_project_m4_v3(float r[3], const float mat[4][4], const float vec[3]) { const float w = fabsf(mul_project_m4_v3_zfac(mat, vec)); mul_v3_m4v3(r, mat, vec); @@ -535,7 +602,7 @@ void mul_v3_project_m4_v3(float r[3], float mat[4][4], const float vec[3]) r[2] /= w; } -void mul_v2_project_m4_v3(float r[2], float mat[4][4], const float vec[3]) +void mul_v2_project_m4_v3(float r[2], const float mat[4][4], const float vec[3]) { const float w = fabsf(mul_project_m4_v3_zfac(mat, vec)); mul_v2_m4v3(r, mat, vec); @@ -544,7 +611,7 @@ void mul_v2_project_m4_v3(float r[2], float mat[4][4], const float vec[3]) r[1] /= w; } -void mul_v4_m4v4(float r[4], float mat[4][4], const float v[4]) +void mul_v4_m4v4(float r[4], const float mat[4][4], const float v[4]) { const float x = v[0]; const float y = v[1]; @@ -556,12 +623,12 @@ void mul_v4_m4v4(float r[4], float mat[4][4], const float v[4]) r[3] = x * mat[0][3] + y * mat[1][3] + z * mat[2][3] + mat[3][3] * v[3]; } -void mul_m4_v4(float mat[4][4], float r[4]) +void mul_m4_v4(const float mat[4][4], float r[4]) { mul_v4_m4v4(r, mat, r); } -void mul_v4d_m4v4d(double r[4], float mat[4][4], double v[4]) +void mul_v4d_m4v4d(double r[4], const float mat[4][4], const double v[4]) { const double x = v[0]; const double y = v[1]; @@ -573,12 +640,21 @@ void mul_v4d_m4v4d(double r[4], float mat[4][4], double v[4]) r[3] = x * (double)mat[0][3] + y * (double)mat[1][3] + z * (double)mat[2][3] + (double)mat[3][3] * v[3]; } -void mul_m4_v4d(float mat[4][4], double r[4]) +void mul_m4_v4d(const float mat[4][4], double r[4]) { mul_v4d_m4v4d(r, mat, r); } -void mul_v3_m3v3(float r[3], float M[3][3], const float a[3]) +void mul_v4_m4v3(float r[4], const float M[4][4], const float v[3]) +{ + /* v has implicit w = 1.0f */ + r[0] = v[0] * M[0][0] + v[1] * M[1][0] + M[2][0] * v[2] + M[3][0]; + r[1] = v[0] * M[0][1] + v[1] * M[1][1] + M[2][1] * v[2] + M[3][1]; + r[2] = v[0] * M[0][2] + v[1] * M[1][2] + M[2][2] * v[2] + M[3][2]; + r[3] = v[0] * M[0][3] + v[1] * M[1][3] + M[2][3] * v[2] + M[3][3]; +} + +void mul_v3_m3v3(float r[3], const float M[3][3], const float a[3]) { BLI_assert(r != a); @@ -587,7 +663,7 @@ void mul_v3_m3v3(float r[3], float M[3][3], const float a[3]) r[2] = M[0][2] * a[0] + M[1][2] * a[1] + M[2][2] * a[2]; } -void mul_v3_m3v3_db(double r[3], double M[3][3], const double a[3]) +void mul_v3_m3v3_db(double r[3], const double M[3][3], const double a[3]) { BLI_assert(r != a); @@ -596,7 +672,7 @@ void mul_v3_m3v3_db(double r[3], double M[3][3], const double a[3]) r[2] = M[0][2] * a[0] + M[1][2] * a[1] + M[2][2] * a[2]; } -void mul_v2_m3v3(float r[2], float M[3][3], const float a[3]) +void mul_v2_m3v3(float r[2], const float M[3][3], const float a[3]) { BLI_assert(r != a); @@ -604,17 +680,17 @@ void mul_v2_m3v3(float r[2], float M[3][3], const float a[3]) r[1] = M[0][1] * a[0] + M[1][1] * a[1] + M[2][1] * a[2]; } -void mul_m3_v3(float M[3][3], float r[3]) +void mul_m3_v3(const float M[3][3], float r[3]) { mul_v3_m3v3(r, M, (const float[3]){UNPACK3(r)}); } -void mul_m3_v3_db(double M[3][3], double r[3]) +void mul_m3_v3_db(const double M[3][3], double r[3]) { mul_v3_m3v3_db(r, M, (const double[3]){UNPACK3(r)}); } -void mul_transposed_m3_v3(float mat[3][3], float vec[3]) +void mul_transposed_m3_v3(const float mat[3][3], float vec[3]) { const float x = vec[0]; const float y = vec[1]; @@ -624,7 +700,7 @@ void mul_transposed_m3_v3(float mat[3][3], float vec[3]) vec[2] = x * mat[2][0] + y * mat[2][1] + mat[2][2] * vec[2]; } -void mul_transposed_mat3_m4_v3(float mat[4][4], float vec[3]) +void mul_transposed_mat3_m4_v3(const float mat[4][4], float vec[3]) { const float x = vec[0]; const float y = vec[1]; @@ -688,7 +764,7 @@ void negate_m4(float m[4][4]) m[i][j] *= -1.0f; } -void mul_m3_v3_double(float mat[3][3], double vec[3]) +void mul_m3_v3_double(const float mat[3][3], double vec[3]) { const double x = vec[0]; const double y = vec[1]; @@ -698,7 +774,7 @@ void mul_m3_v3_double(float mat[3][3], double vec[3]) vec[2] = x * (double)mat[0][2] + y * (double)mat[1][2] + (double)mat[2][2] * vec[2]; } -void add_m3_m3m3(float m1[3][3], float m2[3][3], float m3[3][3]) +void add_m3_m3m3(float m1[3][3], const float m2[3][3], const float m3[3][3]) { int i, j; @@ -707,7 +783,7 @@ void add_m3_m3m3(float m1[3][3], float m2[3][3], float m3[3][3]) m1[i][j] = m2[i][j] + m3[i][j]; } -void add_m4_m4m4(float m1[4][4], float m2[4][4], float m3[4][4]) +void add_m4_m4m4(float m1[4][4], const float m2[4][4], const float m3[4][4]) { int i, j; @@ -716,7 +792,7 @@ void add_m4_m4m4(float m1[4][4], float m2[4][4], float m3[4][4]) m1[i][j] = m2[i][j] + m3[i][j]; } -void sub_m3_m3m3(float m1[3][3], float m2[3][3], float m3[3][3]) +void sub_m3_m3m3(float m1[3][3], const float m2[3][3], const float m3[3][3]) { int i, j; @@ -725,7 +801,7 @@ void sub_m3_m3m3(float m1[3][3], float m2[3][3], float m3[3][3]) m1[i][j] = m2[i][j] - m3[i][j]; } -void sub_m4_m4m4(float m1[4][4], float m2[4][4], float m3[4][4]) +void sub_m4_m4m4(float m1[4][4], const float m2[4][4], const float m3[4][4]) { int i, j; @@ -734,7 +810,7 @@ void sub_m4_m4m4(float m1[4][4], float m2[4][4], float m3[4][4]) m1[i][j] = m2[i][j] - m3[i][j]; } -float determinant_m3_array(float m[3][3]) +float determinant_m3_array(const float m[3][3]) { return (m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]) - m[1][0] * (m[0][1] * m[2][2] - m[0][2] * m[2][1]) + @@ -750,7 +826,7 @@ bool invert_m3_ex(float m[3][3], const float epsilon) return success; } -bool invert_m3_m3_ex(float m1[3][3], float m2[3][3], const float epsilon) +bool invert_m3_m3_ex(float m1[3][3], const float m2[3][3], const float epsilon) { float det; int a, b; @@ -786,7 +862,7 @@ bool invert_m3(float m[3][3]) return success; } -bool invert_m3_m3(float m1[3][3], float m2[3][3]) +bool invert_m3_m3(float m1[3][3], const float m2[3][3]) { float det; int a, b; @@ -821,74 +897,11 @@ bool invert_m4(float m[4][4]) return success; } -/* - * invertmat - - * computes the inverse of mat and puts it in inverse. Returns - * true on success (i.e. can always find a pivot) and false on failure. - * Uses Gaussian Elimination with partial (maximal column) pivoting. - * - * Mark Segal - 1992 - */ - -bool invert_m4_m4(float inverse[4][4], float mat[4][4]) +bool invert_m4_m4(float inverse[4][4], const float mat[4][4]) { - int i, j, k; - double temp; - float tempmat[4][4]; - float max; - int maxj; - - BLI_assert(inverse != mat); - - /* 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 = fabsf(tempmat[i][i]); - maxj = i; - for (j = i + 1; j < 4; j++) { - if (fabsf(tempmat[j][i]) > max) { - max = fabsf(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]); - } - } - - if (UNLIKELY(tempmat[i][i] == 0.0f)) { - return false; /* No non-zero pivot */ - } - temp = (double)tempmat[i][i]; - for (k = 0; k < 4; k++) { - tempmat[i][k] = (float)((double)tempmat[i][k] / temp); - inverse[i][k] = (float)((double)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] -= (float)((double)tempmat[i][k] * temp); - inverse[j][k] -= (float)((double)inverse[i][k] * temp); - } - } - } - } - return true; + /* Use optimized matrix inverse from Eigen, since performance + * impact of this function is significant in complex rigs. */ + return EIG_invert_m4_m4(inverse, mat); } /****************************** Linear Algebra *******************************/ @@ -908,7 +921,7 @@ void transpose_m3(float mat[3][3]) mat[2][1] = t; } -void transpose_m3_m3(float rmat[3][3], float mat[3][3]) +void transpose_m3_m3(float rmat[3][3], const float mat[3][3]) { BLI_assert(rmat != mat); @@ -924,7 +937,7 @@ void transpose_m3_m3(float rmat[3][3], float mat[3][3]) } /* seems obscure but in-fact a common operation */ -void transpose_m3_m4(float rmat[3][3], float mat[4][4]) +void transpose_m3_m4(float rmat[3][3], const float mat[4][4]) { BLI_assert(&rmat[0][0] != &mat[0][0]); @@ -965,7 +978,7 @@ void transpose_m4(float mat[4][4]) mat[3][2] = t; } -void transpose_m4_m4(float rmat[4][4], float mat[4][4]) +void transpose_m4_m4(float rmat[4][4], const float mat[4][4]) { BLI_assert(rmat != mat); @@ -987,7 +1000,8 @@ void transpose_m4_m4(float rmat[4][4], float mat[4][4]) rmat[3][3] = mat[3][3]; } -int compare_m4m4(float mat1[4][4], float mat2[4][4], float limit) +/* TODO: return bool */ +int compare_m4m4(const float mat1[4][4], const float mat2[4][4], float limit) { if (compare_v4v4(mat1[0], mat2[0], limit)) if (compare_v4v4(mat1[1], mat2[1], limit)) @@ -1165,7 +1179,7 @@ void orthogonalize_m4(float mat[4][4], int axis) mul_v3_fl(mat[2], size[2]); } -bool is_orthogonal_m3(float m[3][3]) +bool is_orthogonal_m3(const float m[3][3]) { int i, j; @@ -1179,7 +1193,7 @@ bool is_orthogonal_m3(float m[3][3]) return true; } -bool is_orthogonal_m4(float m[4][4]) +bool is_orthogonal_m4(const float m[4][4]) { int i, j; @@ -1194,7 +1208,7 @@ bool is_orthogonal_m4(float m[4][4]) return true; } -bool is_orthonormal_m3(float m[3][3]) +bool is_orthonormal_m3(const float m[3][3]) { if (is_orthogonal_m3(m)) { int i; @@ -1209,7 +1223,7 @@ bool is_orthonormal_m3(float m[3][3]) return false; } -bool is_orthonormal_m4(float m[4][4]) +bool is_orthonormal_m4(const float m[4][4]) { if (is_orthogonal_m4(m)) { int i; @@ -1224,7 +1238,7 @@ bool is_orthonormal_m4(float m[4][4]) return false; } -bool is_uniform_scaled_m3(float m[3][3]) +bool is_uniform_scaled_m3(const float m[3][3]) { const float eps = 1e-7f; float t[3][3]; @@ -1252,7 +1266,7 @@ bool is_uniform_scaled_m3(float m[3][3]) return false; } -bool is_uniform_scaled_m4(float m[4][4]) +bool is_uniform_scaled_m4(const float m[4][4]) { float t[3][3]; copy_m3_m4(t, m); @@ -1274,14 +1288,14 @@ void normalize_m3(float mat[3][3]) } } -void normalize_m3_m3_ex(float rmat[3][3], float mat[3][3], float r_scale[3]) +void normalize_m3_m3_ex(float rmat[3][3], const float mat[3][3], float r_scale[3]) { int i; for (i = 0; i < 3; i++) { r_scale[i] = normalize_v3_v3(rmat[i], mat[i]); } } -void normalize_m3_m3(float rmat[3][3], float mat[3][3]) +void normalize_m3_m3(float rmat[3][3], const float mat[3][3]) { int i; for (i = 0; i < 3; i++) { @@ -1310,7 +1324,7 @@ void normalize_m4(float mat[4][4]) } } -void normalize_m4_m4_ex(float rmat[4][4], float mat[4][4], float r_scale[3]) +void normalize_m4_m4_ex(float rmat[4][4], const float mat[4][4], float r_scale[3]) { int i; for (i = 0; i < 3; i++) { @@ -1319,7 +1333,7 @@ void normalize_m4_m4_ex(float rmat[4][4], float mat[4][4], float r_scale[3]) } copy_v4_v4(rmat[3], mat[3]); } -void normalize_m4_m4(float rmat[4][4], float mat[4][4]) +void normalize_m4_m4(float rmat[4][4], const float mat[4][4]) { int i; for (i = 0; i < 3; i++) { @@ -1329,7 +1343,7 @@ void normalize_m4_m4(float rmat[4][4], float mat[4][4]) copy_v4_v4(rmat[3], mat[3]); } -void adjoint_m2_m2(float m1[2][2], float m[2][2]) +void adjoint_m2_m2(float m1[2][2], const float m[2][2]) { BLI_assert(m1 != m); m1[0][0] = m[1][1]; @@ -1338,7 +1352,7 @@ void adjoint_m2_m2(float m1[2][2], float m[2][2]) m1[1][1] = m[0][0]; } -void adjoint_m3_m3(float m1[3][3], float m[3][3]) +void adjoint_m3_m3(float m1[3][3], const float m[3][3]) { BLI_assert(m1 != m); m1[0][0] = m[1][1] * m[2][2] - m[1][2] * m[2][1]; @@ -1354,7 +1368,7 @@ void adjoint_m3_m3(float m1[3][3], float m[3][3]) m1[2][2] = m[0][0] * m[1][1] - m[0][1] * m[1][0]; } -void adjoint_m4_m4(float out[4][4], float in[4][4]) /* out = ADJ(in) */ +void adjoint_m4_m4(float out[4][4], const float in[4][4]) /* out = ADJ(in) */ { float a1, a2, a3, a4, b1, b2, b3, b4; float c1, c2, c3, c4, d1, d2, d3, d4; @@ -1420,7 +1434,7 @@ float determinant_m3(float a1, float a2, float a3, return ans; } -float determinant_m4(float m[4][4]) +float determinant_m4(const float m[4][4]) { float ans; float a1, a2, a3, a4, b1, b2, b3, b4, c1, c2, c3, c4, d1, d2, d3, d4; @@ -1488,14 +1502,14 @@ void size_to_mat4(float mat[4][4], const float size[3]) mat[3][3] = 1.0f; } -void mat3_to_size(float size[3], float mat[3][3]) +void mat3_to_size(float size[3], const float mat[3][3]) { size[0] = len_v3(mat[0]); size[1] = len_v3(mat[1]); size[2] = len_v3(mat[2]); } -void mat4_to_size(float size[3], float mat[4][4]) +void mat4_to_size(float size[3], const float mat[4][4]) { size[0] = len_v3(mat[0]); size[1] = len_v3(mat[1]); @@ -1505,7 +1519,7 @@ void mat4_to_size(float size[3], float mat[4][4]) /* this gets the average scale of a matrix, only use when your scaling * data that has no idea of scale axis, examples are bone-envelope-radius * and curve radius */ -float mat3_to_scale(float mat[3][3]) +float mat3_to_scale(const float mat[3][3]) { /* unit length vector */ float unit_vec[3]; @@ -1514,7 +1528,7 @@ float mat3_to_scale(float mat[3][3]) return len_v3(unit_vec); } -float mat4_to_scale(float mat[4][4]) +float mat4_to_scale(const float mat[4][4]) { /* unit length vector */ float unit_vec[3]; @@ -1524,7 +1538,7 @@ float mat4_to_scale(float mat[4][4]) } /** Return 2D scale (in XY plane) of given mat4. */ -float mat4_to_xy_scale(float M[4][4]) +float mat4_to_xy_scale(const float M[4][4]) { /* unit length vector in xy plane */ float unit_vec[3] = {(float)M_SQRT1_2, (float)M_SQRT1_2, 0.0f}; @@ -1532,7 +1546,7 @@ float mat4_to_xy_scale(float M[4][4]) return len_v3(unit_vec); } -void mat3_to_rot_size(float rot[3][3], float size[3], float mat3[3][3]) +void mat3_to_rot_size(float rot[3][3], float size[3], const float mat3[3][3]) { /* keep rot as a 3x3 matrix, the caller can convert into a quat or euler */ size[0] = normalize_v3_v3(rot[0], mat3[0]); @@ -1544,7 +1558,7 @@ void mat3_to_rot_size(float rot[3][3], float size[3], float mat3[3][3]) } } -void mat4_to_loc_rot_size(float loc[3], float rot[3][3], float size[3], float wmat[4][4]) +void mat4_to_loc_rot_size(float loc[3], float rot[3][3], float size[3], const float wmat[4][4]) { float mat3[3][3]; /* wmat -> 3x3 */ @@ -1555,7 +1569,7 @@ void mat4_to_loc_rot_size(float loc[3], float rot[3][3], float size[3], float wm copy_v3_v3(loc, wmat[3]); } -void mat4_to_loc_quat(float loc[3], float quat[4], float wmat[4][4]) +void mat4_to_loc_quat(float loc[3], float quat[4], const float wmat[4][4]) { float mat3[3][3]; float mat3_n[3][3]; /* normalized mat3 */ @@ -1573,7 +1587,7 @@ void mat4_to_loc_quat(float loc[3], float quat[4], float wmat[4][4]) copy_v3_v3(loc, wmat[3]); } -void mat4_decompose(float loc[3], float quat[4], float size[3], float wmat[4][4]) +void mat4_decompose(float loc[3], float quat[4], float size[3], const float wmat[4][4]) { float rot[3][3]; mat4_to_loc_rot_size(loc, rot, size, wmat); @@ -1590,7 +1604,7 @@ void mat4_decompose(float loc[3], float quat[4], float size[3], float wmat[4][4] * See https://en.wikipedia.org/wiki/Polar_decomposition for more. */ #ifndef MATH_STANDALONE -void mat3_polar_decompose(float mat3[3][3], float r_U[3][3], float r_P[3][3]) +void mat3_polar_decompose(const float mat3[3][3], float r_U[3][3], float r_P[3][3]) { /* From svd decomposition (M = WSV*), we have: * U = WV* @@ -1701,7 +1715,7 @@ void transform_pivot_set_m4(float mat[4][4], const float pivot[3]) mul_m4_m4m4(mat, mat, tmat); } -void blend_m3_m3m3(float out[3][3], float dst[3][3], float src[3][3], const float srcweight) +void blend_m3_m3m3(float out[3][3], const float dst[3][3], const float src[3][3], const float srcweight) { float srot[3][3], drot[3][3]; float squat[4], dquat[4], fquat[4]; @@ -1724,7 +1738,7 @@ void blend_m3_m3m3(float out[3][3], float dst[3][3], float src[3][3], const floa mul_m3_m3m3(out, rmat, smat); } -void blend_m4_m4m4(float out[4][4], float dst[4][4], float src[4][4], const float srcweight) +void blend_m4_m4m4(float out[4][4], const float dst[4][4], const float src[4][4], const float srcweight) { float sloc[3], dloc[3], floc[3]; float srot[3][3], drot[3][3]; @@ -1762,7 +1776,7 @@ void blend_m4_m4m4(float out[4][4], float dst[4][4], float src[4][4], const floa * \param B: Input matrix which is totally effective with `t = 1.0`. * \param t: Interpolation factor. */ -void interp_m3_m3m3(float R[3][3], float A[3][3], float B[3][3], const float t) +void interp_m3_m3m3(float R[3][3], const float A[3][3], const float B[3][3], const float t) { /* 'Rotation' component ('U' part of polar decomposition, the closest orthogonal matrix to M3 rot/scale * transformation matrix), spherically interpolated. */ @@ -1797,7 +1811,7 @@ void interp_m3_m3m3(float R[3][3], float A[3][3], float B[3][3], const float t) * \param B: Input matrix which is totally effective with `t = 1.0`. * \param t: Interpolation factor. */ -void interp_m4_m4m4(float R[4][4], float A[4][4], float B[4][4], const float t) +void interp_m4_m4m4(float R[4][4], const float A[4][4], const float B[4][4], const float t) { float A3[3][3], B3[3][3], R3[3][3]; @@ -1818,27 +1832,27 @@ void interp_m4_m4m4(float R[4][4], float A[4][4], float B[4][4], const float t) } #endif /* MATH_STANDALONE */ -bool is_negative_m3(float mat[3][3]) +bool is_negative_m3(const float mat[3][3]) { float vec[3]; cross_v3_v3v3(vec, mat[0], mat[1]); return (dot_v3v3(vec, mat[2]) < 0.0f); } -bool is_negative_m4(float mat[4][4]) +bool is_negative_m4(const float mat[4][4]) { float vec[3]; cross_v3_v3v3(vec, mat[0], mat[1]); return (dot_v3v3(vec, mat[2]) < 0.0f); } -bool is_zero_m3(float mat[3][3]) +bool is_zero_m3(const float mat[3][3]) { return (is_zero_v3(mat[0]) && is_zero_v3(mat[1]) && is_zero_v3(mat[2])); } -bool is_zero_m4(float mat[4][4]) +bool is_zero_m4(const float mat[4][4]) { return (is_zero_v4(mat[0]) && is_zero_v4(mat[1]) && @@ -1846,14 +1860,14 @@ bool is_zero_m4(float mat[4][4]) is_zero_v4(mat[3])); } -bool equals_m3m3(float mat1[3][3], float mat2[3][3]) +bool equals_m3m3(const float mat1[3][3], const float mat2[3][3]) { return (equals_v3v3(mat1[0], mat2[0]) && equals_v3v3(mat1[1], mat2[1]) && equals_v3v3(mat1[2], mat2[2])); } -bool equals_m4m4(float mat1[4][4], float mat2[4][4]) +bool equals_m4m4(const float mat1[4][4], const float mat2[4][4]) { return (equals_v4v4(mat1[0], mat2[0]) && equals_v4v4(mat1[1], mat2[1]) && @@ -1944,7 +1958,7 @@ void loc_axisangle_size_to_mat4(float mat[4][4], const float loc[3], const float /*********************************** Other ***********************************/ -void print_m3(const char *str, float m[3][3]) +void print_m3(const char *str, const float m[3][3]) { printf("%s\n", str); printf("%f %f %f\n", m[0][0], m[1][0], m[2][0]); @@ -1953,7 +1967,7 @@ void print_m3(const char *str, float m[3][3]) printf("\n"); } -void print_m4(const char *str, float m[4][4]) +void print_m4(const char *str, const float m[4][4]) { printf("%s\n", str); printf("%f %f %f %f\n", m[0][0], m[1][0], m[2][0], m[3][0]); @@ -2407,7 +2421,7 @@ void svd_m4(float U[4][4], float s[4], float V[4][4], float A_[4][4]) } } -void pseudoinverse_m4_m4(float Ainv[4][4], float A_[4][4], float epsilon) +void pseudoinverse_m4_m4(float Ainv[4][4], const float A_[4][4], float epsilon) { /* compute Moore-Penrose pseudo inverse of matrix, singular values * below epsilon are ignored for stability (truncated SVD) */ @@ -2428,7 +2442,7 @@ void pseudoinverse_m4_m4(float Ainv[4][4], float A_[4][4], float epsilon) mul_m4_series(Ainv, U, Wm, V); } -void pseudoinverse_m3_m3(float Ainv[3][3], float A[3][3], float epsilon) +void pseudoinverse_m3_m3(float Ainv[3][3], const float A[3][3], float epsilon) { /* try regular inverse when possible, otherwise fall back to slow svd */ if (!invert_m3_m3(Ainv, A)) { @@ -2440,14 +2454,14 @@ void pseudoinverse_m3_m3(float Ainv[3][3], float A[3][3], float epsilon) } } -bool has_zero_axis_m4(float matrix[4][4]) +bool has_zero_axis_m4(const float matrix[4][4]) { return len_squared_v3(matrix[0]) < FLT_EPSILON || len_squared_v3(matrix[1]) < FLT_EPSILON || len_squared_v3(matrix[2]) < FLT_EPSILON; } -void invert_m4_m4_safe(float Ainv[4][4], float A[4][4]) +void invert_m4_m4_safe(float Ainv[4][4], const float A[4][4]) { if (!invert_m4_m4(Ainv, A)) { float Atemp[4][4]; @@ -2495,7 +2509,7 @@ void invert_m4_m4_safe(float Ainv[4][4], float A[4][4]) * this defines a transform matrix TM such that (x', y', z') = TM * (x, y, z) * where (x', y', z') are the coordinates of P' in target space such that it keeps (X, Y, Z) coordinates in global space. */ -void BLI_space_transform_from_matrices(SpaceTransform *data, float local[4][4], float target[4][4]) +void BLI_space_transform_from_matrices(SpaceTransform *data, const float local[4][4], const float target[4][4]) { float itarget[4][4]; invert_m4_m4(itarget, target); @@ -2513,7 +2527,7 @@ void BLI_space_transform_from_matrices(SpaceTransform *data, float local[4][4], * this defines a transform matrix TM such that (X', Y', Z') = TM * (X, Y, Z) * where (X', Y', Z') are the coordinates of p' in global space such that it keeps (x, y, z) coordinates in target space. */ -void BLI_space_transform_global_from_matrices(SpaceTransform *data, float local[4][4], float target[4][4]) +void BLI_space_transform_global_from_matrices(SpaceTransform *data, const float local[4][4], const float target[4][4]) { float ilocal[4][4]; invert_m4_m4(ilocal, local); diff --git a/source/blender/blenlib/intern/math_vector_inline.c b/source/blender/blenlib/intern/math_vector_inline.c index 715e2e65c96..189b94a6f13 100644 --- a/source/blender/blenlib/intern/math_vector_inline.c +++ b/source/blender/blenlib/intern/math_vector_inline.c @@ -491,7 +491,7 @@ MINLINE void mul_v2_v2_ccw(float r[2], const float mat[2], const float vec[2]) * return co_4d[3]; * \endcode */ -MINLINE float mul_project_m4_v3_zfac(float mat[4][4], const float co[3]) +MINLINE float mul_project_m4_v3_zfac(const float mat[4][4], const float co[3]) { return (mat[0][3] * co[0]) + (mat[1][3] * co[1]) + @@ -501,15 +501,15 @@ MINLINE float mul_project_m4_v3_zfac(float mat[4][4], const float co[3]) /** * Has the effect of #mul_m3_v3(), on a single axis. */ -MINLINE float dot_m3_v3_row_x(float M[3][3], const float a[3]) +MINLINE float dot_m3_v3_row_x(const float M[3][3], const float a[3]) { return M[0][0] * a[0] + M[1][0] * a[1] + M[2][0] * a[2]; } -MINLINE float dot_m3_v3_row_y(float M[3][3], const float a[3]) +MINLINE float dot_m3_v3_row_y(const float M[3][3], const float a[3]) { return M[0][1] * a[0] + M[1][1] * a[1] + M[2][1] * a[2]; } -MINLINE float dot_m3_v3_row_z(float M[3][3], const float a[3]) +MINLINE float dot_m3_v3_row_z(const float M[3][3], const float a[3]) { return M[0][2] * a[0] + M[1][2] * a[1] + M[2][2] * a[2]; } @@ -518,15 +518,15 @@ MINLINE float dot_m3_v3_row_z(float M[3][3], const float a[3]) * Has the effect of #mul_mat3_m4_v3(), on a single axis. * (no adding translation) */ -MINLINE float dot_m4_v3_row_x(float M[4][4], const float a[3]) +MINLINE float dot_m4_v3_row_x(const float M[4][4], const float a[3]) { return M[0][0] * a[0] + M[1][0] * a[1] + M[2][0] * a[2]; } -MINLINE float dot_m4_v3_row_y(float M[4][4], const float a[3]) +MINLINE float dot_m4_v3_row_y(const float M[4][4], const float a[3]) { return M[0][1] * a[0] + M[1][1] * a[1] + M[2][1] * a[2]; } -MINLINE float dot_m4_v3_row_z(float M[4][4], const float a[3]) +MINLINE float dot_m4_v3_row_z(const float M[4][4], const float a[3]) { return M[0][2] * a[0] + M[1][2] * a[1] + M[2][2] * a[2]; } diff --git a/source/blender/blenlib/intern/rand.c b/source/blender/blenlib/intern/rand.c index 110757ac3c0..9e56ce6b2cf 100644 --- a/source/blender/blenlib/intern/rand.c +++ b/source/blender/blenlib/intern/rand.c @@ -265,29 +265,6 @@ void BLI_rng_skip(RNG *rng, int n) /***/ -/* initialize with some non-zero seed */ -static RNG theBLI_rng = {611330372042337130}; - -void BLI_srandom(unsigned int seed) -{ - BLI_rng_srandom(&theBLI_rng, seed); -} - -int BLI_rand(void) -{ - return BLI_rng_get_int(&theBLI_rng); -} - -float BLI_frand(void) -{ - return BLI_rng_get_float(&theBLI_rng); -} - -void BLI_frand_unit_v3(float v[3]) -{ - BLI_rng_get_float_unit_v3(&theBLI_rng, v); -} - float BLI_hash_frand(unsigned int seed) { RNG rng; @@ -386,6 +363,8 @@ void BLI_halton_1D(unsigned int prime, double offset, int n, double *r) { const double invprime = 1.0 / (double)prime; + *r = 0.0; + for (int s = 0; s < n; s++) { *r = halton_ex(invprime, &offset); } @@ -395,6 +374,8 @@ void BLI_halton_2D(unsigned int prime[2], double offset[2], int n, double *r) { const double invprimes[2] = {1.0 / (double)prime[0], 1.0 / (double)prime[1]}; + r[0] = r[1] = 0.0; + for (int s = 0; s < n; s++) { for (int i = 0; i < 2; i++) { r[i] = halton_ex(invprimes[i], &offset[i]); @@ -406,6 +387,8 @@ void BLI_halton_3D(unsigned int prime[3], double offset[3], int n, double *r) { const double invprimes[3] = {1.0 / (double)prime[0], 1.0 / (double)prime[1], 1.0 / (double)prime[2]}; + r[0] = r[1] = r[2] = 0.0; + for (int s = 0; s < n; s++) { for (int i = 0; i < 3; i++) { r[i] = halton_ex(invprimes[i], &offset[i]); diff --git a/source/blender/blenlib/intern/scanfill_utils.c b/source/blender/blenlib/intern/scanfill_utils.c index 3b6ab99ae86..ef1df479f13 100644 --- a/source/blender/blenlib/intern/scanfill_utils.c +++ b/source/blender/blenlib/intern/scanfill_utils.c @@ -95,20 +95,6 @@ void BLI_scanfill_obj_dump(ScanFillContext *sf_ctx) } #endif -#if 0 -void BLI_scanfill_view3d_dump(ScanFillContext *sf_ctx) -{ - ScanFillEdge *eed; - - bl_debug_draw_quad_clear(); - bl_debug_color_set(0x0000ff); - - for (eed = sf_ctx->filledgebase.first; eed; eed = eed->next) { - bl_debug_draw_edge_add(eed->v1->co, eed->v2->co); - } -} -#endif - static ListBase *edge_isect_ls_ensure(GHash *isect_hash, ScanFillEdge *eed) { ListBase *e_ls; |