并行加速实战 双边滤波器

最新推荐文章于 2024-07-30 09:20:48 发布
最新推荐文章于 2024-07-30 09:20:48 发布
阅读量2.3k 收藏 4
点赞数
分类专栏： C/C++编程 SIMD
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/billbliss/article/details/78438097
版权

之前分析了 二维中值滤波器的并行加速

由于二维中值滤波器是控制密集型的滤波器（排序操作），所以SSE加速不太明显


这次选用了计算密集型的双边滤波器

针对双边滤波器在5*5的滤波核下的运算速度做优化和分析

以下会有主区域、全图、主循环、完整（初始化+主循环）的概念


1.     由于双边滤波的滤波半径为2+1，所以不能忽略图像四周边界的区域了。

所以，以下会对主区域和全图滤波做一个预算时间的对比。

 

2.     在快速算法中还做了查找表优化，所以滤波函数是个有状态滤波器，算法需要初始化。

以下会对主循环和包含初始化的运算时间做个对比


总结先写：

发现同样的版本 使用浮点比整形还要快。浮点的SSE并行并没有提速4倍。

目前分析是因为编译器自动做了优化，将浮点运算的速度提升了2-3倍。以至于比整形的版本还略快一点。

 

运算时间

旧版整型 主区域 主循环

8.766ms

新版整型 全图   主循环

6.324ms

Opencv  浮点 全图 主循环

5.713ms

新版    浮点 全图 主循环

5.301ms

SSE 浮点 全图 主循环

4.778ms

omp浮点 全图 主循环

1.527ms

SSE+omp浮点 全图 主循环

1.355ms


并行算法优化分析

1.     整型双边滤波

 

运算时间

旧版整型 主区域 主循环

8.766ms

新版整型 全图   主循环

6.324ms

 

2.     整型、浮点型双边滤波

 

运算时间

新版  整型 全图 主循环

6.324ms

新版  浮点 全图 主循环

5.301ms

opencv浮点 全图 主循环

5.713ms

 

3.     浮点型SSE双边滤波主区域、全图耗时

 

 

主循环

完整

浮点型

主区域

5.002ms

5.078ms

全图

5.198ms

5.301ms

浮点型SSE

主区域

4.658ms

4.764ms

全图

4.778ms

5.076ms

 

4.     浮点型SSE omp双边滤波

 

 

主循环

浮点型

主区域

5.002ms

全图

5.198ms

浮点型SSE

主区域

4.658ms

全图

4.778ms

浮点型omp

主区域

1.434ms

全图

1.527ms

浮点型SSE+omp

主区域

1.285ms

全图

1.355ms

 

 

 

一．重构了关于整型优化的双边滤波器

原先版本是对矩形区域做滤波的，现在改成了圆形区域。减少了近一半的计算量。

旧版的整型双边滤波主区域主循环耗时为8.766ms

新版的整型双边滤波全图主循环耗时为6.324ms

 

二．设计了浮点型优化的双边滤波器

浮点型双边滤波主区域主循环耗时为5.002ms

浮点型双边滤波主区域完整耗时为5.078ms

浮点型双边滤波全图主循环耗时为5.198ms

浮点型双边滤波全图完整耗时为5.301ms

Opencv浮点型双边滤波全图完整耗时为5.713ms

 

这里可以看出全图运算大概比主区域多耗时0.2ms

算法初始化耗时0.1ms

 

三．增加了SSE加速优化的双边滤波器

浮点型SSE加速双边滤波主区域主循环耗时为4.658ms

浮点型SSE加速双边滤波主区域完整耗时为4.764ms

浮点型SSE加速双边滤波全图主循环耗时为4.778ms

浮点型SSE加速双边滤波全图完整耗时为5.076ms

 

这里可以看出全图运算大概比主区域多耗时0.1ms

算法初始化耗时0.2ms

 

四．增加了omp加速优化的双边滤波

浮点型omp加速双边滤波主区域主循环耗时为1.434ms

浮点型omp加速双边滤波全图主循环耗时为1.527ms

 

这里可以看出全图运算大概比主区域多耗时0.1ms

 

五．增加了SSE+omp加速优化的双边滤波

浮点型SSE+omp加速双边滤波主区域主循环耗时为1.285ms

浮点型SSE+omp加速双边滤波全图主循环耗时为1.355ms

 

这里可以看出全图运算大概比主区域多耗时0.1ms

 

以下是具体的运算优化耗时

l  opencv full 5.713ms

l  mainBody mainLoop 5.002ms

l  mainBody 5.078ms

l  Full mainLoop5.198ms

l  Full 5.301ms

l  mainBody sse mainLoop 4.658ms

l  mainBody sse 4.764ms

l  Full sse mainLoop 4.778ms

l  Full sse 5.076ms

l  mainBody omp mainLoop 1.434ms

l  Full omp mainLoop 1.527ms

l  mainBody sse_omp mainLoop 1.285ms

l  Full sse_omp mainLoop 1.355ms

l  Int Full mainLoop 6.324ms

l  INT old version mainBody mainLoop 8.766ms

 

 具体代码如下

宏定义

因为整型的双边会有求和溢出风险，所以这里限制了滤波直径为11/半径5

#define MALLOC               malloc
#define FREE(p)              if (p != NULL) { free(p); p = NULL;}

#define ALIGN16              __declspec(align(16))
#define ALIGN_MALLOC16(n)    _aligned_malloc(n, 16)
#define ALIGN_MALLOC32(n)    _aligned_malloc(n, 32)
#define ALIGN_MALLOC64(n)    _aligned_malloc(n, 64)
#define ALIGN_MALLOC128(n)   _aligned_malloc(n, 128)
#define ALIGN_FREE(p)        if (p != NULL) { _aligned_free(p); p = NULL;}

#define BF_INT_BITS      (10)
#define BF_INT_SCALE     (1 << BF_INT_BITS)
#define BF_INT_SHIFT     ((S32)(((BF_INT_BITS + 1) << 1) - 31 + 16) + 7) //more 7 bits[>11*11]
#define BF_INT_BITS2     ((S32)((BF_INT_BITS << 1) - BF_INT_SHIFT))
#define BF_INT_SCALE2    (1 << BF_INT_BITS2)

#define BF_BUF_LEN	 (1024)
#define EDGEPRES_R_MAX   (5)
类的定义 
typedef class edgePresFiltMain
{
public:

    edgePresFiltMain(S32 width, S32 height, S32 radius, F32 sigmaVal, F32 sigmaSp);
    ~edgePresFiltMain();

    void edgePreserve_mainBody(U16 *src, U16 *dst);

    void edgePreserve_mainBody_omp(U16 *src, U16 *dst);

    void edgePreserve_mainBody_sse(U16 *src, U16 *dst);

    void edgePreserve_mainBody_sse_omp(U16 *src, U16 *dst);

    void *hdl;

}edgePresFiltMain_;

typedef class edgePresFilt
{
public:

    edgePresFilt(S32 width, S32 height, S32 radius, F32 sigmaVal, F32 sigmaSp);
    ~edgePresFilt();

    void edgePreserve(U16 *src, U16 *dst);

    void edgePreserve_omp(U16 *src, U16 *dst);

    void edgePreserve_sse(U16 *src, U16 *dst);

    void edgePreserve_sse_omp(U16 *src, U16 *dst);

    void *hdl;

}edgePresFilt_;

typedef class edgePresFiltInt
{
public:

    edgePresFiltInt(S32 width, S32 height, S32 radius, F32 sigmaVal, F32 sigmaSp);
    ~edgePresFiltInt();

    void edgePreserve(U16 *src, U16 *dst);

    void edgePreserve_omp(U16 *src, U16 *dst);

    void *hdl;

}edgePresFiltInt_;
纯C函数声明 
// no smooth on the border
void edgePreserve_mainBody(U16 *src, U16 *dst, S32 width, S32 height,
    S32 radius, F32 sigmaVal, F32 sigmaSp);

// smooth on the border
void edgePreserve(U16 *src, U16 *dst, S32 width, S32 height,
    S32 radius, F32 sigmaVal, F32 sigmaSp);

// no smooth on the border
void edgePreserve_mainBody_sse(U16 *src, U16 *dst, S32 width, S32 height,
    S32 radius, F32 sigmaVal, F32 sigmaSp);

// smooth on the border
void edgePreserve_sse(U16 *src, U16 *dst, S32 width, S32 height,
    S32 radius, F32 sigmaVal, F32 sigmaSp);

// smooth on the border
void edgePreserveInt(U16 *src, U16 *dst, S32 width, S32 height,
    S32 radius, F32 sigmaVal, F32 sigmaSp);
边界处理 
static void borderReflect(U16 *src, S32 width, S32 height, U16 *dst, S32 radius)
{
    S32 i = 0;
    S32 j = 0;
    S32 itmp1 = radius - 1;
    S32 itmp2 = -1 - radius;
    S32 width2 = width + (radius << 1);
    U16 *psrc = src;
    U16 *pdst = dst + width2 * radius;

    for (i = 0; i < height; i++)
    {
        for (j = 0; j < radius; j++)
        {
            pdst[j] = psrc[itmp1 - j];
        }

        memcpy(pdst + radius, psrc, sizeof(U16) * width);

        psrc += width;
        pdst += width2;

        for (j = -1; j >= -radius; j--)
        {
            pdst[j] = psrc[itmp2 - j];
        }
    }

    psrc = pdst - width2;

    for (i = 0; i < radius; i++)
    {
        memcpy(pdst, psrc, sizeof(U16) * width2);
        psrc -= width2;
        pdst += width2;
    }

    psrc = dst + width2 * radius;
    pdst = dst + width2 * (radius - 1);

    for (i = 0; i < radius; i++)
    {
        memcpy(pdst, psrc, sizeof(U16) * width2);
        psrc += width2;
        pdst -= width2;
    }
}

static void borderCopy(U16 *src, U16 *dst, S32 width, S32 height, S32 radius)
{
    S32 i = 0;
    S32 j = 0;
    S32 xend = height - (radius << 1);
    U16 *psrc = src;
    U16 *pdst = dst;

    memcpy(pdst, psrc, sizeof(U16) * (width * radius - radius));

    psrc += radius * width;
    pdst += radius * width;

    for (i = 0; i < xend; i++)
    {
        memcpy(pdst - radius, psrc - radius, sizeof(U16) * (radius << 1));
        psrc += width;
        pdst += width;
    }

    memcpy(pdst - radius, psrc - radius, sizeof(U16) * (width * radius + radius));
}


static void edgePreserve_LUT(S32 radius, S32 width, F32 sigmaVal, F32 sigmaSp,
    S32 *spOfs, F32 *spWt, F32 *valWt)
{
    S32 i = 0;
    S32 j = 0;
    S32 k = 0;
    F32 sigmaValCoeff = -0.5f / (sigmaVal * sigmaVal);
    F32 sigmaSpCoeff = -0.5f / (sigmaSp * sigmaSp);

    for (i = -radius; i <= radius; i++)
    {
        for (j = -radius; j <= radius; j++)
        {
            if (sqrtf((i * i + j * j) * 1.f) <= radius)
            {
                if ((i == 0) && (j == 0))
                {
                    continue;
                }

                spOfs[k] = i * width + j;
                spWt[k] = expf((i * i + j * j) * sigmaSpCoeff);
                k++;
            }
        }
    }

    for (i = 0; i < BF_BUF_LEN - 1; i++)
    {
        valWt[i] = expf(i * i * sigmaValCoeff);
    }
    valWt[i] = 0.f;
}


// Smooth main body of src to the dst img
// src img is bigger than dst img
static void edgePreserve_mainBody_process(U16 *src, S32 srcWidth, S32 srcHeight, S32 srcStep,
    U16 *dst, S32 dstStep, S32 radius, S32 *spOfs, F32 *spWt, F32 *valWt, S32 maxk)
{
    const S32 xEnd = srcWidth - radius;
    const S32 dstHeight = srcHeight - (radius << 1);

    S32 i = 0;
    S32 j = 0;
    S32 k = 0;
    S32 n = 0;
    U16 val0 = 0;
    U16 val = 0;
    S32 tmp = 0;
    F32 w = 0.f;
    F32 sum = 0.f;
    F32 wsum = 0.f;

    U16 *psrc = src + radius * srcStep;
    U16 *pdst = dst;

    for (i = 0; i < dstHeight; i++)
    {
        for (j = radius; j < xEnd; j++)
        {
            val0 = psrc[j];
            sum = 0.f;
            wsum = 0.f;

            for (k = 0; k < maxk; k++)
            {
                val = psrc[j + spOfs[k]];
                tmp = val - val0;
                tmp = ABS(tmp);
                tmp = MIN(tmp, BF_BUF_LEN - 1);
                w = spWt[k] * valWt[tmp];
                sum += val * w;
                wsum += w;
            }

            pdst[j - radius] = (U16)((sum + val0) / (wsum + 1.f));
        }

        psrc += srcStep;
        pdst += dstStep;
    }
}

// Smooth main body of src to the dst img
// src img is bigger than dst img
static void edgePreserve_mainBody_omp_process(U16 *src, S32 srcWidth, S32 srcHeight, S32 srcStep,
    U16 *dst, S32 dstStep, S32 radius, S32 *spOfs, F32 *spWt, F32 *valWt, S32 maxk)
{
    const S32 xEnd = srcWidth - radius;
    const S32 dstHeight = srcHeight - (radius << 1);

    S32 i = 0;

#pragma omp parallel for
    for (i = 0; i < dstHeight; i++)
    {
        S32 j = 0;
        S32 k = 0;
        U16 val0 = 0;
        U16 val = 0;
        S32 tmp = 0;
        F32 w = 0.f;
        F32 sum = 0.f;
        F32 wsum = 0.f;

        U16 *psrc = src + (radius + i) * srcStep;
        U16 *pdst = dst + i * dstStep;

        for (j = radius; j < xEnd; j++)
        {
            val0 = psrc[j];
            sum = 0.f;
            wsum = 0.f;

            for (k = 0; k < maxk; k++)
            {
                val = psrc[j + spOfs[k]];
                tmp = val - val0;
                tmp = ABS(tmp);
                tmp = MIN(tmp, BF_BUF_LEN - 1);
                w = spWt[k] * valWt[tmp];
                sum += val * w;
                wsum += w;
            }

            pdst[j - radius] = (U16)((sum + val0) / (wsum + 1.f));
        }
    }
}

// Smooth main body of src to the dst img
// src img is bigger than dst img
static void edgePreserve_mainBody_sse_process(U16 *src, S32 srcWidth, S32 srcHeight, S32 srcStep,
    U16 *dst, S32 dstStep, S32 radius, S32 *spOfs, F32 *spWt, F32 *valWt, S32 maxk)
{
    const S32 xEnd = srcWidth - radius;
    const S32 dstHeight = srcHeight - (radius << 1);
    const U32 ALIGN16 bufSignMask[] = { 0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff };
    const F32 ALIGN16 bufLutLen[] = { BF_BUF_LEN - 1, BF_BUF_LEN - 1, BF_BUF_LEN - 1, BF_BUF_LEN - 1 };

    S32 ALIGN16 buf[4];
    F32 ALIGN16 bufSum[4];
    S32 i = 0;
    S32 j = 0;
    S32 k = 0;
    F32 val0 = 0.f;
    F32 sum = 0.f;
    F32 wsum = 0.f;

    U16 *psrc = src + radius * srcStep;
    U16 *pdst = dst;

    __m128 _val;
    __m128 _val0;
    __m128 _idx;
    __m128 _psum;
    __m128 _sw;
    __m128 _cw;
    __m128 _w;

    const __m128 _signMask = _mm_load_ps((const float*)bufSignMask);
    const __m128 _lutLen = _mm_load_ps((const float*)bufLutLen);

    for (i = 0; i < dstHeight; i++)
    {
        for (j = radius; j < xEnd; j++)
        {
            val0 = psrc[j] * 1.0f;
            _psum = _mm_setzero_ps();
            _val0 = _mm_set1_ps(val0);

            for (k = 0; k <= maxk - 4; k += 4)
            {
                _val = _mm_set_ps(psrc[j + spOfs[k + 3]], psrc[j + spOfs[k + 2]],
                    psrc[j + spOfs[k + 1]], psrc[j + spOfs[k]]);
                //                 _sw = _mm_loadu_ps(spWt + k);
                _sw = _mm_load_ps(spWt + k);
                _idx = _mm_and_ps(_signMask, _mm_sub_ps(_val, _val0));
                _mm_store_si128((__m128i*)buf, _mm_cvtps_epi32(_mm_min_ps(_idx, _lutLen)));
                _cw = _mm_set_ps(valWt[buf[3]], valWt[buf[2]], valWt[buf[1]], valWt[buf[0]]);
                _w = _mm_mul_ps(_cw, _sw);
                _val = _mm_mul_ps(_w, _val);
                _sw = _mm_hadd_ps(_w, _val);
                _sw = _mm_hadd_ps(_sw, _sw);
                _psum = _mm_add_ps(_sw, _psum);
            }

            _mm_storel_pi((__m64*)bufSum, _psum);
            sum = bufSum[1] + val0;
            wsum = bufSum[0] + 1.f;

            pdst[j - radius] = (U16)(sum / wsum);
        }

        psrc += srcStep;
        pdst += dstStep;
    }
}

// Smooth main body of src to the dst img
// src img is bigger than dst img
static void edgePreserve_mainBody_sse_omp_process(U16 *src, S32 srcWidth, S32 srcHeight, S32 srcStep,
    U16 *dst, S32 dstStep, S32 radius, S32 *spOfs, F32 *spWt, F32 *valWt, S32 maxk)
{
    const S32 xEnd = srcWidth - radius;
    const S32 dstHeight = srcHeight - (radius << 1);
    const U32 ALIGN16 bufSignMask[] = { 0x7fffffff, 0x7fffffff, 0x7fffffff, 0x7fffffff };
    const F32 ALIGN16 bufLutLen[] = { BF_BUF_LEN - 1, BF_BUF_LEN - 1, BF_BUF_LEN - 1, BF_BUF_LEN - 1 };
    const __m128 _signMask = _mm_load_ps((const float*)bufSignMask);
    const __m128 _lutLen = _mm_load_ps((const float*)bufLutLen);

    S32 i = 0;

#pragma omp parallel for
    for (i = 0; i < dstHeight; i++)
    {
        S32 ALIGN16 buf[4];
        F32 ALIGN16 bufSum[4];
        S32 j = 0;
        S32 k = 0;
        F32 val0 = 0.f;
        F32 sum = 0.f;
        F32 wsum = 0.f;

        U16 *psrc = src + (radius + i) * srcStep;
        U16 *pdst = dst + i * dstStep;

        __m128 _val;
        __m128 _val0;
        __m128 _idx;
        __m128 _psum;
        __m128 _sw;
        __m128 _cw;
        __m128 _w;

        for (j = radius; j < xEnd; j++)
        {
            val0 = psrc[j] * 1.0f;
            _psum = _mm_setzero_ps();
            _val0 = _mm_set1_ps(val0);

            for (k = 0; k <= maxk - 4; k += 4)
            {
                _val = _mm_set_ps(psrc[j + spOfs[k + 3]], psrc[j + spOfs[k + 2]],
                    psrc[j + spOfs[k + 1]], psrc[j + spOfs[k]]);
                //                 _sw = _mm_loadu_ps(spWt + k);
                _sw = _mm_load_ps(spWt + k);
                _idx = _mm_and_ps(_signMask, _mm_sub_ps(_val, _val0));
                _mm_store_si128((__m128i*)buf, _mm_cvtps_epi32(_mm_min_ps(_idx, _lutLen)));
                _cw = _mm_set_ps(valWt[buf[3]], valWt[buf[2]], valWt[buf[1]], valWt[buf[0]]);
                _w = _mm_mul_ps(_cw, _sw);
                _val = _mm_mul_ps(_w, _val);
                _sw = _mm_hadd_ps(_w, _val);
                _sw = _mm_hadd_ps(_sw, _sw);
                _psum = _mm_add_ps(_sw, _psum);
            }

            _mm_storel_pi((__m64*)bufSum, _psum);
            sum = bufSum[1] + val0;
            wsum = bufSum[0] + 1.f;

            pdst[j - radius] = (U16)(sum / wsum);
        }
    }
}

// no smooth on the border
void edgePreserve_mainBody(U16 *src, U16 *dst, S32 width, S32 height,
    S32 radius, F32 sigmaVal, F32 sigmaSp)
{
    S32 i = 0;
    S32 j = 0;
    S32 maxk = -1;    //exclude the center

    S32 *spOfs = NULL;
    F32 *spWt = NULL;
    F32 *valWt = NULL;

    for (i = -radius; i <= radius; i++)
    {
        for (j = -radius; j <= radius; j++)
        {
            if (sqrtf((i * i + j * j) * 1.f) <= radius)
            {
                maxk++;
            }
        }
    }

    spWt = (F32 *)ALIGN_MALLOC16(sizeof(F32) * maxk);
    spOfs = (S32 *)MALLOC(sizeof(S32) * maxk);
    valWt = (F32 *)MALLOC(sizeof(F32) * BF_BUF_LEN);

    edgePreserve_LUT(radius, width, sigmaVal, sigmaSp, spOfs, spWt, valWt);

    borderCopy(src, dst, width, height, radius);

    edgePreserve_mainBody_process(src, width, height, width,
        dst + width * radius + radius, width,
        radius, spOfs, spWt, valWt, maxk);

    ALIGN_FREE(spWt);
    FREE(spOfs);
    FREE(valWt);
}

// smooth on the border
void edgePreserve(U16 *src, U16 *dst, S32 width, S32 height,
    S32 radius, F32 sigmaVal, F32 sigmaSp)
{
    S32 i = 0;
    S32 j = 0;
    S32 maxk = -1;    //exclude the center
    S32 width2 = width + (radius << 1);
    S32 height2 = height + (radius << 1);

    S32 *spOfs = NULL;
    F32 *spWt = NULL;
    F32 *valWt = NULL;

    U16 *buf = (U16 *)MALLOC(sizeof(U16) * width2 * height2);

    borderReflect(src, width, height, buf, radius);

    for (i = -radius; i <= radius; i++)
    {
        for (j = -radius; j <= radius; j++)
        {
            if (sqrtf((i * i + j * j) * 1.f) <= radius)
            {
                maxk++;
            }
        }
    }

    spWt = (F32 *)ALIGN_MALLOC16(sizeof(F32) * maxk);
    spOfs = (S32 *)MALLOC(sizeof(S32) * maxk);
    valWt = (F32 *)MALLOC(sizeof(F32) * BF_BUF_LEN);

    edgePreserve_LUT(radius, width2, sigmaVal, sigmaSp, spOfs, spWt, valWt);

    edgePreserve_mainBody_process(buf, width2, height2, width2,
        dst, width, radius, spOfs, spWt, valWt, maxk);

    ALIGN_FREE(spWt);
    FREE(spOfs);
    FREE(valWt);
    FREE(buf);
}

// no smooth on the border
void edgePreserve_mainBody_sse(U16 *src, U16 *dst, S32 width, S32 height,
    S32 radius, F32 sigmaVal, F32 sigmaSp)
{
    S32 i = 0;
    S32 j = 0;
    S32 maxk = -1;    //exclude the center

    S32 *spOfs = NULL;
    F32 *spWt = NULL;
    F32 *valWt = NULL;

    for (i = -radius; i <= radius; i++)
    {
        for (j = -radius; j <= radius; j++)
        {
            if (sqrtf((i * i + j * j) * 1.f) <= radius)
            {
                maxk++;
            }
        }
    }

    spWt = (F32 *)ALIGN_MALLOC16(sizeof(F32) * maxk);
    spOfs = (S32 *)MALLOC(sizeof(S32) * maxk);
    valWt = (F32 *)MALLOC(sizeof(F32) * BF_BUF_LEN);

    edgePreserve_LUT(radius, width, sigmaVal, sigmaSp, spOfs, spWt, valWt);

    borderCopy(src, dst, width, height, radius);

    //     edgePreserve_noBorder_sse_mainloop(src, dst, width, height, radius, spOfs, spWt, valWt, maxk);

    edgePreserve_mainBody_sse_process(src, width, height, width,
        dst + width * radius + radius, width,
        radius, spOfs, spWt, valWt, maxk);

    ALIGN_FREE(spWt);
    FREE(spOfs);
    FREE(valWt);
}

// smooth on the border
void edgePreserve_sse(U16 *src, U16 *dst, S32 width, S32 height,
    S32 radius, F32 sigmaVal, F32 sigmaSp)
{
    S32 i = 0;
    S32 j = 0;
    S32 maxk = -1;    //exclude the center
    S32 width2 = width + (radius << 1);
    S32 height2 = height + (radius << 1);

    S32 *spOfs = NULL;
    F32 *spWt = NULL;
    F32 *valWt = NULL;

    U16 *buf = (U16 *)MALLOC(sizeof(U16) * width2 * height2);

    borderReflect(src, width, height, buf, radius);

    for (i = -radius; i <= radius; i++)
    {
        for (j = -radius; j <= radius; j++)
        {
            if (sqrtf((i * i + j * j) * 1.f) <= radius)
            {
                maxk++;
            }
        }
    }

    spWt = (F32 *)ALIGN_MALLOC16(sizeof(F32) * maxk);
    spOfs = (S32 *)MALLOC(sizeof(S32) * maxk);
    valWt = (F32 *)MALLOC(sizeof(F32) * BF_BUF_LEN);

    edgePreserve_LUT(radius, width2, sigmaVal, sigmaSp, spOfs, spWt, valWt);

    //     edgePreserve_sse_mainloop(buf, dst, width, height, radius, spOfs, spWt, valWt, maxk);

    edgePreserve_mainBody_sse_process(buf, width2, height2, width2,
        dst, width, radius, spOfs, spWt, valWt, maxk);


    ALIGN_FREE(spWt);
    FREE(spOfs);
    FREE(valWt);
    FREE(buf);
}


typedef struct EDGE_PRES_FILT_HDL
{
    S32 maxk;
    S32 width;
    S32 height;
    S32 radius;
    S32 *spOfs;
    F32 *spWt;
    F32 *valWt;
    U16 *buf;
}EDGE_PRES_FILT_HDL_;

edgePresFiltMain::edgePresFiltMain(S32 width, S32 height, S32 radius, F32 sigmaVal, F32 sigmaSp)
{
    S32 i = 0;
    S32 j = 0;
    S32 maxk = -1;    //exclude the center
    EDGE_PRES_FILT_HDL *pHdl = NULL;

    pHdl = new EDGE_PRES_FILT_HDL;
    hdl = pHdl;

    for (i = -radius; i <= radius; i++)
    {
        for (j = -radius; j <= radius; j++)
        {
            if (sqrtf((i * i + j * j) * 1.f) <= radius)
            {
                maxk++;
            }
        }
    }

    pHdl->maxk = maxk;
    pHdl->width = width;
    pHdl->height = height;
    pHdl->radius = radius;
    pHdl->spWt = (F32 *)ALIGN_MALLOC16(sizeof(F32) * maxk);
    pHdl->spOfs = (S32 *)MALLOC(sizeof(S32) * maxk);
    pHdl->valWt = (F32 *)MALLOC(sizeof(F32) * BF_BUF_LEN);

    edgePreserve_LUT(radius, width, sigmaVal, sigmaSp, pHdl->spOfs, pHdl->spWt, pHdl->valWt);
}

edgePresFiltMain::~edgePresFiltMain()
{
    EDGE_PRES_FILT_HDL *pHdl = (EDGE_PRES_FILT_HDL *)hdl;

    ALIGN_FREE(pHdl->spWt);
    FREE(pHdl->spOfs);
    FREE(pHdl->valWt);

    delete hdl;
}

void edgePresFiltMain::edgePreserve_mainBody(U16 *src, U16 *dst)
{
    EDGE_PRES_FILT_HDL *pHdl = (EDGE_PRES_FILT_HDL *)hdl;
    S32 width = pHdl->width;
    S32 height = pHdl->height;
    S32 radius = pHdl->radius;

    borderCopy(src, dst, width, height, radius);

    edgePreserve_mainBody_process(src, width, height, width,
        dst + width * radius + radius, width,
        radius, pHdl->spOfs, pHdl->spWt, pHdl->valWt, pHdl->maxk);
}

void edgePresFiltMain::edgePreserve_mainBody_omp(U16 *src, U16 *dst)
{
    EDGE_PRES_FILT_HDL *pHdl = (EDGE_PRES_FILT_HDL *)hdl;
    S32 width = pHdl->width;
    S32 height = pHdl->height;
    S32 radius = pHdl->radius;

    borderCopy(src, dst, width, height, radius);

    edgePreserve_mainBody_omp_process(src, width, height, width,
        dst + width * radius + radius, width,
        radius, pHdl->spOfs, pHdl->spWt, pHdl->valWt, pHdl->maxk);
}

void edgePresFiltMain::edgePreserve_mainBody_sse(U16 *src, U16 *dst)
{
    EDGE_PRES_FILT_HDL *pHdl = (EDGE_PRES_FILT_HDL *)hdl;
    S32 width = pHdl->width;
    S32 height = pHdl->height;
    S32 radius = pHdl->radius;

    borderCopy(src, dst, width, height, radius);

    edgePreserve_mainBody_sse_process(src, width, height, width,
        dst + width * radius + radius, width,
        radius, pHdl->spOfs, pHdl->spWt, pHdl->valWt, pHdl->maxk);
}

void edgePresFiltMain::edgePreserve_mainBody_sse_omp(U16 *src, U16 *dst)
{
    EDGE_PRES_FILT_HDL *pHdl = (EDGE_PRES_FILT_HDL *)hdl;
    S32 width = pHdl->width;
    S32 height = pHdl->height;
    S32 radius = pHdl->radius;

    borderCopy(src, dst, width, height, radius);

    edgePreserve_mainBody_sse_omp_process(src, width, height, width,
        dst + width * radius + radius, width,
        radius, pHdl->spOfs, pHdl->spWt, pHdl->valWt, pHdl->maxk);
}

edgePresFilt::edgePresFilt(S32 width, S32 height, S32 radius, F32 sigmaVal, F32 sigmaSp)
{
    S32 i = 0;
    S32 j = 0;
    S32 maxk = -1;    //exclude the center
    S32 width2 = width + (radius << 1);
    S32 height2 = height + (radius << 1);
    EDGE_PRES_FILT_HDL *pHdl = NULL;

    pHdl = new EDGE_PRES_FILT_HDL;
    hdl = pHdl;

    for (i = -radius; i <= radius; i++)
    {
        for (j = -radius; j <= radius; j++)
        {
            if (sqrtf((i * i + j * j) * 1.f) <= radius)
            {
                maxk++;
            }
        }
    }

    pHdl->maxk = maxk;
    pHdl->width = width;
    pHdl->height = height;
    pHdl->radius = radius;
    pHdl->spWt = (F32 *)ALIGN_MALLOC16(sizeof(F32) * maxk);
    pHdl->spOfs = (S32 *)MALLOC(sizeof(S32) * maxk);
    pHdl->valWt = (F32 *)MALLOC(sizeof(F32) * BF_BUF_LEN);
    pHdl->buf = (U16 *)MALLOC(sizeof(U16) * width2 * height2);

    edgePreserve_LUT(radius, width2, sigmaVal, sigmaSp, pHdl->spOfs, pHdl->spWt, pHdl->valWt);
}

edgePresFilt::~edgePresFilt()
{
    EDGE_PRES_FILT_HDL *pHdl = (EDGE_PRES_FILT_HDL *)hdl;

    ALIGN_FREE(pHdl->spWt);
    FREE(pHdl->spOfs);
    FREE(pHdl->valWt);
    FREE(pHdl->buf);

    delete hdl;
}

void edgePresFilt::edgePreserve(U16 *src, U16 *dst)
{
    EDGE_PRES_FILT_HDL *pHdl = (EDGE_PRES_FILT_HDL *)hdl;
    S32 width = pHdl->width;
    S32 height = pHdl->height;
    S32 radius = pHdl->radius;
    S32 width2 = width + (radius << 1);
    S32 height2 = height + (radius << 1);

    borderReflect(src, width, height, pHdl->buf, radius);

    edgePreserve_mainBody_process(pHdl->buf, width2, height2, width2,
        dst, width, radius, pHdl->spOfs, pHdl->spWt, pHdl->valWt, pHdl->maxk);
}

void edgePresFilt::edgePreserve_omp(U16 *src, U16 *dst)
{
    EDGE_PRES_FILT_HDL *pHdl = (EDGE_PRES_FILT_HDL *)hdl;
    S32 width = pHdl->width;
    S32 height = pHdl->height;
    S32 radius = pHdl->radius;
    S32 width2 = width + (radius << 1);
    S32 height2 = height + (radius << 1);

    borderReflect(src, width, height, pHdl->buf, radius);

    edgePreserve_mainBody_omp_process(pHdl->buf, width2, height2, width2,
        dst, width, radius, pHdl->spOfs, pHdl->spWt, pHdl->valWt, pHdl->maxk);
}

void edgePresFilt::edgePreserve_sse(U16 *src, U16 *dst)
{
    EDGE_PRES_FILT_HDL *pHdl = (EDGE_PRES_FILT_HDL *)hdl;
    S32 width = pHdl->width;
    S32 height = pHdl->height;
    S32 radius = pHdl->radius;
    S32 width2 = width + (radius << 1);
    S32 height2 = height + (radius << 1);

    borderReflect(src, width, height, pHdl->buf, radius);

    edgePreserve_mainBody_sse_process(pHdl->buf, width2, height2, width2,
        dst, width, radius, pHdl->spOfs, pHdl->spWt, pHdl->valWt, pHdl->maxk);
}

void edgePresFilt::edgePreserve_sse_omp(U16 *src, U16 *dst)
{
    EDGE_PRES_FILT_HDL *pHdl = (EDGE_PRES_FILT_HDL *)hdl;
    S32 width = pHdl->width;
    S32 height = pHdl->height;
    S32 radius = pHdl->radius;
    S32 width2 = width + (radius << 1);
    S32 height2 = height + (radius << 1);

    borderReflect(src, width, height, pHdl->buf, radius);

    edgePreserve_mainBody_sse_omp_process(pHdl->buf, width2, height2, width2,
        dst, width, radius, pHdl->spOfs, pHdl->spWt, pHdl->valWt, pHdl->maxk);
}

//

static void edgePreserveInt_LUT(S32 radius, S32 width, F32 sigmaVal, F32 sigmaSp,
    S32 *spOfs, S32 *spWt, S32 *valWt)
{
    S32 i = 0;
    S32 j = 0;
    S32 k = 0;
    F32 sigmaValCoeff = -0.5f / (sigmaVal * sigmaVal);
    F32 sigmaSpCoeff = -0.5f / (sigmaSp * sigmaSp);

    for (i = -radius; i <= radius; i++)
    {
        for (j = -radius; j <= radius; j++)
        {
            if (sqrtf((i * i + j * j) * 1.f) <= radius)
            {
                if ((i == 0) && (j == 0))
                {
                    continue;
                }

                spOfs[k] = i * width + j;
                spWt[k] = (S32)(expf((i * i + j * j) * sigmaSpCoeff) * BF_INT_SCALE);
                k++;
            }
        }
    }

    for (i = 0; i < BF_BUF_LEN - 1; i++)
    {
        valWt[i] = (S32)(expf(i * i * sigmaValCoeff) * BF_INT_SCALE);
    }
    valWt[i] = 0;
}

// Smooth main body of src to the dst img
// src img is bigger than dst img
static void edgePreserveInt_mainBody_process(U16 *src, S32 srcWidth, S32 srcHeight, S32 srcStep,
    U16 *dst, S32 dstStep, S32 radius, S32 *spOfs, S32 *spWt, S32 *valWt, S32 maxk)
{
    const S32 xEnd = srcWidth - radius;
    const S32 dstHeight = srcHeight - (radius << 1);

    S32 i = 0;
    S32 j = 0;
    S32 k = 0;
    S32 n = 0;
    U16 val0 = 0;
    U16 val = 0;
    S32 tmp = 0;
    S32 w = 0;
    S32 sum = 0;
    S32 wsum = 0;

    U16 *psrc = src + radius * srcStep;
    U16 *pdst = dst;

    for (i = 0; i < dstHeight; i++)
    {
        for (j = radius; j < xEnd; j++)
        {
            val0 = psrc[j];
            sum = 0;
            wsum = 0;

            for (k = 0; k < maxk; k++)
            {
                val = psrc[j + spOfs[k]];
                tmp = val - val0;
                tmp = ABS(tmp);
                tmp = MIN(tmp, BF_BUF_LEN - 1);
                w = spWt[k] * valWt[tmp];
                w >>= BF_INT_SHIFT;
                wsum += w;
                sum += (val * w);
            }

            pdst[j - radius] = (U16)((sum + (val0 << BF_INT_BITS2)) / (wsum + BF_INT_SCALE2));
        }

        psrc += srcStep;
        pdst += dstStep;
    }
}

// smooth on the border
void edgePreserveInt(U16 *src, U16 *dst, S32 width, S32 height,
    S32 radius_, F32 sigmaVal, F32 sigmaSp)
{
    S32 i = 0;
    S32 j = 0;
    S32 maxk = -1;    //exclude the center
    S32 radius = GP_MIN(radius_, GP_EDGEPRES_R_MAX);
    S32 width2 = width + (radius << 1);
    S32 height2 = height + (radius << 1);

    S32 *spOfs = NULL;
    S32 *spWt = NULL;
    S32 *valWt = NULL;

    U16 *buf = (U16 *)MALLOC(sizeof(U16) * width2 * height2);

    borderReflect(src, width, height, buf, radius);

    for (i = -radius; i <= radius; i++)
    {
        for (j = -radius; j <= radius; j++)
        {
            if (sqrtf((i * i + j * j) * 1.f) <= radius)
            {
                maxk++;
            }
        }
    }

    spWt = (S32 *)ALIGN_MALLOC16(sizeof(S32) * maxk);
    spOfs = (S32 *)MALLOC(sizeof(S32) * maxk);
    valWt = (S32 *)MALLOC(sizeof(S32) * BF_BUF_LEN);

    edgePreserveInt_LUT(radius, width2, sigmaVal, sigmaSp, spOfs, spWt, valWt);

    edgePreserveInt_mainBody_process(buf, width2, height2, width2,
        dst, width, radius, spOfs, spWt, valWt, maxk);

    ALIGN_FREE(spWt);
    FREE(spOfs);
    FREE(valWt);
    FREE(buf);
}

typedef struct EDGE_PRES_FILT_INT_HDL
{
    S32 maxk;
    S32 width;
    S32 height;
    S32 radius;
    S32 *spOfs;
    S32 *spWt;
    S32 *valWt;
    U16 *buf;
}EDGE_PRES_FILT_INT_HDL_;

edgePresFiltInt::edgePresFiltInt(S32 width, S32 height, S32 radius_, F32 sigmaVal, F32 sigmaSp)
{
    S32 i = 0;
    S32 j = 0;
    S32 maxk = -1;    //exclude the center
    S32 radius = GP_MIN(radius_, GP_EDGEPRES_R_MAX);
    S32 width2 = width + (radius << 1);
    S32 height2 = height + (radius << 1);
    EDGE_PRES_FILT_INT_HDL *pHdl = NULL;

    pHdl = new EDGE_PRES_FILT_INT_HDL;
    hdl = pHdl;

    for (i = -radius; i <= radius; i++)
    {
        for (j = -radius; j <= radius; j++)
        {
            if (sqrtf((i * i + j * j) * 1.f) <= radius)
            {
                maxk++;
            }
        }
    }

    pHdl->maxk = maxk;
    pHdl->width = width;
    pHdl->height = height;
    pHdl->radius = radius;
    pHdl->spWt = (S32 *)ALIGN_MALLOC16(sizeof(S32) * maxk);
    pHdl->spOfs = (S32 *)MALLOC(sizeof(S32) * maxk);
    pHdl->valWt = (S32 *)MALLOC(sizeof(S32) * BF_BUF_LEN);
    pHdl->buf = (U16 *)MALLOC(sizeof(U16) * width2 * height2);

    edgePreserveInt_LUT(radius, width2, sigmaVal, sigmaSp, pHdl->spOfs, pHdl->spWt, pHdl->valWt);
}

edgePresFiltInt::~edgePresFiltInt()
{
    EDGE_PRES_FILT_INT_HDL *pHdl = (EDGE_PRES_FILT_INT_HDL *)hdl;

    ALIGN_FREE(pHdl->spWt);
    FREE(pHdl->spOfs);
    FREE(pHdl->valWt);
    FREE(pHdl->buf);

    delete hdl;
}

void edgePresFiltInt::edgePreserve(U16 *src, U16 *dst)
{
    EDGE_PRES_FILT_INT_HDL *pHdl = (EDGE_PRES_FILT_INT_HDL *)hdl;
    S32 width = pHdl->width;
    S32 height = pHdl->height;
    S32 radius = pHdl->radius;
    S32 width2 = width + (radius << 1);
    S32 height2 = height + (radius << 1);

    borderReflect(src, width, height, pHdl->buf, radius);

    edgePreserveInt_mainBody_process(pHdl->buf, width2, height2, width2,
        dst, width, radius, pHdl->spOfs, pHdl->spWt, pHdl->valWt, pHdl->maxk);
}





billbliss
关注
  • 0
    点赞
  • 踩
  • 4
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
双边滤波快速算法-图像处理方向
05-10
双边滤波快速算法-图像处理方向 介绍了双边滤波的概念-很精练形象 以及在图像处理中的应用-edge preserving filtering 介绍了现有的多种加速算法
增维型快速双边滤波-A Fast Approximation of the Bilateral Filter using a Signal Processing Approach
Nick Blog
12-17 2017
1. 简介 双边滤波非常有用，但速度很慢，因为它是非线性的，传统的加速算法例如在FFT之后执行卷积，是不适用的。本文提出了对双边滤波的新解释，即高维卷积，然后是两个非线性操作。其基本思想就是将非线性的双边滤波改成可分离的线性操作和非线性操作。换句话说，原来的双边滤波在图像不同位置应用不同的权重，也就是位移改变卷积，他们通过增加一个维度，也就是将灰度值作为一个新的维度，将双边滤波表达成3D空间中的线性位移不变卷积，最后再执行非线性的归一化操作。 2. 原理 2.1 公式推导 （1） 首先将原始双边滤波公式等
图像处理之双边滤波
yannkaii的博客
06-08 393
值域核通过计算模板窗口中每个位置对应的灰度值与模板窗口中心对应的灰度值之间的差值来分配该位置的权重，如果该位置对应灰度与中心坐标对应灰度差异较小，则权重较大，反之则权值较小，这也是双边滤波能够在滤除噪声同时保留边缘的关键所在。空间域核与高斯滤波核类似，主要通过模板窗口中各位置与模板窗口中心的欧式距离给模板窗口中每一个位置计算分配一个权值，权值分布为模板窗口中心大边缘小（具体权值分布模型可参考高斯滤波器）。上式中为模板窗口中心位置灰度，为模板窗口中其他位置灰度；从滤波前后的效果来看，噪声被去除，边缘被保留。
快速双边滤波在High-Dynamic Range（HDR)图像显示中的应用
big tiger的博客
03-05 3356
快速双边滤波在High-Dynamic Range（HDR)图像显示中的应用。                分类：            图像处理2013-11-15 20:38628人阅读评论(3)收藏举报 bilateral filter快速双边滤波Hdr      高动态范围图像（High-Dynamic Range，简称HDR），相比普通的图像，可以提供更多的动态范围和
PCL快速双边滤波（fast_bilateral、fast_bilateral_omp）
最新发布
xhm01291212的博客
07-30 324
的速度较慢。快速双边滤波使用提升双边滤波的速度。（1）对进行。（2）同双边滤波一样，也是先计算，得到的空间领域内的所有点。（3）在该领域范围内，计算操作点和空间邻域内所有点的强度差异是否在高斯权重的范围内，如果不在高斯权重的范围内，那么使用高斯权重平滑强度。（4）得到最终结果。
快速双边滤波在High-Dynamic Range（HDR)图像显示中的应用。
凌风探梅的专栏
04-15 2439
快速双边滤波在High-Dynamic Range（HDR)图像显示中的应用。 2014年10月21日 ⁄ 综合 ⁄ 共 4617字 ⁄ 字号 小 中 大 ⁄ 评论关闭      高动态范围图像（High-Dynamic Range，简称HDR），相比普通的图像，可以提供更多的动态范围和图像细节，根据不同的曝光时间的LDR(Low-Dynamic
双边滤波器（bilateral filter）
12-20
- **并行计算**：利用GPU的并行计算能力，加速滤波过程。 - **多尺度策略**：在多个尺度上应用双边滤波，以适应不同场景的需求。 **5. 限制与挑战** 虽然双边滤波器在很多方面表现出色，但也存在一些限制，如处理...
双边滤波器C++ 代码
02-23
1. **并行处理**：利用多核CPU的并行计算能力，如OpenMP，来加速滤波过程。 2. **内存管理**：合理使用内存，避免不必要的数据拷贝。 3. **参数调整**：根据实际图像内容和需求，适当调整滤波器的直径和标准差，以...
双边滤波器
06-15
- **并行处理与硬件加速**：为了进一步提升双边滤波器的性能，可以通过利用并行处理平台（如GPU）和可重构硬件（如FPGA）来加速计算过程。 **4. 应用领域** 双边滤波器因其独特的特性，在多个领域都有广泛应用，...
OpenACC并行编程实战
07-22
资源名称： OpenACC并行编程实战资源截图： 资源太大，传百度网盘了，链接在附件中，有需要的同学自取。
Fast Bilateral Filt​er:快速 O (1) 双边滤波器-matlab开发
06-01
以下论文中描述的快速 O(1) 双边滤波器的 MATLAB 实现： [1] KN Chaudhury、D. Sage 和 M. Unser，“使用三角范围内核的快速 O(1) 双边滤波”，IEEE Trans。 图像处理，卷。 20，没有。 2011 年 11 月。 [2] KN Chaudhury，“用于双边滤波和非局部均值的可移动 O(1) 算法的加速”，arXiv：1203.5128v1。
Fast Bilateral Filtering for the Display of High-Dynamic-Range Image
02-27
提出一种快速双边滤波，并用其来作图像的色调映射。
一些滤波器（boxfilter，快速双边滤波器，快速引导滤波器和 permutohedral 双边滤波器）C++
07-12
双边滤波器、permutohedral双边滤波器和传播滤波器 箱式过滤器 快速导向过滤器 快速双边滤波器 多面体双边滤波器 传播过滤器 更多详情、使用方法，请下载后阅读README.md文件
快速双边滤波的美颜算法
02-25
本文实现了基于快速双边滤波得美颜算法。详细信息可以查看博客：https://blog.csdn.net/zhangpengzp/article/details/87884584
快速双边滤波matlab代码
08-20
双边滤波的逐段线性近似，参考文献Fast Bilateral Filtering for the Display of High-Dynamic-Range Images
3D点云处理：双边滤波器&快速双边滤波器（附源码）
lddx_123456的博客
09-22 854
3D点云处理：双边滤波器&快速双边滤波器。
c语言双边滤波算法,快速双边滤波 附完整C代码
weixin_34569101的博客
05-23 1281
很早之前写过《双边滤波算法的简易实现bilateralFilter》。当时学习参考的代码来自cuda的样例。相关代码可以参阅:https://github.com/johng12/cudaSamples/tree/master/cudaSamples/3_Imaging/bilateralFilter由于算法逻辑非常清晰，就不多解释了。需要补课的，请移步《o(1)复杂度之双边滤波算法的原理、流程、...
技术的真相 | 基于双边滤波的磨皮算法及优化
旷视开发者技术社区
09-06 1516
一．背景简介现在视频类应用非常火热，直播、美妆、医美应用层出不穷。用户们在使用这类应用时都希望自己在屏幕上的样子美美的，皮肤细腻光滑。本文就介绍一种实现简单、效果很好的磨皮算法以及对它的优...
图像处理（六）递归双边滤波磨皮
Bob的专栏
10-27 3639
递归双边滤波是双边滤波的一种加速算法，加速比非常大，就像高斯模糊快速算法一样，加速起来，速度跟原算法相比，提高了十几倍。双边滤波的加速算法很多，文献都比较老旧，我这边主要讲一下比较新的算法：递归双边滤波，对应的paper为：《Recursive Bilateral Filtering》 这个算法比我另外一篇博文将的双指数滤波快一点，两篇文献的基本思想一样： 递归公式跟双指数的递归公
基于fpga的双边滤波器设计
07-29
1. 实时性能：FPGA的并行计算能力使得双边滤波器设计可以在实时应用中使用，如实时图像处理、视频传输等。 2. 灵活性：FPGA的可编程性使得双边滤波器设计可以根据实际需求进行优化和修改。 3. 高性能：FPGA的计算...
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值