Merge pull request #52 from ludv1x/edgefilter

Edgefilter

modules/ximgproc/CMakeLists.txt

@@ -0,0 +1,4 @@
+set(the_description "Extended image processing module. It includes edge-aware filters and etc.")
+ocv_define_module(ximgproc opencv_imgproc opencv_core opencv_highgui)
+
+target_link_libraries(opencv_ximgproc)

modules/ximgproc/doc/edge_aware_filters.rst

@@ -0,0 +1,259 @@
+.. highlight:: cpp
+
+Domain Transform filter
+====================================
+
+This section describes interface for Domain Transform filter.
+For more details about this filter see [Gastal11]_ and References_.
+
+DTFilter
+------------------------------------
+.. ocv:class:: DTFilter : public Algorithm
+
+Interface for realizations of Domain Transform filter.
+
+createDTFilter
+------------------------------------
+Factory method, create instance of :ocv:class:`DTFilter` and produce initialization routines.
+
+.. ocv:function:: Ptr<DTFilter> createDTFilter(InputArray guide, double sigmaSpatial, double sigmaColor, int mode = DTF_NC, int numIters = 3)
+
+.. ocv:pyfunction:: cv2.createDTFilter(guide, sigmaSpatial, sigmaColor[, mode[, numIters]]) -> instance
+
+    :param guide: guided image (used to build transformed distance, which describes edge structure of guided image).
+    :param sigmaSpatial: :math:`{\sigma}_H` parameter in the original article, it's similar to the sigma in the coordinate space into :ocv:func:`bilateralFilter`.
+    :param sigmaColor: :math:`{\sigma}_r` parameter in the original article, it's similar to the sigma in the color space into :ocv:func:`bilateralFilter`.
+    :param mode: one form three modes ``DTF_NC``, ``DTF_RF`` and ``DTF_IC`` which corresponds to three modes for filtering 2D signals in the article.
+    :param numIters: optional number of iterations used for filtering, 3 is quite enough.
+
+For more details about Domain Transform filter parameters, see the original article [Gastal11]_ and `Domain Transform filter homepage <http://www.inf.ufrgs.br/~eslgastal/DomainTransform/>`_.
+
+DTFilter::filter
+------------------------------------
+Produce domain transform filtering operation on source image.
+
+.. ocv:function:: void DTFilter::filter(InputArray src, OutputArray dst, int dDepth = -1)
+
+.. ocv:pyfunction:: cv2.DTFilter.filter(src, dst[, dDepth]) -> None
+
+    :param src: filtering image with unsigned 8-bit or floating-point 32-bit depth and up to 4 channels.
+    :param dst: destination image.
+    :param dDepth: optional depth of the output image. ``dDepth`` can be set to -1, which will be equivalent to ``src.depth()``.
+    
+dtFilter
+------------------------------------
+Simple one-line Domain Transform filter call.
+If you have multiple images to filter with the same guided image then use :ocv:class:`DTFilter` interface to avoid extra computations on initialization stage.
+
+.. ocv:function:: void dtFilter(InputArray guide, InputArray src, OutputArray dst, double sigmaSpatial, double sigmaColor, int mode = DTF_NC, int numIters = 3)
+
+.. ocv:pyfunction:: cv2.dtFilter(guide, src, sigmaSpatial, sigmaColor[, mode[, numIters]]) -> None
+
+    :param guide: guided image (also called as joint image) with unsigned 8-bit or floating-point 32-bit depth and up to 4 channels.
+    :param src: filtering image with unsigned 8-bit or floating-point 32-bit depth and up to 4 channels.
+    :param sigmaSpatial: :math:`{\sigma}_H` parameter in the original article, it's similar to the sigma in the coordinate space into :ocv:func:`bilateralFilter`.
+    :param sigmaColor: :math:`{\sigma}_r` parameter in the original article, it's similar to the sigma in the color space into :ocv:func:`bilateralFilter`.
+    :param mode: one form three modes ``DTF_NC``, ``DTF_RF`` and ``DTF_IC`` which corresponds to three modes for filtering 2D signals in the article.
+    :param numIters: optional number of iterations used for filtering, 3 is quite enough.
+    
+.. seealso:: :ocv:func:`bilateralFilter`, :ocv:func:`guidedFilter`, :ocv:func:`amFilter`
+
+Guided Filter
+====================================
+
+This section describes interface for Guided Filter.
+For more details about this filter see [Kaiming10]_ and References_.
+
+GuidedFilter
+------------------------------------
+.. ocv:class:: GuidedFilter : public Algorithm
+
+Interface for realizations of Guided Filter.
+
+createGuidedFilter
+------------------------------------
+Factory method, create instance of :ocv:class:`GuidedFilter` and produce initialization routines.
+
+.. ocv:function:: Ptr<GuidedFilter> createGuidedFilter(InputArray guide, int radius, double eps)
+
+.. ocv:pyfunction:: cv2.createGuidedFilter(guide, radius, eps) -> instance
+
+    :param guide: guided image (or array of images) with up to 3 channels, if it have more then 3 channels then only first 3 channels will be used.
+    :param radius: radius of Guided Filter.
+    :param eps: regularization term of Guided Filter. :math:`{eps}^2` is similar to the sigma in the color space into :ocv:func:`bilateralFilter`.
+
+For more details about Guided Filter parameters, see the original article [Kaiming10]_.
+
+GuidedFilter::filter
+------------------------------------
+Apply Guided Filter to the filtering image.
+
+.. ocv:function:: void GuidedFilter::filter(InputArray src, OutputArray dst, int dDepth = -1)
+
+.. ocv:pyfunction:: cv2.GuidedFilter.filter(src, dst[, dDepth]) -> None
+
+    :param src: filtering image with any numbers of channels.
+    :param dst: output image.
+    :param dDepth: optional depth of the output image. ``dDepth`` can be set to -1, which will be equivalent to ``src.depth()``.
+    
+guidedFilter
+------------------------------------
+Simple one-line Guided Filter call.
+If you have multiple images to filter with the same guided image then use :ocv:class:`GuidedFilter` interface to avoid extra computations on initialization stage.
+
+.. ocv:function:: void guidedFilter(InputArray guide, InputArray src, OutputArray dst, int radius, double eps, int dDepth = -1)
+
+.. ocv:pyfunction:: cv2.guidedFilter(guide, src, dst, radius, eps, [, dDepth]) -> None
+
+    :param guide: guided image (or array of images) with up to 3 channels, if it have more then 3 channels then only first 3 channels will be used.
+    :param src: filtering image with any numbers of channels.
+    :param dst: output image.
+    :param radius: radius of Guided Filter.
+    :param eps: regularization term of Guided Filter. :math:`{eps}^2` is similar to the sigma in the color space into :ocv:func:`bilateralFilter`.
+    :param dDepth: optional depth of the output image.
+    
+.. seealso:: :ocv:func:`bilateralFilter`, :ocv:func:`dtFilter`, :ocv:func:`amFilter`
+
+Adaptive Manifold Filter
+====================================
+
+This section describes interface for Adaptive Manifold Filter.
+
+For more details about this filter see [Gastal12]_ and References_.
+
+AdaptiveManifoldFilter
+------------------------------------
+.. ocv:class:: AdaptiveManifoldFilter : public Algorithm
+    
+    Interface for Adaptive Manifold Filter realizations.
+    
+    Below listed optional parameters which may be set up with :ocv:func:`Algorithm::set` function.
+    
+    .. ocv:member:: double sigma_s = 16.0
+    
+        Spatial standard deviation.
+        
+    .. ocv:member:: double sigma_r = 0.2
+    
+        Color space standard deviation.
+        
+    .. ocv:member:: int tree_height = -1
+    
+        Height of the manifold tree (default = -1 : automatically computed).
+    
+    .. ocv:member:: int num_pca_iterations = 1
+    
+        Number of iterations to computed the eigenvector.
+    
+    .. ocv:member:: bool adjust_outliers = false
+    
+        Specify adjust outliers using Eq. 9 or not.
+        
+    .. ocv:member:: bool use_RNG = true
+    
+        Specify use random number generator to compute eigenvector or not.
+
+createAMFilter
+------------------------------------
+Factory method, create instance of :ocv:class:`AdaptiveManifoldFilter` and produce some initialization routines.
+
+.. ocv:function:: Ptr<AdaptiveManifoldFilter> createAMFilter(double sigma_s, double sigma_r, bool adjust_outliers = false)
+
+.. ocv:pyfunction:: cv2.createAMFilter(sigma_s, sigma_r, adjust_outliers) -> instance
+
+    :param sigma_s: spatial standard deviation.
+    :param sigma_r: color space standard deviation, it is similar to the sigma in the color space into :ocv:func:`bilateralFilter`.
+    :param adjust_outliers: optional, specify perform outliers adjust operation or not, (Eq. 9) in the original paper.
+
+For more details about Adaptive Manifold Filter parameters, see the original article [Gastal12]_.
+
+.. note::
+    Joint images with `CV_8U` and `CV_16U` depth converted to images with `CV_32F` depth and [0; 1] color range before processing.
+    Hence color space sigma `sigma_r` must be in [0; 1] range, unlike same sigmas in :ocv:func:`bilateralFilter` and :ocv:func:`dtFilter` functions.
+
+AdaptiveManifoldFilter::filter
+------------------------------------
+Apply high-dimensional filtering using adaptive manifolds.
+
+.. ocv:function:: void AdaptiveManifoldFilter::filter(InputArray src, OutputArray dst, InputArray joint = noArray())
+
+.. ocv:pyfunction:: cv2.AdaptiveManifoldFilter.filter(src, dst[, joint]) -> None
+
+    :param src: filtering image with any numbers of channels.
+    :param dst: output image.
+    :param joint: optional joint (also called as guided) image with any numbers of channels.
+    
+amFilter
+------------------------------------
+Simple one-line Adaptive Manifold Filter call.
+
+.. ocv:function:: void amFilter(InputArray joint, InputArray src, OutputArray dst, double sigma_s, double sigma_r, bool adjust_outliers = false)
+
+.. ocv:pyfunction:: cv2.amFilter(joint, src, dst, sigma_s, sigma_r, [, adjust_outliers]) -> None
+
+    :param joint: joint (also  called as guided) image or array of images with any numbers of channels.
+    :param src: filtering image with any numbers of channels.
+    :param dst: output image.
+    :param sigma_s: spatial standard deviation.
+    :param sigma_r: color space standard deviation, it is similar to the sigma in the color space into :ocv:func:`bilateralFilter`.
+    :param adjust_outliers: optional, specify perform outliers adjust operation or not, (Eq. 9) in the original paper.
+    
+.. note::
+    Joint images with `CV_8U` and `CV_16U` depth converted to images with `CV_32F` depth and [0; 1] color range before processing.
+    Hence color space sigma `sigma_r` must be in [0; 1] range, unlike same sigmas in :ocv:func:`bilateralFilter` and :ocv:func:`dtFilter` functions.
+    
+.. seealso:: :ocv:func:`bilateralFilter`, :ocv:func:`dtFilter`, :ocv:func:`guidedFilter`
+
+Joint Bilateral Filter
+====================================
+
+jointBilateralFilter
+------------------------------------
+Applies the joint bilateral filter to an image.
+
+.. ocv:function:: void jointBilateralFilter(InputArray joint, InputArray src, OutputArray dst, int d, double sigmaColor, double sigmaSpace, int borderType = BORDER_DEFAULT)
+
+.. ocv:pyfunction:: cv2.jointBilateralFilter(joint, src, dst, d, sigmaColor, sigmaSpace, [, borderType]) -> None
+
+    :param joint: Joint 8-bit or floating-point, 1-channel or 3-channel image.
+
+    :param src: Source 8-bit or floating-point, 1-channel or 3-channel image with the same depth as joint image.
+
+    :param dst: Destination image of the same size and type as  ``src`` .
+
+    :param d: Diameter of each pixel neighborhood that is used during filtering. If it is non-positive, it is computed from  ``sigmaSpace`` .
+
+    :param sigmaColor: Filter sigma in the color space. A larger value of the parameter means that farther colors within the pixel neighborhood (see  ``sigmaSpace`` ) will be mixed together, resulting in larger areas of semi-equal color.
+
+    :param sigmaSpace: Filter sigma in the coordinate space. A larger value of the parameter means that farther pixels will influence each other as long as their colors are close enough (see  ``sigmaColor`` ). When  ``d>0`` , it specifies the neighborhood size regardless of  ``sigmaSpace`` . Otherwise,  ``d``  is proportional to  ``sigmaSpace`` .
+
+.. note:: :ocv:func:`bilateralFilter` and :ocv:func:`jointBilateralFilter` use L1 norm to compute difference between colors.
+    
+.. seealso:: :ocv:func:`bilateralFilter`, :ocv:func:`amFilter`
+
+References
+==========
+
+..  [Gastal11] E. Gastal and M. Oliveira, "Domain Transform for Edge-Aware Image and Video Processing", Proceedings of SIGGRAPH, 2011, vol. 30, pp. 69:1 - 69:12.
+
+    The paper is available `online <http://www.inf.ufrgs.br/~eslgastal/DomainTransform/>`__.
+
+
+..  [Gastal12] E. Gastal and M. Oliveira, "Adaptive manifolds for real-time high-dimensional filtering," Proceedings of SIGGRAPH, 2012, vol. 31, pp. 33:1 - 33:13.
+
+    The paper is available `online <http://inf.ufrgs.br/~eslgastal/AdaptiveManifolds/>`__.
+
+    
+..  [Kaiming10] Kaiming He et. al., "Guided Image Filtering," ECCV 2010, pp. 1 - 14.
+
+    The paper is available `online <http://research.microsoft.com/en-us/um/people/kahe/eccv10/>`__.
+    
+    
+.. [Tomasi98] Carlo Tomasi and Roberto Manduchi, “Bilateral filtering for gray and color images,” in Computer Vision, 1998. Sixth International Conference on . IEEE, 1998, pp. 839– 846.
+
+    The paper is available `online <https://www.cs.duke.edu/~tomasi/papers/tomasi/tomasiIccv98.pdf>`__.
+    
+
+..  [Ziyang13] Ziyang Ma et al., "Constant Time Weighted Median Filtering for Stereo Matching and Beyond," ICCV, 2013, pp. 49 - 56.
+
+    The paper is available `online <http://www.cv-foundation.org/openaccess/content_iccv_2013/papers/Ma_Constant_Time_Weighted_2013_ICCV_paper.pdf>`__.

modules/ximgproc/doc/ximgproc.rst

@@ -0,0 +1,10 @@
+********************************************
+ximgproc. Extended image processing module.
+********************************************
+
+.. highlight:: cpp
+
+.. toctree::
+    :maxdepth: 2
+
+    edge_aware_filters

modules/ximgproc/include/opencv2/ximgproc.hpp

@@ -0,0 +1,41 @@
+/*
+ * Software License Agreement (BSD License)
+ *
+ *  Copyright (c) 2009, Willow Garage, Inc.
+ *  All rights reserved.
+ *
+ *  Redistribution and use in source and binary forms, with or without
+ *  modification, are permitted provided that the following conditions
+ *  are met:
+ *
+ *   * Redistributions of source code must retain the above copyright
+ *     notice, this list of conditions and the following disclaimer.
+ *   * Redistributions 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.
+ *   * Neither the name of Willow Garage, Inc. nor the names of its
+ *     contributors may 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
+ *  COPYRIGHT OWNER 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.
+ *
+ */
+
+#ifndef __OPENCV_XIMGPROC_HPP__
+#define __OPENCV_XIMGPROC_HPP__
+
+#include "ximgproc/edge_filter.hpp"
+
+#endif

modules/ximgproc/include/opencv2/ximgproc/edge_filter.hpp

@@ -0,0 +1,127 @@
+/*
+ * Software License Agreement (BSD License)
+ *
+ *  Copyright (c) 2009, Willow Garage, Inc.
+ *  All rights reserved.
+ *
+ *  Redistribution and use in source and binary forms, with or without
+ *  modification, are permitted provided that the following conditions
+ *  are met:
+ *
+ *   * Redistributions of source code must retain the above copyright
+ *     notice, this list of conditions and the following disclaimer.
+ *   * Redistributions 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.
+ *   * Neither the name of Willow Garage, Inc. nor the names of its
+ *     contributors may 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
+ *  COPYRIGHT OWNER 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.
+ *
+ */
+
+#ifndef __OPENCV_EDGEFILTER_HPP__
+#define __OPENCV_EDGEFILTER_HPP__
+#ifdef __cplusplus
+
+#include <opencv2/core.hpp>
+
+namespace cv
+{
+namespace ximgproc
+{
+
+enum EdgeAwareFiltersList
+{
+    DTF_NC,
+    DTF_IC,
+    DTF_RF,
+
+    GUIDED_FILTER,
+    AM_FILTER
+};
+
+
+/*Interface for DT filters*/
+class CV_EXPORTS DTFilter : public Algorithm
+{
+public:
+
+    virtual void filter(InputArray src, OutputArray dst, int dDepth = -1) = 0;
+};
+
+typedef Ptr<DTFilter> DTFilterPtr;
+
+/*Fabric function for DT filters*/
+CV_EXPORTS
+Ptr<DTFilter> createDTFilter(InputArray guide, double sigmaSpatial, double sigmaColor, int mode = DTF_NC, int numIters = 3);
+
+/*One-line DT filter call*/
+CV_EXPORTS
+void dtFilter(InputArray guide, InputArray src, OutputArray dst, double sigmaSpatial, double sigmaColor, int mode = DTF_NC, int numIters = 3);
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+/*Interface for Guided Filter*/
+class CV_EXPORTS GuidedFilter : public Algorithm
+{
+public:
+
+    virtual void filter(InputArray src, OutputArray dst, int dDepth = -1) = 0;
+};
+
+/*Fabric function for Guided Filter*/
+CV_EXPORTS Ptr<GuidedFilter> createGuidedFilter(InputArray guide, int radius, double eps);
+
+/*One-line Guided Filter call*/
+CV_EXPORTS void guidedFilter(InputArray guide, InputArray src, OutputArray dst, int radius, double eps, int dDepth = -1);
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+class CV_EXPORTS AdaptiveManifoldFilter : public Algorithm
+{
+public:
+    /**
+        * @brief Apply High-dimensional filtering using adaptive manifolds
+        * @param src       Input image to be filtered.
+        * @param dst       Adaptive-manifold filter response.
+        * @param src_joint Image for joint filtering (optional).
+        */
+    virtual void filter(InputArray src, OutputArray dst, InputArray joint = noArray()) = 0;
+
+    virtual void collectGarbage() = 0;
+
+    static Ptr<AdaptiveManifoldFilter> create();
+};
+
+//Fabric function for AM filter algorithm
+CV_EXPORTS Ptr<AdaptiveManifoldFilter> createAMFilter(double sigma_s, double sigma_r, bool adjust_outliers = false);
+
+//One-line Adaptive Manifold filter call
+CV_EXPORTS void amFilter(InputArray joint, InputArray src, OutputArray dst, double sigma_s, double sigma_r, bool adjust_outliers = false);
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+CV_EXPORTS
+void jointBilateralFilter(InputArray joint, InputArray src, OutputArray dst, int d, double sigmaColor, double sigmaSpace, int borderType = BORDER_DEFAULT);
+
+}
+}
+#endif
+#endif

modules/ximgproc/perf/perf_adaptive_manifold.cpp

@@ -0,0 +1,57 @@
+#include "perf_precomp.hpp"
+
+namespace cvtest
+{
+
+using std::tr1::tuple;
+using std::tr1::get;
+using namespace perf;
+using namespace testing;
+using namespace cv;
+using namespace cv::ximgproc;
+
+
+typedef tuple<bool, Size, int, int, MatType> AMPerfTestParam;
+typedef TestBaseWithParam<AMPerfTestParam> AdaptiveManifoldPerfTest;
+    
+PERF_TEST_P( AdaptiveManifoldPerfTest, perf,
+    Combine(
+    Values(true, false),            //adjust_outliers flag
+    Values(sz1080p, sz720p),        //size
+    Values(1, 3, 8),                //joint channels num
+    Values(1, 3),                   //source channels num
+    Values(CV_8U, CV_32F)           //source and joint depth
+    )
+)
+{
+    AMPerfTestParam params = GetParam();
+    bool adjustOutliers = get<0>(params);
+    Size sz             = get<1>(params);
+    int jointCnNum      = get<2>(params);
+    int srcCnNum        = get<3>(params);
+    int depth           = get<4>(params);
+
+    Mat joint(sz, CV_MAKE_TYPE(depth, jointCnNum));
+    Mat src(sz, CV_MAKE_TYPE(depth, srcCnNum));
+    Mat dst(sz, CV_MAKE_TYPE(depth, srcCnNum));
+    
+    cv::setNumThreads(cv::getNumberOfCPUs());
+
+    declare.in(joint, src, WARMUP_RNG).out(dst).tbb_threads(cv::getNumberOfCPUs());
+
+    double sigma_s = 16;
+    double sigma_r = 0.5;
+    TEST_CYCLE_N(3)
+    {
+        Mat res;
+        amFilter(joint, src, res, sigma_s, sigma_r, adjustOutliers);
+
+        //at 5th cycle sigma_s will be five times more and tree depth will be 5
+        sigma_s *= 1.38;
+        sigma_r /= 1.38;
+    }
+
+    SANITY_CHECK(dst);
+}
+
+}

modules/ximgproc/perf/perf_domain_transform.cpp

@@ -0,0 +1,53 @@
+#include "perf_precomp.hpp"
+
+namespace cvtest
+{
+
+using std::tr1::tuple;
+using std::tr1::get;
+using namespace perf;
+using namespace testing;
+using namespace cv;
+using namespace cv::ximgproc;
+
+CV_ENUM(GuideMatType, CV_8UC1, CV_8UC3, CV_32FC1, CV_32FC3) //reduced set
+CV_ENUM(SourceMatType, CV_8UC1, CV_8UC2, CV_8UC3, CV_8UC4, CV_32FC1, CV_32FC2, CV_32FC3, CV_32FC4) //full supported set
+CV_ENUM(DTFMode, DTF_NC, DTF_IC, DTF_RF)
+typedef tuple<GuideMatType, SourceMatType, Size, double, double, DTFMode> DTTestParams;
+
+typedef TestBaseWithParam<DTTestParams> DomainTransformTest;
+
+PERF_TEST_P( DomainTransformTest, perf,
+             Combine(
+                      GuideMatType::all(),
+                      SourceMatType::all(),
+                      Values(szVGA, sz720p),
+                      Values(10.0, 80.0),
+                      Values(30.0, 50.0),
+                      DTFMode::all()
+                    )
+           )
+{
+    int guideType       = get<0>(GetParam());
+    int srcType         = get<1>(GetParam());
+    Size size           = get<2>(GetParam());
+    double sigmaSpatial = get<3>(GetParam());
+    double sigmaColor   = get<4>(GetParam());
+    int dtfType         = get<5>(GetParam());
+
+    Mat guide(size, guideType);
+    Mat src(size, srcType);
+    Mat dst(size, srcType);
+
+    declare.in(guide, src, WARMUP_RNG).out(dst).tbb_threads(cv::getNumberOfCPUs());
+
+    cv::setNumThreads(cv::getNumberOfCPUs());
+    TEST_CYCLE_N(5)
+    {
+        dtFilter(guide, src, dst, sigmaSpatial, sigmaColor, dtfType);
+    }
+    
+    SANITY_CHECK(dst);
+}
+
+}

modules/ximgproc/perf/perf_guided_filter.cpp

@@ -0,0 +1,45 @@
+#include "perf_precomp.hpp"
+
+namespace cvtest
+{
+
+using std::tr1::tuple;
+using std::tr1::get;
+using namespace perf;
+using namespace testing;
+using namespace cv;
+using namespace cv::ximgproc;
+
+CV_ENUM(GuideTypes, CV_8UC1, CV_8UC2, CV_8UC3, CV_32FC1, CV_32FC2, CV_32FC3);
+CV_ENUM(SrcTypes, CV_8UC1, CV_8UC3, CV_32FC1, CV_32FC3);
+typedef tuple<GuideTypes, SrcTypes, Size> GFParams;
+    
+typedef TestBaseWithParam<GFParams> GuidedFilterPerfTest;
+    
+PERF_TEST_P( GuidedFilterPerfTest, perf, Combine(GuideTypes::all(), SrcTypes::all(), Values(sz1080p, sz2K)) )
+{
+    RNG rng(0);
+
+    GFParams params = GetParam();
+    int guideType   = get<0>(params);
+    int srcType     = get<1>(params);
+    Size sz         = get<2>(params);
+
+    Mat guide(sz, guideType);
+    Mat src(sz, srcType);
+    Mat dst(sz, srcType);
+
+    declare.in(guide, src, WARMUP_RNG).out(dst).tbb_threads(cv::getNumberOfCPUs());
+
+    cv::setNumThreads(cv::getNumberOfCPUs());
+    TEST_CYCLE_N(3)
+    {
+        int radius = rng.uniform(5, 30);
+        double eps = rng.uniform(0.1, 1e5);
+        guidedFilter(guide, src, dst, radius, eps);
+    }
+
+    SANITY_CHECK(dst);
+}
+
+}

modules/ximgproc/perf/perf_main.cpp

@@ -0,0 +1,3 @@
+#include "perf_precomp.hpp"
+
+CV_PERF_TEST_MAIN(edgefilter)

modules/ximgproc/perf/perf_precomp.hpp

@@ -0,0 +1,17 @@
+#ifdef __GNUC__
+#  pragma GCC diagnostic ignored "-Wmissing-declarations"
+#  if defined __clang__ || defined __APPLE__
+#    pragma GCC diagnostic ignored "-Wmissing-prototypes"
+#    pragma GCC diagnostic ignored "-Wextra"
+#  endif
+#endif
+
+#ifndef __OPENCV_PERF_PRECOMP_HPP__
+#define __OPENCV_PERF_PRECOMP_HPP__
+
+#include <opencv2/ts.hpp>
+#include <opencv2/imgproc.hpp>
+#include <opencv2/highgui.hpp>
+#include <opencv2/ximgproc.hpp>
+
+#endif

modules/ximgproc/perf/pref_joint_bilateral_filter.cpp

@@ -0,0 +1,49 @@
+#include "perf_precomp.hpp"
+
+namespace cvtest
+{
+
+using std::tr1::tuple;
+using std::tr1::get;
+using namespace perf;
+using namespace testing;
+using namespace cv;
+using namespace cv::ximgproc;
+
+typedef tuple<double, Size, MatType, int, int> JBFTestParam;
+typedef TestBaseWithParam<JBFTestParam> JointBilateralFilterTest;
+
+PERF_TEST_P(JointBilateralFilterTest, perf, 
+    Combine(
+    Values(2.0, 4.0, 6.0, 10.0),
+    SZ_TYPICAL,
+    Values(CV_8U, CV_32F),
+    Values(1, 3),
+    Values(1, 3))
+)
+{
+    JBFTestParam params = GetParam();
+    double sigmaS   = get<0>(params);
+    Size sz         = get<1>(params);
+    int depth       = get<2>(params);
+    int jCn         = get<3>(params);
+    int srcCn       = get<4>(params);
+
+    Mat joint(sz, CV_MAKE_TYPE(depth, jCn));
+    Mat src(sz, CV_MAKE_TYPE(depth, srcCn));
+    Mat dst(sz, src.type());
+
+    cv::setNumThreads(cv::getNumberOfCPUs());
+    declare.in(joint, src, WARMUP_RNG).out(dst).tbb_threads(cv::getNumberOfCPUs());
+
+    RNG rnd(cvRound(10*sigmaS) + sz.height + depth + jCn + srcCn);
+    double sigmaC = rnd.uniform(1.0, 255.0);
+
+    TEST_CYCLE_N(1)
+    {
+        jointBilateralFilter(joint, src, dst, 0, sigmaC, sigmaS);
+    }
+    
+    SANITY_CHECK(dst);
+}
+}

modules/ximgproc/samples/CMakeLists.txt

@@ -0,0 +1,9 @@
+cmake_minimum_required(VERSION 2.8)
+project(live_demo)
+find_package(OpenCV 3.0 REQUIRED)
+
+set(SOURCES live_demo.cpp)
+
+include_directories(${OpenCV_INCLUDE_DIRS})
+add_executable(live_demo ${SOURCES} ${HEADERS})
+target_link_libraries(live_demo ${OpenCV_LIBS})

modules/ximgproc/samples/live_demo.cpp

@@ -0,0 +1,195 @@
+#include <opencv2/core.hpp>
+#include <opencv2/core/utility.hpp>
+#include <opencv2/highgui.hpp>
+#include <opencv2/imgproc.hpp>
+#include <opencv2/ximgproc.hpp>
+using namespace cv;
+using namespace cv::ximgproc;
+
+#include <iostream>
+using namespace std;
+
+typedef void(*FilteringOperation)(const Mat& src, Mat& dst);
+//current mode (filtering operation example)
+FilteringOperation g_filterOp = NULL;
+
+//list of filtering operations
+void filterDoNothing(const Mat& frame, Mat& dst);
+void filterBlurring(const Mat& frame, Mat& dst);
+void filterStylize(const Mat& frame, Mat& dst);
+void filterDetailEnhancement(const Mat& frame8u, Mat& dst);
+
+//common sliders for every mode
+int g_sigmaColor = 25;
+int g_sigmaSpatial = 10;
+
+//for Stylizing mode
+int g_edgesGamma = 100;
+
+//for Details Enhancement mode
+int g_contrastBase = 100;
+int g_detailsLevel = 100;
+
+int g_numberOfCPUs = cv::getNumberOfCPUs();
+
+//We will use two callbacks to change parameters
+void changeModeCallback(int state, void *filter);
+void changeNumberOfCpuCallback(int count, void*);
+
+void splitScreen(const Mat& rawFrame, Mat& outputFrame, Mat& srcFrame, Mat& processedFrame);
+
+//trivial filter
+void filterDoNothing(const Mat& frame, Mat& dst)
+{
+    frame.copyTo(dst);
+}
+
+//simple edge-aware blurring
+void filterBlurring(const Mat& frame, Mat& dst)
+{
+    dtFilter(frame, frame, dst, g_sigmaSpatial, g_sigmaColor, DTF_RF);
+}
+
+//stylizing filter
+void filterStylize(const Mat& frame, Mat& dst)
+{
+    //blur frame
+    Mat filtered;
+    dtFilter(frame, frame, filtered, g_sigmaSpatial, g_sigmaColor, DTF_NC);
+
+    //compute grayscale blurred frame
+    Mat filteredGray;
+    cvtColor(filtered, filteredGray, COLOR_BGR2GRAY);
+
+    //find gradients of blurred image
+    Mat gradX, gradY;
+    Sobel(filteredGray, gradX, CV_32F, 1, 0, 3, 1.0/255);
+    Sobel(filteredGray, gradY, CV_32F, 0, 1, 3, 1.0/255);
+
+    //compute magnitude of gradient and fit it accordingly the gamma parameter
+    Mat gradMagnitude;
+    magnitude(gradX, gradY, gradMagnitude);
+    cv::pow(gradMagnitude, g_edgesGamma/100.0, gradMagnitude);
+
+    //multiply a blurred frame to the value inversely proportional to the magnitude
+    Mat multiplier = 1.0/(1.0 + gradMagnitude);
+    cvtColor(multiplier, multiplier, COLOR_GRAY2BGR);
+    multiply(filtered, multiplier, dst, 1, dst.type());
+}
+
+void filterDetailEnhancement(const Mat& frame8u, Mat& dst)
+{
+    Mat frame;
+    frame8u.convertTo(frame, CV_32F, 1.0/255);
+
+    //Decompose image to 3 Lab channels
+    Mat frameLab, frameLabCn[3];
+    cvtColor(frame, frameLab, COLOR_BGR2Lab);
+    split(frameLab, frameLabCn);
+
+    //Generate progressively smoother versions of the lightness channel
+    Mat layer0 = frameLabCn[0]; //first channel is original lightness
+    Mat layer1, layer2;
+    dtFilter(layer0, layer0, layer1, g_sigmaSpatial, g_sigmaColor, DTF_IC);
+    dtFilter(layer1, layer1, layer2, 2*g_sigmaSpatial, g_sigmaColor, DTF_IC);
+
+    //Compute detail layers
+    Mat detailLayer1 = layer0 - layer1;
+    Mat detailLayer2 = layer1 - layer2;
+
+    double cBase = g_contrastBase / 100.0;
+    double cDetails1 = g_detailsLevel / 100.0;
+    double cDetails2 = 2.0 - g_detailsLevel / 100.0;
+
+    //Generate lightness
+    double meanLigtness = mean(frameLabCn[0])[0];
+    frameLabCn[0]  = cBase*(layer2 - meanLigtness) + meanLigtness; //fit contrast of base (most blurred) layer
+    frameLabCn[0] += cDetails1*detailLayer1; //add weighted sum of detail layers to new lightness
+    frameLabCn[0] += cDetails2*detailLayer2; //
+
+    //Update new lightness
+    merge(frameLabCn, 3, frameLab);
+    cvtColor(frameLab, frame, COLOR_Lab2BGR);
+    frame.convertTo(dst, CV_8U, 255);
+}
+
+void changeModeCallback(int state, void *filter)
+{
+    if (state == 1)
+        g_filterOp = (FilteringOperation) filter;
+}
+
+void changeNumberOfCpuCallback(int count, void*)
+{
+    count = std::max(1, count);
+    cv::setNumThreads(count);
+    g_numberOfCPUs = count;
+}
+
+//divide screen on two parts: srcFrame and processed Frame
+void splitScreen(const Mat& rawFrame, Mat& outputFrame, Mat& srcFrame, Mat& processedFrame)
+{
+    int h = rawFrame.rows;
+    int w = rawFrame.cols;
+    int cn = rawFrame.channels();
+
+    outputFrame.create(h, 2 * w, CV_MAKE_TYPE(CV_8U, cn));
+    srcFrame = outputFrame(Range::all(), Range(0, w));
+    processedFrame = outputFrame(Range::all(), Range(w, 2 * w));
+    rawFrame.convertTo(srcFrame, srcFrame.type());
+}
+
+int main()
+{
+    VideoCapture cap(0);
+    if (!cap.isOpened())
+    {
+        cerr << "Capture device was not found" << endl;
+        return -1;
+    }
+
+    namedWindow("Demo");
+    displayOverlay("Demo", "Press Ctrl+P to show property window", 5000);
+    
+    //Thread trackbar
+    cv::setNumThreads(g_numberOfCPUs); //speedup filtering
+    createTrackbar("Threads", NULL, &g_numberOfCPUs, cv::getNumberOfCPUs(), changeNumberOfCpuCallback);
+
+    //Buttons to choose different modes
+    createButton("Mode Details Enhancement", changeModeCallback, (void*)filterDetailEnhancement, QT_RADIOBOX, true);
+    createButton("Mode Stylizing",           changeModeCallback, (void*)filterStylize,           QT_RADIOBOX, false);
+    createButton("Mode Blurring",            changeModeCallback, (void*)filterBlurring,          QT_RADIOBOX, false);
+    createButton("Mode DoNothing",           changeModeCallback, (void*)filterDoNothing,         QT_RADIOBOX, false);
+
+    //sliders for Details Enhancement mode
+    g_filterOp = filterDetailEnhancement; //set Details Enhancement as default filter
+    createTrackbar("Detail contrast", NULL, &g_contrastBase, 200);
+    createTrackbar("Detail level" , NULL, &g_detailsLevel, 200);
+    
+    //sliders for Stylizing mode
+    createTrackbar("Style gamma", NULL, &g_edgesGamma, 300);
+
+    //sliders for every mode
+    createTrackbar("Sigma Spatial", NULL, &g_sigmaSpatial, 200);
+    createTrackbar("Sigma Color"  , NULL, &g_sigmaColor, 200);
+
+    Mat rawFrame, outputFrame;
+    Mat srcFrame, processedFrame;
+
+    for (;;)
+    {
+        do
+        {
+            cap >> rawFrame;
+        } while (rawFrame.empty());
+
+        splitScreen(rawFrame, outputFrame, srcFrame, processedFrame);
+        g_filterOp(srcFrame, processedFrame);
+
+        imshow("Demo", outputFrame);
+
+        if (waitKey(1) == 27) break;
+    }
+
+    return 0;
+}

modules/ximgproc/src/adaptive_manifold_filter_n.cpp

@@ -0,0 +1,774 @@
+#include "precomp.hpp"
+#include "edgeaware_filters_common.hpp"
+#include <cmath>
+#include <cstring>
+#include <limits>
+
+namespace
+{
+
+using std::numeric_limits;
+using std::vector;
+using namespace cv;
+using namespace cv::ximgproc;
+using namespace cv::ximgproc::intrinsics;
+
+#ifndef SQR
+#define SQR(x) ((x)*(x))
+#endif
+
+void computeEigenVector(const Mat1f& X, const Mat1b& mask, Mat1f& dst, int num_pca_iterations, const Mat1f& rand_vec);
+
+inline double Log2(double n)
+{
+    return log(n) / log(2.0);
+}
+
+inline double floor_to_power_of_two(double r)
+{
+    return pow(2.0, floor(Log2(r)));
+}
+
+inline int computeManifoldTreeHeight(double sigma_s, double sigma_r)
+{
+    const double Hs = floor(Log2(sigma_s)) - 1.0;
+    const double Lr = 1.0 - sigma_r;
+    return max(2, static_cast<int>(ceil(Hs * Lr)));
+}
+
+static void splitChannels(InputArrayOfArrays src, vector<Mat>& dst)
+{
+    CV_Assert(src.isMat() || src.isUMat() || src.isMatVector() || src.isUMatVector());
+
+    if ( src.isMat() || src.isUMat() )
+    {
+        split(src, dst);
+    }
+    else
+    {
+        Size sz;
+        int depth, totalCnNum;
+
+        checkSameSizeAndDepth(src, sz, depth);
+        totalCnNum = getTotalNumberOfChannels(src);
+
+        dst.resize(totalCnNum);
+        vector<int> fromTo(2 * totalCnNum);
+        for (int i = 0; i < totalCnNum; i++)
+        {
+            fromTo[i * 2 + 0] = i;
+            fromTo[i * 2 + 1] = i;
+
+            dst[i].create(sz, CV_MAKE_TYPE(depth, 1));
+        }
+
+        mixChannels(src, dst, fromTo);
+    }
+}
+
+class AdaptiveManifoldFilterN : public AdaptiveManifoldFilter
+{
+public:
+    AlgorithmInfo* info() const;
+
+    AdaptiveManifoldFilterN();
+
+    void filter(InputArray src, OutputArray dst, InputArray joint);
+
+    void collectGarbage();
+
+protected:
+
+    bool adjust_outliers_;
+    double sigma_s_;
+    double sigma_r_;
+    int tree_height_;
+    int num_pca_iterations_;
+    bool useRNG;
+
+private:
+    
+    Size srcSize;
+    Size smallSize;
+    int jointCnNum;
+    int srcCnNum;
+
+    vector<Mat> jointCn;
+    vector<Mat> srcCn;
+
+    vector<Mat> etaFull;
+
+    vector<Mat> sum_w_ki_Psi_blur_;
+    Mat sum_w_ki_Psi_blur_0_;    
+    
+    Mat w_k;
+    Mat Psi_splat_0_small;
+    vector<Mat> Psi_splat_small;
+
+    Mat1f minDistToManifoldSquared;
+    
+    int curTreeHeight;
+    float sigma_r_over_sqrt_2;
+
+    RNG rnd;
+
+private: /*inline functions*/
+
+    double getNormalizer(int depth)
+    {
+        double normalizer = 1.0;
+
+        if (depth == CV_8U)
+            normalizer = 1.0 / 0xFF;
+        else if (depth == CV_16U)
+            normalizer = 1.0 / 0xFFFF;
+
+        return normalizer;
+    }
+
+    double getResizeRatio()
+    {
+        double df = min(sigma_s_ / 4.0, 256.0 * sigma_r_);
+        df = floor_to_power_of_two(df);
+        df = max(1.0, df);
+        return df;
+    }
+
+    Size getSmallSize()
+    {
+        double df = getResizeRatio();
+        return Size( cvRound(srcSize.width * (1.0/df)), cvRound(srcSize.height*(1.0/df)) ) ;
+    }
+
+    void downsample(InputArray src, OutputArray dst)
+    {
+        if (src.isMatVector())
+        {
+            vector<Mat>& srcv = *static_cast< vector<Mat>* >(src.getObj());
+            vector<Mat>& dstv = *static_cast< vector<Mat>* >(dst.getObj());
+            dstv.resize(srcv.size());
+            for (int i = 0; i < (int)srcv.size(); i++)
+                downsample(srcv[i], dstv[i]);
+        }
+        else
+        {
+            double df = getResizeRatio();
+            CV_DbgAssert(src.empty() || src.size() == srcSize);
+            resize(src, dst, Size(), 1.0 / df, 1.0 / df, INTER_LINEAR);
+            CV_DbgAssert(dst.size() == smallSize);
+        }
+    }
+
+    void upsample(InputArray src, OutputArray dst)
+    {
+        if (src.isMatVector())
+        {
+            vector<Mat>& srcv = *static_cast< vector<Mat>* >(src.getObj());
+            vector<Mat>& dstv = *static_cast< vector<Mat>* >(dst.getObj());
+            dstv.resize(srcv.size());
+            for (int i = 0; i < (int)srcv.size(); i++)
+                upsample(srcv[i], dstv[i]);
+        }
+        else
+        {
+            CV_DbgAssert(src.empty() || src.size() == smallSize);
+            resize(src, dst, srcSize, 0, 0);
+        }
+    }
+
+private:
+
+    void initBuffers(InputArray src_, InputArray joint_);
+
+    void initSrcAndJoint(InputArray src_, InputArray joint_);
+
+    void buildManifoldsAndPerformFiltering(vector<Mat>& eta, Mat1b& cluster, int treeLevel);
+
+    void gatherResult(InputArray src_, OutputArray dst_);
+
+    void compute_w_k(vector<Mat>& etak, Mat& dst, float sigma, int curTreeLevel);
+
+    void computeClusters(Mat1b& cluster, Mat1b& cluster_minus, Mat1b& cluster_plus);
+
+    void computeEta(Mat& teta, Mat1b& cluster, vector<Mat>& etaDst);
+
+
+    static void h_filter(const Mat1f& src, Mat& dst, float sigma);
+
+    static void RFFilterPass(vector<Mat>& joint, vector<Mat>& Psi_splat, Mat& Psi_splat_0, vector<Mat>& Psi_splat_dst, Mat& Psi_splat_0_dst, float ss, float sr);
+
+    static void computeDTHor(vector<Mat>& srcCn, Mat& dst, float ss, float sr);
+
+    static void computeDTVer(vector<Mat>& srcCn, Mat& dst, float ss, float sr);
+};
+
+CV_INIT_ALGORITHM(AdaptiveManifoldFilterN, "AdaptiveManifoldFilter",
+    obj.info()->addParam(obj, "sigma_s", obj.sigma_s_, false, 0, 0, "Filter spatial standard deviation");
+    obj.info()->addParam(obj, "sigma_r", obj.sigma_r_, false, 0, 0, "Filter range standard deviation");
+    obj.info()->addParam(obj, "tree_height", obj.tree_height_, false, 0, 0, "Height of the manifold tree (default = -1 : automatically computed)");
+    obj.info()->addParam(obj, "num_pca_iterations", obj.num_pca_iterations_, false, 0, 0, "Number of iterations to computed the eigenvector v1");
+    obj.info()->addParam(obj, "adjust_outliers", obj.adjust_outliers_, false, 0, 0, "Specify adjust outliers using Eq. 9 or not");
+    obj.info()->addParam(obj, "use_RNG", obj.useRNG, false, 0, 0, "Specify use random to compute eigenvector or not");)
+
+AdaptiveManifoldFilterN::AdaptiveManifoldFilterN()
+{
+    sigma_s_ = 16.0;
+    sigma_r_ = 0.2;
+    tree_height_ = -1;
+    num_pca_iterations_ = 1;
+    adjust_outliers_ = false;
+    useRNG = true;
+}
+
+void AdaptiveManifoldFilterN::initBuffers(InputArray src_, InputArray joint_)
+{
+    initSrcAndJoint(src_, joint_);
+
+    jointCn.resize(jointCnNum);
+    Psi_splat_small.resize(jointCnNum);
+    for (int i = 0; i < jointCnNum; i++)
+    {
+        //jointCn[i].create(srcSize, CV_32FC1);
+        Psi_splat_small[i].create(smallSize, CV_32FC1);
+    }
+
+    srcCn.resize(srcCnNum);
+    sum_w_ki_Psi_blur_.resize(srcCnNum);
+    for (int i = 0; i < srcCnNum; i++)
+    {
+        //srcCn[i].create(srcSize, CV_32FC1);
+        sum_w_ki_Psi_blur_[i] = Mat::zeros(srcSize, CV_32FC1);
+    }
+
+    sum_w_ki_Psi_blur_0_ = Mat::zeros(srcSize, CV_32FC1);
+    w_k.create(srcSize, CV_32FC1);
+    Psi_splat_0_small.create(smallSize, CV_32FC1);
+    
+    if (adjust_outliers_)
+        minDistToManifoldSquared.create(srcSize);
+}
+
+void AdaptiveManifoldFilterN::initSrcAndJoint(InputArray src_, InputArray joint_)
+{
+    srcSize = src_.size();
+    smallSize = getSmallSize();
+    srcCnNum = src_.channels();
+
+    split(src_, srcCn);
+    if (src_.depth() != CV_32F)
+    {
+        for (int i = 0; i < srcCnNum; i++)
+            srcCn[i].convertTo(srcCn[i], CV_32F);
+    }
+
+    if (joint_.empty() || joint_.getObj() == src_.getObj())
+    {
+        jointCnNum = srcCnNum;
+
+        if (src_.depth() == CV_32F)
+        {
+            jointCn = srcCn;
+        }
+        else
+        {
+            jointCn.resize(jointCnNum);
+            for (int i = 0; i < jointCnNum; i++)
+                srcCn[i].convertTo(jointCn[i], CV_32F, getNormalizer(src_.depth()));
+        }
+    }
+    else
+    {
+        splitChannels(joint_, jointCn);
+        
+        jointCnNum = (int)jointCn.size();
+        int jointDepth = jointCn[0].depth();
+        Size jointSize = jointCn[0].size();
+
+        CV_Assert( jointSize == srcSize && (jointDepth == CV_8U || jointDepth == CV_16U || jointDepth == CV_32F) );
+
+        if (jointDepth != CV_32F)
+        {
+            for (int i = 0; i < jointCnNum; i++)
+                jointCn[i].convertTo(jointCn[i], CV_32F, getNormalizer(jointDepth));
+        }
+    }
+}
+
+void AdaptiveManifoldFilterN::filter(InputArray src, OutputArray dst, InputArray joint)
+{
+    CV_Assert(sigma_s_ >= 1 && (sigma_r_ > 0 && sigma_r_ <= 1));
+    num_pca_iterations_ = std::max(1, num_pca_iterations_);
+
+    initBuffers(src, joint);
+
+    curTreeHeight = tree_height_ <= 0 ? computeManifoldTreeHeight(sigma_s_, sigma_r_) : tree_height_;
+
+    sigma_r_over_sqrt_2 = (float) (sigma_r_ / sqrt(2.0));
+
+    const double seedCoef = jointCn[0].at<float>(srcSize.height/2, srcSize.width/2);
+    const uint64 baseCoef = numeric_limits<uint64>::max() / 0xFFFF;
+    rnd.state = static_cast<int64>(baseCoef*seedCoef);
+    
+    Mat1b cluster0(srcSize, 0xFF);
+    vector<Mat> eta0(jointCnNum);
+    for (int i = 0; i < jointCnNum; i++)
+        h_filter(jointCn[i], eta0[i], (float)sigma_s_);
+
+    buildManifoldsAndPerformFiltering(eta0, cluster0, 1);
+
+    gatherResult(src, dst);
+}
+
+void AdaptiveManifoldFilterN::gatherResult(InputArray src_, OutputArray dst_)
+{
+    int dDepth = src_.depth();
+    vector<Mat> dstCn(srcCnNum);
+
+    if (!adjust_outliers_)
+    {
+        for (int i = 0; i < srcCnNum; i++)
+            divide(sum_w_ki_Psi_blur_[i], sum_w_ki_Psi_blur_0_, dstCn[i], 1.0, dDepth);
+
+        merge(dstCn, dst_);
+    }
+    else
+    {
+        Mat1f& alpha = minDistToManifoldSquared;
+        double sigmaMember = -0.5 / (sigma_r_*sigma_r_);
+        multiply(minDistToManifoldSquared, sigmaMember, minDistToManifoldSquared);
+        cv::exp(minDistToManifoldSquared, alpha);
+
+        for (int i = 0; i < srcCnNum; i++)
+        {
+            Mat& f = srcCn[i];
+            Mat& g = dstCn[i];
+
+            divide(sum_w_ki_Psi_blur_[i], sum_w_ki_Psi_blur_0_, g);
+
+            subtract(g, f, g);
+            multiply(alpha, g, g);
+            add(g, f, g);
+
+            g.convertTo(g, dDepth);
+        }
+
+        merge(dstCn, dst_);
+    }
+}
+
+void AdaptiveManifoldFilterN::buildManifoldsAndPerformFiltering(vector<Mat>& eta, Mat1b& cluster, int treeLevel)
+{
+    CV_DbgAssert((int)eta.size() == jointCnNum);
+
+    //splatting
+    Size etaSize = eta[0].size();
+    CV_DbgAssert(etaSize == srcSize || etaSize == smallSize);
+
+    if (etaSize == srcSize)
+    {
+        compute_w_k(eta, w_k, sigma_r_over_sqrt_2, treeLevel);
+        etaFull = eta;
+        downsample(eta, eta);
+    }
+    else
+    {
+        upsample(eta, etaFull);
+        compute_w_k(etaFull, w_k, sigma_r_over_sqrt_2, treeLevel);
+    }
+    
+    //blurring
+    Psi_splat_small.resize(srcCnNum);
+    for (int si = 0; si < srcCnNum; si++)
+    {
+        Mat tmp;
+        multiply(srcCn[si], w_k, tmp);
+        downsample(tmp, Psi_splat_small[si]);
+    }
+    downsample(w_k, Psi_splat_0_small);
+
+    vector<Mat>& Psi_splat_small_blur = Psi_splat_small;
+    Mat& Psi_splat_0_small_blur = Psi_splat_0_small;
+
+    float rf_ss = (float)(sigma_s_ / getResizeRatio());
+    float rf_sr = (float)(sigma_r_over_sqrt_2);
+    RFFilterPass(eta, Psi_splat_small, Psi_splat_0_small, Psi_splat_small_blur, Psi_splat_0_small_blur, rf_ss, rf_sr);
+
+    //slicing
+    {
+        Mat tmp;
+        for (int i = 0; i < srcCnNum; i++)
+        {
+            upsample(Psi_splat_small_blur[i], tmp);
+            multiply(tmp, w_k, tmp);
+            add(sum_w_ki_Psi_blur_[i], tmp, sum_w_ki_Psi_blur_[i]);
+        }
+        upsample(Psi_splat_0_small_blur, tmp);
+        multiply(tmp, w_k, tmp);
+        add(sum_w_ki_Psi_blur_0_, tmp, sum_w_ki_Psi_blur_0_);
+    }
+
+    //build new manifolds
+    if (treeLevel < curTreeHeight)
+    {
+        Mat1b cluster_minus, cluster_plus;
+
+        computeClusters(cluster, cluster_minus, cluster_plus);
+
+        vector<Mat> eta_minus(jointCnNum), eta_plus(jointCnNum);
+        {
+            Mat1f teta = 1.0 - w_k;
+            computeEta(teta, cluster_minus, eta_minus);
+            computeEta(teta, cluster_plus, eta_plus);
+        }
+
+        //free memory to continue deep recursion
+        eta.clear();
+        cluster.release();
+
+        buildManifoldsAndPerformFiltering(eta_minus, cluster_minus, treeLevel + 1);
+        buildManifoldsAndPerformFiltering(eta_plus, cluster_plus, treeLevel + 1);
+    }
+}
+
+void AdaptiveManifoldFilterN::collectGarbage()
+{
+    srcCn.clear();
+    jointCn.clear();
+    etaFull.clear();
+    sum_w_ki_Psi_blur_.clear();
+    Psi_splat_small.clear();
+
+    sum_w_ki_Psi_blur_0_.release();
+    w_k.release();
+    Psi_splat_0_small.release();
+    minDistToManifoldSquared.release();
+}
+
+void AdaptiveManifoldFilterN::h_filter(const Mat1f& src, Mat& dst, float sigma)
+{
+    CV_DbgAssert(src.depth() == CV_32F);
+
+    const float a = exp(-sqrt(2.0f) / sigma);
+
+    dst.create(src.size(), CV_32FC1);
+
+    for (int y = 0; y < src.rows; ++y)
+    {
+        const float* src_row = src[y];
+        float* dst_row = dst.ptr<float>(y);
+
+        dst_row[0] = src_row[0];
+        for (int x = 1; x < src.cols; ++x)
+        {
+            dst_row[x] = src_row[x] + a * (dst_row[x - 1] - src_row[x]);
+        }
+        for (int x = src.cols - 2; x >= 0; --x)
+        {
+            dst_row[x] = dst_row[x] + a * (dst_row[x + 1] - dst_row[x]);
+        }
+    }
+
+    for (int y = 1; y < src.rows; ++y)
+    {
+        float* dst_cur_row = dst.ptr<float>(y);
+        float* dst_prev_row = dst.ptr<float>(y-1);
+
+        rf_vert_row_pass(dst_cur_row, dst_prev_row, a, src.cols);
+    }
+    for (int y = src.rows - 2; y >= 0; --y)
+    {
+        float* dst_cur_row = dst.ptr<float>(y);
+        float* dst_prev_row = dst.ptr<float>(y+1);
+
+        rf_vert_row_pass(dst_cur_row, dst_prev_row, a, src.cols);
+    }
+}
+
+void AdaptiveManifoldFilterN::compute_w_k(vector<Mat>& etak, Mat& dst, float sigma, int curTreeLevel)
+{
+    CV_DbgAssert((int)etak.size() == jointCnNum);
+
+    dst.create(srcSize, CV_32FC1);
+    float argConst = -0.5f / (sigma*sigma);
+
+    for (int i = 0; i < srcSize.height; i++)
+    {
+        float *dstRow = dst.ptr<float>(i);
+
+        for (int cn = 0; cn < jointCnNum; cn++)
+        {
+            float *eta_kCnRow = etak[cn].ptr<float>(i);
+            float *jointCnRow = jointCn[cn].ptr<float>(i);
+
+            if (cn == 0)
+            {
+                sqr_dif(dstRow, eta_kCnRow, jointCnRow, srcSize.width);
+            }
+            else
+            {
+                add_sqr_dif(dstRow, eta_kCnRow, jointCnRow, srcSize.width);
+            }
+        }
+
+        if (adjust_outliers_)
+        {
+            float *minDistRow = minDistToManifoldSquared.ptr<float>(i);
+
+            if (curTreeLevel != 1)
+            {
+                min_(minDistRow, minDistRow, dstRow, srcSize.width);
+            }
+            else
+            {
+                std::memcpy(minDistRow, dstRow, srcSize.width*sizeof(float));
+            }
+        }
+
+        mul(dstRow, dstRow, argConst, srcSize.width);
+        //Exp_32f(dstRow, dstRow, srcSize.width);
+    }
+
+    cv::exp(dst, dst);
+}
+
+void AdaptiveManifoldFilterN::computeDTHor(vector<Mat>& srcCn, Mat& dst, float sigma_s, float sigma_r)
+{
+    int cnNum = (int)srcCn.size();
+    int h = srcCn[0].rows;
+    int w = srcCn[0].cols;
+
+    float sigmaRatioSqr = (float) SQR(sigma_s / sigma_r);
+    float lnAlpha       = (float) (-sqrt(2.0) / sigma_s);
+
+    dst.create(h, w-1, CV_32F);
+
+    for (int i = 0; i < h; i++)
+    {
+        float *dstRow = dst.ptr<float>(i);
+
+        for (int cn = 0; cn < cnNum; cn++)
+        {
+            float *curCnRow = srcCn[cn].ptr<float>(i);
+
+            if (cn == 0)
+                sqr_dif(dstRow, curCnRow, curCnRow + 1, w - 1);
+            else
+                add_sqr_dif(dstRow, curCnRow, curCnRow + 1, w - 1);
+        }
+        
+        mad(dstRow, dstRow, sigmaRatioSqr, 1.0f, w - 1);
+        sqrt_(dstRow, dstRow, w - 1);
+        mul(dstRow, dstRow, lnAlpha, w - 1);
+        //Exp_32f(dstRow, dstRow, w - 1);
+    }
+    
+    cv::exp(dst, dst);
+}
+
+void AdaptiveManifoldFilterN::computeDTVer(vector<Mat>& srcCn, Mat& dst, float sigma_s, float sigma_r)
+{
+    int cnNum = (int)srcCn.size();
+    int h = srcCn[0].rows;
+    int w = srcCn[0].cols;
+
+    dst.create(h-1, w, CV_32F);
+
+    float sigmaRatioSqr = (float) SQR(sigma_s / sigma_r);
+    float lnAlpha       = (float) (-sqrt(2.0) / sigma_s);
+
+    for (int i = 0; i < h-1; i++)
+    {
+        float *dstRow = dst.ptr<float>(i);
+
+        for (int cn = 0; cn < cnNum; cn++)
+        {
+            float *srcRow1 = srcCn[cn].ptr<float>(i);
+            float *srcRow2 = srcCn[cn].ptr<float>(i+1);
+
+            if (cn == 0)
+                sqr_dif(dstRow, srcRow1, srcRow2, w);
+            else
+                add_sqr_dif(dstRow, srcRow1, srcRow2, w);
+        }
+
+        mad(dstRow, dstRow, sigmaRatioSqr, 1.0f, w);
+        sqrt_(dstRow, dstRow, w);
+
+        mul(dstRow, dstRow, lnAlpha, w);
+        //Exp_32f(dstRow, dstRow, w);
+    }
+
+    cv::exp(dst, dst);
+}
+
+void AdaptiveManifoldFilterN::RFFilterPass(vector<Mat>& joint, vector<Mat>& Psi_splat, Mat& Psi_splat_0, vector<Mat>& Psi_splat_dst, Mat& Psi_splat_0_dst, float ss, float sr)
+{
+    int h = joint[0].rows;
+    int w = joint[0].cols;
+    int cnNum = (int)Psi_splat.size();
+
+    Mat adth, adtv;
+    computeDTHor(joint, adth, ss, sr);
+    computeDTVer(joint, adtv, ss, sr);
+
+    Psi_splat_0_dst.create(h, w, CV_32FC1);
+    Psi_splat_dst.resize(cnNum);
+    for (int cn = 0; cn < cnNum; cn++)
+        Psi_splat_dst[cn].create(h, w, CV_32FC1);
+
+    Ptr<DTFilter> dtf = createDTFilterRF(adth, adtv, ss, sr, 1);
+    for (int cn = 0; cn < cnNum; cn++)
+        dtf->filter(Psi_splat[cn], Psi_splat_dst[cn]);
+    dtf->filter(Psi_splat_0, Psi_splat_0_dst);
+}
+
+void AdaptiveManifoldFilterN::computeClusters(Mat1b& cluster, Mat1b& cluster_minus, Mat1b& cluster_plus)
+{
+    Mat difEtaSrc;
+    {
+        vector<Mat> eta_difCn(jointCnNum);
+        for (int i = 0; i < jointCnNum; i++)
+            subtract(jointCn[i], etaFull[i], eta_difCn[i]);
+
+        merge(eta_difCn, difEtaSrc);
+        difEtaSrc = difEtaSrc.reshape(1, (int)difEtaSrc.total());
+        CV_DbgAssert(difEtaSrc.cols == jointCnNum);
+    }
+
+    Mat1f initVec(1, jointCnNum);
+    if (useRNG)
+    {
+        rnd.fill(initVec, RNG::UNIFORM, -0.5, 0.5);
+    }
+    else
+    {
+        for (int i = 0; i < (int)initVec.total(); i++)
+            initVec(0, i) = (i % 2 == 0) ? 0.5f : -0.5f;
+    }
+
+    Mat1f eigenVec(1, jointCnNum);
+    computeEigenVector(difEtaSrc, cluster, eigenVec, num_pca_iterations_, initVec);
+
+    Mat1f difOreientation;
+    gemm(difEtaSrc, eigenVec, 1, noArray(), 0, difOreientation, GEMM_2_T);
+    difOreientation = difOreientation.reshape(1, srcSize.height);
+    CV_DbgAssert(difOreientation.size() == srcSize);
+
+    compare(difOreientation, 0, cluster_minus, CMP_LT);
+    bitwise_and(cluster_minus, cluster, cluster_minus);
+
+    compare(difOreientation, 0, cluster_plus, CMP_GE);
+    bitwise_and(cluster_plus, cluster, cluster_plus);
+}
+
+void AdaptiveManifoldFilterN::computeEta(Mat& teta, Mat1b& cluster, vector<Mat>& etaDst)
+{
+    CV_DbgAssert(teta.size() == srcSize && cluster.size() == srcSize);
+
+    Mat1f tetaMasked = Mat1f::zeros(srcSize);
+    teta.copyTo(tetaMasked, cluster);
+    
+    float sigma_s = (float)(sigma_s_ / getResizeRatio());
+
+    Mat1f tetaMaskedBlur;
+    downsample(tetaMasked, tetaMaskedBlur);
+    h_filter(tetaMaskedBlur, tetaMaskedBlur, sigma_s);
+
+    Mat mul;
+    etaDst.resize(jointCnNum);
+    for (int i = 0; i < jointCnNum; i++)
+    {
+        multiply(tetaMasked, jointCn[i], mul);
+        downsample(mul, etaDst[i]);
+        h_filter(etaDst[i], etaDst[i], sigma_s);
+        divide(etaDst[i], tetaMaskedBlur, etaDst[i]);
+    }
+}
+
+void computeEigenVector(const Mat1f& X, const Mat1b& mask, Mat1f& dst, int num_pca_iterations, const Mat1f& rand_vec)
+{
+    CV_DbgAssert( X.cols == rand_vec.cols );
+    CV_DbgAssert( X.rows == mask.size().area() );
+    CV_DbgAssert( rand_vec.rows == 1 );
+
+    dst.create(rand_vec.size());
+    rand_vec.copyTo(dst);
+
+    Mat1f t(X.size());
+
+    float* dst_row = dst[0];
+
+    for (int i = 0; i < num_pca_iterations; ++i)
+    {
+        t.setTo(Scalar::all(0));
+
+        for (int y = 0, ind = 0; y < mask.rows; ++y)
+        {
+            const uchar* mask_row = mask[y];
+
+            for (int x = 0; x < mask.cols; ++x, ++ind)
+            {
+                if (mask_row[x])
+                {
+                    const float* X_row = X[ind];
+                    float* t_row = t[ind];
+
+                    float dots = 0.0;
+                    for (int c = 0; c < X.cols; ++c)
+                        dots += dst_row[c] * X_row[c];
+
+                    for (int c = 0; c < X.cols; ++c)
+                        t_row[c] = dots * X_row[c];
+                }
+            }
+        }
+
+        dst.setTo(0.0);
+        for (int k = 0; k < X.rows; ++k)
+        {
+            const float* t_row = t[k];
+
+            for (int c = 0; c < X.cols; ++c)
+            {
+                dst_row[c] += t_row[c];
+            }
+        }
+    }
+
+    double n = norm(dst);
+    divide(dst, n, dst);
+}
+}
+
+
+namespace cv
+{
+namespace ximgproc
+{
+
+Ptr<AdaptiveManifoldFilter> AdaptiveManifoldFilter::create()
+{
+    return Ptr<AdaptiveManifoldFilter>(new AdaptiveManifoldFilterN());
+}
+
+CV_EXPORTS_W
+Ptr<AdaptiveManifoldFilter> createAMFilter(double sigma_s, double sigma_r, bool adjust_outliers)
+{
+    Ptr<AdaptiveManifoldFilter> amf(new AdaptiveManifoldFilterN());
+    
+    amf->set("sigma_s", sigma_s);
+    amf->set("sigma_r", sigma_r);
+    amf->set("adjust_outliers", adjust_outliers);
+
+    return amf;
+}
+
+CV_EXPORTS_W
+void amFilter(InputArray joint, InputArray src, OutputArray dst, double sigma_s, double sigma_r, bool adjust_outliers)
+{
+    Ptr<AdaptiveManifoldFilter> amf = createAMFilter(sigma_s, sigma_r, adjust_outliers);
+    amf->filter(src, dst, joint);
+}
+
+}
+}

modules/ximgproc/src/domain_transform.cpp

@@ -0,0 +1,24 @@
+#include "precomp.hpp"
+#include "dtfilter_cpu.hpp"
+
+namespace cv
+{
+namespace ximgproc
+{
+
+CV_EXPORTS_W
+Ptr<DTFilter> createDTFilter(InputArray guide, double sigmaSpatial, double sigmaColor, int mode, int numIters)
+{
+    return Ptr<DTFilter>(DTFilterCPU::create(guide, sigmaSpatial, sigmaColor, mode, numIters));
+}
+
+CV_EXPORTS_W
+void dtFilter(InputArray guide, InputArray src, OutputArray dst, double sigmaSpatial, double sigmaColor, int mode, int numIters)
+{
+    Ptr<DTFilterCPU> dtf = DTFilterCPU::create(guide, sigmaSpatial, sigmaColor, mode, numIters);
+    dtf->setSingleFilterCall(true);
+    dtf->filter(src, dst);
+}
+
+}
+}

modules/ximgproc/src/dtfilter_cpu.cpp

@@ -0,0 +1,177 @@
+#include "precomp.hpp"
+#include "dtfilter_cpu.hpp"
+
+namespace cv
+{
+namespace ximgproc
+{
+
+typedef Vec<uchar, 1> Vec1b;
+typedef Vec<float, 1> Vec1f;
+
+Ptr<DTFilterCPU> DTFilterCPU::create(InputArray guide, double sigmaSpatial, double sigmaColor, int mode, int numIters)
+{
+    Ptr<DTFilterCPU> dtf(new DTFilterCPU());
+    dtf->init(guide, sigmaSpatial, sigmaColor, mode, numIters);
+    return dtf;
+}
+
+Ptr<DTFilterCPU> DTFilterCPU::createRF(InputArray adistHor, InputArray adistVert, double sigmaSpatial, double sigmaColor, int numIters /*= 3*/)
+{
+    Mat adh = adistHor.getMat();
+    Mat adv = adistVert.getMat();
+    CV_Assert(adh.type() == CV_32FC1 && adv.type() == CV_32FC1 && adh.rows == adv.rows + 1 && adh.cols == adv.cols - 1);
+
+    Ptr<DTFilterCPU> dtf(new DTFilterCPU());
+    dtf->release();
+    dtf->mode = DTF_RF;
+    dtf->numIters = std::max(1, numIters);
+
+    dtf->h = adh.rows;
+    dtf->w = adh.cols + 1;
+
+    dtf->sigmaSpatial = std::max(1.0f, (float)sigmaSpatial);
+    dtf->sigmaColor = std::max(0.01f, (float)sigmaColor);
+
+    dtf->a0distHor = adh;
+    dtf->a0distVert = adv;
+
+    return dtf;
+}
+
+void DTFilterCPU::init(InputArray guide_, double sigmaSpatial_, double sigmaColor_, int mode_, int numIters_)
+{
+    Mat guide = guide_.getMat();
+
+    int cn = guide.channels();
+    int depth = guide.depth();
+
+    CV_Assert(cn <= 4);
+    CV_Assert((depth == CV_8U || depth == CV_32F) && !guide.empty());
+
+    #define CREATE_DTF(Vect) init_<Vect>(guide, sigmaSpatial_, sigmaColor_, mode_, numIters_);
+
+    if (cn == 1)
+    {
+        if (depth == CV_8U)
+            CREATE_DTF(Vec1b);
+        if (depth == CV_32F)
+            CREATE_DTF(Vec1f);
+    }
+    else if (cn == 2)
+    {
+        if (depth == CV_8U)
+            CREATE_DTF(Vec2b);
+        if (depth == CV_32F)
+            CREATE_DTF(Vec2f);
+    }
+    else if (cn == 3)
+    {
+        if (depth == CV_8U)
+            CREATE_DTF(Vec3b);
+        if (depth == CV_32F)
+            CREATE_DTF(Vec3f);
+    }
+    else if (cn == 4)
+    {
+        if (depth == CV_8U)
+            CREATE_DTF(Vec4b);
+        if (depth == CV_32F)
+            CREATE_DTF(Vec4f);
+    }
+
+    #undef CREATE_DTF
+}
+
+void DTFilterCPU::filter(InputArray src_, OutputArray dst_, int dDepth)
+{
+    Mat src = src_.getMat();
+    dst_.create(src.size(), src.type());
+    Mat& dst = dst_.getMatRef();
+
+    int cn = src.channels();
+    int depth = src.depth();
+
+    CV_Assert(cn <= 4 && (depth == CV_8U || depth == CV_32F));
+
+    if (cn == 1)
+    {
+        if (depth == CV_8U)
+            filter_<Vec1b>(src, dst, dDepth);
+        if (depth == CV_32F)
+            filter_<Vec1f>(src, dst, dDepth);
+    }
+    else if (cn == 2)
+    {
+        if (depth == CV_8U)
+            filter_<Vec2b>(src, dst, dDepth);
+        if (depth == CV_32F)
+            filter_<Vec2f>(src, dst, dDepth);
+    }
+    else if (cn == 3)
+    {
+        if (depth == CV_8U)
+            filter_<Vec3b>(src, dst, dDepth);
+        if (depth == CV_32F)
+            filter_<Vec3f>(src, dst, dDepth);
+    }
+    else if (cn == 4)
+    {
+        if (depth == CV_8U)
+            filter_<Vec4b>(src, dst, dDepth);
+        if (depth == CV_32F)
+            filter_<Vec4f>(src, dst, dDepth);
+    }
+}
+
+void DTFilterCPU::setSingleFilterCall(bool value)
+{
+    singleFilterCall = value;
+}
+
+void DTFilterCPU::release()
+{
+    if (mode == -1) return;
+
+    idistHor.release();
+    idistVert.release();
+
+    distHor.release();
+    distVert.release();
+
+    a0distHor.release();
+    a0distVert.release();
+
+    adistHor.release();
+    adistVert.release();
+}
+
+Mat DTFilterCPU::getWExtendedMat(int h, int w, int type, int brdleft /*= 0*/, int brdRight /*= 0*/, int cacheAlign /*= 0*/)
+{
+    int wrapperCols = w + brdleft + brdRight;
+    if (cacheAlign > 0)
+        wrapperCols += ((wrapperCols + cacheAlign-1) / cacheAlign)*cacheAlign;
+    Mat mat(h, wrapperCols, type);
+    return mat(Range::all(), Range(brdleft, w + brdleft));
+}
+
+
+Range DTFilterCPU::getWorkRangeByThread(const Range& itemsRange, const Range& rangeThread, int declaredNumThreads)
+{
+    if (declaredNumThreads <= 0)
+        declaredNumThreads = cv::getNumThreads();
+
+    int chunk   = itemsRange.size() / declaredNumThreads;
+    int start   = itemsRange.start + chunk * rangeThread.start;
+    int end     = itemsRange.start + ((rangeThread.end >= declaredNumThreads) ? itemsRange.size() : chunk * rangeThread.end);
+
+    return Range(start, end);
+}
+
+Range DTFilterCPU::getWorkRangeByThread(int items, const Range& rangeThread, int declaredNumThreads)
+{
+    return getWorkRangeByThread(Range(0, items), rangeThread, declaredNumThreads);
+}
+
+}
+}

modules/ximgproc/src/dtfilter_cpu.hpp

@@ -0,0 +1,258 @@
+#ifndef __OPENCV_DTFILTER_HPP__
+#define __OPENCV_DTFILTER_HPP__
+#include "precomp.hpp"
+
+#ifdef _MSC_VER
+#pragma warning(disable: 4512)
+#pragma warning(disable: 4127)
+#endif
+
+#define CV_GET_NUM_THREAD_WORKS_PROPERLY
+#undef CV_GET_NUM_THREAD_WORKS_PROPERLY
+
+namespace cv
+{
+namespace ximgproc
+{
+
+class DTFilterCPU : public DTFilter
+{
+public: /*Non-template methods*/
+
+    static Ptr<DTFilterCPU> create(InputArray guide, double sigmaSpatial, double sigmaColor, int mode = DTF_NC, int numIters = 3);
+
+    static Ptr<DTFilterCPU> createRF(InputArray adistHor, InputArray adistVert, double sigmaSpatial, double sigmaColor, int numIters = 3);
+
+    void filter(InputArray src, OutputArray dst, int dDepth = -1);
+
+    void setSingleFilterCall(bool value);
+
+public: /*Template methods*/
+
+    /*Use this static methods instead of constructor*/
+    template<typename GuideVec>
+    static DTFilterCPU* create_p_(const Mat& guide, double sigmaSpatial, double sigmaColor, int mode = DTF_NC, int numIters = 3);
+
+    template<typename GuideVec>
+    static DTFilterCPU create_(const Mat& guide, double sigmaSpatial, double sigmaColor, int mode = DTF_NC, int numIters = 3);
+
+    template<typename GuideVec>
+    void init_(Mat& guide, double sigmaSpatial, double sigmaColor, int mode = DTF_NC, int numIters = 3);
+
+    template<typename SrcVec>
+    void filter_(const Mat& src, Mat& dst, int dDepth = -1);
+
+protected: /*Typedefs declarations*/
+
+    typedef float                   IDistType;
+    typedef Vec<IDistType, 1>       IDistVec;
+
+    typedef float                   DistType;
+    typedef Vec<DistType, 1>        DistVec;
+
+    typedef float                   WorkType;
+
+public: /*Members declarations*/
+
+    int h, w, mode;
+    float sigmaSpatial, sigmaColor;
+
+    bool singleFilterCall;
+    int numFilterCalls;
+
+    Mat idistHor, idistVert;
+    Mat distHor, distVert;
+
+    Mat a0distHor, a0distVert;
+    Mat adistHor, adistVert;
+    int numIters;
+
+protected: /*Functions declarations*/
+
+    DTFilterCPU() : mode(-1), singleFilterCall(false), numFilterCalls(0) {}
+
+    void init(InputArray guide, double sigmaSpatial, double sigmaColor, int mode = DTF_NC, int numIters = 3);
+
+    void release();
+
+    template<typename GuideVec>
+    inline IDistType getTransformedDistance(const GuideVec &l, const GuideVec &r)
+    {
+        return (IDistType)(1.0f + sigmaSpatial / sigmaColor * norm1<IDistType>(l, r));
+    }
+
+    inline double getIterSigmaH(int iterNum)
+    {
+        return sigmaSpatial * std::pow(2.0, numIters - iterNum) / sqrt(std::pow(4.0, numIters) - 1);
+    }
+
+    inline IDistType getIterRadius(int iterNum)
+    {
+        return (IDistType)(3.0*getIterSigmaH(iterNum));
+    }
+
+    inline float getIterAlpha(int iterNum)
+    {
+        return (float)std::exp(-std::sqrt(2.0 / 3.0) / getIterSigmaH(iterNum));
+    }
+
+protected: /*Wrappers for parallelization*/
+
+    template <typename WorkVec>
+    struct FilterNC_horPass : public ParallelLoopBody
+    {
+        Mat &src, &idist, &dst;
+        float radius;
+
+        FilterNC_horPass(Mat& src_, Mat& idist_, Mat& dst_);
+        void operator() (const Range& range) const;
+    };
+
+    template <typename WorkVec>
+    struct FilterIC_horPass : public ParallelLoopBody
+    {
+        Mat &src, &idist, &dist, &dst, isrcBuf;
+        float radius;
+
+        FilterIC_horPass(Mat& src_, Mat& idist_, Mat& dist_, Mat& dst_);
+        void operator() (const Range& range) const;
+    };
+
+    template <typename WorkVec>
+    struct FilterRF_horPass : public ParallelLoopBody
+    {
+        Mat &res, &alphaD;
+        int iteration;
+
+        FilterRF_horPass(Mat& res_, Mat& alphaD_, int iteration_);
+        void operator() (const Range& range) const;
+        Range getRange() const { return Range(0, res.rows); }
+    };
+
+    template <typename WorkVec>
+    struct FilterRF_vertPass : public ParallelLoopBody
+    {
+        Mat &res, &alphaD;
+        int iteration;
+
+        FilterRF_vertPass(Mat& res_, Mat& alphaD_, int iteration_);
+        void operator() (const Range& range) const;
+        #ifdef CV_GET_NUM_THREAD_WORKS_PROPERLY
+        Range getRange() const { return Range(0, cv::getNumThreads()); }
+        #else
+        Range getRange() const { return Range(0, res.cols); }
+        #endif
+    };
+
+    template <typename GuideVec>
+    struct ComputeIDTHor_ParBody: public ParallelLoopBody
+    {
+        DTFilterCPU &dtf;
+        Mat &guide, &dst;
+
+        ComputeIDTHor_ParBody(DTFilterCPU& dtf_, Mat& guide_, Mat& dst_);
+        void operator() (const Range& range) const;
+        Range getRange() { return Range(0, guide.rows); }
+    };
+
+    template <typename GuideVec>
+    struct ComputeDTandIDTHor_ParBody : public ParallelLoopBody
+    {
+        DTFilterCPU &dtf;
+        Mat &guide, &dist, &idist;
+        IDistType maxRadius;
+
+        ComputeDTandIDTHor_ParBody(DTFilterCPU& dtf_, Mat& guide_, Mat& dist_, Mat& idist_);
+        void operator() (const Range& range) const;
+        Range getRange() { return Range(0, guide.rows); }
+    };
+
+    template <typename GuideVec>
+    struct ComputeA0DTHor_ParBody : public ParallelLoopBody
+    {
+        DTFilterCPU &dtf;
+        Mat &guide;
+        float lna;
+
+        ComputeA0DTHor_ParBody(DTFilterCPU& dtf_, Mat& guide_);
+        void operator() (const Range& range) const;
+        Range getRange() { return Range(0, guide.rows); }
+        ~ComputeA0DTHor_ParBody();
+    };
+
+    template <typename GuideVec>
+    struct ComputeA0DTVert_ParBody : public ParallelLoopBody
+    {
+        DTFilterCPU &dtf;
+        Mat &guide;
+        float lna;
+
+        ComputeA0DTVert_ParBody(DTFilterCPU& dtf_, Mat& guide_);
+        void operator() (const Range& range) const;
+        Range getRange() const { return Range(0, guide.rows - 1); }
+        ~ComputeA0DTVert_ParBody();
+    };
+
+protected: /*Auxiliary implementation functions*/
+
+    static Range getWorkRangeByThread(const Range& itemsRange, const Range& rangeThread, int maxThreads = 0);
+    static Range getWorkRangeByThread(int items, const Range& rangeThread, int maxThreads = 0);
+
+    template<typename SrcVec>
+    static void prepareSrcImg_IC(const Mat& src, Mat& inner, Mat& outer);
+
+    static Mat getWExtendedMat(int h, int w, int type, int brdleft = 0, int brdRight = 0, int cacheAlign = 0);
+
+    template<typename SrcVec, typename SrcWorkVec>
+    static void integrateSparseRow(const SrcVec *src, const float *dist, SrcWorkVec *dst, int cols);
+
+    template<typename SrcVec, typename SrcWorkVec>
+    static void integrateRow(const SrcVec *src, SrcWorkVec *dst, int cols);
+
+    inline static int getLeftBound(IDistType *idist, int pos, IDistType searchValue)
+    {
+        while (idist[pos] < searchValue)
+            pos++;
+        return pos;
+    }
+
+    inline static int getRightBound(IDistType *idist, int pos, IDistType searchValue)
+    {
+        while (idist[pos + 1] < searchValue)
+            pos++;
+        return pos;
+    }
+
+    template <typename T, typename T1, typename T2, int n>
+    inline static T norm1(const cv::Vec<T1, n>& v1, const cv::Vec<T2, n>& v2)
+    {
+        T sum = (T) 0;
+        for (int i = 0; i < n; i++)
+            sum += std::abs( (T)v1[i] - (T)v2[i] );
+        return sum;
+    }
+};
+
+/*One-line template wrappers for DT call*/
+
+template<typename GuideVec, typename SrcVec>
+void domainTransformFilter( const Mat_<GuideVec>& guide,
+                            const Mat_<SrcVec>& source,
+                            Mat& dst,
+                            double sigmaSpatial, double sigmaColor,
+                            int mode = DTF_NC, int numPasses = 3
+                          );
+
+template<typename GuideVec, typename SrcVec>
+void domainTransformFilter( const Mat& guide,
+                            const Mat& source,
+                            Mat& dst,
+                            double sigmaSpatial, double sigmaColor,
+                            int mode = DTF_NC, int numPasses = 3
+                          );
+}
+}
+
+#include "dtfilter_cpu.inl.hpp"
+
+#endif

modules/ximgproc/src/dtfilter_cpu.inl.hpp

@@ -0,0 +1,588 @@
+#ifndef __OPENCV_DTFILTER_INL_HPP__
+#define __OPENCV_DTFILTER_INL_HPP__
+#include "precomp.hpp"
+#include "edgeaware_filters_common.hpp"
+#include <limits>
+
+namespace cv
+{
+namespace ximgproc
+{
+
+using namespace cv::ximgproc::intrinsics;
+
+#define NC_USE_INTEGRAL_SRC
+//#undef NC_USE_INTEGRAL_SRC
+
+template<typename GuideVec>
+DTFilterCPU DTFilterCPU::create_(const Mat& guide, double sigmaSpatial, double sigmaColor, int mode, int numIters)
+{
+    DTFilterCPU dtf;
+    dtf.init_<GuideVec>(guide, sigmaSpatial, sigmaColor, mode, numIters);
+    return dtf;
+}
+
+template<typename GuideVec>
+DTFilterCPU* DTFilterCPU::create_p_(const Mat& guide, double sigmaSpatial, double sigmaColor, int mode, int numIters)
+{
+    DTFilterCPU* dtf = new DTFilterCPU();
+    dtf->init_<GuideVec>(guide, sigmaSpatial, sigmaColor, mode, numIters);
+    return dtf;
+}
+
+template<typename GuideVec>
+void DTFilterCPU::init_(Mat& guide, double sigmaSpatial_, double sigmaColor_, int mode_, int numIters_)
+{
+    CV_Assert(guide.type() == cv::DataType<GuideVec>::type);
+
+    this->release();
+
+    h = guide.rows;
+    w = guide.cols;
+
+    sigmaSpatial = std::max(1.0f, (float)sigmaSpatial_);
+    sigmaColor   = std::max(0.01f, (float)sigmaColor_);
+
+    mode = mode_;
+    numIters = std::max(1, numIters_);
+
+    if (mode == DTF_NC)
+    {
+        {
+            ComputeIDTHor_ParBody<GuideVec> horBody(*this, guide, idistHor);
+            parallel_for_(horBody.getRange(), horBody);
+        }
+        {
+            Mat guideT = guide.t();
+            ComputeIDTHor_ParBody<GuideVec> horBody(*this, guideT, idistVert);
+            parallel_for_(horBody.getRange(), horBody);
+        }
+    }
+    else if (mode == DTF_IC)
+    {
+        {
+            ComputeDTandIDTHor_ParBody<GuideVec> horBody(*this, guide, distHor, idistHor);
+            parallel_for_(horBody.getRange(), horBody);
+        }
+        {
+            Mat guideT = guide.t();
+            ComputeDTandIDTHor_ParBody<GuideVec> horBody(*this, guideT, distVert, idistVert);
+            parallel_for_(horBody.getRange(), horBody);
+        }
+    }
+    else if (mode == DTF_RF)
+    {
+        ComputeA0DTHor_ParBody<GuideVec> horBody(*this, guide);
+        ComputeA0DTVert_ParBody<GuideVec> vertBody(*this, guide);
+        parallel_for_(horBody.getRange(), horBody);
+        parallel_for_(vertBody.getRange(), vertBody);
+    }
+    else
+    {
+        CV_Error(Error::StsBadFlag, "Incorrect DT filter mode");
+    }
+}
+
+template <typename SrcVec>
+void DTFilterCPU::filter_(const Mat& src, Mat& dst, int dDepth)
+{
+    typedef typename DataType<Vec<WorkType, SrcVec::channels> >::vec_type WorkVec;
+    CV_Assert( src.type() == SrcVec::type );
+    if ( src.cols != w || src.rows != h )
+    {
+        CV_Error(Error::StsBadSize, "Size of filtering image must be equal to size of guide image");
+    }
+
+    if (singleFilterCall)
+    {
+        CV_Assert(numFilterCalls == 0);
+    }
+    numFilterCalls++;
+
+    Mat res;
+    if (dDepth == -1) dDepth = src.depth();
+    
+    //small optimization to avoid extra copying of data
+    bool useDstAsRes = (dDepth == WorkVec::depth && (mode == DTF_NC || mode == DTF_RF));
+    if (useDstAsRes)
+    {
+        dst.create(h, w, WorkVec::type);
+        res = dst;
+    }
+
+    if (mode == DTF_NC)
+    {
+        Mat resT(src.cols, src.rows, WorkVec::type);
+        src.convertTo(res, WorkVec::type);
+
+        FilterNC_horPass<WorkVec> horParBody(res, idistHor, resT);
+        FilterNC_horPass<WorkVec> vertParBody(resT, idistVert, res);
+
+        for (int iter = 1; iter <= numIters; iter++)
+        {
+            horParBody.radius = vertParBody.radius = getIterRadius(iter);
+
+            parallel_for_(Range(0, res.rows), horParBody);
+            parallel_for_(Range(0, resT.rows), vertParBody);
+        }
+    }
+    else if (mode == DTF_IC)
+    {
+        Mat resT;
+        prepareSrcImg_IC<WorkVec>(src, res, resT);
+
+        FilterIC_horPass<WorkVec> horParBody(res, idistHor, distHor, resT);
+        FilterIC_horPass<WorkVec> vertParBody(resT, idistVert, distVert, res);
+
+        for (int iter = 1; iter <= numIters; iter++)
+        {
+            horParBody.radius = vertParBody.radius = getIterRadius(iter);
+
+            parallel_for_(Range(0, res.rows), horParBody);
+            parallel_for_(Range(0, resT.rows), vertParBody);
+        }
+    }
+    else if (mode == DTF_RF)
+    {
+        src.convertTo(res, WorkVec::type);
+
+        for (int iter = 1; iter <= numIters; iter++)
+        {
+            if (!singleFilterCall && iter == 2)
+            {
+                a0distHor.copyTo(adistHor);
+                a0distVert.copyTo(adistVert);
+            }
+
+            bool useA0DT = (singleFilterCall || iter == 1);
+            Mat& a0dHor  = (useA0DT) ? a0distHor : adistHor;
+            Mat& a0dVert = (useA0DT) ? a0distVert : adistVert;
+            
+            FilterRF_horPass<WorkVec> horParBody(res, a0dHor, iter);
+            FilterRF_vertPass<WorkVec> vertParBody(res, a0dVert, iter);
+            parallel_for_(horParBody.getRange(), horParBody);
+            parallel_for_(vertParBody.getRange(), vertParBody);
+        }
+    }
+
+    if (!useDstAsRes)
+    {
+        res.convertTo(dst, dDepth);
+    }
+}
+
+template<typename SrcVec, typename SrcWorkVec>
+void DTFilterCPU::integrateRow(const SrcVec *src, SrcWorkVec *dst, int cols)
+{
+    SrcWorkVec sum = SrcWorkVec::all(0);
+    dst[0] = sum;
+
+    for (int j = 0; j < cols; j++)
+    {
+        sum += SrcWorkVec(src[j]);
+        dst[j + 1] = sum;
+    }
+}
+
+
+template<typename SrcVec, typename SrcWorkVec>
+void DTFilterCPU::integrateSparseRow(const SrcVec *src, const float *dist, SrcWorkVec *dst, int cols)
+{
+    SrcWorkVec sum = SrcWorkVec::all(0);
+    dst[0] = sum;
+
+    for (int j = 0; j < cols-1; j++)
+    {
+        sum += dist[j] * 0.5f * (SrcWorkVec(src[j]) + SrcWorkVec(src[j+1]));
+        dst[j + 1] = sum;
+    }
+}
+
+template<typename WorkVec>
+void DTFilterCPU::prepareSrcImg_IC(const Mat& src, Mat& dst, Mat& dstT)
+{
+    Mat dstOut(src.rows, src.cols + 2, WorkVec::type);
+    Mat dstOutT(src.cols, src.rows + 2, WorkVec::type);
+
+    dst = dstOut(Range::all(), Range(1, src.cols+1));
+    dstT = dstOutT(Range::all(), Range(1, src.rows+1));
+
+    src.convertTo(dst, WorkVec::type);
+
+    WorkVec *line;
+    int ri = dstOut.cols - 1;
+    for (int i = 0; i < src.rows; i++)
+    {
+        line        = dstOut.ptr<WorkVec>(i);
+        line[0]     = line[1];
+        line[ri]    = line[ri - 1];
+    }
+
+    WorkVec *topLine = dst.ptr<WorkVec>(0);
+    WorkVec *bottomLine = dst.ptr<WorkVec>(dst.rows - 1);
+    ri = dstOutT.cols - 1;
+    for (int i = 0; i < src.cols; i++)
+    {
+        line        = dstOutT.ptr<WorkVec>(i);
+        line[0]     = topLine[i];
+        line[ri]    = bottomLine[i]; 
+    }
+}
+
+
+template <typename WorkVec>
+DTFilterCPU::FilterNC_horPass<WorkVec>::FilterNC_horPass(Mat& src_, Mat& idist_, Mat& dst_)
+: src(src_), idist(idist_), dst(dst_), radius(1.0f)
+{
+    CV_DbgAssert(src.type() == WorkVec::type && dst.type() == WorkVec::type && dst.rows == src.cols && dst.cols == src.rows);
+}
+
+template <typename WorkVec>
+void DTFilterCPU::FilterNC_horPass<WorkVec>::operator()(const Range& range) const
+{
+    #ifdef NC_USE_INTEGRAL_SRC
+    std::vector<WorkVec> isrcBuf(src.cols + 1);
+    WorkVec *isrcLine = &isrcBuf[0];
+    #endif
+
+    for (int i = range.start; i < range.end; i++)
+    {
+        const WorkVec   *srcLine    = src.ptr<WorkVec>(i);
+        IDistType       *idistLine  = idist.ptr<IDistType>(i);
+        int leftBound = 0, rightBound = 0;
+        WorkVec sum;
+
+        #ifdef NC_USE_INTEGRAL_SRC
+        integrateRow(srcLine, isrcLine, src.cols);
+        #else
+        sum = srcLine[0];
+        #endif
+
+        for (int j = 0; j < src.cols; j++)
+        {
+            IDistType curVal = idistLine[j];
+            #ifdef NC_USE_INTEGRAL_SRC
+            leftBound  = getLeftBound(idistLine, leftBound, curVal - radius);
+            rightBound = getRightBound(idistLine, rightBound, curVal + radius);
+            sum = (isrcLine[rightBound + 1] - isrcLine[leftBound]);
+            #else
+            while (idistLine[leftBound] < curVal - radius)
+            {
+                sum -= srcLine[leftBound];
+                leftBound++;
+            }
+
+            while (idistLine[rightBound + 1] < curVal + radius)
+            {
+                rightBound++;
+                sum += srcLine[rightBound];
+            }
+            #endif
+
+            dst.at<WorkVec>(j, i) = sum / (float)(rightBound + 1 - leftBound);
+        }
+    }
+}
+
+template <typename WorkVec>
+DTFilterCPU::FilterIC_horPass<WorkVec>::FilterIC_horPass(Mat& src_, Mat& idist_, Mat& dist_, Mat& dst_)
+: src(src_), idist(idist_), dist(dist_), dst(dst_), radius(1.0f)
+{
+    CV_DbgAssert(src.type() == WorkVec::type && dst.type() == WorkVec::type && dst.rows == src.cols && dst.cols == src.rows);
+
+    #ifdef CV_GET_NUM_THREAD_WORKS_PROPERLY
+    isrcBuf.create(cv::getNumThreads(), src.cols + 1, WorkVec::type);
+    #else
+    isrcBuf.create(src.rows, src.cols + 1, WorkVec::type);
+    #endif
+}
+
+template <typename WorkVec>
+void DTFilterCPU::FilterIC_horPass<WorkVec>::operator()(const Range& range) const
+{
+    #ifdef CV_GET_NUM_THREAD_WORKS_PROPERLY
+    WorkVec *isrcLine = const_cast<WorkVec*>( isrcBuf.ptr<WorkVec>(cv::getThreadNum()) );
+    #else
+    WorkVec *isrcLine = const_cast<WorkVec*>( isrcBuf.ptr<WorkVec>(range.start) );
+    #endif
+
+    for (int i = range.start; i < range.end; i++)
+    {
+        WorkVec   *srcLine      = src.ptr<WorkVec>(i);
+        DistType  *distLine     = dist.ptr<DistType>(i);
+        IDistType *idistLine    = idist.ptr<IDistType>(i);
+
+        integrateSparseRow(srcLine, distLine, isrcLine, src.cols);
+
+        int leftBound = 0, rightBound = 0;
+        WorkVec sumL, sumR, sumC;
+
+        srcLine[-1] = srcLine[0];
+        srcLine[src.cols] = srcLine[src.cols - 1];
+
+        for (int j = 0; j < src.cols; j++)
+        {
+            IDistType curVal = idistLine[j];
+            IDistType valueLeft = curVal - radius;
+            IDistType valueRight = curVal + radius;
+
+            leftBound = getLeftBound(idistLine, leftBound, valueLeft);
+            rightBound = getRightBound(idistLine, rightBound, valueRight);
+
+            float areaL = idistLine[leftBound] - valueLeft;
+            float areaR = valueRight - idistLine[rightBound];
+            float dl = areaL / distLine[leftBound - 1];
+            float dr = areaR / distLine[rightBound];
+
+            sumL = 0.5f*areaL*(dl*srcLine[leftBound - 1] + (2.0f - dl)*srcLine[leftBound]);
+            sumR = 0.5f*areaR*((2.0f - dr)*srcLine[rightBound] + dr*srcLine[rightBound + 1]);
+            sumC = isrcLine[rightBound] - isrcLine[leftBound];
+
+            dst.at<WorkVec>(j, i) = (sumL + sumC + sumR) / (2.0f * radius);
+        }
+    }
+}
+
+
+template <typename WorkVec>
+DTFilterCPU::FilterRF_horPass<WorkVec>::FilterRF_horPass(Mat& res_, Mat& alphaD_, int iteration_)
+: res(res_), alphaD(alphaD_), iteration(iteration_)
+{
+    CV_DbgAssert(res.type() == WorkVec::type);
+    CV_DbgAssert(res.type() == WorkVec::type && res.size() == res.size());
+}
+
+
+template <typename WorkVec>
+void DTFilterCPU::FilterRF_horPass<WorkVec>::operator()(const Range& range) const
+{
+    for (int i = range.start; i < range.end; i++)
+    {
+        WorkVec     *dstLine = res.ptr<WorkVec>(i);
+        DistType    *adLine  = alphaD.ptr<DistType>(i);
+        int j;
+
+        if (iteration > 1)
+        {
+            for (j = res.cols - 2; j >= 0; j--)
+                adLine[j] *= adLine[j];
+        }
+
+        for (j = 1; j < res.cols; j++)
+        {
+            dstLine[j] += adLine[j-1] * (dstLine[j-1] - dstLine[j]);
+        }
+
+        for (j = res.cols - 2; j >= 0; j--)
+        {
+            dstLine[j] += adLine[j] * (dstLine[j+1] - dstLine[j]);
+        }
+    }
+}
+
+
+template <typename WorkVec>
+DTFilterCPU::FilterRF_vertPass<WorkVec>::FilterRF_vertPass(Mat& res_, Mat& alphaD_, int iteration_)
+: res(res_), alphaD(alphaD_), iteration(iteration_)
+{
+    CV_DbgAssert(res.type() == WorkVec::type);
+    CV_DbgAssert(res.type() == WorkVec::type && res.size() == res.size());
+}
+
+
+template <typename WorkVec>
+void DTFilterCPU::FilterRF_vertPass<WorkVec>::operator()(const Range& range) const
+{
+    #ifdef CV_GET_NUM_THREAD_WORKS_PROPERLY
+    Range rcols = getWorkRangeByThread(res.cols, range);
+    #else
+    Range rcols = range;
+    #endif
+
+    for (int i = 1; i < res.rows; i++)
+    {
+        WorkVec     *curRow  = res.ptr<WorkVec>(i);
+        WorkVec     *prevRow = res.ptr<WorkVec>(i - 1);
+        DistType    *adRow   = alphaD.ptr<DistType>(i - 1);
+
+        if (iteration > 1)
+        {
+            for (int j = rcols.start; j < rcols.end; j++)
+                adRow[j] *= adRow[j];
+        }
+
+        for (int j = rcols.start; j < rcols.end; j++)
+        {
+            curRow[j] += adRow[j] * (prevRow[j] - curRow[j]);
+        }
+    }
+
+    for (int i = res.rows - 2; i >= 0; i--)
+    {
+        WorkVec     *prevRow = res.ptr<WorkVec>(i + 1);
+        WorkVec     *curRow  = res.ptr<WorkVec>(i);
+        DistType    *adRow   = alphaD.ptr<DistType>(i);
+
+        for (int j = rcols.start; j < rcols.end; j++)
+        {
+            curRow[j] += adRow[j] * (prevRow[j] - curRow[j]);
+        }
+    }
+}
+
+template <typename GuideVec>
+DTFilterCPU::ComputeIDTHor_ParBody<GuideVec>::ComputeIDTHor_ParBody(DTFilterCPU& dtf_, Mat& guide_, Mat& dst_)
+: dtf(dtf_), guide(guide_), dst(dst_)
+{
+    dst.create(guide.rows, guide.cols + 1, IDistVec::type);
+}
+
+template <typename GuideVec>
+void DTFilterCPU::ComputeIDTHor_ParBody<GuideVec>::operator()(const Range& range) const
+{
+    for (int i = range.start; i < range.end; i++)
+    {
+        const GuideVec *guideLine   = guide.ptr<GuideVec>(i);
+        IDistType *idistLine        = dst.ptr<IDistType>(i);
+
+        IDistType curDist   = (IDistType)0;
+        idistLine[0]        = (IDistType)0;
+
+        for (int j = 1; j < guide.cols; j++)
+        {
+            curDist += dtf.getTransformedDistance(guideLine[j-1], guideLine[j]);
+            idistLine[j] = curDist;
+        }
+        idistLine[guide.cols] = std::numeric_limits<IDistType>::max();
+    }
+}
+
+template <typename GuideVec>
+DTFilterCPU::ComputeDTandIDTHor_ParBody<GuideVec>::ComputeDTandIDTHor_ParBody(DTFilterCPU& dtf_, Mat& guide_, Mat& dist_, Mat& idist_)
+: dtf(dtf_), guide(guide_), dist(dist_), idist(idist_)
+{
+    dist = getWExtendedMat(guide.rows, guide.cols, IDistVec::type, 1, 1);
+    idist = getWExtendedMat(guide.rows, guide.cols + 1, IDistVec::type);
+    maxRadius = dtf.getIterRadius(1);
+}
+
+template <typename GuideVec>
+void DTFilterCPU::ComputeDTandIDTHor_ParBody<GuideVec>::operator()(const Range& range) const
+{
+    for (int i = range.start; i < range.end; i++)
+    {
+        const GuideVec  *guideLine  = guide.ptr<GuideVec>(i);
+              DistType  *distLine   = dist.ptr<DistType>(i);
+              IDistType *idistLine  = idist.ptr<IDistType>(i);
+
+        DistType  curDist;
+        IDistType curIDist = (IDistType)0;
+        int j;
+                
+        distLine[-1] = maxRadius;
+        //idistLine[-1] = curIDist - maxRadius;
+        idistLine[0] = curIDist;
+        for (j = 0; j < guide.cols-1; j++)
+        {
+            curDist = (DistType) dtf.getTransformedDistance(guideLine[j], guideLine[j + 1]);
+            curIDist += curDist;
+
+            distLine[j] = curDist;
+            idistLine[j + 1] = curIDist;
+        }
+        idistLine[j + 1] = curIDist + maxRadius;
+        distLine[j] = maxRadius;
+    }
+}
+
+template <typename GuideVec>
+DTFilterCPU::ComputeA0DTHor_ParBody<GuideVec>::ComputeA0DTHor_ParBody(DTFilterCPU& dtf_, Mat& guide_)
+: dtf(dtf_), guide(guide_)
+{
+    dtf.a0distHor.create(guide.rows, guide.cols - 1, DistVec::type);
+    lna = std::log(dtf.getIterAlpha(1));
+}
+
+template <typename GuideVec>
+void DTFilterCPU::ComputeA0DTHor_ParBody<GuideVec>::operator()(const Range& range) const
+{
+    for (int i = range.start; i < range.end; i++)
+    {
+        const GuideVec  *guideLine  = guide.ptr<GuideVec>(i);
+              DistType  *dstLine    = dtf.a0distHor.ptr<DistType>(i);
+
+        for (int j = 0; j < guide.cols - 1; j++)
+        {
+            DistType d = (DistType)dtf.getTransformedDistance(guideLine[j], guideLine[j + 1]);
+            dstLine[j] = lna*d;
+        }
+    }
+}
+
+template <typename GuideVec>
+DTFilterCPU::ComputeA0DTHor_ParBody<GuideVec>::~ComputeA0DTHor_ParBody()
+{
+    cv::exp(dtf.a0distHor, dtf.a0distHor);
+}
+
+template <typename GuideVec>
+DTFilterCPU::ComputeA0DTVert_ParBody<GuideVec>::ComputeA0DTVert_ParBody(DTFilterCPU& dtf_, Mat& guide_)
+: dtf(dtf_), guide(guide_)
+{
+    dtf.a0distVert.create(guide.rows - 1, guide.cols, DistVec::type);
+    lna = std::log(dtf.getIterAlpha(1));
+}
+
+template <typename GuideVec>
+void DTFilterCPU::ComputeA0DTVert_ParBody<GuideVec>::operator()(const Range& range) const
+{
+    for (int i = range.start; i < range.end; i++)
+    {
+        DistType *dstLine = dtf.a0distVert.ptr<DistType>(i);
+        GuideVec *guideRow1 = guide.ptr<GuideVec>(i);
+        GuideVec *guideRow2 = guide.ptr<GuideVec>(i+1);
+
+        for (int j = 0; j < guide.cols; j++)
+        {
+            DistType d = (DistType)dtf.getTransformedDistance(guideRow1[j], guideRow2[j]);
+            dstLine[j] = lna*d;
+        }
+    }
+}
+
+template <typename GuideVec>
+DTFilterCPU::ComputeA0DTVert_ParBody<GuideVec>::~ComputeA0DTVert_ParBody()
+{
+    cv::exp(dtf.a0distVert, dtf.a0distVert);
+}
+
+
+template<typename GuideVec, typename SrcVec>
+void domainTransformFilter( const Mat_<GuideVec>& guide,
+                            const Mat_<SrcVec>& source,
+                            Mat& dst,
+                            double sigmaSpatial, double sigmaColor,
+                            int mode, int numPasses
+                            )
+{
+    DTFilterCPU *dtf = DTFilterCPU::create_p_<GuideVec>(guide, sigmaSpatial, sigmaColor, mode, numPasses);
+    dtf->filter_<SrcVec>(source, dst);
+    delete dtf;
+}
+
+template<typename GuideVec, typename SrcVec>
+void domainTransformFilter( const Mat& guide,
+                            const Mat& source,
+                            Mat& dst,
+                            double sigmaSpatial, double sigmaColor,
+                            int mode, int numPasses
+                          )
+{
+    DTFilterCPU *dtf = DTFilterCPU::create_p_<GuideVec>(guide, sigmaSpatial, sigmaColor, mode, numPasses);
+    dtf->filter_<SrcVec>(source, dst);
+    delete dtf;
+}
+
+}
+}
+#endif

modules/ximgproc/src/edgeaware_filters_common.cpp

@@ -0,0 +1,515 @@
+#include "precomp.hpp"
+#include "edgeaware_filters_common.hpp"
+#include "dtfilter_cpu.hpp"
+
+#include <opencv2/core/cvdef.h>
+#include <opencv2/core/utility.hpp>
+#include <cmath>
+using namespace std;
+
+#ifndef SQR
+#define SQR(x) ((x)*(x))
+#endif
+
+#if defined(CV_SSE)
+static volatile bool CPU_SUPPORT_SSE1 = cv::checkHardwareSupport(CV_CPU_SSE);
+#endif
+
+#ifdef CV_SSE2
+static volatile bool CPU_SUPPORT_SSE2 = cv::checkHardwareSupport(CV_CPU_SSE2);
+#endif
+
+namespace cv
+{
+namespace ximgproc
+{
+
+Ptr<DTFilter> createDTFilterRF(InputArray adistHor, InputArray adistVert, double sigmaSpatial, double sigmaColor, int numIters)
+{
+    return Ptr<DTFilter>(DTFilterCPU::createRF(adistHor, adistVert, sigmaSpatial, sigmaColor, numIters));
+}
+
+int getTotalNumberOfChannels(InputArrayOfArrays src)
+{
+    CV_Assert(src.isMat() || src.isUMat() || src.isMatVector() || src.isUMatVector());
+
+    if (src.isMat() || src.isUMat())
+    {
+        return src.channels();
+    }
+    else if (src.isMatVector())
+    {
+        int cnNum = 0;
+        const vector<Mat>& srcv = *static_cast<const vector<Mat>*>(src.getObj());
+        for (unsigned i = 0; i < srcv.size(); i++)
+            cnNum += srcv[i].channels();
+        return cnNum;
+    }
+    else if (src.isUMatVector())
+    {
+        int cnNum = 0;
+        const vector<UMat>& srcv = *static_cast<const vector<UMat>*>(src.getObj());
+        for (unsigned i = 0; i < srcv.size(); i++)
+            cnNum += srcv[i].channels();
+        return cnNum;
+    }
+    else
+    {
+        return 0;
+    }
+}
+
+void checkSameSizeAndDepth(InputArrayOfArrays src, Size &sz, int &depth)
+{
+    CV_Assert(src.isMat() || src.isUMat() || src.isMatVector() || src.isUMatVector());
+
+    if (src.isMat() || src.isUMat())
+    {
+        CV_Assert(!src.empty());
+        sz = src.size();
+        depth = src.depth();
+    }
+    else if (src.isMatVector())
+    {
+        const vector<Mat>& srcv = *static_cast<const vector<Mat>*>(src.getObj());
+        CV_Assert(srcv.size() > 0);
+        for (unsigned i = 0; i < srcv.size(); i++)
+        {
+            CV_Assert(srcv[i].depth() == srcv[0].depth());
+            CV_Assert(srcv[i].size() == srcv[0].size());
+        }
+        sz = srcv[0].size();
+        depth = srcv[0].depth();
+    }
+    else if (src.isUMatVector())
+    {
+        const vector<UMat>& srcv = *static_cast<const vector<UMat>*>(src.getObj());
+        CV_Assert(srcv.size() > 0);
+        for (unsigned i = 0; i < srcv.size(); i++)
+        {
+            CV_Assert(srcv[i].depth() == srcv[0].depth());
+            CV_Assert(srcv[i].size() == srcv[0].size());
+        }
+        sz = srcv[0].size();
+        depth = srcv[0].depth();
+    }
+}
+
+namespace intrinsics
+{
+
+inline float getFloatSignBit()
+{
+    union
+    {
+        int signInt;
+        float signFloat;
+    };
+    signInt = 0x80000000;
+
+    return signFloat;
+}
+
+void add_(register float *dst, register float *src1, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a, b;
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_loadu_ps(src1 + j);
+            b = _mm_loadu_ps(dst + j);
+            b = _mm_add_ps(b, a);
+            _mm_storeu_ps(dst + j, b);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        dst[j] += src1[j];
+}
+
+void mul(register float *dst, register float *src1, register float *src2, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a, b;
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_loadu_ps(src1 + j);
+            b = _mm_loadu_ps(src2 + j);
+            b = _mm_mul_ps(a, b);
+            _mm_storeu_ps(dst + j, b);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        dst[j] = src1[j] * src2[j];
+}
+
+void mul(register float *dst, register float *src1, float src2, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a, b;
+        b = _mm_set_ps1(src2);
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_loadu_ps(src1 + j);
+            a = _mm_mul_ps(a, b);
+            _mm_storeu_ps(dst + j, a);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        dst[j] = src1[j]*src2;
+}
+
+void mad(register float *dst, register float *src1, float alpha, float beta, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a, b, c;
+        a = _mm_set_ps1(alpha);
+        b = _mm_set_ps1(beta);
+        for (; j < w - 3; j += 4)
+        {
+            c = _mm_loadu_ps(src1 + j);
+            c = _mm_mul_ps(c, a);
+            c = _mm_add_ps(c, b);
+            _mm_storeu_ps(dst + j, c);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        dst[j] = alpha*src1[j] + beta;
+}
+
+void sqr_(register float *dst, register float *src1, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a;
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_loadu_ps(src1 + j);
+            a = _mm_mul_ps(a, a);
+            _mm_storeu_ps(dst + j, a);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        dst[j] = src1[j] * src1[j];
+}
+
+void sqr_dif(register float *dst, register float *src1, register float *src2, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 d;
+        for (; j < w - 3; j += 4)
+        {
+            d = _mm_sub_ps(_mm_loadu_ps(src1 + j), _mm_loadu_ps(src2 + j));
+            d = _mm_mul_ps(d, d);
+            _mm_storeu_ps(dst + j, d);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        dst[j] = (src1[j] - src2[j])*(src1[j] - src2[j]);
+}
+
+void add_mul(register float *dst, register float *src1, register float *src2, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a, b, c;
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_loadu_ps(src1 + j);
+            b = _mm_loadu_ps(src2 + j);
+            b = _mm_mul_ps(b, a);
+            c = _mm_loadu_ps(dst + j);
+            c = _mm_add_ps(c, b);
+            _mm_storeu_ps(dst + j, c);
+        }
+    }
+#endif
+    for (; j < w; j++)
+    {
+        dst[j] += src1[j] * src2[j];
+    }
+}
+
+void add_sqr(register float *dst, register float *src1, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a, c;
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_loadu_ps(src1 + j);
+            a = _mm_mul_ps(a, a);
+            c = _mm_loadu_ps(dst + j);
+            c = _mm_add_ps(c, a);
+            _mm_storeu_ps(dst + j, c);
+        }
+    }
+#endif
+    for (; j < w; j++)
+    {
+        dst[j] += src1[j] * src1[j];
+    }
+}
+
+void add_sqr_dif(register float *dst, register float *src1, register float *src2, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a, d;
+        for (; j < w - 3; j += 4)
+        {
+            d = _mm_sub_ps(_mm_loadu_ps(src1 + j), _mm_loadu_ps(src2 + j));
+            d = _mm_mul_ps(d, d);
+            a = _mm_loadu_ps(dst + j);
+            a = _mm_add_ps(a, d);
+            _mm_storeu_ps(dst + j, a);
+        }
+    }
+#endif
+    for (; j < w; j++)
+    {
+        dst[j] += (src1[j] - src2[j])*(src1[j] - src2[j]);
+    }
+}
+
+void sub_mul(register float *dst, register float *src1, register float *src2, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a, b, c;
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_loadu_ps(src1 + j);
+            b = _mm_loadu_ps(src2 + j);
+            b = _mm_mul_ps(b, a);
+            c = _mm_loadu_ps(dst + j);
+            c = _mm_sub_ps(c, b);
+            _mm_storeu_ps(dst + j, c);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        dst[j] -= src1[j] * src2[j];
+}
+
+void sub_mad(register float *dst, register float *src1, register float *src2, float c0, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a, b, c;
+        __m128 cnst = _mm_set_ps1(c0);
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_loadu_ps(src1 + j);
+            b = _mm_loadu_ps(src2 + j);
+            b = _mm_mul_ps(b, a);
+            c = _mm_loadu_ps(dst + j);
+            c = _mm_sub_ps(c, cnst);
+            c = _mm_sub_ps(c, b);
+            _mm_storeu_ps(dst + j, c);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        dst[j] -= src1[j] * src2[j] + c0;
+}
+
+void det_2x2(register float *dst, register float *a00, register float *a01, register float *a10, register float *a11, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a, b;
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_mul_ps(_mm_loadu_ps(a00 + j), _mm_loadu_ps(a11 + j));
+            b = _mm_mul_ps(_mm_loadu_ps(a01 + j), _mm_loadu_ps(a10 + j));
+            a = _mm_sub_ps(a, b);
+            _mm_storeu_ps(dst + j, a);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        dst[j] = a00[j]*a11[j] - a01[j]*a10[j];
+}
+
+void div_det_2x2(register float *a00, register float *a01, register float *a11, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        const __m128 SIGN_MASK = _mm_set_ps1(getFloatSignBit());
+
+        __m128 a, b, _a00, _a01, _a11;
+        for (; j < w - 3; j += 4)
+        {
+            _a00 = _mm_loadu_ps(a00 + j);
+            _a11 = _mm_loadu_ps(a11 + j);
+            a = _mm_mul_ps(_a00, _a11);
+
+            _a01 = _mm_loadu_ps(a01 + j);
+            _a01 = _mm_xor_ps(_a01, SIGN_MASK);
+            b = _mm_mul_ps(_a01, _a01);
+            
+            a = _mm_sub_ps(a, b);
+            
+            _a01 = _mm_div_ps(_a01, a);
+            _a00 = _mm_div_ps(_a00, a);
+            _a11 = _mm_div_ps(_a11, a);
+            
+            _mm_storeu_ps(a01 + j, _a01);
+            _mm_storeu_ps(a00 + j, _a00);
+            _mm_storeu_ps(a11 + j, _a11);
+        }
+    }
+#endif
+    for (; j < w; j++)
+    {
+        float det = a00[j] * a11[j] - a01[j] * a01[j];
+        a00[j] /= det;
+        a11[j] /= det;
+        a01[j] /= -det;
+    }
+}
+
+void div_1x(register float *a1, register float *b1, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 _a1, _b1;
+        for (; j < w - 3; j += 4)
+        {
+            _b1 = _mm_loadu_ps(b1 + j);
+            _a1 = _mm_loadu_ps(a1 + j);
+            _mm_storeu_ps(a1 + j, _mm_div_ps(_a1, _b1));
+        }
+    }
+#endif
+    for (; j < w; j++)
+    {
+        a1[j] /= b1[j];
+    }
+}
+
+void inv_self(register float *src, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a;
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_rcp_ps(_mm_loadu_ps(src + j));
+            _mm_storeu_ps(src + j, a);
+        }
+    }
+#endif
+    for (; j < w; j++)
+    {
+        src[j] = 1.0f / src[j];
+    }
+}
+
+void sqrt_(register float *dst, register float *src, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a;
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_sqrt_ps(_mm_loadu_ps(src + j));
+            _mm_storeu_ps(dst + j, a);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        dst[j] = sqrt(src[j]);
+}
+
+void min_(register float *dst, register float *src1, register float *src2, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 a, b;
+        for (; j < w - 3; j += 4)
+        {
+            a = _mm_loadu_ps(src1 + j);
+            b = _mm_loadu_ps(src2 + j);
+            b = _mm_min_ps(b, a);
+
+            _mm_storeu_ps(dst + j, b);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        dst[j] = std::min(src1[j], src2[j]);
+}
+
+void rf_vert_row_pass(register float *curRow, register float *prevRow, float alphaVal, int w)
+{
+    register int j = 0;
+#ifdef CV_SSE
+    if (CPU_SUPPORT_SSE1)
+    {
+        __m128 cur, prev, res;
+        __m128 alpha = _mm_set_ps1(alphaVal);
+        for (; j < w - 3; j += 4)
+        {
+            cur = _mm_loadu_ps(curRow + j);
+            prev = _mm_loadu_ps(prevRow + j);
+
+            res = _mm_mul_ps(alpha, _mm_sub_ps(prev, cur));
+            res = _mm_add_ps(res, cur);
+            _mm_storeu_ps(curRow + j, res);
+        }
+    }
+#endif
+    for (; j < w; j++)
+        curRow[j] += alphaVal*(prevRow[j] - curRow[j]);
+}
+
+} //end of cv::ximgproc::intrinsics
+
+} //end of cv::ximgproc 
+} //end of cv 

modules/ximgproc/src/edgeaware_filters_common.hpp

@@ -0,0 +1,61 @@
+#ifndef __EDGEAWAREFILTERS_COMMON_HPP__
+#define __EDGEAWAREFILTERS_COMMON_HPP__
+#ifdef __cplusplus
+
+namespace cv
+{
+namespace ximgproc
+{
+
+Ptr<DTFilter> createDTFilterRF(InputArray adistHor, InputArray adistVert, double sigmaSpatial, double sigmaColor, int numIters);
+
+int getTotalNumberOfChannels(InputArrayOfArrays src);
+
+void checkSameSizeAndDepth(InputArrayOfArrays src, Size &sz, int &depth);
+
+namespace intrinsics
+{  
+    void add_(register float *dst, register float *src1, int w);
+
+    void mul(register float *dst, register float *src1, register float *src2, int w);
+
+    void mul(register float *dst, register float *src1, float src2, int w);
+
+    //dst = alpha*src + beta
+    void mad(register float *dst, register float *src1, float alpha, float beta, int w);
+
+    void add_mul(register float *dst, register float *src1, register float *src2, int w);
+
+    void sub_mul(register float *dst, register float *src1, register float *src2, int w);
+
+    void sub_mad(register float *dst, register float *src1, register float *src2, float c0, int w);
+
+    void det_2x2(register float *dst, register float *a00, register float *a01, register float *a10, register float *a11, int w);
+
+    void div_det_2x2(register float *a00, register float *a01, register float *a11, int w);
+
+    void div_1x(register float *a1, register float *b1, int w);
+
+    void inv_self(register float *src, int w);
+
+    
+    void sqr_(register float *dst, register float *src1, int w);
+
+    void sqrt_(register float *dst, register float *src, int w);
+
+    void sqr_dif(register float *dst, register float *src1, register float *src2, int w);
+
+    void add_sqr_dif(register float *dst, register float *src1, register float *src2, int w);
+
+    void add_sqr(register float *dst, register float *src1, int w);
+
+    void min_(register float *dst, register float *src1, register float *src2, int w);
+
+    void rf_vert_row_pass(register float *curRow, register float *prevRow, float alphaVal, int w);
+}
+
+}
+}
+
+#endif
+#endif

modules/ximgproc/src/guided_filter.cpp

@@ -0,0 +1,755 @@
+#include "precomp.hpp"
+#include "edgeaware_filters_common.hpp"
+#include <vector>
+
+#ifdef _MSC_VER
+#   pragma warning(disable: 4512)
+#endif
+
+namespace cv
+{
+namespace ximgproc
+{
+
+using std::vector;
+using namespace cv::ximgproc::intrinsics;
+
+template <typename T>
+struct SymArray2D
+{
+    vector<T> vec;
+    int sz;
+
+    SymArray2D()
+    {
+        sz = 0;
+    }
+
+    void create(int sz_)
+    {
+        CV_DbgAssert(sz_ > 0);
+        sz = sz_;
+        vec.resize(total());
+    }
+
+    inline T& operator()(int i, int j)
+    {
+        CV_DbgAssert(i >= 0 && i < sz && j >= 0 && j < sz);
+        if (i < j) std::swap(i, j);
+        return vec[i*(i+1)/2 + j];
+    }
+
+    inline T& operator()(int i)
+    {
+        return vec[i];
+    }
+
+    int total() const
+    {
+        return sz*(sz + 1)/2;
+    }
+
+    void release()
+    {
+        vec.clear();
+        sz = 0;
+    }
+};
+
+
+template <typename XMat>
+static void splitFirstNChannels(InputArrayOfArrays src, vector<XMat>& dst, int maxDstCn)
+{
+    CV_Assert(src.isMat() || src.isUMat() || src.isMatVector() || src.isUMatVector());
+
+    if ( (src.isMat() || src.isUMat()) && src.channels() == maxDstCn )
+    {
+        split(src, dst);
+    }
+    else
+    {
+        Size sz;
+        int depth, totalCnNum;
+        
+        checkSameSizeAndDepth(src, sz, depth);
+        totalCnNum = std::min(maxDstCn, getTotalNumberOfChannels(src));
+        
+        dst.resize(totalCnNum);
+        vector<int> fromTo(2*totalCnNum);
+        for (int i = 0; i < totalCnNum; i++)
+        {
+            fromTo[i*2 + 0] = i;
+            fromTo[i*2 + 1] = i;
+
+            dst[i].create(sz, CV_MAKE_TYPE(depth, 1));
+        }
+
+        mixChannels(src, dst, fromTo);
+    }
+}
+
+class GuidedFilterImpl : public GuidedFilter
+{
+public:
+    
+    static Ptr<GuidedFilterImpl> create(InputArray guide, int radius, double eps);
+
+    void filter(InputArray src, OutputArray dst, int dDepth = -1);
+
+protected:
+
+    int radius;
+    double eps;
+    int h, w;
+
+    vector<Mat> guideCn;
+    vector<Mat> guideCnMean;
+
+    SymArray2D<Mat> covarsInv;
+
+    int gCnNum;
+
+protected:
+
+    GuidedFilterImpl() {}
+    
+    void init(InputArray guide, int radius, double eps);
+
+    void computeCovGuide(SymArray2D<Mat>& covars);
+
+    void computeCovGuideAndSrc(vector<Mat>& srcCn, vector<Mat>& srcCnMean, vector<vector<Mat> >& cov);
+
+    void getWalkPattern(int eid, int &cn1, int &cn2);
+
+    inline void meanFilter(Mat& src, Mat& dst)
+    {
+        boxFilter(src, dst, CV_32F, Size(2 * radius + 1, 2 * radius + 1), cv::Point(-1, -1), true, BORDER_REFLECT);
+    }
+
+    inline void convertToWorkType(Mat& src, Mat& dst)
+    {
+        src.convertTo(dst, CV_32F);
+    }
+
+private: /*Routines to parallelize boxFilter and convertTo*/
+    
+    typedef void (GuidedFilterImpl::*TransformFunc)(Mat& src, Mat& dst);
+
+    struct GFTransform_ParBody : public ParallelLoopBody
+    {
+        GuidedFilterImpl &gf;
+        mutable vector<Mat*> src;
+        mutable vector<Mat*> dst;
+        TransformFunc func;
+
+        GFTransform_ParBody(GuidedFilterImpl& gf_, vector<Mat>& srcv, vector<Mat>& dstv, TransformFunc func_);
+        GFTransform_ParBody(GuidedFilterImpl& gf_, vector<vector<Mat> >& srcvv, vector<vector<Mat> >& dstvv, TransformFunc func_);
+
+        void operator () (const Range& range) const;
+
+        Range getRange() const
+        {
+            return Range(0, (int)src.size());
+        }
+    };
+
+    template<typename V>
+    void parConvertToWorkType(V &src, V &dst)
+    {
+        GFTransform_ParBody pb(*this, src, dst, &GuidedFilterImpl::convertToWorkType);
+        parallel_for_(pb.getRange(), pb);
+    }
+
+    template<typename V>
+    void parMeanFilter(V &src, V &dst)
+    {
+        GFTransform_ParBody pb(*this, src, dst, &GuidedFilterImpl::meanFilter);
+        parallel_for_(pb.getRange(), pb);
+    }
+
+private: /*Parallel body classes*/
+
+    inline void runParBody(const ParallelLoopBody& pb)
+    {
+        parallel_for_(Range(0, h), pb);
+    }
+
+    struct MulChannelsGuide_ParBody : public ParallelLoopBody
+    {
+        GuidedFilterImpl &gf;
+        SymArray2D<Mat> &covars;
+
+        MulChannelsGuide_ParBody(GuidedFilterImpl& gf_, SymArray2D<Mat>& covars_)
+            : gf(gf_), covars(covars_) {}
+
+        void operator () (const Range& range) const;
+    };
+
+    struct ComputeCovGuideFromChannelsMul_ParBody : public ParallelLoopBody
+    {
+        GuidedFilterImpl &gf;
+        SymArray2D<Mat> &covars;
+
+        ComputeCovGuideFromChannelsMul_ParBody(GuidedFilterImpl& gf_, SymArray2D<Mat>& covars_)
+            : gf(gf_), covars(covars_) {}
+
+        void operator () (const Range& range) const;
+    };
+
+    struct ComputeCovGuideInv_ParBody : public ParallelLoopBody
+    {
+        GuidedFilterImpl &gf;
+        SymArray2D<Mat> &covars;
+
+        ComputeCovGuideInv_ParBody(GuidedFilterImpl& gf_, SymArray2D<Mat>& covars_);
+
+        void operator () (const Range& range) const;
+    };
+
+
+    struct MulChannelsGuideAndSrc_ParBody : public ParallelLoopBody
+    {
+        GuidedFilterImpl &gf;
+        vector<vector<Mat> > &cov;
+        vector<Mat> &srcCn;
+
+        MulChannelsGuideAndSrc_ParBody(GuidedFilterImpl& gf_, vector<Mat>& srcCn_, vector<vector<Mat> >& cov_)
+            : gf(gf_), cov(cov_), srcCn(srcCn_) {}
+
+        void operator () (const Range& range) const;
+    };
+
+    struct ComputeCovFromSrcChannelsMul_ParBody : public ParallelLoopBody
+    {
+        GuidedFilterImpl &gf;
+        vector<vector<Mat> > &cov;
+        vector<Mat> &srcCnMean;
+
+        ComputeCovFromSrcChannelsMul_ParBody(GuidedFilterImpl& gf_, vector<Mat>& srcCnMean_, vector<vector<Mat> >& cov_)
+            : gf(gf_), cov(cov_), srcCnMean(srcCnMean_) {}
+
+        void operator () (const Range& range) const;
+    };
+
+    struct ComputeAlpha_ParBody : public ParallelLoopBody
+    {
+        GuidedFilterImpl &gf;
+        vector<vector<Mat> > &alpha;
+        vector<vector<Mat> > &covSrc;
+
+        ComputeAlpha_ParBody(GuidedFilterImpl& gf_, vector<vector<Mat> >& alpha_, vector<vector<Mat> >& covSrc_)
+            : gf(gf_), alpha(alpha_), covSrc(covSrc_) {}
+
+        void operator () (const Range& range) const;
+    };
+
+    struct ComputeBeta_ParBody : public ParallelLoopBody
+    {
+        GuidedFilterImpl &gf;
+        vector<vector<Mat> > &alpha;
+        vector<Mat> &srcCnMean;
+        vector<Mat> &beta;
+
+        ComputeBeta_ParBody(GuidedFilterImpl& gf_, vector<vector<Mat> >& alpha_, vector<Mat>& srcCnMean_, vector<Mat>& beta_)
+            : gf(gf_), alpha(alpha_), srcCnMean(srcCnMean_), beta(beta_) {}
+
+        void operator () (const Range& range) const;
+    };
+
+    struct ApplyTransform_ParBody : public ParallelLoopBody
+    {
+        GuidedFilterImpl &gf;
+        vector<vector<Mat> > &alpha;
+        vector<Mat> &beta;
+
+        ApplyTransform_ParBody(GuidedFilterImpl& gf_, vector<vector<Mat> >& alpha_, vector<Mat>& beta_)
+            : gf(gf_),  alpha(alpha_), beta(beta_) {}
+
+        void operator () (const Range& range) const;
+    };
+};
+
+void GuidedFilterImpl::MulChannelsGuide_ParBody::operator()(const Range& range) const
+{
+    int total = covars.total();
+
+    for (int i = range.start; i < range.end; i++)
+    {
+        int c1, c2;
+        float *cov, *guide1, *guide2;
+
+        for (int k = 0; k < total; k++)
+        {
+            gf.getWalkPattern(k, c1, c2);
+
+            guide1 = gf.guideCn[c1].ptr<float>(i);
+            guide2 = gf.guideCn[c2].ptr<float>(i);
+            cov = covars(c1, c2).ptr<float>(i);
+
+            mul(cov, guide1, guide2, gf.w);
+        }
+    }
+}
+
+void GuidedFilterImpl::ComputeCovGuideFromChannelsMul_ParBody::operator()(const Range& range) const
+{
+    int total = covars.total();
+    float diagSummand = (float)(gf.eps);
+
+    for (int i = range.start; i < range.end; i++)
+    {
+        int c1, c2;
+        float *cov, *guide1, *guide2;
+
+        for (int k = 0; k < total; k++)
+        {
+            gf.getWalkPattern(k, c1, c2);
+
+            guide1 = gf.guideCnMean[c1].ptr<float>(i);
+            guide2 = gf.guideCnMean[c2].ptr<float>(i);
+            cov = covars(c1, c2).ptr<float>(i);
+
+            if (c1 != c2)
+            {
+                sub_mul(cov, guide1, guide2, gf.w);
+            }
+            else
+            {
+                sub_mad(cov, guide1, guide2, -diagSummand, gf.w);
+            }
+        }
+    }
+}
+
+GuidedFilterImpl::ComputeCovGuideInv_ParBody::ComputeCovGuideInv_ParBody(GuidedFilterImpl& gf_, SymArray2D<Mat>& covars_) 
+    : gf(gf_), covars(covars_)
+{
+    gf.covarsInv.create(gf.gCnNum);
+
+    if (gf.gCnNum == 3)
+    {
+        for (int k = 0; k < 2; k++)
+            for (int l = 0; l < 3; l++)
+                gf.covarsInv(k, l).create(gf.h, gf.w, CV_32FC1);
+
+        ////trick to avoid memory allocation
+        gf.covarsInv(2, 0).create(gf.h, gf.w, CV_32FC1);
+        gf.covarsInv(2, 1) = covars(2, 1);
+        gf.covarsInv(2, 2) = covars(2, 2);
+
+        return;
+    }
+
+    if (gf.gCnNum == 2)
+    {
+        gf.covarsInv(0, 0) = covars(1, 1);
+        gf.covarsInv(0, 1) = covars(0, 1);
+        gf.covarsInv(1, 1) = covars(0, 0);
+        return;
+    }
+
+    if (gf.gCnNum == 1)
+    {
+        gf.covarsInv(0, 0) = covars(0, 0);
+        return;
+    }
+}
+
+void GuidedFilterImpl::ComputeCovGuideInv_ParBody::operator()(const Range& range) const
+{
+    if (gf.gCnNum == 3)
+    {
+        vector<float> covarsDet(gf.w);
+        float *det = &covarsDet[0];
+
+        for (int i = range.start; i < range.end; i++)
+        {
+            for (int k = 0; k < 3; k++)
+                for (int l = 0; l <= k; l++)
+                {
+                    float *dst = gf.covarsInv(k, l).ptr<float>(i);
+
+                    float *a00 = covars((k + 1) % 3, (l + 1) % 3).ptr<float>(i);
+                    float *a01 = covars((k + 1) % 3, (l + 2) % 3).ptr<float>(i);
+                    float *a10 = covars((k + 2) % 3, (l + 1) % 3).ptr<float>(i);
+                    float *a11 = covars((k + 2) % 3, (l + 2) % 3).ptr<float>(i);
+
+                    det_2x2(dst, a00, a01, a10, a11, gf.w);
+                }
+
+            for (int k = 0; k < 3; k++)
+            {
+                register float *a = covars(k, 0).ptr<float>(i);
+                register float *ac = gf.covarsInv(k, 0).ptr<float>(i);
+
+                if (k == 0)
+                    mul(det, a, ac, gf.w);
+                else
+                    add_mul(det, a, ac, gf.w);
+            }
+
+            for (int k = 0; k < gf.covarsInv.total(); k += 1)
+            {
+                div_1x(gf.covarsInv(k).ptr<float>(i), det, gf.w);
+            }
+        }
+        return;
+    }
+
+    if (gf.gCnNum == 2)
+    {
+        for (int i = range.start; i < range.end; i++)
+        {
+            float *a00 = gf.covarsInv(0, 0).ptr<float>(i);
+            float *a10 = gf.covarsInv(1, 0).ptr<float>(i);
+            float *a11 = gf.covarsInv(1, 1).ptr<float>(i);
+
+            div_det_2x2(a00, a10, a11, gf.w);
+        }
+        return;
+    }
+
+    if (gf.gCnNum == 1)
+    {
+        //divide(1.0, covars(0, 0)(range, Range::all()), gf.covarsInv(0, 0)(range, Range::all()));
+        //return;
+
+        for (int i = range.start; i < range.end; i++)
+        {
+            float *res = covars(0, 0).ptr<float>(i);
+            inv_self(res, gf.w);
+        }
+        return;
+    }
+}
+
+void GuidedFilterImpl::MulChannelsGuideAndSrc_ParBody::operator()(const Range& range) const
+{
+    int srcCnNum = (int)srcCn.size();
+
+    for (int i = range.start; i < range.end; i++)
+    {
+        for (int si = 0; si < srcCnNum; si++)
+        {
+            int step    = (si % 2) * 2 - 1;
+            int start   = (si % 2) ? 0 : gf.gCnNum - 1;
+            int end     = (si % 2) ? gf.gCnNum : -1;
+
+            float *srcLine = srcCn[si].ptr<float>(i);
+
+            for (int gi = start; gi != end; gi += step)
+            {
+                float *guideLine = gf.guideCn[gi].ptr<float>(i);
+                float *dstLine = cov[si][gi].ptr<float>(i);
+
+                mul(dstLine, srcLine, guideLine, gf.w);
+            }
+        }
+    }
+}
+
+void GuidedFilterImpl::ComputeCovFromSrcChannelsMul_ParBody::operator()(const Range& range) const
+{
+    int srcCnNum = (int)srcCnMean.size();
+
+    for (int i = range.start; i < range.end; i++)
+    {
+        for (int si = 0; si < srcCnNum; si++)
+        {
+            int step    = (si % 2) * 2 - 1;
+            int start   = (si % 2) ? 0 : gf.gCnNum - 1;
+            int end     = (si % 2) ? gf.gCnNum : -1;
+
+            register float *srcMeanLine = srcCnMean[si].ptr<float>(i);
+
+            for (int gi = start; gi != end; gi += step)
+            {
+                float *guideMeanLine = gf.guideCnMean[gi].ptr<float>(i);
+                float *covLine = cov[si][gi].ptr<float>(i);
+
+                sub_mul(covLine, srcMeanLine, guideMeanLine, gf.w);
+            }
+        }
+    }
+}
+
+void GuidedFilterImpl::ComputeAlpha_ParBody::operator()(const Range& range) const
+{
+    int srcCnNum = (int)covSrc.size();
+
+    for (int i = range.start; i < range.end; i++)
+    {
+        for (int si = 0; si < srcCnNum; si++)
+        {
+            for (int gi = 0; gi < gf.gCnNum; gi++)
+            {
+                float *y, *A, *dstAlpha;
+
+                dstAlpha = alpha[si][gi].ptr<float>(i);
+                for (int k = 0; k < gf.gCnNum; k++)
+                {
+                    y = covSrc[si][k].ptr<float>(i);
+                    A = gf.covarsInv(gi, k).ptr<float>(i);
+
+                    if (k == 0)
+                    {
+                        mul(dstAlpha, A, y, gf.w);
+                    }
+                    else
+                    {
+                        add_mul(dstAlpha, A, y, gf.w);
+                    }
+                }
+            }
+        }
+    }
+}
+
+void GuidedFilterImpl::ComputeBeta_ParBody::operator()(const Range& range) const
+{
+    int srcCnNum = (int)srcCnMean.size();
+    CV_DbgAssert(&srcCnMean == &beta);
+
+    for (int i = range.start; i < range.end; i++)
+    {
+        float *_g[4];
+        for (int gi = 0; gi < gf.gCnNum; gi++)
+            _g[gi] = gf.guideCnMean[gi].ptr<float>(i);
+
+        float *betaDst, *g, *a;
+        for (int si = 0; si < srcCnNum; si++)
+        {
+            betaDst = beta[si].ptr<float>(i);
+            for (int gi = 0; gi < gf.gCnNum; gi++)
+            {
+                a = alpha[si][gi].ptr<float>(i);
+                g = _g[gi];
+
+                sub_mul(betaDst, a, g, gf.w);
+            }
+        }
+    }
+}
+
+void GuidedFilterImpl::ApplyTransform_ParBody::operator()(const Range& range) const
+{
+    int srcCnNum = (int)alpha.size();
+
+    for (int i = range.start; i < range.end; i++)
+    {
+        float *_g[4];
+        for (int gi = 0; gi < gf.gCnNum; gi++)
+            _g[gi] = gf.guideCn[gi].ptr<float>(i);
+
+        float *betaDst, *g, *a;
+        for (int si = 0; si < srcCnNum; si++)
+        {
+            betaDst = beta[si].ptr<float>(i);
+            for (int gi = 0; gi < gf.gCnNum; gi++)
+            {
+                a = alpha[si][gi].ptr<float>(i);
+                g = _g[gi];
+
+                add_mul(betaDst, a, g, gf.w);
+            }
+        }
+    }
+}
+
+GuidedFilterImpl::GFTransform_ParBody::GFTransform_ParBody(GuidedFilterImpl& gf_, vector<Mat>& srcv, vector<Mat>& dstv, TransformFunc func_)
+    : gf(gf_), func(func_)
+{
+    CV_DbgAssert(srcv.size() == dstv.size());
+    src.resize(srcv.size());
+    dst.resize(srcv.size());
+
+    for (int i = 0; i < (int)srcv.size(); i++)
+    {
+        src[i] = &srcv[i];
+        dst[i] = &dstv[i];
+    }
+}
+
+GuidedFilterImpl::GFTransform_ParBody::GFTransform_ParBody(GuidedFilterImpl& gf_, vector<vector<Mat> >& srcvv, vector<vector<Mat> >& dstvv, TransformFunc func_)
+    : gf(gf_), func(func_)
+{
+    CV_DbgAssert(srcvv.size() == dstvv.size());
+    int n = (int)srcvv.size();
+    int total = 0;
+
+    for (int i = 0; i < n; i++)
+    {
+        CV_DbgAssert(srcvv[i].size() == dstvv[i].size());
+        total += (int)srcvv[i].size();
+    }
+
+    src.resize(total);
+    dst.resize(total);
+
+    int k = 0;
+    for (int i = 0; i < n; i++)
+    {
+        for (int j = 0; j < (int)srcvv[i].size(); j++)
+        {
+            src[k] = &srcvv[i][j];
+            dst[k] = &dstvv[i][j];
+            k++;
+        }
+    }
+}
+
+void GuidedFilterImpl::GFTransform_ParBody::operator()(const Range& range) const
+{
+    for (int i = range.start; i < range.end; i++)
+    {
+        (gf.*func)(*src[i], *dst[i]);
+    }
+}
+
+void GuidedFilterImpl::getWalkPattern(int eid, int &cn1, int &cn2)
+{
+    static int wdata[] = {
+        0, -1, -1, -1, -1, -1,
+        0, -1, -1, -1, -1, -1,
+
+        0,  0,  1, -1, -1, -1,
+        0,  1,  1, -1, -1, -1,
+
+        0,  0,  0,  2,  1,  1, 
+        0,  1,  2,  2,  2,  1,
+    };
+
+    cn1 = wdata[6 * 2 * (gCnNum-1) + eid];
+    cn2 = wdata[6 * 2 * (gCnNum-1) + 6 + eid];
+}
+
+Ptr<GuidedFilterImpl> GuidedFilterImpl::create(InputArray guide, int radius, double eps)
+{
+    GuidedFilterImpl *gf = new GuidedFilterImpl();
+    gf->init(guide, radius, eps);
+    return Ptr<GuidedFilterImpl>(gf);
+}
+
+void GuidedFilterImpl::init(InputArray guide, int radius_, double eps_)
+{
+    CV_Assert( !guide.empty() && radius_ >= 0 && eps_ >= 0 );
+    CV_Assert( (guide.depth() == CV_32F || guide.depth() == CV_8U || guide.depth() == CV_16U) && (guide.channels() <= 3) );
+
+    radius = radius_;
+    eps = eps_;
+
+    splitFirstNChannels(guide, guideCn, 3);
+    gCnNum = (int)guideCn.size();
+    h = guideCn[0].rows;
+    w = guideCn[0].cols;
+
+    guideCnMean.resize(gCnNum);
+    parConvertToWorkType(guideCn, guideCn);
+    parMeanFilter(guideCn, guideCnMean);
+    
+    SymArray2D<Mat> covars;
+    computeCovGuide(covars);
+    runParBody(ComputeCovGuideInv_ParBody(*this, covars));
+    covars.release();
+}
+
+void GuidedFilterImpl::computeCovGuide(SymArray2D<Mat>& covars)
+{
+    covars.create(gCnNum);
+    for (int i = 0; i < covars.total(); i++)
+        covars(i).create(h, w, CV_32FC1);
+
+    runParBody(MulChannelsGuide_ParBody(*this, covars));
+
+    parMeanFilter(covars.vec, covars.vec);
+
+    runParBody(ComputeCovGuideFromChannelsMul_ParBody(*this, covars));
+}
+
+void GuidedFilterImpl::filter(InputArray src, OutputArray dst, int dDepth /*= -1*/)
+{
+    CV_Assert( !src.empty() && (src.depth() == CV_32F || src.depth() == CV_8U) );
+    if (src.rows() != h || src.cols() != w)
+    {
+        CV_Error(Error::StsBadSize, "Size of filtering image must be equal to size of guide image");
+        return;
+    }
+
+    if (dDepth == -1) dDepth = src.depth();
+    int srcCnNum = src.channels();
+
+    vector<Mat> srcCn(srcCnNum);
+    vector<Mat>& srcCnMean = srcCn;
+    split(src, srcCn);
+
+    if (src.depth() != CV_32F)
+    {
+        parConvertToWorkType(srcCn, srcCn);
+    }
+
+    vector<vector<Mat> > covSrcGuide(srcCnNum);
+    computeCovGuideAndSrc(srcCn, srcCnMean, covSrcGuide);
+
+    vector<vector<Mat> > alpha(srcCnNum);
+    for (int si = 0; si < srcCnNum; si++)
+    {
+        alpha[si].resize(gCnNum);
+        for (int gi = 0; gi < gCnNum; gi++)
+            alpha[si][gi].create(h, w, CV_32FC1);
+    }
+    runParBody(ComputeAlpha_ParBody(*this, alpha, covSrcGuide));
+    covSrcGuide.clear();
+
+    vector<Mat>& beta = srcCnMean;
+    runParBody(ComputeBeta_ParBody(*this, alpha, srcCnMean, beta));
+
+    parMeanFilter(beta, beta);
+    parMeanFilter(alpha, alpha);
+
+    runParBody(ApplyTransform_ParBody(*this, alpha, beta));
+    if (dDepth != CV_32F)
+    {
+        for (int i = 0; i < srcCnNum; i++)
+            beta[i].convertTo(beta[i], dDepth);
+    }
+    merge(beta, dst);
+}
+
+void GuidedFilterImpl::computeCovGuideAndSrc(vector<Mat>& srcCn, vector<Mat>& srcCnMean, vector<vector<Mat> >& cov)
+{
+    int srcCnNum = (int)srcCn.size();
+
+    cov.resize(srcCnNum);
+    for (int i = 0; i < srcCnNum; i++)
+    {
+        cov[i].resize(gCnNum);
+        for (int j = 0; j < gCnNum; j++)
+            cov[i][j].create(h, w, CV_32FC1);
+    }
+
+    runParBody(MulChannelsGuideAndSrc_ParBody(*this, srcCn, cov));
+
+    parMeanFilter(srcCn, srcCnMean);
+    parMeanFilter(cov, cov);
+
+    runParBody(ComputeCovFromSrcChannelsMul_ParBody(*this, srcCnMean, cov));
+}
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+CV_EXPORTS_W
+Ptr<GuidedFilter> createGuidedFilter(InputArray guide, int radius, double eps)
+{
+    return Ptr<GuidedFilter>(GuidedFilterImpl::create(guide, radius, eps));
+}
+
+CV_EXPORTS_W
+void guidedFilter(InputArray guide, InputArray src, OutputArray dst, int radius, double eps, int dDepth)
+{
+    Ptr<GuidedFilter> gf = createGuidedFilter(guide, radius, eps);
+    gf->filter(src, dst, dDepth);
+}
+
+}
+}

modules/ximgproc/src/joint_bilateral_filter.cpp

@@ -0,0 +1,366 @@
+#include "precomp.hpp"
+#include <climits>
+#include <iostream>
+using namespace std;
+
+#ifdef _MSC_VER
+#   pragma warning(disable: 4512)
+#endif
+
+namespace cv
+{
+namespace ximgproc
+{
+
+typedef Vec<float, 1> Vec1f;
+typedef Vec<uchar, 1> Vec1b;
+
+#ifndef SQR
+#define SQR(a) ((a)*(a))
+#endif
+
+void jointBilateralFilter_32f(Mat& joint, Mat& src, Mat& dst, int radius, double sigmaColor, double sigmaSpace, int borderType);
+
+void jointBilateralFilter_8u(Mat& joint, Mat& src, Mat& dst, int radius, double sigmaColor, double sigmaSpace, int borderType);
+
+template<typename JointVec, typename SrcVec>
+class JointBilateralFilter_32f : public ParallelLoopBody
+{
+    Mat &joint, &src;
+    Mat &dst;
+    int radius, maxk;
+    float scaleIndex;
+    int *spaceOfs;
+    float *spaceWeights, *expLUT;
+
+public:
+
+    JointBilateralFilter_32f(Mat& joint_, Mat& src_, Mat& dst_, int radius_,
+        int maxk_, float scaleIndex_, int *spaceOfs_, float *spaceWeights_, float *expLUT_)
+        :
+        joint(joint_), src(src_), dst(dst_), radius(radius_), maxk(maxk_), 
+        scaleIndex(scaleIndex_), spaceOfs(spaceOfs_), spaceWeights(spaceWeights_), expLUT(expLUT_)
+    {
+        CV_DbgAssert(joint.type() == JointVec::type && src.type() == dst.type() && src.type() == SrcVec::type);
+        CV_DbgAssert(joint.rows == src.rows && src.rows == dst.rows + 2*radius);
+        CV_DbgAssert(joint.cols == src.cols && src.cols == dst.cols + 2*radius);
+    }
+
+    void operator () (const Range& range) const
+    {
+        for (int i = radius + range.start; i < radius + range.end; i++)
+        {
+            for (int j = radius; j < src.cols - radius; j++)
+            {
+                JointVec *jointCenterPixPtr = joint.ptr<JointVec>(i) + j;
+                SrcVec *srcCenterPixPtr = src.ptr<SrcVec>(i) + j;
+
+                JointVec jointPix0 = *jointCenterPixPtr;
+                SrcVec srcSum = SrcVec::all(0.0f);
+                float wSum = 0.0f;
+
+                for (int k = 0; k < maxk; k++)
+                {
+                    float *jointPix = reinterpret_cast<float*>(jointCenterPixPtr + spaceOfs[k]);
+                    float alpha = 0.0f;
+
+                    for (int cn = 0; cn < JointVec::channels; cn++)
+                        alpha += std::abs(jointPix0[cn] - jointPix[cn]);
+                    alpha *= scaleIndex;
+                    int idx = (int)(alpha);
+                    alpha -= idx;
+                    float weight = spaceWeights[k] * (expLUT[idx] + alpha*(expLUT[idx + 1] - expLUT[idx]));
+
+                    float *srcPix = reinterpret_cast<float*>(srcCenterPixPtr + spaceOfs[k]);
+                    for (int cn = 0; cn < SrcVec::channels; cn++)
+                        srcSum[cn] += weight*srcPix[cn];
+                    wSum += weight;
+                }
+
+                dst.at<SrcVec>(i - radius, j - radius) = srcSum / wSum;
+            }
+        }
+    }
+};
+
+void jointBilateralFilter_32f(Mat& joint, Mat& src, Mat& dst, int radius, double sigmaColor, double sigmaSpace, int borderType)
+{
+    CV_DbgAssert(joint.depth() == CV_32F && src.depth() == CV_32F);
+    
+    int d = 2*radius + 1;
+    int jCn = joint.channels();
+    const int kExpNumBinsPerChannel = 1 << 12;
+    double minValJoint, maxValJoint;
+
+    minMaxLoc(joint, &minValJoint, &maxValJoint);
+    if (abs(maxValJoint - minValJoint) < FLT_EPSILON)
+    {
+        //TODO: make circle pattern instead of square
+        GaussianBlur(src, dst, Size(d, d), sigmaSpace, 0, borderType);
+        return;
+    }
+    float colorRange = (float)(maxValJoint - minValJoint) * jCn;
+    colorRange = std::max(0.01f, colorRange);
+
+    int kExpNumBins = kExpNumBinsPerChannel * jCn;
+    vector<float> expLUTv(kExpNumBins + 2);
+    float *expLUT = &expLUTv[0];
+    float scaleIndex = kExpNumBins/colorRange;
+
+    double gaussColorCoeff = -0.5 / (sigmaColor*sigmaColor);
+    double gaussSpaceCoeff = -0.5 / (sigmaSpace*sigmaSpace);
+
+    for (int i = 0; i < kExpNumBins + 2; i++)
+    {
+        double val = i / scaleIndex;
+        expLUT[i] = (float) std::exp(val * val * gaussColorCoeff);
+    }
+
+    Mat jointTemp, srcTemp;
+    copyMakeBorder(joint, jointTemp, radius, radius, radius, radius, borderType);
+    copyMakeBorder(src, srcTemp, radius, radius, radius, radius, borderType);
+    size_t srcElemStep = srcTemp.step / srcTemp.elemSize();
+    size_t jElemStep = jointTemp.step / jointTemp.elemSize();
+    CV_Assert(srcElemStep == jElemStep);
+
+    vector<float> spaceWeightsv(d*d);
+    vector<int> spaceOfsJointv(d*d);
+    float *spaceWeights = &spaceWeightsv[0];
+    int *spaceOfsJoint = &spaceOfsJointv[0];
+
+    int maxk = 0;
+    for (int i = -radius; i <= radius; i++)
+    {
+        for (int j = -radius; j <= radius; j++)
+        {
+            double r2 = i*i + j*j;
+            if (r2 > SQR(radius))
+                continue;
+
+            spaceWeights[maxk] = (float) std::exp(r2 * gaussSpaceCoeff);
+            spaceOfsJoint[maxk] = (int) (i*jElemStep + j);
+            maxk++;
+        }
+    }
+
+    Range range(0, joint.rows);
+    if (joint.type() == CV_32FC1)
+    {
+        if (src.type() == CV_32FC1)
+        {
+            parallel_for_(range, JointBilateralFilter_32f<Vec1f, Vec1f>(jointTemp, srcTemp, dst, radius, maxk, scaleIndex, spaceOfsJoint, spaceWeights, expLUT));
+        }
+        if (src.type() == CV_32FC3)
+        {
+            parallel_for_(range, JointBilateralFilter_32f<Vec1f, Vec3f>(jointTemp, srcTemp, dst, radius, maxk, scaleIndex, spaceOfsJoint, spaceWeights, expLUT));
+        }
+    }
+
+    if (joint.type() == CV_32FC3)
+    {
+        if (src.type() == CV_32FC1)
+        {
+            parallel_for_(range, JointBilateralFilter_32f<Vec3f, Vec1f>(jointTemp, srcTemp, dst, radius, maxk, scaleIndex, spaceOfsJoint, spaceWeights, expLUT));
+        }
+        if (src.type() == CV_32FC3)
+        {
+            parallel_for_(range, JointBilateralFilter_32f<Vec3f, Vec3f>(jointTemp, srcTemp, dst, radius, maxk, scaleIndex, spaceOfsJoint, spaceWeights, expLUT));
+        }
+    }
+}
+
+template<typename JointVec, typename SrcVec>
+class JointBilateralFilter_8u : public ParallelLoopBody
+{
+    Mat &joint, &src;
+    Mat &dst;
+    int radius, maxk;
+    float scaleIndex;
+    int *spaceOfs;
+    float *spaceWeights, *expLUT;
+
+public:
+
+    JointBilateralFilter_8u(Mat& joint_, Mat& src_, Mat& dst_, int radius_,
+        int maxk_, int *spaceOfs_, float *spaceWeights_, float *expLUT_)
+        :
+        joint(joint_), src(src_), dst(dst_), radius(radius_), maxk(maxk_),
+        spaceOfs(spaceOfs_), spaceWeights(spaceWeights_), expLUT(expLUT_)
+    {
+        CV_DbgAssert(joint.type() == JointVec::type && src.type() == dst.type() && src.type() == SrcVec::type);
+        CV_DbgAssert(joint.rows == src.rows && src.rows == dst.rows + 2 * radius);
+        CV_DbgAssert(joint.cols == src.cols && src.cols == dst.cols + 2 * radius);
+    }
+
+    void operator () (const Range& range) const
+    {
+        typedef Vec<int, JointVec::channels> JointVeci;
+        typedef Vec<float, SrcVec::channels> SrcVecf;
+
+        for (int i = radius + range.start; i < radius + range.end; i++)
+        {
+            for (int j = radius; j < src.cols - radius; j++)
+            {
+                JointVec *jointCenterPixPtr = joint.ptr<JointVec>(i) + j;
+                SrcVec *srcCenterPixPtr = src.ptr<SrcVec>(i) + j;
+
+                JointVeci jointPix0 = JointVeci(*jointCenterPixPtr);
+                SrcVecf srcSum = SrcVecf::all(0.0f);
+                float wSum = 0.0f;
+
+                for (int k = 0; k < maxk; k++)
+                {
+                    uchar *jointPix = reinterpret_cast<uchar*>(jointCenterPixPtr + spaceOfs[k]);
+                    int alpha = 0;
+                    for (int cn = 0; cn < JointVec::channels; cn++)
+                        alpha += std::abs(jointPix0[cn] - (int)jointPix[cn]);
+                    
+                    float weight = spaceWeights[k] * expLUT[alpha];
+
+                    uchar *srcPix = reinterpret_cast<uchar*>(srcCenterPixPtr + spaceOfs[k]);
+                    for (int cn = 0; cn < SrcVec::channels; cn++)
+                        srcSum[cn] += weight*srcPix[cn];
+                    wSum += weight;
+                }
+
+                dst.at<SrcVec>(i - radius, j - radius) = SrcVec(srcSum / wSum);
+            }
+        }
+    }
+};
+
+void jointBilateralFilter_8u(Mat& joint, Mat& src, Mat& dst, int radius, double sigmaColor, double sigmaSpace, int borderType)
+{
+    CV_DbgAssert(joint.depth() == CV_8U && src.depth() == CV_8U);
+
+    int d = 2 * radius + 1;
+    int jCn = joint.channels();
+
+    double gaussColorCoeff = -0.5 / (sigmaColor*sigmaColor);
+    double gaussSpaceCoeff = -0.5 / (sigmaSpace*sigmaSpace);
+
+    vector<float> expLUTv(jCn*0xFF);
+    float *expLUT = &expLUTv[0];
+
+    for (int i = 0; i < (int)expLUTv.size(); i++)
+    {
+        expLUT[i] = (float)std::exp(i * i * gaussColorCoeff);
+    }
+
+    Mat jointTemp, srcTemp;
+    copyMakeBorder(joint, jointTemp, radius, radius, radius, radius, borderType);
+    copyMakeBorder(src, srcTemp, radius, radius, radius, radius, borderType);
+    size_t srcElemStep = srcTemp.step / srcTemp.elemSize();
+    size_t jElemStep = jointTemp.step / jointTemp.elemSize();
+    CV_Assert(srcElemStep == jElemStep);
+
+    vector<float> spaceWeightsv(d*d);
+    vector<int> spaceOfsJointv(d*d);
+    float *spaceWeights = &spaceWeightsv[0];
+    int *spaceOfsJoint = &spaceOfsJointv[0];
+
+    int maxk = 0;
+    for (int i = -radius; i <= radius; i++)
+    {
+        for (int j = -radius; j <= radius; j++)
+        {
+            double r2 = i*i + j*j;
+            if (r2 > SQR(radius))
+                continue;
+
+            spaceWeights[maxk] = (float)std::exp(r2 * gaussSpaceCoeff);
+            spaceOfsJoint[maxk] = (int)(i*jElemStep + j);
+            maxk++;
+        }
+    }
+
+    Range range(0, src.rows);
+    if (joint.type() == CV_8UC1)
+    {
+        if (src.type() == CV_8UC1)
+        {
+            parallel_for_(range, JointBilateralFilter_8u<Vec1b, Vec1b>(jointTemp, srcTemp, dst, radius, maxk, spaceOfsJoint, spaceWeights, expLUT));
+        }
+        if (src.type() == CV_8UC3)
+        {
+            parallel_for_(range, JointBilateralFilter_8u<Vec1b, Vec3b>(jointTemp, srcTemp, dst, radius, maxk, spaceOfsJoint, spaceWeights, expLUT));
+        }
+    }
+
+    if (joint.type() == CV_8UC3)
+    {
+        if (src.type() == CV_8UC1)
+        {
+            parallel_for_(range, JointBilateralFilter_8u<Vec3b, Vec1b>(jointTemp, srcTemp, dst, radius, maxk, spaceOfsJoint, spaceWeights, expLUT));
+        }
+        if (src.type() == CV_8UC3)
+        {
+            parallel_for_(range, JointBilateralFilter_8u<Vec3b, Vec3b>(jointTemp, srcTemp, dst, radius, maxk, spaceOfsJoint, spaceWeights, expLUT));
+        }
+    }
+}
+
+void jointBilateralFilter(InputArray joint_, InputArray src_, OutputArray dst_, int d, double sigmaColor, double sigmaSpace, int borderType)
+{
+    CV_Assert(!src_.empty());
+
+    if (joint_.empty())
+    {
+        bilateralFilter(src_, dst_, d, sigmaColor, sigmaSpace, borderType);
+        return;
+    }
+
+    Mat src = src_.getMat();
+    Mat joint = joint_.getMat();
+
+    if (src.data == joint.data)
+    {
+        bilateralFilter(src_, dst_, d, sigmaColor, sigmaSpace, borderType);
+        return;
+    }
+
+    CV_Assert(src.size() == joint.size());
+    CV_Assert(src.depth() == joint.depth() && (src.depth() == CV_8U || src.depth() == CV_32F) );
+
+    if (sigmaColor <= 0)
+        sigmaColor = 1;
+    if (sigmaSpace <= 0)
+        sigmaSpace = 1;
+
+    int radius;
+    if (d <= 0)
+        radius = cvRound(sigmaSpace*1.5);
+    else
+        radius = d / 2;
+    radius = std::max(radius, 1);
+
+    dst_.create(src.size(), src.type());
+    Mat dst = dst_.getMat();
+
+    if (dst.data == joint.data)
+        joint = joint.clone();
+    if (dst.data == src.data)
+        src = src.clone();
+
+    int jointCnNum = joint.channels();
+    int srcCnNum = src.channels();
+
+    if ( (srcCnNum == 1 || srcCnNum == 3) && (jointCnNum == 1 || jointCnNum == 3) )
+    {
+        if (joint.depth() == CV_8U)
+        {
+            jointBilateralFilter_8u(joint, src, dst, radius, sigmaColor, sigmaSpace, borderType);
+        }
+        else
+        {
+            jointBilateralFilter_32f(joint, src, dst, radius, sigmaColor, sigmaSpace, borderType);
+        }
+    }
+    else
+    {
+        CV_Error(Error::BadNumChannels, "Unsupported number of channels");
+    }
+}
+
+}
+}

modules/ximgproc/src/precomp.hpp

@@ -0,0 +1,15 @@
+#ifndef _OPENCV_EDGEFILTER_PRECOMP_HPP_
+#define _OPENCV_EDGEFILTER_PRECOMP_HPP_
+
+#include <opencv2/core.hpp>
+#include <opencv2/core/ocl.hpp>
+#include <opencv2/core/base.hpp>
+#include <opencv2/core/utility.hpp>
+#include <opencv2/core/cvdef.h>
+#include <opencv2/core/core_c.h>
+#include <opencv2/core/private.hpp>
+#include <opencv2/imgproc.hpp>
+
+#include <opencv2/ximgproc.hpp>
+
+#endif

modules/ximgproc/test/test_adaptive_manifold.cpp

@@ -0,0 +1,174 @@
+#include "test_precomp.hpp"
+
+namespace cvtest
+{
+
+using namespace std;
+using namespace std::tr1;
+using namespace testing;
+using namespace cv;
+using namespace cv::ximgproc;
+
+#ifndef SQR
+#define SQR(x) ((x)*(x))
+#endif
+
+static string getOpenCVExtraDir()
+{
+    return cvtest::TS::ptr()->get_data_path();
+}
+
+static void checkSimilarity(InputArray res, InputArray ref)
+{
+    double normInf = cvtest::norm(res, ref, NORM_INF);
+    double normL2 = cvtest::norm(res, ref, NORM_L2) / res.total();
+
+    EXPECT_LE(normInf, 1);
+    EXPECT_LE(normL2, 1.0 / 64);
+}
+
+TEST(AdaptiveManifoldTest, SplatSurfaceAccuracy)
+{
+    RNG rnd(0);
+
+    for (int i = 0; i < 10; i++)
+    {
+        Size sz(rnd.uniform(512, 1024), rnd.uniform(512, 1024));
+
+        int guideCn = rnd.uniform(1, 8);
+        Mat guide(sz, CV_MAKE_TYPE(CV_32F, guideCn));
+        randu(guide, 0, 1);
+
+        Scalar surfaceValue;
+        int srcCn = rnd.uniform(1, 4);
+        rnd.fill(surfaceValue, RNG::UNIFORM, 0, 255);
+        Mat src(sz, CV_MAKE_TYPE(CV_8U, srcCn), surfaceValue);
+
+        double sigma_s = rnd.uniform(1.0, 50.0);
+        double sigma_r = rnd.uniform(0.1, 0.9);
+
+        Mat res;
+        amFilter(guide, src, res, sigma_s, sigma_r, false);
+
+        double normInf = cvtest::norm(src, res, NORM_INF);
+        EXPECT_EQ(normInf, 0);
+    }
+}
+
+TEST(AdaptiveManifoldTest, AuthorsReferenceAccuracy)
+{
+    String srcImgPath = "cv/edgefilter/kodim23.png";
+    
+    String refPaths[] =
+    {
+        "cv/edgefilter/amf/kodim23_amf_ss5_sr0.3_ref.png",
+        "cv/edgefilter/amf/kodim23_amf_ss30_sr0.1_ref.png",
+        "cv/edgefilter/amf/kodim23_amf_ss50_sr0.3_ref.png"
+    };
+
+    pair<double, double> refParams[] = 
+    {
+        make_pair(5.0, 0.3),
+        make_pair(30.0, 0.1),
+        make_pair(50.0, 0.3)
+    };
+
+    String refOutliersPaths[] = 
+    {
+        "cv/edgefilter/amf/kodim23_amf_ss5_sr0.1_outliers_ref.png",
+        "cv/edgefilter/amf/kodim23_amf_ss15_sr0.3_outliers_ref.png",
+        "cv/edgefilter/amf/kodim23_amf_ss50_sr0.5_outliers_ref.png"
+    };
+
+    pair<double, double> refOutliersParams[] =
+    {
+        make_pair(5.0, 0.1),
+        make_pair(15.0, 0.3),
+        make_pair(50.0, 0.5),
+    };
+
+    Mat srcImg = imread(getOpenCVExtraDir() + srcImgPath);
+    ASSERT_TRUE(!srcImg.empty());
+
+    for (int i = 0; i < 3; i++)
+    {
+        Mat refRes = imread(getOpenCVExtraDir() + refPaths[i]);
+        double sigma_s = refParams[i].first;
+        double sigma_r = refParams[i].second;
+        ASSERT_TRUE(!refRes.empty());
+
+        Mat res;
+        Ptr<AdaptiveManifoldFilter> amf = createAMFilter(sigma_s, sigma_r, false);
+        amf->setBool("use_RNG", false);
+        amf->filter(srcImg, res, srcImg);
+        amf->collectGarbage();
+
+        checkSimilarity(res, refRes);
+    }
+
+    for (int i = 0; i < 3; i++)
+    {
+        Mat refRes = imread(getOpenCVExtraDir() + refOutliersPaths[i]);
+        double sigma_s = refOutliersParams[i].first;
+        double sigma_r = refOutliersParams[i].second;
+        ASSERT_TRUE(!refRes.empty());
+
+        Mat res;
+        Ptr<AdaptiveManifoldFilter> amf = createAMFilter(sigma_s, sigma_r, true);
+        amf->setBool("use_RNG", false);
+        amf->filter(srcImg, res, srcImg);
+        amf->collectGarbage();
+
+        checkSimilarity(res, refRes);
+    }
+}
+
+typedef tuple<string, string> AMRefTestParams;
+typedef TestWithParam<AMRefTestParams> AdaptiveManifoldRefImplTest;
+
+Ptr<AdaptiveManifoldFilter> createAMFilterRefImpl(double sigma_s, double sigma_r, bool adjust_outliers = false);
+
+TEST_P(AdaptiveManifoldRefImplTest, RefImplAccuracy)
+{
+    AMRefTestParams params = GetParam();
+
+    string guideFileName = get<0>(params);
+    string srcFileName = get<1>(params);
+
+    Mat guide = imread(getOpenCVExtraDir() + guideFileName);
+    Mat src = imread(getOpenCVExtraDir() + srcFileName);
+    ASSERT_TRUE(!guide.empty() && !src.empty());
+
+    int seed = 10 * (int)guideFileName.length() + (int)srcFileName.length();
+    RNG rnd(seed);
+
+    //inconsistent downsample/upsample operations in reference implementation
+    Size dstSize((guide.cols + 15) & ~15, (guide.rows + 15) & ~15);
+
+    resize(guide, guide, dstSize);
+    resize(src, src, dstSize);
+
+    for (int iter = 0; iter < 6; iter++)
+    {
+        double sigma_s = rnd.uniform(1.0, 50.0);
+        double sigma_r = rnd.uniform(0.1, 0.9);
+        bool adjust_outliers = (iter % 2 == 0);
+
+        Mat res;
+        amFilter(guide, src, res, sigma_s, sigma_r, adjust_outliers);
+
+        Mat resRef;
+        Ptr<AdaptiveManifoldFilter> amf = createAMFilterRefImpl(sigma_s, sigma_r, adjust_outliers);
+        amf->filter(src, resRef, guide);
+
+        checkSimilarity(res, resRef);
+    }
+}
+
+INSTANTIATE_TEST_CASE_P(TypicalSet, AdaptiveManifoldRefImplTest, 
+    Combine(
+    Values("cv/shared/lena.png", "cv/edgefilter/kodim23.png", "cv/npr/test4.png"),
+    Values("cv/shared/lena.png", "cv/edgefilter/kodim23.png", "cv/npr/test4.png")
+));
+
+}

modules/ximgproc/test/test_adaptive_manifold_ref_impl.cpp

@@ -0,0 +1,948 @@
+/*
+ * The MIT License(MIT)
+ * 
+ * Copyright(c) 2013 Vladislav Vinogradov
+ * 
+ * Permission is hereby granted, free of charge, to any person obtaining a copy of
+ * this software and associated documentation files(the "Software"), to deal in
+ * the Software without restriction, including without limitation the rights to
+ * use, copy, modify, merge, publish, distribute, sublicense, and / or sell copies of
+ * the Software, and to permit persons to whom the Software is furnished to do so,
+ * subject to the following conditions :
+ * 
+ * The above copyright notice and this permission notice shall be included in all
+ * copies or substantial portions of the Software.
+ * 
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
+ * FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.IN NO EVENT SHALL THE AUTHORS OR
+ * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
+ * IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+ * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
+ */
+
+#include "test_precomp.hpp"
+#include <opencv2/core/private.hpp>
+#include <cmath>
+
+namespace
+{
+
+    using std::numeric_limits;
+    using namespace cv;
+    using namespace cv::ximgproc;
+
+    struct Buf
+    {
+        Mat_<Point3f> eta_1;
+        Mat_<uchar> cluster_1;
+
+        Mat_<Point3f> tilde_dst;
+        Mat_<float> alpha;
+        Mat_<Point3f> diff;
+        Mat_<Point3f> dst;
+
+        Mat_<float> V;
+
+        Mat_<Point3f> dIcdx;
+        Mat_<Point3f> dIcdy;
+        Mat_<float> dIdx;
+        Mat_<float> dIdy;
+        Mat_<float> dHdx;
+        Mat_<float> dVdy;
+
+        Mat_<float> t;
+
+        Mat_<float> theta_masked;
+        Mat_<Point3f> mul;
+        Mat_<Point3f> numerator;
+        Mat_<float> denominator;
+        Mat_<Point3f> numerator_filtered;
+        Mat_<float> denominator_filtered;
+
+        Mat_<Point3f> X;
+        Mat_<Point3f> eta_k_small;
+        Mat_<Point3f> eta_k_big;
+        Mat_<Point3f> X_squared;
+        Mat_<float> pixel_dist_to_manifold_squared;
+        Mat_<float> gaussian_distance_weights;
+        Mat_<Point3f> Psi_splat;
+        Mat_<Vec4f> Psi_splat_joined;
+        Mat_<Vec4f> Psi_splat_joined_resized;
+        Mat_<Vec4f> blurred_projected_values;
+        Mat_<Point3f> w_ki_Psi_blur;
+        Mat_<float> w_ki_Psi_blur_0;
+        Mat_<Point3f> w_ki_Psi_blur_resized;
+        Mat_<float> w_ki_Psi_blur_0_resized;
+        Mat_<float> rand_vec;
+        Mat_<float> v1;
+        Mat_<float> Nx_v1_mult;
+        Mat_<float> theta;
+
+        std::vector<Mat_<Point3f> > eta_minus;
+        std::vector<Mat_<uchar> > cluster_minus;
+        std::vector<Mat_<Point3f> > eta_plus;
+        std::vector<Mat_<uchar> > cluster_plus;
+
+        void release();
+    };
+
+    void Buf::release()
+    {
+        eta_1.release();
+        cluster_1.release();
+
+        tilde_dst.release();
+        alpha.release();
+        diff.release();
+        dst.release();
+
+        V.release();
+
+        dIcdx.release();
+        dIcdy.release();
+        dIdx.release();
+        dIdy.release();
+        dHdx.release();
+        dVdy.release();
+
+        t.release();
+
+        theta_masked.release();
+        mul.release();
+        numerator.release();
+        denominator.release();
+        numerator_filtered.release();
+        denominator_filtered.release();
+
+        X.release();
+        eta_k_small.release();
+        eta_k_big.release();
+        X_squared.release();
+        pixel_dist_to_manifold_squared.release();
+        gaussian_distance_weights.release();
+         Psi_splat.release();
+        Psi_splat_joined.release();
+        Psi_splat_joined_resized.release();
+        blurred_projected_values.release();
+        w_ki_Psi_blur.release();
+        w_ki_Psi_blur_0.release();
+        w_ki_Psi_blur_resized.release();
+        w_ki_Psi_blur_0_resized.release();
+        rand_vec.release();
+        v1.release();
+        Nx_v1_mult.release();
+        theta.release();
+
+        eta_minus.clear();
+        cluster_minus.clear();
+        eta_plus.clear();
+        cluster_plus.clear();
+    }
+
+    class AdaptiveManifoldFilterRefImpl : public AdaptiveManifoldFilter
+    {
+    public:
+        AlgorithmInfo* info() const;
+
+        AdaptiveManifoldFilterRefImpl();
+
+        void filter(InputArray src, OutputArray dst, InputArray joint);
+
+        void collectGarbage();
+
+    protected:
+        bool adjust_outliers_;
+        double sigma_s_;
+        double sigma_r_;
+        int tree_height_;
+        int num_pca_iterations_;
+        bool useRNG;
+
+    private:
+        void buildManifoldsAndPerformFiltering(const Mat_<Point3f>& eta_k, const Mat_<uchar>& cluster_k, int current_tree_level);
+
+        Buf buf_;
+
+        Mat_<Point3f> src_f_;
+        Mat_<Point3f> src_joint_f_;
+
+        Mat_<Point3f> sum_w_ki_Psi_blur_;
+        Mat_<float> sum_w_ki_Psi_blur_0_;
+
+        Mat_<float> min_pixel_dist_to_manifold_squared_;
+
+        RNG rng_;
+
+        int cur_tree_height_;
+        float sigma_r_over_sqrt_2_;
+    };
+
+    AdaptiveManifoldFilterRefImpl::AdaptiveManifoldFilterRefImpl()
+    {
+        sigma_s_ = 16.0;
+        sigma_r_ = 0.2;
+        tree_height_ = -1;
+        num_pca_iterations_ = 1;
+        adjust_outliers_ = false;
+        useRNG = true;
+    }
+
+    void AdaptiveManifoldFilterRefImpl::collectGarbage()
+    {
+        buf_.release();
+
+        src_f_.release();
+        src_joint_f_.release();
+
+        sum_w_ki_Psi_blur_.release();
+        sum_w_ki_Psi_blur_0_.release();
+
+        min_pixel_dist_to_manifold_squared_.release();
+    }
+
+    CV_INIT_ALGORITHM(AdaptiveManifoldFilterRefImpl, "AdaptiveManifoldFilterRefImpl",
+                      obj.info()->addParam(obj, "sigma_s", obj.sigma_s_, false, 0, 0, "Filter spatial standard deviation");
+                      obj.info()->addParam(obj, "sigma_r", obj.sigma_r_, false, 0, 0, "Filter range standard deviation");
+                      obj.info()->addParam(obj, "tree_height", obj.tree_height_, false, 0, 0, "Height of the manifold tree (default = -1 : automatically computed)");
+                      obj.info()->addParam(obj, "num_pca_iterations", obj.num_pca_iterations_, false, 0, 0, "Number of iterations to computed the eigenvector v1");
+                      obj.info()->addParam(obj, "adjust_outliers", obj.adjust_outliers_, false, 0, 0, "Specify supress outliers using Eq. 9 or not");
+                      obj.info()->addParam(obj, "use_RNG", obj.useRNG, false, 0, 0, "Specify use random to compute eigenvector or not");)
+
+    inline double Log2(double n)
+    {
+        return log(n) / log(2.0);
+    }
+
+    inline int computeManifoldTreeHeight(double sigma_s, double sigma_r)
+    {
+        const double Hs = floor(Log2(sigma_s)) - 1.0;
+        const double Lr = 1.0 - sigma_r;
+        return max(2, static_cast<int>(ceil(Hs * Lr)));
+    }
+
+    void ensureSizeIsEnough(int rows, int cols, int type, Mat& m)
+    {
+        if (m.empty() || m.type() != type || m.data != m.datastart)
+            m.create(rows, cols, type);
+        else
+        {
+            const size_t esz = m.elemSize();
+            const ptrdiff_t delta2 = m.dataend - m.datastart;
+
+            const size_t minstep = m.cols * esz;
+
+            Size wholeSize;
+            wholeSize.height = std::max(static_cast<int>((delta2 - minstep) / m.step + 1), m.rows);
+            wholeSize.width = std::max(static_cast<int>((delta2 - m.step * (wholeSize.height - 1)) / esz), m.cols);
+
+            if (wholeSize.height < rows || wholeSize.width < cols)
+                m.create(rows, cols, type);
+            else
+            {
+                m.cols = cols;
+                m.rows = rows;
+            }
+        }
+    }
+
+    inline void ensureSizeIsEnough(Size size, int type, Mat& m)
+    {
+        ensureSizeIsEnough(size.height, size.width, type, m);
+    }
+
+    template <typename T>
+    inline void ensureSizeIsEnough(int rows, int cols, Mat_<T>& m)
+    {
+        if (m.empty() || m.data != m.datastart)
+            m.create(rows, cols);
+        else
+        {
+            const size_t esz = m.elemSize();
+            const ptrdiff_t delta2 = m.dataend - m.datastart;
+
+            const size_t minstep = m.cols * esz;
+
+            Size wholeSize;
+            wholeSize.height = std::max(static_cast<int>((delta2 - minstep) / m.step + 1), m.rows);
+            wholeSize.width = std::max(static_cast<int>((delta2 - m.step * (wholeSize.height - 1)) / esz), m.cols);
+
+            if (wholeSize.height < rows || wholeSize.width < cols)
+                m.create(rows, cols);
+            else
+            {
+                m.cols = cols;
+                m.rows = rows;
+            }
+        }
+    }
+
+    template <typename T>
+    inline void ensureSizeIsEnough(Size size, Mat_<T>& m)
+    {
+        ensureSizeIsEnough(size.height, size.width, m);
+    }
+
+    template <typename T>
+    void h_filter(const Mat_<T>& src, Mat_<T>& dst, float sigma)
+    {
+        CV_DbgAssert( src.depth() == CV_32F );
+
+        const float a = exp(-sqrt(2.0f) / sigma);
+
+        ensureSizeIsEnough(src.size(), dst);
+
+        for (int y = 0; y < src.rows; ++y)
+        {
+            const T* src_row = src[y];
+            T* dst_row = dst[y];
+
+            dst_row[0] = src_row[0];
+            for (int x = 1; x < src.cols; ++x)
+            {
+                //dst_row[x] = src_row[x] + a * (src_row[x - 1] - src_row[x]);
+                dst_row[x] = src_row[x] + a * (dst_row[x - 1] - src_row[x]); //!!!
+            }
+            for (int x = src.cols - 2; x >= 0; --x)
+            {
+                dst_row[x] = dst_row[x] + a * (dst_row[x + 1] - dst_row[x]);
+            }
+        }
+
+        for (int y = 1; y < src.rows; ++y)
+        {
+            T* dst_cur_row = dst[y];
+            T* dst_prev_row = dst[y - 1];
+
+            for (int x = 0; x < src.cols; ++x)
+            {
+                dst_cur_row[x] = dst_cur_row[x] + a * (dst_prev_row[x] - dst_cur_row[x]);
+            }
+        }
+        for (int y = src.rows - 2; y >= 0; --y)
+        {
+            T* dst_cur_row = dst[y];
+            T* dst_prev_row = dst[y + 1];
+
+            for (int x = 0; x < src.cols; ++x)
+            {
+                dst_cur_row[x] = dst_cur_row[x] + a * (dst_prev_row[x] - dst_cur_row[x]);
+            }
+        }
+    }
+
+    template <typename T>
+    void rdivide(const Mat_<T>& a, const Mat_<float>& b, Mat_<T>& dst)
+    {
+        CV_DbgAssert( a.depth() == CV_32F );
+        CV_DbgAssert( a.size() == b.size() );
+
+        ensureSizeIsEnough(a.size(), dst);
+        dst.setTo(0);
+
+        for (int y = 0; y < a.rows; ++y)
+        {
+            const T* a_row = a[y];
+            const float* b_row = b[y];
+            T* dst_row = dst[y];
+
+            for (int x = 0; x < a.cols; ++x)
+            {
+                //if (b_row[x] > numeric_limits<float>::epsilon())
+                    dst_row[x] = a_row[x] * (1.0f / b_row[x]);
+            }
+        }
+    }
+
+    template <typename T>
+    void times(const Mat_<T>& a, const Mat_<float>& b, Mat_<T>& dst)
+    {
+        CV_DbgAssert( a.depth() == CV_32F );
+        CV_DbgAssert( a.size() == b.size() );
+
+        ensureSizeIsEnough(a.size(), dst);
+
+        for (int y = 0; y < a.rows; ++y)
+        {
+            const T* a_row = a[y];
+            const float* b_row = b[y];
+            T* dst_row = dst[y];
+
+            for (int x = 0; x < a.cols; ++x)
+            {
+                dst_row[x] = a_row[x] * b_row[x];
+            }
+        }
+    }
+
+    void AdaptiveManifoldFilterRefImpl::filter(InputArray _src, OutputArray _dst, InputArray _joint)
+    {
+        const Mat src = _src.getMat();
+        const Mat src_joint = _joint.getMat();
+
+        const Size srcSize = src.size();
+
+        CV_Assert( src.type() == CV_8UC3 );
+        CV_Assert( src_joint.empty() || (src_joint.type() == src.type() && src_joint.size() == srcSize) );
+
+        ensureSizeIsEnough(srcSize, src_f_);
+        src.convertTo(src_f_, src_f_.type(), 1.0 / 255.0);
+
+        ensureSizeIsEnough(srcSize, sum_w_ki_Psi_blur_);
+        sum_w_ki_Psi_blur_.setTo(Scalar::all(0));
+
+        ensureSizeIsEnough(srcSize, sum_w_ki_Psi_blur_0_);
+        sum_w_ki_Psi_blur_0_.setTo(Scalar::all(0));
+
+        ensureSizeIsEnough(srcSize, min_pixel_dist_to_manifold_squared_);
+        min_pixel_dist_to_manifold_squared_.setTo(Scalar::all(numeric_limits<float>::max()));
+
+        // If the tree_height was not specified, compute it using Eq. (10) of our paper.
+        cur_tree_height_ = tree_height_ > 0 ? tree_height_ : computeManifoldTreeHeight(sigma_s_, sigma_r_);
+
+        // If no joint signal was specified, use the original signal
+        ensureSizeIsEnough(srcSize, src_joint_f_);
+        if (src_joint.empty())
+            src_f_.copyTo(src_joint_f_);
+        else
+            src_joint.convertTo(src_joint_f_, src_joint_f_.type(), 1.0 / 255.0);
+
+        // Use the center pixel as seed to random number generation.
+        const double seedCoef = src_joint_f_(src_joint_f_.rows / 2, src_joint_f_.cols / 2).x;
+        const uint64 baseCoef = numeric_limits<uint64>::max() / 0xFFFF;
+        rng_.state = static_cast<uint64>(baseCoef*seedCoef);
+
+        // Dividing the covariance matrix by 2 is equivalent to dividing the standard deviations by sqrt(2).
+        sigma_r_over_sqrt_2_ = static_cast<float>(sigma_r_ / sqrt(2.0));
+
+        // Algorithm 1, Step 1: compute the first manifold by low-pass filtering.
+        h_filter(src_joint_f_, buf_.eta_1, static_cast<float>(sigma_s_));
+
+        ensureSizeIsEnough(srcSize, buf_.cluster_1);
+        buf_.cluster_1.setTo(Scalar::all(1));
+
+        buf_.eta_minus.resize(cur_tree_height_);
+        buf_.cluster_minus.resize(cur_tree_height_);
+        buf_.eta_plus.resize(cur_tree_height_);
+        buf_.cluster_plus.resize(cur_tree_height_);
+        buildManifoldsAndPerformFiltering(buf_.eta_1, buf_.cluster_1, 1);
+
+        // Compute the filter response by normalized convolution -- Eq. (4)
+        rdivide(sum_w_ki_Psi_blur_, sum_w_ki_Psi_blur_0_, buf_.tilde_dst);
+
+        if (!adjust_outliers_)
+        {
+            buf_.tilde_dst.convertTo(_dst, CV_8U, 255.0);
+        }
+        else
+        {
+            // Adjust the filter response for outlier pixels -- Eq. (10)
+            ensureSizeIsEnough(srcSize, buf_.alpha);
+            exp(min_pixel_dist_to_manifold_squared_ * (-0.5 / sigma_r_ / sigma_r_), buf_.alpha);
+
+            ensureSizeIsEnough(srcSize, buf_.diff);
+            subtract(buf_.tilde_dst, src_f_, buf_.diff);
+            times(buf_.diff, buf_.alpha, buf_.diff);
+
+            ensureSizeIsEnough(srcSize, buf_.dst);
+            add(src_f_, buf_.diff, buf_.dst);
+
+            buf_.dst.convertTo(_dst, CV_8U, 255.0);
+        }
+    }
+
+    inline double floor_to_power_of_two(double r)
+    {
+        return pow(2.0, floor(Log2(r)));
+    }
+
+    void channelsSum(const Mat_<Point3f>& src, Mat_<float>& dst)
+    {
+        ensureSizeIsEnough(src.size(), dst);
+
+        for (int y = 0; y < src.rows; ++y)
+        {
+            const Point3f* src_row = src[y];
+            float* dst_row = dst[y];
+
+            for (int x = 0; x < src.cols; ++x)
+            {
+                const Point3f src_val = src_row[x];
+                dst_row[x] = src_val.x + src_val.y + src_val.z;
+            }
+        }
+    }
+
+    void phi(const Mat_<float>& src, Mat_<float>& dst, float sigma)
+    {
+        ensureSizeIsEnough(src.size(), dst);
+
+        for (int y = 0; y < dst.rows; ++y)
+        {
+            const float* src_row = src[y];
+            float* dst_row = dst[y];
+
+            for (int x = 0; x < dst.cols; ++x)
+            {
+                dst_row[x] = exp(-0.5f * src_row[x] / sigma / sigma);
+            }
+        }
+    }
+
+    void catCn(const Mat_<Point3f>& a, const Mat_<float>& b, Mat_<Vec4f>& dst)
+    {
+        ensureSizeIsEnough(a.size(), dst);
+
+        for (int y = 0; y < a.rows; ++y)
+        {
+            const Point3f* a_row = a[y];
+            const float* b_row = b[y];
+            Vec4f* dst_row = dst[y];
+
+            for (int x = 0; x < a.cols; ++x)
+            {
+                const Point3f a_val = a_row[x];
+                const float b_val = b_row[x];
+                dst_row[x] = Vec4f(a_val.x, a_val.y, a_val.z, b_val);
+            }
+        }
+    }
+
+    void diffY(const Mat_<Point3f>& src, Mat_<Point3f>& dst)
+    {
+        ensureSizeIsEnough(src.rows - 1, src.cols, dst);
+
+        for (int y = 0; y < src.rows - 1; ++y)
+        {
+            const Point3f* src_cur_row = src[y];
+            const Point3f* src_next_row = src[y + 1];
+            Point3f* dst_row = dst[y];
+
+            for (int x = 0; x < src.cols; ++x)
+            {
+                dst_row[x] = src_next_row[x] - src_cur_row[x];
+            }
+        }
+    }
+
+    void diffX(const Mat_<Point3f>& src, Mat_<Point3f>& dst)
+    {
+        ensureSizeIsEnough(src.rows, src.cols - 1, dst);
+
+        for (int y = 0; y < src.rows; ++y)
+        {
+            const Point3f* src_row = src[y];
+            Point3f* dst_row = dst[y];
+
+            for (int x = 0; x < src.cols - 1; ++x)
+            {
+                dst_row[x] = src_row[x + 1] - src_row[x];
+            }
+        }
+    }
+
+    void TransformedDomainRecursiveFilter(const Mat_<Vec4f>& I, const Mat_<float>& DH, const Mat_<float>& DV, Mat_<Vec4f>& dst, float sigma, Buf& buf)
+    {
+        CV_DbgAssert( I.size() == DH.size() );
+
+        const float a = exp(-sqrt(2.0f) / sigma);
+
+        ensureSizeIsEnough(I.size(), dst);
+        I.copyTo(dst);
+
+        ensureSizeIsEnough(DH.size(), buf.V);
+
+        for (int y = 0; y < DH.rows; ++y)
+        {
+            const float* D_row = DH[y];
+            float* V_row = buf.V[y];
+
+            for (int x = 0; x < DH.cols; ++x)
+            {
+                V_row[x] = pow(a, D_row[x]);
+            }
+        }
+        for (int y = 0; y < I.rows; ++y)
+        {
+            const float* V_row = buf.V[y];
+            Vec4f* dst_row = dst[y];
+
+            for (int x = 1; x < I.cols; ++x)
+            {
+                Vec4f dst_cur_val = dst_row[x];
+                const Vec4f dst_prev_val = dst_row[x - 1];
+                const float V_val = V_row[x];
+
+                dst_cur_val[0] += V_val * (dst_prev_val[0] - dst_cur_val[0]);
+                dst_cur_val[1] += V_val * (dst_prev_val[1] - dst_cur_val[1]);
+                dst_cur_val[2] += V_val * (dst_prev_val[2] - dst_cur_val[2]);
+                dst_cur_val[3] += V_val * (dst_prev_val[3] - dst_cur_val[3]);
+
+                dst_row[x] = dst_cur_val;
+            }
+            for (int x = I.cols - 2; x >= 0; --x)
+            {
+                Vec4f dst_cur_val = dst_row[x];
+                const Vec4f dst_prev_val = dst_row[x + 1];
+                //const float V_val = V_row[x];
+                const float V_val = V_row[x+1];
+
+                dst_cur_val[0] += V_val * (dst_prev_val[0] - dst_cur_val[0]);
+                dst_cur_val[1] += V_val * (dst_prev_val[1] - dst_cur_val[1]);
+                dst_cur_val[2] += V_val * (dst_prev_val[2] - dst_cur_val[2]);
+                dst_cur_val[3] += V_val * (dst_prev_val[3] - dst_cur_val[3]);
+
+                dst_row[x] = dst_cur_val;
+            }
+        }
+
+        for (int y = 0; y < DV.rows; ++y)
+        {
+            const float* D_row = DV[y];
+            float* V_row = buf.V[y];
+
+            for (int x = 0; x < DV.cols; ++x)
+            {
+                V_row[x] = pow(a, D_row[x]);
+            }
+        }
+        for (int y = 1; y < I.rows; ++y)
+        {
+            const float* V_row = buf.V[y];
+            Vec4f* dst_cur_row = dst[y];
+            Vec4f* dst_prev_row = dst[y - 1];
+
+            for (int x = 0; x < I.cols; ++x)
+            {
+                Vec4f dst_cur_val = dst_cur_row[x];
+                const Vec4f dst_prev_val = dst_prev_row[x];
+                const float V_val = V_row[x];
+
+                dst_cur_val[0] += V_val * (dst_prev_val[0] - dst_cur_val[0]);
+                dst_cur_val[1] += V_val * (dst_prev_val[1] - dst_cur_val[1]);
+                dst_cur_val[2] += V_val * (dst_prev_val[2] - dst_cur_val[2]);
+                dst_cur_val[3] += V_val * (dst_prev_val[3] - dst_cur_val[3]);
+
+                dst_cur_row[x] = dst_cur_val;
+            }
+        }
+        for (int y = I.rows - 2; y >= 0; --y)
+        {
+            //const float* V_row = buf.V[y];
+            const float* V_row = buf.V[y + 1];
+            Vec4f* dst_cur_row = dst[y];
+            Vec4f* dst_prev_row = dst[y + 1];
+
+            for (int x = 0; x < I.cols; ++x)
+            {
+                Vec4f dst_cur_val = dst_cur_row[x];
+                const Vec4f dst_prev_val = dst_prev_row[x];
+                const float V_val = V_row[x];
+
+                dst_cur_val[0] += V_val * (dst_prev_val[0] - dst_cur_val[0]);
+                dst_cur_val[1] += V_val * (dst_prev_val[1] - dst_cur_val[1]);
+                dst_cur_val[2] += V_val * (dst_prev_val[2] - dst_cur_val[2]);
+                dst_cur_val[3] += V_val * (dst_prev_val[3] - dst_cur_val[3]);
+
+                dst_cur_row[x] = dst_cur_val;
+            }
+        }
+    }
+
+    void RF_filter(const Mat_<Vec4f>& src, const Mat_<Point3f>& src_joint, Mat_<Vec4f>& dst, float sigma_s, float sigma_r, Buf& buf)
+    {
+        CV_DbgAssert( src_joint.size() == src.size() );
+
+        diffX(src_joint, buf.dIcdx);
+        diffY(src_joint, buf.dIcdy);
+
+        ensureSizeIsEnough(src.size(), buf.dIdx);
+        buf.dIdx.setTo(Scalar::all(0));
+        for (int y = 0; y < src.rows; ++y)
+        {
+            const Point3f* dIcdx_row = buf.dIcdx[y];
+            float* dIdx_row = buf.dIdx[y];
+
+            for (int x = 1; x < src.cols; ++x)
+            {
+                const Point3f val = dIcdx_row[x - 1];
+                dIdx_row[x] = val.dot(val);
+            }
+        }
+
+        ensureSizeIsEnough(src.size(), buf.dIdy);
+        buf.dIdy.setTo(Scalar::all(0));
+        for (int y = 1; y < src.rows; ++y)
+        {
+            const Point3f* dIcdy_row = buf.dIcdy[y - 1];
+            float* dIdy_row = buf.dIdy[y];
+
+            for (int x = 0; x < src.cols; ++x)
+            {
+                const Point3f val = dIcdy_row[x];
+                dIdy_row[x] = val.dot(val);
+            }
+        }
+
+        ensureSizeIsEnough(buf.dIdx.size(), buf.dHdx);
+        buf.dIdx.convertTo(buf.dHdx, buf.dHdx.type(), (sigma_s / sigma_r) * (sigma_s / sigma_r), (sigma_s / sigma_s) * (sigma_s / sigma_s));
+        sqrt(buf.dHdx, buf.dHdx);
+
+        ensureSizeIsEnough(buf.dIdy.size(), buf.dVdy);
+        buf.dIdy.convertTo(buf.dVdy, buf.dVdy.type(), (sigma_s / sigma_r) * (sigma_s / sigma_r), (sigma_s / sigma_s) * (sigma_s / sigma_s));
+        sqrt(buf.dVdy, buf.dVdy);
+
+        ensureSizeIsEnough(src.size(), dst);
+        src.copyTo(dst);
+        TransformedDomainRecursiveFilter(src, buf.dHdx, buf.dVdy, dst, sigma_s, buf);
+    }
+
+    void split_3_1(const Mat_<Vec4f>& src, Mat_<Point3f>& dst1, Mat_<float>& dst2)
+    {
+        ensureSizeIsEnough(src.size(), dst1);
+        ensureSizeIsEnough(src.size(), dst2);
+
+        for (int y = 0; y < src.rows; ++y)
+        {
+            const Vec4f* src_row = src[y];
+            Point3f* dst1_row = dst1[y];
+            float* dst2_row = dst2[y];
+
+            for (int x = 0; x < src.cols; ++x)
+            {
+                Vec4f val = src_row[x];
+                dst1_row[x] = Point3f(val[0], val[1], val[2]);
+                dst2_row[x] = val[3];
+            }
+        }
+    }
+
+    void computeEigenVector(const Mat_<float>& X, const Mat_<uchar>& mask, Mat_<float>& dst, int num_pca_iterations, const Mat_<float>& rand_vec, Buf& buf)
+    {
+        CV_DbgAssert( X.cols == rand_vec.cols );
+        CV_DbgAssert( X.rows == mask.size().area() );
+        CV_DbgAssert( rand_vec.rows == 1 );
+
+        ensureSizeIsEnough(rand_vec.size(), dst);
+        rand_vec.copyTo(dst);
+
+        ensureSizeIsEnough(X.size(), buf.t);
+
+        float* dst_row = dst[0];
+
+        for (int i = 0; i < num_pca_iterations; ++i)
+        {
+            buf.t.setTo(Scalar::all(0));
+
+            for (int y = 0, ind = 0; y < mask.rows; ++y)
+            {
+                const uchar* mask_row = mask[y];
+
+                for (int x = 0; x < mask.cols; ++x, ++ind)
+                {
+                    if (mask_row[x])
+                    {
+                        const float* X_row = X[ind];
+                        float* t_row = buf.t[ind];
+
+                        float dots = 0.0;
+                        for (int c = 0; c < X.cols; ++c)
+                            dots += dst_row[c] * X_row[c];
+
+                        for (int c = 0; c < X.cols; ++c)
+                            t_row[c] = dots * X_row[c];
+                    }
+                }
+            }
+
+            dst.setTo(0.0);
+            for (int k = 0; k < X.rows; ++k)
+            {
+                const float* t_row = buf.t[k];
+
+                for (int c = 0; c < X.cols; ++c)
+                {
+                    dst_row[c] += t_row[c];
+                }
+            }
+        }
+
+        double n = norm(dst);
+        divide(dst, n, dst);
+    }
+
+    void calcEta(const Mat_<Point3f>& src_joint_f, const Mat_<float>& theta, const Mat_<uchar>& cluster, Mat_<Point3f>& dst, float sigma_s, float df, Buf& buf)
+    {
+        ensureSizeIsEnough(theta.size(), buf.theta_masked);
+        buf.theta_masked.setTo(Scalar::all(0));
+        theta.copyTo(buf.theta_masked, cluster);
+
+        times(src_joint_f, buf.theta_masked, buf.mul);
+
+        const Size nsz = Size(saturate_cast<int>(buf.mul.cols * (1.0 / df)), saturate_cast<int>(buf.mul.rows * (1.0 / df)));
+
+        ensureSizeIsEnough(nsz, buf.numerator);
+        resize(buf.mul, buf.numerator, Size(), 1.0 / df, 1.0 / df);
+
+        ensureSizeIsEnough(nsz, buf.denominator);
+        resize(buf.theta_masked, buf.denominator, Size(), 1.0 / df, 1.0 / df);        
+        h_filter(buf.numerator, buf.numerator_filtered, sigma_s / df);
+        h_filter(buf.denominator, buf.denominator_filtered, sigma_s / df);
+
+
+        rdivide(buf.numerator_filtered, buf.denominator_filtered, dst);
+    }
+
+    void AdaptiveManifoldFilterRefImpl::buildManifoldsAndPerformFiltering(const Mat_<Point3f>& eta_k, const Mat_<uchar>& cluster_k, int current_tree_level)
+    {
+        // Compute downsampling factor
+
+        double df = min(sigma_s_ / 4.0, 256.0 * sigma_r_);
+        df = floor_to_power_of_two(df);
+        df = max(1.0, df);
+
+        // Splatting: project the pixel values onto the current manifold eta_k
+
+        if (eta_k.rows == src_joint_f_.rows)
+        {
+            ensureSizeIsEnough(src_joint_f_.size(), buf_.X);
+            subtract(src_joint_f_, eta_k, buf_.X);
+
+            const Size nsz = Size(saturate_cast<int>(eta_k.cols * (1.0 / df)), saturate_cast<int>(eta_k.rows * (1.0 / df)));
+            ensureSizeIsEnough(nsz, buf_.eta_k_small);
+            resize(eta_k, buf_.eta_k_small, Size(), 1.0 / df, 1.0 / df);
+        }
+        else
+        {
+            ensureSizeIsEnough(eta_k.size(), buf_.eta_k_small);
+            eta_k.copyTo(buf_.eta_k_small);
+
+            ensureSizeIsEnough(src_joint_f_.size(), buf_.eta_k_big);
+            resize(eta_k, buf_.eta_k_big, src_joint_f_.size());
+
+            ensureSizeIsEnough(src_joint_f_.size(), buf_.X);
+            subtract(src_joint_f_, buf_.eta_k_big, buf_.X);
+        }
+
+        // Project pixel colors onto the manifold -- Eq. (3), Eq. (5)
+
+        ensureSizeIsEnough(buf_.X.size(), buf_.X_squared);
+        multiply(buf_.X, buf_.X, buf_.X_squared);
+
+        channelsSum(buf_.X_squared, buf_.pixel_dist_to_manifold_squared);
+
+        phi(buf_.pixel_dist_to_manifold_squared, buf_.gaussian_distance_weights, sigma_r_over_sqrt_2_);
+
+        times(src_f_, buf_.gaussian_distance_weights, buf_.Psi_splat);
+
+        const Mat_<float>& Psi_splat_0 = buf_.gaussian_distance_weights;
+
+        // Save min distance to later perform adjustment of outliers -- Eq. (10)
+
+        if (adjust_outliers_)
+        {
+            cv::min(_InputArray(min_pixel_dist_to_manifold_squared_), _InputArray(buf_.pixel_dist_to_manifold_squared), _OutputArray(min_pixel_dist_to_manifold_squared_));
+        }
+
+        // Blurring: perform filtering over the current manifold eta_k
+
+        catCn(buf_.Psi_splat, Psi_splat_0, buf_.Psi_splat_joined);
+
+        ensureSizeIsEnough(buf_.eta_k_small.size(), buf_.Psi_splat_joined_resized);
+        resize(buf_.Psi_splat_joined, buf_.Psi_splat_joined_resized, buf_.eta_k_small.size());
+
+        RF_filter(buf_.Psi_splat_joined_resized, buf_.eta_k_small, buf_.blurred_projected_values, static_cast<float>(sigma_s_ / df), sigma_r_over_sqrt_2_, buf_);
+
+        split_3_1(buf_.blurred_projected_values, buf_.w_ki_Psi_blur, buf_.w_ki_Psi_blur_0);
+
+        // Slicing: gather blurred values from the manifold
+
+        // Since we perform splatting and slicing at the same points over the manifolds,
+        // the interpolation weights are equal to the gaussian weights used for splatting.
+        const Mat_<float>& w_ki = buf_.gaussian_distance_weights;
+
+        ensureSizeIsEnough(src_f_.size(), buf_.w_ki_Psi_blur_resized);
+        resize(buf_.w_ki_Psi_blur, buf_.w_ki_Psi_blur_resized, src_f_.size());
+        times(buf_.w_ki_Psi_blur_resized, w_ki, buf_.w_ki_Psi_blur_resized);
+        add(sum_w_ki_Psi_blur_, buf_.w_ki_Psi_blur_resized, sum_w_ki_Psi_blur_);
+
+        ensureSizeIsEnough(src_f_.size(), buf_.w_ki_Psi_blur_0_resized);
+        resize(buf_.w_ki_Psi_blur_0, buf_.w_ki_Psi_blur_0_resized, src_f_.size());
+        times(buf_.w_ki_Psi_blur_0_resized, w_ki, buf_.w_ki_Psi_blur_0_resized);
+        add(sum_w_ki_Psi_blur_0_, buf_.w_ki_Psi_blur_0_resized, sum_w_ki_Psi_blur_0_);
+
+        // Compute two new manifolds eta_minus and eta_plus
+
+        if (current_tree_level < cur_tree_height_)
+        {
+            // Algorithm 1, Step 2: compute the eigenvector v1
+            const Mat_<float> nX(src_joint_f_.size().area(), 3, (float*) buf_.X.data);
+
+            ensureSizeIsEnough(1, nX.cols, buf_.rand_vec);
+            if (useRNG)
+            {
+                rng_.fill(buf_.rand_vec, RNG::UNIFORM, -0.5, 0.5);
+            }
+            else
+            {
+                for (int i = 0; i < (int)buf_.rand_vec.total(); i++)
+                    buf_.rand_vec(0, i) = (i % 2 == 0) ? 0.5f : -0.5f;
+            }
+
+            computeEigenVector(nX, cluster_k, buf_.v1, num_pca_iterations_, buf_.rand_vec, buf_);
+
+            // Algorithm 1, Step 3: Segment pixels into two clusters -- Eq. (6)
+
+            ensureSizeIsEnough(nX.rows, buf_.v1.rows, buf_.Nx_v1_mult);
+            gemm(nX, buf_.v1, 1.0, noArray(), 0.0, buf_.Nx_v1_mult, GEMM_2_T);
+
+            const Mat_<float> dot(src_joint_f_.rows, src_joint_f_.cols, (float*) buf_.Nx_v1_mult.data);
+
+            Mat_<uchar>& cluster_minus = buf_.cluster_minus[current_tree_level];
+            ensureSizeIsEnough(dot.size(), cluster_minus);
+            compare(dot, 0, cluster_minus, CMP_LT);
+            bitwise_and(cluster_minus, cluster_k, cluster_minus);
+
+            Mat_<uchar>& cluster_plus = buf_.cluster_plus[current_tree_level];
+            ensureSizeIsEnough(dot.size(), cluster_plus);
+            //compare(dot, 0, cluster_plus, CMP_GT);
+            compare(dot, 0, cluster_plus, CMP_GE);
+            bitwise_and(cluster_plus, cluster_k, cluster_plus);
+
+            // Algorithm 1, Step 4: Compute new manifolds by weighted low-pass filtering -- Eq. (7-8)
+
+            ensureSizeIsEnough(w_ki.size(), buf_.theta);
+            buf_.theta.setTo(Scalar::all(1.0));
+            subtract(buf_.theta, w_ki, buf_.theta);
+
+            Mat_<Point3f>& eta_minus = buf_.eta_minus[current_tree_level];
+            calcEta(src_joint_f_, buf_.theta, cluster_minus, eta_minus, (float)sigma_s_, (float)df, buf_);
+
+            Mat_<Point3f>& eta_plus = buf_.eta_plus[current_tree_level];
+            calcEta(src_joint_f_, buf_.theta, cluster_plus, eta_plus, (float)sigma_s_, (float)df, buf_);
+
+            // Algorithm 1, Step 5: recursively build more manifolds.
+
+            buildManifoldsAndPerformFiltering(eta_minus, cluster_minus, current_tree_level + 1);
+            buildManifoldsAndPerformFiltering(eta_plus, cluster_plus, current_tree_level + 1);
+        }
+    }
+}
+
+namespace cvtest
+{
+
+using namespace cv::ximgproc;
+
+Ptr<AdaptiveManifoldFilter> createAMFilterRefImpl(double sigma_s, double sigma_r, bool adjust_outliers)
+{
+    Ptr<AdaptiveManifoldFilter> amf(new AdaptiveManifoldFilterRefImpl());
+    
+    amf->set("sigma_s", sigma_s);
+    amf->set("sigma_r", sigma_r);
+    amf->set("adjust_outliers", adjust_outliers);
+
+    return amf;
+}
+
+}

modules/ximgproc/test/test_domain_transform.cpp

@@ -0,0 +1,220 @@
+#include "test_precomp.hpp"
+
+namespace cvtest
+{
+
+using namespace std;
+using namespace std::tr1;
+using namespace testing;
+using namespace perf;
+using namespace cv;
+using namespace cv::ximgproc;
+
+static string getOpenCVExtraDir()
+{
+    return cvtest::TS::ptr()->get_data_path();
+}
+
+CV_ENUM(SupportedTypes, CV_8UC1, CV_8UC2, CV_8UC3, CV_8UC4, CV_32FC1, CV_32FC2, CV_32FC3, CV_32FC4);
+CV_ENUM(ModeType, DTF_NC, DTF_IC, DTF_RF)
+typedef tuple<Size, ModeType, SupportedTypes, SupportedTypes> DTParams;
+
+Mat convertTypeAndSize(Mat src, int dstType, Size dstSize)
+{
+    Mat dst;
+    CV_Assert(src.channels() == 3);
+
+    int dstChannels = CV_MAT_CN(dstType);
+
+    if (dstChannels == 1)
+    {
+        cvtColor(src, dst, COLOR_BGR2GRAY);
+    }
+    else if (dstChannels == 2)
+    {
+        Mat srcCn[3];
+        split(src, srcCn);
+        merge(srcCn, 2, dst);
+    }
+    else if (dstChannels == 3)
+    {
+        dst = src.clone();
+    }
+    else if (dstChannels == 4)
+    {
+        Mat srcCn[4];
+        split(src, srcCn);
+        srcCn[3] = srcCn[0].clone();
+        merge(srcCn, 4, dst);
+    }
+
+    dst.convertTo(dst, dstType);
+    resize(dst, dst, dstSize);
+
+    return dst;
+}
+
+TEST(DomainTransformTest, SplatSurfaceAccuracy)
+{
+    static int dtModes[] = {DTF_NC, DTF_RF, DTF_IC};
+    RNG rnd(0);
+
+    for (int i = 0; i < 15; i++)
+    {
+        Size sz(rnd.uniform(512, 1024), rnd.uniform(512, 1024));
+
+        int guideCn = rnd.uniform(1, 4);
+        Mat guide(sz, CV_MAKE_TYPE(CV_32F, guideCn));
+        randu(guide, 0, 255);
+
+        Scalar surfaceValue;
+        int srcCn = rnd.uniform(1, 4);
+        rnd.fill(surfaceValue, RNG::UNIFORM, 0, 255);
+        Mat src(sz, CV_MAKE_TYPE(CV_8U, srcCn), surfaceValue);
+
+        double sigma_s = rnd.uniform(1.0, 100.0);
+        double sigma_r = rnd.uniform(1.0, 100.0);
+        int mode = dtModes[i%3];
+
+        Mat res;
+        dtFilter(guide, src, res, sigma_s, sigma_r, mode, 1);
+
+        double normL1 = cvtest::norm(src, res, NORM_L1)/src.total()/src.channels();
+        EXPECT_LE(normL1, 1.0/64);
+    }
+}
+
+typedef TestWithParam<DTParams> DomainTransformTest;
+TEST_P(DomainTransformTest, MultiThreadReproducibility)
+{
+    if (cv::getNumberOfCPUs() == 1)
+        return;
+
+    double MAX_DIF = 1.0;
+    double MAX_MEAN_DIF = 1.0 / 256.0;
+    int loopsCount = 2;
+    RNG rng(0);
+
+    DTParams params = GetParam();
+    Size size = get<0>(params);
+    int mode = get<1>(params);
+    int guideType = get<2>(params);
+    int srcType = get<3>(params);
+
+    Mat original = imread(getOpenCVExtraDir() + "cv/edgefilter/statue.png");
+    Mat guide = convertTypeAndSize(original, guideType, size);
+    Mat src = convertTypeAndSize(original, srcType, size);
+
+    for (int iter = 0; iter <= loopsCount; iter++)
+    {
+        double ss = rng.uniform(0.0, 100.0);
+        double sc = rng.uniform(0.0, 100.0);
+
+        cv::setNumThreads(cv::getNumberOfCPUs());
+        Mat resMultithread;
+        dtFilter(guide, src, resMultithread, ss, sc, mode);
+
+        cv::setNumThreads(1);
+        Mat resSingleThread;
+        dtFilter(guide, src, resSingleThread, ss, sc, mode);
+
+        EXPECT_LE(cv::norm(resSingleThread, resMultithread, NORM_INF), MAX_DIF);
+        EXPECT_LE(cv::norm(resSingleThread, resMultithread, NORM_L1), MAX_MEAN_DIF*src.total());
+    }
+}
+
+INSTANTIATE_TEST_CASE_P(FullSet, DomainTransformTest,
+    Combine(Values(szODD, szQVGA), ModeType::all(), SupportedTypes::all(), SupportedTypes::all())
+);
+
+template<typename SrcVec>
+Mat getChessMat1px(Size sz, double whiteIntensity = 255)
+{
+    typedef typename DataType<SrcVec>::channel_type SrcType;
+
+    Mat dst(sz, DataType<SrcVec>::type);
+
+    SrcVec black = SrcVec::all(0);
+    SrcVec white = SrcVec::all((SrcType)whiteIntensity);
+
+    for (int i = 0; i < dst.rows; i++)
+        for (int j = 0; j < dst.cols; j++)
+            dst.at<SrcVec>(i, j) = ((i + j) % 2) ? white : black;
+
+    return dst;
+}
+
+TEST(DomainTransformTest, ChessBoard_NC_accuracy)
+{
+    RNG rng(0);
+    double MAX_DIF = 1;
+    Size sz = szVGA;
+    double ss = 80;
+    double sc = 60;
+
+    Mat srcb = randomMat(rng, sz, CV_8UC4, 0, 255, true);
+    Mat srcf = randomMat(rng, sz, CV_32FC4, 0, 255, true);
+    Mat chessb = getChessMat1px<Vec3b>(sz);
+
+    Mat dstb, dstf;
+    dtFilter(chessb, srcb.clone(), dstb, ss, sc, DTF_NC);
+    dtFilter(chessb, srcf.clone(), dstf, ss, sc, DTF_NC);
+
+    EXPECT_LE(cv::norm(srcb, dstb, NORM_INF), MAX_DIF);
+    EXPECT_LE(cv::norm(srcf, dstf, NORM_INF), MAX_DIF);
+}
+
+TEST(DomainTransformTest, BoxFilter_NC_accuracy)
+{
+    double MAX_DIF = 1;
+    int radius = 5;
+    double sc = 1.0;
+    double ss = 1.01*radius / sqrt(3.0);
+
+    Mat src = imread(getOpenCVExtraDir() + "cv/edgefilter/statue.png");
+    ASSERT_TRUE(!src.empty());
+
+    Mat1b guide(src.size(), 200);
+    Mat res_dt, res_box;
+
+    blur(src, res_box, Size(2 * radius + 1, 2 * radius + 1));
+    dtFilter(guide, src, res_dt, ss, sc, DTF_NC, 1);
+
+    EXPECT_LE(cv::norm(res_dt, res_box, NORM_L2), MAX_DIF*src.total());
+}
+
+TEST(DomainTransformTest, AuthorReferenceAccuracy)
+{
+    string dir = getOpenCVExtraDir() + "cv/edgefilter";
+    double ss = 30;
+    double sc = 0.2 * 255;
+
+    Mat src = imread(dir + "/statue.png");
+    Mat ref_NC = imread(dir + "/dt/authors_statue_NC_ss30_sc0.2.png");
+    Mat ref_IC = imread(dir + "/dt/authors_statue_IC_ss30_sc0.2.png");
+    Mat ref_RF = imread(dir + "/dt/authors_statue_RF_ss30_sc0.2.png");
+
+    ASSERT_FALSE(src.empty());
+    ASSERT_FALSE(ref_NC.empty());
+    ASSERT_FALSE(ref_IC.empty());
+    ASSERT_FALSE(ref_RF.empty());
+
+    cv::setNumThreads(cv::getNumberOfCPUs());
+    Mat res_NC, res_IC, res_RF;
+    dtFilter(src, src, res_NC, ss, sc, DTF_NC);
+    dtFilter(src, src, res_IC, ss, sc, DTF_IC);
+    dtFilter(src, src, res_RF, ss, sc, DTF_RF);
+
+    double totalMaxError = 1.0/64.0*src.total();
+    
+    EXPECT_LE(cvtest::norm(res_NC, ref_NC, NORM_L2), totalMaxError);
+    EXPECT_LE(cvtest::norm(res_NC, ref_NC, NORM_INF), 1);
+
+    EXPECT_LE(cvtest::norm(res_IC, ref_IC, NORM_L2), totalMaxError);
+    EXPECT_LE(cvtest::norm(res_IC, ref_IC, NORM_INF), 1);
+
+    EXPECT_LE(cvtest::norm(res_RF, ref_RF, NORM_L2), totalMaxError);
+    EXPECT_LE(cvtest::norm(res_IC, ref_IC, NORM_INF), 1);
+}
+
+}

modules/ximgproc/test/test_guided_filter.cpp

@@ -0,0 +1,362 @@
+#include "test_precomp.hpp"
+
+namespace cvtest
+{
+
+using namespace std;
+using namespace std::tr1;
+using namespace testing;
+using namespace cv;
+using namespace cv::ximgproc;
+
+#ifndef SQR
+#define SQR(x) ((x)*(x))
+#endif
+
+static string getOpenCVExtraDir()
+{
+    return cvtest::TS::ptr()->get_data_path();
+}
+
+static Mat convertTypeAndSize(Mat src, int dstType, Size dstSize)
+{
+    Mat dst;
+    int srcCnNum = src.channels();
+    int dstCnNum = CV_MAT_CN(dstType);
+    CV_Assert(srcCnNum == 3);
+
+    if (srcCnNum == dstCnNum)
+    {
+        src.copyTo(dst);
+    }
+    else if (dstCnNum == 1 && srcCnNum == 3)
+    {
+        cvtColor(src, dst, COLOR_BGR2GRAY);
+    }
+    else if (dstCnNum == 1 && srcCnNum == 4)
+    {
+        cvtColor(src, dst, COLOR_BGRA2GRAY);
+    }
+    else
+    {
+        vector<Mat> srcCn;
+        split(src, srcCn);
+        srcCn.resize(dstCnNum);
+
+        uint64 seed = 10000 * src.rows + 1000 * src.cols + 100 * dstSize.height + 10 * dstSize.width + dstType;
+        RNG rnd(seed);
+
+        for (int i = srcCnNum; i < dstCnNum; i++)
+        {
+            Mat& donor = srcCn[i % srcCnNum];
+
+            double minVal, maxVal;
+            minMaxLoc(donor, &minVal, &maxVal);
+
+            Mat randItem(src.size(), CV_MAKE_TYPE(src.depth(), 1));
+            randn(randItem, 0, (maxVal - minVal) / 100);
+
+            add(donor, randItem, srcCn[i]);
+        }
+
+        merge(srcCn, dst);
+    }
+
+    dst.convertTo(dst, dstType);
+    resize(dst, dst, dstSize);
+
+    return dst;
+}
+
+class GuidedFilterRefImpl : public GuidedFilter
+{
+    int height, width, rad, chNum;
+    Mat det;
+    Mat *channels, *exps, **vars, **A;
+    double eps;
+
+    void meanFilter(const Mat &src, Mat & dst);
+
+    void computeCovGuide();
+
+    void computeCovGuideInv();
+
+    void applyTransform(int cNum, Mat *Ichannels, Mat *beta, Mat **alpha, int dDepth);
+
+    void computeCovGuideAndSrc(int cNum, Mat **vars_I, Mat *Ichannels, Mat *exp_I);
+
+    void computeBeta(int cNum, Mat *beta, Mat *exp_I, Mat **alpha);
+
+    void computeAlpha(int cNum, Mat **alpha, Mat **vars_I);
+
+public:
+
+    GuidedFilterRefImpl(InputArray guide_, int rad, double eps);
+
+    void filter(InputArray src, OutputArray dst, int dDepth = -1);
+
+    ~GuidedFilterRefImpl();
+};
+
+void GuidedFilterRefImpl::meanFilter(const Mat &src, Mat & dst)
+{
+    boxFilter(src, dst, CV_32F, Size(2 * rad + 1, 2 * rad + 1), Point(-1, -1), true, BORDER_REFLECT);
+}
+
+GuidedFilterRefImpl::GuidedFilterRefImpl(InputArray _guide, int _rad, double _eps) :
+  height(_guide.rows()), width(_guide.cols()), rad(_rad), chNum(_guide.channels()), eps(_eps)
+{
+    Mat guide = _guide.getMat();
+    CV_Assert(chNum > 0 && chNum <= 3);
+
+    channels = new Mat[chNum];
+    exps     = new Mat[chNum];
+
+    A    = new Mat *[chNum];
+    vars = new Mat *[chNum];
+    for (int i = 0; i < chNum; ++i)
+    {
+        A[i]    = new Mat[chNum];
+        vars[i] = new Mat[chNum];
+    }
+
+    split(guide, channels);
+    for (int i = 0; i < chNum; ++i)
+    {
+        channels[i].convertTo(channels[i], CV_32F);
+        meanFilter(channels[i], exps[i]);
+    }
+
+    computeCovGuide();
+
+    computeCovGuideInv();
+}
+
+void GuidedFilterRefImpl::computeCovGuide()
+{
+    static const int pY[] = { 0, 0, 1, 0, 1, 2 };
+    static const int pX[] = { 0, 1, 1, 2, 2, 2 };
+
+    int numOfIterations = (SQR(chNum) - chNum) / 2 + chNum;
+    for (int k = 0; k < numOfIterations; ++k)
+    {
+        int i = pY[k], j = pX[k];
+
+        vars[i][j] = channels[i].mul(channels[j]);
+        meanFilter(vars[i][j], vars[i][j]);
+        vars[i][j] -= exps[i].mul(exps[j]);
+
+        if (i == j)
+            vars[i][j] += eps * Mat::ones(height, width, CV_32F);
+        else
+            vars[j][i] = vars[i][j];
+    }
+}
+
+void GuidedFilterRefImpl::computeCovGuideInv()
+{
+    static const int pY[] = { 0, 0, 1, 0, 1, 2 };
+    static const int pX[] = { 0, 1, 1, 2, 2, 2 };
+
+    int numOfIterations = (SQR(chNum) - chNum) / 2 + chNum;
+    if (chNum == 3)
+    {
+        for (int k = 0; k < numOfIterations; ++k){
+            int i = pY[k], i1 = (pY[k] + 1) % 3, i2 = (pY[k] + 2) % 3;
+            int j = pX[k], j1 = (pX[k] + 1) % 3, j2 = (pX[k] + 2) % 3;
+
+            A[i][j] = vars[i1][j1].mul(vars[i2][j2])
+                - vars[i1][j2].mul(vars[i2][j1]);
+        }
+    }
+    else if (chNum == 2)
+    {
+        A[0][0] = vars[1][1];
+        A[1][1] = vars[0][0];
+        A[0][1] = -vars[0][1];
+    }
+    else if (chNum == 1)
+        A[0][0] = Mat::ones(height, width, CV_32F);
+
+    for (int i = 0; i < chNum; ++i)
+        for (int j = 0; j < i; ++j)
+            A[i][j] = A[j][i];
+
+    det = vars[0][0].mul(A[0][0]);
+    for (int k = 0; k < chNum - 1; ++k)
+        det += vars[0][k + 1].mul(A[0][k + 1]);
+}
+
+GuidedFilterRefImpl::~GuidedFilterRefImpl(){
+    delete [] channels;
+    delete [] exps;
+
+    for (int i = 0; i < chNum; ++i)
+    {
+        delete [] A[i];
+        delete [] vars[i];
+    }
+
+    delete [] A;
+    delete [] vars;
+}
+
+void GuidedFilterRefImpl::filter(InputArray src_, OutputArray dst_, int dDepth)
+{
+    if (dDepth == -1) dDepth = src_.depth();
+    dst_.create(height, width, src_.type());
+    Mat src = src_.getMat();
+    Mat dst = dst_.getMat();
+    int cNum = src.channels();
+
+    CV_Assert(height == src.rows && width == src.cols);
+
+    Mat *Ichannels, *exp_I, **vars_I, **alpha, *beta;
+    Ichannels = new Mat[cNum];
+    exp_I     = new Mat[cNum];
+    beta      = new Mat[cNum];
+
+    vars_I = new Mat *[chNum];
+    alpha  = new Mat *[chNum];
+    for (int i = 0; i < chNum; ++i){
+        vars_I[i] = new Mat[cNum];
+        alpha[i]  = new Mat[cNum];
+    }
+
+    split(src, Ichannels);
+    for (int i = 0; i < cNum; ++i)
+    {
+        Ichannels[i].convertTo(Ichannels[i], CV_32F);
+        meanFilter(Ichannels[i], exp_I[i]);
+    }
+
+    computeCovGuideAndSrc(cNum, vars_I, Ichannels, exp_I);
+
+    computeAlpha(cNum, alpha, vars_I);
+
+    computeBeta(cNum, beta, exp_I, alpha);
+
+    for (int i = 0; i < chNum + 1; ++i)
+        for (int j = 0; j < cNum; ++j)
+            if (i < chNum)
+                meanFilter(alpha[i][j], alpha[i][j]);
+            else
+                meanFilter(beta[j], beta[j]);
+
+    applyTransform(cNum, Ichannels, beta, alpha, dDepth);
+    merge(Ichannels, cNum, dst);
+
+    delete [] Ichannels;
+    delete [] exp_I;
+    delete [] beta;
+
+    for (int i = 0; i < chNum; ++i)
+    {
+        delete [] vars_I[i];
+        delete [] alpha[i];
+    }
+    delete [] vars_I;
+    delete [] alpha;
+}
+
+void GuidedFilterRefImpl::computeAlpha(int cNum, Mat **alpha, Mat **vars_I)
+{
+    for (int i = 0; i < chNum; ++i)
+        for (int j = 0; j < cNum; ++j)
+        {
+            alpha[i][j] = vars_I[0][j].mul(A[i][0]);
+            for (int k = 1; k < chNum; ++k)
+                alpha[i][j] += vars_I[k][j].mul(A[i][k]);
+            alpha[i][j] /= det;
+        }
+}
+
+void GuidedFilterRefImpl::computeBeta(int cNum, Mat *beta, Mat *exp_I, Mat **alpha)
+{
+    for (int i = 0; i < cNum; ++i)
+    {
+        beta[i] = exp_I[i];
+        for (int j = 0; j < chNum; ++j)
+            beta[i] -= alpha[j][i].mul(exps[j]);
+    }
+}
+
+void GuidedFilterRefImpl::computeCovGuideAndSrc(int cNum, Mat **vars_I, Mat *Ichannels, Mat *exp_I)
+{
+    for (int i = 0; i < chNum; ++i)
+        for (int j = 0; j < cNum; ++j)
+        {
+            vars_I[i][j] = channels[i].mul(Ichannels[j]);
+            meanFilter(vars_I[i][j], vars_I[i][j]);
+            vars_I[i][j] -= exp_I[j].mul(exps[i]);
+        }
+}
+
+void GuidedFilterRefImpl::applyTransform(int cNum, Mat *Ichannels, Mat *beta, Mat **alpha, int dDepth)
+{
+    for (int i = 0; i < cNum; ++i)
+    {
+        Ichannels[i] = beta[i];
+        for (int j = 0; j < chNum; ++j)
+            Ichannels[i] += alpha[j][i].mul(channels[j]);
+        Ichannels[i].convertTo(Ichannels[i], dDepth);
+    }
+}
+
+typedef tuple<int, int, string, string> GFParams;
+typedef TestWithParam<GFParams> GuidedFilterTest;
+
+TEST_P(GuidedFilterTest, accuracy)
+{
+    GFParams params = GetParam();
+
+    int guideCnNum = get<0>(params);
+    int srcCnNum = get<1>(params);
+
+    string guideFileName = get<2>(params);
+    string srcFileName = get<3>(params);
+
+    int seed = 100 * guideCnNum + 50 * srcCnNum + 5*(int)guideFileName.length() + (int)srcFileName.length();
+    RNG rng(seed);
+
+    Mat guide = imread(getOpenCVExtraDir() + guideFileName);
+    Mat src = imread(getOpenCVExtraDir() + srcFileName);
+    ASSERT_TRUE(!guide.empty() && !src.empty());
+
+    Size dstSize(guide.cols + 1 + rng.uniform(0, 3), guide.rows);
+
+    guide = convertTypeAndSize(guide, CV_MAKE_TYPE(guide.depth(), guideCnNum), dstSize);
+    src = convertTypeAndSize(src, CV_MAKE_TYPE(src.depth(), srcCnNum), dstSize);
+
+    for (int iter = 0; iter < 3; iter++)
+    {
+        int radius = rng.uniform(0, 50);
+        double eps = rng.uniform(0.0, SQR(255.0));
+
+        cv::setNumThreads(cv::getNumberOfCPUs());
+        Mat res;
+        Ptr<GuidedFilter> gf = createGuidedFilter(guide, radius, eps);
+        gf->filter(src, res);
+
+        cv::setNumThreads(1);
+        Mat resRef;
+        Ptr<GuidedFilter> gfRef(new GuidedFilterRefImpl(guide, radius, eps));
+        gfRef->filter(src, resRef);
+
+        double normInf = cv::norm(res, resRef, NORM_INF);
+        double normL2 = cv::norm(res, resRef, NORM_L2) / guide.total();
+
+        EXPECT_LE(normInf, 1.0);
+        EXPECT_LE(normL2, 1.0/64.0);
+    }
+}
+
+INSTANTIATE_TEST_CASE_P(TypicalSet, GuidedFilterTest, 
+    Combine(
+    Values(1, 2, 3),
+    Values(1, 2, 3), 
+    Values("cv/shared/lena.png", "cv/shared/baboon.png", "cv/npr/test2.png"),
+    Values("cv/shared/lena.png", "cv/shared/baboon.png", "cv/npr/test2.png")
+));
+
+}

modules/ximgproc/test/test_joint_bilateral_filter.cpp

@@ -0,0 +1,256 @@
+#include "test_precomp.hpp"
+
+namespace cvtest
+{
+
+using namespace std;
+using namespace std::tr1;
+using namespace testing;
+using namespace cv;
+using namespace cv::ximgproc;
+
+static std::string getOpenCVExtraDir()
+{
+    return cvtest::TS::ptr()->get_data_path();
+}
+
+static void checkSimilarity(InputArray src, InputArray ref)
+{
+    double normInf = cvtest::norm(src, ref, NORM_INF);
+    double normL2 = cvtest::norm(src, ref, NORM_L2) / (src.total()*src.channels());
+
+    EXPECT_LE(normInf, 1.0);
+    EXPECT_LE(normL2, 1.0 / 16);
+}
+
+static Mat convertTypeAndSize(Mat src, int dstType, Size dstSize)
+{
+    Mat dst;
+    int srcCnNum = src.channels();
+    int dstCnNum = CV_MAT_CN(dstType);
+
+    if (srcCnNum == dstCnNum)
+    {
+        src.copyTo(dst);
+    }
+    else if (srcCnNum == 3 && dstCnNum == 1)
+    {
+        cvtColor(src, dst, COLOR_BGR2GRAY);
+    }
+    else if (srcCnNum == 1 && dstCnNum == 3)
+    {
+        cvtColor(src, dst, COLOR_GRAY2BGR);
+    }
+    else
+    {
+        CV_Error(Error::BadNumChannels, "Bad num channels in src");
+    }
+
+    dst.convertTo(dst, dstType);
+    resize(dst, dst, dstSize);
+
+    return dst;
+}
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+void jointBilateralFilterNaive(InputArray joint, InputArray src, OutputArray dst, int d, double sigmaColor, double sigmaSpace, int borderType = BORDER_DEFAULT);
+
+typedef Vec<float, 1> Vec1f;
+typedef Vec<uchar, 1> Vec1b;
+
+#ifndef SQR
+#define SQR(a) ((a)*(a))
+#endif
+
+template<typename T, int cn>
+float normL1Sqr(const Vec<T, cn>& a, const Vec<T, cn>& b)
+{
+    float res = 0.0f;
+    for (int i = 0; i < cn; i++)
+        res += std::abs((float)a[i] - (float)b[i]);
+    return res*res;
+}
+
+template<typename JointVec, typename SrcVec>
+void jointBilateralFilterNaive_(InputArray joint_, InputArray src_, OutputArray dst_, int d, double sigmaColor, double sigmaSpace, int borderType)
+{
+    CV_Assert(joint_.size() == src_.size());
+    CV_Assert(joint_.type() == JointVec::type && src_.type() == SrcVec::type);
+    typedef Vec<float, SrcVec::channels> SrcVecf;
+
+    if (sigmaColor <= 0)
+        sigmaColor = 1;
+    if (sigmaSpace <= 0)
+        sigmaSpace = 1;
+
+    int radius;
+    if (d <= 0)
+        radius = cvRound(sigmaSpace*1.5);
+    else
+        radius = d / 2;
+    radius = std::max(radius, 1);
+    d = 2 * radius + 1;
+
+    dst_.create(src_.size(), src_.type());
+    Mat_<SrcVec> dst = dst_.getMat();
+    Mat_<JointVec> jointExt;
+    Mat_<SrcVec> srcExt;
+    copyMakeBorder(src_, srcExt, radius, radius, radius, radius, borderType);
+    copyMakeBorder(joint_, jointExt, radius, radius, radius, radius, borderType);
+
+    float colorGaussCoef = (float)(-0.5 / (sigmaColor*sigmaColor));
+    float spaceGaussCoef = (float)(-0.5 / (sigmaSpace*sigmaSpace));
+
+    for (int i = radius; i < srcExt.rows - radius; i++)
+    {
+        for (int j = radius; j < srcExt.cols - radius; j++)
+        {
+            JointVec joint0 = jointExt(i, j);
+            SrcVecf sum = SrcVecf::all(0.0f);
+            float sumWeights = 0.0f;
+
+            for (int k = -radius; k <= radius; k++)
+            {
+                for (int l = -radius; l <= radius; l++)
+                {
+                    float spatDistSqr = (float)(k*k + l*l);
+                    if (spatDistSqr > SQR(radius)) continue;
+
+                    float colorDistSqr = normL1Sqr(joint0, jointExt(i + k, j + l));
+
+                    float weight = std::exp(spatDistSqr*spaceGaussCoef + colorDistSqr*colorGaussCoef);
+
+                    sum += weight*SrcVecf(srcExt(i + k, j + l));
+                    sumWeights += weight;
+                }
+            }
+
+            dst(i - radius, j - radius) = sum / sumWeights;
+        }
+    }
+}
+
+void jointBilateralFilterNaive(InputArray joint, InputArray src, OutputArray dst, int d, double sigmaColor, double sigmaSpace, int borderType)
+{
+    CV_Assert(src.size() == joint.size() && src.depth() == joint.depth());
+    CV_Assert(src.type() == CV_32FC1 || src.type() == CV_32FC3 || src.type() == CV_8UC1 || src.type() == CV_8UC3);
+    CV_Assert(joint.type() == CV_32FC1 || joint.type() == CV_32FC3 || joint.type() == CV_8UC1 || joint.type() == CV_8UC3);
+
+    int jointType = joint.type();
+    int srcType = src.type();
+
+    #define JBF_naive(VecJoint, VecSrc) jointBilateralFilterNaive_<VecJoint, VecSrc>(joint, src, dst, d, sigmaColor, sigmaSpace, borderType);
+    if (jointType == CV_8UC1)
+    {
+        if (srcType == CV_8UC1)  JBF_naive(Vec1b, Vec1b);
+        if (srcType == CV_8UC3)  JBF_naive(Vec1b, Vec3b);
+    }
+    if (jointType == CV_8UC3)
+    {
+        if (srcType == CV_8UC1)  JBF_naive(Vec3b, Vec1b);
+        if (srcType == CV_8UC3)  JBF_naive(Vec3b, Vec3b);
+    }
+    if (jointType == CV_32FC1)
+    {
+        if (srcType == CV_32FC1) JBF_naive(Vec1f, Vec1f);
+        if (srcType == CV_32FC3) JBF_naive(Vec1f, Vec3f);
+    }
+    if (jointType == CV_32FC3)
+    {
+        if (srcType == CV_32FC1) JBF_naive(Vec3f, Vec1f);
+        if (srcType == CV_32FC3) JBF_naive(Vec3f, Vec3f);
+    }
+    #undef JBF_naive
+}
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+typedef tuple<double, string, string, int, int, int> JBFTestParam;
+typedef TestWithParam<JBFTestParam> JointBilateralFilterTest_NaiveRef;
+
+TEST_P(JointBilateralFilterTest_NaiveRef, Accuracy)
+{
+    JBFTestParam param = GetParam();
+    double sigmaS       = get<0>(param);
+    string jointPath    = get<1>(param);
+    string srcPath      = get<2>(param);
+    int depth           = get<3>(param);
+    int jCn             = get<4>(param);
+    int srcCn           = get<5>(param);
+    int jointType       = CV_MAKE_TYPE(depth, jCn);
+    int srcType         = CV_MAKE_TYPE(depth, srcCn);
+
+    Mat joint = imread(getOpenCVExtraDir() + jointPath);
+    Mat src = imread(getOpenCVExtraDir() + srcPath);
+    ASSERT_TRUE(!joint.empty() && !src.empty());
+
+    joint = convertTypeAndSize(joint, jointType, joint.size());
+    src = convertTypeAndSize(src, srcType, joint.size());
+
+    RNG rnd(cvRound(10*sigmaS) + jointType + srcType + jointPath.length() + srcPath.length() + joint.rows + joint.cols);
+    double sigmaC = rnd.uniform(0, 255);
+
+    Mat resNaive;
+    jointBilateralFilterNaive(joint, src, resNaive, 0, sigmaC, sigmaS);
+
+    cv::setNumThreads(cv::getNumberOfCPUs());
+    Mat res;
+    jointBilateralFilter(joint, src, res, 0, sigmaC, sigmaS);
+    cv::setNumThreads(1);
+
+    checkSimilarity(res, resNaive);
+}
+
+INSTANTIATE_TEST_CASE_P(Set2, JointBilateralFilterTest_NaiveRef,
+    Combine(
+    Values(4.0, 6.0, 8.0),
+    Values("/cv/shared/airplane.png", "/cv/shared/fruits.png"),
+    Values("/cv/shared/airplane.png", "/cv/shared/lena.png", "/cv/shared/fruits.png"),
+    Values(CV_8U, CV_32F),
+    Values(1, 3),
+    Values(1, 3))
+);
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+typedef tuple<double, string, int> BFTestParam;
+typedef TestWithParam<BFTestParam> JointBilateralFilterTest_BilateralRef;
+
+TEST_P(JointBilateralFilterTest_BilateralRef, Accuracy)
+{
+    BFTestParam param   = GetParam();
+    double sigmaS       = get<0>(param);
+    string srcPath      = get<1>(param);
+    int srcType         = get<2>(param);
+
+    Mat src = imread(getOpenCVExtraDir() + srcPath);
+    ASSERT_TRUE(!src.empty());
+    src = convertTypeAndSize(src, srcType, src.size());
+
+    RNG rnd(cvRound(10*sigmaS) + srcPath.length() + srcType + src.rows);
+    double sigmaC = rnd.uniform(0.0, 255.0);
+
+    cv::setNumThreads(cv::getNumberOfCPUs());
+
+    Mat resRef;
+    bilateralFilter(src, resRef, 0, sigmaC, sigmaS);
+
+    Mat res, joint = src.clone();
+    jointBilateralFilter(joint, src, res, 0, sigmaC, sigmaS);
+
+    checkSimilarity(res, resRef);
+}
+
+INSTANTIATE_TEST_CASE_P(Set1, JointBilateralFilterTest_BilateralRef,
+    Combine(
+    Values(4.0, 6.0, 8.0),
+    Values("/cv/shared/pic2.png", "/cv/shared/lena.png", "cv/shared/box_in_scene.png"),
+    Values(CV_8UC1, CV_8UC3, CV_32FC1, CV_32FC3)
+    )
+);
+
+}

modules/ximgproc/test/test_main.cpp

@@ -0,0 +1,3 @@
+#include "test_precomp.hpp"
+
+CV_TEST_MAIN("")

modules/ximgproc/test/test_precomp.hpp

@@ -0,0 +1,20 @@
+#ifdef __GNUC__
+#  pragma GCC diagnostic ignored "-Wmissing-declarations"
+#  if defined __clang__ || defined __APPLE__
+#    pragma GCC diagnostic ignored "-Wmissing-prototypes"
+#    pragma GCC diagnostic ignored "-Wextra"
+#  endif
+#endif
+
+#ifndef __OPENCV_TEST_PRECOMP_HPP__
+#define __OPENCV_TEST_PRECOMP_HPP__
+
+#include <opencv2/ts.hpp>
+#include <opencv2/ts/ts_perf.hpp>
+#include <opencv2/core.hpp>
+#include <opencv2/core/utility.hpp>
+#include <opencv2/imgproc.hpp>
+#include <opencv2/highgui.hpp>
+#include <opencv2/ximgproc.hpp>
+
+#endif