Использование цепей Маркова Монте-Карло в сочетании с моделью Изинга для очистки зашумленного двоичного изображения

В этой статье я продемонстрирую использование цепей Маркова Монте-Карло для удаления шума из двоичного изображения.

Марковская цепь Монте-Карло или сокращенно MCMC относится к классу методов, используемых для оценки распределения вероятностей путем выборки из него. Различные методы, составляющие MCMC, отличаются друг от друга на основе метода, используемого для отбора образцов. Некоторые из наиболее известных методов MCMC - это Метрополис - Гастингс, выборка Гиббса и гамильтониан Монте-Карло. Техника, которую я буду использовать, - это выборка Гиббса.

Выборка Гиббса - это метод выборки из многомерного распределения, при условии, что все остальные переменные остаются неизменными.

Например, если есть распределение только с двумя переменными x1 и x2, выборка будет следующей:

Он начинается с точки 1, а затем, когда необходимо выполнить выборку второй точки, он выглядит следующим образом:

P (x2 | x1), он ищет следующий образец в той же строке, что и x1 (сохраняя его постоянным)

После чего выбирается третья точка следующим образом:

P (x1 | x2), он ищет следующий образец в той же строке, что и x2 (сохраняя его постоянным)

И так далее, процесс выборки продолжается для заданного количества точек, позволяя ему пересечь все пространство.

Модель Изинга - это математическая модель, соответствующая квадратной решетке, используемой для моделирования фазовых переходов. Каждый элемент решетки может существовать в двух дискретных состояниях и может быть представлен как +1 и -1. Каждый элемент оказывает влияние на все свои соседние элементы и пытается достичь состояния равновесия, когда все элементы находятся в одном и том же состоянии.

Приложение:

Бинарное изображение можно представить в виде решетки, где каждый пиксель представляет один элемент. Состояние пикселя может быть представлено как 1 или -1 в зависимости от цвета пикселя. Изображение можно представить состоящим из двух слоев: нижележащий слой представляет истинное изображение, а верхний слой представляет шум. Говорят, что гауссов шум накладывается на изображение, причем в некоторых местах он совпадает с реальным изображением, а в некоторых местах принимает противоположные значения.

К вышеупомянутому слою, состоящему из шума, применяется модель Изинга. На шум влияет каждый из его соседей в зависимости от того, насколько тесно они связаны друг с другом, представленного краевым потенциалом. Чем ближе они связаны друг с другом, тем больше стараются находиться в одном состоянии. Формула для краевого потенциала имеет вид:

exp (J, X ₐ , X ₙ )

Здесь J - константа связи, которая обозначает, насколько тесно соседи связаны друг с другом. Xₐ представляет рассматриваемый пиксель, а Xₙ представляет наблюдаемые значения его соседей.

Модель наблюдения Гаусса используется для моделирования связи между шумом и ранее существовавшими пикселями. Считается, что шум является функцией фактического значения нижележащего пикселя и стандартного отклонения от него. Он может быть представлен:

N(Yₐ | Xₐ, σ ²)

Здесь Yₐ представляет собой наблюдаемое значение пикселя, Xₐ - рассматриваемый пиксель, а σ² - стандартное отклонение.

Именно взаимодействие между этими двумя силами определяет конечную ценность пикселя. Константа связи пытается сохранить все соседние пиксели (шум) в одном и том же состоянии или цвете. Гауссовская модель пытается сопоставить фактический цвет пикселя с наблюдаемым цветом (шумом). Воздействием этих двух сил можно управлять, изменяя значения J и σ².

Применяемая выборка Гиббса будет производить выборку каждого пикселя, обусловленного всеми его соседями и лежащим в основе истинным значением. Затем он исправит это значение, перейдет к следующему элементу и повторит ту же операцию. Одна итерация будет завершена, когда она завершит обход всей решетки. В зависимости от количества итераций конечное качество изображения может отличаться.

Xₐ рассматриваемый пиксель может принимать значения +1 или -1. Вышеприведенное выражение дает вероятность того, что Xₐ примет значения +1 или -1.

Числитель дает вероятность того, что Xₐ будет, например, +1, и проверяет его связь с наблюдаемым значением пикселя Yₐ и с наблюдаемыми значениями его соседей Xₙ.

Оно делится на сумму вероятности того, что каждое из них принимает значения +1 и -1, чтобы получить вероятность того, что фактическое значение Xₐ будет +1.

Код:

Цель кода - восстановить исходное изображение из поврежденного.

import numpy as np
import cv2 
import random 
import scipy
from scipy.spatial import distance
from scipy.stats import multivariate_normal
import pandas as pd
from PIL import Image
data = Image.open('noise_img.png')
image = np.asarray(data).astype(np.float32)

Изображение было импортировано и сохранено в виде двумерного массива со значениями оттенков серого пикселей. Это изображение было преобразовано в массив для удобства работы.

Затем он был преобразован в модель ising путем замены всех 0 (соответствующих черному) на -1 и всех 255 (соответствующих белому) на +1. Массив был дополнен нулями со всех сторон, чтобы упростить итерацию по решетке.

#Convert image values to ising model
for i in range(len(image)):
    for j in range(len(image[0])):
        if image[i,j,:] == 255:
            image[i,j,:] = 1
        else:
            image[i,j,:] = -1
#Create array to perform operations on
ising = np.zeros((len(image)+2,len(image[0])+2))
for i in range(len(image)):
    for j in range(len(image[0])):
        ising[i+1,j+1] = image[i,j,:]

Первая строка и столбец - это отступы. Остальные строки и столбец соответствуют цвету в один пиксель каждая. Номер строки и номер столбца соответствуют положению пикселя, а значение ячейки соответствует цвету.

Перед началом отбора образцов по Гиббсу было установлено значение силы сцепления.

#Coupling strength
J=4
#Gibbs sampling 
for n in range(3):
    for i in range(1,len(ising[0])-1):
        for j in range(1,len(ising)-1):
            pot = []
            for x in [-1, 1]:
                edge_pot = np.exp(J*ising[j-1,i]*x) * np.exp(J*ising[j,i-1]*x) * np.exp(J*ising[j+1,i]*x) * np.exp(J*ising[j,i+1]*x)
                pot.append(edge_pot)
            prob1 = multivariate_normal.pdf(image[j-1,i-1,:], mean = 1, cov = 1)*pot[1]/(multivariate_normal.pdf(image[j-1,i-1,:], mean = 1, cov = 1)*pot[1] + multivariate_normal.pdf(image[j-1,i-1,:], mean = -1, cov = 1)*pot[0]) 
            if np.random.uniform() <= prob1:
                ising[j,i] = 1
            else:
                ising[j,i] = -1

Все значения в исходном массиве, за исключением нулей, повторяются. Каждая точка выбирается путем проверки вероятности того, что исходный цвет равен +1 или -1. Для обоих случаев рассчитывается краевой потенциал с соседями. Наряду с краевым потенциалом вероятность наблюдаемого значения вычисляется относительно истинного значения, равного +1 или -1.

На основе двух вышеупомянутых значений вероятность того, что пиксель будет иметь значение +1, вычисляется путем деления его на вероятность того, что пиксель равен +1 и -1.

Затем это значение сравнивается со значением, случайно полученным из стандартного нормального распределения. Если вероятность того, что исходное значение будет +1, выше, тогда значение будет установлено на +1, иначе оно будет установлено на -1.

Когда выбирается следующий элемент, он принимает новое значение указанного выше элемента при вычислении краевого потенциала. Этот код проходит через каждый элемент всей решетки трижды.

#Retrieving the final array
final = np.zeros((len(image),len(image[0])))
final = ising[1:len(ising)-1,1:len(ising[0])-1]
#Converting it back to image
for i in range(len(final[0])):
    for j in range(len(final)):
        if final[j,i] == 1:
            final[j,i] = 255
        else:
            final[j,i] = 0

Затем отступы удаляются, а затем значения преобразуются обратно в значения в градациях серого, а затем визуализируются.

Наконец, у нас есть изображение с шумоподавлением.

Вы также можете связаться со мной в LinkedIn.