From fc020e4418f0260a38f59573624058f227ed47c8 Mon Sep 17 00:00:00 2001
From: Manuel Grizonnet <manuel.grizonnet@orfeo-toolbox.org>
Date: Sat, 7 May 2011 14:53:00 +0200
Subject: [PATCH] ENH:add Gamma Function
---
Code/Common/otbGamma.h | 107 +++++++++++++++++++++++
Testing/Code/Common/CMakeLists.txt | 15 +++-
Testing/Code/Common/otbCommonTests14.cxx | 30 +++++++
Testing/Code/Common/otbGammaTest.cxx | 37 ++++++++
4 files changed, 188 insertions(+), 1 deletion(-)
create mode 100644 Code/Common/otbGamma.h
create mode 100644 Testing/Code/Common/otbCommonTests14.cxx
create mode 100644 Testing/Code/Common/otbGammaTest.cxx
diff --git a/Code/Common/otbGamma.h b/Code/Common/otbGamma.h
new file mode 100644
index 0000000000..cbc0497960
--- /dev/null
+++ b/Code/Common/otbGamma.h
@@ -0,0 +1,107 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+#ifndef __otbGamma_h
+#define __otbGamma_h
+
+#ifdef _MSC_VER
+#pragma warning ( disable : 4290 )
+#endif
+
+#include "otbMath.h"
+
+namespace otb
+{
+/** \class Gamma
+ * \brief Gamma function (necessary for Fisher distribution
+ *
+ */
+class Gamma {
+public:
+ double gamma(double x)
+ {
+ int i,k,m;
+ double ga,gr,r,z;
+
+ static double g[] = {
+ 1.0,
+ 0.5772156649015329,
+ -0.6558780715202538,
+ -0.420026350340952e-1,
+ 0.1665386113822915,
+ -0.421977345555443e-1,
+ -0.9621971527877e-2,
+ 0.7218943246663e-2,
+ -0.11651675918591e-2,
+ -0.2152416741149e-3,
+ 0.1280502823882e-3,
+ -0.201348547807e-4,
+ -0.12504934821e-5,
+ 0.1133027232e-5,
+ -0.2056338417e-6,
+ 0.6116095e-8,
+ 0.50020075e-8,
+ -0.11812746e-8,
+ 0.1043427e-9,
+ 0.77823e-11,
+ -0.36968e-11,
+ 0.51e-12,
+ -0.206e-13,
+ -0.54e-14,
+ 0.14e-14};
+
+ if (x > 171.0) return 1e308; // This value is an overflow flag.
+ if (x == (int)x) {
+ if (x > 0.0) {
+ ga = 1.0; // use factorial
+ for (i=2;i<x;i++) {
+ ga *= i;
+ }
+ }
+ else
+ ga = 1e308;
+ }
+ else {
+ if (fabs(x) > 1.0) {
+ z = fabs(x);
+ m = (int)z;
+ r = 1.0;
+ for (k=1;k<=m;k++) {
+ r *= (z-k);
+ }
+ z -= m;
+ }
+ else
+ z = x;
+ gr = g[24];
+ for (k=23;k>=0;k--) {
+ gr = gr*z+g[k];
+ }
+ ga = 1.0/(gr*z);
+ if (fabs(x) > 1.0) {
+ ga *= r;
+ if (x < 0.0) {
+ ga = -M_PI/(x*ga*sin(M_PI*x));
+ }
+ }
+ }
+ return ga;
+ }
+};
+} // end namespace otb
+
+#endif
diff --git a/Testing/Code/Common/CMakeLists.txt b/Testing/Code/Common/CMakeLists.txt
index d7ef54353b..887f0a9f29 100644
--- a/Testing/Code/Common/CMakeLists.txt
+++ b/Testing/Code/Common/CMakeLists.txt
@@ -34,6 +34,7 @@ ENDIF(OTB_USE_PQXX)
SET(COMMON_TESTS11 ${CXX_TEST_PATH}/otbCommonTests11)
SET(COMMON_TESTS12 ${CXX_TEST_PATH}/otbCommonTests12)
SET(COMMON_TESTS13 ${CXX_TEST_PATH}/otbCommonTests13)
+SET(COMMON_TESTS14 ${CXX_TEST_PATH}/otbCommonTests14)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbCommonTests1 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -1043,7 +1044,13 @@ ADD_TEST(coTvRAMDrivenTiledStreamingManager ${COMMON_TESTS13}
${TEMP}/coTvRAMDrivenTiledStreamingManager.txt
)
-
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ otbCommonTests14 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# Gamma function tests #
+ADD_TEST(coTvGamma ${COMMON_TESTS14}
+ otbGammaTest
+)
# ------- Fichiers sources CXX -----------------------------------
SET(BasicCommon_SRCS1
@@ -1209,6 +1216,11 @@ otbPipelineMemoryPrintCalculatorTest.cxx
otbStreamingManager.cxx
)
+SET(BasicCommon_SRCS14
+otbCommonTests14.cxx
+otbGammaTest.cxx
+)
+
OTB_ADD_EXECUTABLE(otbCommonTests1 "${BasicCommon_SRCS1}" "OTBIO;OTBTesting")
OTB_ADD_EXECUTABLE(otbCommonTests2 "${BasicCommon_SRCS2}" "OTBIO;OTBTesting")
@@ -1231,5 +1243,6 @@ TARGET_LINK_LIBRARIES(otbCommonTests11 OTBIO OTBTesting itkvcl)
OTB_ADD_EXECUTABLE(otbCommonTests12 "${BasicCommon_SRCS12}" "OTBIO;OTBTesting")
OTB_ADD_EXECUTABLE(otbCommonTests13 "${BasicCommon_SRCS13}" "OTBIO;OTBTesting")
+OTB_ADD_EXECUTABLE(otbCommonTests14 "${BasicCommon_SRCS14}" "OTBIO;OTBTesting")
ENDIF( NOT OTB_DISABLE_CXX_TESTING AND BUILD_TESTING )
diff --git a/Testing/Code/Common/otbCommonTests14.cxx b/Testing/Code/Common/otbCommonTests14.cxx
new file mode 100644
index 0000000000..a26f15d108
--- /dev/null
+++ b/Testing/Code/Common/otbCommonTests14.cxx
@@ -0,0 +1,30 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+
+// this file defines the otbCommonTest for the test driver
+// and all it expects is that you have a function called RegisterTests
+#if defined(_MSC_VER)
+#pragma warning ( disable : 4786 )
+#endif
+
+#include "otbTestMain.h"
+
+void RegisterTests()
+{
+ REGISTER_TEST(otbGammaTest);
+}
diff --git a/Testing/Code/Common/otbGammaTest.cxx b/Testing/Code/Common/otbGammaTest.cxx
new file mode 100644
index 0000000000..627e039add
--- /dev/null
+++ b/Testing/Code/Common/otbGammaTest.cxx
@@ -0,0 +1,37 @@
+/*=========================================================================
+
+ Program: ORFEO Toolbox
+ Language: C++
+ Date: $Date$
+ Version: $Revision$
+
+
+ Copyright (c) Centre National d'Etudes Spatiales. All rights reserved.
+ See OTBCopyright.txt for details.
+
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notices for more information.
+
+=========================================================================*/
+
+#include "otbGamma.h"
+#include "otbMacro.h"
+#include <iostream>
+
+
+int otbGammaTest(int argc, char * argv[])
+{
+ typedef otb::Gamma GammaFunctionType;
+
+ //Instantiation
+ GammaFunctionType * gam = new GammaFunctionType();
+ const double epsilon = 0.0000000001;
+
+ if ( vcl_abs(gam->gamma(1) - 1) > epsilon ) return EXIT_FAILURE;
+ if ( vcl_abs(gam->gamma(0.5) - 1.77245385091) > epsilon ) return EXIT_FAILURE;
+ if ( vcl_abs(gam->gamma(4) - 6) > epsilon ) return EXIT_FAILURE;
+
+ return EXIT_SUCCESS;
+}
--
GitLab