Commit fc020e44 authored by Manuel Grizonnet's avatar Manuel Grizonnet

ENH:add Gamma Function

parent 36026775
/*=========================================================================
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
......@@ -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 )
/*=========================================================================
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);
}
/*=========================================================================
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;
}
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment