Line data    Source code

       1              : // This file is part of the Mamico project. For conditions of distribution
       2              : // and use, please see the copyright notice in Mamico's main folder, or at
       3              : // www5.in.tum.de/mamico
       4              : 
       5              : #pragma once
       6              : 
       7              : #include <cmath>
       8              : #include <string>
       9              : #include <vector>
      10              : 
      11              : // #define DEBUG_GAUSS
      12              : #include "coupling/filtering/interfaces/FilterInterface.h"
      13              : 
      14              : namespace coupling {
      15              : namespace filtering {
      16              : template <unsigned int dim, coupling::indexing::IndexTrait... scope> class Gauss;
      17              : 
      18              : // cf. member variable in coupling::Gauss for more details
      19              : enum GaussExtrapolationStrategy { NONE, MIRROR, REFLECT };
      20              : } // namespace filtering
      21              : } // namespace coupling
      22              : 
      23              : // Define kernel radius. e.g. radius = 1 means kernel size of 3
      24              : #define GAUSS_KERNEL_RADIUS 1
      25              : 
      26              : /*
      27              :  * Gaussian filter.
      28              :  * Operates in one dimension: If you wish to use a multidimensional gaussian
      29              :  * filter, simply chain multiple instances of this filter in one FilterSequence.
      30              :  *
      31              :  * @author Felix Maurer
      32              :  */
      33              : template <unsigned int dim, coupling::indexing::IndexTrait... scope> class coupling::filtering::Gauss : public coupling::filtering::FilterInterface<dim> {
      34              :   using coupling::filtering::FilterInterface<dim>::_inputCells;
      35              :   using coupling::filtering::FilterInterface<dim>::_outputCells;
      36              :   using coupling::filtering::FilterInterface<dim>::_scalarAccessFunctionPairs;
      37              :   using coupling::filtering::FilterInterface<dim>::_vectorAccessFunctionPairs;
      38              : 
      39              :   using ScalarIndex = coupling::indexing::CellIndex<dim, scope..., coupling::indexing::IndexTrait::md2macro, coupling::indexing::IndexTrait::noGhost>;
      40              :   using VectorIndex = coupling::indexing::CellIndex<dim, coupling::indexing::IndexTrait::vector, scope..., coupling::indexing::IndexTrait::md2macro,
      41              :                                                     coupling::indexing::IndexTrait::noGhost>;
      42              : 
      43              : public:
      44            0 :   Gauss(const std::vector<coupling::datastructures::CouplingCell<dim>*>& inputCellVector,
      45              :         const std::vector<coupling::datastructures::CouplingCell<dim>*>& outputCellVector, const std::array<bool, 7> filteredValues, unsigned int dimension,
      46              :         int sigma, const char* extrapolationStrategy)
      47            0 :       : coupling::filtering::FilterInterface<dim>(inputCellVector, outputCellVector, filteredValues, "GAUSS"), _dim(dimension), _sigma(sigma),
      48            0 :         _kernel(generateKernel()) {
      49              :     // TODO
      50              :     if (GAUSS_KERNEL_RADIUS != 1)
      51              :       throw std::runtime_error("ERROR: GAUSS: Kernel radius != 1 currently not supported.");
      52              : 
      53            0 :     if (extrapolationStrategy == nullptr || std::strcmp(extrapolationStrategy, "none") == 0)
      54            0 :       _extrapolationStrategy = NONE;
      55            0 :     else if (std::strcmp(extrapolationStrategy, "mirror") == 0)
      56            0 :       _extrapolationStrategy = MIRROR;
      57            0 :     else if (std::strcmp(extrapolationStrategy, "reflect") == 0)
      58            0 :       _extrapolationStrategy = REFLECT;
      59              :     else {
      60            0 :       std::cout << "Extrapolation strategy: " << extrapolationStrategy << std::endl;
      61            0 :       throw std::runtime_error("ERROR: GAUSS: Unknown extrapolation strategy.");
      62              :     }
      63              : 
      64              : #ifdef DEBUG_GAUSS
      65              :     std::cout << "           GAUSS (Dim: " << _dim << "): Created Gaussian filter." << std::endl;
      66              :     if (_extrapolationStrategy == NONE)
      67              :       std::cout << "         It will not use extrapolation." << std::endl;
      68              :     if (_extrapolationStrategy == MIRROR)
      69              :       std::cout << "         It will use mirroring extrapolation." << std::endl;
      70              :     if (_extrapolationStrategy == REFLECT)
      71              :       std::cout << "         It will use reflecting extrapolation." << std::endl;
      72              : #endif
      73            0 :   }
      74              : 
      75            0 :   ~Gauss() {
      76              : #ifdef DEBUG_GAUSS
      77              :     std::cout << "           GAUSS (Dim: " << _dim << "): Gaussian filter deconstructed" << std::endl;
      78              : #endif
      79            0 :   }
      80              : 
      81              :   void operator()();
      82              : 
      83              : private:
      84              :   std::array<double, 1 + 2 * GAUSS_KERNEL_RADIUS> generateKernel();
      85              : 
      86              :   constexpr double gaussianDensityFunction(int x);
      87              : 
      88              :   /*
      89              :    * Returns the index of the cell cell that's above the cell at index on the
      90              :    * d-axis If no such index exists, index (the first parameter) is returned.
      91              :    *
      92              :    * Index is assumed to be in terms of the MD2Macro domain, i.e. (0,..0) is the
      93              :    * lowest cell that gets sent from MD to Macro.
      94              :    */
      95              : 
      96              :   VectorIndex getIndexAbove(const VectorIndex index, unsigned int d);
      97              : 
      98              :   /*
      99              :    * Returns the index of the cell that's below the cell at index on the d-axis
     100              :    * If no such index exists, index (the first parameter) is returned.
     101              :    *
     102              :    * Index is assumed to be in terms of the MD2Macro domain, i.e. (0,..0) is the
     103              :    * lowest cell that gets sent from MD to Macro.
     104              :    */
     105              : 
     106              :   VectorIndex getIndexBelow(const VectorIndex index, unsigned int d);
     107              : 
     108              :   // on which axis this filter operates. 0 <= _dim <= dim
     109              :   const unsigned int _dim;
     110              : 
     111              :   // standard deviation used
     112              :   const double _sigma;
     113              : 
     114              :   std::array<double, 1 + 2 * GAUSS_KERNEL_RADIUS> _kernel;
     115              : 
     116              :   /**
     117              :    * Determines how to apply filter to border cells:
     118              :    * NONE = only use existing cells and normalize their weight accordingly
     119              :    * MIRROR = b | a b c d | c
     120              :    * REFLECT = a | a b c d | d
     121              :    *
     122              :    * The last two are congruent to SciPy's gaussian filter's respective
     123              :    * extrapolation modes
     124              :    */
     125              :   coupling::filtering::GaussExtrapolationStrategy _extrapolationStrategy;
     126              : };
     127              : 
     128              : // include implementation of header
     129              : #include "Gauss.cpph"