filtering.hpp

Go to the documentation of this file.

/*****************************************************************************
    NumeRe: Framework fuer Numerische Rechnungen
    Copyright (C) 2020  Erik Haenel et al.
 
    This program is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.
 
    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.
 
    You should have received a copy of the GNU General Public License
    along with this program.  If not, see <http://www.gnu.org/licenses/>.
******************************************************************************/
 
#ifndef FILTERING_HPP
#define FILTERING_HPP
 
#include <utility>
#include <vector>
#include <cmath>
#include <queue>
#include "../io/file.hpp"
#include "../utils/tools.hpp"
#include "../ParserLib/muParserDef.h"
 
//void showMatrix(const vector<vector<double> >&);
 
namespace NumeRe
{
    struct FilterSettings
    {
        enum FilterType
        {
            FILTER_NONE,
            FILTER_WEIGHTED_LINEAR,
            FILTER_GAUSSIAN,
            FILTER_SAVITZKY_GOLAY
        };
 
        FilterType type;
        size_t row;
        size_t col;
        double alpha;
 
        FilterSettings(FilterType _type = FILTER_NONE, size_t _row = 1u, size_t _col = 1u, double _alpha = NAN) : type(_type), row(std::max(_row, 1u)), col(std::max(_col, 1u)), alpha(_alpha)
        {
            //
        }
    };
 
 
    inline double pow2(double val)
    {
        return val * val;
    }
 
 
    inline double pow3(double val)
    {
        return val * val * val;
    }
 
 
    inline double pow4(double val)
    {
        return val * val * val * val;
    }
 
 
    using FilterBuffer = std::queue<mu::value_type>;
    using FilterBuffer2D = std::queue<std::vector<mu::value_type>>;
 
 
    class Filter
    {
        protected:
            FilterSettings::FilterType m_type;
            std::pair<size_t, size_t> m_windowSize;
            bool m_isConvolution;
            FilterBuffer m_buffer;
            FilterBuffer2D m_buffer2D;
 
        public:
            Filter(size_t row, size_t col) : m_type(FilterSettings::FILTER_NONE), m_isConvolution(false)
            {
                // This expression avoids that someone tries to
                // create a filter with a zero dimension.
                m_windowSize = std::make_pair(std::max(row, 1u), std::max(col, 1u));
            }
 
            virtual ~Filter() {}
 
            virtual double operator()(size_t i, size_t j) const = 0;
 
            virtual mu::value_type apply(size_t i, size_t j, const mu::value_type& val) const = 0;
 
            bool isConvolution() const
            {
                return m_isConvolution;
            }
 
            FilterSettings::FilterType getType() const
            {
                return m_type;
            }
 
            std::pair<size_t,size_t> getWindowSize() const
            {
                return m_windowSize;
            }
 
            FilterBuffer& getBuffer()
            {
                return m_buffer;
            }
 
            FilterBuffer2D& get2DBuffer()
            {
                return m_buffer2D;
            }
    };
 
 
    class WeightedLinearFilter : public Filter
    {
        private:
            bool is2D;
            std::vector<std::vector<double> > m_filterKernel;
 
            void createKernel()
            {
                if (!is2D)
                    is2D = m_windowSize.first > 1 && m_windowSize.second > 1;
 
                if (!(m_windowSize.first % 2))
                    m_windowSize.first++;
 
                m_filterKernel = std::vector<std::vector<double> >(m_windowSize.first, std::vector<double>(m_windowSize.second, 0.0));
 
                double mean_row = m_windowSize.first > 1 ? 0.5 : 0.0;
                double mean_col = m_windowSize.second > 1 ? 0.5 : 0.0;
                double sum = 0;
 
                // Calculate the filter values
                for (size_t i = 0; i < m_windowSize.first; i++)
                {
                    for (size_t j = 0; j < m_windowSize.second; j++)
                    {
                        if (sqrt(pow2(i/((double)std::max(1u, m_windowSize.first-1))-mean_row) + pow2(j/((double)std::max(1u, m_windowSize.second-1))-mean_col)) <= 0.5)
                        {
                            m_filterKernel[i][j] = fabs(sqrt(pow2(i/((double)std::max(1u, m_windowSize.first-1))-mean_row) + pow2(j/((double)std::max(1u, m_windowSize.second-1))-mean_col)) - 0.5);
                            sum += m_filterKernel[i][j];
                        }
                    }
                }
 
                for (size_t i = 0; i < m_windowSize.first; i++)
                {
                    for (size_t j = 0; j < m_windowSize.second; j++)
                    {
                        m_filterKernel[i][j] /= sum;
                    }
                }
            }
 
        public:
            WeightedLinearFilter(size_t row, size_t col, bool force2D = false) : Filter(row, col)
            {
                m_type = FilterSettings::FILTER_WEIGHTED_LINEAR;
                m_isConvolution = true;
                is2D = force2D;
 
                createKernel();
            }
 
            virtual ~WeightedLinearFilter() override
            {
                m_filterKernel.clear();
            }
 
            virtual double operator()(size_t i, size_t j) const override
            {
                if (i >= m_windowSize.first || j >= m_windowSize.second)
                    return NAN;
 
                return m_filterKernel[i][j];
            }
 
            virtual mu::value_type apply(size_t i, size_t j, const mu::value_type& val) const override
            {
                if (i >= m_windowSize.first || j >= m_windowSize.second)
                    return NAN;
 
                return m_filterKernel[i][j]*val;
            }
    };
 
 
    class GaussianFilter : public Filter
    {
        private:
            std::vector<std::vector<double> > m_filterKernel;
 
            void createKernel(double sigma)
            {
                m_filterKernel = std::vector<std::vector<double> >(m_windowSize.first, std::vector<double>(m_windowSize.second, 0.0));
 
                double sum = 0.0;
                double mean_row = (m_windowSize.first - 1) * 0.5;
                double mean_col = (m_windowSize.second - 1) * 0.5;
 
                for (size_t i = 0; i < m_windowSize.first; i++)
                {
                    for (size_t j = 0; j < m_windowSize.second; j++)
                    {
                        m_filterKernel[i][j] = exp(-(pow2(i-mean_row) + pow2(j-mean_col)) / (2*pow2(sigma))) / (2*M_PI*pow2(sigma));
                        sum += m_filterKernel[i][j];
                    }
                }
 
                for (size_t i = 0; i < m_windowSize.first; i++)
                {
                    for (size_t j = 0; j < m_windowSize.second; j++)
                    {
                        m_filterKernel[i][j] /= sum;
                    }
                }
            }
 
        public:
            GaussianFilter(size_t row, size_t col, double sigma) : Filter(row, col)
            {
                m_type = FilterSettings::FILTER_GAUSSIAN;
                m_isConvolution = true;
 
                createKernel(sigma);
            }
 
            virtual ~GaussianFilter() override
            {
                m_filterKernel.clear();
            }
 
            virtual double operator()(size_t i, size_t j) const override
            {
                if (i >= m_windowSize.first || j >= m_windowSize.second)
                    return NAN;
 
                return m_filterKernel[i][j];
            }
 
            virtual mu::value_type apply(size_t i, size_t j, const mu::value_type& val) const override
            {
                if (i >= m_windowSize.first || j >= m_windowSize.second)
                    return NAN;
 
                if (mu::isnan(val))
                    return 0.0;
 
                // Gaussian is symmetric, therefore the convolution does not
                // need to be inverted
                return m_filterKernel[i][j] * val;
            }
    };
 
 
    class SavitzkyGolayFilter : public Filter
    {
        private:
            std::vector<std::vector<double>> m_filterKernel;
 
            long long int findColumn(const FileView& _view)
            {
                std::vector<size_t> vWindowSizes;
 
                // Decode all contained window sizes
                for (long long int j = 0; j < _view.getCols(); j++)
                {
                    vWindowSizes.push_back(StrToInt(_view.getColumnHead(j).substr(0, _view.getColumnHead(j).find('x'))));
                }
 
                // Window size is already smaller than the first
                // available window size
                if (m_windowSize.first < vWindowSizes.front())
                {
                    m_windowSize.first = vWindowSizes.front();
                    m_windowSize.second = m_windowSize.first;
                    m_filterKernel = std::vector<std::vector<double> >(m_windowSize.first, std::vector<double>(m_windowSize.second, 0.0));
 
                    return 0;
                }
 
                for (long long int j = 0; j < vWindowSizes.size(); j++)
                {
                    // Found a perfect match?
                    if (m_windowSize.first == vWindowSizes[j])
                    {
                        m_windowSize.second = m_windowSize.first;
                        m_filterKernel = std::vector<std::vector<double> >(m_windowSize.first, std::vector<double>(m_windowSize.second, 0.0));
 
                        return j;
                    }
 
                    // Is there a nearest match? (We assume
                    // ascending order of the window sizes)
                    if (m_windowSize.first > vWindowSizes[j] && j+1 < vWindowSizes.size() && m_windowSize.first < vWindowSizes[j+1])
                    {
                        // If the following fits better, increment the index
                        if (m_windowSize.first - vWindowSizes[j] > vWindowSizes[j+1] - m_windowSize.first)
                            j++;
 
                        m_windowSize.first = vWindowSizes[j];
                        m_windowSize.second = m_windowSize.first;
                        m_filterKernel = std::vector<std::vector<double> >(m_windowSize.first, std::vector<double>(m_windowSize.second, 0.0));
 
                        return j;
                    }
                    else if (m_windowSize.first > vWindowSizes[j] && j+1 == vWindowSizes.size())
                    {
                        m_windowSize.first = vWindowSizes[j];
                        m_windowSize.second = m_windowSize.first;
                        m_filterKernel = std::vector<std::vector<double> >(m_windowSize.first, std::vector<double>(m_windowSize.second, 0.0));
 
                        return j;
                    }
                }
 
                return 0;
            }
 
            void createKernel()
            {
                // Ensure odd numbers
                if (!(m_windowSize.first % 2))
                    m_windowSize.first++;
 
                // Create a 2D kernel
                if (m_windowSize.second > 1)
                {
                    // Create a file instance
                    GenericFile* _file = getFileByType("<>/params/savitzky_golay_coeffs_2D.dat");
 
                    if (!_file)
                        return;
 
                    // Read the contents and assign it to a view
                    try
                    {
                        _file->read();
                    }
                    catch (...)
                    {
                        delete _file;
                        throw;
                    }
 
                    FileView _view(_file);
 
                    // Find the possible window sizes
                    long long int j = findColumn(_view);
 
                    // First element in the column is the
                    // central element in the matrix
                    m_filterKernel[m_windowSize.first/2][m_windowSize.first/2] = _view.getElement(m_windowSize.first*m_windowSize.first/2, j).real();
 
                    for (long long int i = 0; i < m_windowSize.first*m_windowSize.first/2; i++)
                    {
                        // left part
                        m_filterKernel[m_windowSize.first - 1 - i % m_windowSize.first][i / m_windowSize.first] = _view.getElement(i, j).real();
                        // right part
                        m_filterKernel[i % m_windowSize.first][m_windowSize.first - i / m_windowSize.first - 1] = _view.getElement(i, j).real();
                        // middle column
                    }
 
                    delete _file;
                    return;
                }
 
                // Create a 1D kernel
                m_filterKernel = std::vector<std::vector<double> >(m_windowSize.first, std::vector<double>(m_windowSize.second, 0.0));
 
                for (size_t i = 0; i < m_windowSize.first; i++)
                {
                    m_filterKernel[i][0] = (3.0*pow2(m_windowSize.first) - 7.0 - 20.0*pow2((int)i - (int)m_windowSize.first/2)) / 4.0 / (m_windowSize.first * (pow2(m_windowSize.first) - 4.0) / 3.0);
                }
            }
 
        public:
            SavitzkyGolayFilter(size_t row, size_t col) : Filter(row, col)
            {
                m_type = FilterSettings::FILTER_SAVITZKY_GOLAY;
                m_isConvolution = true;
 
                createKernel();
            }
 
            virtual ~SavitzkyGolayFilter() override
            {
                m_filterKernel.clear();
            }
 
            virtual double operator()(size_t i, size_t j) const override
            {
                if (i >= m_windowSize.first || j >= m_windowSize.second)
                    return NAN;
 
                return m_filterKernel[i][j];
            }
 
            virtual mu::value_type apply(size_t i, size_t j, const mu::value_type& val) const override
            {
                if (i >= m_windowSize.first || j >= m_windowSize.second)
                    return NAN;
 
                if (mu::isnan(val))
                    return 0.0;
 
                return m_filterKernel[i][j] * val;
            }
    };
 
 
    class SavitzkyGolayDiffFilter : public Filter
    {
        private:
            std::vector<double> m_filterKernel;
 
            void createKernel(size_t nthDerivative)
            {
                // Ensure odd numbers
                if (!(m_windowSize.first % 2))
                    m_windowSize.first++;
 
                // Create a 1D kernel
                m_filterKernel = std::vector<double>(m_windowSize.first, 0.0);
                double m = m_windowSize.first;
 
                for (size_t i = 0; i < m_windowSize.first; i++)
                {
                    double I = (int)i - floor(m/2.0);
 
                    if (nthDerivative == 1)
                        m_filterKernel[i] = 15.0 * (5.0*(3.0*pow4(m) - 18.0*pow2(m) + 31.0)*I - 28.0*(3.0*pow2(m)-7.0)*pow3(I))
                                             / (m * (pow2(m) - 1.0)*(3.0*pow4(m) - 39.0*pow2(m) + 108.0));
                    else if (nthDerivative == 2)
                        m_filterKernel[i] = 30.0 * (12.0*m*pow2(I) - m * (pow2(m) - 1.0))
                                             / (pow2(m) * (pow2(m) - 1.0) * (pow2(m) - 4.0));
                    else if (nthDerivative == 3)
                        m_filterKernel[i] = 6 * 420.0 * (-(3.0*pow2(m) - 7.0)*I + 20.0*pow3(I))
                                             / (m * (pow2(m) - 1.0)*(3.0*pow4(m) - 39.0*pow2(m) + 108.0));
                    else
                        m_filterKernel[i] = 15.0 * (5.0*(3.0*pow4(m) - 18.0*pow2(m) + 31.0)*I - 28.0*(3.0*pow2(m)-7.0)*pow3(I))
                                             / (m * (pow2(m) - 1.0)*(3.0*pow4(m) - 39.0*pow2(m) + 108.0));
                }
            }
 
        public:
            SavitzkyGolayDiffFilter(size_t row, size_t nthDerivative) : Filter(row, 1)
            {
                m_type = FilterSettings::FILTER_SAVITZKY_GOLAY;
                m_isConvolution = true;
 
                createKernel(nthDerivative);
            }
 
            virtual ~SavitzkyGolayDiffFilter() override
            {
                m_filterKernel.clear();
            }
 
            virtual double operator()(size_t i, size_t j) const override
            {
                if (i >= m_windowSize.first)
                    return NAN;
 
                return m_filterKernel[i];
            }
 
            virtual mu::value_type apply(size_t i, size_t j, const mu::value_type& val) const override
            {
                if (i >= m_windowSize.first)
                    return NAN;
 
                if (mu::isnan(val))
                    return 0.0;
 
                return m_filterKernel[i] * val;
            }
    };
 
 
    inline Filter* createFilter(const FilterSettings& _settings)
    {
        switch (_settings.type)
        {
            case FilterSettings::FILTER_NONE:
                return nullptr;
            case FilterSettings::FILTER_WEIGHTED_LINEAR:
                return new WeightedLinearFilter(_settings.row, _settings.col);
            case FilterSettings::FILTER_GAUSSIAN:
                return new GaussianFilter(_settings.row, _settings.col, (std::max(_settings.row, _settings.col)-1)/(2*_settings.alpha));
            case FilterSettings::FILTER_SAVITZKY_GOLAY:
                return new SavitzkyGolayFilter(_settings.row, _settings.col);
        }
 
        return nullptr;
    }
 
 
    class RetouchRegion : public Filter
    {
        private:
            std::vector<mu::value_type> m_left;
            std::vector<mu::value_type> m_right;
            std::vector<mu::value_type> m_top;
            std::vector<mu::value_type> m_bottom;
            bool is2D;
            std::vector<std::vector<double> > m_filterKernel;
            mu::value_type m_fallback;
            bool m_invertedKernel;
 
            void createKernel()
            {
                m_filterKernel = std::vector<std::vector<double> >(m_windowSize.first, std::vector<double>(m_windowSize.second, 0.0));
 
                if (!is2D)
                    is2D = m_windowSize.first > 1 && m_windowSize.second > 1;
 
                double mean_row = m_windowSize.first > 1 ? 0.5 : 0.0;
                double mean_col = m_windowSize.second > 1 ? 0.5 : 0.0;
 
                // Calculate the filter values
                for (size_t i = 0; i < m_windowSize.first; i++)
                {
                    for (size_t j = 0; j < m_windowSize.second; j++)
                    {
                        if (sqrt(pow2(i/((double)std::max(1u, m_windowSize.first-1))-mean_row) + pow2(j/((double)std::max(1u, m_windowSize.second-1))-mean_col)) <= 0.5)
                        {
                            m_filterKernel[i][j] = fabs(sqrt(pow2(i/((double)std::max(1u, m_windowSize.first-1))-mean_row) + pow2(j/((double)std::max(1u, m_windowSize.second-1))-mean_col)) - 0.5);
                        }
                    }
                }
            }
 
            mu::value_type left(size_t i, size_t j) const
            {
                if (!is2D)
                    return validize(m_left.front());
 
                return validize(m_left[i+1]);
            }
 
            mu::value_type right(size_t i, size_t j) const
            {
                if (!is2D)
                    return validize(m_right.front());
 
                return validize(m_right[i+1]);
            }
 
            mu::value_type top(size_t i, size_t j) const
            {
                if (!is2D)
                    return 0.0;
 
                return validize(m_top[j+1]);
            }
 
            mu::value_type bottom(size_t i, size_t j) const
            {
                if (!is2D)
                    return 0.0;
 
                return validize(m_bottom[j+1]);
            }
 
            mu::value_type topleft(size_t i, size_t j) const
            {
                if (i >= j)
                    return validize(m_left[i-j]);
 
                return validize(m_top[j-i]);
            }
 
            mu::value_type topright(size_t i, size_t j) const
            {
                if (i+j <= m_windowSize.second-1)
                    return validize(m_top[i+j+2]); // size(m_top) + 2 == m_order
 
                return validize(m_right[i+j-m_windowSize.second+1]);
            }
 
            mu::value_type bottomleft(size_t i, size_t j) const
            {
                if (i+j <= m_windowSize.first-1)
                    return validize(m_left[i+j+2]); // size(m_left) + 2 == m_order
 
                return validize(m_bottom[i+j-m_windowSize.first+1]);
            }
 
            mu::value_type bottomright(size_t i, size_t j) const
            {
                if (i >= j)
                    return validize(m_bottom[m_windowSize.second-(i-j)+1]); // size(m_bottom) + 2 == m_order
 
                return validize(m_right[m_windowSize.first-(j-i)+1]);
            }
 
            mu::value_type validize(mu::value_type val) const
            {
                if (mu::isnan(val))
                    return m_fallback;
 
                return val;
            }
 
        public:
            RetouchRegion(size_t _row, size_t _col, const mu::value_type& _dMedian) : Filter(_row, _col)
            {
                m_type = FilterSettings::FILTER_WEIGHTED_LINEAR;
                m_isConvolution = false;
                m_fallback = _dMedian;
                m_invertedKernel = true;
                is2D = true;
 
                createKernel();
            }
 
            virtual ~RetouchRegion() override
            {
                m_filterKernel.clear();
            }
 
            virtual double operator()(size_t i, size_t j) const override
            {
                if (i >= m_windowSize.first || j >= m_windowSize.second)
                    return NAN;
 
                return m_filterKernel[i][j];
            }
 
            virtual mu::value_type apply(size_t i, size_t j, const mu::value_type& val) const override
            {
                if (i >= m_windowSize.first || j >= m_windowSize.second)
                    return NAN;
 
                if (!is2D)
                {
                    if (m_invertedKernel)
                        return m_filterKernel[i][j]*val + (1-m_filterKernel[i][j])*(left(i,j) + (right(i,j) - left(i,j))/(m_windowSize.first + 1.0)*(i + 1.0));
 
                    return (1-m_filterKernel[i][j])*val + m_filterKernel[i][j]*(left(i,j) + (right(i,j) - left(i,j))/(m_windowSize.first + 1.0)*(i + 1.0));
                }
 
                // cross hair: summarize the linearily interpolated values along the rows and cols at the desired position
                // Summarize implies that the value is not averaged yet
                mu::value_type dAverage = top(i,j) + (bottom(i,j) - top(i,j)) / (m_windowSize.first + 1.0) * (i+1.0)
                                        + left(i,j) + (right(i,j) - left(i,j)) / (m_windowSize.second + 1.0) * (j+1.0);
 
                // Additional weighting because are the nearest neighbours
                dAverage /= 2.0;
 
                // Calculate along columns
                // Find the diagonal neighbours and interpolate the value
                if (i >= j)
                    dAverage += topleft(i,j) + (bottomright(i,j) - topleft(i,j)) / (m_windowSize.second - fabs(i-j) + 1.0) * (j+1.0);
                else
                    dAverage += topleft(i,j) + (bottomright(i,j) - topleft(i,j)) / (m_windowSize.first - fabs(i-j) + 1.0) * (i+1.0);
 
                // calculate along rows
                // Find the diagonal neighbours and interpolate the value
                if (i + j <= m_windowSize.first + 1)
                    dAverage += bottomleft(i,j) + (topright(i,j) - bottomleft(i,j)) / (i+j+2.0) * (j+1.0);
                else
                    dAverage += bottomleft(i,j) + (topright(i,j) - bottomleft(i,j)) / (double)(m_windowSize.first + m_windowSize.second - i - j) * (double)(m_windowSize.first-i);
 
                // Restore the desired average
                dAverage /= 6.0;
 
                if (m_invertedKernel)
                    return (1-m_filterKernel[i][j])*dAverage + m_filterKernel[i][j]*val;
 
                return (1-m_filterKernel[i][j])*val + m_filterKernel[i][j]*dAverage;
            }
 
            void setBoundaries(const std::vector<mu::value_type>& left, const std::vector<mu::value_type>& right, const std::vector<mu::value_type>& top = std::vector<mu::value_type>(), const std::vector<mu::value_type>& bottom = std::vector<mu::value_type>())
            {
                m_left = left;
                m_right = right;
                m_top = top;
                m_bottom = bottom;
            }
 
            mu::value_type retouch(size_t i, size_t j, const mu::value_type& val, const mu::value_type& med)
            {
                if (mu::isnan(val) && !mu::isnan(med))
                    return 0.5*(apply(i, j, m_fallback) + med);
                else if (mu::isnan(val) && mu::isnan(med))
                    return apply(i, j, m_fallback);
 
                return val;
            }
    };
}
 
 
 
#endif // FILTERING_HPP