Moved the square-root-function to `src/utils_math.c'.

[libopano.git] / src / utils_math.c

diff --git a/src/utils_math.c b/src/utils_math.c
new file mode 100644 (file)
index 0000000..09b424f
--- /dev/null
+++ b/src/utils_math.c
@@ -0,0 +1,141 @@
+/*
+ *  Copyright (C) 1998,1999  Helmut Dersch <der@fh-furtwangen.de>
+ *  Copyright (C) 2007  Florian octo Forster <octo at verplant.org>
+ *
+ *  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
+ **/
+
+#include <stdlib.h>
+#include <stdio.h>
+#include <math.h>
+
+#include "utils_math.h"
+
+// Lookup Tables for Trig-functions and interpolator
+
+#define ATAN_TABLE_SIZE 4096
+#define SQRT_TABLE_SIZE 4096
+
+static int *atan_table = NULL;
+static int *sqrt_table = NULL;
+
+static int setup_sqrt_table (void)
+{
+  double z;
+  int i;
+
+  sqrt_table = (int *) malloc (SQRT_TABLE_SIZE * sizeof (int));
+  if (sqrt_table == NULL)
+    return (-1);
+
+  /* Calculate the table entries. */
+  for (i = 0; i < SQRT_TABLE_SIZE; i++)
+  {
+    z = ((double) i) / ((double) (SQRT_TABLE_SIZE - 1));
+    sqrt_table[i] = (int) (sqrt( 1.0 + z*z ) * 4096.0 + .5);
+  }
+
+  return (0);
+} /* int setup_sqrt_table */
+
+/*
+ * Calculate
+ *   256 * sqrt(x1^2 + x2^2)
+ * efficiently
+ */
+int utils_sqrt2 (int x1, int x2)
+{
+  int ret;
+
+  if (sqrt_table == NULL)
+    if (setup_sqrt_table () != 0)
+      return (-1);
+
+  if ((x1 == 0) && (x2 == 0))
+  {
+    ret = 0;
+  }
+  else if (x1 > x2)
+  {
+    ret = x1 * sqrt_table[(SQRT_TABLE_SIZE - 1) * x2 / x1] / 16;
+  }
+  else /* (x2 < x1) */
+  {
+    ret = x2 * sqrt_table[(SQRT_TABLE_SIZE - 1) * x1 / x2] / 16;
+  }
+
+  return (ret);
+} /* int utils_sqrt2 */
+
+static int setup_atan_table (void)
+{
+  int i;
+  double z;
+
+  atan_table = (int *) malloc (ATAN_TABLE_SIZE * sizeof (int));
+  if (atan_table == NULL)
+    return (-1);
+
+  for (i = 0; i < (ATAN_TABLE_SIZE - 1); i++)
+  {
+    z = ((double) i) / ((double) (ATAN_TABLE_SIZE - 1));
+    atan_table[i] = atan ( z / (1.0 - z ) ) * ATAN_TABLE_SIZE * 256;
+  }
+  atan_table[ATAN_TABLE_SIZE - 1] = M_PI_4 * ATAN_TABLE_SIZE * 256;
+
+  return 0;
+} /* int setup_atan_table */
+
+int utils_atan2 (int y, int x)
+{
+  if (atan_table == NULL)
+    if (setup_atan_table () != 0)
+      return (-1);
+
+  // return atan2(y,x) * 256*ATAN_TABLE_SIZE;
+  if( x > 0 )
+  {
+    if( y > 0 )
+    {
+      return  atan_table[(int)( ATAN_TABLE_SIZE * y / ( x + y ))];
+    }
+    else
+    {
+      return -atan_table[ (int)(ATAN_TABLE_SIZE * (-y) / ( x - y ))];
+    }
+  }
+
+  if( x == 0 )
+  {
+    if( y > 0 )
+      return  (int)(256*ATAN_TABLE_SIZE * M_PI_2);
+    else
+      return  -(int)(256*ATAN_TABLE_SIZE * M_PI_2);
+  }
+
+  if( y < 0 )
+  {
+    return  atan_table[(int)( ATAN_TABLE_SIZE * y / ( x + y ))] - (int)(M_PI*256*ATAN_TABLE_SIZE);
+  }
+  else
+  {
+    return -atan_table[ (int)(ATAN_TABLE_SIZE * (-y) / ( x - y ))] + (int)(M_PI*256*ATAN_TABLE_SIZE);
+  }
+
+}
+
+/*
+ * vim: set tabstop=8 softtabstop=2 shiftwidth=2 :
+ */