Bayes,HMM,CRF & Gibbs Distribution在图像降噪中的应用

  • 1 问题描述
  • 2 数学模型
    • 2.1Bayes推断模型
    • 2.2 HMM模型
    • 2.3 CRF模型
    • 2.4 Gibbs分布
  • 3 快速算法
  • 4 代码与结果分析

1 问题描述

假设原始图像 X m × n \mathbf{X}_{m\times n} Xm×n​,第 i i i行,第 j j j列的灰度为 x i , j x_{i,j} xi,j​,经过高斯噪声后
Y = X + N \mathbf{Y}=\mathbf{X}+\mathbf{N} Y=X+N 其中 N \mathbf{N} N是 m × n m \times n m×n噪声矩阵,第 i i i行,第 j j j列的灰度为 N i , j N_{i,j} Ni,j​, N i j ∼ N ( 0 , σ 2 ) N_{ij} \sim \mathcal{N}(0,\sigma^2) Nij​∼N(0,σ2)。例如,原始图像 I \mathbf{I} I为

添加噪声后得到观测 Y \mathbf{Y} Y

问题是,如何利用像素之间的相互关系,基于 Y \mathbf{Y} Y 恢复 X \mathbf{X} X

2 数学模型

2.1Bayes推断模型

将降噪问题建模为基于观测 Y \mathbf{Y} Y推断 X \mathbf{X} X问题,基于最大后验概率准则(MAP:Maximum A Posteriori):
X ^ = arg ⁡ max ⁡ X ℓ = p ( X ∣ Y ) \hat{\mathbf{X}}=\arg \max_\mathbf{X} \ell=p(\mathbf{X}|\mathbf{Y}) X^=argXmax​ℓ=p(X∣Y)根据Bayes公式:
p ( X ∣ Y ) = p ( Y ∣ X ) p ( X ) p ( Y ) p(\mathbf{X}|\mathbf{Y}) = \frac{p(\mathbf{Y}|\mathbf{X})p(\mathbf{X})}{p(\mathbf{Y})} p(X∣Y)=p(Y)p(Y∣X)p(X)​
优化问题等价于:
X ^ = arg ⁡ max ⁡ X ℓ ′ = p ( Y ∣ X ) p ( X ) \hat{\mathbf{X}}=\arg \max_\mathbf{X} \ell'=p(\mathbf{Y}|\mathbf{X})p(\mathbf{X}) X^=argXmax​ℓ′=p(Y∣X)p(X)其中似然函数 p ( Y ∣ X ) p(\mathbf{Y}|\mathbf{X}) p(Y∣X)由噪声模型给出(假设各像素的噪声相互独立):
p ( Y ∣ X ) = 1 2 π σ exp ⁡ { − ∑ i = 1 m ∑ j = 1 n ( y i , j − x i , j ) 2 2 σ 2 } p(\mathbf{Y}|\mathbf{X})=\frac{1}{\sqrt{2\pi}\sigma}\exp{\left\{-\sum_{i=1}^m\sum_{j=1}^n\frac{(y_{i,j}-x_{i,j})^2}{2\sigma^2}\right\}} p(Y∣X)=2π ​σ1​exp{−i=1∑m​j=1∑n​2σ2(yi,j​−xi,j​)2​}问题的主要困难有两处:

  1. 如何对 X \mathbf{X} X的先验分布 p ( X ) p(\mathbf{X}) p(X)进行准确建模?
  2. X \mathbf{X} X的搜索空间非常庞大,如何高效求解上述优化模型?

2.2 HMM模型

先验分布 p ( X ) p(\mathbf{X}) p(X)表示在观测前,对于图像 X \mathbf{X} X的了解程度。根据概率链式法则计算先验概率 X \mathbf{X} X:
p ( X ) = p ( x n ∣ x n − 1 , x n − 2 , . . . , x 1 ) ⋯ p ( x 3 ∣ x 2 , x 1 ) p ( x 2 ∣ x 1 ) p ( x 1 ) p(\mathbf{X})=p(x_n|x_{n-1},x_{n-2},...,x_1 )\cdots p(x_3|x_2,x_1)p(x_2|x_1)p(x_1) p(X)=p(xn​∣xn−1​,xn−2​,...,x1​)⋯p(x3​∣x2​,x1​)p(x2​∣x1​)p(x1​)上述展开是很困难的,为了简化模型,可以将概率模型进行简化,认为 X \mathbf{X} X的各个分量按顺序形成一条马尔科夫链:

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

x0
x1
x2
x3

此时有
p ( X ) = p ( x n ∣ x n − 1 ) ⋯ p ( x 3 ∣ x 2 ) p ( x 2 ∣ x 1 ) p ( x 1 ) p(\mathbf{X})=p(x_n|x_{n-1})\cdots p(x_3|x_2)p(x_2|x_1)p(x_1) p(X)=p(xn​∣xn−1​)⋯p(x3​∣x2​)p(x2​∣x1​)p(x1​) 这就是隐马尔科夫模型 HMM(HMM: Hiden Markov Model)所采用的隐状态概率模型。

HMM 是Bayes推断的一种特例。HMM模型对Bayes推断中的 p ( X ) p(\mathbf{X}) p(X)进行简化假设。HMM假设隐状态 X \mathbf{X} X的各个随机变量 x i x_i xi​是条马尔科夫链,即当前状态 x i x_i xi​只与上一个状态 x i − 1 x_{i-1} xi−1​有关,而与如何到达上一个状态的过程无关。
HMM的目标是基于观测 Y \mathbf{Y} Y推断 X \mathbf{X} X,先验知识是观测模型 Y = X + N \mathbf{Y}=\mathbf{X}+\mathbf{N} Y=X+N,以及状态转移概率模型 p ( x i ∣ x i − 1 ) p(x_i|x_{i-1}) p(xi​∣xi−1​),通常HMM采用Bayes方法,即MAP构建目标函数,采用Viterbbi算法求解。Viterbbi算法的基本原理是动态规划(DP: Dynamic Programming)方法。
在语音识别中,常用HMM模型,将音素(phoneme)作为 x i x_i xi​,构建相邻两个音素之间的概率关系,构成一条马尔科夫链,运用HMM的方法估计 X \mathbf{X} X。
为了进一步贴近事实,也有部分研究对HMM进行了扩展,构建高阶HMM,例如三音素HMM模型,即认为当前音素与前三个因素的状态有关,构建三阶HMM,后面做法与单阶HMM相同。
更多HMM内容参考这里。

2.3 CRF模型

HMM只能表示链状概率转移关系,而有些场景下,链状关系无法准确描述不同随机变量的概率依存关系,这就需要构建MRF(Markov Random Feild)模型了。MRF是链状Markov过程的扩展。MRF中常用农场举例说明。

农场的一个点上种植的作物这个随机变量 x i x_i xi​,仅与临近的地种植的作物品种的随机变量 x j x_j xj​有关,而且大概率与临近种植作物相同,例如
p ( x i = 莴 苣 ∣ x j = 莴 苣 ) = 0.8 , p ( x i = 非 莴 苣 ∣ p j = 莴 苣 ) = 0.2 ) p(x_i=莴苣|x_j=莴苣)=0.8, p(x_i=非莴苣|p_j=莴苣)=0.2) p(xi​=莴苣∣xj​=莴苣)=0.8,p(xi​=非莴苣∣pj​=莴苣)=0.2)
如果每个隐变量还有一个观测 Y \mathbf{Y} Y,我们需要根据观测推断隐变量 X \mathbf{X} X的概率分布,即估计后验概率 p ( X ∣ Y ) p(\mathbf{X}|\mathbf{Y}) p(X∣Y),则构成条件随机场(CRF:Conditional Random Field)
用概率图表示就是:

图中一个蓝色的圈 x i j x_{ij} xij​,指的是一个像素的真实状态,对应绿色的圈是该真实状态的观测灰度值 y i , j y_{i,j} yi,j​。每个蓝色圈对应的像素的真实灰度值与邻近的4个像素的颜色相关,并且被观测到绿色的结果 y i j y_{ij} yij​

区别于Bayes网络,MRF是一个无向图,适合于解释没有明确因果关系的网络关系。
针对图像而言,先验知识是图像大概率是连续的,光滑的,即图像像素的灰度大概率与临近像素的灰度接近。为了简单,可以假设每个像素的灰度只与4邻域或者8邻域像素的灰度有关,不妨设:
p ( x i ˉ , j ˉ ∣ x i , j ) = 1 2 π σ x exp ⁡ { − ( x i ˉ , j ˉ − x i , j ) 2 2 σ x 2 } p(x_{\bar{i},\bar{j}}|x_{i,j})=\frac{1}{\sqrt{2\pi}\sigma_x}\exp{\left\{\frac{-(x_{\bar{i},\bar{j}}-x_{i,j})^2}{2\sigma_x^2}\right\}} p(xiˉ,jˉ​​∣xi,j​)=2π ​σx​1​exp{2σx2​−(xiˉ,jˉ​​−xi,j​)2​}其中 x i ˉ , j ˉ x_{\bar{i},\bar{j}} xiˉ,jˉ​​表示第 i i i行第 j j j列的某个相邻像素的灰度值,例如4邻域中, ( i ˉ , j ˉ ) ∈ { ( i − 1 , j ) , ( i + 1 , j ) , ( i , j − 1 ) , ( i , j + 1 ) } (\bar{i},\bar{j})\in\{(i-1,j),(i+1,j),(i,j-1),(i,j+1)\} (iˉ,jˉ​)∈{(i−1,j),(i+1,j),(i,j−1),(i,j+1)}。 σ v \sigma_v σv​是相邻像素灰度差异的标准差,该数值可以通过统计大量自然图像估计得到。
基于上述模型,可以得到一个局部的MRF:

该局部模型的Bayes推断模型为:
p ( y i , j ∣ x i , j ) = 1 2 π σ exp ⁡ { − ( y i , j − x i , j ) 2 2 σ 2 } p ( x i ˉ , j ˉ ∣ x i , j ) = 1 2 π σ v exp ⁡ { − ( x i ˉ , j ˉ − x i , j ) 2 2 σ v 2 } p(y_{i,j}|x_{i,j})=\frac{1}{\sqrt{2\pi}\sigma}\exp{\left\{-\frac{(y_{i,j}-x_{i,j})^2}{2\sigma^2}\right\}}\\ p(x_{\bar{i},\bar{j}}|x_{i,j})=\frac{1}{\sqrt{2\pi}\sigma_v}\exp{\left\{-\frac{(x_{\bar{i},\bar{j}}-x_{i,j})^2}{2\sigma_v^2}\right\}} p(yi,j​∣xi,j​)=2π ​σ1​exp{−2σ2(yi,j​−xi,j​)2​}p(xiˉ,jˉ​​∣xi,j​)=2π ​σv​1​exp{−2σv2​(xiˉ,jˉ​​−xi,j​)2​}其中前式为观测模型,后式为状态关系模型。
若用Bayes推断的方式估计最优 X \mathbf{X} X,需要计算 p ( Y ∣ X ) p(\mathbf{Y}|\mathbf{X} ) p(Y∣X)和 p ( X ) p(\mathbf{X}) p(X)。其中,假设观测噪声独立,则有:
p ( Y ∣ X ) = ∑ i = 1 m ∑ j = 1 n p ( y i , j ∣ x i , j ) p(\mathbf{Y}|\mathbf{X})=\sum_{i=1}^m \sum_{j=1}^n p(y_{i,j}|x_{i,j}) p(Y∣X)=i=1∑m​j=1∑n​p(yi,j​∣xi,j​)现在新的问题是需要根据局部条件概率 p ( x i ˉ , j ˉ ∣ x i , j ) p(x_{\bar{i},\bar{j}}|x_{i,j}) p(xiˉ,jˉ​​∣xi,j​),估计联合概率分布 p ( X ) p(\mathbf{X}) p(X)。也就是说,需要基于概率图中每条边的概率转移关系,推断联合概率分布。传统方法是基于链式法则展开,但展开式会非常复杂。这里就需要引入Gibbs分布的概念。

2.4 Gibbs分布

MRF研究的两个重要问题是:

  1. 如何根据概率图中边的概率关系 p ( x i ∣ x j ) p(x_i|x_j) p(xi​∣xj​)推断所有随机变量的联合概率 p ( X ) p(\mathbf{X}) p(X);
  2. 如何根据概率图的概率关系生成服从 p ( X ) p(\mathbf{X}) p(X)分布的随机样本,主要用于马尔科夫蒙特卡洛(MCMC: Markov Chain Monte Carlo),该问题可以使用metropolis-hasting算法(低维随机数生成有效)和它的特例Gibbs采样算法(主要针对高维随机数生成)解决。

针对第二种方法,在后续文中讨论,这里主要利用MRF研究的第一个问题解决图像降噪问题。
Hammersley Clifford Theorem给出了马尔科夫随机场的联合概率 p ( X ) p(\mathbf{X}) p(X)的计算方法,即MRF的联合概率 p ( X ) p(\mathbf{X}) p(X)和Gibbs分布是一致的。也就是说,可以使用Gibbs分布估计 p ( X ) p(\mathbf{X}) p(X)。

Hammersley Clifford Theorem指出:

  1. Gibbs分布一定满足由node separation导致的条件独立性
  2. 马尔科夫随机场的概率分布一定可以表示成最大团(Maximum Clique)上的非负函数乘积形式
    证明方法可以看这里

概率图中的”团”指的是概率图中的一个全连接子图,即“团”内的任意两点都有边连接。“最大团”指的是不被任何其他“团”包含的团。在图像降噪问题的例子中,没有两个以上的全连接子图,最大团是包含两个节点一条边的团,所以图像降噪例子的各个像素灰度的联合分布可以表示成每条边上一个非负函数乘积的形式:
p ( X ∣ Y ) = 1 z ∏ i = 1 m ∏ j = 1 m exp ⁡ F i , j ( X , Y ) ) p(\mathbf{X}|\mathbf{Y})=\frac{1}{z}\prod_{i=1}^m \prod_{j=1}^m \exp{F_{i,j}(\mathbf{X,Y}))} p(X∣Y)=z1​i=1∏m​j=1∏m​expFi,j​(X,Y))其中 z z z是归一化系数, F i , j ( X ) F_{i,j}(\mathbf{X}) Fi,j​(X)表示与像素 x i , j x_{i,j} xi,j​有关的势函数,它包含两部分,一部分是像素 x i , j x_{i,j} xi,j​与观测 y i , j y_{i,j} yi,j​的势函数 F o ( x i , j , y i , j ) F_o(x_{i,j},y_{i,j}) Fo​(xi,j​,yi,j​),另一部分是像素 x i , j x_{i,j} xi,j​与临近像素 x i − 1 , j x_{i-1,j} xi−1,j​和 x i , j − 1 x_{i,j-1} xi,j−1​的势函数 F g ( x i , j , x i − 1 , j , x i , j − 1 ) F_g(x_{i,j},x_{i-1,j},x_{i,j-1}) Fg​(xi,j​,xi−1,j​,xi,j−1​),并且
F i , j ( X ) = F o ( x i , j , y i , j ) F g ( x i , j , x i − 1 , j , x i , j − 1 ) F_{i,j}(\mathbf{X})=F_o(x_{i,j},y_{i,j})F_g(x_{i,j},x_{i-1,j},x_{i,j-1}) Fi,j​(X)=Fo​(xi,j​,yi,j​)Fg​(xi,j​,xi−1,j​,xi,j−1​)为了避免重复计算,只计算每个像素与左侧和下侧像素的势函数 F g ( x i , j , x − i , j ) F_g(x_{i,j},x_{-i,j}) Fg​(xi,j​,x−i,j​),为了避免第一行和第一列计算异常,定义 x i , j = 0 x_{i,j}=0 xi,j​=0,if i = 0 i=0 i=0 or j = 0 j=0 j=0)。
在高斯噪声假设和相邻像素灰度差异服从高斯分布的假设下,定义势函数为:
F o ( x i , j , y i , j ) ≜ − 1 2 σ 2 ( x i , j − y i , j ) 2 F g ( x i , j , x i − 1 , j , x i , j − 1 ) ≜ − 1 2 σ v 2 ( x i , j − x i − 1 , j ) 2 + ( x i , j − x i , j − 1 ) 2 F_o(x_{i,j},y_{i,j}) \triangleq -\frac{1}{2 \sigma^2}(x_{i,j}-y_{i,j})^2 \\ F_g(x_{i,j},x_{i-1,j},x_{i,j-1}) \triangleq -\frac{1}{2\sigma_v^2} (x_{i,j}-x_{i-1,j})^2+(x_{i,j}-x_{i,j-1})^2 Fo​(xi,j​,yi,j​)≜−2σ21​(xi,j​−yi,j​)2Fg​(xi,j​,xi−1,j​,xi,j−1​)≜−2σv2​1​(xi,j​−xi−1,j​)2+(xi,j​−xi,j−1​)2
根据最大后验概率(MAP)准则,对后验概率函数取对数,去除常数项后:
X ^ = arg ⁡ min ⁡ ℓ ′ ′ = ∑ i = 1 m ∑ j = 1 n { λ o ( x i , j − y i , j ) 2 + λ g [ ( x i , j − x i − 1 , j ) 2 + ( x i , j − x i , j − 1 ) 2 ] } \hat{\mathbf{X}} = \arg \min \ell'' = \sum_{i=1}^m \sum_{j=1}^n \left\{ \lambda_o (x_{i,j}-y_{i,j})^2 + \lambda _g \left[ (x_{i,j}-x_{i-1,j})^2+(x_{i,j}-x_{i,j-1})^2 \right ] \right\} X^=argminℓ′′=i=1∑m​j=1∑n​{λo​(xi,j​−yi,j​)2+λg​[(xi,j​−xi−1,j​)2+(xi,j​−xi,j−1​)2]}其中 λ o = 1 2 σ 2 \lambda_o=\frac{1}{2\sigma^2} λo​=2σ21​表示观测的权重, λ g = 1 2 σ v 2 \lambda_g=\frac{1}{2\sigma_v^2} λg​=2σv2​1​表示图像先验知识的权重, x 0 , j ≜ x 1 , j , x i , 0 ≜ x i , 1 x_{0,j}\triangleq x_{1,j},x_{i,0} \triangleq x_{i,1} x0,j​≜x1,j​,xi,0​≜xi,1​。事实上,对比Bayes推断框架,前者对应观测似然函数 p ( Y ∣ X ) p(\mathbf{Y|X}) p(Y∣X),后者对应先验知识似然函数 p ( X ) p(\mathbf{X}) p(X)。

经过上述过程,构建起了推断 X \mathbf{X} X的数学模型。在求解该优化问题时,就需要遍历所有的 X \mathbf{X} X,找出最大后验的解,此时搜索规模是 25 6 m × n 256^{m\times n} 256m×n,显然是无法接受的,需要设计有效的搜索算法。

3 快速算法

优化目标是搜索最佳的 x i , j x_{i,j} xi,j​,使两类距离(观测与估计距离&相邻像素之间距离)最小。若能每步迭代都减小一点距离,并且有个好的初值,就能以比较大的概率收敛。考虑以下迭代方法:
x ^ i , j ( t + 1 ) = arg ⁡ min ⁡ x i , j ( t + 1 ) ℓ i , j ( t + 1 ) = λ o ( x i , j − y i , j ) 2 + λ g [ ( x i , j − x i − 1 , j ( t ) ) 2 + ( x i , j − x i , j − 1 ( t ) ) 2 ] (1) \hat{x}_{i,j}^{(t+1)}=\arg \min_{x_{i,j}^{(t+1)}} \ell_{i,j} ^{(t+1)}= \lambda_o (x_{i,j}-y_{i,j})^2 + \lambda _g \left[ (x_{i,j}-x_{i-1,j}^{(t)})^2+(x_{i,j}-x_{i,j-1}^{(t)})^2 \right ] \tag{1} x^i,j(t+1)​=argxi,j(t+1)​min​ℓi,j(t+1)​=λo​(xi,j​−yi,j​)2+λg​[(xi,j​−xi−1,j(t)​)2+(xi,j​−xi,j−1(t)​)2](1) 直观意义就是优化第 i i i行第 j j j列的灰度值 x i , j x_{i,j} xi,j​,使其与观测 y i , j y_{i,j} yi,j​和临近像素灰度值 x i − 1 , j , x i , j − 1 x_{i-1,j},x_{i,j-1} xi−1,j​,xi,j−1​尽可能接近。基于Gibbs抽样的相关结论,不难证明,该迭代过程是收敛于MAP估计的(具体证明方法在后续文章给出)。
一个好的初值能够加速收敛过程,注意到真实图像 X \mathbf{X} X应该接近观测 Y \mathbf{Y} Y,所以使用 Y \mathbf{Y} Y作为 X \mathbf{X} X的初值是合理的。
综上得到算法流程:

Created with Raphaël 2.3.0 开始 初始化X=Y 计算目标函数值 终止条件? 结束 逐像素更新 yes no

其中逐像素更新模块为求解上述优化模型(1),迭代终止条件为图像灰度值不再发生变化。

4 代码与结果分析

代码如下,代码部分参考了这里:

% close all;
clear all;
figure;
I=rgb2gray(imread('lena.jpg'));
imwrite(I,'lena_gray.jpg');
I1=double(I);
subplot(2,2,1)
imagesc(I1);
title('原图像');
colormap('gray')
J1 = double(imnoise(I,'gaussian',0,0.001));
subplot(2,2,2)
imagesc(J1);
title('噪声图')
colormap('gray')
imwrite(uint8(J1),'lena_noised.jpg');Y = J1;
[m,n]=size(Y);
X = Y;
h=0;beta=0.05;eta=1;
z = [-1,0,1];
while 1tot = 0;for i=2:1:m-1for j=2:1:n-1temp = X(i,j);sq = X(i-1:i+1,j-1:j+1);sq1 = sq;sq1(2,2) = sq(2,2)-1;   %根据定义计算势函数E1 = beta*((sq1(2,2)-sq1(1,2))^2+(sq1(2,2)-sq1(2,1))^2+(sq1(2,2)-sq1(2,3))^2+(sq1(2,2)-sq1(3,2))^2)+eta*(sq1(2,2)-Y(i,j))^2/9;sq2 = sq;sq2(2,2) = sq(2,2);   %根据定义计算势函数E2 = beta*((sq2(2,2)-sq2(1,2))^2+(sq2(2,2)-sq2(2,1))^2+(sq2(2,2)-sq1(2,3))^2+(sq2(2,2)-sq2(3,2))^2)+eta*(sq2(2,2)-Y(i,j))^2/9;sq3 = sq;sq3(2,2) = sq(2,2)+1;   %根据定义计算势函数E3 = beta*((sq3(2,2)-sq3(1,2))^2+(sq3(2,2)-sq3(2,1))^2+(sq3(2,2)-sq3(2,3))^2+(sq3(2,2)-sq3(3,2))^2)+eta*(sq3(2,2)-Y(i,j))^2/9;[~,index] = min([E1,E2,E3]);X(i,j) = X(i,j) + z(index);if temp~=X(i,j)tot=tot+1;endendendif tot<1break;end
endJ2=X;
subplot(2,2,3)
imagesc(J2);
colormap('gray')
title('mrf降噪结果')subplot(2,2,4)
imagesc(J2-J1);
colormap('gray')
title('差异')
imwrite(uint8(J2),'lena_denoised.jpg');

结果

原始图

加噪声图

降噪图

从结果看出,噪声得到一定程度的遏制,但是图像锐度发生下降,这是因为代价函数总是趋向于让图像变得光滑导致的。如果把更准确的先验知识应用代价函数中,则可以获得更理想的效果。
例如为了避免图像变化剧烈的正常边缘被模糊,将先验知识中相邻像素的势函数修改为
F g ( x i , j , x i − 1 , j , x i , j − 1 ) ≜ − 1 2 σ v k ∣ ( x i , j − x i − 1 , j ) ∣ k + ∣ ( x i , j − x i , j − 1 ) ∣ k F_g(x_{i,j},x_{i-1,j},x_{i,j-1}) \triangleq -\frac{1}{2\sigma_v^k} |(x_{i,j}-x_{i-1,j})|^k+|(x_{i,j}-x_{i,j-1})|^k Fg​(xi,j​,xi−1,j​,xi,j−1​)≜−2σvk​1​∣(xi,j​−xi−1,j​)∣k+∣(xi,j​−xi,j−1​)∣k
k k k越小,则表示越忽视过大差异的势,若边缘模糊,则说明代价函数过多关注了相邻像素的光滑性,需要将k减小,以提高图像的锐度,例如当 k = 1.5 k=1.5 k=1.5时,图像锐度得到了提高:


上述三张图分别时含噪图, k = 2 k=2 k=2 和 k = 1.5 k=1.5 k=1.5,可以看出 k k k的下降有利于保护原图中高反差的边缘(例如胳膊和背景,帽子亮部和背景),但一些低反差的细节可能会损失(例如发丝细节)。实际中可以通过机器学习的方式获得最佳的参数估计。

Bayes,HMM,MRF Gibbs Distribution在图像降噪中的应用相关推荐

  1. 【知识星球】图像降噪模型和数据集内容开启更新,经典问题永垂不朽!

    欢迎大家来到<知识星球>专栏,这里是网络结构1000变小专题,今天介绍的是我们知识星球图像降噪模型和数据集相关专题上线. 作者&编辑 | 言有三 1 图像降噪问题 图像降噪是深度学 ...

  2. 可复现的图像降噪算法总结——超赞整理

    点击上方"小白学视觉",选择加"星标"或"置顶" 重磅干货,第一时间送达 本文转自:AI算法与图像处理 图像降噪,是最简单也是最基础的图像处 ...

  3. 视觉进阶 | 用于图像降噪的卷积自编码器

    点击上方"小白学视觉",选择加"星标"或"置顶" 重磅干货,第一时间送达 本文转自:磐创AI 作者|Dataman 编译|Arno 来源|A ...

  4. LIVE 预告 | CVPR 2021 预讲 · 旷视专场,覆盖目标检测、蒸馏、图像降噪、人体姿态估计等...

    CVPR 2021 大会将于6月19日至 25日线上举行.为促进同行之间的交流与合作,智源社区近期举办了系列CVPR 2021预讲报告,其中实验室系列将汇聚国内顶尖高校和企业实验室的研究人员为大家分享 ...

  5. 「每周CV论文」深度学习图像降噪应该阅读哪些文章

    图像降噪是图像处理领域中非常传统和经典的问题,今天给大家推荐学习该领域值得读的文章. 作者&编辑 | 言有三 1 基本CNN结构 图像去噪模型的输出是无噪声的图像,与输入图像大小相同,所以可以 ...

  6. 「技术综述」一文道尽传统图像降噪方法

    https://www.toutiao.com/a6713171323893318151/ 作者 | 黄小邪/言有三 编辑 | 黄小邪/言有三 图像预处理算法的好坏直接关系到后续图像处理的效果,如图像 ...

  7. 自适应图像降噪滤波器的设计与实现

    原文:http://blog.csdn.net/baimafujinji/article/details/73302911 有一个说法:好奇心和懒惰是推动人类发明创造的两大动力!所以今天我又再次展现了 ...

  8. 图像降噪算法——低秩聚类:WNNM算法

    图像降噪算法--低秩聚类:WNNM算法 图像降噪算法--低秩聚类:WNNM算法 1. 基本原理 2. matlab代码 3. 结论 图像降噪算法--低秩聚类:WNNM算法 同样是为了完善自己知识版图的 ...

  9. 图像降噪算法——稀疏表达:K-SVD算法

    图像降噪算法--稀疏表达:K-SVD算法 图像降噪算法--稀疏表达:K-SVD算法 1. 基本原理 2. python代码 3. 结论 图像降噪算法--稀疏表达:K-SVD算法 为了完善下自己降噪算法 ...

最新文章

  1. LA 3353 最优巴士线路设计
  2. 【转载】(EM算法)The EM Algorithm
  3. python核心编程第六章练习6-12
  4. c++ 返回智能指针_C++核心指南(17) I.11 禁止使用指针(T*)或引用(T)来转移所有权...
  5. CENTOS 7.0 安装discuz ,搭 mysql +php+apache 环境
  6. zoj1022 Parallel Expectations(DP)
  7. css常用选择器选择器
  8. 龙芯2k1000-pmon(5)- pmon无法修改环境变量的问题
  9. 数位板和sai2安装使用
  10. ET框架-03 ET框架-Demo工程的编译与运行
  11. jquery怎么读(jquery怎么读音英语)
  12. Java实现批量修改文件名
  13. qconshanghai2017
  14. 橙瓜码字多端同步、十份云储存本地实时备份,最放心的码字软件
  15. 关于系统前端开发的那些事
  16. PostgreSQL忘记密码
  17. Day11-软件测试设计之银行储蓄系统
  18. 【Python】random.randint()用法
  19. Vue全家桶系列之Vuex(一)
  20. pytorch: 给tensor删除或者添加维度为1的维度(squeeze和unsqueeze)

热门文章

  1. 世界主要货币英文缩写
  2. php 策略模式 理解
  3. 【重要】一部手机失窃而揭露的黑色产业链—完整修订版
  4. 【数据结构】顺序表的创建、插入、删除、合并
  5. 火爆的人工智能项目都在这里了|Gitee项目推荐
  6. 手机如何借用笔记本网络上网
  7. 计算机公式与函数试题,计算机国考样题EXCEL之公式与函数的应用一
  8. Android开发——如何设计开发一款Android App
  9. redis中的AKF理论和CAP理论详解
  10. opencv 学习笔记(七) 灰度变换