文章目录

  • 图像恢复的一般运算过程
  • 什么是「最小均方差滤波」
  • 实现步骤
  • 实现代码
  • 最后的结果

图像恢复的一般运算过程

我们从前几章的基本理论出发，退化信号恢复成原始信号的步骤，可以概括成两步基本公式。对于空间来说，有如下公式

{ g ( x , y ) = h ( x , y ) ⊙ f ( x , y ) + η ( x , y ) f ^ ( x , y ) = ω ( x , y ) ∗ g ( x , y ) − η ( x , y ) \left \{ \begin{matrix} g(x, y) = h(x, y) \odot f(x, y) + \eta(x, y) \\ \hat f(x, y) = \omega(x, y) * g(x, y) - \eta(x, y) \end{matrix} \right. {g(x,y)=h(x,y)⊙f(x,y)+η(x,y)f^​(x,y)=ω(x,y)∗g(x,y)−η(x,y)​

退化的过程，是由退化函数 h ( x , y ) h(x, y) h(x,y) 与原图像 f ( x , y ) f(x, y) f(x,y) 通过卷积，然后加上噪声函数 η ( x , y ) \eta(x, y) η(x,y) 最终得到退化图像 g ( x , y ) g(x, y) g(x,y)；而恢复方法则是通过评估或建模恢复函数 ω ( x , y ) \omega (x, y) ω(x,y) ，与退化后的图像 g ( x , y ) g(x, y) g(x,y) 做运算后得到。

由于针对不同退化效果的还原运算方式的不同，这决定了 ω ( x , y ) \omega(x, y) ω(x,y) 不一定是 h ( x , y ) h(x, y) h(x,y) 的反运算，而且很难有通用的技巧可言。相比之下，图像恢复过程在频域上则简单许多，它的还原过程被表述如下：

{ G ( u , v ) = H ( u , v ) ∗ F ( u , v ) + D ( u , v ) F ^ ( u , v ) = G ( u , v ) ∗ H ( u , v ) − 1 − D ( u , v ) \left \{ \begin{matrix} G(u, v) = H(u, v) * F(u, v) + D(u, v) \\ \hat F(u, v) = G(u, v) * H(u, v)^{-1} - D(u, v) \end{matrix} \right. {G(u,v)=H(u,v)∗F(u,v)+D(u,v)F^(u,v)=G(u,v)∗H(u,v)−1−D(u,v)​

希腊字母 Eta 的大写是 H \Eta H，为了跟已有的 H 做区别，这里改成了 D。可以看到，在频域下，只要得到退化函数的反函数，就可以直接得到还原图像 F ^ \hat F F^。另外稍微说一下，无论对于空间还是频域来说，评估图像的噪音函数并不困难，所以以上过程又经常被人省略了噪音的影响。

事实上，在做图像退化与恢复的研究时，经常会省略噪音函数的影响。

在上一章 《数字图像学笔记——15. 图像退化与复原（使用逆滤波复原图像）》 通过逆滤波函数已经演示了为什么在知道 H ( u , v ) H(u, v) H(u,v) 的情况下，可以直接还原出图像 F ^ ( u , v ) \hat F(u, v) F^(u,v)。

现在我们需要来考虑一个更一般的问题，如果我们不知道导致图像退化的函数 H ( u , v ) H(u, v) H(u,v) ，我们又应该如何还原图像呢。在这样一个假设前提下，有人提出了一类被称为的自适应的滤波方法，而「最小均方差滤波」就是诸多自适应滤波的其中一种。

什么是「最小均方差滤波」

最小均方滤波器（Least Mean Square Filter，或LMS Filter）是一类可通过最小化误差信号（error signal）的均方值（mean square）而修正滤波器系数，以模拟所需理想滤波器的自适应滤波器1。

在冈萨雷斯版《数字图像处理》关于这一节中，老爷子并没有给出详细的证明或者推导过程，这对于第一次接触这个概念的朋友显然造成了困扰，而且这个问题显然在他最新版的著作里也没有得到修复。

所以，作为这本书的学习笔记，我只好先从信号学的角度出发，解释什么是最小均方差滤波。

补充信息
关于自适应算法，能够考证到的文献应该来自于1959年，Widrow等人的论文「Adaptive pattern recognition system」。
有需求的朋友可以顺着这个线索自行查阅相关论文。

从信号学的角度看，「最小均方差」算法属于一种被称为「梯度下降」技术的应用扩展。

补充信息
关于什么是梯度下降，有需要的朋友可以参考我这几篇文章内容。
深度学习知识总结—— 1.1. 什么是梯度下降
深度学习知识总结—— 1.2.梯度下降算法实现

它的评价指标，就是用均方差评价恢复的信号 F ^ ( u , v ) \hat F(u, v) F^(u,v) 与原始信号 F ( u , v ) F(u, v) F(u,v) 的离散程度。目标是找到一个能接受，趋近于 0 的误差。

所以，一般均方差公式尽管写为

M S E = 1 n ∑ i = 1 K ( y i − y ˉ ) 2 MSE = \frac{1}{n} \sum_{i=1}^{K} (y_i - \bar y)^2 MSE=n1​i=1∑K​(yi​−yˉ​)2

由于对于评价指标来说 1 / n 1/n 1/n 其实没什么必要性，所以你在信号学或者这本《数字图像处理》里，看到的就是下面这个样子

E { e 2 ( n ) } = { F ^ ( n ) − F ( n ) } 2 E\{ e^2(n) \} = \{\hat F(n) - F(n) \} ^2 E{e2(n)}={F^(n)−F(n)}2

和之前的内容稍微有点出入的是，我们这里需要考虑算法迭代 n n n。什么意思呢，在第一次运算后，我们的算法会给出一个 H ( u , v ) 1 − 1 H(u, v)_1^{-1} H(u,v)1−1​ ，我们于是可以得出 F ^ ( 1 ) \hat F(1) F^(1)，它会与原始图像计算得到一个 E { e 2 ( 1 ) } E\{ e^2(1) \} E{e2(1)}；然后这个步骤会重复几次，最终得到一个 E { e 2 ( n ) } E\{ e^2(n) \} E{e2(n)}，而这时的误差应该最小，得到我们理想的结果。

#mermaid-svg-pGQeQ4y0vVFPtLFp .label{font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family);fill:#333;color:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .label text{fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .node rect,#mermaid-svg-pGQeQ4y0vVFPtLFp .node circle,#mermaid-svg-pGQeQ4y0vVFPtLFp .node ellipse,#mermaid-svg-pGQeQ4y0vVFPtLFp .node polygon,#mermaid-svg-pGQeQ4y0vVFPtLFp .node path{fill:#ECECFF;stroke:#9370db;stroke-width:1px}#mermaid-svg-pGQeQ4y0vVFPtLFp .node .label{text-align:center;fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .node.clickable{cursor:pointer}#mermaid-svg-pGQeQ4y0vVFPtLFp .arrowheadPath{fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .edgePath .path{stroke:#333;stroke-width:1.5px}#mermaid-svg-pGQeQ4y0vVFPtLFp .flowchart-link{stroke:#333;fill:none}#mermaid-svg-pGQeQ4y0vVFPtLFp .edgeLabel{background-color:#e8e8e8;text-align:center}#mermaid-svg-pGQeQ4y0vVFPtLFp .edgeLabel rect{opacity:0.9}#mermaid-svg-pGQeQ4y0vVFPtLFp .edgeLabel span{color:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .cluster rect{fill:#ffffde;stroke:#aa3;stroke-width:1px}#mermaid-svg-pGQeQ4y0vVFPtLFp .cluster text{fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp div.mermaidTooltip{position:absolute;text-align:center;max-width:200px;padding:2px;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family);font-size:12px;background:#ffffde;border:1px solid #aa3;border-radius:2px;pointer-events:none;z-index:100}#mermaid-svg-pGQeQ4y0vVFPtLFp .actor{stroke:#ccf;fill:#ECECFF}#mermaid-svg-pGQeQ4y0vVFPtLFp text.actor>tspan{fill:#000;stroke:none}#mermaid-svg-pGQeQ4y0vVFPtLFp .actor-line{stroke:grey}#mermaid-svg-pGQeQ4y0vVFPtLFp .messageLine0{stroke-width:1.5;stroke-dasharray:none;stroke:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .messageLine1{stroke-width:1.5;stroke-dasharray:2, 2;stroke:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp #arrowhead path{fill:#333;stroke:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .sequenceNumber{fill:#fff}#mermaid-svg-pGQeQ4y0vVFPtLFp #sequencenumber{fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp #crosshead path{fill:#333;stroke:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .messageText{fill:#333;stroke:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .labelBox{stroke:#ccf;fill:#ECECFF}#mermaid-svg-pGQeQ4y0vVFPtLFp .labelText,#mermaid-svg-pGQeQ4y0vVFPtLFp .labelText>tspan{fill:#000;stroke:none}#mermaid-svg-pGQeQ4y0vVFPtLFp .loopText,#mermaid-svg-pGQeQ4y0vVFPtLFp .loopText>tspan{fill:#000;stroke:none}#mermaid-svg-pGQeQ4y0vVFPtLFp .loopLine{stroke-width:2px;stroke-dasharray:2, 2;stroke:#ccf;fill:#ccf}#mermaid-svg-pGQeQ4y0vVFPtLFp .note{stroke:#aa3;fill:#fff5ad}#mermaid-svg-pGQeQ4y0vVFPtLFp .noteText,#mermaid-svg-pGQeQ4y0vVFPtLFp .noteText>tspan{fill:#000;stroke:none}#mermaid-svg-pGQeQ4y0vVFPtLFp .activation0{fill:#f4f4f4;stroke:#666}#mermaid-svg-pGQeQ4y0vVFPtLFp .activation1{fill:#f4f4f4;stroke:#666}#mermaid-svg-pGQeQ4y0vVFPtLFp .activation2{fill:#f4f4f4;stroke:#666}#mermaid-svg-pGQeQ4y0vVFPtLFp .mermaid-main-font{font-family:"trebuchet ms", verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-pGQeQ4y0vVFPtLFp .section{stroke:none;opacity:0.2}#mermaid-svg-pGQeQ4y0vVFPtLFp .section0{fill:rgba(102,102,255,0.49)}#mermaid-svg-pGQeQ4y0vVFPtLFp .section2{fill:#fff400}#mermaid-svg-pGQeQ4y0vVFPtLFp .section1,#mermaid-svg-pGQeQ4y0vVFPtLFp .section3{fill:#fff;opacity:0.2}#mermaid-svg-pGQeQ4y0vVFPtLFp .sectionTitle0{fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .sectionTitle1{fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .sectionTitle2{fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .sectionTitle3{fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .sectionTitle{text-anchor:start;font-size:11px;text-height:14px;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-pGQeQ4y0vVFPtLFp .grid .tick{stroke:#d3d3d3;opacity:0.8;shape-rendering:crispEdges}#mermaid-svg-pGQeQ4y0vVFPtLFp .grid .tick text{font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-pGQeQ4y0vVFPtLFp .grid path{stroke-width:0}#mermaid-svg-pGQeQ4y0vVFPtLFp .today{fill:none;stroke:red;stroke-width:2px}#mermaid-svg-pGQeQ4y0vVFPtLFp .task{stroke-width:2}#mermaid-svg-pGQeQ4y0vVFPtLFp .taskText{text-anchor:middle;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-pGQeQ4y0vVFPtLFp .taskText:not([font-size]){font-size:11px}#mermaid-svg-pGQeQ4y0vVFPtLFp .taskTextOutsideRight{fill:#000;text-anchor:start;font-size:11px;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-pGQeQ4y0vVFPtLFp .taskTextOutsideLeft{fill:#000;text-anchor:end;font-size:11px}#mermaid-svg-pGQeQ4y0vVFPtLFp .task.clickable{cursor:pointer}#mermaid-svg-pGQeQ4y0vVFPtLFp .taskText.clickable{cursor:pointer;fill:#003163 !important;font-weight:bold}#mermaid-svg-pGQeQ4y0vVFPtLFp .taskTextOutsideLeft.clickable{cursor:pointer;fill:#003163 !important;font-weight:bold}#mermaid-svg-pGQeQ4y0vVFPtLFp .taskTextOutsideRight.clickable{cursor:pointer;fill:#003163 !important;font-weight:bold}#mermaid-svg-pGQeQ4y0vVFPtLFp .taskText0,#mermaid-svg-pGQeQ4y0vVFPtLFp .taskText1,#mermaid-svg-pGQeQ4y0vVFPtLFp .taskText2,#mermaid-svg-pGQeQ4y0vVFPtLFp .taskText3{fill:#fff}#mermaid-svg-pGQeQ4y0vVFPtLFp .task0,#mermaid-svg-pGQeQ4y0vVFPtLFp .task1,#mermaid-svg-pGQeQ4y0vVFPtLFp .task2,#mermaid-svg-pGQeQ4y0vVFPtLFp .task3{fill:#8a90dd;stroke:#534fbc}#mermaid-svg-pGQeQ4y0vVFPtLFp .taskTextOutside0,#mermaid-svg-pGQeQ4y0vVFPtLFp .taskTextOutside2{fill:#000}#mermaid-svg-pGQeQ4y0vVFPtLFp .taskTextOutside1,#mermaid-svg-pGQeQ4y0vVFPtLFp .taskTextOutside3{fill:#000}#mermaid-svg-pGQeQ4y0vVFPtLFp .active0,#mermaid-svg-pGQeQ4y0vVFPtLFp .active1,#mermaid-svg-pGQeQ4y0vVFPtLFp .active2,#mermaid-svg-pGQeQ4y0vVFPtLFp .active3{fill:#bfc7ff;stroke:#534fbc}#mermaid-svg-pGQeQ4y0vVFPtLFp .activeText0,#mermaid-svg-pGQeQ4y0vVFPtLFp .activeText1,#mermaid-svg-pGQeQ4y0vVFPtLFp .activeText2,#mermaid-svg-pGQeQ4y0vVFPtLFp .activeText3{fill:#000 !important}#mermaid-svg-pGQeQ4y0vVFPtLFp .done0,#mermaid-svg-pGQeQ4y0vVFPtLFp .done1,#mermaid-svg-pGQeQ4y0vVFPtLFp .done2,#mermaid-svg-pGQeQ4y0vVFPtLFp .done3{stroke:grey;fill:#d3d3d3;stroke-width:2}#mermaid-svg-pGQeQ4y0vVFPtLFp .doneText0,#mermaid-svg-pGQeQ4y0vVFPtLFp .doneText1,#mermaid-svg-pGQeQ4y0vVFPtLFp .doneText2,#mermaid-svg-pGQeQ4y0vVFPtLFp .doneText3{fill:#000 !important}#mermaid-svg-pGQeQ4y0vVFPtLFp .crit0,#mermaid-svg-pGQeQ4y0vVFPtLFp .crit1,#mermaid-svg-pGQeQ4y0vVFPtLFp .crit2,#mermaid-svg-pGQeQ4y0vVFPtLFp .crit3{stroke:#f88;fill:red;stroke-width:2}#mermaid-svg-pGQeQ4y0vVFPtLFp .activeCrit0,#mermaid-svg-pGQeQ4y0vVFPtLFp .activeCrit1,#mermaid-svg-pGQeQ4y0vVFPtLFp .activeCrit2,#mermaid-svg-pGQeQ4y0vVFPtLFp .activeCrit3{stroke:#f88;fill:#bfc7ff;stroke-width:2}#mermaid-svg-pGQeQ4y0vVFPtLFp .doneCrit0,#mermaid-svg-pGQeQ4y0vVFPtLFp .doneCrit1,#mermaid-svg-pGQeQ4y0vVFPtLFp .doneCrit2,#mermaid-svg-pGQeQ4y0vVFPtLFp .doneCrit3{stroke:#f88;fill:#d3d3d3;stroke-width:2;cursor:pointer;shape-rendering:crispEdges}#mermaid-svg-pGQeQ4y0vVFPtLFp .milestone{transform:rotate(45deg) scale(0.8, 0.8)}#mermaid-svg-pGQeQ4y0vVFPtLFp .milestoneText{font-style:italic}#mermaid-svg-pGQeQ4y0vVFPtLFp .doneCritText0,#mermaid-svg-pGQeQ4y0vVFPtLFp .doneCritText1,#mermaid-svg-pGQeQ4y0vVFPtLFp .doneCritText2,#mermaid-svg-pGQeQ4y0vVFPtLFp .doneCritText3{fill:#000 !important}#mermaid-svg-pGQeQ4y0vVFPtLFp .activeCritText0,#mermaid-svg-pGQeQ4y0vVFPtLFp .activeCritText1,#mermaid-svg-pGQeQ4y0vVFPtLFp .activeCritText2,#mermaid-svg-pGQeQ4y0vVFPtLFp .activeCritText3{fill:#000 !important}#mermaid-svg-pGQeQ4y0vVFPtLFp .titleText{text-anchor:middle;font-size:18px;fill:#000;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-pGQeQ4y0vVFPtLFp g.classGroup text{fill:#9370db;stroke:none;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family);font-size:10px}#mermaid-svg-pGQeQ4y0vVFPtLFp g.classGroup text .title{font-weight:bolder}#mermaid-svg-pGQeQ4y0vVFPtLFp g.clickable{cursor:pointer}#mermaid-svg-pGQeQ4y0vVFPtLFp g.classGroup rect{fill:#ECECFF;stroke:#9370db}#mermaid-svg-pGQeQ4y0vVFPtLFp g.classGroup line{stroke:#9370db;stroke-width:1}#mermaid-svg-pGQeQ4y0vVFPtLFp .classLabel .box{stroke:none;stroke-width:0;fill:#ECECFF;opacity:0.5}#mermaid-svg-pGQeQ4y0vVFPtLFp .classLabel .label{fill:#9370db;font-size:10px}#mermaid-svg-pGQeQ4y0vVFPtLFp .relation{stroke:#9370db;stroke-width:1;fill:none}#mermaid-svg-pGQeQ4y0vVFPtLFp .dashed-line{stroke-dasharray:3}#mermaid-svg-pGQeQ4y0vVFPtLFp #compositionStart{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-pGQeQ4y0vVFPtLFp #compositionEnd{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-pGQeQ4y0vVFPtLFp #aggregationStart{fill:#ECECFF;stroke:#9370db;stroke-width:1}#mermaid-svg-pGQeQ4y0vVFPtLFp #aggregationEnd{fill:#ECECFF;stroke:#9370db;stroke-width:1}#mermaid-svg-pGQeQ4y0vVFPtLFp #dependencyStart{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-pGQeQ4y0vVFPtLFp #dependencyEnd{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-pGQeQ4y0vVFPtLFp #extensionStart{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-pGQeQ4y0vVFPtLFp #extensionEnd{fill:#9370db;stroke:#9370db;stroke-width:1}#mermaid-svg-pGQeQ4y0vVFPtLFp .commit-id,#mermaid-svg-pGQeQ4y0vVFPtLFp .commit-msg,#mermaid-svg-pGQeQ4y0vVFPtLFp .branch-label{fill:lightgrey;color:lightgrey;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-pGQeQ4y0vVFPtLFp .pieTitleText{text-anchor:middle;font-size:25px;fill:#000;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-pGQeQ4y0vVFPtLFp .slice{font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-pGQeQ4y0vVFPtLFp g.stateGroup text{fill:#9370db;stroke:none;font-size:10px;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-pGQeQ4y0vVFPtLFp g.stateGroup text{fill:#9370db;fill:#333;stroke:none;font-size:10px}#mermaid-svg-pGQeQ4y0vVFPtLFp g.statediagram-cluster .cluster-label text{fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp g.stateGroup .state-title{font-weight:bolder;fill:#000}#mermaid-svg-pGQeQ4y0vVFPtLFp g.stateGroup rect{fill:#ECECFF;stroke:#9370db}#mermaid-svg-pGQeQ4y0vVFPtLFp g.stateGroup line{stroke:#9370db;stroke-width:1}#mermaid-svg-pGQeQ4y0vVFPtLFp .transition{stroke:#9370db;stroke-width:1;fill:none}#mermaid-svg-pGQeQ4y0vVFPtLFp .stateGroup .composit{fill:white;border-bottom:1px}#mermaid-svg-pGQeQ4y0vVFPtLFp .stateGroup .alt-composit{fill:#e0e0e0;border-bottom:1px}#mermaid-svg-pGQeQ4y0vVFPtLFp .state-note{stroke:#aa3;fill:#fff5ad}#mermaid-svg-pGQeQ4y0vVFPtLFp .state-note text{fill:black;stroke:none;font-size:10px}#mermaid-svg-pGQeQ4y0vVFPtLFp .stateLabel .box{stroke:none;stroke-width:0;fill:#ECECFF;opacity:0.7}#mermaid-svg-pGQeQ4y0vVFPtLFp .edgeLabel text{fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .stateLabel text{fill:#000;font-size:10px;font-weight:bold;font-family:'trebuchet ms', verdana, arial;font-family:var(--mermaid-font-family)}#mermaid-svg-pGQeQ4y0vVFPtLFp .node circle.state-start{fill:black;stroke:black}#mermaid-svg-pGQeQ4y0vVFPtLFp .node circle.state-end{fill:black;stroke:white;stroke-width:1.5}#mermaid-svg-pGQeQ4y0vVFPtLFp #statediagram-barbEnd{fill:#9370db}#mermaid-svg-pGQeQ4y0vVFPtLFp .statediagram-cluster rect{fill:#ECECFF;stroke:#9370db;stroke-width:1px}#mermaid-svg-pGQeQ4y0vVFPtLFp .statediagram-cluster rect.outer{rx:5px;ry:5px}#mermaid-svg-pGQeQ4y0vVFPtLFp .statediagram-state .divider{stroke:#9370db}#mermaid-svg-pGQeQ4y0vVFPtLFp .statediagram-state .title-state{rx:5px;ry:5px}#mermaid-svg-pGQeQ4y0vVFPtLFp .statediagram-cluster.statediagram-cluster .inner{fill:white}#mermaid-svg-pGQeQ4y0vVFPtLFp .statediagram-cluster.statediagram-cluster-alt .inner{fill:#e0e0e0}#mermaid-svg-pGQeQ4y0vVFPtLFp .statediagram-cluster .inner{rx:0;ry:0}#mermaid-svg-pGQeQ4y0vVFPtLFp .statediagram-state rect.basic{rx:5px;ry:5px}#mermaid-svg-pGQeQ4y0vVFPtLFp .statediagram-state rect.divider{stroke-dasharray:10,10;fill:#efefef}#mermaid-svg-pGQeQ4y0vVFPtLFp .note-edge{stroke-dasharray:5}#mermaid-svg-pGQeQ4y0vVFPtLFp .statediagram-note rect{fill:#fff5ad;stroke:#aa3;stroke-width:1px;rx:0;ry:0}:root{--mermaid-font-family: '"trebuchet ms", verdana, arial';--mermaid-font-family: "Comic Sans MS", "Comic Sans", cursive}#mermaid-svg-pGQeQ4y0vVFPtLFp .error-icon{fill:#522}#mermaid-svg-pGQeQ4y0vVFPtLFp .error-text{fill:#522;stroke:#522}#mermaid-svg-pGQeQ4y0vVFPtLFp .edge-thickness-normal{stroke-width:2px}#mermaid-svg-pGQeQ4y0vVFPtLFp .edge-thickness-thick{stroke-width:3.5px}#mermaid-svg-pGQeQ4y0vVFPtLFp .edge-pattern-solid{stroke-dasharray:0}#mermaid-svg-pGQeQ4y0vVFPtLFp .edge-pattern-dashed{stroke-dasharray:3}#mermaid-svg-pGQeQ4y0vVFPtLFp .edge-pattern-dotted{stroke-dasharray:2}#mermaid-svg-pGQeQ4y0vVFPtLFp .marker{fill:#333}#mermaid-svg-pGQeQ4y0vVFPtLFp .marker.cross{stroke:#333}:root { --mermaid-font-family: "trebuchet ms", verdana, arial;} #mermaid-svg-pGQeQ4y0vVFPtLFp {color: rgba(0, 0, 0, 0.75);font: ;}

Start
初始化参数
生成H*
生成F*
MSE
输出误差E
更新参数
Stop

听起来很简单，实现起来也不复杂，我们接下来看看算法的实现细节。

实现步骤

LMS 滤波器背后的基本思想是通过以收敛到最佳滤波器的方式更新滤波器权重 R − 1 P R^{-1} P R−1P ，采取渐进地逼近策略，也就是梯度下降算法的根本思想。

对于我们来说 F ^ ( u , v ) = G ( u , v ) ∗ H ( u , v ) − 1 − D ( u , v ) \hat F(u, v) = G(u, v) * H(u, v)^{-1} - D(u, v) F^(u,v)=G(u,v)∗H(u,v)−1−D(u,v) 其实可以看作矩阵的线性计算， H ( u , v ) − 1 H(u, v)^{-1} H(u,v)−1 可以看作权重 ω ( u , v ) \omega(u, v) ω(u,v)，输入参数为 G ( u , v ) G(u, v) G(u,v)，输出参数为 F ^ ( u , v ) \hat F(u, v) F^(u,v)。噪音 D ( u , v ) D(u, v) D(u,v) 可以当作偏移量 b b b，于是可以把运算过程表达为在矩阵中的基本线性计算

y ^ = ω x + b \hat y = \omega x + b y^​=ωx+b

然后我们引入权重更新公式

ω ^ = ω − λ ∇ ( e 2 ) \hat \omega = \omega - \lambda \nabla (e^2) ω^=ω−λ∇(e2)

补充信息
梯度下降的权重更新公式推导比较麻烦，我在这篇文章里就不再赘述，有需要的朋友可以看看我这篇文章的内容
深度学习知识总结—— 2. 计算图与反向传播

由于梯度在计算过程中可能为负，也可能为正。也就是说，如果 MSE 梯度为正，则意味着将相同的权重用于进一步迭代，误差将继续正增加，这意味着我们需要减少权重。同样，如果梯度为负，我们需要增加权重。

为了方便计算，对于 e 2 e^2 e2 的求导运算，我们可以换成对梯度的差分计算。

补充信息
关于梯度的差分运算，有需要的朋友可以看看我这篇文章的内容
浅谈矢量场 —— 1. 梯度、散度与拉普拉斯算子

然后，我们把以上信息组合后，就可以得到 「最小均方差滤波」的实现方法。

  1. 先将原始图像从 f ( x , y ) f(x, y) f(x,y) 转化为 F ( u , v ) F(u,v) F(u,v)；
  2. 再讲退化图像从 g ( x , y ) g(x, y) g(x,y) 转化为 G ( u , v ) G(u, v) G(u,v)；
  3. 初始化一个大小和 G ( u , v ) G(u, v) G(u,v) 一样的权重矩阵 W ( u , v ) W(u, v) W(u,v)；
  4. 设定更新率 λ \lambda λ，通常是小于 < 1 < 1 <1 的浮点数；
  5. 由权重矩阵 W W W 和原始图像 F F F 计算得出一个退化图像 G G G；
  6. 计算方差 M S E MSE MSE；
  7. 更新权重矩阵 W W W，重复第5步；

看起来有点复杂，但是基本逻辑很简单。

实现代码

首先，我们要创建网络权重矩阵

def create_parameters(dft):# generate empty matrixweight = np.ones_like(dft, dtype=np.float32)# return to callerreturn weight

然后创建一个用来计算均方差的函数，它主要帮助我们判断算法收敛情况。

def compute_loss(dft_predicated, dft_target):# compute mse of original and recovery imagesnew_errors = (dft_predicated - dft_target)**2# return sumreturn np.sum(new_errors)

然后一个通过导数更新参数的模块。为了方便理解，我们把导数部分以绝对值的形式表现出来，这样对于我们来说，接下来只需要关心在不同方向上，梯度的增减方向了。

def backward(weights, dft_predicated, dft_target, lambs):# update value chainupdated_val = np.absolute(lambs * (dft_predicated - dft_target))width, height, channel = dft_predicated.shapefor w in range(width):for h in range(height):for c in range(channel):if dft_target[w, h, c] >= dft_predicated[w, h, c] >= 0:if weights[w, h, c] >= 0:weights[w, h, c] = weights[w, h, c] + updated_val[w, h, c]else:weights[w, h, c] = weights[w, h, c] - updated_val[w, h, c]continueelif dft_predicated[w, h, c] >= dft_target[w, h, c] >= 0:if weights[w, h, c] >= 0:weights[w, h, c] = weights[w, h, c] - updated_val[w, h, c]else:weights[w, h, c] = weights[w, h, c] + updated_val[w, h, c]continueelif 0 >= dft_target[w, h, c] >= dft_predicated[w, h, c]:if weights[w, h, c] >= 0:weights[w, h, c] = weights[w, h, c] - updated_val[w, h, c]else:weights[w, h, c] = weights[w, h, c] + updated_val[w, h, c]continueelif 0 >= dft_predicated[w, h, c] >= dft_target[w, h, c]:if weights[w, h, c] >= 0:weights[w, h, c] = weights[w, h, c] + updated_val[w, h, c]else:weights[w, h, c] = weights[w, h, c] - updated_val[w, h, c]continueelse:# dft_target[w, h, c] > 0 >= dft_predicated[w, h, c]:# dft_predicated[w, h, c] > 0 >= dft_target[w, h, c]:if weights[w, h, c] >= 0:weights[w, h, c] = weights[w, h, c] - updated_val[w, h, c]else:weights[w, h, c] = weights[w, h, c] + updated_val[w, h, c]continue# return the updated weightsreturn weights

然后一个前向运算

def forward(weights, dft_origin):return weights * dft_origin  # predicated_y = weight * 1

这里稍微注意一点，我们引入一个标准化模块，因为DFT转化后数值特别大，这会导致梯度爆炸，所以需要适度的压缩数据。压缩算法有很多种，比如指数或逆指数；但是对于我们来说我们更关心全局情况，而且没有数据特征上的偏好，所以传统的 normalization 方法比较符合我们的需要。

def normalize_data(dft1, dft2):max_dft = np.maximum(dft1, dft2)max_data = np.max(max_dft.reshape(-1))min_dft = np.minimum(dft1, dft2)min_data = np.min(min_dft.reshape(-1))length = max_data - min_datadft1 = dft1 / length * 100.0dft2 = dft2 / length * 100.0return dft1, dft2

然后就是通过梯度下降把需要的网络权重算出来，这样可以得到我们需要的滤波器。

def gradient_descent(img, distort):# to compress the data to avoid gradient explosiondft_origin, dft_distort = normalize_data(img_to_dft(img), img_to_dft(distort))# lambda ratelambs = 0.01# create weightweight = create_parameters(dft_origin)# loop and descent gradients for DFT weightsfor i in range(500):# forwarddft_predicated = forward(weight, dft_origin)# compute lossmse = compute_loss(dft_predicated, dft_distort)# backward computationweight = backward(weight, dft_predicated, dft_distort, lambs)# print debug messageprint("mse: %.4f" % mse)return weightdef display_result(diagrams):plot = DiagramPlotter()plot.append_image(diagrams[0], "Original")plot.append_image(diagrams[1], "Violent k=0.0015")plot.append_image(diagrams[2], "Recovery from Violent with k=-0.0010")plot.append_image(diagrams[2], "Recovery from LMS Filter")plot.show(1, 4)def recovery_image(img, distorted, disrecv):# copy for backupdistorted_backup = distorted.copy()# obtain weightweight = gradient_descent(img, distorted_backup)# use LMS filter to recovery distorted imagedft = img_to_dft(distorted)dft = dft / weightimg_ds = dft_to_img(dft)# show resultdisplay_result((img, distorted, disrecv, img_ds))

以上代码并非最优，其实可优化的地方很多，不过主要的业务逻辑就这么多，我们来看看最后的效果。

最后的结果


最右是LMS过滤器经过500次迭代，我的MSE只收敛到约 130.7263左右，尽管还可以继续，但作为示例就没必要继续花费时间了，可以看到效果恢复确实不错。参数我是根据一些经验调整的，你如果有兴趣可以修改代码看看，是不是能更快收敛到最佳效果。

理论上来说，只要梯度持续下降，那么可以达到一个全局最优的情况，应该还可以得到几乎接近原图的效果。作为自适应滤波器，LMS 基本上也能应付绝大多数场景需要，除了那种数据信息被完全摧毁的。

当然，这种算法也有问题，想要它找出最合适的模型，需要花费挺久的时间，而且需要有标签数据进行比对。

  1. https://zh.wikipedia.org/wiki/%E6%9C%80%E5%B0%8F%E5%9D%87%E6%96%B9%E6%BB%A4%E6%B3%A2%E5%99%A8 ↩︎

数字图像学笔记 —— 16. 图像退化与复原（自适应滤波之「最小均方差滤波」）相关推荐

  1. 数字图像学笔记——13. 图像退化与复原（退化函数的评估方法：观察法、实验法、数学建模法与湍流导致的退化）

    在对受到多种原因影响的图像进行复原时,我们经常需要先行评估对图像质量产生影响的退化函数,有时甚至需要尝试建模.通过这些手段,能够最大程度上恢复图像上的噪音,并重建高清的图像细节. 文章目录 线性位置不 ...

  2. 数字图像学笔记——14. 图像退化与复原（线性退化）

    文章目录 运动导致的退化(线性退化) 水平运动导致的退化 垂直运动导致的退化 运动导致的退化(线性退化) 在上一章 <数字图像学笔记--13. 图像退化与复原(退化函数的评估方法:观察法.实验法 ...

  3. 数字图像处理——图像退化与复原

    图像退化与复原 内容简介 1.图像退化与复原的原理 1.1 图像退化的数学模型 1.2 图像退化的原理 1.3 图像复原的原理 2. 图像去噪 2.1 噪声模型 2.2 噪声参数的估计 2.3 针对噪 ...

  4. 数字图像学笔记——7. 噪音生成（泊松噪音生成方法）

    文章目录 泊松噪音 Knuth算法 散列生成算法 生成泊松噪音的图像 泊松噪音 Knuth算法 首先,回顾泊松分布的函数: P(x=k)=e−λλkk!P(x=k) = \frac{e^{- \lam ...

  5. 数字图像学笔记——3.彩色转黑白

    文章目录 一些说明 关于示例代码 关于依赖环境 关于教材 灰度图.亮度图(Gray Image) 彩色图转灰度图 一般亮度转换(luminosity method) 亮度优先转换(luminosity ...

  6. 数字图像学笔记——4. 直方图计算、线性变换、对数变换、Gamma变换

    文章目录 灰度直方图(Gray Histogram) 直方图的计算方法 简单的图像转换方法 线性变换 / 图像翻转(Image Nagatives) 对数变换(Log Transformation) ...

  7. 数字图像学笔记——10. 频域与傅里叶分析方法

    频域滤波技术,目前主要使用的有两种类型,一种是傅立叶变换技术,还有一种是小波分析.基本逻辑就是把原始信号映射到频率空间中,使得在时域空间无法处理的信号,得以在另外一种空间体系下能够被较有效的处理. 目 ...

  8. matlab 图像退化,数字图像退化与复原系统设计doc完整版(MATLAB).doc

    数字图像退化与复原系统设计doc完整版(MATLAB) 题一 数字图像退化与复原系统设计 摘要 题一针对数字图像退化与复原系统设计,利用MATLAB软件设计出了较为合理的用户界面并基本上实现了系统的功 ...

  9. 图像退化 / 复原处理的模型

    文章目录 一.概述 二.图像退化/复原处理的模型 一.概述 图形复原技术的主要目的是以预先确定的目标来改善图像 图像复原试图利用退化现象的某种先验知识来复原一幅退化的图像,因而复原技术是面向退化模型的 ...

最新文章

  1. python模块之pickle
  2. python是全栈_Python全栈之路-3-字符串
  3. 有理数加减乘除 计算机应用带答案,列50道有理数的混合运算（加减乘除）包括答案 初一的...
  4. maven 如何看jar是否被修改_如何在线修改jar文件
  5. 【收集】EJB3.0的各应用服务器提供的JNDI接口
  6. 【HDU6667】Roundgod and Milk Tea【贪心】
  7. 【专访】PP租车孙览江：与有梦想的人一拍即合，PM都有改变世界的小情怀
  8. windows+mysql+解压版_Windows操作系统安装MySQL解压版
  9. C/C++ unsigned char*类型
  10. 由“深”至“广”，探索2022音视频技术的无限可能
  11. linux mysql cron_定时MySQL动作－Linux下用Cron现定时执行脚本
  12. PyTorch | torch.randperm()使用方法
  13. python一个函数调用另一个函数的返回值_python-调用另一个函数后立即从函数返回...
  14. [数据库]oracle导出数据库
  15. php post 视频教程,PHP教程：POST数据的三种方法
  16. 在ASP.NET 2.0中实现URL重写
  17. JNPF.java前后端分离框架，SpringBoot+SpringCloud开发微服务平台
  18. html看汉字选拼音小游戏
  19. 电子烟能破壳类四大天王“和大天壹”新物种么？
  20. 当计算机遇上经济学：如何量化你的投资并获得第一桶金

热门文章

  1. 开机运行到window启动画面马上蓝屏重启
  2. 基于三菱PLC的饮料售货机控制系统设计
  3. 高通电源管理qpnp-vm-bms驱动-电量计
  4. Java岗面试必问！mysql视频教程百度云
  5. 异步同步通信数据帧格式
  6. 以太网数据帧详细解析 逐字节分析
  7. STM32F103 CAN通讯实操
  8. navicat for mysql 注册码
  9. 数据库设计时的一些注意事项
  10. python基础（条件及循环）