Merge pull request #733 from bpawlik:bm3d

modules/xphoto/include/opencv2/xphoto.hpp

@@ -49,4 +49,5 @@
 #include "xphoto/inpainting.hpp"
 #include "xphoto/white_balance.hpp"
 #include "xphoto/dct_image_denoising.hpp"
+#include "xphoto/bm3d_image_denoising.hpp"
 #endif

modules/xphoto/include/opencv2/xphoto/bm3d_image_denoising.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009-2011, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_BM3D_IMAGE_DENOISING_HPP__
+#define __OPENCV_BM3D_IMAGE_DENOISING_HPP__
+
+/** @file
+@date Jul 19, 2016
+@author Bartek Pawlik
+*/
+
+#include <opencv2/core.hpp>
+
+namespace cv
+{
+    namespace xphoto
+    {
+        //! @addtogroup xphoto
+        //! @{
+
+        //! BM3D transform types
+        enum TransformTypes
+        {
+            /** Un-normalized Haar transform */
+            HAAR = 0
+        };
+
+        //! BM3D algorithm steps
+        enum Bm3dSteps
+        {
+            /** Execute all steps of the algorithm */
+            BM3D_STEPALL = 0,
+            /** Execute only first step of the algorithm */
+            BM3D_STEP1 = 1,
+            /** Execute only second step of the algorithm */
+            BM3D_STEP2 = 2
+        };
+
+        /** @brief Performs image denoising using the Block-Matching and 3D-filtering algorithm
+        <http://www.cs.tut.fi/~foi/GCF-BM3D/BM3D_TIP_2007.pdf> with several computational
+        optimizations. Noise expected to be a gaussian white noise.
+
+        @param src Input 8-bit or 16-bit 1-channel image.
+        @param dstStep1 Output image of the first step of BM3D with the same size and type as src.
+        @param dstStep2 Output image of the second step of BM3D with the same size and type as src.
+        @param h Parameter regulating filter strength. Big h value perfectly removes noise but also
+        removes image details, smaller h value preserves details but also preserves some noise.
+        @param templateWindowSize Size in pixels of the template patch that is used for block-matching.
+        Should be power of 2.
+        @param searchWindowSize Size in pixels of the window that is used to perform block-matching.
+        Affect performance linearly: greater searchWindowsSize - greater denoising time.
+        Must be larger than templateWindowSize.
+        @param blockMatchingStep1 Block matching threshold for the first step of BM3D (hard thresholding),
+        i.e. maximum distance for which two blocks are considered similar.
+        Value expressed in euclidean distance.
+        @param blockMatchingStep2 Block matching threshold for the second step of BM3D (Wiener filtering),
+        i.e. maximum distance for which two blocks are considered similar.
+        Value expressed in euclidean distance.
+        @param groupSize Maximum size of the 3D group for collaborative filtering.
+        @param slidingStep Sliding step to process every next reference block.
+        @param beta Kaiser window parameter that affects the sidelobe attenuation of the transform of the
+        window. Kaiser window is used in order to reduce border effects. To prevent usage of the window,
+        set beta to zero.
+        @param normType Norm used to calculate distance between blocks. L2 is slower than L1
+        but yields more accurate results.
+        @param step Step of BM3D to be executed. Possible variants are: step 1, step 2, both steps.
+        @param transformType Type of the orthogonal transform used in collaborative filtering step.
+        Currently only Haar transform is supported.
+
+        This function expected to be applied to grayscale images. Advanced usage of this function
+        can be manual denoising of colored image in different colorspaces.
+
+        @sa
+        fastNlMeansDenoising
+        */
+        CV_EXPORTS_W void bm3dDenoising(
+            InputArray src,
+            InputOutputArray dstStep1,
+            OutputArray dstStep2,
+            float h = 1,
+            int templateWindowSize = 4,
+            int searchWindowSize = 16,
+            int blockMatchingStep1 = 2500,
+            int blockMatchingStep2 = 400,
+            int groupSize = 8,
+            int slidingStep = 1,
+            float beta = 2.0f,
+            int normType = cv::NORM_L2,
+            int step = cv::xphoto::BM3D_STEPALL,
+            int transformType = cv::xphoto::HAAR);
+
+        /** @brief Performs image denoising using the Block-Matching and 3D-filtering algorithm
+        <http://www.cs.tut.fi/~foi/GCF-BM3D/BM3D_TIP_2007.pdf> with several computational
+        optimizations. Noise expected to be a gaussian white noise.
+
+        @param src Input 8-bit or 16-bit 1-channel image.
+        @param dst Output image with the same size and type as src.
+        @param h Parameter regulating filter strength. Big h value perfectly removes noise but also
+        removes image details, smaller h value preserves details but also preserves some noise.
+        @param templateWindowSize Size in pixels of the template patch that is used for block-matching.
+        Should be power of 2.
+        @param searchWindowSize Size in pixels of the window that is used to perform block-matching.
+        Affect performance linearly: greater searchWindowsSize - greater denoising time.
+        Must be larger than templateWindowSize.
+        @param blockMatchingStep1 Block matching threshold for the first step of BM3D (hard thresholding),
+        i.e. maximum distance for which two blocks are considered similar.
+        Value expressed in euclidean distance.
+        @param blockMatchingStep2 Block matching threshold for the second step of BM3D (Wiener filtering),
+        i.e. maximum distance for which two blocks are considered similar.
+        Value expressed in euclidean distance.
+        @param groupSize Maximum size of the 3D group for collaborative filtering.
+        @param slidingStep Sliding step to process every next reference block.
+        @param beta Kaiser window parameter that affects the sidelobe attenuation of the transform of the
+        window. Kaiser window is used in order to reduce border effects. To prevent usage of the window,
+        set beta to zero.
+        @param normType Norm used to calculate distance between blocks. L2 is slower than L1
+        but yields more accurate results.
+        @param step Step of BM3D to be executed. Allowed are only BM3D_STEP1 and BM3D_STEPALL.
+        BM3D_STEP2 is not allowed as it requires basic estimate to be present.
+        @param transformType Type of the orthogonal transform used in collaborative filtering step.
+        Currently only Haar transform is supported.
+
+        This function expected to be applied to grayscale images. Advanced usage of this function
+        can be manual denoising of colored image in different colorspaces.
+
+        @sa
+        fastNlMeansDenoising
+        */
+        CV_EXPORTS_W void bm3dDenoising(
+            InputArray src,
+            OutputArray dst,
+            float h = 1,
+            int templateWindowSize = 4,
+            int searchWindowSize = 16,
+            int blockMatchingStep1 = 2500,
+            int blockMatchingStep2 = 400,
+            int groupSize = 8,
+            int slidingStep = 1,
+            float beta = 2.0f,
+            int normType = cv::NORM_L2,
+            int step = cv::xphoto::BM3D_STEPALL,
+            int transformType = cv::xphoto::HAAR);
+        //! @}
+    }
+}
+
+#endif // __OPENCV_BM3D_IMAGE_DENOISING_HPP__

modules/xphoto/samples/bm3d_image_denoising.cpp

+#include "opencv2/xphoto.hpp"
+#include "opencv2/highgui.hpp"
+
+const char* keys =
+{
+    "{i || input image name}"
+    "{o || output image name}"
+    "{sigma || expected noise standard deviation}"
+    "{tw |4| template window size}"
+    "{sw |16| search window size}"
+};
+
+int main(int argc, const char** argv)
+{
+    bool printHelp = (argc == 1);
+    printHelp = printHelp || (argc == 2 && std::string(argv[1]) == "--help");
+    printHelp = printHelp || (argc == 2 && std::string(argv[1]) == "-h");
+
+    if (printHelp)
+    {
+        printf("\nThis sample demonstrates BM3D image denoising\n"
+            "Call:\n"
+            "    bm3d_image_denoising -i=<string> -sigma=<double> -tw=<int> -sw=<int> [-o=<string>]\n\n");
+        return 0;
+    }
+
+    cv::CommandLineParser parser(argc, argv, keys);
+    if (!parser.check())
+    {
+        parser.printErrors();
+        return -1;
+    }
+
+    std::string inFilename = parser.get<std::string>("i");
+    std::string outFilename = parser.get<std::string>("o");
+
+    cv::Mat src = cv::imread(inFilename, cv::IMREAD_GRAYSCALE);
+    if (src.empty())
+    {
+        printf("Cannot read image file: %s\n", inFilename.c_str());
+        return -1;
+    }
+
+    float sigma = parser.get<float>("sigma");
+    if (sigma == 0.0)
+        sigma = 15.0;
+
+    int templateWindowSize = parser.get<int>("tw");
+    if (templateWindowSize == 0)
+        templateWindowSize = 4;
+
+    int searchWindowSize = parser.get<int>("sw");
+    if (searchWindowSize == 0)
+        searchWindowSize = 16;
+
+    cv::Mat res(src.size(), src.type());
+    cv::xphoto::bm3dDenoising(src, res, sigma, templateWindowSize, searchWindowSize);
+
+    if (outFilename.empty())
+    {
+        cv::namedWindow("input image", cv::WINDOW_NORMAL);
+        cv::imshow("input image", src);
+        cv::namedWindow("denoising result", cv::WINDOW_NORMAL);
+        cv::imshow("denoising result", res);
+        cv::waitKey(0);
+    }
+    else
+    {
+        cv::imwrite(outFilename, res);
+    }
+
+    return 0;
+}

modules/xphoto/src/bm3d_denoising_invoker_commons.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective icvers.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_BM3D_DENOISING_INVOKER_COMMONS_HPP__
+#define __OPENCV_BM3D_DENOISING_INVOKER_COMMONS_HPP__
+
+#include "bm3d_denoising_invoker_structs.hpp"
+
+// std::isnan is a part of C++11 and it is not supported in MSVS2010/2012
+#if defined _MSC_VER && _MSC_VER < 1800 /* MSVC 2013 */
+#include <float.h>
+namespace std {
+    template <typename T> bool isnan(T value) { return _isnan(value) != 0; }
+}
+#endif
+
+namespace cv
+{
+namespace xphoto
+{
+
+// Returns largest power of 2 smaller than the input value
+inline int getLargestPowerOf2SmallerThan(unsigned x)
+{
+    x = x | (x >> 1);
+    x = x | (x >> 2);
+    x = x | (x >> 4);
+    x = x | (x >> 8);
+    x = x | (x >> 16);
+    return x - (x >> 1);
+}
+
+// Returns true if x is a power of 2. Otherwise false.
+inline bool isPowerOf2(int x)
+{
+    return (x > 0) && !(x & (x - 1));
+}
+
+
+template <typename T>
+inline static void shrink(T &val, T &nonZeroCount, const T &threshold)
+{
+    if (std::abs(val) < threshold)
+        val = 0;
+    else
+        ++nonZeroCount;
+}
+
+template <typename T>
+inline static void hardThreshold2D(T *dst, T *thrMap, const int &templateWindowSizeSq)
+{
+    for (int i = 1; i < templateWindowSizeSq; ++i)
+    {
+        if (std::abs(dst[i] < thrMap[i]))
+            dst[i] = 0;
+    }
+}
+
+template <int N, typename T, typename DT, typename CT>
+inline static T HardThreshold(BlockMatch<T, DT, CT> *z, const int &n, T *&thrMap)
+{
+    T nonZeroCount = 0;
+
+    for (int i = 0; i < N; ++i)
+        shrink(z[i][n], nonZeroCount, *thrMap++);
+
+    return nonZeroCount;
+}
+
+template <typename T, typename DT, typename CT>
+inline static T HardThreshold(BlockMatch<T, DT, CT> *z, const int &n, T *&thrMap, const int &N)
+{
+    T nonZeroCount = 0;
+
+    for (int i = 0; i < N; ++i)
+        shrink(z[i][n], nonZeroCount, *thrMap++);
+
+    return nonZeroCount;
+}
+
+template <int N, typename T, typename DT, typename CT>
+inline static int WienerFiltering(BlockMatch<T, DT, CT> *zSrc, BlockMatch<T, DT, CT> *zBasic, const int &n, T *&thrMap)
+{
+    int wienerCoeffs = 0;
+
+    for (int i = 0; i < N; ++i)
+    {
+        // Possible optimization point here to get rid of floats and casts
+        int basicSq = zBasic[i][n] * zBasic[i][n];
+        int sigmaSq = *thrMap * *thrMap;
+        int denom = basicSq + sigmaSq;
+        float wie = (denom == 0) ? 1.0f : ((float)basicSq / (float)denom);
+
+        zBasic[i][n] = (T)(zSrc[i][n] * wie);
+        wienerCoeffs += (int)wie;
+        ++thrMap;
+    }
+
+    return wienerCoeffs;
+}
+
+template <typename T, typename DT, typename CT>
+inline static int WienerFiltering(BlockMatch<T, DT, CT> *zSrc, BlockMatch<T, DT, CT> *zBasic, const int &n, T *&thrMap, const unsigned &N)
+{
+    int wienerCoeffs = 0;
+
+    for (unsigned i = 0; i < N; ++i)
+    {
+        // Possible optimization point here to get rid of floats and casts
+        int basicSq = zBasic[i][n] * zBasic[i][n];
+        int sigmaSq = *thrMap * *thrMap;
+        int denom = basicSq + sigmaSq;
+        float wie = (denom == 0) ? 1.0f : ((float)basicSq / (float)denom);
+
+        zBasic[i][n] = (T)(zSrc[i][n] * wie);
+        wienerCoeffs += (int)wie;
+        ++thrMap;
+    }
+
+    return wienerCoeffs;
+}
+
+
+}  // namespace xphoto
+}  // namespace cv
+
+#endif

modules/xphoto/src/bm3d_denoising_invoker_step1.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective icvers.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_BM3D_DENOISING_INVOKER_STEP1_HPP__
+#define __OPENCV_BM3D_DENOISING_INVOKER_STEP1_HPP__
+
+#include "bm3d_denoising_invoker_commons.hpp"
+#include "bm3d_denoising_transforms.hpp"
+#include "kaiser_window.hpp"
+
+namespace cv
+{
+namespace xphoto
+{
+
+template <typename T, typename D, typename WT, typename TT, typename TC>
+struct Bm3dDenoisingInvokerStep1 : public ParallelLoopBody
+{
+public:
+    Bm3dDenoisingInvokerStep1(
+        const Mat& src,
+        Mat& dst,
+        const int &templateWindowSize,
+        const int &searchWindowSize,
+        const float &h,
+        const int &hBM,
+        const int &groupSize,
+        const int &slidingStep,
+        const float &beta);
+
+    virtual ~Bm3dDenoisingInvokerStep1();
+    void operator() (const Range& range) const;
+
+private:
+    // Unimplemented operator in order to satisfy compiler warning.
+    void operator= (const Bm3dDenoisingInvokerStep1&);
+
+    void calcDistSumsForFirstElementInRow(
+        int i,
+        Array2d<int>& distSums,
+        Array3d<int>& colDistSums,
+        Array3d<int>& lastColDistSums,
+        BlockMatch<TT, int, TT> *bm,
+        int &elementSize) const;
+
+    void calcDistSumsForAllElementsInFirstRow(
+        int i,
+        int j,
+        int firstColNum,
+        Array2d<int>& distSums,
+        Array3d<int>& colDistSums,
+        Array3d<int>& lastColDistSums,
+        BlockMatch<TT, int, TT> *bm,
+        int &elementSize) const;
+
+    // Image containers
+    const Mat& src_;
+    Mat& dst_;
+    Mat srcExtended_;
+
+    // Border size of the extended src and basic images
+    int borderSize_;
+
+    // Template and window size
+    int templateWindowSize_;
+    int searchWindowSize_;
+
+    // Half template and window size
+    int halfTemplateWindowSize_;
+    int halfSearchWindowSize_;
+
+    // Squared template and window size
+    int templateWindowSizeSq_;
+    int searchWindowSizeSq_;
+
+    // Block matching threshold
+    int hBM_;
+
+    // Maximum size of 3D group
+    int groupSize_;
+
+    // Sliding step
+    const int slidingStep_;
+
+    // Threshold map
+    TT *thrMap_;
+
+    // Kaiser window
+    float *kaiser_;
+};
+
+template <typename T, typename D, typename WT, typename TT, typename TC>
+Bm3dDenoisingInvokerStep1<T, D, WT, TT, TC>::Bm3dDenoisingInvokerStep1(
+    const Mat& src,
+    Mat& dst,
+    const int &templateWindowSize,
+    const int &searchWindowSize,
+    const float &h,
+    const int &hBM,
+    const int &groupSize,
+    const int &slidingStep,
+    const float &beta) :
+    src_(src), dst_(dst), groupSize_(groupSize), slidingStep_(slidingStep), thrMap_(NULL), kaiser_(NULL)
+{
+    groupSize_ = getLargestPowerOf2SmallerThan(groupSize);
+    CV_Assert(groupSize > 0);
+
+    halfTemplateWindowSize_ = templateWindowSize >> 1;
+    halfSearchWindowSize_ = searchWindowSize >> 1;
+    templateWindowSize_ = templateWindowSize;
+    searchWindowSize_ = searchWindowSize;
+    templateWindowSizeSq_ = templateWindowSize_ * templateWindowSize_;
+    searchWindowSizeSq_ = searchWindowSize_ * searchWindowSize_;
+
+    // Extend image to avoid border problem
+    borderSize_ = halfSearchWindowSize_ + halfTemplateWindowSize_;
+    copyMakeBorder(src_, srcExtended_, borderSize_, borderSize_, borderSize_, borderSize_, BORDER_DEFAULT);
+
+    // Calculate block matching threshold
+    hBM_ = D::template calcBlockMatchingThreshold<int>(hBM, templateWindowSizeSq_);
+
+    // Select transforms depending on the template size
+    TC::RegisterTransforms2D(templateWindowSize_);
+
+    // Precompute threshold map
+    TC::calcThresholdMap3D(thrMap_, h, templateWindowSize_, groupSize_);
+
+    // Generate kaiser window
+    calcKaiserWindow2D(kaiser_, templateWindowSize_, beta);
+}
+
+template<typename T, typename D, typename WT, typename TT, typename TC>
+inline Bm3dDenoisingInvokerStep1<T, D, WT, TT, TC>::~Bm3dDenoisingInvokerStep1()
+{
+    delete[] thrMap_;
+    delete[] kaiser_;
+}
+
+template <typename T, typename D, typename WT, typename TT, typename TC>
+void Bm3dDenoisingInvokerStep1<T, D, WT, TT, TC>::operator() (const Range& range) const
+{
+    const int size = (range.size() + 2 * borderSize_) * srcExtended_.cols;
+    std::vector<WT> weightedSum(size, 0.0);
+    std::vector<WT> weights(size, 0.0);
+    int row_from = range.start;
+    int row_to = range.end - 1;
+
+    // Local vars for faster processing
+    const int blockSize = templateWindowSize_;
+    const int blockSizeSq = templateWindowSizeSq_;
+    const int halfBlockSize = halfTemplateWindowSize_;
+    const int searchWindowSize = searchWindowSize_;
+    const int searchWindowSizeSq = searchWindowSizeSq_;
+    const TT halfSearchWindowSize = (TT)halfSearchWindowSize_;
+    const int hBM = hBM_;
+    const int groupSize = groupSize_;
+
+    const int step = srcExtended_.cols;
+    const int dstStep = srcExtended_.cols;
+    const int weiStep = srcExtended_.cols;
+    const int dstcstep = dstStep - blockSize;
+    const int weicstep = weiStep - blockSize;
+
+    // Buffer to store 3D group
+    BlockMatch<TT, int, TT> *bm = new BlockMatch<TT, int, TT>[searchWindowSizeSq];
+    for (int i = 0; i < searchWindowSizeSq; ++i)
+        bm[i].init(blockSizeSq);
+
+    // First element in a group is always the reference patch. Hence distance is 0.
+    bm[0](0, halfSearchWindowSize, halfSearchWindowSize);
+
+    // Sums of columns and rows for current pixel
+    Array2d<int> distSums(searchWindowSize, searchWindowSize);
+
+    // Sums of columns for current pixel (for lazy calc optimization)
+    Array3d<int> colDistSums(blockSize, searchWindowSize, searchWindowSize);
+
+    // Last elements of column sum (for each element in a row)
+    Array3d<int> lastColDistSums(src_.cols, searchWindowSize, searchWindowSize);
+
+    int firstColNum = -1;
+    for (int j = row_from, jj = 0; j <= row_to; j += slidingStep_, jj += slidingStep_)
+    {
+        for (int i = 0; i < src_.cols; i += slidingStep_)
+        {
+            const T *currentPixel = srcExtended_.ptr<T>(0) + step*j + i;
+            int elementSize = 1;
+
+            // Calculate distSums using moving average filter approach.
+            if (i == 0)
+            {
+                // Calculate distSums for the first element in a row
+                calcDistSumsForFirstElementInRow(j, distSums, colDistSums, lastColDistSums, bm, elementSize);
+                firstColNum = 0;
+            }
+            else
+            {
+                if (j == row_from)
+                {
+                    // Calculate distSums for all elements in the first row
+                    calcDistSumsForAllElementsInFirstRow(
+                        j, i, firstColNum, distSums, colDistSums, lastColDistSums, bm, elementSize);
+                }
+                else
+                {
+                    const int start_bx = blockSize + i - 1;
+                    const int start_by = j - 1;
+                    const int ax = halfSearchWindowSize + start_bx;
+                    const int ay = halfSearchWindowSize + start_by;
+
+                    const T a_up = srcExtended_.at<T>(ay, ax);
+                    const T a_down = srcExtended_.at<T>(ay + blockSize, ax);
+
+                    for (TT y = 0; y < searchWindowSize; y++)
+                    {
+                        int *distSumsRow = distSums.row_ptr(y);
+                        int *colDistSumsRow = colDistSums.row_ptr(firstColNum, y);
+                        int *lastColDistSumsRow = lastColDistSums.row_ptr(i, y);
+
+                        const T *b_up_ptr = srcExtended_.ptr<T>(start_by + y);
+                        const T *b_down_ptr = srcExtended_.ptr<T>(start_by + y + blockSize);
+
+                        for (TT x = 0; x < searchWindowSize; x++)
+                        {
+                            // Remove from current pixel sum column sum with index "firstColNum"
+                            distSumsRow[x] -= colDistSumsRow[x];
+
+                            const int bx = start_bx + x;
+                            colDistSumsRow[x] = lastColDistSumsRow[x] +
+                                D::template calcUpDownDist<T>(a_up, a_down, b_up_ptr[bx], b_down_ptr[bx]);
+
+                            distSumsRow[x] += colDistSumsRow[x];
+                            lastColDistSumsRow[x] = colDistSumsRow[x];
+
+                            if (x == halfSearchWindowSize && y == halfSearchWindowSize)
+                                continue;
+
+                            // Save the distance, coordinate and increase the counter
+                            if (distSumsRow[x] < hBM)
+                                bm[elementSize++](distSumsRow[x], x, y);
+                        }
+                    }
+                }
+
+                firstColNum = (firstColNum + 1) % blockSize;
+            }
+
+            // Sort bm by distance (first element is already sorted)
+            std::sort(bm + 1, bm + elementSize);
+
+            // Find the nearest power of 2 and cap the group size from the top
+            elementSize = getLargestPowerOf2SmallerThan(elementSize);
+            if (elementSize > groupSize)
+                elementSize = groupSize;
+
+            // Transform 2D patches
+            for (int n = 0; n < elementSize; ++n)
+            {
+                const T *candidatePatch = currentPixel + step * bm[n].coord_y + bm[n].coord_x;
+                TC::forwardTransform2D(candidatePatch, bm[n].data(), step, blockSize);
+            }
+
+            // Transform and shrink 1D columns
+            TT sumNonZero = 0;
+            TT *thrMapPtr1D = thrMap_ + (elementSize - 1) * blockSizeSq;
+            switch (elementSize)
+            {
+            case 16:
+                for (int n = 0; n < blockSizeSq; n++)
+                {
+                    TC::forwardTransform16(bm, n);
+                    sumNonZero += HardThreshold<16>(bm, n, thrMapPtr1D);
+                    TC::inverseTransform16(bm, n);
+                }
+                break;
+            case 8:
+                for (int n = 0; n < blockSizeSq; n++)
+                {
+                    TC::forwardTransform8(bm, n);
+                    sumNonZero += HardThreshold<8>(bm, n, thrMapPtr1D);
+                    TC::inverseTransform8(bm, n);
+                }
+                break;
+            case 4:
+                for (int n = 0; n < blockSizeSq; n++)
+                {
+                    TC::forwardTransform4(bm, n);
+                    sumNonZero += HardThreshold<4>(bm, n, thrMapPtr1D);
+                    TC::inverseTransform4(bm, n);
+                }
+                break;
+            case 2:
+                for (int n = 0; n < blockSizeSq; n++)
+                {
+                    TC::forwardTransform2(bm, n);
+                    TC::forwardTransform2(bm, n);
+                    sumNonZero += HardThreshold<2>(bm, n, thrMapPtr1D);
+                    TC::inverseTransform2(bm, n);
+                }
+                break;
+            case 1:
+                {
+                    TT *block = bm[0].data();
+                    for (int n = 0; n < blockSizeSq; n++)
+                        shrink(block[n], sumNonZero, *thrMapPtr1D++);
+                }
+                break;
+            default:
+                for (int n = 0; n < blockSizeSq; n++)
+                {
+                    TC::forwardTransformN(bm, n, elementSize);
+                    sumNonZero += HardThreshold(bm, n, thrMapPtr1D, elementSize);
+                    TC::inverseTransformN(bm, n, elementSize);
+                }
+            }
+
+            // Inverse 2D transform
+            for (int n = 0; n < elementSize; ++n)
+                TC::inverseTransform2D(bm[n].data(), blockSize);
+
+            // Aggregate the results (increase sumNonZero to avoid division by zero)
+            float weight = 1.0f / (float)(++sumNonZero);
+
+            // Scale weight by element size
+            weight *= elementSize;
+            weight /= groupSize;
+
+            // Put patches back to their original positions
+            WT *dstPtr = weightedSum.data() + jj * dstStep + i;
+            WT *weiPtr = weights.data() + jj * dstStep + i;
+            const float *kaiser = kaiser_;
+
+            for (int l = 0; l < elementSize; ++l)
+            {
+                const TT *block = bm[l].data();
+                int offset = bm[l].coord_y * dstStep + bm[l].coord_x;
+                WT *d = dstPtr + offset;
+                WT *dw = weiPtr + offset;
+
+                for (int n = 0; n < blockSize; ++n)
+                {
+                    for (int m = 0; m < blockSize; ++m)
+                    {
+                        unsigned idx = n * blockSize + m;
+                        *d += kaiser[idx] * block[idx] * weight;
+                        *dw += kaiser[idx] * weight;
+                        ++d, ++dw;
+                    }
+                    d += dstcstep;
+                    dw += weicstep;
+                }
+            }
+        } // i
+    } // j
+
+    // Cleanup
+    for (int i = 0; i < searchWindowSizeSq; ++i)
+        bm[i].release();
+    delete[] bm;
+
+    // Divide accumulation buffer by the corresponding weights
+    for (int i = row_from, ii = 0; i <= row_to; ++i, ++ii)
+    {
+        T *d = dst_.ptr<T>(i);
+        float *dE = weightedSum.data() + (ii + halfSearchWindowSize + halfBlockSize) * dstStep + halfSearchWindowSize;
+        float *dw = weights.data() + (ii + halfSearchWindowSize + halfBlockSize) * dstStep + halfSearchWindowSize;
+        for (int j = 0; j < dst_.cols; ++j)
+            d[j] = cv::saturate_cast<T>(dE[j + halfBlockSize] / dw[j + halfBlockSize]);
+    }
+}
+
+template <typename T, typename D, typename WT, typename TT, typename TC>
+inline void Bm3dDenoisingInvokerStep1<T, D, WT, TT, TC>::calcDistSumsForFirstElementInRow(
+    int i,
+    Array2d<int>& distSums,
+    Array3d<int>& colDistSums,
+    Array3d<int>& lastColDistSums,
+    BlockMatch<TT, int, TT> *bm,
+    int &elementSize) const
+{
+    int j = 0;
+    const int hBM = hBM_;
+    const int blockSize = templateWindowSize_;
+    const int searchWindowSize = searchWindowSize_;
+    const TT halfSearchWindowSize = (TT)halfSearchWindowSize_;
+    const int ay = halfSearchWindowSize + i;
+    const int ax = halfSearchWindowSize + j;
+
+    for (TT y = 0; y < searchWindowSize; ++y)
+    {
+        for (TT x = 0; x < searchWindowSize; ++x)
+        {
+            // Zeroize arrays
+            distSums[y][x] = 0;
+            for (int tx = 0; tx < blockSize; tx++)
+                colDistSums[tx][y][x] = 0;
+
+            int start_y = i + y;
+            int start_x = j + x;
+
+            for (int ty = 0; ty < blockSize; ty++)
+                for (int tx = 0; tx < blockSize; tx++)
+                {
+                    int dist = D::template calcDist<T>(
+                        srcExtended_,
+                        ay + ty,
+                        ax + tx,
+                        start_y + ty,
+                        start_x + tx);
+
+                    distSums[y][x] += dist;
+                    colDistSums[tx][y][x] += dist;
+                }
+
+            lastColDistSums[j][y][x] = colDistSums[blockSize - 1][y][x];
+
+            if (x == halfSearchWindowSize && y == halfSearchWindowSize)
+                continue;
+
+            if (distSums[y][x] < hBM)
+                bm[elementSize++](distSums[y][x], x, y);
+        }
+    }
+}
+
+template <typename T, typename D, typename WT, typename TT, typename TC>
+inline void Bm3dDenoisingInvokerStep1<T, D, WT, TT, TC>::calcDistSumsForAllElementsInFirstRow(
+    int i,
+    int j,
+    int firstColNum,
+    Array2d<int>& distSums,
+    Array3d<int>& colDistSums,
+    Array3d<int>& lastColDistSums,
+    BlockMatch<TT, int, TT> *bm,
+    int &elementSize) const
+{
+    const int hBM = hBM_;
+    const int blockSize = templateWindowSize_;
+    const int searchWindowSize = searchWindowSize_;
+    const TT halfSearchWindowSize = (TT)halfSearchWindowSize_;
+
+    const int bx_start = blockSize - 1 + j;
+    const int ax = halfSearchWindowSize + bx_start;
+    const int ay = halfSearchWindowSize + i;
+
+    for (TT y = 0; y < searchWindowSize; ++y)
+    {
+        for (TT x = 0; x < searchWindowSize; ++x)
+        {
+            distSums[y][x] -= colDistSums[firstColNum][y][x];
+
+            colDistSums[firstColNum][y][x] = 0;
+            int by = i + y;
+            int bx = bx_start + x;
+
+            for (int ty = 0; ty < blockSize; ty++)
+                colDistSums[firstColNum][y][x] += D::template calcDist<T>(
+                    srcExtended_,
+                    ay + ty,
+                    ax,
+                    by + ty,
+                    bx);
+
+            distSums[y][x] += colDistSums[firstColNum][y][x];
+            lastColDistSums[j][y][x] = colDistSums[firstColNum][y][x];
+
+            if (x == halfSearchWindowSize && y == halfSearchWindowSize)
+                continue;
+
+            if (distSums[y][x] < hBM)
+                bm[elementSize++](distSums[y][x], x, y);
+        }
+    }
+}
+
+}  // namespace xphoto
+}  // namespace cv
+
+#endif
\ No newline at end of file

modules/xphoto/src/bm3d_denoising_invoker_step2.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective icvers.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_BM3D_DENOISING_INVOKER_STEP2_HPP__
+#define __OPENCV_BM3D_DENOISING_INVOKER_STEP2_HPP__
+
+#include "bm3d_denoising_invoker_commons.hpp"
+#include "bm3d_denoising_transforms.hpp"
+#include "kaiser_window.hpp"
+
+namespace cv
+{
+namespace xphoto
+{
+
+template <typename T, typename D, typename WT, typename TT, typename TC>
+struct Bm3dDenoisingInvokerStep2 : public ParallelLoopBody
+{
+public:
+    Bm3dDenoisingInvokerStep2(
+        const Mat& src,
+        const Mat& basic,
+        Mat& dst,
+        const int &templateWindowSize,
+        const int &searchWindowSize,
+        const float &h,
+        const int &hBM,
+        const int &groupSize,
+        const int &slidingStep,
+        const float &beta);
+
+    virtual ~Bm3dDenoisingInvokerStep2();
+    void operator() (const Range& range) const;
+
+private:
+    // Unimplemented operator in order to satisfy compiler warning.
+    void operator= (const Bm3dDenoisingInvokerStep2&);
+
+    void calcDistSumsForFirstElementInRow(
+        int i,
+        Array2d<int>& distSums,
+        Array3d<int>& colDistSums,
+        Array3d<int>& lastColDistSums,
+        BlockMatch<TT, int, TT> *bm,
+        int &elementSize) const;
+
+    void calcDistSumsForAllElementsInFirstRow(
+        int i,
+        int j,
+        int firstColNum,
+        Array2d<int>& distSums,
+        Array3d<int>& colDistSums,
+        Array3d<int>& lastColDistSums,
+        BlockMatch<TT, int, TT> *bm,
+        int &elementSize) const;
+
+    // Image containers
+    const Mat& src_;
+    const Mat& basic_;
+    Mat& dst_;
+    Mat srcExtended_;
+    Mat basicExtended_;
+
+    // Border size of the extended src and basic images
+    int borderSize_;
+
+    // Template and window size
+    int templateWindowSize_;
+    int searchWindowSize_;
+
+    // Half template and window size
+    int halfTemplateWindowSize_;
+    int halfSearchWindowSize_;
+
+    // Squared template and window size
+    int templateWindowSizeSq_;
+    int searchWindowSizeSq_;
+
+    // Block matching threshold
+    int hBM_;
+
+    // Maximum size of 3D group
+    int groupSize_;
+
+    // Sliding step
+    const int slidingStep_;
+
+    // Threshold map
+    TT *thrMap_;
+
+    // Kaiser window
+    float *kaiser_;
+};
+
+template <typename T, typename D, typename WT, typename TT, typename TC>
+Bm3dDenoisingInvokerStep2<T, D, WT, TT, TC>::Bm3dDenoisingInvokerStep2(
+    const Mat& src,
+    const Mat& basic,
+    Mat& dst,
+    const int &templateWindowSize,
+    const int &searchWindowSize,
+    const float &h,
+    const int &hBM,
+    const int &groupSize,
+    const int &slidingStep,
+    const float &beta) :
+    src_(src), basic_(basic), dst_(dst), groupSize_(groupSize), slidingStep_(slidingStep), thrMap_(NULL), kaiser_(NULL)
+{
+    groupSize_ = getLargestPowerOf2SmallerThan(groupSize);
+    CV_Assert(groupSize > 0);
+
+    halfTemplateWindowSize_ = templateWindowSize >> 1;
+    halfSearchWindowSize_ = searchWindowSize >> 1;
+    templateWindowSize_ = templateWindowSize;
+    searchWindowSize_ = searchWindowSize;
+    templateWindowSizeSq_ = templateWindowSize_ * templateWindowSize_;
+    searchWindowSizeSq_ = searchWindowSize_ * searchWindowSize_;
+
+    // Extend image to avoid border problem
+    borderSize_ = halfSearchWindowSize_ + halfTemplateWindowSize_;
+    copyMakeBorder(src_, srcExtended_, borderSize_, borderSize_, borderSize_, borderSize_, BORDER_DEFAULT);
+    copyMakeBorder(basic_, basicExtended_, borderSize_, borderSize_, borderSize_, borderSize_, BORDER_DEFAULT);
+
+    // Calculate block matching threshold
+    hBM_ = D::template calcBlockMatchingThreshold<int>(hBM, templateWindowSizeSq_);
+
+    // Select transforms depending on the template size
+    TC::RegisterTransforms2D(templateWindowSize_);
+
+    // Precompute threshold map
+    TC::calcThresholdMap3D(thrMap_, h, templateWindowSize_, groupSize_);
+
+    // Generate kaiser window
+    calcKaiserWindow2D(kaiser_, templateWindowSize_, beta);
+}
+
+template<typename T, typename D, typename WT, typename TT, typename TC>
+inline Bm3dDenoisingInvokerStep2<T, D, WT, TT, TC>::~Bm3dDenoisingInvokerStep2()
+{
+    delete[] thrMap_;
+    delete[] kaiser_;
+}
+
+template <typename T, typename D, typename WT, typename TT, typename TC>
+void Bm3dDenoisingInvokerStep2<T, D, WT, TT, TC>::operator() (const Range& range) const
+{
+    const int size = (range.size() + 2 * borderSize_) * srcExtended_.cols;
+    std::vector<WT> weightedSum(size, 0.0);
+    std::vector<WT> weights(size, 0.0);
+    int row_from = range.start;
+    int row_to = range.end - 1;
+
+    // Local vars for faster processing
+    const int blockSize = templateWindowSize_;
+    const int blockSizeSq = templateWindowSizeSq_;
+    const int halfBlockSize = halfTemplateWindowSize_;
+    const int searchWindowSize = searchWindowSize_;
+    const int searchWindowSizeSq = searchWindowSizeSq_;
+    const TT halfSearchWindowSize = (TT)halfSearchWindowSize_;
+    const int hBM = hBM_;
+    const int groupSize = groupSize_;
+
+    const int step = srcExtended_.cols;
+    const int dstStep = srcExtended_.cols;
+    const int weiStep = srcExtended_.cols;
+    const int dstcstep = dstStep - blockSize;
+    const int weicstep = weiStep - blockSize;
+
+    // Buffer to store 3D group
+    BlockMatch<TT, int, TT> *bmBasic = new BlockMatch<TT, int, TT>[searchWindowSizeSq];
+    BlockMatch<TT, int, TT> *bmSrc = new BlockMatch<TT, int, TT>[searchWindowSizeSq];
+    for (int i = 0; i < searchWindowSizeSq; ++i)
+    {
+        bmBasic[i].init(blockSizeSq);
+        bmSrc[i].init(blockSizeSq);
+    }
+
+    // First element in a group is always the reference patch. Hence distance is 0.
+    bmBasic[0](0, halfSearchWindowSize, halfSearchWindowSize);
+    bmSrc[0](0, halfSearchWindowSize, halfSearchWindowSize);
+
+    // Sums of columns and rows for current pixel
+    Array2d<int> distSums(searchWindowSize, searchWindowSize);
+
+    // Sums of columns for current pixel (for lazy calc optimization)
+    Array3d<int> colDistSums(blockSize, searchWindowSize, searchWindowSize);
+
+    // Last elements of column sum (for each element in a row)
+    Array3d<int> lastColDistSums(src_.cols, searchWindowSize, searchWindowSize);
+
+    int firstColNum = -1;
+    for (int j = row_from, jj = 0; j <= row_to; j += slidingStep_, jj += slidingStep_)
+    {
+        for (int i = 0; i < src_.cols; i += slidingStep_)
+        {
+            const T *currentPixelSrc = srcExtended_.ptr<T>(0) + step*j + i;
+            const T *currentPixelBasic = basicExtended_.ptr<T>(0) + step*j + i;
+
+            int elementSize = 1;
+
+            // Calculate distSums using moving average filter approach.
+            if (i == 0)
+            {
+                // Calculate distSums for the first element in a row
+                calcDistSumsForFirstElementInRow(j, distSums, colDistSums, lastColDistSums, bmBasic, elementSize);
+                firstColNum = 0;
+            }
+            else
+            {
+                if (j == row_from)
+                {
+                    // Calculate distSums for all elements in the first row
+                    calcDistSumsForAllElementsInFirstRow(
+                        j, i, firstColNum, distSums, colDistSums, lastColDistSums, bmBasic, elementSize);
+                }
+                else
+                {
+                    const int start_bx = blockSize + i - 1;
+                    const int start_by = j - 1;
+                    const int ax = halfSearchWindowSize + start_bx;
+                    const int ay = halfSearchWindowSize + start_by;
+
+                    const T a_up = basicExtended_.at<T>(ay, ax);
+                    const T a_down = basicExtended_.at<T>(ay + blockSize, ax);
+
+                    for (TT y = 0; y < searchWindowSize; y++)
+                    {
+                        int *distSumsRow = distSums.row_ptr(y);
+                        int *colDistSumsRow = colDistSums.row_ptr(firstColNum, y);
+                        int *lastColDistSumsRow = lastColDistSums.row_ptr(i, y);
+
+                        const T *b_up_ptr = basicExtended_.ptr<T>(start_by + y);
+                        const T *b_down_ptr = basicExtended_.ptr<T>(start_by + y + blockSize);
+
+                        for (TT x = 0; x < searchWindowSize; x++)
+                        {
+                            // Remove from current pixel sum column sum with index "firstColNum"
+                            distSumsRow[x] -= colDistSumsRow[x];
+
+                            const int bx = start_bx + x;
+                            colDistSumsRow[x] = lastColDistSumsRow[x] +
+                                D::template calcUpDownDist<T>(a_up, a_down, b_up_ptr[bx], b_down_ptr[bx]);
+
+                            distSumsRow[x] += colDistSumsRow[x];
+                            lastColDistSumsRow[x] = colDistSumsRow[x];
+
+                            if (x == halfSearchWindowSize && y == halfSearchWindowSize)
+                                continue;
+
+                            // Save the distance, coordinate and increase the counter
+                            if (distSumsRow[x] < hBM)
+                                bmBasic[elementSize++](distSumsRow[x], x, y);
+                        }
+                    }
+                }
+
+                firstColNum = (firstColNum + 1) % blockSize;
+            }
+
+            // Sort bmBasic by distance (first element is already sorted)
+            std::sort(bmBasic + 1, bmBasic + elementSize);
+
+            // Find the nearest power of 2 and cap the group size from the top
+            elementSize = getLargestPowerOf2SmallerThan(elementSize);
+            if (elementSize > groupSize)
+                elementSize = groupSize;
+
+            // Transform 2D patches
+            for (int n = 0; n < elementSize; ++n)
+            {
+                const T *candidatePatchSrc = currentPixelSrc + step * bmBasic[n].coord_y + bmBasic[n].coord_x;
+                const T *candidatePatchBasic = currentPixelBasic + step * bmBasic[n].coord_y + bmBasic[n].coord_x;
+                TC::forwardTransform2D(candidatePatchSrc, bmSrc[n].data(), step, blockSize);
+                TC::forwardTransform2D(candidatePatchBasic, bmBasic[n].data(), step, blockSize);
+            }
+
+            // Transform and shrink 1D columns
+            int wienerCoefficients = 0;
+            TT *thrMapPtr1D = thrMap_ + (elementSize - 1) * blockSizeSq;
+            switch (elementSize)
+            {
+            case 16:
+                for (int n = 0; n < blockSizeSq; n++)
+                {
+                    TC::forwardTransform16(bmSrc, n);
+                    TC::forwardTransform16(bmBasic, n);
+                    wienerCoefficients += WienerFiltering<16>(bmSrc, bmBasic, n, thrMapPtr1D);
+                    TC::inverseTransform16(bmBasic, n);
+                }
+                break;
+            case 8:
+                for (int n = 0; n < blockSizeSq; n++)
+                {
+                    TC::forwardTransform8(bmSrc, n);
+                    TC::forwardTransform8(bmBasic, n);
+                    wienerCoefficients += WienerFiltering<8>(bmSrc, bmBasic, n, thrMapPtr1D);
+                    TC::inverseTransform8(bmBasic, n);
+                }
+                break;
+            case 4:
+                for (int n = 0; n < blockSizeSq; n++)
+                {
+                    TC::forwardTransform4(bmSrc, n);
+                    TC::forwardTransform4(bmBasic, n);
+                    wienerCoefficients += WienerFiltering<4>(bmSrc, bmBasic, n, thrMapPtr1D);
+                    TC::inverseTransform4(bmBasic, n);
+                }
+                break;
+            case 2:
+                for (int n = 0; n < blockSizeSq; n++)
+                {
+                    TC::forwardTransform2(bmSrc, n);
+                    TC::forwardTransform2(bmBasic, n);
+                    wienerCoefficients += WienerFiltering<2>(bmSrc, bmBasic, n, thrMapPtr1D);
+                    TC::inverseTransform2(bmBasic, n);
+                }
+                break;
+            case 1:
+            {
+                for (int n = 0; n < blockSizeSq; n++)
+                    wienerCoefficients += WienerFiltering<1>(bmSrc, bmBasic, n, thrMapPtr1D);
+            }
+            break;
+            default:
+                for (int n = 0; n < blockSizeSq; n++)
+                {
+                    TC::forwardTransformN(bmSrc, n, elementSize);
+                    TC::forwardTransformN(bmBasic, n, elementSize);
+                    wienerCoefficients += WienerFiltering(bmSrc, bmBasic, n, thrMapPtr1D, elementSize);
+                    TC::inverseTransformN(bmBasic, n, elementSize);
+                }
+            }
+
+            // Inverse 2D transform
+            for (int n = 0; n < elementSize; ++n)
+                TC::inverseTransform2D(bmBasic[n].data(), blockSize);
+
+            // Aggregate the results (increase sumNonZero to avoid division by zero)
+            float weight = 1.0f / (float)(++wienerCoefficients);
+
+            // Scale weight by element size
+            weight *= elementSize;
+            weight /= groupSize;
+
+            // Put patches back to their original positions
+            WT *dstPtr = weightedSum.data() + jj * dstStep + i;
+            WT *weiPtr = weights.data() + jj * dstStep + i;
+            const float *kaiser = kaiser_;
+
+            for (int l = 0; l < elementSize; ++l)
+            {
+                const TT *block = bmBasic[l].data();
+                int offset = bmBasic[l].coord_y * dstStep + bmBasic[l].coord_x;
+                WT *d = dstPtr + offset;
+                WT *dw = weiPtr + offset;
+
+                for (int n = 0; n < blockSize; ++n)
+                {
+                    for (int m = 0; m < blockSize; ++m)
+                    {
+                        unsigned idx = n * blockSize + m;
+                        *d += kaiser[idx] * block[idx] * weight;
+                        *dw += kaiser[idx] * weight;
+                        ++d, ++dw;
+                    }
+                    d += dstcstep;
+                    dw += weicstep;
+                }
+            }
+        } // i
+    } // j
+
+    // Cleanup
+    for (int i = 0; i < searchWindowSizeSq; ++i)
+    {
+        bmBasic[i].release();
+        bmSrc[i].release();
+    }
+
+    delete[] bmSrc;
+    delete[] bmBasic;
+
+    // Divide accumulation buffer by the corresponding weights
+    for (int i = row_from, ii = 0; i <= row_to; ++i, ++ii)
+    {
+        T *d = dst_.ptr<T>(i);
+        float *dE = weightedSum.data() + (ii + halfSearchWindowSize + halfBlockSize) * dstStep + halfSearchWindowSize;
+        float *dw = weights.data() + (ii + halfSearchWindowSize + halfBlockSize) * dstStep + halfSearchWindowSize;
+        for (int j = 0; j < dst_.cols; ++j)
+            d[j] = cv::saturate_cast<T>(dE[j + halfBlockSize] / dw[j + halfBlockSize]);
+    }
+}
+
+
+template <typename T, typename D, typename WT, typename TT, typename TC>
+inline void Bm3dDenoisingInvokerStep2<T, D, WT, TT, TC>::calcDistSumsForFirstElementInRow(
+    int i,
+    Array2d<int>& distSums,
+    Array3d<int>& colDistSums,
+    Array3d<int>& lastColDistSums,
+    BlockMatch<TT, int, TT> *bm,
+    int &elementSize) const
+{
+    int j = 0;
+    const int hBM = hBM_;
+    const int blockSize = templateWindowSize_;
+    const int searchWindowSize = searchWindowSize_;
+    const TT halfSearchWindowSize = (TT)halfSearchWindowSize_;
+    const int ay = halfSearchWindowSize + i;
+    const int ax = halfSearchWindowSize + j;
+
+    for (TT y = 0; y < searchWindowSize; ++y)
+    {
+        for (TT x = 0; x < searchWindowSize; ++x)
+        {
+            // Zeroize arrays
+            distSums[y][x] = 0;
+            for (int tx = 0; tx < blockSize; tx++)
+                colDistSums[tx][y][x] = 0;
+
+            int start_y = i + y;
+            int start_x = j + x;
+
+            for (int ty = 0; ty < blockSize; ty++)
+                for (int tx = 0; tx < blockSize; tx++)
+                {
+                    int dist = D::template calcDist<T>(
+                        basicExtended_,
+                        ay + ty,
+                        ax + tx,
+                        start_y + ty,
+                        start_x + tx);
+
+                    distSums[y][x] += dist;
+                    colDistSums[tx][y][x] += dist;
+                }
+
+            lastColDistSums[j][y][x] = colDistSums[blockSize - 1][y][x];
+
+            if (x == halfSearchWindowSize && y == halfSearchWindowSize)
+                continue;
+
+            if (distSums[y][x] < hBM)
+                bm[elementSize++](distSums[y][x], x, y);
+        }
+    }
+}
+
+template <typename T, typename D, typename WT, typename TT, typename TC>
+inline void Bm3dDenoisingInvokerStep2<T, D, WT, TT, TC>::calcDistSumsForAllElementsInFirstRow(
+    int i,
+    int j,
+    int firstColNum,
+    Array2d<int>& distSums,
+    Array3d<int>& colDistSums,
+    Array3d<int>& lastColDistSums,
+    BlockMatch<TT, int, TT> *bm,
+    int &elementSize) const
+{
+    const int hBM = hBM_;
+    const int blockSize = templateWindowSize_;
+    const int searchWindowSize = searchWindowSize_;
+    const TT halfSearchWindowSize = (TT)halfSearchWindowSize_;
+
+    const int bx_start = blockSize - 1 + j;
+    const int ax = halfSearchWindowSize + bx_start;
+    const int ay = halfSearchWindowSize + i;
+
+    for (TT y = 0; y < searchWindowSize; ++y)
+    {
+        for (TT x = 0; x < searchWindowSize; ++x)
+        {
+            distSums[y][x] -= colDistSums[firstColNum][y][x];
+
+            colDistSums[firstColNum][y][x] = 0;
+            int by = i + y;
+            int bx = bx_start + x;
+
+            for (int ty = 0; ty < blockSize; ty++)
+                colDistSums[firstColNum][y][x] += D::template calcDist<T>(
+                    basicExtended_,
+                    ay + ty,
+                    ax,
+                    by + ty,
+                    bx);
+
+            distSums[y][x] += colDistSums[firstColNum][y][x];
+            lastColDistSums[j][y][x] = colDistSums[firstColNum][y][x];
+
+            if (x == halfSearchWindowSize && y == halfSearchWindowSize)
+                continue;
+
+            if (distSums[y][x] < hBM)
+                bm[elementSize++](distSums[y][x], x, y);
+        }
+    }
+}
+
+}  // namespace xphoto
+}  // namespace cv
+
+#endif
\ No newline at end of file

modules/xphoto/src/bm3d_denoising_invoker_structs.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective icvers.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_BM3D_DENOISING_INVOKER_STRUCTS_HPP__
+#define __OPENCV_BM3D_DENOISING_INVOKER_STRUCTS_HPP__
+
+namespace cv
+{
+namespace xphoto
+{
+
+template <typename T, typename DT, typename CT>
+class BlockMatch
+{
+public:
+    // Data accessor
+    T* data()
+    {
+        return data_;
+    }
+
+    // Const version of data accessor
+    const T* data() const
+    {
+        return data_;
+    }
+
+    // Allocate memory for data
+    void init(const int &blockSizeSq)
+    {
+        data_ = new T[blockSizeSq];
+    }
+
+    // Release data memory
+    void release()
+    {
+        delete[] data_;
+    }
+
+    // Overloaded operator for convenient assignment
+    void operator()(const DT &_dist, const CT &_coord_x, const CT &_coord_y)
+    {
+        dist = _dist;
+        coord_x = _coord_x;
+        coord_y = _coord_y;
+    }
+
+    // Overloaded array subscript operator
+    T& operator[](const std::size_t &idx)
+    {
+        return data_[idx];
+    };
+
+    // Overloaded const array subscript operator
+    const T& operator[](const std::size_t &idx) const
+    {
+        return data_[idx];
+    };
+
+    // Overloaded comparison operator for sorting
+    bool operator<(const BlockMatch& right) const
+    {
+        return dist < right.dist;
+    }
+
+    // Block matching distance
+    DT dist;
+
+    // Relative coordinates to the current search window
+    CT coord_x;
+    CT coord_y;
+
+private:
+    // Pointer to the pixel values of the block
+    T *data_;
+};
+
+class DistAbs
+{
+    template <typename T>
+    struct calcDist_
+    {
+        static inline int f(const T &a, const T &b)
+        {
+            return std::abs(a - b);
+        }
+    };
+
+    template <typename ET>
+    struct calcDist_<Vec<ET, 2> >
+    {
+        static inline int f(const Vec<ET, 2> a, const Vec<ET, 2> b)
+        {
+            return std::abs((int)(a[0] - b[0])) + std::abs((int)(a[1] - b[1]));
+        }
+    };
+
+    template <typename ET>
+    struct calcDist_<Vec<ET, 3> >
+    {
+        static inline int f(const Vec<ET, 3> a, const Vec<ET, 3> b)
+        {
+            return
+                std::abs((int)(a[0] - b[0])) +
+                std::abs((int)(a[1] - b[1])) +
+                std::abs((int)(a[2] - b[2]));
+        }
+    };
+
+    template <typename ET>
+    struct calcDist_<Vec<ET, 4> >
+    {
+        static inline int f(const Vec<ET, 4> a, const Vec<ET, 4> b)
+        {
+            return
+                std::abs((int)(a[0] - b[0])) +
+                std::abs((int)(a[1] - b[1])) +
+                std::abs((int)(a[2] - b[2])) +
+                std::abs((int)(a[3] - b[3]));
+        }
+    };
+
+public:
+    template <typename T>
+    static inline int calcDist(const T &a, const T &b)
+    {
+        return calcDist_<T>::f(a, b);
+    }
+
+    template <typename T>
+    static inline int calcDist(const Mat& m, int i1, int j1, int i2, int j2)
+    {
+        const T a = m.at<T>(i1, j1);
+        const T b = m.at<T>(i2, j2);
+        return calcDist<T>(a, b);
+    }
+
+    template <typename T>
+    static inline int calcUpDownDist(T a_up, T a_down, T b_up, T b_down)
+    {
+        return calcDist<T>(a_down, b_down) - calcDist<T>(a_up, b_up);
+    };
+
+    template <typename T>
+    static inline T calcBlockMatchingThreshold(const T &blockMatchThrL2, const T &blockSizeSq)
+    {
+        return (T)(std::sqrt((double)blockMatchThrL2) * blockSizeSq);
+    }
+};
+
+class DistSquared
+{
+    template <typename T>
+    struct calcDist_
+    {
+        static inline int f(const T &a, const T &b)
+        {
+            return (a - b) * (a - b);
+        }
+    };
+
+    template <typename ET>
+    struct calcDist_<Vec<ET, 2> >
+    {
+        static inline int f(const Vec<ET, 2> a, const Vec<ET, 2> b)
+        {
+            return (int)(a[0] - b[0])*(int)(a[0] - b[0]) + (int)(a[1] - b[1])*(int)(a[1] - b[1]);
+        }
+    };
+
+    template <typename ET>
+    struct calcDist_<Vec<ET, 3> >
+    {
+        static inline int f(const Vec<ET, 3> a, const Vec<ET, 3> b)
+        {
+            return
+                (int)(a[0] - b[0])*(int)(a[0] - b[0]) +
+                (int)(a[1] - b[1])*(int)(a[1] - b[1]) +
+                (int)(a[2] - b[2])*(int)(a[2] - b[2]);
+        }
+    };
+
+    template <typename ET>
+    struct calcDist_<Vec<ET, 4> >
+    {
+        static inline int f(const Vec<ET, 4> a, const Vec<ET, 4> b)
+        {
+            return
+                (int)(a[0] - b[0])*(int)(a[0] - b[0]) +
+                (int)(a[1] - b[1])*(int)(a[1] - b[1]) +
+                (int)(a[2] - b[2])*(int)(a[2] - b[2]) +
+                (int)(a[3] - b[3])*(int)(a[3] - b[3]);
+        }
+    };
+
+    template <typename T> struct calcUpDownDist_
+    {
+        static inline int f(T a_up, T a_down, T b_up, T b_down)
+        {
+            int A = a_down - b_down;
+            int B = a_up - b_up;
+            return (A - B)*(A + B);
+        }
+    };
+
+    template <typename ET, int n> struct calcUpDownDist_<Vec<ET, n> >
+    {
+    private:
+        typedef Vec<ET, n> T;
+    public:
+        static inline int f(T a_up, T a_down, T b_up, T b_down)
+        {
+            return calcDist<T>(a_down, b_down) - calcDist<T>(a_up, b_up);
+        }
+    };
+
+public:
+    template <typename T>
+    static inline int calcDist(const T &a, const T &b)
+    {
+        return calcDist_<T>::f(a, b);
+    }
+
+    template <typename T>
+    static inline int calcDist(const Mat& m, int i1, int j1, int i2, int j2)
+    {
+        const T a = m.at<T>(i1, j1);
+        const T b = m.at<T>(i2, j2);
+        return calcDist<T>(a, b);
+    }
+
+    template <typename T>
+    static inline int calcUpDownDist(T a_up, T a_down, T b_up, T b_down)
+    {
+        return calcUpDownDist_<T>::f(a_up, a_down, b_up, b_down);
+    };
+
+    template <typename T>
+    static inline T calcBlockMatchingThreshold(const T &blockMatchThrL2, const T &blockSizeSq)
+    {
+        return blockMatchThrL2 * blockSizeSq;
+    }
+};
+
+template <class T>
+struct Array2d
+{
+    T* a;
+    int n1, n2;
+    bool needToDeallocArray;
+
+    Array2d(const Array2d& array2d) :
+        a(array2d.a), n1(array2d.n1), n2(array2d.n2), needToDeallocArray(false)
+    {
+        if (array2d.needToDeallocArray)
+        {
+            CV_Error(Error::BadDataPtr, "Copy constructor for self allocating arrays not supported");
+        }
+    }
+
+    Array2d(T* _a, int _n1, int _n2) :
+        a(_a), n1(_n1), n2(_n2), needToDeallocArray(false)
+    {
+    }
+
+    Array2d(int _n1, int _n2) :
+        n1(_n1), n2(_n2), needToDeallocArray(true)
+    {
+        a = new T[n1*n2];
+    }
+
+    ~Array2d()
+    {
+        if (needToDeallocArray)
+            delete[] a;
+    }
+
+    T* operator [] (int i)
+    {
+        return a + i*n2;
+    }
+
+    inline T* row_ptr(int i)
+    {
+        return (*this)[i];
+    }
+};
+
+template <class T>
+struct Array3d
+{
+    T* a;
+    int n1, n2, n3;
+    bool needToDeallocArray;
+
+    Array3d(T* _a, int _n1, int _n2, int _n3) :
+        a(_a), n1(_n1), n2(_n2), n3(_n3), needToDeallocArray(false)
+    {
+    }
+
+    Array3d(int _n1, int _n2, int _n3) :
+        n1(_n1), n2(_n2), n3(_n3), needToDeallocArray(true)
+    {
+        a = new T[n1*n2*n3];
+    }
+
+    ~Array3d()
+    {
+        if (needToDeallocArray)
+            delete[] a;
+    }
+
+    Array2d<T> operator [] (int i)
+    {
+        Array2d<T> array2d(a + i*n2*n3, n2, n3);
+        return array2d;
+    }
+
+    inline T* row_ptr(int i1, int i2)
+    {
+        return a + i1*n2*n3 + i2*n3;
+    }
+};
+
+}  // namespace xphoto
+}  // namespace cv
+
+#endif
\ No newline at end of file

modules/xphoto/src/bm3d_denoising_transforms.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective icvers.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_BM3D_DENOISING_TRANSFORMS_HPP__
+#define __OPENCV_BM3D_DENOISING_TRANSFORMS_HPP__
+
+#include "bm3d_denoising_transforms_haar.hpp"
+
+namespace cv
+{
+namespace xphoto
+{
+
+// Following class contains interface of the tranform domain functions.
+template <typename T, typename TT>
+class Transform
+{
+  public:
+    // 2D transforms
+    typedef void(*Forward2D)(const T *ptr, TT *dst, const int &step, const int blockSize);
+    typedef void(*Inverse2D)(TT *src, const int blockSize);
+
+    // 1D transforms
+    typedef void(*Forward1D)(BlockMatch<TT, int, TT> *z, const int &n, const unsigned &N);
+    typedef void(*Inverse1D)(BlockMatch<TT, int, TT> *z, const int &n, const unsigned &N);
+
+    // Specialized 1D transforms
+    typedef void(*Forward1Ds)(BlockMatch<TT, int, TT> *z, const int &n);
+    typedef void(*Inverse1Ds)(BlockMatch<TT, int, TT> *z, const int &n);
+};
+
+}  // namespace xphoto
+}  // namespace cv
+
+#endif
\ No newline at end of file

modules/xphoto/src/bm3d_denoising_transforms_1D.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective icvers.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_BM3D_DENOISING_TRANSFORMS_1D_HPP__
+#define __OPENCV_BM3D_DENOISING_TRANSFORMS_1D_HPP__
+
+namespace cv
+{
+namespace xphoto
+{
+
+class HaarTransform1D
+{
+    static void CalculateIndicesN(unsigned *diffIndices, const unsigned &size, const unsigned &N)
+    {
+        unsigned diffIdx = 1;
+        unsigned diffAllIdx = 0;
+        for (unsigned i = 1; i <= N; i <<= 1)
+        {
+            diffAllIdx += (i >> 1);
+            for (unsigned j = 0; j < (i >> 1); ++j)
+                diffIndices[diffIdx++] = size - (--diffAllIdx);
+            diffAllIdx += i;
+        }
+    }
+
+public:
+    /// 1D forward transformations of array of arbitrary size
+    template <typename T, typename DT, typename CT>
+    inline static void ForwardTransformN(BlockMatch<T, DT, CT> *src, const int &n, const unsigned &N)
+    {
+        const unsigned size = N + (N << 1) - 2;
+        T *dstX = new T[size];
+
+        // Fill dstX with source values
+        for (unsigned i = 0; i < N; ++i)
+            dstX[i] = src[i][n];
+
+        unsigned idx = 0, dstIdx = N;
+        for (unsigned i = N; i > 1; i >>= 1)
+        {
+            // Get sums
+            for (unsigned j = 0; j < (i >> 1); ++j)
+                dstX[dstIdx++] = (dstX[idx + 2 * j] + dstX[idx + j * 2 + 1] + 1) >> 1;
+
+            // Get diffs
+            for (unsigned j = 0; j < (i >> 1); ++j)
+                dstX[dstIdx++] = dstX[idx + 2 * j] - dstX[idx + j * 2 + 1];
+
+            idx = dstIdx - i;
+        }
+
+        // Calculate indices in the destination matrix.
+        unsigned *diffIndices = new unsigned[N];
+        CalculateIndicesN(diffIndices, size, N);
+
+        // Fill in destination matrix
+        src[0][n] = dstX[size - 2];
+        for (unsigned i = 1; i < N; ++i)
+            src[i][n] = dstX[diffIndices[i]];
+
+        delete[] dstX;
+        delete[] diffIndices;
+    }
+
+    /// 1D inverse transformation of array of arbitrary size
+    template <typename T, typename DT, typename CT>
+    inline static void InverseTransformN(BlockMatch<T, DT, CT> *src, const int &n, const unsigned &N)
+    {
+        const unsigned dstSize = (N << 1) - 2;
+        T *dstX = new T[dstSize];
+        T *srcX = new T[N];
+
+        // Fill srcX with source values
+        srcX[0] = src[0][n] * 2;
+        for (unsigned i = 1; i < N; ++i)
+            srcX[i] = src[i][n];
+
+        // Take care of first two elements
+        dstX[0] = srcX[0] + srcX[1];
+        dstX[1] = srcX[0] - srcX[1];
+
+        unsigned idx = 0, dstIdx = 2;
+        for (unsigned i = 4; i < N; i <<= 1)
+        {
+            for (unsigned j = 0; j < (i >> 1); ++j)
+            {
+                dstX[dstIdx++] = dstX[idx + j] + srcX[idx + 2 + j];
+                dstX[dstIdx++] = dstX[idx + j] - srcX[idx + 2 + j];
+            }
+            idx += (i >> 1);
+        }
+
+        // Handle the last X elements
+        dstIdx = 0;
+        for (unsigned j = 0; j < (N >> 1); ++j)
+        {
+            src[dstIdx++][n] = (dstX[idx + j] + srcX[idx + 2 + j]) >> 1;
+            src[dstIdx++][n] = (dstX[idx + j] - srcX[idx + 2 + j]) >> 1;
+        }
+
+        delete[] srcX;
+        delete[] dstX;
+    }
+
+    /// 1D forward transformations of fixed array size: 2, 4, 8 and 16
+
+    template <typename T, typename DT, typename CT>
+    inline static void ForwardTransform2(BlockMatch<T, DT, CT> *z, const int &n)
+    {
+        T sum = (z[0][n] + z[1][n] + 1) >> 1;
+        T dif = z[0][n] - z[1][n];
+
+        z[0][n] = sum;
+        z[1][n] = dif;
+    }
+
+    template <typename T, typename DT, typename CT>
+    inline static void ForwardTransform4(BlockMatch<T, DT, CT> *z, const int &n)
+    {
+        T sum0 = (z[0][n] + z[1][n] + 1) >> 1;
+        T sum1 = (z[2][n] + z[3][n] + 1) >> 1;
+        T dif0 = z[0][n] - z[1][n];
+        T dif1 = z[2][n] - z[3][n];
+
+        T sum00 = (sum0 + sum1 + 1) >> 1;
+        T dif00 = sum0 - sum1;
+
+        z[0][n] = sum00;
+        z[1][n] = dif00;
+        z[2][n] = dif0;
+        z[3][n] = dif1;
+    }
+
+    template <typename T, typename DT, typename CT>
+    inline static void ForwardTransform8(BlockMatch<T, DT, CT> *z, const int &n)
+    {
+        T sum0 = (z[0][n] + z[1][n] + 1) >> 1;
+        T sum1 = (z[2][n] + z[3][n] + 1) >> 1;
+        T sum2 = (z[4][n] + z[5][n] + 1) >> 1;
+        T sum3 = (z[6][n] + z[7][n] + 1) >> 1;
+        T dif0 = z[0][n] - z[1][n];
+        T dif1 = z[2][n] - z[3][n];
+        T dif2 = z[4][n] - z[5][n];
+        T dif3 = z[6][n] - z[7][n];
+
+        T sum00 = (sum0 + sum1 + 1) >> 1;
+        T sum11 = (sum2 + sum3 + 1) >> 1;
+        T dif00 = sum0 - sum1;
+        T dif11 = sum2 - sum3;
+
+        T sum000 = (sum00 + sum11 + 1) >> 1;
+        T dif000 = sum00 - sum11;
+
+        z[0][n] = sum000;
+        z[1][n] = dif000;
+        z[2][n] = dif00;
+        z[3][n] = dif11;
+        z[4][n] = dif0;
+        z[5][n] = dif1;
+        z[6][n] = dif2;
+        z[7][n] = dif3;
+    }
+
+    template <typename T, typename DT, typename CT>
+    inline static void ForwardTransform16(BlockMatch<T, DT, CT> *z, const int &n)
+    {
+        T sum0 = (z[0][n] + z[1][n] + 1) >> 1;
+        T sum1 = (z[2][n] + z[3][n] + 1) >> 1;
+        T sum2 = (z[4][n] + z[5][n] + 1) >> 1;
+        T sum3 = (z[6][n] + z[7][n] + 1) >> 1;
+        T sum4 = (z[8][n] + z[9][n] + 1) >> 1;
+        T sum5 = (z[10][n] + z[11][n] + 1) >> 1;
+        T sum6 = (z[12][n] + z[13][n] + 1) >> 1;
+        T sum7 = (z[14][n] + z[15][n] + 1) >> 1;
+        T dif0 = z[0][n] - z[1][n];
+        T dif1 = z[2][n] - z[3][n];
+        T dif2 = z[4][n] - z[5][n];
+        T dif3 = z[6][n] - z[7][n];
+        T dif4 = z[8][n] - z[9][n];
+        T dif5 = z[10][n] - z[11][n];
+        T dif6 = z[12][n] - z[13][n];
+        T dif7 = z[14][n] - z[15][n];
+
+        T sum00 = (sum0 + sum1 + 1) >> 1;
+        T sum11 = (sum2 + sum3 + 1) >> 1;
+        T sum22 = (sum4 + sum5 + 1) >> 1;
+        T sum33 = (sum6 + sum7 + 1) >> 1;
+        T dif00 = sum0 - sum1;
+        T dif11 = sum2 - sum3;
+        T dif22 = sum4 - sum5;
+        T dif33 = sum6 - sum7;
+
+        T sum000 = (sum00 + sum11 + 1) >> 1;
+        T sum111 = (sum22 + sum33 + 1) >> 1;
+        T dif000 = sum00 - sum11;
+        T dif111 = sum22 - sum33;
+
+        T sum0000 = (sum000 + sum111 + 1) >> 1;
+        T dif0000 = dif000 - dif111;
+
+        z[0][n] = sum0000;
+        z[1][n] = dif0000;
+        z[2][n] = dif000;
+        z[3][n] = dif111;
+        z[4][n] = dif00;
+        z[5][n] = dif11;
+        z[6][n] = dif22;
+        z[7][n] = dif33;
+        z[8][n] = dif0;
+        z[9][n] = dif1;
+        z[10][n] = dif2;
+        z[11][n] = dif3;
+        z[12][n] = dif4;
+        z[13][n] = dif5;
+        z[14][n] = dif6;
+        z[15][n] = dif7;
+    }
+
+    /// 1D inverse transformations of fixed array size: 2, 4, 8 and 16
+
+    template <typename T, typename DT, typename CT>
+    inline static void InverseTransform2(BlockMatch<T, DT, CT> *src, const int &n)
+    {
+        T src0 = src[0][n] * 2;
+        T src1 = src[1][n];
+
+        src[0][n] = (src0 + src1) >> 1;
+        src[1][n] = (src0 - src1) >> 1;
+    }
+
+    template <typename T, typename DT, typename CT>
+    inline static void InverseTransform4(BlockMatch<T, DT, CT> *src, const int &n)
+    {
+        T src0 = src[0][n] * 2;
+        T src1 = src[1][n];
+        T src2 = src[2][n];
+        T src3 = src[3][n];
+
+        T sum0 = src0 + src1;
+        T dif0 = src0 - src1;
+
+        src[0][n] = (sum0 + src2) >> 1;
+        src[1][n] = (sum0 - src2) >> 1;
+        src[2][n] = (dif0 + src3) >> 1;
+        src[3][n] = (dif0 - src3) >> 1;
+    }
+
+    template <typename T, typename DT, typename CT>
+    inline static void InverseTransform8(BlockMatch<T, DT, CT> *src, const int &n)
+    {
+        T src0 = src[0][n] * 2;
+        T src1 = src[1][n];
+        T src2 = src[2][n];
+        T src3 = src[3][n];
+        T src4 = src[4][n];
+        T src5 = src[5][n];
+        T src6 = src[6][n];
+        T src7 = src[7][n];
+
+        T sum0 = src0 + src1;
+        T dif0 = src0 - src1;
+
+        T sum00 = sum0 + src2;
+        T dif00 = sum0 - src2;
+        T sum11 = dif0 + src3;
+        T dif11 = dif0 - src3;
+
+        src[0][n] = (sum00 + src4) >> 1;
+        src[1][n] = (sum00 - src4) >> 1;
+        src[2][n] = (dif00 + src5) >> 1;
+        src[3][n] = (dif00 - src5) >> 1;
+        src[4][n] = (sum11 + src6) >> 1;
+        src[5][n] = (sum11 - src6) >> 1;
+        src[6][n] = (dif11 + src7) >> 1;
+        src[7][n] = (dif11 - src7) >> 1;
+    }
+
+    template <typename T, typename DT, typename CT>
+    inline static void InverseTransform16(BlockMatch<T, DT, CT> *src, const int &n)
+    {
+        T src0 = src[0][n] * 2;
+        T src1 = src[1][n];
+        T src2 = src[2][n];
+        T src3 = src[3][n];
+        T src4 = src[4][n];
+        T src5 = src[5][n];
+        T src6 = src[6][n];
+        T src7 = src[7][n];
+        T src8 = src[8][n];
+        T src9 = src[9][n];
+        T src10 = src[10][n];
+        T src11 = src[11][n];
+        T src12 = src[12][n];
+        T src13 = src[13][n];
+        T src14 = src[14][n];
+        T src15 = src[15][n];
+
+        T sum0 = src0 + src1;
+        T dif0 = src0 - src1;
+
+        T sum00 = sum0 + src2;
+        T dif00 = sum0 - src2;
+        T sum11 = dif0 + src3;
+        T dif11 = dif0 - src3;
+
+        T sum000 = sum00 + src4;
+        T dif000 = sum00 - src4;
+        T sum111 = dif00 + src5;
+        T dif111 = dif00 - src5;
+        T sum222 = sum11 + src6;
+        T dif222 = sum11 - src6;
+        T sum333 = dif11 + src7;
+        T dif333 = dif11 - src7;
+
+        src[0][n] = (sum000 + src8) >> 1;
+        src[1][n] = (sum000 - src8) >> 1;
+        src[2][n] = (dif000 + src9) >> 1;
+        src[3][n] = (dif000 - src9) >> 1;
+        src[4][n] = (sum111 + src10) >> 1;
+        src[5][n] = (sum111 - src10) >> 1;
+        src[6][n] = (dif111 + src11) >> 1;
+        src[7][n] = (dif111 - src11) >> 1;
+        src[8][n] = (sum222 + src12) >> 1;
+        src[9][n] = (sum222 - src12) >> 1;
+        src[10][n] = (dif222 + src13) >> 1;
+        src[11][n] = (dif222 - src13) >> 1;
+        src[12][n] = (sum333 + src14) >> 1;
+        src[13][n] = (sum333 - src14) >> 1;
+        src[14][n] = (dif333 + src15) >> 1;
+        src[15][n] = (dif333 - src15) >> 1;
+    }
+};
+
+}  // namespace xphoto
+}  // namespace cv
+
+#endif
\ No newline at end of file

modules/xphoto/src/bm3d_denoising_transforms_2D.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective icvers.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_BM3D_DENOISING_TRANSFORMS_2D_HPP__
+#define __OPENCV_BM3D_DENOISING_TRANSFORMS_2D_HPP__
+
+namespace cv
+{
+namespace xphoto
+{
+
+class HaarTransform2D
+{
+    template <int X>
+    static void CalculateIndices(unsigned *diffIndices, const unsigned &size)
+    {
+        unsigned diffIdx = 1;
+        unsigned diffAllIdx = 0;
+        for (unsigned i = 1; i <= X; i <<= 1)
+        {
+            diffAllIdx += (i >> 1);
+            for (unsigned j = 0; j < (i >> 1); ++j)
+                diffIndices[diffIdx++] = size - (--diffAllIdx);
+            diffAllIdx += i;
+        }
+    }
+
+public:
+    /// Transforms for 2D block of arbitrary size
+    template <typename T, typename TT, int X, int N>
+    inline static void ForwardTransformX(const T *src, TT *dst, const int &step)
+    {
+        const unsigned size = X + (X << 1) - 2;
+        TT dstX[size];
+
+        // Fill dstX with source values
+        for (unsigned i = 0; i < X; ++i)
+            dstX[i] = *(src + i * step);
+
+        unsigned idx = 0, dstIdx = X;
+        for (unsigned i = X; i > 1; i >>= 1)
+        {
+            // Get sums
+            for (unsigned j = 0; j < (i >> 1); ++j)
+                dstX[dstIdx++] = (dstX[idx + 2 * j] + dstX[idx + j * 2 + 1] + 1) >> 1;
+
+            // Get diffs
+            for (unsigned j = 0; j < (i >> 1); ++j)
+                dstX[dstIdx++] = dstX[idx + 2 * j] - dstX[idx + j * 2 + 1];
+
+            idx = dstIdx - i;
+        }
+
+        // Calculate indices in the destination matrix.
+        unsigned diffIndices[X];
+        CalculateIndices<X>(diffIndices, size);
+
+        // Fill in destination matrix
+        dst[0] = dstX[size - 2];
+        for (int i = 1; i < X; ++i)
+            dst[i * N] = dstX[diffIndices[i]];
+    }
+
+    template <typename T, typename TT, int X>
+    inline static void ForwardTransformXxX(const T *ptr, TT *dst, const int &step, const int /*blockSize*/)
+    {
+        TT temp[X * X];
+
+        // Transform columns first
+        for (unsigned i = 0; i < X; ++i)
+            ForwardTransformX<T, TT, X, X>(ptr + i, temp + i, step);
+
+        // Then transform rows
+        for (unsigned i = 0; i < X; ++i)
+            ForwardTransformX<TT, TT, X, 1>(temp + i * X, dst + i * X, 1);
+    }
+
+    template <typename T, int X, int N>
+    inline static void InverseTransformX(T *src, T *dst)
+    {
+        const unsigned dstSize = (X << 1) - 2;
+        T dstX[dstSize];
+        T srcX[X];
+
+        // Fill srcX with source values
+        srcX[0] = src[0] * 2;
+        for (int i = 1; i < X; ++i)
+            srcX[i] = src[i * N];
+
+        // Take care of first two elements
+        dstX[0] = srcX[0] + srcX[1];
+        dstX[1] = srcX[0] - srcX[1];
+
+        unsigned idx = 0, dstIdx = 2;
+        for (int i = 4; i < X; i <<= 1)
+        {
+            for (int j = 0; j < (i >> 1); ++j)
+            {
+                dstX[dstIdx++] = dstX[idx + j] + srcX[idx + 2 + j];
+                dstX[dstIdx++] = dstX[idx + j] - srcX[idx + 2 + j];
+            }
+            idx += (i >> 1);
+        }
+
+        // Handle the last X elements
+        dstIdx = 0;
+        for (int j = 0; j < (X >> 1); ++j)
+        {
+            dst[dstIdx++ * N] = (dstX[idx + j] + srcX[idx + 2 + j]) >> 1;
+            dst[dstIdx++ * N] = (dstX[idx + j] - srcX[idx + 2 + j]) >> 1;
+        }
+    }
+
+    template <typename T, int X>
+    inline static void InverseTransformXxX(T *src, const int /*blockSize*/)
+    {
+        T temp[X * X];
+
+        // Invert columns first
+        for (int i = 0; i < X; ++i)
+            InverseTransformX<T, X, X>(src + i, temp + i);
+
+        // Then invert rows
+        for (int i = 0; i < X; ++i)
+            InverseTransformX<T, X, 1>(temp + i * X, src + i * X);
+    }
+
+    // Same as above but X and N are arguments, not template parameters.
+
+    static void CalculateIndices(int *diffIndices, const int &size, const int &X)
+    {
+        int diffIdx = 1;
+        int diffAllIdx = 0;
+        for (int i = 1; i <= X; i <<= 1)
+        {
+            diffAllIdx += (i >> 1);
+            for (int j = 0; j < (i >> 1); ++j)
+                diffIndices[diffIdx++] = size - (--diffAllIdx);
+            diffAllIdx += i;
+        }
+    }
+
+    template <typename T, typename TT>
+    inline static void ForwardTransformX(const T *src, TT *dst, const int &step, const int &X, const int &N)
+    {
+        const int size = X + (X << 1) - 2;
+        TT *dstX = new TT[size];
+
+        // Fill dstX with source values
+        for (int i = 0; i < X; ++i)
+            dstX[i] = *(src + i * step);
+
+        int idx = 0, dstIdx = X;
+        for (int i = X; i > 1; i >>= 1)
+        {
+            // Get sums
+            for (int j = 0; j < (i >> 1); ++j)
+                dstX[dstIdx++] = (dstX[idx + 2 * j] + dstX[idx + j * 2 + 1] + 1) >> 1;
+
+            // Get diffs
+            for (int j = 0; j < (i >> 1); ++j)
+                dstX[dstIdx++] = dstX[idx + 2 * j] - dstX[idx + j * 2 + 1];
+
+            idx = dstIdx - i;
+        }
+
+        // Calculate indices in the destination matrix.
+        int *diffIndices = new int[X];
+        CalculateIndices(diffIndices, size, X);
+
+        // Fill in destination matrix
+        dst[0] = dstX[size - 2];
+        for (int i = 1; i < X; ++i)
+            dst[i * N] = dstX[diffIndices[i]];
+
+        delete[] diffIndices;
+        delete[] dstX;
+    }
+
+    template <typename T, typename TT>
+    inline static void ForwardTransformXxX(const T *ptr, TT *dst, const int &step, const int X)
+    {
+        TT *temp = new TT[X * X];
+
+        // Transform columns first
+        for (int i = 0; i < X; ++i)
+            ForwardTransformX<T, TT>(ptr + i, temp + i, step, X, X);
+
+        // Then transform rows
+        for (int i = 0; i < X; ++i)
+            ForwardTransformX<TT, TT>(temp + i * X, dst + i * X, 1, X, 1);
+
+        delete[] temp;
+    }
+
+    template <typename T>
+    inline static void InverseTransformX(T *src, T *dst, const int &X, const int &N)
+    {
+        const unsigned dstSize = (X << 1) - 2;
+        T *dstX = new T[dstSize];
+        T *srcX = new T[X];
+
+        // Fill srcX with source values
+        srcX[0] = src[0] * 2;
+        for (int i = 1; i < X; ++i)
+            srcX[i] = src[i * N];
+
+        // Take care of first two elements
+        dstX[0] = srcX[0] + srcX[1];
+        dstX[1] = srcX[0] - srcX[1];
+
+        unsigned idx = 0, dstIdx = 2;
+        for (int i = 4; i < X; i <<= 1)
+        {
+            for (int j = 0; j < (i >> 1); ++j)
+            {
+                dstX[dstIdx++] = dstX[idx + j] + srcX[idx + 2 + j];
+                dstX[dstIdx++] = dstX[idx + j] - srcX[idx + 2 + j];
+            }
+            idx += (i >> 1);
+        }
+
+        // Handle the last X elements
+        dstIdx = 0;
+        for (int j = 0; j < (X >> 1); ++j)
+        {
+            dst[dstIdx++ * N] = (dstX[idx + j] + srcX[idx + 2 + j]) >> 1;
+            dst[dstIdx++ * N] = (dstX[idx + j] - srcX[idx + 2 + j]) >> 1;
+        }
+
+        delete[] dstX;
+        delete[] srcX;
+    }
+
+    template <typename T>
+    inline static void InverseTransformXxX(T *src, const int X)
+    {
+        T *temp = new T[X * X];
+
+        // Invert columns first
+        for (int i = 0; i < X; ++i)
+            InverseTransformX<T>(src + i, temp + i, X, X);
+
+        // Then invert rows
+        for (int i = 0; i < X; ++i)
+            InverseTransformX<T>(temp + i * X, src + i * X, X, 1);
+
+        delete[] temp;
+    }
+
+    /// Transforms for 2x2 2D block
+
+    template <typename T, typename TT, int N>
+    inline static void ForwardTransform2(const T *src, TT *dst, const int &step)
+    {
+        const T *src0 = src;
+        const T *src1 = src + 1 * step;
+
+        dst[0 * N] = (*src0 + *src1 + 1) >> 1;
+        dst[1 * N] = *src0 - *src1;
+    }
+
+    template <typename T, typename TT>
+    inline static void ForwardTransform2x2(const T *ptr, TT *dst, const int &step, const int /*blockSize*/)
+    {
+        TT temp[4];
+
+        // Transform columns first
+        for (int i = 0; i < 2; ++i)
+            ForwardTransform2<T, TT, 2>(ptr + i, temp + i, step);
+
+        // Then transform rows
+        for (int i = 0; i < 2; ++i)
+            ForwardTransform2<TT, TT, 1>(temp + i * 2, dst + i * 2, 1);
+    }
+
+    template <typename TT, int N>
+    inline static void InverseTransform2(TT *src, TT *dst)
+    {
+        TT src0 = src[0 * N] * 2;
+        TT src1 = src[1 * N];
+
+        dst[0 * N] = (src0 + src1) >> 1;
+        dst[1 * N] = (src0 - src1) >> 1;
+    }
+
+    template <typename T>
+    inline static void InverseTransform2x2(T *src, const int /*blockSize*/)
+    {
+        T temp[4];
+
+        // Invert columns first
+        for (int i = 0; i < 2; ++i)
+            InverseTransform2<T, 2>(src + i, temp + i);
+
+        // Then invert rows
+        for (int i = 0; i < 2; ++i)
+            InverseTransform2<T, 1>(temp + i * 2, src + i * 2);
+    }
+
+    /// Transforms for 4x4 2D block
+
+    template <typename T, typename TT, int N>
+    inline static void ForwardTransform4(const T *src, TT *dst, const int &step)
+    {
+        const T *src0 = src;
+        const T *src1 = src + 1 * step;
+        const T *src2 = src + 2 * step;
+        const T *src3 = src + 3 * step;
+
+        TT sum0 = (*src0 + *src1 + 1) >> 1;
+        TT sum1 = (*src2 + *src3 + 1) >> 1;
+        TT dif0 = *src0 - *src1;
+        TT dif1 = *src2 - *src3;
+
+        TT sum00 = (sum0 + sum1 + 1) >> 1;
+        TT dif00 = sum0 - sum1;
+
+        dst[0 * N] = sum00;
+        dst[1 * N] = dif00;
+        dst[2 * N] = dif0;
+        dst[3 * N] = dif1;
+    }
+
+    template <typename T, typename TT>
+    inline static void ForwardTransform4x4(const T *ptr, TT *dst, const int &step, const int /*blockSize*/)
+    {
+        TT temp[16];
+
+        // Transform columns first
+        for (int i = 0; i < 4; ++i)
+            ForwardTransform4<T, TT, 4>(ptr + i, temp + i, step);
+
+        // Then transform rows
+        for (int i = 0; i < 4; ++i)
+            ForwardTransform4<TT, TT, 1>(temp + i * 4, dst + i * 4, 1);
+    }
+
+    template <typename TT, int N>
+    inline static void InverseTransform4(TT *src, TT *dst)
+    {
+        TT src0 = src[0 * N] * 2;
+        TT src1 = src[1 * N];
+        TT src2 = src[2 * N];
+        TT src3 = src[3 * N];
+
+        TT sum0 = src0 + src1;
+        TT dif0 = src0 - src1;
+
+        dst[0 * N] = (sum0 + src2) >> 1;
+        dst[1 * N] = (sum0 - src2) >> 1;
+        dst[2 * N] = (dif0 + src3) >> 1;
+        dst[3 * N] = (dif0 - src3) >> 1;
+    }
+
+    template <typename T>
+    inline static void InverseTransform4x4(T *src, const int /*blockSize*/)
+    {
+        T temp[16];
+
+        // Invert columns first
+        for (int i = 0; i < 4; ++i)
+            InverseTransform4<T, 4>(src + i, temp + i);
+
+        // Then invert rows
+        for (int i = 0; i < 4; ++i)
+            InverseTransform4<T, 1>(temp + i * 4, src + i * 4);
+    }
+
+    /// Transforms for 8x8 2D block
+
+    template <typename T, typename TT, int N>
+    inline static void ForwardTransform8(const T *src, TT *dst, const int &step)
+    {
+        const T *src0 = src;
+        const T *src1 = src + 1 * step;
+        const T *src2 = src + 2 * step;
+        const T *src3 = src + 3 * step;
+        const T *src4 = src + 4 * step;
+        const T *src5 = src + 5 * step;
+        const T *src6 = src + 6 * step;
+        const T *src7 = src + 7 * step;
+
+        TT sum0 = (*src0 + *src1 + 1) >> 1;
+        TT sum1 = (*src2 + *src3 + 1) >> 1;
+        TT sum2 = (*src4 + *src5 + 1) >> 1;
+        TT sum3 = (*src6 + *src7 + 1) >> 1;
+        TT dif0 = *src0 - *src1;
+        TT dif1 = *src2 - *src3;
+        TT dif2 = *src4 - *src5;
+        TT dif3 = *src6 - *src7;
+
+        TT sum00 = (sum0 + sum1 + 1) >> 1;
+        TT sum11 = (sum2 + sum3 + 1) >> 1;
+        TT dif00 = sum0 - sum1;
+        TT dif11 = sum2 - sum3;
+
+        TT sum000 = (sum00 + sum11 + 1) >> 1;
+        TT dif000 = sum00 - sum11;
+
+        dst[0 * N] = sum000;
+        dst[1 * N] = dif000;
+        dst[2 * N] = dif00;
+        dst[3 * N] = dif11;
+        dst[4 * N] = dif0;
+        dst[5 * N] = dif1;
+        dst[6 * N] = dif2;
+        dst[7 * N] = dif3;
+    }
+
+    template <typename T, typename TT>
+    inline static void ForwardTransform8x8(const T *ptr, TT *dst, const int &step, const int /*blockSize*/)
+    {
+        TT temp[64];
+
+        // Transform columns first
+        for (int i = 0; i < 8; ++i)
+            ForwardTransform8<T, TT, 8>(ptr + i, temp + i, step);
+
+        // Then transform rows
+        for (int i = 0; i < 8; ++i)
+            ForwardTransform8<TT, TT, 1>(temp + i * 8, dst + i * 8, 1);
+    }
+
+    template <typename T, int N>
+    inline static void InverseTransform8(T *src, T *dst)
+    {
+        T src0 = src[0] * 2;
+        T src1 = src[1 * N];
+        T src2 = src[2 * N];
+        T src3 = src[3 * N];
+        T src4 = src[4 * N];
+        T src5 = src[5 * N];
+        T src6 = src[6 * N];
+        T src7 = src[7 * N];
+
+        T sum0 = src0 + src1;
+        T dif0 = src0 - src1;
+
+        T sum00 = sum0 + src2;
+        T dif00 = sum0 - src2;
+        T sum11 = dif0 + src3;
+        T dif11 = dif0 - src3;
+
+        dst[0 * N] = (sum00 + src4) >> 1;
+        dst[1 * N] = (sum00 - src4) >> 1;
+        dst[2 * N] = (dif00 + src5) >> 1;
+        dst[3 * N] = (dif00 - src5) >> 1;
+        dst[4 * N] = (sum11 + src6) >> 1;
+        dst[5 * N] = (sum11 - src6) >> 1;
+        dst[6 * N] = (dif11 + src7) >> 1;
+        dst[7 * N] = (dif11 - src7) >> 1;
+    }
+
+    template <typename T>
+    inline static void InverseTransform8x8(T *src, const int /*blockSize*/)
+    {
+        T temp[64];
+
+        // Invert columns first
+        for (int i = 0; i < 8; ++i)
+            InverseTransform8<T, 8>(src + i, temp + i);
+
+        // Then invert rows
+        for (int i = 0; i < 8; ++i)
+            InverseTransform8<T, 1>(temp + i * 8, src + i * 8);
+    }
+};
+
+}  // namespace xphoto
+}  // namespace cv
+
+#endif
\ No newline at end of file

modules/xphoto/src/bm3d_denoising_transforms_haar.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective icvers.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_BM3D_DENOISING_TRANSFORMS_HAAR_HPP__
+#define __OPENCV_BM3D_DENOISING_TRANSFORMS_HAAR_HPP__
+
+#include "bm3d_denoising_transforms_1D.hpp"
+#include "bm3d_denoising_transforms_2D.hpp"
+
+namespace cv
+{
+namespace xphoto
+{
+
+// Forward declaration
+template <typename T, typename TT>
+class Transform;
+
+template <typename T, typename TT>
+class HaarTransform
+{
+    static void calcCoefficients1D(cv::Mat &coeff1D, const int &numberOfElements)
+    {
+        // Generate base array and initialize with zeros
+        cv::Mat baseArr = cv::Mat::zeros(numberOfElements, numberOfElements, CV_32FC1);
+
+        // Calculate base array coefficients.
+        int currentRow = 0;
+        for (int i = numberOfElements; i > 0; i /= 2)
+        {
+            for (int k = 0, sign = -1; k < numberOfElements; ++k)
+            {
+                // Alternate sign every i-th element
+                if (k % i == 0)
+                    sign *= -1;
+
+                // Move to the next row every 2*i-th element
+                if (k != 0 && (k % (2 * i) == 0))
+                    ++currentRow;
+
+                baseArr.at<float>(currentRow, k) = sign * 1.0f / i;
+            }
+            ++currentRow;
+        }
+
+        // Square each elements of the base array
+        float *ptr = baseArr.ptr<float>(0);
+        for (unsigned i = 0; i < baseArr.total(); ++i)
+            ptr[i] = ptr[i] * ptr[i];
+
+        // Multiply baseArray with 1D vector of ones
+        cv::Mat unitaryArr = cv::Mat::ones(numberOfElements, 1, CV_32FC1);
+        coeff1D = baseArr * unitaryArr;
+    }
+
+    // Method to generate threshold coefficients for 1D transform depending on the number of elements.
+    static void fillHaarCoefficients1D(float *thrCoeff1D, int &idx, const int &numberOfElements)
+    {
+        cv::Mat coeff1D;
+        calcCoefficients1D(coeff1D, numberOfElements);
+
+        // Square root the array to get standard deviation
+        float *ptr = coeff1D.ptr<float>(0);
+        for (unsigned i = 0; i < coeff1D.total(); ++i)
+        {
+            ptr[i] = std::sqrt(ptr[i]);
+            thrCoeff1D[idx++] = ptr[i];
+        }
+    }
+
+    // Method to generate threshold coefficients for 2D transform depending on the number of elements.
+    static void fillHaarCoefficients2D(float *thrCoeff2D, const int &templateWindowSize)
+    {
+        cv::Mat coeff1D;
+        calcCoefficients1D(coeff1D, templateWindowSize);
+
+        // Calculate 2D array
+        cv::Mat coeff1Dt;
+        cv::transpose(coeff1D, coeff1Dt);
+        cv::Mat coeff2D = coeff1D * coeff1Dt;
+
+        // Square root the array to get standard deviation
+        float *ptr = coeff2D.ptr<float>(0);
+        for (unsigned i = 0; i < coeff2D.total(); ++i)
+            thrCoeff2D[i] = std::sqrt(ptr[i]);
+    }
+
+public:
+    // Method to calculate 1D threshold map based on the maximum number of elements
+    // Allocates memory for the output array.
+    static void calcThresholdMap1D(float *&thrMap1D, const int &numberOfElements)
+    {
+        CV_Assert(numberOfElements > 0);
+
+        // Allocate memory for the array
+        const int arrSize = (numberOfElements << 1) - 1;
+        if (thrMap1D == NULL)
+            thrMap1D = new float[arrSize];
+
+        for (int i = 1, idx = 0; i <= numberOfElements; i *= 2)
+            fillHaarCoefficients1D(thrMap1D, idx, i);
+    }
+
+    // Method to calculate 2D threshold map based on the maximum number of elements
+    // Allocates memory for the output array.
+    static void calcThresholdMap2D(float *&thrMap2D, const int &templateWindowSize)
+    {
+        // Allocate memory for the array
+        if (thrMap2D == NULL)
+            thrMap2D = new float[templateWindowSize * templateWindowSize];
+
+        fillHaarCoefficients2D(thrMap2D, templateWindowSize);
+    }
+
+    // Method to calculate 3D threshold map based on the maximum number of elements.
+    // Allocates memory for the output array.
+    static void calcThresholdMap3D(
+        TT *&outThrMap1D,
+        const float &hardThr1D,
+        const int &templateWindowSize,
+        const int &groupSize)
+    {
+        const int templateWindowSizeSq = templateWindowSize * templateWindowSize;
+
+        // Allocate memory for the output array
+        if (outThrMap1D == NULL)
+            outThrMap1D = new TT[templateWindowSizeSq * ((groupSize << 1) - 1)];
+
+        // Generate 1D coefficients map
+        float *thrMap1D = NULL;
+        calcThresholdMap1D(thrMap1D, groupSize);
+
+        // Generate 2D coefficients map
+        float *thrMap2D = NULL;
+        calcThresholdMap2D(thrMap2D, templateWindowSize);
+
+        // Generate 3D threshold map
+        TT *thrMapPtr1D = outThrMap1D;
+        for (int i = 1, ii = 0; i <= groupSize; ++ii, i *= 2)
+        {
+            float coeff = (i == 1) ? 1.0f : std::sqrt(2.0f * std::log((float)i));
+            for (int jj = 0; jj < templateWindowSizeSq; ++jj)
+            {
+                for (int ii1 = 0; ii1 < (1 << ii); ++ii1)
+                {
+                    int indexIn1D = (1 << ii) - 1 + ii1;
+                    int indexIn2D = jj;
+                    int thr = static_cast<int>(thrMap1D[indexIn1D] * thrMap2D[indexIn2D] * hardThr1D * coeff);
+
+                    // Set DC component to zero
+                    if (jj == 0 && ii1 == 0)
+                        thr = 0;
+
+                    *thrMapPtr1D++ = cv::saturate_cast<TT>(thr);
+                }
+            }
+        }
+
+        delete[] thrMap1D;
+        delete[] thrMap2D;
+    }
+
+    // Method that registers 2D transform calls
+    static void RegisterTransforms2D(const int &templateWindowSize)
+    {
+        // Check if template window size is a power of two
+        if (!isPowerOf2(templateWindowSize))
+            CV_Error(Error::StsBadArg, "Unsupported template size! Template size must be power of two!");
+
+        switch (templateWindowSize)
+        {
+        case 2:
+            forwardTransform2D = HaarTransform2D::ForwardTransform2x2<T, TT>;
+            inverseTransform2D = HaarTransform2D::InverseTransform2x2<TT>;
+            break;
+        case 4:
+            forwardTransform2D = HaarTransform2D::ForwardTransform4x4<T, TT>;
+            inverseTransform2D = HaarTransform2D::InverseTransform4x4<TT>;
+            break;
+        case 8:
+            forwardTransform2D = HaarTransform2D::ForwardTransform8x8<T, TT>;
+            inverseTransform2D = HaarTransform2D::InverseTransform8x8<TT>;
+            break;
+        case 16:
+            forwardTransform2D = HaarTransform2D::ForwardTransformXxX<T, TT, 16>;
+            inverseTransform2D = HaarTransform2D::InverseTransformXxX<TT, 16>;
+            break;
+        case 32:
+            forwardTransform2D = HaarTransform2D::ForwardTransformXxX<T, TT, 32>;
+            inverseTransform2D = HaarTransform2D::InverseTransformXxX<TT, 32>;
+            break;
+        case 64:
+            forwardTransform2D = HaarTransform2D::ForwardTransformXxX<T, TT, 64>;
+            inverseTransform2D = HaarTransform2D::InverseTransformXxX<TT, 64>;
+            break;
+        default:
+            forwardTransform2D = HaarTransform2D::ForwardTransformXxX<T, TT>;
+            inverseTransform2D = HaarTransform2D::InverseTransformXxX<TT>;
+        }
+    }
+
+    // 2D transform pointers
+    static typename Transform<T, TT>::Forward2D forwardTransform2D;
+    static typename Transform<T, TT>::Inverse2D inverseTransform2D;
+
+    // 1D transform pointers
+    static typename Transform<T, TT>::Forward1D forwardTransformN;
+    static typename Transform<T, TT>::Inverse1D inverseTransformN;
+
+    // Specialized 1D forward transform pointers
+    static typename Transform<T, TT>::Forward1Ds forwardTransform2;
+    static typename Transform<T, TT>::Forward1Ds forwardTransform4;
+    static typename Transform<T, TT>::Forward1Ds forwardTransform8;
+    static typename Transform<T, TT>::Forward1Ds forwardTransform16;
+
+    // Specialized 1D inverse transform pointers
+    static typename Transform<T, TT>::Inverse1Ds inverseTransform2;
+    static typename Transform<T, TT>::Inverse1Ds inverseTransform4;
+    static typename Transform<T, TT>::Inverse1Ds inverseTransform8;
+    static typename Transform<T, TT>::Inverse1Ds inverseTransform16;
+};
+
+/// Explicit static members initialization
+
+#define INITIALIZE_HAAR_TRANSFORM(type, member, value)              \
+template <typename T, typename TT>                                  \
+typename Transform<T, TT>::type HaarTransform<T, TT>::member = value;
+
+// 2D transforms
+INITIALIZE_HAAR_TRANSFORM(Forward2D, forwardTransform2D, NULL)
+INITIALIZE_HAAR_TRANSFORM(Inverse2D, inverseTransform2D, NULL)
+
+// 1D transforms
+INITIALIZE_HAAR_TRANSFORM(Forward1D, forwardTransformN, HaarTransform1D::ForwardTransformN)
+INITIALIZE_HAAR_TRANSFORM(Inverse1D, inverseTransformN, HaarTransform1D::InverseTransformN)
+
+// Specialized 1D forward transforms
+INITIALIZE_HAAR_TRANSFORM(Forward1Ds, forwardTransform2, HaarTransform1D::ForwardTransform2)
+INITIALIZE_HAAR_TRANSFORM(Forward1Ds, forwardTransform4, HaarTransform1D::ForwardTransform4)
+INITIALIZE_HAAR_TRANSFORM(Forward1Ds, forwardTransform8, HaarTransform1D::ForwardTransform8)
+INITIALIZE_HAAR_TRANSFORM(Forward1Ds, forwardTransform16, HaarTransform1D::ForwardTransform16)
+
+// Specialized 1D inverse transforms
+INITIALIZE_HAAR_TRANSFORM(Inverse1Ds, inverseTransform2, HaarTransform1D::InverseTransform2)
+INITIALIZE_HAAR_TRANSFORM(Inverse1Ds, inverseTransform4, HaarTransform1D::InverseTransform4)
+INITIALIZE_HAAR_TRANSFORM(Inverse1Ds, inverseTransform8, HaarTransform1D::InverseTransform8)
+INITIALIZE_HAAR_TRANSFORM(Inverse1Ds, inverseTransform16, HaarTransform1D::InverseTransform16)
+
+}  // namespace xphoto
+}  // namespace cv
+
+#endif
\ No newline at end of file

modules/xphoto/src/bm3d_image_denoising.cpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective icvers.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "opencv2/xphoto.hpp"
+#include "opencv2/core.hpp"
+
+#ifdef OPENCV_ENABLE_NONFREE
+
+#include "bm3d_denoising_invoker_step1.hpp"
+#include "bm3d_denoising_invoker_step2.hpp"
+#include "bm3d_denoising_transforms.hpp"
+
+#endif
+
+namespace cv
+{
+namespace xphoto
+{
+
+#ifdef OPENCV_ENABLE_NONFREE
+
+template<typename ST, typename D, typename TT>
+static void bm3dDenoising_(
+    const Mat& src,
+    Mat& basic,
+    Mat& dst,
+    const float& h,
+    const int &templateWindowSize,
+    const int &searchWindowSize,
+    const int &hBMStep1,
+    const int &hBMStep2,
+    const int &groupSize,
+    const int &slidingStep,
+    const float &beta,
+    const int &step)
+{
+    double granularity = (double)std::max(1., (double)src.total() / (1 << 16));
+
+    switch (CV_MAT_CN(src.type())) {
+    case 1:
+        if (step == BM3D_STEP1 || step == BM3D_STEPALL)
+        {
+            parallel_for_(cv::Range(0, src.rows),
+                Bm3dDenoisingInvokerStep1<ST, D, float, TT, HaarTransform<ST, TT> >(
+                    src,
+                    basic,
+                    templateWindowSize,
+                    searchWindowSize,
+                    h,
+                    hBMStep1,
+                    groupSize,
+                    slidingStep,
+                    beta),
+                granularity);
+        }
+        if (step == BM3D_STEP2 || step == BM3D_STEPALL)
+        {
+            parallel_for_(cv::Range(0, src.rows),
+                Bm3dDenoisingInvokerStep2<ST, D, float, TT, HaarTransform<ST, TT> >(
+                    src,
+                    basic,
+                    dst,
+                    templateWindowSize,
+                    searchWindowSize,
+                    h,
+                    hBMStep2,
+                    groupSize,
+                    slidingStep,
+                    beta),
+                granularity);
+        }
+        break;
+    default:
+        CV_Error(Error::StsBadArg,
+            "Unsupported number of channels! Only 1 channel is supported at the moment.");
+    }
+}
+
+void bm3dDenoising(
+    InputArray _src,
+    InputOutputArray _basic,
+    OutputArray _dst,
+    float h,
+    int templateWindowSize,
+    int searchWindowSize,
+    int blockMatchingStep1,
+    int blockMatchingStep2,
+    int groupSize,
+    int slidingStep,
+    float beta,
+    int normType,
+    int step,
+    int transformType)
+{
+    int type = _src.type(), depth = CV_MAT_DEPTH(type), cn = CV_MAT_CN(type);
+    CV_Assert(1 == cn);
+    CV_Assert(HAAR == transformType);
+    CV_Assert(searchWindowSize > templateWindowSize);
+    CV_Assert(slidingStep > 0 && slidingStep < templateWindowSize);
+
+    Size srcSize = _src.size();
+
+    switch (step)
+    {
+    case BM3D_STEP1:
+        _basic.create(srcSize, type);
+        break;
+    case BM3D_STEP2:
+        CV_Assert(type == _basic.type());
+        _dst.create(srcSize, type);
+        break;
+    case BM3D_STEPALL:
+        _dst.create(srcSize, type);
+        break;
+    default:
+        CV_Error(Error::StsBadArg, "Unsupported BM3D step!");
+    }
+
+    Mat src = _src.getMat();
+    Mat basic = _basic.getMat().empty() ? Mat(srcSize, type) : _basic.getMat();
+    Mat dst = _dst.getMat();
+
+    switch (normType) {
+    case cv::NORM_L2:
+        switch (depth) {
+        case CV_8U:
+            bm3dDenoising_<uchar, DistSquared, short>(
+                src,
+                basic,
+                dst,
+                h,
+                templateWindowSize,
+                searchWindowSize,
+                blockMatchingStep1,
+                blockMatchingStep2,
+                groupSize,
+                slidingStep,
+                beta,
+                step);
+            break;
+        default:
+            CV_Error(Error::StsBadArg,
+                "Unsupported depth! Only CV_8U is supported for NORM_L2");
+        }
+        break;
+    case cv::NORM_L1:
+        switch (depth) {
+        case CV_8U:
+            bm3dDenoising_<uchar, DistAbs, short>(
+                src,
+                basic,
+                dst,
+                h,
+                templateWindowSize,
+                searchWindowSize,
+                blockMatchingStep1,
+                blockMatchingStep2,
+                groupSize,
+                slidingStep,
+                beta,
+                step);
+            break;
+        case CV_16U:
+            bm3dDenoising_<ushort, DistAbs, int>(
+                src,
+                basic,
+                dst,
+                h,
+                templateWindowSize,
+                searchWindowSize,
+                blockMatchingStep1,
+                blockMatchingStep2,
+                groupSize,
+                slidingStep,
+                beta,
+                step);
+            break;
+        default:
+            CV_Error(Error::StsBadArg,
+                "Unsupported depth! Only CV_8U and CV_16U are supported for NORM_L1");
+        }
+        break;
+    default:
+        CV_Error(Error::StsBadArg,
+            "Unsupported norm type! Only NORM_L2 and NORM_L1 are supported");
+    }
+}
+
+void bm3dDenoising(
+    InputArray _src,
+    OutputArray _dst,
+    float h,
+    int templateWindowSize,
+    int searchWindowSize,
+    int blockMatchingStep1,
+    int blockMatchingStep2,
+    int groupSize,
+    int slidingStep,
+    float beta,
+    int normType,
+    int step,
+    int transformType)
+{
+    if (step == BM3D_STEP2)
+        CV_Error(Error::StsBadArg,
+            "Unsupported step type! To use BM3D_STEP2 one need to provide basic image.");
+
+    Mat basic;
+
+    bm3dDenoising(
+        _src,
+        basic,
+        _dst,
+        h,
+        templateWindowSize,
+        searchWindowSize,
+        blockMatchingStep1,
+        blockMatchingStep2,
+        groupSize,
+        slidingStep,
+        beta,
+        normType,
+        step,
+        transformType);
+
+    if (step == BM3D_STEP1)
+        _dst.assign(basic);
+}
+
+#else
+
+void bm3dDenoising(
+    InputArray _src,
+    InputOutputArray _basic,
+    OutputArray _dst,
+    float h,
+    int templateWindowSize,
+    int searchWindowSize,
+    int blockMatchingStep1,
+    int blockMatchingStep2,
+    int groupSize,
+    int slidingStep,
+    float beta,
+    int normType,
+    int step,
+    int transformType)
+{
+    // Empty implementation
+
+    CV_UNUSED(_src);
+    CV_UNUSED(_basic);
+    CV_UNUSED(_dst);
+    CV_UNUSED(h);
+    CV_UNUSED(templateWindowSize);
+    CV_UNUSED(searchWindowSize);
+    CV_UNUSED(blockMatchingStep1);
+    CV_UNUSED(blockMatchingStep2);
+    CV_UNUSED(groupSize);
+    CV_UNUSED(slidingStep);
+    CV_UNUSED(beta);
+    CV_UNUSED(normType);
+    CV_UNUSED(step);
+    CV_UNUSED(transformType);
+
+    CV_Error(Error::StsNotImplemented,
+        "This algorithm is patented and is excluded in this configuration;"
+        "Set OPENCV_ENABLE_NONFREE CMake option and rebuild the library");
+}
+
+void bm3dDenoising(
+    InputArray _src,
+    OutputArray _dst,
+    float h,
+    int templateWindowSize,
+    int searchWindowSize,
+    int blockMatchingStep1,
+    int blockMatchingStep2,
+    int groupSize,
+    int slidingStep,
+    float beta,
+    int normType,
+    int step,
+    int transformType)
+{
+    // Empty implementation
+
+    CV_UNUSED(_src);
+    CV_UNUSED(_dst);
+    CV_UNUSED(h);
+    CV_UNUSED(templateWindowSize);
+    CV_UNUSED(searchWindowSize);
+    CV_UNUSED(blockMatchingStep1);
+    CV_UNUSED(blockMatchingStep2);
+    CV_UNUSED(groupSize);
+    CV_UNUSED(slidingStep);
+    CV_UNUSED(beta);
+    CV_UNUSED(normType);
+    CV_UNUSED(step);
+    CV_UNUSED(transformType);
+
+    CV_Error(Error::StsNotImplemented,
+        "This algorithm is patented and is excluded in this configuration;"
+        "Set OPENCV_ENABLE_NONFREE CMake option and rebuild the library");
+}
+
+#endif
+
+}  // namespace xphoto
+}  // namespace cv
\ No newline at end of file

modules/xphoto/src/kaiser_window.hpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                        Intel License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000, Intel Corporation, all rights reserved.
+// Third party copyrights are property of their respective icvers.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of Intel Corporation may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#ifndef __OPENCV_BM3D_DENOISING_KAISER_WINDOW_HPP__
+#define __OPENCV_BM3D_DENOISING_KAISER_WINDOW_HPP__
+
+#include "opencv2/core.hpp"
+#include <cmath>
+
+namespace cv
+{
+namespace xphoto
+{
+
+static int factorial(int n)
+{
+    if (n == 0)
+        return 1;
+
+    int val = 1;
+    for (int idx = 1; idx <= n; ++idx)
+        val *= idx;
+
+    return val;
+}
+
+template <int MAX_ITER>
+static float bessel0(const float &x)
+{
+    float sum = 0.0f;
+
+    for (int m = 0; m < MAX_ITER; ++m)
+    {
+        float factM = (float)factorial(m);
+        float inc = std::pow(1.0f / factM * std::pow(x * 0.5f, (float)m), 2.0f);
+        sum += inc;
+
+        if ((inc / sum) < 0.001F)
+            break;
+    }
+
+    return sum;
+}
+
+#define MAX_ITER_BESSEL 100
+
+static void calcKaiserWindow1D(cv::Mat &dst, const int N, const float beta)
+{
+    if (dst.empty())
+        dst.create(cv::Size(1, N), CV_32FC1);
+
+    CV_Assert(dst.total() == (size_t)N);
+    CV_Assert(dst.type() == CV_32FC1);
+    CV_Assert(N > 0);
+
+    float *p = dst.ptr<float>(0);
+    for (int i = 0; i < N; ++i)
+    {
+        float b = beta * std::sqrt(1.0f - std::pow(2.0f * i / (N - 1.0f) - 1.0f, 2.0f));
+        p[i] = bessel0<MAX_ITER_BESSEL>(b) / bessel0<MAX_ITER_BESSEL>(beta);
+    }
+}
+
+static void calcKaiserWindow2D(float *&kaiser, const int N, const float beta)
+{
+    if (kaiser == NULL)
+        kaiser = new float[N * N];
+
+    if (beta == 0.0f)
+    {
+        for (int i = 0; i < N * N; ++i)
+            kaiser[i] = 1.0f;
+        return;
+    }
+
+    cv::Mat kaiser1D;
+    calcKaiserWindow1D(kaiser1D, N, beta);
+
+    cv::Mat kaiser1Dt;
+    cv::transpose(kaiser1D, kaiser1Dt);
+
+    cv::Mat kaiser2D = kaiser1D * kaiser1Dt;
+    float *p = kaiser2D.ptr<float>(0);
+    for (unsigned i = 0; i < kaiser2D.total(); ++i)
+        kaiser[i] = p[i];
+}
+
+}  // namespace xphoto
+}  // namespace cv
+
+#endif
\ No newline at end of file

modules/xphoto/test/test_denoise_bm3d.cpp

+/*M///////////////////////////////////////////////////////////////////////////////////////
+//
+//  IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING.
+//
+//  By downloading, copying, installing or using the software you agree to this license.
+//  If you do not agree to this license, do not download, install,
+//  copy or use the software.
+//
+//
+//                           License Agreement
+//                For Open Source Computer Vision Library
+//
+// Copyright (C) 2000-2008, Intel Corporation, all rights reserved.
+// Copyright (C) 2009, Willow Garage Inc., all rights reserved.
+// Third party copyrights are property of their respective owners.
+//
+// Redistribution and use in source and binary forms, with or without modification,
+// are permitted provided that the following conditions are met:
+//
+//   * Redistribution's of source code must retain the above copyright notice,
+//     this list of conditions and the following disclaimer.
+//
+//   * Redistribution's in binary form must reproduce the above copyright notice,
+//     this list of conditions and the following disclaimer in the documentation
+//     and/or other materials provided with the distribution.
+//
+//   * The name of the copyright holders may not be used to endorse or promote products
+//     derived from this software without specific prior written permission.
+//
+// This software is provided by the copyright holders and contributors "as is" and
+// any express or implied warranties, including, but not limited to, the implied
+// warranties of merchantability and fitness for a particular purpose are disclaimed.
+// In no event shall the Intel Corporation or contributors be liable for any direct,
+// indirect, incidental, special, exemplary, or consequential damages
+// (including, but not limited to, procurement of substitute goods or services;
+// loss of use, data, or profits; or business interruption) however caused
+// and on any theory of liability, whether in contract, strict liability,
+// or tort (including negligence or otherwise) arising in any way out of
+// the use of this software, even if advised of the possibility of such damage.
+//
+//M*/
+
+#include "test_precomp.hpp"
+#include <string>
+
+//#define DUMP_RESULTS
+//#define TEST_TRANSFORMS
+
+#ifdef TEST_TRANSFORMS
+#include "..\..\xphoto\src\bm3d_denoising_invoker_commons.hpp"
+#include "..\..\xphoto\src\bm3d_denoising_transforms.hpp"
+#include "..\..\xphoto\src\kaiser_window.hpp"
+using namespace cv::xphoto;
+#endif
+
+#ifdef DUMP_RESULTS
+#  define DUMP(image, path) imwrite(path, image)
+#else
+#  define DUMP(image, path)
+#endif
+
+#ifdef OPENCV_ENABLE_NONFREE
+
+namespace cvtest
+{
+    TEST(xphoto_DenoisingBm3dGrayscale, regression_L2)
+    {
+        std::string folder = std::string(cvtest::TS::ptr()->get_data_path()) + "cv/xphoto/bm3d_image_denoising/";
+        std::string original_path = folder + "lena_noised_gaussian_sigma=10.png";
+        std::string expected_path = folder + "lena_noised_denoised_bm3d_wiener_grayscale_l2_tw=4_sw=16_h=10_bm=400.png";
+
+        cv::Mat original = cv::imread(original_path, cv::IMREAD_GRAYSCALE);
+        cv::Mat expected = cv::imread(expected_path, cv::IMREAD_GRAYSCALE);
+
+        ASSERT_FALSE(original.empty()) << "Could not load input image " << original_path;
+        ASSERT_FALSE(expected.empty()) << "Could not load reference image " << expected_path;
+
+        // BM3D: two different calls doing exactly the same thing
+        cv::Mat result, resultSec;
+        cv::xphoto::bm3dDenoising(original, cv::Mat(), resultSec, 10, 4, 16, 2500, 400, 8, 1, 0.0f, cv::NORM_L2, cv::xphoto::BM3D_STEPALL);
+        cv::xphoto::bm3dDenoising(original, result, 10, 4, 16, 2500, 400, 8, 1, 0.0f, cv::NORM_L2, cv::xphoto::BM3D_STEPALL);
+
+        DUMP(result, expected_path + ".res.png");
+
+        ASSERT_EQ(cvtest::norm(result, resultSec, cv::NORM_L2), 0);
+        ASSERT_LT(cvtest::norm(result, expected, cv::NORM_L2), 200);
+    }
+
+    TEST(xphoto_DenoisingBm3dGrayscale, regression_L2_separate)
+    {
+        std::string folder = std::string(cvtest::TS::ptr()->get_data_path()) + "cv/xphoto/bm3d_image_denoising/";
+        std::string original_path = folder + "lena_noised_gaussian_sigma=10.png";
+        std::string expected_basic_path = folder + "lena_noised_denoised_bm3d_grayscale_l2_tw=4_sw=16_h=10_bm=2500.png";
+        std::string expected_path = folder + "lena_noised_denoised_bm3d_wiener_grayscale_l2_tw=4_sw=16_h=10_bm=400.png";
+
+        cv::Mat original = cv::imread(original_path, cv::IMREAD_GRAYSCALE);
+        cv::Mat expected_basic = cv::imread(expected_basic_path, cv::IMREAD_GRAYSCALE);
+        cv::Mat expected = cv::imread(expected_path, cv::IMREAD_GRAYSCALE);
+
+        ASSERT_FALSE(original.empty()) << "Could not load input image " << original_path;
+        ASSERT_FALSE(expected_basic.empty()) << "Could not load reference image " << expected_basic_path;
+        ASSERT_FALSE(expected.empty()) << "Could not load input image " << expected_path;
+
+        cv::Mat basic, result;
+
+        // BM3D step 1
+        cv::xphoto::bm3dDenoising(original, basic, 10, 4, 16, 2500, -1, 8, 1, 0.0f, cv::NORM_L2, cv::xphoto::BM3D_STEP1);
+        ASSERT_LT(cvtest::norm(basic, expected_basic, cv::NORM_L2), 200);
+        DUMP(basic, expected_basic_path + ".res.basic.png");
+
+        // BM3D step 2
+        cv::xphoto::bm3dDenoising(original, basic, result, 10, 4, 16, 2500, 400, 8, 1, 0.0f, cv::NORM_L2, cv::xphoto::BM3D_STEP2);
+        ASSERT_LT(cvtest::norm(basic, expected_basic, cv::NORM_L2), 200);
+        DUMP(basic, expected_basic_path + ".res.basic2.png");
+
+        DUMP(result, expected_path + ".res.png");
+
+        ASSERT_LT(cvtest::norm(result, expected, cv::NORM_L2), 200);
+    }
+
+    TEST(xphoto_DenoisingBm3dGrayscale, regression_L1)
+    {
+        std::string folder = std::string(cvtest::TS::ptr()->get_data_path()) + "cv/xphoto/bm3d_image_denoising/";
+        std::string original_path = folder + "lena_noised_gaussian_sigma=10.png";
+        std::string expected_path = folder + "lena_noised_denoised_bm3d_grayscale_l1_tw=4_sw=16_h=10_bm=2500.png";
+
+        cv::Mat original = cv::imread(original_path, cv::IMREAD_GRAYSCALE);
+        cv::Mat expected = cv::imread(expected_path, cv::IMREAD_GRAYSCALE);
+
+        ASSERT_FALSE(original.empty()) << "Could not load input image " << original_path;
+        ASSERT_FALSE(expected.empty()) << "Could not load reference image " << expected_path;
+
+        cv::Mat result;
+        cv::xphoto::bm3dDenoising(original, result, 10, 4, 16, 2500, -1, 8, 1, 0.0f, cv::NORM_L1, cv::xphoto::BM3D_STEP1);
+
+        DUMP(result, expected_path + ".res.png");
+
+        ASSERT_LT(cvtest::norm(result, expected, cv::NORM_L2), 200);
+    }
+
+    TEST(xphoto_DenoisingBm3dGrayscale, regression_L2_8x8)
+    {
+        std::string folder = std::string(cvtest::TS::ptr()->get_data_path()) + "cv/xphoto/bm3d_image_denoising/";
+        std::string original_path = folder + "lena_noised_gaussian_sigma=10.png";
+        std::string expected_path = folder + "lena_noised_denoised_bm3d_grayscale_l2_tw=8_sw=16_h=10_bm=2500.png";
+
+        cv::Mat original = cv::imread(original_path, cv::IMREAD_GRAYSCALE);
+        cv::Mat expected = cv::imread(expected_path, cv::IMREAD_GRAYSCALE);
+
+        ASSERT_FALSE(original.empty()) << "Could not load input image " << original_path;
+        ASSERT_FALSE(expected.empty()) << "Could not load reference image " << expected_path;
+
+        cv::Mat result;
+        cv::xphoto::bm3dDenoising(original, result, 10, 8, 16, 2500, -1, 8, 1, 0.0f, cv::NORM_L2, cv::xphoto::BM3D_STEP1);
+
+        DUMP(result, expected_path + ".res.png");
+
+        ASSERT_LT(cvtest::norm(result, expected, cv::NORM_L2), 200);
+    }
+
+#ifdef TEST_TRANSFORMS
+
+    TEST(xphoto_DenoisingBm3dKaiserWindow, regression_4)
+    {
+        float beta = 2.0f;
+        int N = 4;
+
+        cv::Mat kaiserWindow;
+        calcKaiserWindow1D(kaiserWindow, N, beta);
+
+        float kaiser4[] = {
+            0.43869004f,
+            0.92432547f,
+            0.92432547f,
+            0.43869004f
+        };
+
+        for (int i = 0; i < N; ++i)
+            ASSERT_FLOAT_EQ(kaiser4[i], kaiserWindow.at<float>(i));
+    }
+
+    TEST(xphoto_DenoisingBm3dKaiserWindow, regression_8)
+    {
+        float beta = 2.0f;
+        int N = 8;
+
+        cv::Mat kaiserWindow;
+        calcKaiserWindow1D(kaiserWindow, N, beta);
+
+        float kaiser8[] = {
+            0.43869004f,
+            0.68134475f,
+            0.87685609f,
+            0.98582518f,
+            0.98582518f,
+            0.87685609f,
+            0.68134463f,
+            0.43869004f
+        };
+
+        for (int i = 0; i < N; ++i)
+            ASSERT_FLOAT_EQ(kaiser8[i], kaiserWindow.at<float>(i));
+    }
+
+    TEST(xphoto_DenoisingBm3dTransforms, regression_2D_generic)
+    {
+        const int templateWindowSize = 8;
+        const int templateWindowSizeSq = templateWindowSize * templateWindowSize;
+
+        uchar src[templateWindowSizeSq];
+        short dst[templateWindowSizeSq];
+        short dstSec[templateWindowSizeSq];
+
+        // Initialize array
+        for (uchar i = 0; i < templateWindowSizeSq; ++i)
+            src[i] = (i % 10) * 10;
+
+        // Use tailored transforms
+        HaarTransform<uchar, short>::RegisterTransforms2D(templateWindowSize);
+        HaarTransform<uchar, short>::forwardTransform2D(src, dst, templateWindowSize, templateWindowSize);
+        HaarTransform<uchar, short>::inverseTransform2D(dst, templateWindowSize);
+
+        // Use generic transforms
+        HaarTransform2D::ForwardTransformXxX<uchar, short, templateWindowSize>(src, dstSec, templateWindowSize, templateWindowSize);
+        HaarTransform2D::InverseTransformXxX<short, templateWindowSize>(dstSec, templateWindowSize);
+
+        for (unsigned i = 0; i < templateWindowSizeSq; ++i)
+            ASSERT_EQ(dst[i], dstSec[i]);
+    }
+
+    TEST(xphoto_DenoisingBm3dTransforms, regression_2D_4x4)
+    {
+        const int templateWindowSize = 4;
+        const int templateWindowSizeSq = templateWindowSize * templateWindowSize;
+
+        uchar src[templateWindowSizeSq];
+        short dst[templateWindowSizeSq];
+
+        // Initialize array
+        for (uchar i = 0; i < templateWindowSizeSq; ++i)
+        {
+            src[i] = i;
+        }
+
+        HaarTransform2D::ForwardTransform4x4(src, dst, templateWindowSize, templateWindowSize);
+        HaarTransform2D::InverseTransform4x4(dst, templateWindowSize);
+
+        for (uchar i = 0; i < templateWindowSizeSq; ++i)
+            ASSERT_EQ(static_cast<short>(src[i]), dst[i]);
+    }
+
+    TEST(xphoto_DenoisingBm3dTransforms, regression_2D_8x8)
+    {
+        const int templateWindowSize = 8;
+        const int templateWindowSizeSq = templateWindowSize * templateWindowSize;
+
+        uchar src[templateWindowSizeSq];
+        short dst[templateWindowSizeSq];
+
+        // Initialize array
+        for (uchar i = 0; i < templateWindowSizeSq; ++i)
+        {
+            src[i] = i;
+        }
+
+        HaarTransform2D::ForwardTransform8x8(src, dst, templateWindowSize, templateWindowSize);
+        HaarTransform2D::InverseTransform8x8(dst, templateWindowSize);
+
+        for (uchar i = 0; i < templateWindowSizeSq; ++i)
+            ASSERT_EQ(static_cast<short>(src[i]), dst[i]);
+    }
+
+    template <typename T, typename DT, typename CT>
+    static void Test1dTransform(
+        T *thrMap,
+        int groupSize,
+        int templateWindowSizeSq,
+        BlockMatch<T, DT, CT> *bm,
+        BlockMatch<T, DT, CT> *bmOrig,
+        int expectedNonZeroCount = -1)
+    {
+        if (expectedNonZeroCount < 0)
+            expectedNonZeroCount = groupSize * templateWindowSizeSq;
+
+        // Test group size
+        short sumNonZero = 0;
+        T *thrMapPtr1D = thrMap + (groupSize - 1) * templateWindowSizeSq;
+        for (int n = 0; n < templateWindowSizeSq; n++)
+        {
+            switch (groupSize)
+            {
+            case 16:
+                HaarTransform1D::ForwardTransform16(bm, n);
+                sumNonZero += HardThreshold<16>(bm, n, thrMapPtr1D);
+                HaarTransform1D::InverseTransform16(bm, n);
+                break;
+            case 8:
+                HaarTransform1D::ForwardTransform8(bm, n);
+                sumNonZero += HardThreshold<8>(bm, n, thrMapPtr1D);
+                HaarTransform1D::InverseTransform8(bm, n);
+                break;
+            case 4:
+                HaarTransform1D::ForwardTransform4(bm, n);
+                sumNonZero += HardThreshold<4>(bm, n, thrMapPtr1D);
+                HaarTransform1D::InverseTransform4(bm, n);
+                break;
+            case 2:
+                HaarTransform1D::ForwardTransform2(bm, n);
+                sumNonZero += HardThreshold<2>(bm, n, thrMapPtr1D);
+                HaarTransform1D::InverseTransform2(bm, n);
+                break;
+            default:
+                HaarTransform1D::ForwardTransformN(bm, n, groupSize);
+                sumNonZero += HardThreshold(bm, n, thrMapPtr1D, groupSize);
+                HaarTransform1D::InverseTransformN(bm, n, groupSize);
+            }
+        }
+
+        // Assert transform
+        if (expectedNonZeroCount == groupSize * templateWindowSizeSq)
+        {
+            for (int i = 0; i < groupSize; ++i)
+                for (int j = 0; j < templateWindowSizeSq; ++j)
+                    ASSERT_EQ(bm[i][j], bmOrig[i][j]);
+        }
+
+        // Assert shrinkage
+        ASSERT_EQ(sumNonZero, expectedNonZeroCount);
+    }
+
+    TEST(xphoto_DenoisingBm3dTransforms, regression_1D_transform)
+    {
+        const int templateWindowSize = 4;
+        const int templateWindowSizeSq = templateWindowSize * templateWindowSize;
+        const int searchWindowSize = 16;
+        const int searchWindowSizeSq = searchWindowSize * searchWindowSize;
+        const float h = 10;
+        int maxGroupSize = 64;
+
+        // Precompute separate maps for transform and shrinkage verification
+        short *thrMapTransform = NULL;
+        short *thrMapShrinkage = NULL;
+        HaarTransform<short, short>::calcThresholdMap3D(thrMapTransform, 0, templateWindowSize, maxGroupSize);
+        HaarTransform<short, short>::calcThresholdMap3D(thrMapShrinkage, h, templateWindowSize, maxGroupSize);
+
+        // Generate some data
+        BlockMatch<short, int, short> *bm = new BlockMatch<short, int, short>[maxGroupSize];
+        BlockMatch<short, int, short> *bmOrig = new BlockMatch<short, int, short>[maxGroupSize];
+        for (int i = 0; i < maxGroupSize; ++i)
+        {
+            bm[i].init(templateWindowSizeSq);
+            bmOrig[i].init(templateWindowSizeSq);
+        }
+
+        for (short i = 0; i < maxGroupSize; ++i)
+        {
+            for (short j = 0; j < templateWindowSizeSq; ++j)
+            {
+                bm[i][j] = (j + 1);
+                bmOrig[i][j] = bm[i][j];
+            }
+        }
+
+        // Verify transforms
+        Test1dTransform<short, int, short>(thrMapTransform, 2, templateWindowSizeSq, bm, bmOrig);
+        Test1dTransform<short, int, short>(thrMapTransform, 4, templateWindowSizeSq, bm, bmOrig);
+        Test1dTransform<short, int, short>(thrMapTransform, 8, templateWindowSizeSq, bm, bmOrig);
+        Test1dTransform<short, int, short>(thrMapTransform, 16, templateWindowSizeSq, bm, bmOrig);
+        Test1dTransform<short, int, short>(thrMapTransform, 32, templateWindowSizeSq, bm, bmOrig);
+        Test1dTransform<short, int, short>(thrMapTransform, 64, templateWindowSizeSq, bm, bmOrig);
+
+        // Verify shrinkage
+        Test1dTransform<short, int, short>(thrMapShrinkage, 2, templateWindowSizeSq, bm, bmOrig, 6);
+        Test1dTransform<short, int, short>(thrMapShrinkage, 4, templateWindowSizeSq, bm, bmOrig, 6);
+        Test1dTransform<short, int, short>(thrMapShrinkage, 8, templateWindowSizeSq, bm, bmOrig, 6);
+        Test1dTransform<short, int, short>(thrMapShrinkage, 16, templateWindowSizeSq, bm, bmOrig, 6);
+        Test1dTransform<short, int, short>(thrMapShrinkage, 32, templateWindowSizeSq, bm, bmOrig, 6);
+        Test1dTransform<short, int, short>(thrMapShrinkage, 64, templateWindowSizeSq, bm, bmOrig, 14);
+    }
+
+    const float sqrt2 = std::sqrt(2.0f);
+
+    TEST(xphoto_DenoisingBm3dTransforms, regression_1D_generate)
+    {
+        const int numberOfElements = 8;
+        const int arrSize = (numberOfElements << 1) - 1;
+        float *thrMap1D = NULL;
+        HaarTransform<short, short>::calcThresholdMap1D(thrMap1D, numberOfElements);
+
+        // Expected array
+        const float kThrMap1D[arrSize] = {
+            1.0f,  // 1 element
+            sqrt2 / 2.0f,    sqrt2, // 2 elements
+            0.5f,            1.0f,            sqrt2,       sqrt2,  // 4 elements
+            sqrt2 / 4.0f,    sqrt2 / 2.0f,    1.0f,        1.0f,  sqrt2, sqrt2, sqrt2, sqrt2  // 8 elements
+        };
+
+        for (int j = 0; j < arrSize; ++j)
+            ASSERT_EQ(thrMap1D[j], kThrMap1D[j]);
+
+        delete[] thrMap1D;
+    }
+
+    TEST(xphoto_DenoisingBm3dTransforms, regression_2D_generate_4x4)
+    {
+        const int templateWindowSize = 4;
+        float *thrMap2D = NULL;
+        HaarTransform<short, short>::calcThresholdMap2D(thrMap2D, templateWindowSize);
+
+        // Expected array
+        const float kThrMap4x4[templateWindowSize * templateWindowSize] = {
+            0.25f,           0.5f,       sqrt2 / 2.0f,    sqrt2 / 2.0f,
+            0.5f,            1.0f,       sqrt2,           sqrt2,
+            sqrt2 / 2.0f,    sqrt2,      2.0f,            2.0f,
+            sqrt2 / 2.0f,    sqrt2,      2.0f,            2.0f
+        };
+
+        for (int j = 0; j < templateWindowSize * templateWindowSize; ++j)
+            ASSERT_EQ(thrMap2D[j], kThrMap4x4[j]);
+
+        delete[] thrMap2D;
+    }
+
+    TEST(xphoto_DenoisingBm3dTransforms, regression_2D_generate_8x8)
+    {
+        const int templateWindowSize = 8;
+        float *thrMap2D = NULL;
+        HaarTransform<short, short>::calcThresholdMap2D(thrMap2D, templateWindowSize);
+
+        // Expected array
+        const float kThrMap8x8[templateWindowSize * templateWindowSize] = {
+            0.125f,       0.25f,        sqrt2 / 4.0f, sqrt2 / 4.0f, 0.5f,  0.5f,  0.5f,  0.5f,
+            0.25f,        0.5f,         sqrt2 / 2.0f, sqrt2 / 2.0f, 1.0f,  1.0f,  1.0f,  1.0f,
+            sqrt2 / 4.0f, sqrt2 / 2.0f, 1.0f,         1.0f,         sqrt2, sqrt2, sqrt2, sqrt2,
+            sqrt2 / 4.0f, sqrt2 / 2.0f, 1.0f,         1.0f,         sqrt2, sqrt2, sqrt2, sqrt2,
+            0.5f,         1.0f,         sqrt2,        sqrt2,        2.0f,  2.0f,  2.0f,  2.0f,
+            0.5f,         1.0f,         sqrt2,        sqrt2,        2.0f,  2.0f,  2.0f,  2.0f,
+            0.5f,         1.0f,         sqrt2,        sqrt2,        2.0f,  2.0f,  2.0f,  2.0f,
+            0.5f,         1.0f,         sqrt2,        sqrt2,        2.0f,  2.0f,  2.0f,  2.0f
+        };
+
+        for (int j = 0; j < templateWindowSize * templateWindowSize; ++j)
+            ASSERT_EQ(thrMap2D[j], kThrMap8x8[j]);
+
+        delete[] thrMap2D;
+    }
+
+    TEST(xphoto_Bm3dDenoising, powerOf2)
+    {
+        ASSERT_EQ(8, getLargestPowerOf2SmallerThan(9));
+        ASSERT_EQ(16, getLargestPowerOf2SmallerThan(21));
+        ASSERT_EQ(4, getLargestPowerOf2SmallerThan(7));
+        ASSERT_EQ(8, getLargestPowerOf2SmallerThan(8));
+        ASSERT_EQ(4, getLargestPowerOf2SmallerThan(5));
+        ASSERT_EQ(4, getLargestPowerOf2SmallerThan(4));
+        ASSERT_EQ(2, getLargestPowerOf2SmallerThan(3));
+        ASSERT_EQ(1, getLargestPowerOf2SmallerThan(1));
+        ASSERT_EQ(0, getLargestPowerOf2SmallerThan(0));
+    }
+
+#endif  // TEST_TRANSFORMS
+
+}
+
+#endif  // OPENCV_ENABLE_NONFREE
\ No newline at end of file