| [opengtl-commits] [238] * merge the C++ written in the OpenGTL library ( and make them C symbols) |
[ Thread Index |
Date Index
| More lists.tuxfamily.org/opengtl-commits Archives
]
Revision: 238
Author: cyrille
Date: 2008-06-26 22:12:48 +0200 (Thu, 26 Jun 2008)
Log Message:
-----------
* merge the C++ written in the OpenGTL library (and make them C symbols)
* add lookup3D functions
* add matrix/vector arithmetic operations
Modified Paths:
--------------
trunk/OpenGTL/OpenCTL/CMakeLists.txt
trunk/OpenGTL/OpenCTL/OpenCTL/CMakeLists.txt
trunk/OpenGTL/OpenCTL/OpenCTL/compiler/Compiler.cpp
Added Paths:
-----------
trunk/OpenGTL/OpenCTL/OpenCTL/StdLibFunctions.cpp
trunk/OpenGTL/OpenCTL/OpenCTL/ctlstdlib.ctl
Removed Paths:
-------------
trunk/OpenGTL/OpenCTL/CtlStdLib/
Modified: trunk/OpenGTL/OpenCTL/CMakeLists.txt
===================================================================
--- trunk/OpenGTL/OpenCTL/CMakeLists.txt 2008-06-26 20:09:43 UTC (rev 237)
+++ trunk/OpenGTL/OpenCTL/CMakeLists.txt 2008-06-26 20:12:48 UTC (rev 238)
@@ -10,7 +10,6 @@
set(OPENCTL_LIB_SOVERSION ${OPENGTL_LIB_SOVERSION})
add_subdirectory(OpenCTL)
-add_subdirectory(CtlStdLib)
add_subdirectory(tools)
# After tool to get CTLC variable
Modified: trunk/OpenGTL/OpenCTL/OpenCTL/CMakeLists.txt
===================================================================
--- trunk/OpenGTL/OpenCTL/OpenCTL/CMakeLists.txt 2008-06-26 20:09:43 UTC (rev 237)
+++ trunk/OpenGTL/OpenCTL/OpenCTL/CMakeLists.txt 2008-06-26 20:12:48 UTC (rev 238)
@@ -8,6 +8,7 @@
ModulesManager.cpp
Program.cpp
# Internal files
+ StdLibFunctions.cpp
compiler/Compiler.cpp
compiler/LexerNG.cpp
compiler/ParserNG.cpp
@@ -22,14 +23,10 @@
add_definitions( -DCOUMPONENT_NAME=\"\\\"OpenCTL\\\"\" )
if (NOT APPLE)
- add_definitions( -D_OPENCTL_CTL_STD_LIB_BUILD_DIR_=\\"${CMAKE_CURRENT_BINARY_DIR}/../CtlStdLib/\\")
- add_definitions( -D_OPENCTL_LIB_INSTALL_=\\"${CMAKE_INSTALL_PREFIX}/lib/\\")
- add_definitions( -D_OPENCTL_CTL_STD_LIB_SRC_DIR_=\\"${CMAKE_CURRENT_SOURCE_DIR}/../CtlStdLib/\\")
+ add_definitions( -D_OPENCTL_CTL_STD_LIB_SRC_DIR_=\\"${CMAKE_CURRENT_SOURCE_DIR}/\\")
add_definitions( -D_OPENCTL_CTL_SHARE_DIR_=\\"${SHARE_INSTALL_DIR}/ctl\\" )
else (NOT APPLE)
- add_definitions( -D_OPENCTL_CTL_STD_LIB_BUILD_DIR_='\"${CMAKE_CURRENT_BINARY_DIR}/../CtlStdLib/\"')
- add_definitions( -D_OPENCTL_LIB_INSTALL_='\"${CMAKE_INSTALL_PREFIX}/lib/\"')
- add_definitions( -D_OPENCTL_CTL_STD_LIB_SRC_DIR_='\"${CMAKE_CURRENT_SOURCE_DIR}/../CtlStdLib/\"')
+ add_definitions( -D_OPENCTL_CTL_STD_LIB_SRC_DIR_='\"${CMAKE_CURRENT_SOURCE_DIR}/\"')
add_definitions( -D_OPENCTL_CTL_SHARE_DIR_='\"${SHARE_INSTALL_DIR}/ctl\"' )
endif(NOT APPLE)
@@ -45,3 +42,5 @@
configure_file("OpenCTL.pc.cmake" "${CMAKE_CURRENT_BINARY_DIR}/OpenCTL.pc" @ONLY)
install(FILES ${CMAKE_CURRENT_BINARY_DIR}/OpenCTL.pc DESTINATION ${LIB_INSTALL_DIR}/pkgconfig)
+
+install(FILES ctlstdlib.ctl DESTINATION ${SHARE_INSTALL_DIR}/ctl )
Copied: trunk/OpenGTL/OpenCTL/OpenCTL/StdLibFunctions.cpp (from rev 183, trunk/OpenGTL/OpenCTL/CtlStdLib/Print.cpp)
===================================================================
--- trunk/OpenGTL/OpenCTL/OpenCTL/StdLibFunctions.cpp (rev 0)
+++ trunk/OpenGTL/OpenCTL/OpenCTL/StdLibFunctions.cpp 2008-06-26 20:12:48 UTC (rev 238)
@@ -0,0 +1,93 @@
+/*
+ * Copyright (c) 2008 Cyrille Berger <cberger@xxxxxxxxxxx>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation;
+ * version 2 of the License.
+ *
+ * This library 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this library; see the file COPYING. If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ * Boston, MA 02110-1301, USA.
+ */
+
+#include "GTLCore/String.h"
+
+#include <iostream>
+#include <cstdarg>
+
+#undef assert
+
+extern "C" {
+
+void assert( bool test)
+{
+ if(not test)
+ {
+ throw GTLCore::String("Assert");
+ }
+}
+
+bool isnan_f (float x)
+{
+ union {float f; unsigned int i;} u;
+ u.f = x;
+ return u.i == 0x7fffffff;
+}
+bool isnan_h (float x)
+{
+ union {float f; unsigned int i;} u;
+ u.f = x;
+ return u.i == 0x7fffffff;
+}
+
+
+void print(int nb, ...)
+{
+// std::cout << nb << " : ";
+ va_list argp;
+ va_start(argp, nb);
+ for( int i = 0; i < nb; ++i)
+ {
+ int type = va_arg(argp, int);
+// std::cout << type << " ";
+ switch(type)
+ {
+ case 0:
+ {
+ std::cout << va_arg(argp, int) << " ";
+ }
+ break;
+ case 1:
+ {
+ std::cout << va_arg(argp, double) << " ";
+ }
+ break;
+ case 2:
+ {
+ if( va_arg(argp, int))
+ std::cout << "true ";
+ else
+ std::cout << "false ";
+ }
+ break;
+ case 3:
+ {
+ std::cout << va_arg(argp, char*) << " ";
+ break;
+ }
+ default:
+ std::cout << "Unknown type (" << type << ")";
+ }
+ }
+ std::cout << std::endl;
+ va_end(argp);
+}
+
+}
Modified: trunk/OpenGTL/OpenCTL/OpenCTL/compiler/Compiler.cpp
===================================================================
--- trunk/OpenGTL/OpenCTL/OpenCTL/compiler/Compiler.cpp 2008-06-26 20:09:43 UTC (rev 237)
+++ trunk/OpenGTL/OpenCTL/OpenCTL/compiler/Compiler.cpp 2008-06-26 20:12:48 UTC (rev 238)
@@ -132,37 +132,6 @@
d->errorMessage = "";
setModuleName( moduleName );
- // Load the standard library
- GTLCore::String Error;
- bool notLoaded;
- if((notLoaded = llvm::sys::DynamicLibrary::LoadLibraryPermanently("libctlstdlib.so", &Error))) {
- OCTL_DEBUG( "Error opening 'libctlstdlib.so': " << Error );
- appendError( GTLCore::ErrorMessage(Error));
- }
-#ifdef _OPENCTL_CTL_STD_LIB_BUILD_DIR_
- if( notLoaded)
- {
- if((notLoaded = llvm::sys::DynamicLibrary::LoadLibraryPermanently((GTLCore::String(_OPENCTL_CTL_STD_LIB_BUILD_DIR_) + "libctlstdlib.so").c_str(), &Error))) {
- OCTL_DEBUG( "Error opening 'libctlstdlib.so': " << Error );
- appendError( GTLCore::ErrorMessage(Error));
- }
- }
-#endif
-#ifdef _OPENCTL_CTL_STD_LIB_BUILD_DIR_
- if(notLoaded)
- {
- if((notLoaded = llvm::sys::DynamicLibrary::LoadLibraryPermanently((GTLCore::String(_OPENCTL_LIB_INSTALL_) + "libctlstdlib.so").c_str(), &Error))) {
- OCTL_DEBUG( "Error opening 'libctlstdlib.so': " << Error );
- appendError( GTLCore::ErrorMessage(Error));
- }
- }
-#endif
-
- if( notLoaded )
- {
- return 0;
- }
-
// Create LLVM module
d->module = moduleData->llvmModule();
d->moduleData = moduleData;
@@ -170,8 +139,9 @@
// Create Standard Library functions
// CtlStdLib functions (except print which get a special treatement)
- createStdLibFunction( "assert", "_Z6assertb", GTLCore::Type::Void, 1, GTLCore::Type::Boolean);
- createStdLibFunction( "isnan_f", "_Z7isnan_ff", GTLCore::Type::Boolean, 1, GTLCore::Type::Float);
+ createStdLibFunction( "assert", "assert", GTLCore::Type::Void, 1, GTLCore::Type::Boolean);
+ createStdLibFunction( "isnan_f", "isnan_f", GTLCore::Type::Boolean, 1, GTLCore::Type::Float);
+ createStdLibFunction( "isnan_h", "isnan_h", GTLCore::Type::Boolean, 1, GTLCore::Type::Float);
// C Math functions
createStdLibFunction( "acos", "acosf", GTLCore::Type::Float, 1, GTLCore::Type::Float);
createStdLibFunction( "asin", "asinf", GTLCore::Type::Float, 1, GTLCore::Type::Float);
Copied: trunk/OpenGTL/OpenCTL/OpenCTL/ctlstdlib.ctl (from rev 226, trunk/OpenGTL/OpenCTL/CtlStdLib/ctlstdlib.ctl)
===================================================================
--- trunk/OpenGTL/OpenCTL/OpenCTL/ctlstdlib.ctl (rev 0)
+++ trunk/OpenGTL/OpenCTL/OpenCTL/ctlstdlib.ctl 2008-06-26 20:12:48 UTC (rev 238)
@@ -0,0 +1,694 @@
+/*
+ * Copyright (c) 2008 Cyrille Berger <cberger@xxxxxxxxxxx>
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation;
+ * version 2 of the License.
+ *
+ * This library 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
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public License
+ * along with this library; see the file COPYING. If not, write to
+ * the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
+ * Boston, MA 02110-1301, USA.
+ */
+
+// Float operation
+
+bool isfinite_f (float x)
+{
+ return x != FLT_POS_INF && x != FLT_NEG_INF && x != FLT_NAN;
+}
+
+bool isfinite_h (half x)
+{
+ return x != HALF_POS_INF && x != HALF_NEG_INF && x != HALF_NAN;
+}
+
+bool isnormal_f (float x)
+{
+ return isfinite_f(x) && x != 0.0;
+}
+bool isnormal_h (half x)
+{
+ return isfinite_h(x) && x != 0.0;
+}
+
+bool isinf_f (float x)
+{
+ return x == FLT_POS_INF || x == FLT_NEG_INF;
+}
+bool isinf_h (half x)
+{
+ return x == HALF_POS_INF || x != HALF_NEG_INF;
+}
+
+float pow10(float x)
+{
+ return pow(10,x);
+}
+
+half pow10_h(float x)
+{
+ return pow(10,x);
+}
+
+float hypot(float x, float y)
+{
+ return sqrt(x*x + y*y);
+}
+
+// misc mathematic function
+
+float min( float f1, float f2) // WARNING not mentioned in the Spec !
+{
+ if( f1 < f2)
+ {
+ return f1;
+ } else {
+ return f2;
+ }
+}
+
+float max( float f1, float f2) // WARNING not mentioned in the Spec !
+{
+ if( f1 > f2)
+ {
+ return f1;
+ } else {
+ return f2;
+ }
+}
+
+// lookup
+
+float lookup1D(float table[], float pMin, float pMax, float p)
+{
+ if( p < pMin ) return table[ 0 ];
+ if( p > pMax ) return table[ table.size - 1 ];
+ float t = (p - pMin) / (pMax - pMin ) * (table.size - 1 );
+ int i = floor( t );
+ float s = t - i;
+ return table[i] * ( 1 - s ) + table[i+1] * s;
+}
+
+float lookupCubic1D (float table[], float pMin, float pMax, float p)
+{
+ if( p < pMin ) return table[ 0 ];
+ if( p > pMax ) return table[ table.size - 1 ];
+
+ float t = (p - pMin) / (pMax - pMin) * (table.size-1);
+ int i = floor (t);
+ float s = t - i;
+ float m0;
+ float m1;
+ if( i > 0 )
+ {
+ m0 = (table[i+1] - table[i-1]) / 2;
+ }
+ if( i < table.size-2 )
+ {
+ m1 = (table[i+2] - table[i]) / 2;
+ }
+ if( i == 0) {
+ m0 = (3 * table[i+1] - table[i] - m1);
+ }
+ if( i == table.size-2 )
+ {
+ m1 = (3 * table[i+1] - table[i] - m0);
+ }
+ return table[i] * (2 * s*s*s - 3 * s*s + 1) + m0 * (s*s*s - 2 * s*s + s) + table[i+1] * (-2 * s*s*s + 3 * s*s) + m1 * (s*s*s - s*s);
+}
+
+float interpolate1D (float table[][2], float p)
+{
+ if( p < table[0][0] ) return table[0][1];
+ if( p > table[table.size-1][0] ) return table[table.size-1][1];
+
+ for( int i = 0; i < table.size - 1; ++i )
+ {
+ if( table[i][0] <= p && p < table[i+1][0] )
+ {
+ float s = (p - table[i][0]) / (table[i+1][0] - table[i][0]);
+ return table[i][1] * ( 1 - s ) + table[i+1][1] * s;
+ }
+ }
+ return 0.0;
+}
+
+float interpolateCubic1D (float table[][2], float p)
+{
+ if( p < table[0][0] ) return table[0][1];
+ if( p > table[table.size-1][0] ) return table[table.size-1][1];
+
+ for( int i = 0; i < table.size - 1; ++i )
+ {
+ if( table[i][0] <= p && p < table[i+1][0] )
+ {
+ float s = (p - table[i][0]) / (table[i+1][0] - table[i][0]);
+ float dx0 = (table[i][0] - table[i-1][0]);
+ float dx1 = (table[i+1][0] - table[i][0]);
+ float dx2 = (table[i+2][0] - table[i+1][0]);
+ float dy0 = (table[i][1] - table[i-1][1]);
+ float dy1 = (table[i+1][1] - table[i][1]);
+ float dy2 = (table[i+2][1] - table[i+1][1]);
+
+ float m0;
+ float m1;
+ if( i > 0 )
+ {
+ m0 = (dy1 + dx1 * dy0 / dx0) / 2;
+ }
+ if( i < table.size-2 )
+ {
+ m1 = (dy1 + dx1 * dy2 / dx2) / 2;
+ }
+ if( i == 0) {
+ m0 = (3 * dy1 - m1) / 2;
+ }
+ if( i == table.size-2 )
+ {
+ m1 = (3 * dy1 - m0) / 2;
+ }
+ return table[i][1] * (2 * s*s*s - 3 * s*s + 1) +
+ m0 * (s*s*s - 2 * s*s + s) +
+ table[i+1][1] * (-2 * s*s*s + 3 * s*s) +
+ m1 * (s*s*s - s*s);
+ }
+ }
+ return 0.0;
+}
+
+float[3] lookup3D_f3
+ (float table[][][][3],
+ float pMin[3],
+ float pMax[3],
+ float p[3])
+{
+ float result[3];
+ int iMax = table.size - 1;
+ int jMax = table[0].size - 1;
+ int kMax = table[0][0].size - 1;
+ float q[3];
+ q[0] = max (pMin[0], min (pMax[0], p[0]));
+ q[1] = max (pMin[1], min (pMax[1], p[1]));
+ q[2] = max (pMin[2], min (pMax[2], p[2]));
+
+ float ti = (p[0] - pMin[0]) / (pMax[0] - pMin[0]) * iMax;
+ int i = floor (ti);
+ float si = ti - i;
+ float tj = (p[1] - pMin[1]) / (pMax[1] - pMin[1]) * jMax;
+ int j = floor (tj);
+ float sj = tj - j;
+ float tk = (p[2] - pMin[2]) / (pMax[2] - pMin[2]) * kMax;
+ int k = floor (tk);
+ float sk = tk - k;
+
+ for( int l = 0; l < 3; ++l)
+ {
+ result[ l ] = ((table[i][j][k][l] * (1-si) + table[i+1][j][k][l] * si) * (1-sj)
+ + (table[i][j+1][k][l] * (1-si) + table[i+1][j+1][k][l] * si) * sj ) * (1-sk)
+ + ((table[i][j][k+1][l] * (1-si) + table[i+1][j][k+1][l] * si) * (1-sj)
+ + (table[i][j+1][k+1][l] * (1-si) + table[i+1][j+1][k+1][l] * si) * sj ) * sk;
+ }
+
+ return result;
+}
+
+void lookup3D_f
+ (float table[][][][3],
+ float pMin[3],
+ float pMax[3],
+ float p0, float p1, float p2,
+ output float q0, output float q1, output float q2)
+{
+ float p[3];
+ p[0] = p0;
+ p[1] = p1;
+ p[2] = p2;
+ float result[3] = lookup3D_f3( table, pMin, pMax, p );
+ q0 = result[0];
+ q1 = result[1];
+ q2 = result[2];
+}
+
+void lookup3D_h
+ (float table[][][][3],
+ float pMin[3],
+ float pMax[3],
+ half p0, half p1, half p2,
+ output half q0, output half q1, output half q2)
+{
+ lookup3D_f( table, pMin, pMax, p0, p1, p2, q0, q1, q2 );
+}
+
+// Color conversion functions
+
+struct Chromaticities
+{
+ float red[2];
+ float green[2];
+ float blue[2];
+ float white[2];
+};
+
+// float[4][4] RGBtoXYZ(Chromaticities c, float Y)
+// {
+// }
+//
+// float[4][4] XYZtoRGB(Chromaticities c, float Y)
+// {
+// }
+//
+// float[3] XYZtoLuv(float XYZ[3], float XYZn[3]);
+// float[3] LuvtoXYZ(float Luv[3], float XYZn[3]);
+//float[3] XYZtoLab(float XYZ[3], float XYZn[3]);
+// float[3] LabtoXYZ(float Lab[3], float XYZn[3]);
+//
+
+// Vectors/Matrix operations
+
+float[3] mult_f_f3(float f, float x[3])
+{
+ float r[3];
+ for( int k = 0; k < 3; ++k)
+ {
+ r[k] = f * x[k];
+ }
+ return r;
+}
+
+float[3] add_f3_f3(float x[3], float y[3])
+{
+ float r[3];
+ for( int k = 0; k < 3; ++k)
+ {
+ r[k] = x[k] + y[k];
+ }
+ return r;
+}
+
+float[3] sub_f3_f3(float x[3], float y[3])
+{
+ float r[3];
+ for( int k = 0; k < 3; ++k)
+ {
+ r[k] = x[k] - y[k];
+ }
+ return r;
+}
+
+float[3] cross_f3_f3(float x[3], float y[3])
+{
+ float r[3];
+ r[2] = x[0] * y[1] - x[1] * y[0];
+ r[0] = x[1] * y[2] - x[2] * y[1];
+ r[1] = x[2] * y[0] - x[0] * y[2];
+ return r;
+}
+
+float dot_f3_f3(float x[3], float y[3])
+{
+ return x[0] * y[0] + x[1] * y[1] + x[2] * y[2];
+}
+
+float length_f3 (float x[3])
+{
+ return sqrt( x[0] * x[0] + x[1] * x[1] + x[2] * x[2] );
+}
+
+float[3][3] mult_f33_f33 (float A[3][3], float B[3][3])
+{
+ float r[3][3];
+ for( int i = 0; i < 3; ++i)
+ {
+ for( int j = 0; j < 3; ++j)
+ {
+ r[i][j] = 0.0;
+ for( int k = 0; k < 3; ++k)
+ {
+ r[i][j] = r[i][j] + A[i][k] * B[k][j];
+ }
+ }
+ }
+ return r;
+}
+
+float[4][4] mult_f44_f44 (float A[4][4], float B[4][4])
+{
+ float r[4][4];
+ for( int i = 0; i < 4; ++i)
+ {
+ for( int j = 0; j < 4; ++j)
+ {
+ r[i][j] = 0.0;
+ for( int k = 0; k < 4; ++k)
+ {
+ r[i][j] = r[i][j] + A[i][k] * B[k][j];
+ }
+ }
+ }
+ return r;
+}
+
+float[3][3] mult_f_f33 (float f, float A[3][3])
+{
+ float r[3][3];
+ for( int i = 0; i < 3; ++i )
+ {
+ for( int j = 0; j < 3; ++j )
+ {
+ r[i][j] = f * A[i][j];
+ }
+ }
+ return r;
+}
+
+float[4][4] mult_f_f44 (float f, float A[4][4])
+{
+ float r[4][4];
+ for( int i = 0; i < 4; ++i )
+ {
+ for( int j = 0; j < 4; ++j )
+ {
+ r[i][j] = f * A[i][j];
+ }
+ }
+ return r;
+}
+
+float[3][3] add_f33_f33(float A[3][3], float B[3][3])
+{
+ float r[3][3];
+ for( int i = 0; i < 3; ++i )
+ {
+ for( int j = 0; j < 3; ++j )
+ {
+ r[i][j] = A[i][j] + B[i][j];
+ }
+ }
+ return r;
+}
+
+float[4][4] add_f44_f44(float A[4][4], float B[4][4])
+{
+ float r[4][4];
+ for( int i = 0; i < 4; ++i )
+ {
+ for( int j = 0; j < 4; ++j )
+ {
+ r[i][j] = A[i][j] + B[i][j];
+ }
+ }
+ return r;
+}
+
+float[3][3] invert_f33 (float A[3][3])
+{
+ float result[3][3];
+ float det = A[0][0] * A[1][1] * A[2][2]
+ + A[0][1] * A[1][2] * A[2][0]
+ + A[0][2] * A[1][0] * A[2][1]
+ - A[2][0] * A[1][1] * A[0][2]
+ - A[2][1] * A[1][2] * A[0][0]
+ - A[2][2] * A[1][0] * A[0][1];
+ if( det != 0.0 )
+ {
+ result[0][0] = A[1][1] * A[2][2] - A[1][2] * A[2][1];
+ result[0][1] = A[2][1] * A[0][2] - A[2][2] * A[0][1];
+ result[0][2] = A[0][1] * A[1][2] - A[0][2] * A[1][1];
+ result[1][0] = A[2][0] * A[1][2] - A[1][0] * A[2][2];
+ result[1][1] = A[0][0] * A[2][2] - A[2][0] * A[0][2];
+ result[1][2] = A[1][0] * A[0][2] - A[0][0] * A[1][2];
+ result[2][0] = A[1][0] * A[2][1] - A[2][0] * A[1][1];
+ result[2][1] = A[2][0] * A[0][1] - A[0][0] * A[2][1];
+ result[2][2] = A[0][0] * A[1][1] - A[1][0] * A[0][1];
+
+ return mult_f_f33( 1.0 / det, result);
+ }
+ result = { { 1, 0, 0 }, { 0, 1, 0 }, { 0, 0, 1 } };
+ return result;
+}
+
+
+//---------- Special copyright notice for function invert_f44 -----------//
+
+///////////////////////////////////////////////////////////////////////////
+//
+// Copyright (c) 2008, Cyrille Berger <cberger@xxxxxxxxxxx>
+// Copyright (c) 2002, Industrial Light & Magic, a division of Lucas
+// Digital Ltd. LLC
+//
+// All rights reserved.
+//
+// Redistribution and use in source and binary forms, with or without
+// modification, are permitted provided that the following conditions are
+// met:
+// * Redistributions of source code must retain the above copyright
+// notice, this list of conditions and the following disclaimer.
+// * Redistributions in binary form must reproduce the above
+// copyright notice, this list of conditions and the following disclaimer
+// in the documentation and/or other materials provided with the
+// distribution.
+// * Neither the name of Industrial Light & Magic nor the names of
+// its contributors may be used to endorse or promote products derived
+// from this software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
+// "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
+// LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
+// A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
+// OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
+// SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
+// LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
+// DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
+// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
+// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+///////////////////////////////////////////////////////////////////////////
+
+float[4][4] invert_f44 (float A[4][4])
+{
+
+
+ if (A[0][3] != 0 || A[1][3] != 0 || A[2][3] != 0 || A[3][3] != 1)
+ {
+ int i;
+ int j;
+ int k;
+ float s[4][4] = {{1,0,0,0},{0,1,0,0},{0,0,1,0},{0,0,0,1}};
+ float t[4][4] = A;
+
+ // Forward elimination
+
+ for (i = 0; i < 3 ; ++i)
+ {
+ int pivot = i;
+
+ float pivotsize = t[i][i];
+
+ if (pivotsize < 0)
+ {
+ pivotsize = -pivotsize;
+ }
+
+ for (j = i + 1; j < 4; ++j)
+ {
+ float tmp = t[j][i];
+
+ if (tmp < 0)
+ tmp = -tmp;
+
+ if (tmp > pivotsize)
+ {
+ pivot = j;
+ pivotsize = tmp;
+ }
+ }
+ if (pivotsize == 0)
+ {
+ s = { { 1, 0, 0, 0 }, { 0, 1, 0, 0 }, { 0, 0, 1, 0 }, { 0, 0, 0, 1 } };
+ return s;
+ }
+
+ if (pivot != i)
+ {
+ for (j = 0; j < 4; ++j)
+ {
+ float tmp = t[i][j];
+ t[i][j] = t[pivot][j];
+ t[pivot][j] = tmp;
+
+ tmp = s[i][j];
+ s[i][j] = s[pivot][j];
+ s[pivot][j] = tmp;
+ }
+ }
+
+ for (j = i + 1; j < 4; ++j)
+ {
+ float f = t[j][i] / t[i][i];
+ for (k = 0; k < 4; ++k)
+ {
+ t[j][k] = t[j][k] - f * t[i][k];
+ s[j][k] = s[j][k] - f * s[i][k];
+ }
+ }
+ }
+ // Backward substitution
+ for (i = 3; i >= 0; --i)
+ {
+ float f = t[i][i];
+ if ( f == 0)
+ {
+ s = { { 1, 0, 0, 0 }, { 0, 1, 0, 0 }, { 0, 0, 1, 0 }, { 0, 0, 0, 1 } };
+ return s;
+ }
+ for (j = 0; j < 4; ++j)
+ {
+ t[i][j] = t[i][j] / f;
+ s[i][j] = s[i][j] / f;
+ }
+
+ for (j = 0; j < i; ++j)
+ {
+ f = t[j][i];
+
+ for (k = 0; k < 4; ++k)
+ {
+ t[j][k] = t[j][k] - f * t[i][k];
+ s[j][k] = s[j][k] - f * s[i][k];
+ }
+ }
+ }
+ return s;
+ } else {
+
+ float result[4][4];
+ result[0][0] = A[1][1] * A[2][2] - A[2][1] * A[1][2];
+ result[0][1] = A[2][1] * A[0][2] - A[0][1] * A[2][2];
+ result[0][2] = A[0][1] * A[1][2] - A[1][1] * A[0][2];
+ result[0][3] = 0;
+
+ result[1][0] = A[2][0] * A[1][2] - A[1][0] * A[2][2];
+ result[1][1] = A[0][0] * A[2][2] - A[2][0] * A[0][2];
+ result[1][2] = A[1][0] * A[0][2] - A[0][0] * A[1][2];
+ result[1][3] = 0;
+
+ result[0][0] = A[1][0] * A[2][1] - A[2][0] * A[1][1];
+ result[0][1] = A[2][0] * A[0][1] - A[0][0] * A[2][1];
+ result[0][2] = A[0][0] * A[1][1] - A[1][0] * A[0][1];
+ result[0][3] = 0;
+
+ result[3][0] = 0;
+ result[3][1] = 0;
+ result[3][2] = 0;
+ result[3][3] = 1;
+
+ float r = A[0][0] * result[0][0] + A[0][1] * result[1][0] + A[0][2] * result[2][0];
+
+ if(fabs(r) >= 1)
+ {
+ for(int i = 0; i < 3; ++i)
+ {
+ for(int j = 0; j < 3; ++j)
+ {
+ result[i][j] = result[i][j] / r;
+ }
+ }
+ }
+ else
+ {
+ float mr = fabs(r) / FLT_MIN;
+
+ for (int i = 0; i < 3; ++i)
+ {
+ for (int j = 0; j < 3; ++j)
+ {
+ if (mr > fabs(result[i][j]))
+ {
+ result[i][j] = result[i][j] / r;
+ }
+ else
+ {
+ result = { { 1, 0, 0, 0 }, { 0, 1, 0, 0 }, { 0, 0, 1, 0 }, { 0, 0, 0, 1 } };
+ return result;
+ }
+ }
+ }
+ }
+ result[3][0] = -A[3][0] * result[0][0] - A[3][1] * result[1][0] - A[3][2] * result[2][0];
+ result[3][1] = -A[3][0] * result[0][1] - A[3][1] * result[1][1] - A[3][2] * result[2][1];
+ result[3][2] = -A[3][0] * result[0][2] - A[3][1] * result[1][2] - A[3][2] * result[2][2];
+ return result;
+ }
+ float result[4][4];
+ return result;
+}
+
+// End of special copyright notice
+
+float[3][3] transpose_f33 (float A[3][3])
+{
+ float r[3][3];
+ for( int i = 0; i < 3; ++i)
+ {
+ for( int j = 0; j < 3; ++j)
+ {
+ r[i][j] = A[j][i];
+ }
+ }
+ return r;
+}
+
+float[4][4] transpose_f44 (float A[4][4])
+{
+ float r[4][4];
+ for( int i = 0; i < 4; ++i)
+ {
+ for( int j = 0; j < 4; ++j)
+ {
+ r[i][j] = A[j][i];
+ }
+ }
+ return r;
+}
+
+float[3] mult_f3_f33 (float x[3], float A[3][3])
+{
+ float r[3];
+ for( int i = 0; i < 3; ++i)
+ {
+ r[i] = 0.0;
+ for( int j = 0; j < 3; ++j)
+ {
+ r[i] = r[i] + x[j] * A[j][i];
+ }
+ }
+ return r;
+}
+
+float[3] mult_f3_f44 (float x[3], float A[4][4])
+{
+ float r[3];
+ for( int i = 0; i < 3; ++i)
+ {
+ r[i] = 0.0;
+ for( int j = 0; j < 3; ++j)
+ {
+ r[i] = r[i] + x[j] * A[j][i];
+ }
+ r[i] = r[i] + A[3][i];
+ }
+ float s = 1.0 / (x[0] * A[0][3] + x[1] * A[1][3] + x[2] * A[2][3] + A[3][3]);
+ for( int k = 0; k < 3; ++k)
+ {
+ r[k] = r[k] * s;
+ }
+ return r;
+}