图像中每个像素$s$都用一个特征向量$\vec { f } _ { s }$表示: $$f = \left\{ \vec { f } _ { s } : s \in S \right\}$$
每一个像素都对应一个类别: $$\omega = \left\{ \omega _ { s } , s \in S \right\}$$
于是,图像分割相当于找到每个像素点对应的类别,即给定特征向量$f$时,求概率最大的$\omega$: $$\hat { \omega } ^ { M A P } = \arg \max _ { \omega \in \Omega } P ( \omega | f )$$
根据贝叶斯定理: $$P ( \omega | f ) = \frac { P ( f | \omega ) P ( \omega ) } { P ( f ) }\propto P ( f | \omega ) P ( \omega )$$
其中$P(\omega)$是先验概率,$P (f | \omega)$是条件概率.
马尔科夫状态的先验概率服从吉布斯分布: $$P ( \omega ) = \frac { 1 } { Z } \exp ( - U ( \omega ) ) = \frac { 1 } { Z } \exp \left( - \sum _ { c \in C } V _ { c } ( \omega ) \right)$$
$$\mathrm { U } ( \omega ) = \sum _ { \mathrm { c } \in \mathrm { C } } \mathrm { V } _ { \mathrm { c } } ( \omega ) = \sum _ { \mathrm { i } \in \mathrm { C } _ { 1 } } \mathrm { V } _ { \mathrm { C } _ { 1 } } \left( \omega _ { \mathrm { i } } \right) + \sum _ { ( \mathrm { i } , \mathrm { j } ) \in \mathrm { C } _ { 2 } } \mathrm { V } _ { \mathrm { C } _ { 2 } } \left( \omega _ { \mathrm { i } } , \omega _ { \mathrm { j } } \right) + \ldots$$假设已知像素点的类别$\omega$,色值$f$会服从高斯分布,即:
$$P \left( f _ { s } | \omega _ { s } \right) = \frac { 1 } { \sqrt { 2 \pi } \sigma _ { \omega _ { s } } } \exp \left( - \frac { \left( f _ { s } - \mu _ { \omega _ { s } } \right) ^ { 2 } } { 2 \sigma _ { \omega _ { s } } ^ { 2 } } \right)$$可以证明,后验概率还是吉布斯分布: $$P ( \omega | f ) = P ( \omega )P \left( f _ { s } | \omega _ { s } \right) =\frac { 1 } { Z } \exp ( - U ( \omega ) )$$
$$U ( \omega ) = \sum _ { s } \left( \log \left( \sqrt { 2 \pi } \sigma _ { \omega _ { s } } \right) + \frac { \left( f _ { s } - \mu _ { \omega _ { s } } \right) ^ { 2 } } { 2 \sigma _ { \omega _ { s } } ^ { 2 } } \right) + \sum _ { s , r } \beta \delta \left( \omega _ { s } , \omega _ { r } \right)$$马尔科夫模型指一件事物的当前状态只与它之前的1个或者n个状态有关，而与再之前的状态没有关系，比如今天天气好坏只与昨天天气有关，而与前天乃至大前天都没有关系。符合这样的一种特性的事物认为其具有马尔科夫性。那么引申到图像领域，就是认为图像中某一点的特征（一般都是像素点灰色、颜色值等）只与其附近的一小块领域有关，而与其他的领域无关。想想也是这样，大多数时候，图像中某一点像素和附近的像素是相关的，附近领域像素是黑的，那它八成也是黑的，很好理解。根据这一性质,像素点仅与邻居相关:
对于一个像素点属于某个类别的概率$P(\omega | f )$,可以理解成像素本身的色值与周围像素类别共同决定的,其中像素本身的色值相当于$U ( \omega )$中的: $$\sum _ { s } \left( \log \left( \sqrt { 2 \pi } \sigma _ { \omega _ { s } } \right) + \frac { \left( f _ { s } - \mu _ { \omega _ { s } } \right) ^ { 2 } } { 2 \sigma _ { \omega _ { s } } ^ { 2 } } \right)$$ 而周围像素类别相当于$U ( \omega )$中的: $$\sum _ { s , r } \beta \delta \left( \omega _ { s } , \omega _ { r } \right)$$ 其中: $$\beta \delta \left( \omega _ { i } , \omega _ { j } \right) = \left\{ \begin{array} { l l } { - \beta } & { \text { if } } & { \omega _ { i } = \omega _ { j } } \\ { + \beta } & { \text { if } } & { \omega _ { i } \neq \omega _ { j } } \end{array} \right.$$
于是最大化后验概率相当于最小化$U ( \omega )$,在物理学中$U ( \omega )$代表着能量:
$$\hat { \omega } ^ { M A P } = \arg \max _ { \omega \in \Omega } P ( \omega | f ) = \arg \min _ { \omega \in \Omega } U ( \omega )$$编程实现¶
import cv2
import numpy as np
import matplotlib.pyplot as plt
from random import randint
from math import *
def isSafe(shape, x, y):
a, b = shape
return x >= 0 and x < a and y >= 0 and y < b
def delta(i, l):
if (i == l):
return -BETA
return BETA
用于将类别转换为色值
def reconstruct(labs):
labels = labs
for i in range(len(labels)):
for j in range(len(labels[0])):
labels[i][j] = (labels[i][j] * 255) / (SEGS - 1)
return labels
def calculateEnergy(img, variances, labels):
energy = 0.0
for i in range(len(img)):
for j in range(len(img[0])):
l = labels[i][j]
energy += log(sqrt(variances[l]))
for (p, q) in NEIGHBORS:
if isSafe(img.shape, i + p, j + q):
energy += (delta(l, labels[i + p][j + q]) / 2.0)
return energy
方差等于平方的均值减去均值的平方
def variance(sums1, squares1, nos1):
return squares1 / nos1 - (sums1 / nos1) ** 2
def initialize(img):
labels = np.zeros(shape=img.shape, dtype=np.uint8)
nos = [0.0] * SEGS
sums = [0.0] * SEGS
squares = [0.0] * SEGS
for i in range(len(img)):
for j in range(len(img[0])):
l = randint(0, SEGS - 1)
sums[l] += img[i][j]
squares[l] += img[i][j] ** 2
nos[l] += 1.0
labels[i][j] = l
return (sums, squares, nos, labels)
imagepath = 'images/scar.jpg'
SEGS = 2 # 切分成多少部分
NEIGHBORS = [(-1, 0), (1, 0), (0, -1), (0, 1)]
BETA = 1
temperature = 4.0
iterations = 1000000
COOLRATE = 0.95
读取图片
original = cv2.imread(imagepath)
origflt = original.astype(float)
转换成灰度图片
随机初始化每个像素点的类别,然后计算下列参数
sums, squares, nos, labels = initialize(img)
variances = [variance(sums[i], squares[i], nos[i]) for i in range(SEGS)]
energy = calculateEnergy(img, variances, labels)
variances, sums, squares, nos, energy
使用模拟退火算法让能量最低
while iterations > 0:
(a, b) = img.shape
change = False
# 随机取一个点
x = randint(0, a - 1)
y = randint(0, b - 1)
# 取出这个点的值,和类别
val = float(img[x][y])
l = labels[x][y]
# 随机取一个label
newl = l
while newl == l:
newl = randint(0, SEGS - 1)
# 计算方差,将原来的类别减去,增加到新的类别
remsums = sums[l] - val
addsums = sums[newl] + val
remsquares = squares[l] - val * val
addsquares = squares[newl] + val * val
remvar = variance(remsums, remsquares, nos[l] - 1)
addvar = variance(addsums, addsquares, nos[newl] + 1)
# 计算色值对类别的影响
newenergy = energy
# 更新类别l的能量
# energy += np.log(np.sqrt(2*np.pi*var))
# energy += ((pixels[i][j]-mean)**2)/(2*var)
newenergy -= log(sqrt(variance(sums[l], squares[l], nos[l]))) * (nos[l])
newenergy += log(sqrt(remvar)) * (nos[l] - 1)
# 更新类别newl的能量
newenergy -= log(sqrt(variance(sums[newl], squares[newl], nos[newl]))) * (nos[newl])
newenergy += log(sqrt(addvar)) * (nos[newl] + 1)
# 如果和周围的点类别相同,减少能量,否则增加能量
for (p, q) in NEIGHBORS:
if isSafe((a, b), x + p, y + q):
newenergy -= delta(l, labels[x + p][y + q])
newenergy += delta(newl, labels[x + p][y + q])
# 如果能量减小则修改类别,否则以一定概率接受修改
if newenergy < energy:
change = True
else:
prob = 1.1
if temperature != 0:
prob = np.exp((energy - newenergy) / temperature)
if prob >= (randint(0, 1000) + 0.0) / 1000:
change = True
if change:
labels[x][y] = newl
energy = newenergy
nos[l] -= 1
sums[l] = remsums
squares[l] = remsquares
nos[newl] += 1
sums[newl] = addsums
squares[newl] = addsquares
temperature *= COOLRATE
iterations -= 1