gmpPiecewiseCurve.h

Go to the documentation of this file.

/************************************************************************
**
** Copyright (C) 2014 by Carlos Augusto Teixera Mendes
** All rights reserved.
**
** This file is part of the "GeMA" software. It's use should respect
** the terms in the license agreement that can be found together
** with this source code.
** It is provided AS IS, with NO WARRANTY OF ANY KIND,
** INCLUDING THE WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR
** A PARTICULAR PURPOSE.
**
************************************************************************/
 
#include <gmTrace.h>
#include <gmLog.h>
#include <assert.h>
#include <cmath>
#include <limits>
#include <gmMathUtils.h>
 
#ifndef _GEMA_PIECEWISECURVE_H_
#define _GEMA_PIECEWISECURVE_H_
 
#ifndef EPSILON
#define EPSILON      2.2204460492503131e-016
#endif
 
namespace GmpPiecewiseCurve
{
        // -----------------------------------------------------------------------------------------------------------
        int findInterval(double value, GmCRMatrix& curve, int col)
        {
                S_TRACE();
                // Check if the given point is located inside the interval defined by the curve
                if (curve(0, col) < curve(curve.n_rows - 1, col))
                {
                        /* Monotonically increasing curve */
                        if (value < curve(0, col))
                                return 1;
                        if (value > curve(curve.n_rows - 1, col))
                                return curve.n_rows - 1;
                        int i = 1;
                        while (curve(i, col) < value)
                                i++;
                        return i;
                }
                else
                {
                        /* Monotonically decreasing curve */
                        if (value > curve(0, col))
                                return 1;
                        if (value < curve(curve.n_rows - 1, col))
                                return curve.n_rows - 1;
                        int i = 1;
                        while (curve(i, col) > value)
                                i++;
                        return i;
                }
 
                return -1;
        }
        
        // -----------------------------------------------------------------------------------------------------------
        double interpolate(double x, GmCRMatrix& curve, const GmLogCategory& logger)
        {
                S_TRACE();
                // Find the interval where the x value is located
                int i = findInterval(x, curve, 0);
                if (i < 1)
                {
                        gmErrorMsg(logger, "GmpPiecewiseCurve: Error 'x' value is out of the curve range.");
                        return std::numeric_limits<double>::quiet_NaN();
                }
                // Interpolate the y value
                double dx = curve(i, 0) - curve(i - 1, 0);
                double dy = curve(i, 1) - curve(i - 1, 1);
                if (std::abs(dx) < EPSILON)
                {
                        gmErrorMsg(logger, "GmpPiecewiseCurve: Error the given curve has repeated x values.");
                        return std::numeric_limits<double>::quiet_NaN();
                }
                double y = (dy / dx) * (x - curve(i - 1, 0)) + curve(i - 1, 1);
                return y;
                }
 
        // -----------------------------------------------------------------------------------------------------------
        double derivativeGivenY(double y, GmCRMatrix& curve, const GmLogCategory& logger)
        {
                S_TRACE();
 
                // Find the interval where the y value is located
                int i = findInterval(y, curve, 1);
                if (i < 1)
                {
                        gmErrorMsg(logger, "GmpPiecewiseCurve: Error 'y' value is out of the curve range.");
                        return std::numeric_limits<double>::quiet_NaN();
                }
 
                // Interpolate the y value
                double dx = curve(i, 0) - curve(i - 1, 0);
                double dy = curve(i, 1) - curve(i - 1, 1);
                if (std::abs(dx) > EPSILON)
                {
                        return dy / dx;
                }
                else
                {
                        return GmMathUtils::sign(dy) / EPSILON;
                }
        }
};
#endif