图像中每个像素$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 )$$

编程实现

In [1]:
import cv2
import numpy as np
import matplotlib.pyplot as plt
from random import randint
from math import *
In [2]:
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

用于将类别转换为色值

In [3]:
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
In [4]:
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

方差等于平方的均值减去均值的平方

In [5]:
def variance(sums1, squares1, nos1):
    return squares1 / nos1 - (sums1 / nos1) ** 2
In [6]:
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)
In [7]:
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

读取图片

In [8]:
original = cv2.imread(imagepath)
origflt = original.astype(float)

转换成灰度图片

In [9]:
img = cv2.cvtColor(original, cv2.COLOR_BGR2GRAY)
plt.imshow(img, cmap='gray')
Out[9]:
<matplotlib.image.AxesImage at 0x7fdc40fc6d68>

随机初始化每个像素点的类别,然后计算下列参数

In [10]:
sums, squares, nos, labels = initialize(img)
variances = [variance(sums[i], squares[i], nos[i]) for i in range(SEGS)]
energy = calculateEnergy(img, variances, labels)
In [11]:
variances, sums, squares, nos, energy
Out[11]:
([6564.583863855842, 6574.225853623524],
 [11399606.0, 11363030.0],
 [2251465872.0, 2239780624.0],
 [73447.0, 73508.0],
 645816.3449991856)

使用模拟退火算法让能量最低

In [12]:
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
In [13]:
plt.imshow(reconstruct(labels), interpolation='nearest', cmap='gray')
Out[13]:
<matplotlib.image.AxesImage at 0x7fdc406fcd30>
In [ ]:
 
posted @ 2019-01-19 10:54:57
评论加载中...

发表评论