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
3 :
4 : #pragma once
5 :
6 : #include "coupling/filtering/filters/Datastructures.h"
7 : #include "coupling/filtering/interfaces/JunctorInterface.h"
8 : #include <array>
9 :
10 : namespace coupling {
11 : namespace filtering {
12 : template <unsigned int dim, coupling::indexing::IndexTrait... scope> class NLM;
13 : }
14 : } // namespace coupling
15 :
16 : /**
17 : * Noise reduction algorithm using non-local means (NLM) method
18 : * See 'Fast Non Local Means Denoising for 3D MR Images' by Coupé et al. 2006
19 : * and 'Non-Local Means Denoising' by Buades et al. 2011.
20 : *
21 : * @author Piet Jarmatz
22 : *
23 : */
24 : template <unsigned int dim, coupling::indexing::IndexTrait... scope> class coupling::filtering::NLM : public coupling::filtering::JunctorInterface<dim, 2, 1> {
25 : public:
26 : using CellIndex_T = coupling::indexing::CellIndex<dim, coupling::indexing::IndexTrait::vector, scope..., coupling::indexing::IndexTrait::md2macro,
27 : coupling::indexing::IndexTrait::noGhost>;
28 :
29 0 : NLM(const std::vector<coupling::datastructures::CouplingCell<dim>*> inputCellVector_unfiltered,
30 : const std::vector<coupling::datastructures::CouplingCell<dim>*> inputCellVector_prefiltered,
31 : const std::vector<coupling::datastructures::CouplingCell<dim>*> outputCellVector, const std::array<bool, 7> filteredValues, int tws, double sigsq,
32 : double sigsq_rel, double hsq, double hsq_rel, int M = 2, int d = 1)
33 : : coupling::filtering::JunctorInterface<dim, 2, 1>({inputCellVector_unfiltered, inputCellVector_prefiltered}, {outputCellVector}, filteredValues, "NLM"),
34 0 : _timeWindowSize(tws), _M(M), _d(d), _sigsq(sigsq), _sigsq_rel(sigsq_rel), _hsq(hsq), _hsq_rel(hsq_rel), _cycleCounter(0), _t(0),
35 0 : _flowfield(CellIndex_T::numberCellsInDomain, tws), _flowfield_prefiltered(CellIndex_T::numberCellsInDomain, tws),
36 0 : _patchfield(CellIndex_T::numberCellsInDomain - tarch::la::Vector<dim, unsigned int>(2), tws), _innerCellIndices() {
37 : // Initialize flowfield
38 0 : for (int t = 0; t < tws; ++t)
39 0 : for (auto idx : CellIndex_T()) {
40 0 : tarch::la::Vector<dim, unsigned int> idxv(idx.get());
41 0 : coupling::filtering::Quantities<dim>& q = _flowfield(idxv, t);
42 0 : q[0] = 1;
43 0 : for (unsigned int d = 1; d <= dim; ++d)
44 0 : q[d] = 0;
45 0 : coupling::filtering::Quantities<dim>& qp = _flowfield_prefiltered(idxv, t);
46 0 : qp[0] = 1;
47 0 : for (unsigned int d = 1; d <= dim; ++d)
48 0 : qp[d] = 0;
49 : }
50 :
51 : // Initialize innerCellIndices
52 0 : auto domainSize = CellIndex_T::numberCellsInDomain;
53 0 : for (auto idx : CellIndex_T()) {
54 0 : tarch::la::Vector<dim, unsigned int> idxv(idx.get());
55 0 : for (unsigned int d = 0; d < dim; d++)
56 0 : if (idxv[d] == 0 || idxv[d] == domainSize[d] - 1)
57 0 : goto continue_loop;
58 0 : _innerCellIndices.push_back(idx);
59 0 : continue_loop:;
60 : }
61 :
62 : #ifdef NLM_DEBUG
63 : std::cout << " NLM: Created NLM instance." << std::endl;
64 : std::cout << " WARNING: Regardless of configuration, NLM always filters "
65 : "macroscopic mass and momentum."
66 : << std::endl;
67 : #endif
68 0 : }
69 :
70 0 : virtual ~NLM() {
71 : #ifdef NLM_DEBUG
72 : std::cout << " NLM: Destroyed NLM instance." << std::endl;
73 : #endif
74 0 : }
75 :
76 : void operator()();
77 :
78 : private:
79 : /**
80 : * Stores new values from inputCellVectors into _flowfield and
81 : * _flowfield_prefiltered for this timestep _t.
82 : */
83 : void save_data();
84 :
85 : /**
86 : * Reconstructs _patchfield(idx,_t) for all inner cell indices idx,
87 : * so that new information from updated flowfield will be included.
88 : */
89 : void update_patchfield();
90 :
91 : /**
92 : * Main loop of NLM algorithm, computes filtering output for one coupling
93 : * cycle.
94 : */
95 : void denoise();
96 :
97 : /**
98 : * Compute which filtering results we need here.
99 : * @param idx_inner Index of this cell on the inner local MD 2 macro domain,
100 : * i.e. the domain where the patchfield is defined
101 : * @return Relative offsets, of all neighboring boundary cells and this cell
102 : * itself
103 : */
104 : std::vector<tarch::la::Vector<dim, int>> compute_boundary_neighbors(const tarch::la::Vector<dim, unsigned int>&);
105 :
106 : /**
107 : * Must be called to indicate that one coupling cycle is finished and timestep
108 : * counters should be incremented.
109 : */
110 : void increment_time();
111 :
112 : const unsigned int _timeWindowSize; // number of snapshots / coupling cycles taken into
113 : // consideration for noise reduction
114 : // unsigned int _timeModulo; // e.g. 2,4,8 ....
115 : const unsigned int _M; // search volume has size (2M+1)^3
116 : const unsigned int _d; // patches have size (2d+1)^4; this makes at least _d
117 : // ghost layer necessary
118 : double _sigsq, _sigsq_rel;
119 : double _hsq, _hsq_rel;
120 :
121 : unsigned int _cycleCounter; // coupling cycle counter, indicates how many data
122 : // snapshots are available already
123 : unsigned int _t; // active temporal index, iterates cyclic between zero and
124 : // _timeWindowSize
125 : coupling::filtering::Flowfield<dim> _flowfield;
126 : coupling::filtering::Flowfield<dim> _flowfield_prefiltered;
127 : coupling::filtering::Patchfield<dim> _patchfield;
128 : std::vector<CellIndex_T> _innerCellIndices;
129 :
130 0 : inline unsigned int posmod(int i, int n) const { return (i % n + n) % n; }
131 : };
132 :
133 : // include implementation of header
134 : #include "NLM.cpph"
|
