PRML Chapitre 8 Algorithme Somme des produits Implémentation Python
Le chapitre 8 décrit le modèle graphique. Un modèle graphique est une méthode d'expression graphique de relations telles que des variables stochastiques et des paramètres de modèle. Les modèles tels que la régression linéaire dans PRML peuvent généralement être représentés à l'aide de modèles graphiques. J'ai implémenté un algorithme de somme de produits qui peut être appliqué à un modèle qui peut être représenté par un modèle graphique et j'ai essayé de supprimer le bruit de l'image. Cette méthode semble être une forme générale de la méthode appliquée au modèle de Markov caché introduit au chapitre 13 de PRML.
Suppression du bruit de l'image
Ici, considérons une image binaire (blanc: 1, noir: -1) dans laquelle les pixels de l'image ont soit une valeur de -1,1.
Étant donné une image avec des valeurs inversées à certains pixels, comme ci-dessous
Restaure une image sans bruit qui ressemble à ceci:
Champ de probabilité de Markov
Comme modèle pour supprimer le bruit d'image, nous utilisons le champ stochastique de Markov (réseau de Markov, modèle graphique non orienté), qui est une sorte de modèle graphique.
La figure ci-dessous montre schématiquement le champ de probabilité de Markov considéré cette fois.
Le nœud ** $ y $ est une variable stochastique qui représente la valeur de pixel de l'image que nous avons observée avec du bruit, et le nœud $ x $ est la variable stochastique qui représente la valeur de pixel de l'image sans bruit que nous voulons restaurer **. Dans une image sans bruit, les pixels adjacents sont considérés comme ayant une forte corrélation, donc même dans ce modèle, les nœuds adjacents sont connectés par une ligne en $ x $. De plus, étant donné que la proportion de bruit est relativement petite, il existe une forte corrélation entre la valeur d'un pixel dans une image sans bruit et la valeur de pixel correspondante dans une image bruyante.
Par exemple, si un pixel dans une image bruyante a une valeur de 1, alors les pixels adjacents à ce pixel sont susceptibles d'être 1, et la valeur de pixel correspondante dans une image bruyante peut également être 1. haute. Cela signifie que les nœuds connectés par ** lignes ont tendance à avoir la même valeur de pixel **.
Exprimant cela comme une formule,
p({\bf x},{\bf y}) = {1\over Z}\exp\left\{-E({\bf x},{\bf y})\right\}
Cependant, en utilisant $ i, j $ comme index pour représenter le nœud,
E({\bf x},{\bf y}) = -\alpha\sum_ix_iy_i-\beta\sum_{Adjacent i et j}x_ix_j
Et $ Z $ est une constante de normalisation.
Étant donné une image avec du bruit
Algorithme de somme de produits
Cet algorithme est une méthode de calcul de la probabilité de la valeur qu'un nœud de variable stochastique sur un modèle graphique peut prendre. La probabilité est calculée en échangeant quelque chose comme une probabilité non normalisée appelée un message des nœuds connectés par une ligne.
Message de $ y_i $ à $ x_i $
m_{y_i\to x_i} = \exp(\alpha x_i y_i)
Cependant, $ y_i $ ici n'est pas une variable stochastique mais une valeur observée.
Message de $ x_j $ à $ x_i $
m_{x_j\to x_i} = \sum_{x_j}\exp(\beta x_ix_j)f(x_j)
Cependant, $ f (x_j) $ est le produit de messages provenant de nœuds autres que le nœud $ i $ adjacent au nœud $ j $. Pour calculer $ f (x_j) $, vous avez besoin d'un message d'un autre nœud au nœud $ i $, mais comme le modèle graphique que nous traitons cette fois-ci est un modèle avec une boucle, initialisez d'abord le message avec une certaine valeur, puis Envoyez les deux messages ci-dessus pour estimer $ {\ bf x} $.
la mise en oeuvre
sum_product.py
import itertools
import numpy as np
from scipy.misc import imread, imsave
ORIGINAL_IMAGE = "qiita_binary.png "
NOISY_IMAGE = "qiita_noise.png "
DENOISED_IMAGE = "qiita_denoised.png "
class Node(object):
def __init__(self):
self.neighbors = []
self.messages = {}
self.prob = None
self.alpha = 10.
self.beta = 5.
def add_neighbor(self, node):
"""
add neighboring node
Parameters
----------
node : Node
neighboring node
"""
self.neighbors.append(node)
def get_neighbors(self):
"""
get neighbor nodes
Returns
-------
neighbors : list
list containing neighbor nodes
"""
return self.neighbors
def init_messeges(self):
"""
initialize messages from neighbor nodes
"""
for neighbor in self.neighbors:
self.messages[neighbor] = np.ones(shape=(2,)) * 0.5
def marginalize(self):
"""
calculate probability
"""
prob = reduce(lambda x, y: x * y, self.messages.values())
self.prob = prob / prob.sum()
def send_message_to(self, node):
"""
calculate message to be sent to the node
Parameters
----------
node : Node
node to send computed message
Returns
-------
message : np.ndarray (2,)
message to be sent to the node
"""
message_from_neighbors = reduce(lambda x, y: x * y, self.messages.values()) / self.messages[node]
F = np.exp(self.beta * (2 * np.eye(2) - 1))
message = F.dot(message_from_neighbors)
node.messages[self] = message / message.sum()
def likelihood(self, value):
"""
calculate likelihood via observation, which is messege to this node
Parameters
----------
value : int
observed value -1 or 1
"""
assert (value == -1) or (value == 1), "{} is not 1 or -1".format(value)
message = np.exp(self.alpha * np.array([-value, value]))
self.messages[self] = message / message.sum()
class MarkovRandomField(object):
def __init__(self):
self.nodes = {}
def add_node(self, location):
"""
add a new node at the location
Parameters
----------
location : tuple
key to access the node
"""
self.nodes[location] = Node()
def get_node(self, location):
"""
get the node at the location
Parameters
----------
location : tuple
key to access the corresponding node
Returns
-------
node : Node
the node at the location
"""
return self.nodes[location]
def add_edge(self, key1, key2):
"""
add edge between nodes corresponding to key1 and key2
Parameters
----------
key1 : tuple
The key to access one of the nodes
key2 : tuple
The key to access the other node.
"""
self.nodes[key1].add_neighbor(self.nodes[key2])
self.nodes[key2].add_neighbor(self.nodes[key1])
def sum_product_algorithm(self, iter_max=10):
"""
Perform sum product algorithm
1. initialize messages
2. send messages from each node to neighboring nodes
3. calculate probabilities using the messages
Parameters
----------
iter_max : int
number of maximum iteration
"""
for node in self.nodes.values():
node.init_messeges()
for i in xrange(iter_max):
print i
for node in self.nodes.values():
for neighbor in node.get_neighbors():
node.send_message_to(neighbor)
for node in self.nodes.values():
node.marginalize()
def denoise(img, n_iter=20):
mrf = MarkovRandomField()
len_x, len_y = img.shape
X = range(len_x)
Y = range(len_y)
for location in itertools.product(X, Y):
mrf.add_node(location)
for x, y in itertools.product(X, Y):
for dx, dy in itertools.permutations(range(2), 2):
try:
mrf.add_edge((x, y), (x + dx, y + dy))
except Exception:
pass
for location in itertools.product(X, Y):
node = mrf.get_node(location)
node.likelihood(img[location])
mrf.sum_product_algorithm(n_iter)
denoised = np.zeros_like(img)
for location in itertools.product(X, Y):
node = mrf.get_node(location)
denoised[location] = 2 * np.argmax(node.prob) - 1
return denoised
def main():
img_original = 2 * (imread(ORIGINAL_IMAGE) / 255).astype(np.int) - 1
img_noise = 2 * (imread(NOISY_IMAGE) / 255).astype(np.int) - 1
img_denoised = denoise(img_noise, 10)
print "error rate before"
print np.sum((img_original != img_noise).astype(np.float)) / img_noise.size
print "error rate after"
print np.sum((img_denoised != img_original).astype(np.float)) / img_noise.size
imsave(DENOISED_IMAGE, (img_denoised + 1) / 2 * 255)
if __name__ == '__main__':
main()
résultat
Image restaurée en supprimant le bruit
Résultat de sortie du terminal
error rate before
0.109477124183
error rate after
0.0316993464052
Le taux d'erreur par rapport à l'image d'origine est réduit.
À la fin
Si le but est de supprimer le bruit, la méthode de coupe graphique semble être la plus précise.
Recommended Posts