Example of DenseCRF with non-RGB data

本笔记本通过一个示例说明如何在非rgb数据上使用DenseCRFs。同时，它将解释基本概念并通过一个示例进行演示，因此即使您正在处理RGB数据，它也可能是有用的，不过也请查看PyDenseCRF's README !

在jupyter notebook上运行https://github.com/lucasb-eyer/pydensecrf/tree/master/examples/Non RGB Example.ipynb

Basic setup

先导入：

import pydensecrf.densecrf as dcrf
from pydensecrf.utils import unary_from_softmax, create_pairwise_bilateral

设置画图风格：

import numpy as np
import matplotlib.pyplot as plt

%matplotlib inline

plt.rcParams['image.interpolation'] = 'nearest' #设置插入风格
plt.rcParams['image.cmap'] = 'gray' #设置颜色风格

Unary Potential一元势

一元势由每个像素的类概率组成。这可以来自任何类型的模型，如随机森林或深度神经网络的softmax。

1）Create unary potential

补充：scipy.stats.multivariate_normal的使用

1.得到x,y的取值

from scipy.stats import multivariate_normal

H, W, NLABELS = , , 

# This creates a gaussian blob...
a = np.mgrid[:H, :W] #(,,) #得到x,y的取值
a,a.shape

返回：

(array([[[  ,   ,   , ...,   ,   ,   ],
         [  ,   ,   , ...,   ,   ,   ],
         [  ,   ,   , ...,   ,   ,   ],
         ...,
         [, , , ..., , , ],
         [, , , ..., , , ],
         [, , , ..., , , ]],

        [[  ,   ,   , ..., , , ],
         [  ,   ,   , ..., , , ],
         [  ,   ,   , ..., , , ],
         ...,
         [  ,   ,   , ..., , , ],
         [  ,   ,   , ..., , , ],
         [  ,   ,   , ..., , , ]]]), (, , ))

2.得到点(x,y)的值

pos = np.stack(a, axis=) #得到(x,y)点的值
pos,pos.shape

返回：

(array([[[  ,   ],
         [  ,   ],
         [  ,   ],
         ...,
         [  , ],
         [  , ],
         [  , ]],

        [[  ,   ],
         [  ,   ],
         [  ,   ],
         ...,
         [  , ],
         [  , ],
         [  , ]],

        [[  ,   ],
         [  ,   ],
         [  ,   ],
         ...,
         [  , ],
         [  , ],
         [  , ]],

        ...,

        [[,   ],
         [,   ],
         [,   ],
         ...,
         [, ],
         [, ],
         [, ]],

        [[,   ],
         [,   ],
         [,   ],
         ...,
         [, ],
         [, ],
         [, ]],

        [[,   ],
         [,   ],
         [,   ],
         ...,
         [, ],
         [, ],
         [, ]]]), (, , ))

3.得到概率值：

#mean=[, ], cov=,12800的单位乘作为协方差矩阵
rv = multivariate_normal([H//2, W//2], (H//4)*(W//4))
probs = rv.pdf(pos)
probs,probs.shape

返回：

(array([[2.01479637e-07, 2.05541767e-07, 2.09669414e-07, ...,
         2.13863243e-07, 2.09669414e-07, 2.05541767e-07],
        [2.04644486e-07, 2.08770423e-07, 2.12962908e-07, ...,
         2.17222613e-07, 2.12962908e-07, 2.08770423e-07],
        [2.07842810e-07, 2.12033230e-07, 2.16291237e-07, ...,
         2.20617517e-07, 2.16291237e-07, 2.12033230e-07],
        ...,
        [2.11074628e-07, 2.15330207e-07, 2.19654423e-07, ...,
         2.24047973e-07, 2.19654423e-07, 2.15330207e-07],
        [2.07842810e-07, 2.12033230e-07, 2.16291237e-07, ...,
         2.20617517e-07, 2.16291237e-07, 2.12033230e-07],
        [2.04644486e-07, 2.08770423e-07, 2.12962908e-07, ...,
         2.17222613e-07, 2.12962908e-07, 2.08770423e-07]]), (, ))

4.将概率值进行处理集中到[0.4,0.6]区间中

先将值集中到[0,1]区间中

# ...which we project into the range [0.4, 0.6]
probs = (probs-probs.min()) / (probs.max()-probs.min())
probs, probs.shape

返回：

(array([[.        , 0.00033208, 0.00066951, ..., 0.00101235, 0.00066951,
         0.00033208],
        [0.00025872, 0.00059602, 0.00093875, ..., 0.00128698, 0.00093875,
         0.00059602],
        [0.00052019, 0.00086275, 0.00121084, ..., 0.00156451, 0.00121084,
         0.00086275],
        ...,
        [0.00078439, 0.00113228, 0.00148578, ..., 0.00184495, 0.00148578,
         0.00113228],
        [0.00052019, 0.00086275, 0.00121084, ..., 0.00156451, 0.00121084,
         0.00086275],
        [0.00025872, 0.00059602, 0.00093875, ..., 0.00128698, 0.00093875,
         0.00059602]]), (, ))

再将其集中到[0.4,0.6]区间中，0.5 + 0.2 * (probs-0.5)等价于0.4+0.2*probs，而此时的probs为[0,1]区间中的值：

probs = 0.5 + 0.2 * (probs-0.5) #此时probs.shape为(, )
probs, probs.shape #将最终的概率结果的值集中在[0.4, 0.6]区间中

返回：

(array([[0.4       , 0.40006642, 0.4001339 , ..., 0.40020247, 0.4001339 ,
         0.40006642],
        [0.40005174, 0.4001192 , 0.40018775, ..., 0.4002574 , 0.40018775,
         0.4001192 ],
        [0.40010404, 0.40017255, 0.40024217, ..., 0.4003129 , 0.40024217,
         0.40017255],
        ...,
        [0.40015688, 0.40022646, 0.40029716, ..., 0.40036899, 0.40029716,
         0.40022646],
        [0.40010404, 0.40017255, 0.40024217, ..., 0.4003129 , 0.40024217,
         0.40017255],
        [0.40005174, 0.4001192 , 0.40018775, ..., 0.4002574 , 0.40018775,
         0.4001192 ]]), (, ))

5.一份数据已经准备好了，要复制生成一个效果相反的数据：

# 第一个维度需要等于类的数量，这里设置两个类，为2
# 让我们有一个"foreground"和一个"background"类
#因此，复制高斯blob，但倒置它去创建与“background”类的概率是相反的“foreground”类。
#np.tile(a,(m,n)):即是把a数组里面的元素复制n次放进一个数组c中，然后再把数组c复制m次放进一个数组b中
#np.tile(a,(,)) #将a数组重复3次形成2行的数组
#np.newaxis为多维数组增加一个轴
#probs[np.newaxis,:,:].shape为(, , )
probs = np.tile(probs[np.newaxis,:,:],(,,))#此时probs.shape为(, , )
probs[,:,:] =  - probs[,:,:] #得到相反的值
probs,probs.shape

返回：

(array([[[0.4       , 0.40006642, 0.4001339 , ..., 0.40020247,
          0.4001339 , 0.40006642],
         [0.40005174, 0.4001192 , 0.40018775, ..., 0.4002574 ,
          0.40018775, 0.4001192 ],
         [0.40010404, 0.40017255, 0.40024217, ..., 0.4003129 ,
          0.40024217, 0.40017255],
         ...,
         [0.40015688, 0.40022646, 0.40029716, ..., 0.40036899,
          0.40029716, 0.40022646],
         [0.40010404, 0.40017255, 0.40024217, ..., 0.4003129 ,
          0.40024217, 0.40017255],
         [0.40005174, 0.4001192 , 0.40018775, ..., 0.4002574 ,
          0.40018775, 0.4001192 ]],

        [[0.6       , 0.59993358, 0.5998661 , ..., 0.59979753,
          0.5998661 , 0.59993358],
         [0.59994826, 0.5998808 , 0.59981225, ..., 0.5997426 ,
          0.59981225, 0.5998808 ],
         [0.59989596, 0.59982745, 0.59975783, ..., 0.5996871 ,
          0.59975783, 0.59982745],
         ...,
         [0.59984312, 0.59977354, 0.59970284, ..., 0.59963101,
          0.59970284, 0.59977354],
         [0.59989596, 0.59982745, 0.59975783, ..., 0.5996871 ,
          0.59975783, 0.59982745],
         [0.59994826, 0.5998808 , 0.59981225, ..., 0.5997426 ,
          0.59981225, 0.5998808 ]]]), (, , ))

6.最终的代码为：

from scipy.stats import multivariate_normal

H, W, NLABELS = , , 

# This creates a gaussian blob...
pos = np.stack(np.mgrid[:H, :W], axis=)

rv = multivariate_normal([H//2, W//2], (H//4)*(W//4))
probs = rv.pdf(pos)

# ...which we project into the range [0.4, 0.6]
probs = (probs-probs.min()) / (probs.max()-probs.min())
probs = 0.5 + 0.2 * (probs-0.5) #此时probs.shape为(, )

# 第一个维度需要等于类的数量
# 让我们有一个"foreground"和一个"background"类
#因此，复制高斯blob，但倒置它去创建与“background”类的概率是相反的“foreground”类。
#np.tile(a,(m,n)):即是把a数组里面的元素复制n次放进一个数组c中，然后再把数组c复制m次放进一个数组b中
#np.tile(a,(,)) #将a数组重复3次形成2行的数组
#np.newaxis为多维数组增加一个轴
#probs[np.newaxis,:,:].shape为(, , )
probs = np.tile(probs[np.newaxis,:,:],(,,))#此时probs.shape为(, , )
probs[,:,:] =  - probs[,:,:] #得到相反的值

# Let's have a look:
plt.figure(figsize=(,))
plt.subplot(,,); plt.imshow(probs[,:,:]); plt.title('Foreground probability'); plt.axis('off'); plt.colorbar();
plt.subplot(,,); plt.imshow(probs[,:,:]); plt.title('Background probability'); plt.axis('off'); plt.colorbar();

返回：

2）Run inference with unary potential 使用一元势运行推理

我们已经可以运行一个只有一元势的DenseCRF。这不是一个好主意，但我们可以做到:

1.unary_from_softmax函数的作用

将softmax类概率转换为一元势(每个节点的NLL)。

即我们之前先对图片使用训练好的网络预测得到最终经过softmax函数得到的分类结果，这里需要将这个结果转成一元势
参数

• sm: numpy.array ,第一个维度是类的softmax的输出，其他所有维度都是flattend。这意味着“sm.shape[0] == n_classes”。
• scale: float,softmax输出的确定性(默认为None)，需要值在(0,1]。如果不为None，则softmax输出被缩放到从[0,scale]概率的范围。
• clip: float,将概率裁剪到的最小值。这是因为一元函数是概率的负对数，而log(0) = inf，所以我们需要把0概率裁剪成正的值。

在这里因为scale=None,clip=None,所以这个函数的作用其实只进行了下面的操作：

-np.log(sm).reshape([num_cls, -]).astype(np.float32)

对值取负对数，并将结果变为(2, 204800)大小，进行了flatten操作，将图片压扁了：

# Inference without pair-wise terms
#设置固定的一元势
from pydensecrf.utils import unary_from_softmax
U = unary_from_softmax(probs)  # note: num classes is first dim，类别在第一维
U,U.shape

返回：

(array([[0.91629076, 0.9161247 , 0.915956  , ..., 0.91564745, 0.9158215 ,
         0.9159928 ],
        [0.51082563, 0.5109363 , 0.5110488 , ..., 0.5112547 , 0.5111386 ,
         0.5110243 ]], dtype=float32), (, ))

此时的值就是一元势了

d = dcrf.DenseCRF2D(W, H, NLABELS)
d.setUnaryEnergy(U)

# Run inference for  iterations，10次迭代
Q_unary = d.inference()

# Q现在是近似后验，我们可以用argmax得到一个MAP估计值。
#np.argmax得到Q_unary列方向上的最大值的索引
map_soln_unary = np.argmax(Q_unary, axis=)#此时形式为(,)
map_soln_unary,map_soln_unary.shape

返回：

(array([, , , ..., , , ]), (,))

# 不幸的是，DenseCRF把所有东西都压扁了，所以把它恢复到图片形式。
map_soln_unary = map_soln_unary.reshape((H,W)) #重新变为(, ))
map_soln_unary,map_soln_unary.shape

返回：

(array([[, , , ..., , , ],
        [, , , ..., , , ],
        [, , , ..., , , ],
        ...,
        [, , , ..., , , ],
        [, , , ..., , , ],
        [, , , ..., , , ]]), (, ))

总的代码为：

# Inference without pair-wise terms
U = unary_from_softmax(probs)  # note: num classes is first dim
d = dcrf.DenseCRF2D(W, H, NLABELS)
d.setUnaryEnergy(U)

# Run inference for  iterations，10次迭代
Q_unary = d.inference()

# Q现在是近似后验，我们可以用argmax得到一个MAP估计值。
#np.argmax得到Q_unary列方向上的最大值的索引
map_soln_unary = np.argmax(Q_unary, axis=)#此时形式为(,)

# 不幸的是，DenseCRF把所有东西都压扁了，所以把它恢复到图片形式。
map_soln_unary = map_soln_unary.reshape((H,W))#重新变为(, ))

# And let's have a look.
plt.imshow(map_soln_unary); plt.axis('off'); plt.title('MAP Solution without pairwise terms');

返回：

因为设置了颜色为gray，实际结果为：

Pairwise terms二元

DenseCRFs的全部意义在于使用某种形式的内容来平滑预测。这是通过“成对”术语来实现的，“成对”术语编码元素之间的关系。

1）Add (non-RGB) pairwise term

例如，在图像处理中，一个流行的成对关系是“双边”(bilateral)关系，它大致表示颜色相似或位置相似的像素可能属于同一个类。

NCHAN= #通道为1，non-RGB

#创造简单的双边图象。
#注意，我们将通道维度放在最后，但我们也可以将它作为第一个维度，只要将chdim参数进一步更改为0。
img = np.zeros((H,W,NCHAN), np.uint8)
img[H//3:2*H//3,W//4:3*W//4,:] = 1 #使得[133:266,128:384,:]内的值为1，下面画为白色

plt.imshow(img[:,:,]); plt.title('Bilateral image'); plt.axis('off'); plt.colorbar();

返回：

2.create_pairwise_bilateral函数的作用

创造成对双边势的Util函数。这适用于所有的图像尺寸。对于2D情况，与“densecrf2 . addpairwisebilateral”函数效果相同。

参数：

• sdims: list or tuple，每个维度的比例因子。这在DenseCRF2D.addPairwiseBilateral中称为“sxy”。
• schan: list or tuple，图像中每个通道的比例因子。这在DenseCRF2D.addPairwiseBilateral中称为“srgb”。
• img: numpy.array，输入的图像
• chdim: int, 可选。这指定了通道维度在图像中的位置。例如，' chdim=2 '用于大小为(240,300,3)的RGB图像，说明其通道3值放在维度2上(维度从0开始)。如果图像没有通道尺寸(例如只有一个通道)，则使用' chdim=-1 '。

#从上图中创建成对的双边术语。
#这两个' s{dims,chan} '参数是模型超参数，分别定义了位置的强度和图像内容的双边参数。
#将通道为1的图变为了通道为3的图
pairwise_energy = create_pairwise_bilateral(sdims=(,), schan=(0.01,), img=img, chdim=) #其shape为(, )

#pairwise_energy现在包含的维度和DenseCRF的特性一样多，在本例中是3:(x,y,channel1)
img_en = pairwise_energy.reshape((-, H, W))  # Reshape just for plotting,为了画图

#下面将三个通道对应的图画出来
plt.figure(figsize=(,))
plt.subplot(,,); plt.imshow(img_en[]); plt.title('Pairwise bilateral [x]'); plt.axis('off'); plt.colorbar();
plt.subplot(,,); plt.imshow(img_en[]); plt.title('Pairwise bilateral [y]'); plt.axis('off'); plt.colorbar();
plt.subplot(,,); plt.imshow(img_en[]); plt.title('Pairwise bilateral [c]'); plt.axis('off'); plt.colorbar();

返回：

2）Run inference of complete DenseCRF 

现在我们可以创建一个包含一元势和成对势的dense CRF ，并对其进行推理以得到最终结果。

d = dcrf.DenseCRF2D(W, H, NLABELS)
d.setUnaryEnergy(U) #设置一元势
#添加二元势
d.addPairwiseEnergy(pairwise_energy, compat=)  # `compat` is the "strength" of this potential.

#这一次，让我们自己分步骤进行推理，这样我们就可以查看中间解，并监视KL-divergence，它表明我们收敛得有多好。
#PyDenseCRF还要求我们跟踪计算所需的两个临时缓冲区，tmp1, tmp2。
Q, tmp1, tmp2 = d.startInference()
for _ in range(): #迭代5次后，查看结果
    d.stepInference(Q, tmp1, tmp2)
kl1 = d.klDivergence(Q) / (H*W)
map_soln1 = np.argmax(Q, axis=).reshape((H,W))

for _ in range():
    d.stepInference(Q, tmp1, tmp2)
kl2 = d.klDivergence(Q) / (H*W)
map_soln2 = np.argmax(Q, axis=).reshape((H,W))

for _ in range():
    d.stepInference(Q, tmp1, tmp2)
kl3 = d.klDivergence(Q) / (H*W)
map_soln3 = np.argmax(Q, axis=).reshape((H,W))

img_en = pairwise_energy.reshape((-, H, W))  # Reshape just for plotting
plt.figure(figsize=(,))
plt.subplot(,,); plt.imshow(map_soln1);
plt.title('MAP Solution with DenseCRF\n(5 steps, KL={:.2f})'.format(kl1)); plt.axis('off');
plt.subplot(,,); plt.imshow(map_soln2);
plt.title('MAP Solution with DenseCRF\n(20 steps, KL={:.2f})'.format(kl2)); plt.axis('off');
plt.subplot(,,); plt.imshow(map_soln3);
plt.title('MAP Solution with DenseCRF\n(50 steps, KL={:.2f})'.format(kl3)); plt.axis('off');

返回：

这个不太好看出效果，下面：

可见50次迭代后效果很好