从相关卷积到图像高斯模糊的快速计算

博客主要介绍了信号相关、线性卷积、圆周卷积的概念及计算方法。相关用于衡量信号相似程度,线性卷积有公式及物理解释,圆周卷积有多种求法。还阐述了利用圆周卷积计算线性卷积的方法,介绍了FFT快速卷积和IIR型高斯滤波器。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

几年前的总结 ,在此记录一下:

一 相关

相关是用来衡量两个信号  "形状" 和 "相位"(起始位置)的相似程度 . 

计算相关性的两个信号必须要有相同起始位置和长度,然后对两个信号逐点乘积加和,当两个信号完全一致时相关性是最高的.

离散信号相关计算定义:

信号中的直流分量通常对信号表达是没有意义的,反而会产生干扰。所以某些相关性计算要先对信号做减直流处理,如数理统计里用协方差来度量相关性: 

 

相关的计算是满足交换律的, "x 与 h的相关性" 等价于 "h与x的相关性" , 在这里我们把 h看成特征模板,x看成未知信号(下文卷积也都按这种约定).从这个角度来看

相关性计算也表达了"模板信号对未知信号加权和"的意思.比如要从一个很长的序列里匹配一段已知信号,相关计算算法可能从任意起始位置截取一段 x(0 ~ n) , 与 h 做加权和,得分最高的位置即是最匹配的信号.

 

二卷积(通常是指线卷积)

卷积与相关类似,按字面理解,先卷(与相关的主要区别),再积,离散卷积公式:

公式如下:

                先把模板h翻卷(按原点对褶) h(i) -> h(-i)

                h的零点对齐(平移) 到x的目标点n  h(-i) -> h(n - i)

                再对两段序列求逐点积,加和

当卷积模板中心对称时,卷积的结果与相关是一致的!!

同样的卷积计算满足交换率,即 (x * h)(n)  <==> (h * x)(n) 

卷积的物理解释:

a  : 假如一幅只有不同黑色程度的水墨画下雨泡了水。  x(p) 表示水墨画在 p 点的含墨量,h(delta) 表示字画泡水后,单位墨向位置 (p + delta)扩散了多少 ( 实际扩散了 x(p) * h(p + delta)  , 线性时不变系统) , 问题求是这幅画最终每一点的墨迹量。  注意, x 的定义域为整个画的二维坐标, h的定义域为与x,y平面原点中心对称的矩形区域。

 

    显然用卷积的定义求解就很容易, 目标点 P1受到模板h定义的所有邻域点影响,对于邻域点 Q1 , 位移了delta 到 P1 , 则 Q1 + delta = P1 ,即受Q1的影响为 x(Q1) * h (P1 - Q1) ,  p1 点收到总的影响为:

b  : 假如小明在t时刻脸被挨打强度为x(t),  受单位强度的力挨打后, 脸鼓起的包随时间变化函数为 h(t) , 那问题来了, 求{0 , ... N}任意时刻小明脸鼓起的包强度。  显然 x(t) 和h(t) 的定义域只有时间正轴,包括零点。

信号里 x 被称为冲激函数, h 为冲激响应,这里冲激函数只会对后面的时刻产生响应。挨打之前脸上就有包那是不负责任的! 同样计算公式如下:

卷积计算的几种方法:             

a 按公式 反褶(与零点)、平移(h 0点到目标位置n) , 逐点乘积 再求和

 

b 多项式方法(竖式计算法则)一种快捷计算方式

X1 = 

X2 =  

 

如图所示,竖式计算时,两个序列要做到齐次项对齐,即定义域零点对齐。

 

三圆周卷积

圆周卷积也称作循环卷积,百度百科里的定义如下图:

 

两个序列x,h (序列下标从0开始) , 作L点的圆周卷积:

        1    首先将x ,h补成L长度的序列(少则补0,多则截断) 

        2    模板序列h逆时针排成一圈,x 顺时针排成一圈 

        3    求 0 点的圆周卷积对应位置相乘加和

        4    求 点n 的圆周圈积,将模板h 顺时针转n下,再逐点乘积求和。

 

也来个定义:

 

假设两个序列长度分别为 M,N , L >= M + N - 1 时 圆周卷积和线性卷积值是一致的!! 

 

圆周卷积求法2:  竖式计算

X = {1 , 4  , 3 } ,  H = {1 , 2 , 3}

方法是用H的每一个元素去乘  X 每一组顺时针循环移位的序列,再求和:

                        1 X    {1 , 4 ,  3}  +  2 X   {3 , 1 , 4}  + 3 X {4 , 3 , 1}  =          1          4            3

                                                                                                                          6          2           8         

                                                                                                                          12        9           3

                                                                                                                    ---------------------------------

                                                                                                                         19        15         14

圆周卷积求法3:  矩阵乘法

X = {1 , 4  , 3 } ,  H = {1 , 2 , 3}

 

                                         0

先将 H周期延拓(1,2,3 , 1,2,3 , 1,2,3),沿0点反转 , 再依次求每组顺时针循环右称的序列

HM = 

    |    1    3     2    |

    |    2    1     3    |

    |    3    2     1    |

F列向量 =    = {19 , 15 , 14}

圆周卷积求法4:  DFT 

经证明圆周卷积与 两个序列先DFT再乘积再逆DFT结果是一致。个人理解圆周卷积只是计算上有DFT 这种快速算法,没有线卷积那么强的物理意义,L = N + M - 1 时二者实现了统一。 L < N + M - 1 时会造成混叠效应。

 

四 利用圆周卷积来计算线性卷积

4.1 模板定义域为正轴数据 

a 将 x & h 扩展成 L = N + M - 1 序列 

b 计算圆周卷积。

c 再截掉 最后  2 *floor(M / 2) = M - 1个值 

 

4.2 模板h定义域包含负值 

这里举一个 h 中心对称的例子,同样 x 序列长度为N , h 长度为M M为奇数, 且h 定义域为 [ -M / 2 , M/ 2]

a 将包含负值的定义域平移到0点变成 [0 , M],这样最终卷积的结果会延时(晚)  floor(M/2) 个点!

b 将 x & h 扩展成 L = N + M - 1 序列 

c 计算圆周圈积。 

4 截掉最前面 floor(M / 2) 和最后 (M / 2) 个点

 

4.3 线卷积matlab 代码 

function [s] = linear_cov(x , h)


N = size(x , 2);
M = size(h , 2);
L = N + M - 1;


% fprintf(2 , 'M = %d , N = %d , L = %d \n' , M ,  N , L);


s = zeros(1 , L);


for i=1:N
    for j=1:M
        s(j + i - 1) = s(j + i - 1) + x(i) * h(j);
    end
end   

测试:

X=[4,5,6,7,8,9,10];
H=[1,2,3];
s=linear_cov(X,H);
%s=conv(H,X);

4.4 圆周卷积matlab 代码

function [s] = circle_cov(x , h , L)


N = size(x , 2);
M = size(h , 2);


if N < L 
    X1 = [x , zeros(1 , L - N)];
else
    X1 = x(1:L);
end


if M < L
    H1 = [h , zeros(1 , L - M)];
else
    H1 = h(1:L);
end


s = zeros(1 , L);


for i=1:L
    for j=1:L
        s(i) = s(i) + X1(j) * H1(1 + mod(i - j , L));
    end
end

测试:

X=[4,5,6,7,8,9,10];
H=[1,2,3];
s=circle_cov(X,H , 9);

 

9点的圆周卷积结果与线卷积一致。

五 FFT 快速卷积

求x, h的卷积 matlab 代码: ( h (x) : x >= 0)

function [s] = fft_cov(x , h)


N = size(x , 2);
M = size(h , 2);
L = N + M - 1;


X = [x , zeros(1 , L - N)];
H = [h , zeros(1 , L - M)];


f_x = fft(X);
f_h = fft(H);


f_s = f_x .* f_h;
s   = ifft(f_s);

如果 h 定义域  x [-a , b] , a < 0 , 卷积输出会延时 a 长度。 对于二维图像卷积适用。

六 IIR 型高斯滤波器

参考:

https://www.cnblogs.com/luo-peng/p/5223910.html

IIR 型高斯滤波器特点是计算复杂度与 核 大小无关有一定参考意义

实现代码

//iir_guass_blur.h

#ifndef __IIR_GUASS_BLUR_H
#define __IIR_GUASS_BLUR_H

#ifdef __cplusplus
extern "C"
{
#endif

long IIR_guass_filter_create(int w , int h , float sigma);

void IIR_guass_filter_do(long  h_flt , const uint8_t * src_plane , int src_pitch , uint8_t *dst_plane , int dst_pitch);

void IIR_guass_filter_destroy(long  h_flt);

long IIR_guass_filterI_create(int w , int h , float sigma);

void IIR_guass_filterI_do(long  h_flt , const uint8_t * src_plane , int src_pitch , uint8_t *dst_plane , int dst_pitch);

void IIR_guass_filterI_destroy(long  h_flt);



#ifdef __cplusplus
}
#endif

#endif
//iir_guass_blur.c

#include <malloc.h>
#include <stdlib.h>
#include <string.h>
#include <stdint.h>
#include <assert.h>
#include <math.h>
#include "iir_guass_blur.h"

/* *** IIR guass blur formula : ***
 *  
 *  forward pass  : 
        w[n]  =  B * x[n] + (b1 * y[n - 1] + b2 * y[n - 2] + b3 * y[n - 3]) / b0;
  
 *  backword pass :
       out[n] =  B * w[n] + (b1 * out[n + 1] + b2 * out[n + 2] + b3 * out[n + 3]) / b0;

 *      B = 1 - (b1 + b2 + b3) / b0; 
 *
*/

typedef struct IIR_guass_coefs
{
    float   B;
    float   b[4];    
}IIR_guass_coefs;

typedef struct IIR_guass_flt_context
{
    IIR_guass_coefs     coeff;
    int                 w;
    int                 h;
    float               sigma;
    float               *w1;
    float               *w2;
    uint8_t             *dst_1pass;
}IIR_guass_flt_context;


static void compute_coefs3(IIR_guass_coefs *c , float sigma)  
{  
    float q, q2, q3;  
    if (sigma >= 2.5){  
        q = 0.98711 * sigma - 0.96330;  
    }  
    else if ((sigma >= 0.5) && (sigma < 2.5)){  
        q = 3.97156 - 4.14554 * (float) sqrtf ((float) 1 - 0.26891 * sigma);  
    }  
    else {  
        q = 0.1147705018520355224609375;  
    }  

    q2 = q * q;  
    q3 = q * q2;  
    c->b[0] = 1/(1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3));  
    c->b[1] = (        (2.44413*q)+(2.85619*q2)+(1.26661 *q3)) *c->b[0];  
    c->b[2] = (                   -((1.4281*q2)+(1.26661 *q3)))*c->b[0];  
    c->b[3] = (                                 (0.422205*q3)) *c->b[0];  
    c->B = 1.0-(c->b[1]+c->b[2]+c->b[3]);  
}  

long IIR_guass_filter_create(int w , int h , float sigma)
{
    IIR_guass_flt_context *ctx;
    int max_len;

    ctx = (IIR_guass_flt_context *)calloc(1 , sizeof(IIR_guass_flt_context));
    assert(ctx != NULL);

    ctx->w      = w;
    ctx->h      = h;
    ctx->sigma  = sigma;

    compute_coefs3(&ctx->coeff , sigma);

    max_len     = w > h ? w + 3 : h + 3;

    ctx->w1     = (float *)malloc(sizeof(float) * max_len);
    ctx->w2     = (float *)malloc(sizeof(float) * max_len);

    assert(ctx->w1 != 0);
    assert(ctx->w2 != 0);

    ctx->dst_1pass = (uint8_t *)malloc(sizeof(uint8_t) * h * w);
    assert(ctx->dst_1pass != 0);
    

    return (long)ctx;
    
}   

void IIR_guass_filter_onepass( IIR_guass_flt_context *ctx , int width , int height , const uint8_t * src_plane , int src_pitch , uint8_t *dst_plane , int dst_pitch)
{
    const uint8_t       *p_src;
    uint8_t             *p_dst;
    int n , h;
    float *w1 , *w2 , *b , B;

	p_src       = src_plane;
    p_dst		= dst_plane;
    w1          = ctx->w1;
    w2          = ctx->w2;
    b           = ctx->coeff.b;
    B           = ctx->coeff.B;

    for(h = 0 ; h < height ; ++h){
        /*forward process*/
        w1[0] = w1[1] = w1[2] = p_src[0];
        for(n = 3 ; n < width + 3 ; ++n){
            w1[n] = B * p_src[n - 3] + b[1] * w1[n - 1] + b[2] * w1[n - 2] + b[3] * w1[n - 3];
        }
        /*backward process*/
        w2[width] = w2[width + 1] = w2[width + 2] = w1[width + 2];
        for(n = width - 1; n >= 0 ; --n){
            w2[n] = B * w1[n + 3]    + b[1] * w2[n + 1] + b[2] * w2[n + 2] + b[3] * w2[n + 3];
        }
        /*store data*/
        for(n = 0 ; n < width ; ++n){
            p_dst[n * dst_pitch] = w2[n];
        }
        p_src += src_pitch;
        p_dst++;
    }
}
	
void IIR_guass_filter_do(long  h_flt , const uint8_t * src_plane , int src_pitch , uint8_t *dst_plane , int dst_pitch)
{
    IIR_guass_flt_context *ctx;

    ctx = (IIR_guass_flt_context *)h_flt;

    /*process horizontal data*/    
	IIR_guass_filter_onepass(ctx , ctx->w , ctx->h , src_plane , src_pitch , ctx->dst_1pass , ctx->h);
	
    /*process vertical data*/ 
	IIR_guass_filter_onepass(ctx , ctx->h , ctx->w , ctx->dst_1pass , ctx->h , dst_plane , dst_pitch);
	
}

void IIR_guass_filter_destroy(long  h_flt)
{
    IIR_guass_flt_context *ctx; 

    ctx = (IIR_guass_flt_context *)h_flt;

    if (ctx->dst_1pass){
        free(ctx->dst_1pass);
        ctx->dst_1pass = 0;
    }

    if (ctx->w1){
        free(ctx->w1);
        ctx->w1 = 0;
    }

    if (ctx->w2){
        free(ctx->w2);
        ctx->w2 = 0;
    }

    free(ctx);
}


typedef struct IIR_guass_coefsI
{
    int32_t   B;
    int32_t   b[4];    
}IIR_guass_coefsI;

typedef struct IIR_guass_fltI_context
{
    IIR_guass_coefsI    coeff;
    int                 w;
    int                 h;
    float               sigma;
    int32_t             *w1;
    int32_t             *w2;
    uint8_t             *dst_1pass;
}IIR_guass_fltI_context;


static void compute_coefs3_I(IIR_guass_coefsI *c , float sigma)  
{  
    double q, q2, q3;  
    if (sigma >= 2.5){  
        q = 0.98711 * sigma - 0.96330;  
    }  
    else if ((sigma >= 0.5) && (sigma < 2.5)){  
        q = 3.97156 - 4.14554 * (float) sqrtf ((float) 1 - 0.26891 * sigma);  
    }  
    else {  
        q = 0.1147705018520355224609375;  
    }  

    q2 = q * q;  
    q3 = q * q2;  
    c->b[0] = 16384 / (1.57825+(2.44413*q)+(1.4281 *q2)+(0.422205*q3));  
    c->b[1] = (        (2.44413*q)+(2.85619*q2)+(1.26661 *q3))  * c->b[0];  
    c->b[2] = (                   -((1.4281*q2)+(1.26661 *q3))) * c->b[0];  
    c->b[3] = (                                 (0.422205*q3))  * c->b[0];  
    c->B = 16384 - (c->b[1]+c->b[2]+c->b[3]);  
}  

long IIR_guass_filterI_create(int w , int h , float sigma)
{
    IIR_guass_fltI_context *ctx;
    int max_len;

    ctx = (IIR_guass_fltI_context *)calloc(1 , sizeof(IIR_guass_fltI_context));
    assert(ctx != NULL);

    ctx->w      = w;
    ctx->h      = h;
    ctx->sigma  = sigma;

    compute_coefs3_I(&ctx->coeff , sigma);

    max_len     = w > h ? w + 3 : h + 3;

    ctx->w1     = (int32_t *)malloc(sizeof(int32_t) * max_len);
    ctx->w2     = (int32_t *)malloc(sizeof(int32_t) * max_len);

    assert(ctx->w1 != 0);
    assert(ctx->w2 != 0);

    ctx->dst_1pass = (uint8_t *)malloc(sizeof(uint8_t) * h * w);
    assert(ctx->dst_1pass != 0);

    return (long)ctx;
    
}   

void IIR_guass_filterI_onepass( IIR_guass_fltI_context *ctx , int width , int height , const uint8_t * src_plane , int src_pitch , uint8_t *dst_plane , int dst_pitch)
{
    const uint8_t       *p_src;
    uint8_t             *p_dst;
    int n , h;
    int32_t *w1 , *w2 , *b , B;

	p_src       = src_plane;
    p_dst		= dst_plane;
    w1          = ctx->w1;
    w2          = ctx->w2;
    b           = ctx->coeff.b;
    B           = ctx->coeff.B;

    for(h = 0 ; h < height ; ++h){
        /*forward process*/
        w1[0] = w1[1] = w1[2] = p_src[0];
        for(n = 3 ; n < width + 3 ; ++n){
            w1[n] = B * p_src[n - 3] + b[1] * w1[n - 1] + b[2] * w1[n - 2] + b[3] * w1[n - 3];
        }
        /*backward process*/
        w2[width] = w2[width + 1] = w2[width + 2] = (w1[width + 2] >> 14);
        for(n = width - 1; n >= 0 ; --n){
            w2[n] = (B * (w1[n + 3] >> 14)) + b[1] * w2[n + 1] + b[2] * w2[n + 2] + b[3] * w2[n + 3];
        }
        /*store data*/
        for(n = 0 ; n < width ; ++n){
            p_dst[n * dst_pitch] = (w2[n] >> 14);
        }
        p_src += src_pitch;
        p_dst++;
    }
}
	
void IIR_guass_filterI_do(long  h_flt , const uint8_t * src_plane , int src_pitch , uint8_t *dst_plane , int dst_pitch)
{
    IIR_guass_fltI_context *ctx;

    ctx = (IIR_guass_fltI_context *)h_flt;

    /*process horizontal data*/    
	IIR_guass_filterI_onepass(ctx , ctx->w , ctx->h , src_plane , src_pitch , ctx->dst_1pass , ctx->h);
	
    /*process vertical data*/ 
	IIR_guass_filterI_onepass(ctx , ctx->h , ctx->w , ctx->dst_1pass , ctx->h , dst_plane , dst_pitch);
	
}

void IIR_guass_filterI_destroy(long  h_flt)
{
    IIR_guass_fltI_context *ctx; 

    ctx = (IIR_guass_fltI_context *)h_flt;

    if (ctx->dst_1pass){
        free(ctx->dst_1pass);
        ctx->dst_1pass = 0;
    }

    if (ctx->w1){
        free(ctx->w1);
        ctx->w1 = 0;
    }

    if (ctx->w2){
        free(ctx->w2);
        ctx->w2 = 0;
    }

    free(ctx);
}


 

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

walletiger

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

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

余额充值