解析opencv的 Box Filter的实现并提出进一步加速的方案。

说明：本文所有算法的涉及到的优化均指在PC上进行的，对于其他构架是否合适未知，请自行试验。

Box Filter，最经典的一种领域操作，在无数的场合中都有着广泛的应用，作为一个很基础的函数，其性能的好坏也直接影响着其他相关函数的性能，最典型莫如现在很好的EPF滤波器：GuideFilter。因此其优化的档次和程度是非常重要的，网络上有很多相关的代码和博客对该算法进行讲解和优化，提出了不少O(1)算法，但所谓的0(1)算法也有优劣之分，0(1)只是表示执行时间和某个参数无关，但本身的耗时还是有区别。比较流行的就是累积直方图法，其实这个是非常不可取的，因为稍大的图就会导致溢出，这里我们解析下 opencv的相关代码，看看他是如何进行优化的。

首先找到opencv的代码的位置，其在\opencv\sources\modules\imgproc\src\smooth.cpp中。

Box Filter 是一种行列可分离的滤波，因此，opencv也是这样处理的，先进行行方向的滤波，得到中间结果，然后再对中间结果进行列方向的处理，得到最终的结果。

opencv 行方向处理的相关代码如下：

template<typename T, typename ST>
struct RowSum :
        public BaseRowFilter
{
    RowSum( int _ksize, int _anchor ) :
        BaseRowFilter()
    {
        ksize = _ksize;
        anchor = _anchor;
    }

    virtual void operator()(const uchar* src, uchar* dst, int width, int cn)
    {
        const T* S = (const T*)src;
        ST* D = (ST*)dst;
        int i = 0, k, ksz_cn = ksize*cn;

        width = (width - 1)*cn;
        for( k = 0; k < cn; k++, S++, D++ )
        {
            ST s = 0;
            for( i = 0; i < ksz_cn; i += cn )
                s += S[i];
            D[0] = s;
            for( i = 0; i < width; i += cn )
            {
                s += S[i + ksz_cn] - S[i];
                D[i+cn] = s;
            }
        }
    }
};

　　这个代码考虑了多个通道以及多种数据类型的情况，为了分析方便我们重写下单通道时的核心代码。

for(Z = 0, Value = RowData[0]; Z < Size; Z++)
    Value += RowData[Z];
LinePD[0] = Value;

for(X = 1; X < Width; X ++)
{
    Value += RowData[X + Size - 1] - RowData[X - 1];
    LinePD[X] = Value;
}

　　上述代码中RowData是指对某一行像素进行扩展后的像素值，其宽度为Width + Radius + Radius, Radius是值滤波的半径， 而代码中Size = 2 * Radius + 1，LinePD即所谓的中间结果，考虑数据类型能表达的范围，必须使用 int类型。

对于每行第一个点很简单，直接用for计算出行方向的指定半径内的累加值。而之后所有点，用前一个累计值加上新加入的点，然后去除掉移出的哪一个点，得到新的累加值。这个方法在很多文献中都有提及，并没有什么新鲜之处，这里不多说。

按照正常的思维，在列方向的处理应该和行方向完全相同，只是处理的方向的不同、处理的数据源的不同以及最后需要对计算结果多一个归一化的过程而已。事实也是如此，有很多算法都是这样做的，但是由于CPU构架的一些原因（主要是cache miss的增加），同样的算法沿列方向处理总是会比沿行方向慢一个档次，解决方式有很多，例如先对中间结果进行转置，然后再按照行方向的规则进行处理，处理完后在将数据转置回去。转置过程是有非常高效的处理方式的，借助于SSE以及Cache优化，能实现比原始两重for循环快4倍以上的效果。还有一种方式就正如opencv中列处理过程所示，这正是下面将要重点描述的。

在opencv的代码中，并没有直接沿列方向处理，而是继续沿着行的方向一行一行处理，我先贴出我自己翻译的有关纯C的代码进行解说：

    for (Y = 0; Y < Size - 1; Y++)            //    注意没有最后一项哦
    {
        int *LinePS = (int *)(Sum->Data + ColPos[Y] * Sum->WidthStep);
        for(X = 0; X < Width * Channel; X++)    ColSum[X] += LinePS[X];
    }

    for (Y = 0; Y < Height; Y++)
    {
        unsigned char* LinePD    = Dest->Data + Y * Dest->WidthStep;
        int *AddPos                = (int*)(Sum->Data + ColPos[Y + Size - 1] * Sum->WidthStep);
        int *SubPos                = (int*)(Sum->Data + ColPos[Y] * Sum->WidthStep);

        for(X = 0; X < Width * Channel; X++)
        {
            Value = ColSum[X] + AddPos[X];
            LinePD[X] = Value * Scale;                    //            / Amount 这样写在PC上会更慢
            ColSum[X] = Value - SubPos[X];
        }
    }

接着我们看opencv的列方向上的代码，这里只贴出他针对uchar类型的图像数据的处理。

template<>
struct ColumnSum<int, uchar> :
        public BaseColumnFilter
{
    ColumnSum( int _ksize, int _anchor, double _scale ) :
        BaseColumnFilter()
    {
        ksize = _ksize;
        anchor = _anchor;
        scale = _scale;
        sumCount = 0;
    }

    virtual void reset() { sumCount = 0; }

    virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
    {
        int i;
        int* SUM;
        bool haveScale = scale != 1;
        double _scale = scale;

        #if CV_SSE2
            bool haveSSE2 = checkHardwareSupport(CV_CPU_SSE2);
        #endif

        if( width != (int)sum.size() )
        {
            sum.resize(width);
            sumCount = 0;
        }

        SUM = &sum[0];
        if( sumCount == 0 )
        {
            memset((void*)SUM, 0, width*sizeof(int));
            for( ; sumCount < ksize - 1; sumCount++, src++ )
            {
                const int* Sp = (const int*)src[0];
                i = 0;
                #if CV_SSE2
                if(haveSSE2)
                {
                    for( ; i <= width-4; i+=4 )
                    {
                        __m128i _sum = _mm_loadu_si128((const __m128i*)(SUM+i));
                        __m128i _sp = _mm_loadu_si128((const __m128i*)(Sp+i));
                        _mm_storeu_si128((__m128i*)(SUM+i),_mm_add_epi32(_sum, _sp));
                    }
                }
                #elif CV_NEON
                for( ; i <= width - 4; i+=4 )
                    vst1q_s32(SUM + i, vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i)));
                #endif
                for( ; i < width; i++ )
                    SUM[i] += Sp[i];
            }
        }
        else
        {
            CV_Assert( sumCount == ksize-1 );
            src += ksize-1;
        }

        for( ; count--; src++ )
        {
            const int* Sp = (const int*)src[0];
            const int* Sm = (const int*)src[1-ksize];
            uchar* D = (uchar*)dst;
            if( haveScale )
            {
                i = 0;
                #if CV_SSE2
                if(haveSSE2)
                {
                    const __m128 scale4 = _mm_set1_ps((float)_scale);
                    for( ; i <= width-8; i+=8 )
                    {
                        __m128i _sm  = _mm_loadu_si128((const __m128i*)(Sm+i));
                        __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));

                        __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
                                                     _mm_loadu_si128((const __m128i*)(Sp+i)));
                        __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
                                                      _mm_loadu_si128((const __m128i*)(Sp+i+4)));

                        __m128i _s0T = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s0)));
                        __m128i _s0T1 = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s01)));

                        _s0T = _mm_packs_epi32(_s0T, _s0T1);

                        _mm_storel_epi64((__m128i*)(D+i), _mm_packus_epi16(_s0T, _s0T));

                        _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
                        _mm_storeu_si128((__m128i*)(SUM+i+4),_mm_sub_epi32(_s01,_sm1));
                    }
                }
                #elif CV_NEON
                float32x4_t v_scale = vdupq_n_f32((float)_scale);
                for( ; i <= width-8; i+=8 )
                {
                    int32x4_t v_s0 = vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i));
                    int32x4_t v_s01 = vaddq_s32(vld1q_s32(SUM + i + 4), vld1q_s32(Sp + i + 4));

                    uint32x4_t v_s0d = cv_vrndq_u32_f32(vmulq_f32(vcvtq_f32_s32(v_s0), v_scale));
                    uint32x4_t v_s01d = cv_vrndq_u32_f32(vmulq_f32(vcvtq_f32_s32(v_s01), v_scale));

                    uint16x8_t v_dst = vcombine_u16(vqmovn_u32(v_s0d), vqmovn_u32(v_s01d));
                    vst1_u8(D + i, vqmovn_u16(v_dst));

                    vst1q_s32(SUM + i, vsubq_s32(v_s0, vld1q_s32(Sm + i)));
                    vst1q_s32(SUM + i + 4, vsubq_s32(v_s01, vld1q_s32(Sm + i + 4)));
                }
                #endif
                for( ; i < width; i++ )
                {
                    int s0 = SUM[i] + Sp[i];
                    D[i] = saturate_cast<uchar>(s0*_scale);
                    SUM[i] = s0 - Sm[i];
                }
            }
            else
            {
                i = 0;
                #if CV_SSE2
                if(haveSSE2)
                {
                    for( ; i <= width-8; i+=8 )
                    {
                        __m128i _sm  = _mm_loadu_si128((const __m128i*)(Sm+i));
                        __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));

                        __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
                                                     _mm_loadu_si128((const __m128i*)(Sp+i)));
                        __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
                                                      _mm_loadu_si128((const __m128i*)(Sp+i+4)));

                        __m128i _s0T = _mm_packs_epi32(_s0, _s01);

                        _mm_storel_epi64((__m128i*)(D+i), _mm_packus_epi16(_s0T, _s0T));

                        _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
                        _mm_storeu_si128((__m128i*)(SUM+i+4),_mm_sub_epi32(_s01,_sm1));
                    }
                }
                #elif CV_NEON
                for( ; i <= width-8; i+=8 )
                {
                    int32x4_t v_s0 = vaddq_s32(vld1q_s32(SUM + i), vld1q_s32(Sp + i));
                    int32x4_t v_s01 = vaddq_s32(vld1q_s32(SUM + i + 4), vld1q_s32(Sp + i + 4));

                    uint16x8_t v_dst = vcombine_u16(vqmovun_s32(v_s0), vqmovun_s32(v_s01));
                    vst1_u8(D + i, vqmovn_u16(v_dst));

                    vst1q_s32(SUM + i, vsubq_s32(v_s0, vld1q_s32(Sm + i)));
                    vst1q_s32(SUM + i + 4, vsubq_s32(v_s01, vld1q_s32(Sm + i + 4)));
                }
                #endif

                for( ; i < width; i++ )
                {
                    int s0 = SUM[i] + Sp[i];
                    D[i] = saturate_cast<uchar>(s0);
                    SUM[i] = s0 - Sm[i];
                }
            }
            dst += dststep;
        }
    }

    double scale;
    int sumCount;
    std::vector<int> sum;
};

　　代码比较多，我稍微精简下（注意，精简后的是不可运行的，只是更清晰的表达了思路）。

virtual void operator()(const uchar** src, uchar* dst, int dststep, int count, int width)
{
    int i;
    int* SUM;
    bool haveScale = scale != 1;
    double _scale = scale;
    if( width != (int)sum.size() )
    {
        sum.resize(width);
        sumCount = 0;
    }

    SUM = &sum[0];
    if( sumCount == 0 )
    {
        memset((void*)SUM, 0, width*sizeof(int));
        for( ; sumCount < ksize - 1; sumCount++, src++ )
        {
            const int* Sp = (const int*)src[0];
            i = 0;

            for( ; i <= width-4; i+=4 )
            {
                __m128i _sum = _mm_loadu_si128((const __m128i*)(SUM+i));
                __m128i _sp = _mm_loadu_si128((const __m128i*)(Sp+i));
                _mm_storeu_si128((__m128i*)(SUM+i),_mm_add_epi32(_sum, _sp));
            }
            for( ; i < width; i++ )
                SUM[i] += Sp[i];
        }
    }
    else
    {
        src += ksize-1;
    }

    for( ; count--; src++ )
    {
        const int* Sp = (const int*)src[0];
        const int* Sm = (const int*)src[1-ksize];
        uchar* D = (uchar*)dst;

        i = 0;
        const __m128 scale4 = _mm_set1_ps((float)_scale);
        for( ; i <= width-8; i+=8 )
        {
            __m128i _sm  = _mm_loadu_si128((const __m128i*)(Sm+i));
            __m128i _sm1  = _mm_loadu_si128((const __m128i*)(Sm+i+4));

            __m128i _s0  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i)),
                                         _mm_loadu_si128((const __m128i*)(Sp+i)));
            __m128i _s01  = _mm_add_epi32(_mm_loadu_si128((const __m128i*)(SUM+i+4)),
                                          _mm_loadu_si128((const __m128i*)(Sp+i+4)));

            __m128i _s0T = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s0)));
            __m128i _s0T1 = _mm_cvtps_epi32(_mm_mul_ps(scale4, _mm_cvtepi32_ps(_s01)));

            _s0T = _mm_packs_epi32(_s0T, _s0T1);

            _mm_storel_epi64((__m128i*)(D+i), _mm_packus_epi16(_s0T, _s0T));

            _mm_storeu_si128((__m128i*)(SUM+i), _mm_sub_epi32(_s0,_sm));
            _mm_storeu_si128((__m128i*)(SUM+i+4),_mm_sub_epi32(_s01,_sm1));
        }
        for( ; i < width; i++ )
        {
            int s0 = SUM[i] + Sp[i];
            D[i] = saturate_cast<uchar>(s0*_scale);
            SUM[i] = s0 - Sm[i];
        }

        dst += dststep;
    }
}

opencv的代码在这里就没有继续优化了，本人认为虽然这里的操作是前后依赖的，全局无法并行化，但是观察这一行代码：

Value += RowData[X + Size - 1] - RowData[X - 1];

其中RowData[X + Size - 1] - RowData[X - 1]; 并不是前后依赖的，是可以并行的，是可以借助于因此如果提前计算出所有的差值，

				


				
时间： 2024-10-05 23:58:09 

		


		


	

	
	

		


		
解析opencv的 Box Filter的实现并提出进一步加速的方案。的相关文章


						

		

			

                解析opencv中Box Filter的实现并提出进一步加速的方案（源码共享）。
			

		

		

									

				说明:本文所有算法的涉及到的优化均指在PC上进行的,对于其他构架是否合适未知,请自行试验. Box Filter,最经典的一种领域操作,在无数的场合中都有着广泛的应用,作为一个很基础的函数,其性能的好坏也直接影响着其他相关函数的性能,最典型莫如现在很好的EPF滤波器:GuideFilter.因此其优化的档次和程度是非常重要的,网络上有很多相关的代码和博客对该算法进行讲解和优化,提出了不少O(1)算法,但所谓的0(1)算法也有优劣之分,0(1)只是表示执行时间和某个参数无关,但本身的耗时还是有区别			

		

	

						

		

			

                用 OpenCV 实现 Guided Filter
			

		

		

									

				最近发现 OpenCV 的 C++ 接口写得和 Matlab 实在是太相似了.我用 OpenCV 2.49 来实现 Guided Filter 的时候几乎是直接照抄作者公开的 Matlab 代码.废话不多直接贴我的代码: cv::Mat guidedFilter(cv::Mat I, cv::Mat p, int r, double eps) { /* % GUIDEDFILTER O(1) time implementation of guided filter. % % - guidance			

		

	

				

		

			

                FilterEngine 类解析——OpenCV图像滤波核心引擎(zz)
			

		

		

									

				<2>FilterEngine 类解析——OpenCV图像滤波核心引擎 FilterEngine类是OpenCV关于图像滤波的主力军类,OpenCV图像滤波功能的核心引擎.各种滤波函数比如blur, GaussianBlur,到头来其实是就是在函数末尾处定义了一个Ptr<FilterEngine>类型的f,然后f->apply( src, dst )了一下而已. 这个类可以把几乎是所有的滤波操作施加到图像上.它包含了所有必要的中间缓存器.有很多和滤波相关的create系函数的			

		

	

				

		

			

                10. Dubbo原理解析-代理之Listener &amp; filter
			

		

		

									

				Listener ExporterListener: dubbo在服务暴露(exporter)以及销毁暴露(unexporter)服务的过程中提供了回调窗口,供用户做业务处理.ProtocolListenerWrapper在暴露过程中构建了监听器链 public class ProtocolListenerWrapper implements Protocol { public <T> Exporter<T> export(Invoker<T>invoker) thro			

		

	

				

		

			

                盒式滤波器Box Filter
			

		

		

									

				其主要功能是:在给定的滑动窗口大小下,对每个窗口内的像素值进行快速相加求和 在模式识别领域,Haar特征是大家非常熟悉的一种图像特征了,它可以应用于许多目标检测的算法中.与Haar相似,图像的局部矩形内像素的和.平方和.均值.方差等特征也可以用类似Haar特征的计算方法来计算.这些特征有时会频繁的在某些算法中使用,因此对它的优化势在必行.Boxfilter就是这样一种优化方法,它可以使复杂度为O(MN)的求和,求方差等运算降低到O(1)或近似于O(1)的复杂度,它的缺点是不支持多尺度. 第一个提			

		

	

				

		

			

                【Dubbo 源码解析】03_Dubbo Protocol&amp;Filter
			

		

		

									

				Protocol & Filter Dubbo 服务暴露和服务引用都是通过的 com.alibaba.dubbo.rpc.Protocol 来实现的.它是一个 SPI 扩展. @SPI("dubbo") public interface Protocol { int getDefaultPort(); @Adaptive <T> Exporter<T> export(Invoker<T> invoker) throws RpcExceptio			

		

	

				

		

			

                超高速导向滤波实现过程纪要（欢迎挑战）
			

		

		

									

				自从何凯明提出导向滤波后,因为其算法的简单性和有效性,该算法得到了广泛的应用,以至于新版的matlab都将其作为标准自带的函数之一了,利用他可以解决的所有的保边滤波器的能解决的问题,比如细节增强.HDR压缩.细节羽化.去雾.风格化,而且由于其保边特性,如果很多传统函数中使用高斯滤波或者均值滤波的地方用他代替,能很好解决一些强边缘的过渡不自然问题,比如retinex.Highlight/shadow等应用中,因此,快速的实现该算法具有很强的适用意义. 本文简要的记录了本人在优化导向滤波实现的过程中			

		

	

				

		

			

                模板匹配
			

		

		

									

				欧式距离 || 图像处理中任意核卷积(matlab中conv2函数)的快速实现 || 解析opencv中Box Filter的实现并提出进一步加速的方案(源码共享) 标准的基于欧式距离的模板匹配优化源码和实现 解析opencv中Box Filter的实现并提出进一步加速的方案(源码共享) 原文地址:https://www.cnblogs.com/LieYanAnYing/p/12404930.html			

		

	

				

		

			

                CVPR论文《100+ Times Faster Weighted Median Filter (WMF)》的实现和解析（附源代码）。
			

		

		

									

				四年前第一次看到<100+ Times FasterWeighted Median Filter (WMF)>一文时,因为他附带了源代码,而且还是CVPR论文,因此,当时也对代码进行了一定的整理和解读,但是当时觉得这个算法虽然对原始速度有不少的提高,但是还是比较慢.因此,没有怎么在意,这几天有几位朋友又提到这篇文章,于是把当时的代码和论文又仔细的研读了一番,对论文的思想和其中的实现也有了一些新的新的,再次做个总结和分享. 这篇文章的官网地址是:http://www.cse.cuhk.edu.h