Fixed a stupid mistake that the order of applying DCT and IDCT is reversed

Fixed a stupid mistake that the order of applying DCT and IDCT is reversed (it should be DCT->filter->IDCT), which leads to over-filtering in fine structures and produces artifacts such as blocking, ringing and desaturation. Now it acctually filters as expected.
Adjust hard-thresholding table to fit the unnormalized DCT, to be equivalent to filtering in orthogonal transform
Default value of hard_thr changed from 2.3 to 2.8 for profile="vn", from 2.2 to 2.7 for other profiles

Some function arguments changed to const reference
Some read-only references in function arguments changed to const
Some read-only member functions changed to const

Cosmetics

Author: mawen1250

README.md

@@ -104,6 +104,7 @@ bm3d.Basic(clip input[, clip ref=input, string profile="fast", float[] sigma=[10
 
 - sigma:<br />
     The strength of denoising, valid range [0, +inf), default [10,10,10].<br />
+    This denoising algorithm is very sensitive to the sigma, so adjust it carefully according to the source.<br />
     Technically, this is the standard deviation of i.i.d. zero mean additive white Gaussian noise in 8 bit scale. BM3D denoising filter is designed based on this noise model, and best fit for attenuating it.<br />
     An array up to 3 elements can be assigned to set different sigma for Y,U,V channels. If less than 3 elements assigned, the last element's value will be assigned to the undefined elements.<br />
     0 will disable the processing for corresponding channel.
@@ -141,7 +142,7 @@ bm3d.Basic(clip input[, clip ref=input, string profile="fast", float[] sigma=[10
 
 - hard_thr:<br />
     The threshold parameter for the hard-thresholding in 3D transformed domain, in 8 bit scale, valid range (0, +inf).<br />
-    Larger results in stronger hard-threshold filtering in frequency domain.<br />
+    Larger value results in stronger hard-threshold filtering in frequency domain.<br />
     Usually, to tweak denoising strength, it's better to adjust "sigma" rather than "hard_thr".
 
 - matrix:<br />
@@ -296,7 +297,7 @@ bm3d.Basic / bm3d.Final / bm3d.VBasic / bm3d.VFinal
 bm3d.VBasic / bm3d.VFinal
 ---------------------------------------------------
 | profile || radius | ps_num | ps_range | ps_step |
---------------------------------------------------
+---------------------------------------------------
 | "fast"  || 1/1    | 2/2    | 4/5      | 1/1/1/1 |
 | "lc"    || 2/2    | 2/2    | 4/5      | 1/1/1/1 |
 | "np"    || 3/3    | 2/2    | 5/6      | 1/1/1/1 |
@@ -310,11 +311,11 @@ bm3d.Basic & bm3d.VBasic / bm3d.Final & bm3d.VFinal
 --------------------------------------------------------------
 | profile || th_mse                              | hard_thr  |
 --------------------------------------------------------------
-| "fast"  || sigma[0]*80+400   / sigma[0]*10+200 | 2.2 / NUL |
-| "lc"    || sigma[0]*80+400   / sigma[0]*10+200 | 2.2 / NUL |
-| "np"    || sigma[0]*80+400   / sigma[0]*10+200 | 2.2 / NUL |
-| "high"  || sigma[0]*80+400   / sigma[0]*10+200 | 2.2 / NUL |
-| "vn"    || sigma[0]*150+1000 / sigma[0]*40+400 | 2.3 / NUL |
+| "fast"  || sigma[0]*80+400   / sigma[0]*10+200 | 2.7 / NUL |
+| "lc"    || sigma[0]*80+400   / sigma[0]*10+200 | 2.7 / NUL |
+| "np"    || sigma[0]*80+400   / sigma[0]*10+200 | 2.7 / NUL |
+| "high"  || sigma[0]*80+400   / sigma[0]*10+200 | 2.7 / NUL |
+| "vn"    || sigma[0]*150+1000 / sigma[0]*40+400 | 2.8 / NUL |
 --------------------------------------------------------------
 ```
 

include/BM3D_Base.h

@@ -230,18 +230,18 @@ class BM3D_Process_Base
         _NewFrame(width, height, dfi == fi);
     }
 
-    void Kernel(FLType *dst, const FLType *src, const FLType *ref);
+    void Kernel(FLType *dst, const FLType *src, const FLType *ref) const;
 
     void Kernel(FLType *dstY, FLType *dstU, FLType *dstV,
         const FLType *srcY, const FLType *srcU, const FLType *srcV,
-        const FLType *refY, const FLType *refU, const FLType *refV);
+        const FLType *refY, const FLType *refU, const FLType *refV) const;
 
-    PosPairCode BlockMatching(const FLType *ref, PCType j, PCType i);
+    PosPairCode BlockMatching(const FLType *ref, PCType j, PCType i) const;
 
     virtual void CollaborativeFilter(int plane,
         FLType *ResNum, FLType *ResDen,
         const FLType *src, const FLType *ref,
-        const PosPairCode &code) = 0;
+        const PosPairCode &code) const = 0;
 };
 
 

include/BM3D_Basic.h

@@ -85,7 +85,7 @@ class BM3D_Basic_Process
     virtual void CollaborativeFilter(int plane,
         FLType *ResNum, FLType *ResDen,
         const FLType *src, const FLType *ref,
-        const PosPairCode &code) override;
+        const PosPairCode &code) const override;
 };
 
 

include/BM3D_Final.h

@@ -85,7 +85,7 @@ class BM3D_Final_Process
     virtual void CollaborativeFilter(int plane,
         FLType *ResNum, FLType *ResDen,
         const FLType *src, const FLType *ref,
-        const PosPairCode &code) override;
+        const PosPairCode &code) const override;
 };
 
 

include/Block.h

@@ -76,15 +76,15 @@ class Block
     }
 
     template < typename _St2, typename _Fn1 >
-    void for_each(_St2 &right, _Fn1 _Func)
+    void for_each(_St2 &data2, _Fn1 _Func)
     {
-        Block_For_each(*this, right, _Func);
+        Block_For_each(*this, data2, _Func);
     }
 
     template < typename _St2, typename _Fn1 >
-    void for_each(_St2 &right, _Fn1 _Func) const
+    void for_each(_St2 &data2, _Fn1 _Func) const
     {
-        Block_For_each(*this, right, _Func);
+        Block_For_each(*this, data2, _Func);
     }
 
     template < typename _Fn1 >
@@ -103,7 +103,7 @@ class Block
     // Default constructor
     Block() {}
 
-    Block(PCType _Height, PCType _Width, PosType pos, bool Init = true, value_type Value = 0)
+    Block(PCType _Height, PCType _Width, const PosType &pos, bool Init = true, value_type Value = 0)
         : Height_(_Height), Width_(_Width), PixelCount_(Height_ * Width_), pos_(pos)
     {
         AlignedMalloc(Data_, size());
@@ -113,13 +113,13 @@ class Block
 
     // Constructor from plane pointer and PosType
     template < typename _St1 >
-    Block(const _St1 *src, PCType src_stride, PCType _Height, PCType _Width, PosType pos)
+    Block(const _St1 *src, PCType src_stride, PCType _Height, PCType _Width, const PosType &pos)
         : Block(_Height, _Width, pos, false)
     {
         From(src, src_stride);
     }
 
-    // Constructor from src
+    // Constructor from src Block
     Block(const _Myt &src, bool Init, value_type Value = 0)
         : Block(src.Height_, src.Width_, src.pos_, Init, Value)
     {}
@@ -275,10 +275,9 @@ class Block
     }
 
     template < typename _St1 >
-    void From(const _St1 *src, PCType src_stride, PosType pos)
+    void From(const _St1 *src, PCType src_stride, const PosType &pos)
     {
         pos_ = pos;
-
         From(src, src_stride);
     }
 
@@ -303,7 +302,7 @@ class Block
     // Accumulate functions
 
     template < typename _St1 >
-    void AddFrom(const _St1 *src, PCType src_stride, PosType pos)
+    void AddFrom(const _St1 *src, PCType src_stride, const PosType &pos)
     {
         auto dstp = data();
         auto srcp = src + pos.y * src_stride + pos.x;
@@ -320,7 +319,7 @@ class Block
     }
 
     template < typename _St1, typename _Gt1 >
-    void AddFrom(const _St1 *src, PCType src_stride, PosType pos, _Gt1 gain)
+    void AddFrom(const _St1 *src, PCType src_stride, const PosType &pos, _Gt1 gain)
     {
         auto dstp = data();
         auto srcp = src + pos.y * src_stride + pos.x;
@@ -403,7 +402,7 @@ class Block
     }
 
     ////////////////////////////////////////////////////////////////
-    // Block matching functions
+    // Single block-matching functions
 
     template < typename _St1 >
     PosPair BlockMatching(const _St1 *src, PCType src_height, PCType src_width, PCType src_stride, _St1 src_range,
@@ -474,6 +473,9 @@ class Block
         return PosPair(static_cast<KeyType>(distMin * distMul), pos);
     }
 
+    ////////////////////////////////////////////////////////////////
+    // Multiple block-matching functions
+
     template < typename _St1 >
     void BlockMatchingMulti(PosPairCode &match_code, const _St1 *src, PCType src_stride, _St1 src_range,
         const PosCode &search_pos, double thMSE) const
@@ -798,7 +800,7 @@ class BlockGroup
 
     // Constructor from plane pointer and Pos3PairCode
     template < typename _St1 >
-    BlockGroup(std::vector<const _St1 *> src, PCType src_stride, const Pos3PairCode &code,
+    BlockGroup(const std::vector<const _St1 *> &src, PCType src_stride, const Pos3PairCode &code,
         PCType _GroupSize = -1, PCType _Height = 16, PCType _Width = 16)
         : Height_(_Height), Width_(_Width)
     {
@@ -984,7 +986,7 @@ class BlockGroup
     // Read/Store functions
 
     template < typename _St1 >
-    _Myt &From(const _St1 *src, PCType src_stride)
+    void From(const _St1 *src, PCType src_stride)
     {
         auto dstp = data();
 
@@ -1002,12 +1004,10 @@ class BlockGroup
                 }
             }
         }
-
-        return *this;
     }
 
     template < typename _St1 >
-    _Myt &From(std::vector<const _St1 *> src, PCType src_stride)
+    void From(const std::vector<const _St1 *> &src, PCType src_stride)
     {
         auto dstp = data();
 
@@ -1025,8 +1025,6 @@ class BlockGroup
                 }
             }
         }
-
-        return *this;
     }
 
     template < typename _Dt1 >
@@ -1051,7 +1049,7 @@ class BlockGroup
     }
 
     template < typename _Dt1 >
-    void To(std::vector<_Dt1 *> dst, PCType dst_stride) const
+    void To(const std::vector<_Dt1 *> &dst, PCType dst_stride) const
     {
         auto srcp = data();
 
@@ -1117,7 +1115,7 @@ class BlockGroup
     }
 
     template < typename _Dt1 >
-    void AddTo(std::vector<_Dt1 *> dst, PCType dst_stride) const
+    void AddTo(const std::vector<_Dt1 *> &dst, PCType dst_stride) const
     {
         auto srcp = data();
 
@@ -1138,7 +1136,7 @@ class BlockGroup
     }
 
     template < typename _Dt1, typename _Gt1 >
-    void AddTo(std::vector<_Dt1 *> dst, PCType dst_stride, _Gt1 gain) const
+    void AddTo(const std::vector<_Dt1 *> &dst, PCType dst_stride, _Gt1 gain) const
     {
         auto srcp = data();
 
@@ -1197,7 +1195,7 @@ class BlockGroup
     }
 
     template < typename _Dt1 >
-    void CountTo(std::vector<_Dt1 *> dst, PCType dst_stride) const
+    void CountTo(const std::vector<_Dt1 *> &dst, PCType dst_stride) const
     {
         for (PCType z = 0; z < GroupSize(); ++z)
         {
@@ -1216,7 +1214,7 @@ class BlockGroup
     }
 
     template < typename _Dt1 >
-    void CountTo(std::vector<_Dt1 *> dst, PCType dst_stride, _Dt1 value) const
+    void CountTo(const std::vector<_Dt1 *> &dst, PCType dst_stride, _Dt1 value) const
     {
         for (PCType z = 0; z < GroupSize(); ++z)
         {
@@ -1239,7 +1237,7 @@ class BlockGroup
 ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
 
 
-template < typename _St1, typename _Fn1 > inline
+template < typename _St1, typename _Fn1 >
 void Block_For_each(_St1 &data, _Fn1 &&_Func)
 {
     auto datap = data.data();
@@ -1250,26 +1248,24 @@ void Block_For_each(_St1 &data, _Fn1 &&_Func)
     }
 }
 
-template < typename _St1, typename _St2, typename _Fn1 > inline
-void Block_For_each(_St1 &left, _St2 &right, _Fn1 &&_Func)
+template < typename _St1, typename _St2, typename _Fn1 >
+void Block_For_each(_St1 &data1, _St2 &data2, _Fn1 &&_Func)
 {
-    const char *FunctionName = "Block_For_each";
-    if (left.size() != right.size())
+    if (data1.size() != data2.size())
     {
-        std::cerr << FunctionName << ": size() of left and right must be the same.\n";
-        exit(EXIT_FAILURE);
+        DEBUG_FAIL("Block_For_each: size() of data1 and data2 must be the same.");
     }
 
-    auto leftp = left.data();
-    auto rightp = right.data();
+    auto data1p = data1.data();
+    auto data2p = data2.data();
 
-    for (auto upper = leftp + left.size(); leftp != upper; ++leftp, ++rightp)
+    for (auto upper = data1p + data1.size(); data1p != upper; ++data1p, ++data2p)
     {
-        _Func(*leftp, *rightp);
+        _Func(*data1p, *data2p);
     }
 }
 
-template < typename _St1, typename _Fn1 > inline
+template < typename _St1, typename _Fn1 >
 void Block_Transform(_St1 &data, _Fn1 &&_Func)
 {
     auto datap = data.data();
@@ -1280,14 +1276,12 @@ void Block_Transform(_St1 &data, _Fn1 &&_Func)
     }
 }
 
-template < typename _Dt1, typename _St1, typename _Fn1 > inline
+template < typename _Dt1, typename _St1, typename _Fn1 >
 void Block_Transform(_Dt1 &dst, const _St1 &src, _Fn1 &&_Func)
 {
-    const char *FunctionName = "Block_Transform";
     if (dst.size() != src.size())
     {
-        std::cerr << FunctionName << ": size() of dst and src must be the same.\n";
-        exit(EXIT_FAILURE);
+        DEBUG_FAIL("Block_Transform: size() of dst and src must be the same.");
     }
 
     auto dstp = dst.data();