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"