精华内容
下载资源
问答
  • 图像复原最小二乘方滤波
  • 数字图像处理作业图像运动模糊&约束最小二乘方滤波MATLAB源码及实验报告
  • 完整代码,可直接运行
  • 约束最小二乘方滤波

    2021-04-24 19:02:15
    最小二乘自适应滤波引 言基于最小均方误差(MMSE)准则的算法, 如最陡下降法、...(002)014 【摘要】介绍了图像退化模型和约束最小二乘滤波器以及平滑约束最小二乘滤 波器,并用 MATLAB7.0 实现约束最小二乘滤波恢复图像...

    最小二乘自适应滤波引 言基于最小均方误差(MMSE)准则的算法, 如最陡下降法、...

    (002)014 【摘要】介绍了图像退化模型和约束最小二乘滤波器以及平滑约束最小二乘滤 波器,并用 MATLAB7.0 实现约束最小二乘滤波恢复图像和平滑约束最小二乘滤......

    最小二乘递推算法和Kalman滤波算法_数学_自然科学_专业资料。第 15 卷第...

    ? (n), P0? , M ?1 ( n)d ( n) (2)前向预测误差滤波器 f M (n) 最小二乘前向预测器是用 n 时刻以前相继的 M 个数据, 对该时刻的 x(n) ......

    第29 卷 第4期 指挥控制与仿真 Command Control & Simulation Vol.29 No.4 2007 年 8 月 文章编号 1673-3819(2007)04-0041-02 Aug.2007 最小二乘滤波在......

    乘法王祖荫地 矿部 岩矿 测试 技术 研究所口 , 北京 有口 摘要 在卡 尔 曼滤 波的应用 中若 将时间因素固定则卡尔曼滤波退化为经典最小二乘法在通 , ,......

    基于时频域滤波及频域广义整体最小二乘辨识的飞机颤振模态参数辨识 [J], 唐炜; 史忠科; 李洪超 2.具有约束条件的最小二乘辨识及其应用 [J], 刘国隆 3.模态......

    针对强干扰环境下,基于线性约束最小方差(LCMV)准则,提出了一种极化域新的自适应滤波算法,采用变极化接收技术,实现对信号的最佳接收,仿真结果也证明了该方法的有效......

    基于MATLAB 的卡尔曼滤波与最小二乘滤波仿真实验设计一、 实验原理: 卡尔曼滤波器是一个最优化自回归数据处理算法, 对于解决很大部分的问题, 他是最优、效率最高......

    正则化参数的选择 4.5 误差平方加权和的更新递归 河南工业大学信息科学与工程学院 在本章中,我们将推广最小二乘的应用,以便推出一种设 计自适应横向滤波器的递归......

    最小二乘算法 赖晓平 【摘要】考虑二维线性相位矩形对称 FIR 滤波器的约束最小二乘设计问题,即在通 带和阻带逼近误差不超过给定值的约束下使逼近误差平方和最小......

    一个投影最小二乘算法,它是一个交替地更新有效约束集及将二次误差无约束极小点(最小二乘解)逐次投影到有效约束边界的迭代过程.通过二维FIR低通圆形滤波器和方 ......

    由已知的 x(n)来估计d (n) , 这时, 横向滤波器的输出是 d(n) 的最小二乘 估计 d?(n) , 即滤波方程为 d?(n) X0,M 1(n)wM (n) 其中, 采用......

    一 、 传统教学方 法与使用g IA 进行仿真实验的 比较 AIB. 统计信号处理课程 中的卡尔曼滤波 以及最d -乘滤波理论 , 在信 息处 理中有极其广泛的应用 ,......

    1. 用矢量空间法描述FTF算法中的4个横向滤波器 (1)最小二乘横向滤波器 wM (n) 设一 M阶横向滤波器的权矢量为 ? ? w ( n ) ? w ( n ), w ( ......

    最小二乘滤波在某型火控系统信息处理中的应用 作者:李金刚;赵勇;郭戈;谷敏 作者...

    一种基于最小二乘和多重渐消因子的改进 kalman 滤波方法 (57)摘要 本发明公开了一种基于最小二乘和多重渐 消因子的改进 kalman 滤波方法,本发明在滤波开 始......

    小型微型计算机系统 JoumalofCllineseComputerSystems 2013年8月第8期V01.34No.82013 小波滤波的移动最小二乘图像变形方法 李洪安,康宝生,张雷 (西北大学信息科学......

    基于最小二乘法的激光雷达数据滤波方法 赵慧珍 程英蕾 屈亚运 【摘要】摘要 最小二乘算法在“等权”的条件下进行,不具有先验知识,不 适用地形复杂的区域。针对......

    目前采用 的较广泛的自适应算法主要有递归最小二乘(RLS) 算法、 Widrow ? Hoff 最小均方( LMS) 算法、 格型滤波器算1 ] 法和无限冲激响应( ) 算法等 [......

    展开全文
  • 约束最小二乘方滤波就是其中一种较好的方法。在维纳滤波那一篇讲过,维纳滤波要求未退化图像和噪声的功率谱必须是已知的,一般这两个功率谱很难估计,尽管用一个常数去估计功率谱比,然而并不老是一个合适的解。约束...

    图像复原,简单讲,就是恢复图像原本的面貌,但因为各类缘由如图像采集过程当中出现的偏差致使获得的数字图像不清晰,不是咱们人眼看到的实物场景那样,所以须要采起技术手段去除图像的不清晰。约束最小二乘方滤波就是其中一种较好的方法。在维纳滤波那一篇讲过,维纳滤波要求未退化图像和噪声的功率谱必须是已知的,一般这两个功率谱很难估计,尽管用一个常数去估计功率谱比,然而并不老是一个合适的解。约束最小二乘方滤波要求噪声的方差和均值,这些参数可经过给定的退化图像计算出来,这是约束最小二乘方滤波的一个重要优势。函数

    一个图像采集系统输入输出的关系能够表示为g(x,y)=H[f(x,y)]+η(x,y),表达为向量-矩阵形式为g=Hf+η,明确地以矩阵形式来表达问题能够简化复原技术的推导。约束最小二乘方滤波的核心是解决退化函数H对噪声的敏感性问题,而减小噪声敏感性问题的一种方法是以平滑度量的最佳复原为基础的,如图像的二阶导数即拉普拉斯变换。因而,咱们找到一个带约束条件的最小准则函数C,定义以下:spa

    cf35f588321ce1c4672f35e557ebc787.png       (1)3d

    其约束条件为code

    7e8c1be746548310cc45cbc96a18fb2b.png             (2)blog

    其中,

    782854f4511a975547046e25b532ceb0.png是欧几里德向量范数,

    fc7ae988581c41d28b1da4af.html是未退化图像的估计。

    fc7ae988581c41d28b1da4af.html是拉普拉斯算子。求函数C的最小值,便获得最好的平滑效果即最佳复原。怎样找到最小值呢,这里咱们采用拉格朗日乘数法,在频率域中,函数C可表示为C=

    fc7ae988581c41d28b1da4af.html,其中P是拉普拉斯算子的傅里叶变换。则频率域中,拉格朗日函数为ci

    L(

    fc7ae988581c41d28b1da4af.html,λ)=

    fc7ae988581c41d28b1da4af.html+

    eb252e4ec4e0dc2f8a19b4b187a9a0ed.png(3)it

    其中N是加性噪声η的傅里叶变换,式(3)对

    fc7ae988581c41d28b1da4af.html求导,获得

    fc7ae988581c41d28b1da4af.html的最小值表达式以下:io

    cffc981347da50d8b7cf43caa5ecfd5d.png      (4)模板

    其中,γ=1/λ,H*为H的复共轭,求导过程以下class

    882b4c8d81787cbed271a8664aabdf4d.png

    约束最小二乘滤波的MATLAB代码以下

    close all;

    clear all;

    clc;

    I = im2double(imread('C:\Program Files\MATLAB\R2013a\bin\work\图像复原\lena.bmp'));

    [hei,wid,~] = size(I);

    subplot(2,2,1),imshow(I);

    title('原图像');

    % 模拟运动模糊.

    LEN = 21;

    THETA = 11;

    PSF = fspecial('motion', LEN, THETA);%产生运动模糊算子,即点扩展函数

    blurred = imfilter(I, PSF, 'conv', 'circular');

    subplot(2,2,2), imshow(blurred); title('模糊图像');

    Pf = psf2otf(PSF,[hei,wid]);%退化函数的FFT

    % 添加加性噪声

    noise_mean = 0;

    noise_var = 0.00001;

    blurred_noisy = imnoise(blurred, 'gaussian',noise_mean, noise_var);

    subplot(2,2,3), imshow(blurred_noisy)

    title('带运动模糊和噪声图像')

    p = [0 -1 0;-1 4 -1;0 -1 0];%拉普拉斯模板

    P = psf2otf(p,[hei,wid]);

    gama = 0.001;

    If = fft2(blurred_noisy);

    numerator = conj(Pf);%conj函数,用于求一个复数的复共轭

    denominator = Pf.^2 + gama*(P.^2);

    deblurred2 = ifft2( numerator.*If./ denominator );%约束最小二乘方滤波在频率域中的表达式

    subplot(2,2,4), imshow(deblurred2)

    title('约束最小二乘方滤波后图像');

    运行效果以下

    a5185ce11c0609e65a193e49815904a1.png

    展开全文
  • 【图像修复】基于约束最小二乘方滤波实现图像复原matlab源码GUI.md
  • 标题约束最小二乘方滤波几何均值滤波 约束最小二乘方滤波 F^(u,v)=[H∗(u,v)∣H(u,v)∣2+γ∣P(u,v)∣2]G(u,v)(5.89)\hat{F}(u,v) = \bigg[\frac{H^*(u,v)}{|H(u,v)|^2 + \gamma |P(u,v)|^2} \bigg]G(u,v) \tag{5.89...

    约束最小二乘方滤波

    F ^ ( u , v ) = [ H ∗ ( u , v ) ∣ H ( u , v ) ∣ 2 + γ ∣ P ( u , v ) ∣ 2 ] G ( u , v ) (5.89) \hat{F}(u,v) = \bigg[\frac{H^*(u,v)}{|H(u,v)|^2 + \gamma |P(u,v)|^2} \bigg]G(u,v) \tag{5.89} F^(u,v)=[H(u,v)2+γP(u,v)2H(u,v)]G(u,v)(5.89)

    P ( u , v ) 为 函 数 p ( x , y ) = [ 0 − 1 0 − 1 4 − 1 0 − 1 0 ] P(u,v) 为函数 p(x, y) = \begin {bmatrix}0 & -1 & 0 \\ -1 & 4 & -1 \\ 0 & -1 & 0 \end {bmatrix} P(u,v)p(x,y)=010141010的傅里叶变换

    我们认为这个函数是一个拉普拉斯核,注意式中在 γ = 0 \gamma = 0 γ=0时会简化为逆滤波。

    函数 P ( u , v ) 和 H ( u , v ) P(u,v)和H(u,v) P(u,v)H(u,v)的大小必须相同。意味着 p ( x , y ) p(x,y) p(x,y)必须嵌入 M M M x N N N 零阵列的中心。为保持 p ( x , y ) p(x,y) p(x,y)的偶对称, M M M N N N必须是偶整数。如果由其获得 H H H的一幅已知退化图像不是偶数维的,则在计算 H H H之前需要酌情删除一行和/或一行,以便可以使用

    2021-21-16 update, Thanks to fans id is ‘weixin_42674476’。

    def get_puv(image):
        h, w = image.shape[:2]
        
        h_pad, w_pad = h - 3, w - 3
        p_xy = np.array([[0, -1, 0],
                     [-1, 4, -1],
                     [0, -1, 0]])
        p_pad = np.pad(p_xy, ((h_pad//2, h_pad - h_pad//2), (w_pad//2, w_pad - w_pad//2)), mode='constant')
        p_uv = np.fft.fft2(p_pad)
        return p_uv
    
    def least_square_filter(image, PSF, eps, gamma=0.01):
        """
        least square filter for image denoise, math: 
        $$\hat{F}(u,v) = \bigg[\frac{H^*(u,v)}{|H(u,v)|^2 + \gamma |P(u,v)|^2} \bigg]G(u,v)$$
        param: image: input image
        param: PSF: input the PSF mask
        param: eps: epsilon
        param: gamma: gamma value for least square filter fuction
        return image after least square filter
        问题:为什么改变gamma值结果也没有变化,sorted 2021-12-16
        """
        fft = np.fft.fft2(image)
        PSF_fft = np.fft.fft2(PSF)
        conj = PSF_fft.conj()
        p_uv = get_puv(image)
        abs_conj = np.abs(PSF_fft) ** 2
        # abs_conj = (PSF_fft * conj).real
        huv = conj / (abs_conj + gamma * (np.abs(p_uv) ** 2))
        result = np.fft.ifft2(fft * huv)
        result = np.abs(np.fft.fftshift(result))
        return result, abs_conj
    
    # 约束最小二乘方滤波
    image = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0526(a)(original_DIP).tif", 0)
    PSF = get_motion_dsf(image.shape[:2], -50, 100)
    blurred = make_blurred(image, PSF, 1e-5)
    
    wiener = wiener_filter(blurred, PSF, 1e-5, K=0.03)     
    least_square, abs_con = least_square_filter(blurred, PSF, 1e-5, gamma=1e-6)
    least_square = np.uint8(normalize(least_square) * 255)
    
    img_diff = image - least_square
    plt.figure(figsize=(15, 5))
    plt.subplot(131), plt.imshow(blurred, 'gray'), plt.title("Motion blurred"), plt.xticks([]), plt.yticks([])
    plt.subplot(132), plt.imshow(wiener, 'gray'), plt.title("Wiener Filter"), plt.xticks([]), plt.yticks([])
    plt.subplot(133), plt.imshow(least_square, 'gray'), plt.title("Least Square Filter"), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    几何均值滤波

    F ^ ( u , v ) = [ H ∗ ( u , v ) ∣ H ( u , v ) ∣ 2 ] α [ H ∗ ( u , v ) ∣ H ( u , v ) ∣ 2 + β [ S η ( u , v ) S f ( u , v ) ] ] 1 − α G ( u , v ) (5.99) \hat{F}(u,v) = \bigg[\frac{H^*(u,v)}{|H(u,v)|^2}\bigg]^\alpha \Bigg[\frac{H^*(u,v)}{|H(u,v)|^2 + \beta \big[\frac{S_{\eta}(u,v)}{S_f(u,v)} \big]}\Bigg]^{1-\alpha} G(u, v) \tag{5.99} F^(u,v)=[H(u,v)2H(u,v)]α[H(u,v)2+β[Sf(u,v)Sη(u,v)]H(u,v)]1αG(u,v)(5.99)

    α 和 β \alpha 和 \beta αβ 是非负的实常数。

    α = 1 \alpha=1 α=1时,几何均值滤波器简化为逆滤波器;当 α = 0 \alpha=0 α=0时,几何均值滤波器变为参数维纳滤波器,参数维纳滤波器在 β = 1 \beta=1 β=1时简化为标准维纳滤波器;当 α = 1 / 2 \alpha=1/2 α=1/2时几何均值滤波器变成幂次相同的两个量的乘积,这就是几何均值的定义。

    β = 1 \beta=1 β=1 α \alpha α大于1/2时,滤波器更的性能更像逆滤波器。
    α \alpha α 小于1/2时,滤波器的性能更像维纳滤波器。
    α \alpha α 等于1/2且 β = 1 \beta=1 β=1时,滤波器通常称为频谱均衡滤波器。

    def geometric_mean_filter(image, PSF, eps, K=1, alpha=1, beta=1):
        """
        geometric mean filter for image denoise, math: 
        $$\hat{F}(u,v) = \bigg[\frac{H^*(u,v)}{|H(u,v)|^2}\bigg]^\alpha 
        \Bigg[\frac{H^*(u,v)}{|H(u,v)|^2 + \beta \big[\frac{S_{\eta}(u,v)}{S_f(u,v)} \big]}\Bigg]^{1-\alpha} G(u, v) $$
        param: image: input image
        param: PSF  : input the PSF mask
        param: eps  : epsilon
        param: K    : K equal to math: \frac{S_{\eta}(u,v)}{S_f(u,v)}
        param: alpha: alpha value for filter fuction
        param: beta : beta value for the filter fuction
        return image after least square filter
        """
        fft = np.fft.fft2(image)
        PSF_fft = np.fft.fft2(PSF)
        conj = PSF_fft.conj()
        abs_square = (PSF_fft * conj).real
        huv = np.power(conj / (abs_square), alpha) * np.power(conj / (abs_square + beta*(K)), 1 - alpha)
        result = np.fft.ifft2(fft * huv)
        result = np.abs(np.fft.fftshift(result))
        return result
    
    #### 几何均值滤波器
    image = cv2.imread("DIP_Figures/DIP3E_Original_Images_CH05/Fig0526(a)(original_DIP).tif", 0)
    PSF = get_motion_dsf(image.shape[:2], -50, 100)
    blurred = make_blurred(image, PSF, 1e-5)
    
    wiener = wiener_filter(blurred, PSF, 1e-5, K=0.03)     
    geometric_mean = geometric_mean_filter(blurred, PSF, 1e-5, K=1, alpha=1/2, beta=0)
    geometric_mean = np.uint8(normalize(geometric_mean) * 255)
    
    plt.figure(figsize=(14, 5))
    plt.subplot(131), plt.imshow(blurred, 'gray'), plt.title("Motion blurred"), plt.xticks([]), plt.yticks([])
    plt.subplot(132), plt.imshow(wiener, 'gray'), plt.title("Wiener Filter"), plt.xticks([]), plt.yticks([])
    plt.subplot(133), plt.imshow(geometric_mean, 'gray'), plt.title("Geometric mean Filter"), plt.xticks([]), plt.yticks([])
    plt.tight_layout()
    plt.show()
    

    在这里插入图片描述

    展开全文
  • 本节为opencv数字图像处理(10):图像复原与重建的第三小节,逆滤波、维纳滤波、约束最小二乘方滤波和几何均值滤波,主要包括:四种滤波复原图像的数学推导以及维纳滤波的C++实现。

    本节为opencv数字图像处理(12):图像复原与重建的第三小节,逆滤波、维纳滤波、约束最小二乘方滤波和几何均值滤波,主要包括:四种滤波复原图像的数学推导以及维纳滤波的C++实现。

    1. 逆滤波器

    emsp; 若退化函数已知或可以得到一个估计,最简单的图像复原方法就是直接做逆滤波,用退化函数除退化图像的傅立叶变换来计算原始图形的傅立叶变换的估计即:
    在这里插入图片描述
      展开计算为:
    在这里插入图片描述
      如果退化噪声为0或很小,噪声就会支配估计值 F ^ ( u , v ) \hat F(u,v) F^(u,v),这时候经常需要限制滤波的频率,使其接近原点。

    1. 最小均方误差滤波/维纳滤波

      这种滤波方法建立在图像和噪声都是随机变量的基础上,目标是找出未污染图像 f ( x , y ) f(x,y) f(x,y)的一个估计 f ^ \hat f f^,使它们之间的均方误差最小,误差度量由下式给出:
    在这里插入图片描述
      假设噪声和图像不相关,二者中有一个有零均值且估计中的灰度级是退化图像中灰度级的线性函数,则误差函数的最小值在频率与中由下式给出:
    在这里插入图片描述
      这个结果就是维纳滤波,方括号中的项组成的滤波器称为最小均方差滤波器或最小二乘误差滤波器。式子中: H ( u , v ) = H(u,v)= H(u,v)=退化函数; H ∗ ( u , v ) H*(u,v) H(u,v)= H ( u , v ) H(u,v) H(u,v)的复共轭; ∣ H ( u , v ) ∣ 2 = H ∗ ( u , v ) H ( u , v ) |H(u,v)|^2=H*(u,v)H(u,v) H(u,v)2=H(u,v)H(u,v) S e t a ( u , v ) = ∣ N ( u , v ) ∣ 2 = S_{eta}(u,v)=|N(u,v)|^2= Seta(u,v)=N(u,v)2=噪声的功率谱。

      如果噪声为0,那么噪声功率谱消失,维纳滤波简化为逆滤波。

      信噪比是一个基于噪声和未退化图像的功率谱为基础的,频率域可用下式近似:
    在这里插入图片描述
      携带低噪声的图像SNR较高,是表征复原算法的性能的一个重要量度。

      均方误差也可以这样描述:
    在这里插入图片描述
      而如果把复原图像考虑为“信号”,复原图像和原图像的差考虑为噪声,则空间域中信噪比可定义如下:
    在这里插入图片描述
       f f f f ^ \hat f f^越接近,这个比值越大。【这个定量与感觉上的图像质量没有很好的必然关系】

      当处理白噪声时 ∣ N ( u , v ) ∣ 2 |N(u,v)|^2 N(u,v)2是一个常数,但是未退化图像和噪声的功率谱通常未知或不能估计,则可用下式近似:
    在这里插入图片描述
      C++实现来自这篇博文

    #include <iostream>
    #include "opencv2/imgproc.hpp"
    #include "opencv2/imgcodecs.hpp"
    #include <opencv2/opencv.hpp>
    
    using namespace cv;
    using namespace std;
    
    void calcPSF(Mat& outputImg, Size filterSize, int len, double theta);
    void fftshift(const Mat& inputImg, Mat& outputImg);
    void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H);
    void calcWnrFilter(const Mat& input_h_PSF, Mat& output_G, double nsr);
    void edgetaper(const Mat& inputImg, Mat& outputImg, double gamma = 5.0, double beta = 0.2);
    
    int LEN = 50;
    int THETA = 360;
    int snr = 8000;
    Mat imgIn;
    Rect roi;
    static void onChange(int pos, void* userInput);
    
    int main(int argc, char* argv[])
    {
    	string strInFileName = "1.JPG";
    
    	imgIn = imread(strInFileName, IMREAD_GRAYSCALE);
    	if (imgIn.empty()) //check whether the image is loaded or not
    	{
    		cout << "ERROR : Image cannot be loaded..!!" << endl;
    		return -1;
    	}
    	imshow("src", imgIn);
    
    	// it needs to process even image only
    	roi = Rect(0, 0, imgIn.cols & -2, imgIn.rows & -2);
    	imgIn = imgIn(roi);
    	cv::namedWindow("inverse");
    
    	createTrackbar("LEN", "inverse", &LEN, 200, onChange, &imgIn);
    	onChange(0, 0);
    	createTrackbar("THETA", "inverse", &THETA, 360, onChange, &imgIn);
    	onChange(0, 0);
    	createTrackbar("snr", "inverse", &snr, 10000, onChange, &imgIn);
    	onChange(0, 0);
    	imshow("inverse", imgIn);
    	cv::waitKey(0);
    
    	return 0;
    }
    
    void calcPSF(Mat& outputImg, Size filterSize, int len, double theta)
    {
    	Mat h(filterSize, CV_32F, Scalar(0));
    	Point point(filterSize.width / 2, filterSize.height / 2);
    	ellipse(h, point, Size(0, cvRound(float(len) / 2.0)), 90.0 - theta,
    		0, 360, Scalar(255), FILLED);
    	Scalar summa = sum(h);
    	outputImg = h / summa[0];
    	Mat tmp;
    	normalize(outputImg, tmp, 1, 0, NORM_MINMAX);
    	imshow("psf", tmp);
    }
    void fftshift(const Mat& inputImg, Mat& outputImg)
    {
    	outputImg = inputImg.clone();
    	int cx = outputImg.cols / 2;
    	int cy = outputImg.rows / 2;
    	Mat q0(outputImg, Rect(0, 0, cx, cy));
    	Mat q1(outputImg, Rect(cx, 0, cx, cy));
    	Mat q2(outputImg, Rect(0, cy, cx, cy));
    	Mat q3(outputImg, Rect(cx, cy, cx, cy));
    	Mat tmp;
    	q0.copyTo(tmp);
    	q3.copyTo(q0);
    	tmp.copyTo(q3);
    	q1.copyTo(tmp);
    	q2.copyTo(q1);
    	tmp.copyTo(q2);
    }
    void filter2DFreq(const Mat& inputImg, Mat& outputImg, const Mat& H)
    {
    	Mat planes[2] = { Mat_<float>(inputImg.clone()), Mat::zeros(inputImg.size(), CV_32F) };
    	Mat complexI;
    	merge(planes, 2, complexI);
    	dft(complexI, complexI, DFT_SCALE);
    	Mat planesH[2] = { Mat_<float>(H.clone()), Mat::zeros(H.size(), CV_32F) };
    	Mat complexH;
    	merge(planesH, 2, complexH);
    	Mat complexIH;
    	mulSpectrums(complexI, complexH, complexIH, 0);
    	idft(complexIH, complexIH);
    	split(complexIH, planes);
    	outputImg = planes[0];
    }
    void calcWnrFilter(const Mat& input_h_PSF, Mat& output_G, double nsr)
    {
    	Mat h_PSF_shifted;
    	fftshift(input_h_PSF, h_PSF_shifted);
    	Mat planes[2] = { Mat_<float>(h_PSF_shifted.clone()), Mat::zeros(h_PSF_shifted.size(), CV_32F) };
    	Mat complexI;
    	merge(planes, 2, complexI);
    	dft(complexI, complexI);
    	split(complexI, planes);
    	Mat denom;
    	pow(abs(planes[0]), 2, denom);
    	denom += nsr;
    	divide(planes[0], denom, output_G);
    }
    void edgetaper(const Mat& inputImg, Mat& outputImg, double gamma, double beta)
    {
    	int Nx = inputImg.cols;
    	int Ny = inputImg.rows;
    	Mat w1(1, Nx, CV_32F, Scalar(0));
    	Mat w2(Ny, 1, CV_32F, Scalar(0));
    	float* p1 = w1.ptr<float>(0);
    	float* p2 = w2.ptr<float>(0);
    	float dx = float(2.0 * CV_PI / Nx);
    	float x = float(-CV_PI);
    	for (int i = 0; i < Nx; i++)
    	{
    		p1[i] = float(0.5 * (tanh((x + gamma / 2) / beta) - tanh((x - gamma / 2) / beta)));
    		x += dx;
    	}
    	float dy = float(2.0 * CV_PI / Ny);
    	float y = float(-CV_PI);
    	for (int i = 0; i < Ny; i++)
    	{
    		p2[i] = float(0.5 * (tanh((y + gamma / 2) / beta) - tanh((y - gamma / 2) / beta)));
    		y += dy;
    	}
    	Mat w = w2 * w1;
    	multiply(inputImg, w, outputImg);
    }
    
    // Trackbar call back function
    static void onChange(int, void* userInput)
    {
    	Mat imgOut;
    	//Hw calculation (start)
    	Mat Hw, h;
    	calcPSF(h, roi.size(), LEN, (double)THETA);
    	calcWnrFilter(h, Hw, 1.0 / double(snr));
    	//Hw calculation (stop)
    	imgIn.convertTo(imgIn, CV_32F);
    	edgetaper(imgIn, imgIn);
    	// filtering (start)
    	filter2DFreq(imgIn(roi), imgOut, Hw);
    	// filtering (stop)
    	imgOut.convertTo(imgOut, CV_8U);
    	normalize(imgOut, imgOut, 0, 255, NORM_MINMAX);
    	//    imwrite("result.jpg", imgOut);
    	imshow("inverse", imgOut);
    }
    

    2. 约束最小二乘方滤波

      虽然维纳滤波中可以根据上式来进行估计,但是很少得到合适的解,我们将图像受噪声污染表达为如下的矩阵形式:
    g = H f + η g=Hf+\eta g=Hf+η
      考虑一个带约束的最小准则函数 C C C,定义如下:
    在这里插入图片描述
      约束为:
    在这里插入图片描述
      在频率域这个最佳化问题的解决方案由如下公式给出:
    在这里插入图片描述
       γ \gamma γ是一个参数,调整以满足上2式, P ( u , v ) P(u,v) P(u,v)是函数:
    在这里插入图片描述
      的傅立叶变换,当 γ = 0 \gamma=0 γ=0时最小二乘方滤波退化为直接逆滤波。

       γ \gamma γ可以交互地手动调整,当然也可以迭代计算,步骤如下。

      定义一个残差向量 r r r为:
    在这里插入图片描述
      由上3式知, F ^ ( u , v ) \hat F(u,v) F^(u,v) γ \gamma γ的函数,所以 r r r是该参数的函数,可以证明:
    在这里插入图片描述
      是 γ \gamma γ的单调递增函数,我们需要调整 γ \gamma γ使得:
    在这里插入图片描述
      若 ∣ ∣ r ∣ ∣ 2 = ∣ ∣ η ∣ ∣ 2 ||r||^2=||\eta||^2 r2=η2,那么可以严格满足上述约束,具体地,寻找一个 γ \gamma γ值得方法:

    • 指定 γ \gamma γ的初始值
    • 计算 ∣ ∣ r ∣ ∣ 2 ||r||^2 r2
    • 满足上1式则停止;否则,若 ∣ ∣ r ∣ ∣ 2 < ∣ ∣ η ∣ ∣ 2 − a ||r||^2<||\eta||^2-a r2<η2a,增大 γ \gamma γ,若 ∣ ∣ r ∣ ∣ 2 = ∣ ∣ η ∣ ∣ 2 + a ||r||^2=||\eta||^2+a r2=η2+a,则减少 γ \gamma γ。然后返回步骤2。重新计算最佳估计 F ^ ( u , v ) \hat F(u,v) F^(u,v)
        上面步骤中,为了计算 ∣ ∣ r ∣ ∣ 2 ||r||^2 r2,由上式3【从此处向上数第三个列出的公式】有:
      在这里插入图片描述
        计算 R ( u , v ) R(u,v) R(u,v)的傅里叶反变换得 r ( x , y ) r(x,y) r(x,y),有:
      在这里插入图片描述
        计算 ∣ ∣ η ∣ ∣ 2 ||\eta||^2 η2的话,首先考虑正负图像的噪声方差:
      在这里插入图片描述
        其中:
      在这里插入图片描述
        为样本均值。注意到上2式【从此处向上数第二个列出的公式】的双重求和即为 ∣ ∣ η ∣ ∣ 2 ||\eta||^2 η2,给出如下表达式:
      在这里插入图片描述
        也就是说,仅仅使用噪声的均值和方差的知识,就可以实现一个最佳复原算法,但必须假设噪声和图像灰度值不相关。但是需要指出,约束最小二乘方意义小的最佳复原并不是视觉效果上最好。

    3. 几何均值滤波

      几何均值滤波是维纳滤波的推广:
    在这里插入图片描述
      当 α = 1 \alpha=1 α=1,退化为直接逆滤波,当 α = 0 \alpha=0 α=0,参数维纳滤波器,参数维纳滤波器在 β = 1 \beta=1 β=1时还原为标准的维纳滤波器。如果 α = 1 / 2 \alpha=1/2 α=1/2,则滤波器变成相同幂次的两个量的积,也就是几何均值的定义。当 β = 1 \beta=1 β=1,随着 α \alpha α减小到 1 / 2 1/2 1/2以下,滤波器性能越来越接近逆滤波器,增大到 1 / 2 1/2 1/2以上,越来越接近维纳滤波器。当 α = 1 / 2 \alpha=1/2 α=1/2 β = 1 \beta=1 β=1时,称为谱均衡滤波器。


    欢迎扫描二维码关注微信公众号 深度学习与数学   [每天获取免费的大数据、AI等相关的学习资源、经典和最新的深度学习相关的论文研读,算法和其他互联网技能的学习,概率论、线性代数等高等数学知识的回顾]
    在这里插入图片描述

    展开全文
  • 约束最小二乘方滤波去模糊

    千次阅读 2017-12-01 17:10:18
    约束最小二乘滤波仅要求噪声方差和均值的知识。 g=Hf+ηg=Hf+\eta 假设g(x,y)的大小为M×NM\times N,g(x,y)g(x,y)第一行的图像元素构成向量gg的第一组NN个元素,第行构成下一组NN个元素,结果向量MN×1MN\times ...
  • 图像去模糊(约束最小二乘方滤波

    万次阅读 多人点赞 2015-08-09 15:44:57
    约束最小二乘方滤波(Constrained Least Squares Filtering,aka Tikhonov filtration,Tikhonov regularization)核心是H对噪声的敏感性问题。减少噪声敏感新问题的一种方法是以平滑度量的最佳复原为基础的,因此...
  • 约束最小二乘方滤波约束最小二乘方图像复原公式推导算法实现步骤matlab 代码实现实验结果总结 约束最小二乘方图像复原 维纳滤波图像复原技术是指除了要求了解关于退化系统的传递函数之外,还需要知道未退化...
  • 文章目录一、维纳滤波二约束最小乘方滤波 一、维纳滤波 对于运动引起的图像模糊,最简单的方法是直接做逆滤波,但是逆滤波对加性噪声特别敏感,使得恢复的图像几乎不可用。最小均方差(维纳)滤波用来去除含有...
  • 约束最小二乘方滤波matlab,具体细节在:https://blog.csdn.net/MARSHCW/article/details/109614832
  • 约束最小二乘方滤波(Constrained Least Squares Filtering,aka Tikhonov filtration,Tikhonov regularization)核心是H对噪声的敏感性问题。减少噪声敏感新问题的一种方法是以平滑度量的最佳复原为基础的,因此我们...
  • 维纳滤波约束最小二乘滤波图像复原自编matlab代码,共有两个文件CLSFilter.m,WienerFilter.m和一张测试图,可直接在R2013b上可以运行,有详细注释,注释里还有参考资料的网页链接,可帮助理解代码。
  • 约束最小二乘方滤波(Constrained Least Squares Filtering,aka Tikhonov filtration,Tikhonov regularization)核心是H对噪声的敏感性问题。减少噪声敏感新问题的一种方法是以平滑度量的最佳复原为基础的,因此...
  • 约束最小二乘滤波复原(a)原始图像 (b)运动造成的图像模糊 (c)Wiener滤波去图像模糊 维纳滤波复原实例 I = imread('lena.jpg'); figure(1),imshow(I); LEN = 21; THETA = 11; PSF = fspecial('motion', LEN, THETA); ...

空空如也

空空如也

1 2 3 4 5 ... 11
收藏数 203
精华内容 81
关键字:

约束最小二乘方滤波