Added joint bilateral filter realization.

modules/ximgproc/doc/edge_aware_filters.rst

@@ -36,7 +36,7 @@ Produce domain transform filtering operation on source image.
 
 .. ocv:pyfunction:: cv2.DTFilter.filter(src, dst[, dDepth]) -> None
 
-    :param src: filtering image with unsigned 8-bit or floating 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 dst: destination image.
     :param dDepth: optional depth of the output image. ``dDepth`` can be set to -1, which will be equivalent to ``src.depth()``.
     
@@ -49,8 +49,8 @@ If you have multiple images to filter with the same guided image then use :ocv:c
 
 .. 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 32-bit depth and up to 4 channels.
-    :param src: filtering image with unsigned 8-bit or floating 32-bit depth and up to 4 channels.
+    :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.
@@ -204,6 +204,33 @@ Simple one-line Adaptive Manifold Filter call.
     
 .. 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
 ==========
 
@@ -220,8 +247,11 @@ References
 ..  [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/include/opencv2/edge_filter.hpp

@@ -116,6 +116,12 @@ CV_EXPORTS_W void amFilter(InputArray joint, InputArray src, OutputArray dst, do
 //////////////////////////////////////////////////////////////////////////
 //////////////////////////////////////////////////////////////////////////
 
+CV_EXPORTS_W
+void jointBilateralFilter(InputArray joint, InputArray src, OutputArray dst, int d, double sigmaColor, double sigmaSpace, int borderType = BORDER_DEFAULT);
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
 /*Interface for Weighted Median Filter*/
 class CV_EXPORTS_W WeightedMedian : public Algorithm
 {

modules/ximgproc/perf/pref_joint_bilateral_filter.cpp

+#include "perf_precomp.hpp"
+
+using namespace std;
+using namespace std::tr1;
+using namespace cv;
+using namespace perf;
+using namespace testing;
+
+namespace cvtest
+{
+
+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(sigmaS + sz.height + depth + jCn + srcCn);
+    double sigmaC = rnd.uniform(0.0, 255.0);
+
+    TEST_CYCLE_N(1)
+    {
+        jointBilateralFilter(joint, src, dst, 0, sigmaC, sigmaS);
+    }
+    
+    SANITY_CHECK(dst);
+}
+}
\ No newline at end of file

modules/ximgproc/src/adaptive_manifold_filter_n.cpp

@@ -643,7 +643,7 @@ void AdaptiveManifoldFilterN::computeClusters(Mat1b& cluster, Mat1b& cluster_min
     }
     else
     {
-        for (int i = 0; i < initVec.total(); i++)
+        for (int i = 0; i < (int)initVec.total(); i++)
             initVec(0, i) = (i % 2 == 0) ? 0.5f : -0.5f;
     }
 

modules/ximgproc/src/joint_bilateral_filter.cpp

+#include "precomp.hpp"
+#include <climits>
+#include <iostream>
+using namespace std;
+
+#ifdef _MSC_VER
+#   pragma warning(disable: 4512)
+#endif
+
+namespace cv
+{
+
+typedef Vec<float, 1> Vec1f;
+typedef Vec<uchar, 1> Vec1b;
+
+#ifndef SQR
+#define SQR(a) ((a)*(a))
+#endif
+
+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;
+
+                register JointVec jointPix0 = *jointCenterPixPtr;
+                register SrcVec srcSum = SrcVec::all(0.0f);
+                register 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;
+
+                register JointVeci jointPix0 = JointVeci(*jointCenterPixPtr);
+                register SrcVecf srcSum = SrcVecf::all(0.0f);
+                register 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() && 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);
+
+    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();
+
+    if (joint.depth() == CV_8U)
+    {
+        jointBilateralFilter_8u(joint, src, dst, radius, sigmaColor, sigmaSpace, borderType);
+    }
+    else
+    {
+        jointBilateralFilter_32f(joint, src, dst, radius, sigmaColor, sigmaSpace, borderType);   
+    }
+    
+}
+
+}
\ No newline at end of file

modules/ximgproc/test/test_adaptive_manifold_ref_impl.cpp

@@ -883,7 +883,7 @@ namespace
             }
             else
             {
-                for (int i = 0; i < buf_.rand_vec.total(); i++)
+                for (int i = 0; i < (int)buf_.rand_vec.total(); i++)
                     buf_.rand_vec(0, i) = (i % 2 == 0) ? 0.5f : -0.5f;
             }
 
@@ -914,10 +914,10 @@ namespace
             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, sigma_s_, df, buf_);
+            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, sigma_s_, df, buf_);
+            calcEta(src_joint_f_, buf_.theta, cluster_plus, eta_plus, (float)sigma_s_, (float)df, buf_);
 
             // Algorithm 1, Step 5: recursively build more manifolds.
 

modules/ximgproc/test/test_joint_bilateral_filter.cpp

+#include "test_precomp.hpp"
+
+using namespace std;
+using namespace std::tr1;
+using namespace cv;
+using namespace testing;
+
+#ifdef _MSC_VER
+#   pragma warning(disable: 4512)
+#endif
+
+namespace cvtest
+{
+
+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;
+}
+
+//////////////////////////////////////////////////////////////////////////
+//////////////////////////////////////////////////////////////////////////
+
+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 = BORDER_DEFAULT)
+{
+    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)
+    )
+);
+
+}
\ No newline at end of file