KrotosIntrepolators.cpp

Go to the documentation of this file.

namespace krotos
{
    float ShapePreservingCubicInterpolator::interpolate(const std::vector<float>& xPositions,
                                                        const std::vector<float>& yPositions, float x)
    {
        jassert(xPositions.size() == yPositions.size());
 
        const int numPoints = static_cast<int>(xPositions.size());
 
        if (numPoints < 2)
        {
            jassertfalse; // Not enough points for interpolation
            return 0.0f;
        }
 
        int index = 0;
 
        // Find the index of the nearest point to the left of x
        for (int i = 0; i < numPoints - 1; ++i)
        {
            if (xPositions[i + 1] >= x)
            {
                index = i;
                break;
            }
        }
 
        const float x0 = xPositions[index];
        const float x1 = xPositions[index + 1];
        const float y0 = yPositions[index];
        const float y1 = yPositions[index + 1];
 
        const float dx = x1 - x0;
        const float dy = y1 - y0;
        const float h0 = (index > 0) ? x0 - xPositions[index - 1] : dx;
        const float h1 = (index < numPoints - 2) ? xPositions[index + 2] - x1 : dx;
 
        const float t = (x - x0) / dx;
        const float t2 = t * t;
 
        const float p = (1.0f - t) * y0 + t * y1;
 
        if (dy * h0 > 0.0f && dy * h1 > 0.0f)
        {
            const float a = (2.0f * dy - h0 * (y1 - y0) / dx) / h0;
            const float b = (2.0f * dy - h1 * (y1 - y0) / dx) / h1;
 
            return t * (1.0f - t) * (a * h0 + b * h1) + p;
        }
        else
        {
            return p;
        }
    }
 
    float ShapePreservingCubicInterpolator::interpolateEigen(const Eigen::VectorXf& xPositions,
                                                             const Eigen::VectorXf& yPositions, float x)
    {
        jassert(xPositions.size() == yPositions.size());
 
        const int numPoints = static_cast<int>(xPositions.size());
 
        if (numPoints < 2)
        {
            jassertfalse; // Not enough points for interpolation
            return 0.0f;
        }
 
        int index = 0;
 
        // Find the index of the nearest point to the left of x
        for (int i = 0; i < numPoints - 1; ++i)
        {
            if (xPositions[i + 1] >= x)
            {
                index = i;
                break;
            }
        }
 
        const float x0 = xPositions[index];
        const float x1 = xPositions[index + 1];
        const float y0 = yPositions[index];
        const float y1 = yPositions[index + 1];
 
        const float dx = x1 - x0;
        const float dy = y1 - y0;
        const float h0 = (index > 0) ? x0 - xPositions[index - 1] : dx;
        const float h1 = (index < numPoints - 2) ? xPositions[index + 2] - x1 : dx;
 
        const float t = (x - x0) / dx;
        const float t2 = t * t;
 
        const float p = (1.0f - t) * y0 + t * y1;
 
        if (dy * h0 > 0.0f && dy * h1 > 0.0f)
        {
            const float a = (2.0f * dy - h0 * (y1 - y0) / dx) / h0;
            const float b = (2.0f * dy - h1 * (y1 - y0) / dx) / h1;
 
            return t * (1.0f - t) * (a * h0 + b * h1) + p;
        }
        else
        {
            return p;
        }
    }
 
    std::vector<std::vector<float>> ShapePreservingCubicInterpolator::interpolate2D(
            const std::vector<float>& xValues, const std::vector<std::vector<float>>& yValues,
            const std::vector<float>& vals)
    {
        // Check if input dimensions are valid
        if (xValues.size() != yValues[0].size() || yValues[0].size() == 0)
        {
            // Invalid input dimensions
            jassertfalse;
        }
 
        std::vector<std::vector<float>> resultMatrix;
 
        for (int row = 0; row < yValues.size(); ++row)
        {
            // Interpolate using Shape Preserving Cubic Interpolation
            std::vector<float> interpolatedValues = interpolateVector(xValues, yValues[row], vals);
 
            resultMatrix.push_back(interpolatedValues);
        }
 
        return resultMatrix;
    }
 
    Eigen::MatrixXf ShapePreservingCubicInterpolator::interpolate2DEigen(const Eigen::VectorXf& xValues,
                                                                         const Eigen::MatrixXf yValues,
                                                                         const Eigen::VectorXf& vals)
    {
        // Check if input dimensions are valid
        if (xValues.size() != yValues.rows() || yValues.rows() == 0)
        {
            // Invalid input dimensions
            jassertfalse;
        }
 
        Eigen::MatrixXf resultMatrix(vals.size(), yValues.cols());
 
        for (int col = 0; col < yValues.cols(); ++col)
        {
            // Interpolate using Shape Preserving Cubic Interpolation
            Eigen::VectorXf interpolatedValues = interpolateVectorEigen(xValues, yValues.col(col), vals);
 
            resultMatrix.col(col) = interpolatedValues;
        }
 
        return resultMatrix;
    }
 
    std::vector<float> ShapePreservingCubicInterpolator::interpolateVector(const std::vector<float>& xPositions,
                                                                           const std::vector<float>& yPositions,
                                                                           const std::vector<float>& xValues)
    {
        jassert(xPositions.size() == yPositions.size());
 
        const int numPoints = static_cast<int>(xPositions.size());
 
        if (numPoints < 2)
        {
            jassertfalse; // Not enough points for interpolation
            return std::vector<float>(xValues.size(), 0.0f);
        }
 
        std::vector<float> interpolatedValues;
        interpolatedValues.reserve(xValues.size());
 
        for (float x : xValues)
        {
            int index = 0;
 
            // Find the index of the nearest point to the left of x
            for (int i = 0; i < numPoints - 1; ++i)
            {
                if (xPositions[i + 1] >= x)
                {
                    index = i;
                    break;
                }
            }
 
            const float x0 = xPositions[index];
            const float x1 = xPositions[index + 1];
            const float y0 = yPositions[index];
            const float y1 = yPositions[index + 1];
 
            const float dx = x1 - x0;
            const float dy = y1 - y0;
            const float h0 = (index > 0) ? x0 - xPositions[index - 1] : dx;
            const float h1 = (index < numPoints - 2) ? xPositions[index + 2] - x1 : dx;
 
            const float t = (x - x0) / dx;
            const float t2 = t * t;
 
            const float p = (1.0f - t) * y0 + t * y1;
 
            if (dy * h0 > 0.0f && dy * h1 > 0.0f)
            {
                const float a = (2.0f * dy - h0 * (y1 - y0) / dx) / h0;
                const float b = (2.0f * dy - h1 * (y1 - y0) / dx) / h1;
 
                interpolatedValues.push_back(t * (1.0f - t) * (a * h0 + b * h1) + p);
            }
            else
            {
                interpolatedValues.push_back(p);
            }
        }
 
        return interpolatedValues;
    }
 
    Eigen::VectorXf ShapePreservingCubicInterpolator::interpolateVectorEigen(const Eigen::VectorXf& xPositions,
                                                                             const Eigen::VectorXf& yPositions,
                                                                             const Eigen::VectorXf& xValues)
    {
        jassert(xPositions.size() == yPositions.size());
 
        const int numPoints = static_cast<int>(xPositions.size());
 
        if (numPoints < 2)
        {
            jassert(false); // Not enough points for interpolation
            return Eigen::VectorXf::Zero(xValues.size());
        }
 
        Eigen::VectorXf interpolatedValues(xValues.size());
 
        for (int k = 0; k < xValues.size(); ++k)
        {
            float x = xValues[k];
            int index = 0;
 
            // Find the index of the nearest point to the left of x
            for (int i = 0; i < numPoints - 1; ++i)
            {
                if (xPositions[i + 1] >= x)
                {
                    index = i;
                    break;
                }
            }
 
            const float x0 = xPositions[index];
            const float x1 = xPositions[index + 1];
            const float y0 = yPositions[index];
            const float y1 = yPositions[index + 1];
 
            const float dx = x1 - x0;
            const float dy = y1 - y0;
            const float h0 = (index > 0) ? x0 - xPositions[index - 1] : dx;
            const float h1 = (index < numPoints - 2) ? xPositions[index + 2] - x1 : dx;
 
            const float t = (x - x0) / dx;
            const float t2 = t * t;
 
            const float p = (1.0f - t) * y0 + t * y1;
 
            if (dy * h0 > 0.0f && dy * h1 > 0.0f)
            {
                const float a = (2.0f * dy - h0 * (y1 - y0) / dx) / h0;
                const float b = (2.0f * dy - h1 * (y1 - y0) / dx) / h1;
 
                interpolatedValues[k] = t * (1.0f - t) * (a * h0 + b * h1) + p;
            }
            else
            {
                interpolatedValues[k] = p;
            }
        }
 
        return interpolatedValues;
    }
 
    //============================================================================================================================================
 
    std::vector<float> LinearInterpolator::interpolate(const std::vector<float>& xValues,
                                                       const std::vector<float>& yValues,
                                                       const std::vector<float>& vals)
    {
        // check if input vectors are of equal size
        if (xValues.size() != yValues.size())
        {
            // handle error: Vectors must be of equal size
            return std::vector<float>(); // return an empty vector to indicate the
            // error
        }
 
        std::vector<float> interpolatedValues;
        int size = static_cast<int>(vals.size());
 
        for (int i = 0; i < size; ++i)
        {
            float valToInterpolateFor = vals[i];
 
            // find the indices of the two nearest x-values in the input vector
            int index1 = 0;
            int index2 = 0;
            int xSize = static_cast<int>(xValues.size());
            for (int j = 0; j < xSize; ++j)
            {
                if (xValues[j] <= valToInterpolateFor)
                    index1 = j;
                else
                {
                    index2 = j;
                    break;
                }
            }
 
            // perform linear interpolation
            float x1 = xValues[index1];
            float x2 = xValues[index2];
            float y1 = yValues[index1];
            float y2 = yValues[index2];
 
            // avoid division by zero
            if (std::isnan(valToInterpolateFor))
            {
                interpolatedValues.push_back(std::numeric_limits<float>::quiet_NaN());
            }
            else if (std::isnan(x1) || std::isnan(x2) || std::isnan(y1) || std::isnan(y2))
            {
                // handle NaN values in xValues or yValues
                interpolatedValues.push_back(std::numeric_limits<float>::quiet_NaN());
            }
            else
            {
                if (x2 == x1)
                {
                    // handle error: Division by zero
                    interpolatedValues.push_back(0.0f); // or push back an appropriate error value
                }
                else
                {
                    float interpolatedValue = y1 + ((valToInterpolateFor - x1) * (y2 - y1)) / (x2 - x1);
                    interpolatedValues.push_back(interpolatedValue);
                }
            }
        }
 
        return interpolatedValues;
    }
 
    Eigen::VectorXf LinearInterpolator::interpolateEigen(const Eigen::VectorXf& xValues, const Eigen::VectorXf& yValues,
                                                         const Eigen::VectorXf& vals)
    {
        // check if input vectors are of equal size
        if (xValues.size() != yValues.size())
        {
            jassertfalse; // Not enough points for interpolation
            return Eigen::VectorXf::Zero(xValues.size());
        }
 
        Eigen::VectorXf interpolatedValues(vals.size());
 
        for (int i = 0; i < vals.size(); ++i)
        {
            float valToInterpolateFor = vals[i];
 
            // find the indices of the two nearest x-values in the input vector
            int index1 = 0;
            int index2 = 0;
            int xSize = static_cast<int>(xValues.size());
            for (int j = 0; j < xSize; ++j)
            {
                if (xValues[j] <= valToInterpolateFor)
                    index1 = j;
                else
                {
                    index2 = j;
                    break;
                }
            }
 
            // perform linear interpolation
            float x1 = xValues[index1];
            float x2 = xValues[index2];
            float y1 = yValues[index1];
            float y2 = yValues[index2];
 
            if (valToInterpolateFor < x1 || valToInterpolateFor > x2)
            {
                // Extrapolation, return NaN
                interpolatedValues[i] = std::numeric_limits<float>::quiet_NaN();
            }
            else
            {
                if (x2 == x1)
                {
                    // handle error: Division by zero
                    interpolatedValues.conservativeResize(interpolatedValues.size() + 1);
                    interpolatedValues(interpolatedValues.size() - 1) =
                            std::numeric_limits<float>::quiet_NaN(); // or push back an appropriate error value
                }
                else
                {
                    float interpolatedValue = y1 + ((valToInterpolateFor - x1) * (y2 - y1)) / (x2 - x1);
                    interpolatedValues[i] = interpolatedValue;
                }
            }
        }
 
        return interpolatedValues;
    }
 
    Eigen::MatrixXf LinearInterpolator::interpolate2DEigen(const Eigen::VectorXf& xValues,
                                                           const Eigen::MatrixXf& yValues, const Eigen::VectorXf& vals)
    {
        // Check if input dimensions are valid
        if (xValues.size() != yValues.rows() || yValues.rows() == 0)
        {
            // Invalid input dimensions
            jassertfalse;
        }
 
        Eigen::MatrixXf resultMatrix(vals.size(), yValues.cols());
 
        for (int col = 0; col < yValues.cols(); ++col)
        {
            // Interpolate using Shape Preserving Cubic Interpolation
            Eigen::VectorXf interpolatedValues = LinearInterpolator::interpolateEigen(xValues, yValues.col(col), vals);
 
            resultMatrix.col(col) = interpolatedValues;
        }
 
        return resultMatrix;
    }
 
    std::vector<std::vector<float>> LinearInterpolator::interpolate2D(const std::vector<float>& xValues,
                                                                      const std::vector<std::vector<float>>& yValues,
                                                                      const std::vector<float>& vals)
    {
        // Check if input dimensions are valid
        if (xValues.size() != yValues[0].size() || yValues[0].size() == 0)
        {
            // Invalid input dimensions
            jassertfalse;
        }
 
        // Convert values in vals to NaN if they are outside the range of xValues
        std::vector<float> sanitizedVals = vals;
        std::transform(sanitizedVals.begin(), sanitizedVals.end(), sanitizedVals.begin(), [&](float val) {
            return (val < xValues.front() || val > xValues.back()) ? std::numeric_limits<float>::quiet_NaN() : val;
        });
 
        std::vector<std::vector<float>> resultMatrix;
 
        for (int row = 0; row < yValues.size(); ++row)
        {
            // Interpolate using Linear Interpolation
            std::vector<float> interpolatedValues =
                    LinearInterpolator::interpolate(xValues, yValues[row], sanitizedVals);
 
            resultMatrix.push_back(interpolatedValues);
        }
 
        return resultMatrix;
    }
 
} // namespace krotos