qwt_grid_data.hpp

/******************************************************************************
 * Qwt Widget Library
 * Copyright (C) 2024   ChenZongYan <czy.t@163.com>
 *****************************************************************************/
#ifndef QWT_GRID_DATA_HPP
#define QWT_GRID_DATA_HPP
 
#include <vector>
#include <algorithm>
#include <limits>
#include <cmath>
#include <cassert>
#include <stdexcept>
 
#include "qwt_math.h"
 
template< typename T,
          typename XContainer    = std::vector< T >,
          typename YContainer    = std::vector< T >,
          typename DataColumn    = std::vector< T >,
          typename DataContainer = std::vector< DataColumn > >
class QwtGridData
{
public:
    using value_type          = T;              
    using x_container_type    = XContainer;     
    using y_container_type    = YContainer;     
    using data_column_type    = DataColumn;     
    using data_container_type = DataContainer;  
    using size_type           = typename DataColumn::size_type;
    enum ResampleMode
    {
        NearestNeighbour,
        BilinearInterpolation,
        BicubicInterpolation
    };
 
    QwtGridData()
        : m_mode(NearestNeighbour), m_xMin(0.0), m_xMax(0.0), m_yMin(0.0), m_yMax(0.0), m_dataMax(0.0), m_dataMin(0.0)
    {
    }
 
    QwtGridData(const x_container_type& xAxis,
                const y_container_type& yAxis,
                const data_container_type& data,
                ResampleMode mode = NearestNeighbour)
        : m_xAxis(xAxis), m_yAxis(yAxis), m_data(data), m_mode(mode)
    {
        validate();
        findValueRange();
        m_xMin = xAxis.front();
        m_xMax = xAxis.back();
        m_yMin = yAxis.front();
        m_yMax = yAxis.back();
    }
 
    void setValue(const x_container_type& xAxis, const y_container_type& yAxis, const data_container_type& data)
    {
        m_xAxis = xAxis;
        m_yAxis = yAxis;
        m_data  = data;
        validate();
        findValueRange();
        m_xMin = xAxis.front();
        m_xMax = xAxis.back();
        m_yMin = yAxis.front();
        m_yMax = yAxis.back();
    }
 
    T operator()(T x, T y) const
    {
        return value(x, y);
    }
 
    T operator[](const std::pair< T, T >& xy) const
    {
        return value(xy.first, xy.second);
    }
    T value(T x, T y) const
    {
        switch (m_mode) {
        case NearestNeighbour:
            return nearestNeighbour(x, y);
        case BilinearInterpolation:
            return bilinearInterpolation(x, y);
        case BicubicInterpolation:
            return bicubicInterpolation(x, y);
        default:
            throw std::runtime_error("Unknown resampling mode.");
        }
    }
 
    void setResampleMode(ResampleMode mode)
    {
        m_mode = mode;
    }
 
    ResampleMode resampleMode() const
    {
        return m_mode;
    }
 
    size_type xSize() const
    {
        return m_xAxis.size();
    }
 
    size_type ySize() const
    {
        return m_yAxis.size();
    }
 
    std::pair< size_type, size_type > valueSize() const
    {
        return std::make_pair(xSize(), ySize());
    }
 
    T atX(size_type ix) const
    {
        return m_xAxis.at(ix);
    }
 
    T atY(size_type iy) const
    {
        return m_yAxis.at(iy);
    }
 
    T atValue(size_type ix, size_type iy) const
    {
        return m_data.at(ix).at(iy);
    }
 
    const x_container_type& xAxis() const
    {
        return m_xAxis;
    }
 
    const y_container_type& yAxis() const
    {
        return m_yAxis;
    }
 
    const data_container_type& data() const
    {
        return m_data;
    }
 
    bool valid() const
    {
        if (m_xAxis.empty() || m_yAxis.empty() || m_data.empty()) {
            return false;
        }
        if (m_xAxis.size() != m_data.size()) {
            return false;
        }
        for (const auto& column : m_data) {
            if (column.size() != m_yAxis.size()) {
                return false;
            }
        }
        if (!std::is_sorted(m_xAxis.begin(), m_xAxis.end())) {
            return false;
        }
        if (!std::is_sorted(m_yAxis.begin(), m_yAxis.end())) {
            return false;
        }
        return true;
    }
 
    void validate()
    {
        if (m_xAxis.empty() || m_yAxis.empty() || m_data.empty()) {
            throw std::invalid_argument("Axes or data cannot be empty.");
        }
        if (m_data.size() != m_xAxis.size()) {
            throw std::invalid_argument("Number of columns in data must match x-axis size.");
        }
        for (const auto& column : m_data) {
            if (column.size() != m_yAxis.size()) {
                throw std::invalid_argument("Number of rows in data must match y-axis size.");
            }
        }
        if (!std::is_sorted(m_xAxis.begin(), m_xAxis.end())) {
            std::sort(m_xAxis.begin(), m_xAxis.end());
        }
        if (!std::is_sorted(m_yAxis.begin(), m_yAxis.end())) {
            std::sort(m_yAxis.begin(), m_yAxis.end());
        }
    }
 
    T xMin() const
    {
        return m_xMin;
    }
    T xMax() const
    {
        return m_xMax;
    }
    T yMin() const
    {
        return m_yMin;
    }
    T yMax() const
    {
        return m_yMax;
    }
    T dataMin() const
    {
        return m_dataMin;
    }
    T dataMax() const
    {
        return m_dataMax;
    }
 
public:
    // static pulic function
    template< typename Container >
    static size_type findClosestIndex(const Container& arr, T val)
    {
        auto it = std::lower_bound(arr.begin(), arr.end(), val);
        if (it == arr.begin())
            return 0;
        if (it == arr.end())
            return arr.size() - 1;
        size_type idx = std::distance(arr.begin(), it);
        return (std::abs(arr[ idx ] - val) < std::abs(arr[ idx - 1 ] - val)) ? idx : idx - 1;
    }
 
    template< typename Container >
    static size_type findLowerIndex(const Container& arr, T val)
    {
        auto it = std::lower_bound(arr.begin(), arr.end(), val);
        if (it == arr.begin())
            return 0;
        return std::distance(arr.begin(), it) - 1;
    }
 
    template< typename V >
    static const V& clamp(const V& value, const V& lo, const V& hi)
    {
        return (value < lo) ? lo : (hi < value) ? hi : value;
    }
 
protected:
    void findValueRange()
    {
        m_dataMin = std::numeric_limits< T >::max();
        m_dataMax = std::numeric_limits< T >::lowest();
        bool found = false;
        for (const auto& column : m_data) {
            for (const auto& val : column) {
                if (qwt_is_nan_or_inf(val))
                    continue;
                m_dataMin = std::min(m_dataMin, val);
                m_dataMax = std::max(m_dataMax, val);
                found = true;
            }
        }
        if (!found) {
            m_dataMin = std::numeric_limits< T >::quiet_NaN();
            m_dataMax = std::numeric_limits< T >::quiet_NaN();
        }
    }
 
    T nearestNeighbour(T x, T y) const
    {
        size_type xIdx = findClosestIndex(m_xAxis, x);
        size_type yIdx = findClosestIndex(m_yAxis, y);
        return m_data[ xIdx ][ yIdx ];
    }
 
    T bilinearInterpolation(T x, T y) const
    {
        size_type x0Idx = findLowerIndex(m_xAxis, x);
        size_type x1Idx = x0Idx + 1;
        size_type y0Idx = findLowerIndex(m_yAxis, y);
        size_type y1Idx = y0Idx + 1;
 
        T x0 = m_xAxis[ x0Idx ], x1 = m_xAxis[ x1Idx ];
        T y0 = m_yAxis[ y0Idx ], y1 = m_yAxis[ y1Idx ];
 
        T f00 = m_data[ x0Idx ][ y0Idx ];
        T f10 = m_data[ x1Idx ][ y0Idx ];
        T f01 = m_data[ x0Idx ][ y1Idx ];
        T f11 = m_data[ x1Idx ][ y1Idx ];
 
        T dx = (x - x0) / (x1 - x0);
        T dy = (y - y0) / (y1 - y0);
 
        return (1 - dx) * (1 - dy) * f00 + dx * (1 - dy) * f10 + (1 - dx) * dy * f01 + dx * dy * f11;
    }
 
    T bicubicInterpolation(T x, T y) const
    {
        // Find surrounding x indices
        auto xIt = std::lower_bound(m_xAxis.begin(), m_xAxis.end(), x);
        size_type x0, x1, x2, x3;
        T xWeight;
 
        if (xIt == m_xAxis.begin()) {
            x0      = 0;
            x1      = 0;
            x2      = 1;
            x3      = 2;
            xWeight = 0.0;
        } else if (xIt == m_xAxis.end()) {
            x0      = m_xAxis.size() - 3;
            x1      = m_xAxis.size() - 2;
            x2      = m_xAxis.size() - 1;
            x3      = m_xAxis.size() - 1;
            xWeight = 1.0;
        } else {
            x0 = xIt - 2 - m_xAxis.begin();
            if (x0 < 0)
                x0 = 0;
            x1 = xIt - 1 - m_xAxis.begin();
            x2 = xIt - m_xAxis.begin();
            x3 = xIt + 1 - m_xAxis.begin();
            if (x3 >= m_xAxis.size())
                x3 = m_xAxis.size() - 1;
            xWeight = static_cast< T >(x - m_xAxis[ x1 ]) / (m_xAxis[ x2 ] - m_xAxis[ x1 ]);
        }
 
        // Find surrounding y indices
        auto yIt = std::lower_bound(m_yAxis.begin(), m_yAxis.end(), y);
        size_type y0, y1, y2, y3;
        T yWeight;
 
        if (yIt == m_yAxis.begin()) {
            y0      = 0;
            y1      = 0;
            y2      = 1;
            y3      = 2;
            yWeight = 0.0;
        } else if (yIt == m_yAxis.end()) {
            y0      = m_yAxis.size() - 3;
            y1      = m_yAxis.size() - 2;
            y2      = m_yAxis.size() - 1;
            y3      = m_yAxis.size() - 1;
            yWeight = 1.0;
        } else {
            y0 = yIt - 2 - m_yAxis.begin();
            if (y0 < 0)
                y0 = 0;
            y1 = yIt - 1 - m_yAxis.begin();
            y2 = yIt - m_yAxis.begin();
            y3 = yIt + 1 - m_yAxis.begin();
            if (y3 >= m_yAxis.size())
                y3 = m_yAxis.size() - 1;
            yWeight = static_cast< T >(y - m_yAxis[ y1 ]) / (m_yAxis[ y2 ] - m_yAxis[ y1 ]);
        }
 
        // Hermite basis functions
        auto h00 = [](T t) { return (1 + 2 * t) * (1 - t) * (1 - t); };
        auto h10 = [](T t) { return t * (1 - t) * (1 - t); };
        auto h01 = [](T t) { return t * t * (3 - 2 * t); };
        auto h11 = [](T t) { return t * t * (t - 1); };
 
        // Interpolate in x direction for each y position
        T values[ 4 ][ 4 ];
        for (int i = 0; i < 4; ++i) {
            // Get the four y values for current x positions
            T v[ 4 ];
            v[ 0 ] = m_data[ x0 ][ y0 + i ];
            v[ 1 ] = m_data[ x1 ][ y0 + i ];
            v[ 2 ] = m_data[ x2 ][ y0 + i ];
            v[ 3 ] = m_data[ x3 ][ y0 + i ];
 
            // Interpolate in x direction
            values[ 0 ][ i ] = h00(xWeight) * v[ 1 ]
                               + h10(xWeight) * (m_xAxis[ x2 ] - m_xAxis[ x1 ])
                                     * ((v[ 2 ] - v[ 1 ]) / (m_xAxis[ x2 ] - m_xAxis[ x1 ])
                                        + (v[ 2 ] - v[ 1 ] - (v[ 1 ] - v[ 0 ]) / (m_xAxis[ x1 ] - m_xAxis[ x0 ]))
                                              / (m_xAxis[ x2 ] - m_xAxis[ x0 ]) * (m_xAxis[ x2 ] - m_xAxis[ x1 ]))
                               + h01(xWeight) * v[ 2 ]
                               + h11(xWeight) * (m_xAxis[ x2 ] - m_xAxis[ x1 ])
                                     * ((v[ 2 ] - v[ 1 ]) / (m_xAxis[ x2 ] - m_xAxis[ x1 ])
                                        + (v[ 3 ] - v[ 2 ] - (v[ 2 ] - v[ 1 ]) / (m_xAxis[ x2 ] - m_xAxis[ x1 ]))
                                              / (m_xAxis[ x3 ] - m_xAxis[ x1 ]) * (m_xAxis[ x2 ] - m_xAxis[ x1 ]));
        }
 
        // Interpolate in y direction
        T v[ 4 ];
        v[ 0 ] = values[ 0 ][ 0 ];
        v[ 1 ] = values[ 0 ][ 1 ];
        v[ 2 ] = values[ 0 ][ 2 ];
        v[ 3 ] = values[ 0 ][ 3 ];
 
        return h00(yWeight) * v[ 1 ]
               + h10(yWeight) * (m_yAxis[ y2 ] - m_yAxis[ y1 ])
                     * ((v[ 2 ] - v[ 1 ]) / (m_yAxis[ y2 ] - m_yAxis[ y1 ])
                        + (v[ 2 ] - v[ 1 ] - (v[ 1 ] - v[ 0 ]) / (m_yAxis[ y1 ] - m_yAxis[ y0 ]))
                              / (m_yAxis[ y2 ] - m_yAxis[ y0 ]) * (m_yAxis[ y2 ] - m_yAxis[ y1 ]))
               + h01(yWeight) * v[ 2 ]
               + h11(yWeight) * (m_yAxis[ y2 ] - m_yAxis[ y1 ])
                     * ((v[ 2 ] - v[ 1 ]) / (m_yAxis[ y2 ] - m_yAxis[ y1 ])
                        + (v[ 3 ] - v[ 2 ] - (v[ 2 ] - v[ 1 ]) / (m_yAxis[ y2 ] - m_yAxis[ y1 ]))
                              / (m_yAxis[ y3 ] - m_yAxis[ y1 ]) * (m_yAxis[ y2 ] - m_yAxis[ y1 ]));
    }
 
private:
    x_container_type m_xAxis;  
    y_container_type m_yAxis;  
    data_container_type m_data;
    ResampleMode m_mode;                                     
    T m_xMin, m_xMax, m_yMin, m_yMax, m_dataMin, m_dataMax;  
};
 
#endif  // QWT_GRID_DATA_HPP