add math

client/CMakeLists.txt

@@ -232,7 +232,7 @@ endif()
 
 # Subdirectories
 if (WITH_CLIENT_CUDA)
-	#add_subdirectory(cycles_cuda)
+	add_subdirectory(cycles_cuda)
 endif()
 
 if (WITH_CLIENT_RENDERENGINE)
@@ -245,6 +245,7 @@ if (WITH_CLIENT_CESNET)
     add_subdirectory(ultragrid)
 endif()
 
+add_subdirectory(blenlib)
 add_subdirectory(cycles)
 add_subdirectory(blender_client)
 

client/blenlib/CMakeLists.txt

+set(INC
+	.
+	../../extern/Eigen3
+    ../../source/blender/blenlib
+	../../source/blender/makesdna
+)
+
+set(SRC
+         math_matrix.cpp
+         math_rotation.cpp
+         math_vector.cpp
+         math_vector_inline.cpp
+		 math_base.cpp
+		 math_base_inline.cpp
+)
+
+set(SRC_HEADERS
+)
+
+if(WITH_OPENMP)
+	add_definitions(-DWITH_OPENMP)
+endif()
+
+add_definitions(-D__SSE2__)
+
+include_directories(${INC})
+
+add_library(blenlib STATIC ${SRC} ${SRC_HEADERS})
+
+install (TARGETS blenlib DESTINATION lib)

client/blenlib/math_base.cpp

+/*
+ * 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) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: some of this file.
+ *
+ * */
+
+/** \file
+ * \ingroup bli
+ */
+
+#include "BLI_math.h"
+
+#include "BLI_strict_flags.h"
+
+int pow_i(int base, int exp)
+{
+  int result = 1;
+  BLI_assert(exp >= 0);
+  while (exp) {
+    if (exp & 1) {
+      result *= base;
+    }
+    exp >>= 1;
+    base *= base;
+  }
+
+  return result;
+}
+
+/* from python 3.1 floatobject.c
+ * ndigits must be between 0 and 21 */
+double double_round(double x, int ndigits)
+{
+  double pow1, pow2, y, z;
+  if (ndigits >= 0) {
+    pow1 = pow(10.0, (double)ndigits);
+    pow2 = 1.0;
+    y = (x * pow1) * pow2;
+    /* if y overflows, then rounded value is exactly x */
+    if (!isfinite(y)) {
+      return x;
+    }
+  }
+  else {
+    pow1 = pow(10.0, (double)-ndigits);
+    pow2 = 1.0; /* unused; silences a gcc compiler warning */
+    y = x / pow1;
+  }
+
+  z = round(y);
+  if (fabs(y - z) == 0.5) {
+    /* halfway between two integers; use round-half-even */
+    z = 2.0 * round(y / 2.0);
+  }
+
+  if (ndigits >= 0) {
+    z = (z / pow2) / pow1;
+  }
+  else {
+    z *= pow1;
+  }
+
+  /* if computation resulted in overflow, raise OverflowError */
+  return z;
+}

client/blenlib/math_base_inline.cpp

+/*
+ * 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) 2001-2002 by NaN Holding BV.
+ * All rights reserved.
+ *
+ * The Original Code is: some of this file.
+ *
+ * */
+
+/** \file
+ * \ingroup bli
+ */
+
+#ifndef __MATH_BASE_INLINE_C__
+#define __MATH_BASE_INLINE_C__
+
+#include <float.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <limits.h>
+
+#ifdef __SSE2__
+#  include <emmintrin.h>
+#endif
+
+#include "BLI_math_base.h"
+
+/* copied from BLI_utildefines.h */
+#ifdef __GNUC__
+#  define UNLIKELY(x) __builtin_expect(!!(x), 0)
+#else
+#  define UNLIKELY(x) (x)
+#endif
+
+//#define BLI_assert
+
+/* powf is really slow for raising to integer powers. */
+MINLINE float pow2f(float x)
+{
+  return x * x;
+}
+MINLINE float pow3f(float x)
+{
+  return pow2f(x) * x;
+}
+MINLINE float pow4f(float x)
+{
+  return pow2f(pow2f(x));
+}
+MINLINE float pow7f(float x)
+{
+  return pow2f(pow3f(x)) * x;
+}
+
+MINLINE float sqrt3f(float f)
+{
+  if (UNLIKELY(f == 0.0f)) {
+    return 0.0f;
+  }
+  else if (UNLIKELY(f < 0.0f)) {
+    return -(float)(exp(log(-f) / 3.0));
+  }
+  else {
+    return (float)(exp(log(f) / 3.0));
+  }
+}
+
+MINLINE double sqrt3d(double d)
+{
+  if (UNLIKELY(d == 0.0)) {
+    return 0.0;
+  }
+  else if (UNLIKELY(d < 0.0)) {
+    return -exp(log(-d) / 3.0);
+  }
+  else {
+    return exp(log(d) / 3.0);
+  }
+}
+
+MINLINE float sqrtf_signed(float f)
+{
+  return (f >= 0.0f) ? sqrtf(f) : -sqrtf(-f);
+}
+
+MINLINE float saacos(float fac)
+{
+  if (UNLIKELY(fac <= -1.0f)) {
+    return (float)M_PI;
+  }
+  else if (UNLIKELY(fac >= 1.0f)) {
+    return 0.0f;
+  }
+  else {
+    return acosf(fac);
+  }
+}
+
+MINLINE float saasin(float fac)
+{
+  if (UNLIKELY(fac <= -1.0f)) {
+    return (float)-M_PI / 2.0f;
+  }
+  else if (UNLIKELY(fac >= 1.0f)) {
+    return (float)M_PI / 2.0f;
+  }
+  else {
+    return asinf(fac);
+  }
+}
+
+MINLINE float sasqrt(float fac)
+{
+  if (UNLIKELY(fac <= 0.0f)) {
+    return 0.0f;
+  }
+  else {
+    return sqrtf(fac);
+  }
+}
+
+MINLINE float saacosf(float fac)
+{
+  if (UNLIKELY(fac <= -1.0f)) {
+    return (float)M_PI;
+  }
+  else if (UNLIKELY(fac >= 1.0f)) {
+    return 0.0f;
+  }
+  else {
+    return acosf(fac);
+  }
+}
+
+MINLINE float saasinf(float fac)
+{
+  if (UNLIKELY(fac <= -1.0f)) {
+    return (float)-M_PI / 2.0f;
+  }
+  else if (UNLIKELY(fac >= 1.0f)) {
+    return (float)M_PI / 2.0f;
+  }
+  else {
+    return asinf(fac);
+  }
+}
+
+MINLINE float sasqrtf(float fac)
+{
+  if (UNLIKELY(fac <= 0.0f)) {
+    return 0.0f;
+  }
+  else {
+    return sqrtf(fac);
+  }
+}
+
+MINLINE float interpf(float target, float origin, float fac)
+{
+  return (fac * target) + (1.0f - fac) * origin;
+}
+
+/* used for zoom values*/
+MINLINE float power_of_2(float val)
+{
+  return (float)pow(2.0, ceil(log((double)val) / M_LN2));
+}
+
+MINLINE int is_power_of_2_i(int n)
+{
+  return (n & (n - 1)) == 0;
+}
+
+MINLINE int power_of_2_max_i(int n)
+{
+  if (is_power_of_2_i(n)) {
+    return n;
+  }
+
+  do {
+    n = n & (n - 1);
+  } while (!is_power_of_2_i(n));
+
+  return n * 2;
+}
+
+MINLINE int power_of_2_min_i(int n)
+{
+  while (!is_power_of_2_i(n)) {
+    n = n & (n - 1);
+  }
+
+  return n;
+}
+
+MINLINE unsigned int power_of_2_max_u(unsigned int x)
+{
+  x -= 1;
+  x |= (x >> 1);
+  x |= (x >> 2);
+  x |= (x >> 4);
+  x |= (x >> 8);
+  x |= (x >> 16);
+  return x + 1;
+}
+
+MINLINE unsigned power_of_2_min_u(unsigned x)
+{
+  x |= (x >> 1);
+  x |= (x >> 2);
+  x |= (x >> 4);
+  x |= (x >> 8);
+  x |= (x >> 16);
+  return x - (x >> 1);
+}
+
+/* rounding and clamping */
+
+#define _round_clamp_fl_impl(arg, ty, min, max) \
+  { \
+    float r = floorf(arg + 0.5f); \
+    if (UNLIKELY(r <= (float)min)) { \
+      return (ty)min; \
+    } \
+    else if (UNLIKELY(r >= (float)max)) { \
+      return (ty)max; \
+    } \
+    else { \
+      return (ty)r; \
+    } \
+  }
+
+#define _round_clamp_db_impl(arg, ty, min, max) \
+  { \
+    double r = floor(arg + 0.5); \
+    if (UNLIKELY(r <= (double)min)) { \
+      return (ty)min; \
+    } \
+    else if (UNLIKELY(r >= (double)max)) { \
+      return (ty)max; \
+    } \
+    else { \
+      return (ty)r; \
+    } \
+  }
+
+#define _round_fl_impl(arg, ty) \
+  { \
+    return (ty)floorf(arg + 0.5f); \
+  }
+#define _round_db_impl(arg, ty) \
+  { \
+    return (ty)floor(arg + 0.5); \
+  }
+
+MINLINE signed char round_fl_to_char(float a){_round_fl_impl(a, signed char)} MINLINE
+    unsigned char round_fl_to_uchar(float a){_round_fl_impl(a, unsigned char)} MINLINE
+    short round_fl_to_short(float a){_round_fl_impl(a, short)} MINLINE
+    unsigned short round_fl_to_ushort(float a){_round_fl_impl(a, unsigned short)} MINLINE
+    int round_fl_to_int(float a){_round_fl_impl(a, int)} MINLINE
+    unsigned int round_fl_to_uint(float a){_round_fl_impl(a, unsigned int)}
+
+MINLINE signed char round_db_to_char(double a){_round_db_impl(a, signed char)} MINLINE
+    unsigned char round_db_to_uchar(double a){_round_db_impl(a, unsigned char)} MINLINE
+    short round_db_to_short(double a){_round_db_impl(a, short)} MINLINE
+    unsigned short round_db_to_ushort(double a){_round_db_impl(a, unsigned short)} MINLINE
+    int round_db_to_int(double a){_round_db_impl(a, int)} MINLINE
+    unsigned int round_db_to_uint(double a)
+{
+  _round_db_impl(a, unsigned int)
+}
+
+#undef _round_fl_impl
+#undef _round_db_impl
+
+MINLINE signed char round_fl_to_char_clamp(float a){
+    _round_clamp_fl_impl(a, signed char, SCHAR_MIN, SCHAR_MAX)} MINLINE
+    unsigned char round_fl_to_uchar_clamp(float a){
+        _round_clamp_fl_impl(a, unsigned char, 0, UCHAR_MAX)} MINLINE
+    short round_fl_to_short_clamp(float a){
+        _round_clamp_fl_impl(a, short, SHRT_MIN, SHRT_MAX)} MINLINE
+    unsigned short round_fl_to_ushort_clamp(float a){
+        _round_clamp_fl_impl(a, unsigned short, 0, USHRT_MAX)} MINLINE
+    int round_fl_to_int_clamp(float a){_round_clamp_fl_impl(a, int, INT_MIN, INT_MAX)} MINLINE
+    unsigned int round_fl_to_uint_clamp(float a){
+        _round_clamp_fl_impl(a, unsigned int, 0, UINT_MAX)}
+
+MINLINE signed char round_db_to_char_clamp(double a){
+    _round_clamp_db_impl(a, signed char, SCHAR_MIN, SCHAR_MAX)} MINLINE
+    unsigned char round_db_to_uchar_clamp(double a){
+        _round_clamp_db_impl(a, unsigned char, 0, UCHAR_MAX)} MINLINE
+    short round_db_to_short_clamp(double a){
+        _round_clamp_db_impl(a, short, SHRT_MIN, SHRT_MAX)} MINLINE
+    unsigned short round_db_to_ushort_clamp(double a){
+        _round_clamp_db_impl(a, unsigned short, 0, USHRT_MAX)} MINLINE
+    int round_db_to_int_clamp(double a){_round_clamp_db_impl(a, int, INT_MIN, INT_MAX)} MINLINE
+    unsigned int round_db_to_uint_clamp(double a)
+{
+  _round_clamp_db_impl(a, unsigned int, 0, UINT_MAX)
+}
+
+#undef _round_clamp_fl_impl
+#undef _round_clamp_db_impl
+
+/* integer division that rounds 0.5 up, particularly useful for color blending
+ * with integers, to avoid gradual darkening when rounding down */
+MINLINE int divide_round_i(int a, int b)
+{
+  return (2 * a + b) / (2 * b);
+}
+
+/**
+ * Integer division that floors negative result.
+ * \note This works like Python's int division.
+ */
+MINLINE int divide_floor_i(int a, int b)
+{
+  int d = a / b;
+  int r = a % b; /* Optimizes into a single division. */
+  return r ? d - ((a < 0) ^ (b < 0)) : d;
+}
+
+/**
+ * modulo that handles negative numbers, works the same as Python's.
+ */
+MINLINE int mod_i(int i, int n)
+{
+  return (i % n + n) % n;
+}
+
+MINLINE float min_ff(float a, float b)
+{
+  return (a < b) ? a : b;
+}
+MINLINE float max_ff(float a, float b)
+{
+  return (a > b) ? a : b;
+}
+
+MINLINE int min_ii(int a, int b)
+{
+  return (a < b) ? a : b;
+}
+MINLINE int max_ii(int a, int b)
+{
+  return (b < a) ? a : b;
+}
+
+MINLINE float min_fff(float a, float b, float c)
+{
+  return min_ff(min_ff(a, b), c);
+}
+MINLINE float max_fff(float a, float b, float c)
+{
+  return max_ff(max_ff(a, b), c);
+}
+
+MINLINE int min_iii(int a, int b, int c)
+{
+  return min_ii(min_ii(a, b), c);
+}
+MINLINE int max_iii(int a, int b, int c)
+{
+  return max_ii(max_ii(a, b), c);
+}
+
+MINLINE float min_ffff(float a, float b, float c, float d)
+{
+  return min_ff(min_fff(a, b, c), d);
+}
+MINLINE float max_ffff(float a, float b, float c, float d)
+{
+  return max_ff(max_fff(a, b, c), d);
+}
+
+MINLINE int min_iiii(int a, int b, int c, int d)
+{
+  return min_ii(min_iii(a, b, c), d);
+}
+MINLINE int max_iiii(int a, int b, int c, int d)
+{
+  return max_ii(max_iii(a, b, c), d);
+}
+
+MINLINE size_t min_zz(size_t a, size_t b)
+{
+  return (a < b) ? a : b;
+}
+MINLINE size_t max_zz(size_t a, size_t b)
+{
+  return (b < a) ? a : b;
+}
+
+MINLINE int clamp_i(int value, int min, int max)
+{
+  return min_ii(max_ii(value, min), max);
+}
+
+MINLINE float clamp_f(float value, float min, float max)
+{
+  if (value > max) {
+    return max;
+  }
+  else if (value < min) {
+    return min;
+  }
+  return value;
+}
+
+MINLINE size_t clamp_z(size_t value, size_t min, size_t max)
+{
+  return min_zz(max_zz(value, min), max);
+}
+
+/**
+ * Almost-equal for IEEE floats, using absolute difference method.
+ *
+ * \param max_diff: the maximum absolute difference.
+ */
+MINLINE int compare_ff(float a, float b, const float max_diff)
+{
+  return fabsf(a - b) <= max_diff;
+}
+
+/**
+ * Almost-equal for IEEE floats, using their integer representation
+ * (mixing ULP and absolute difference methods).
+ *
+ * \param max_diff: is the maximum absolute difference (allows to take care of the near-zero area,
+ * where relative difference methods cannot really work).
+ * \param max_ulps: is the 'maximum number of floats + 1'
+ * allowed between \a a and \a b to consider them equal.
+ *
+ * \see https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
+ */
+MINLINE int compare_ff_relative(float a, float b, const float max_diff, const int max_ulps)
+{
+  union {
+    float f;
+    int i;
+  } ua, ub;
+
+  BLI_assert(sizeof(float) == sizeof(int));
+  BLI_assert(max_ulps < (1 << 22));
+
+  if (fabsf(a - b) <= max_diff) {
+    return 1;
+  }
+
+  ua.f = a;
+  ub.f = b;
+
+  /* Important to compare sign from integers, since (-0.0f < 0) is false
+   * (though this shall not be an issue in common cases)... */
+  return ((ua.i < 0) != (ub.i < 0)) ? 0 : (abs(ua.i - ub.i) <= max_ulps) ? 1 : 0;
+}
+
+MINLINE float signf(float f)
+{
+  return (f < 0.f) ? -1.f : 1.f;
+}
+
+MINLINE int signum_i_ex(float a, float eps)
+{
+  if (a > eps) {
+    return 1;
+  }
+  if (a < -eps) {
+    return -1;