PRML 8.3.3に記載の通り、ノイズを含む観測画像をy、未知の元の画像をxとするマルコフ確率場において、反復条件付きモード(ICM)によって、yが与えられたときp(x|y)を極大化するxを求めます。図8.30と同様に、元の画像、ノイズを含む画像、ICMによって復元された画像とピクセルの復元率を表示します。

library(biOps)
frame()
set.seed(0)
par(mfcol=c(2, 3))
par(mar=c(0, 0, 1, 0))
data(logo)
source <- imagedata(logo[,,3])
NI <- nrow(source)
NJ <- ncol(source)
s <- ifelse(source < 128, -1, 1)
plot(imagedata((s + 1) / 2 * 255))
title("source", cex.main=1.0)

# add noise
y <- s * ((runif(length(s)) > 0.1) * 2 - 1)
plot(imagedata((y + 1) / 2 * 255))
title("y", cex.main=1.0)

doplot <- function(beta, eta, h) {
	x <- y
	repeat {
		energyDiff <- 0
		for (i in (1:NI)) {
			for (j in (1:NJ)) {
				energy <- function(xx) {
					result <- h * xx
					if (i > 1) {
						result <- result - beta * xx * x[i - 1,j]
					}
					if (i < NI) {
						result <- result - beta * xx * x[i + 1,j]
					}
					if (j > 1) {
						result <- result - beta * xx * x[i,j - 1]
					}
					if (j < NJ) {
						result <- result - beta * xx * x[i,j + 1]
					}
					result <- result - eta * xx * y[i,j]
					result
				}
				values <- c(x[i,j], x[i,j] * -1)
				energies <- energy(values)
				minEnergyIndex <- which(energies == min(energies))
				x[i,j] <- values[minEnergyIndex]
				energyDiff <- energyDiff + energies[minEnergyIndex] - energies[1]
			}
		}
		cat("difference of total energy");print(energyDiff)
		if (energyDiff == 0) {
			break
		}
	}
	restored = round(sum(s == x) / length(x) * 100, 1)
	plot(imagedata((x + 1) / 2 * 255))
	title(paste0("x (beta=", beta, " eta=", eta, " h=", h, ") ", restored, "%"), cex.main=1.0)
}

doplot(1.0, 2.1, 0)
doplot(1.0, 2.1, 1.5)
doplot(1.0, 2.1, -1.5)
doplot(5.0, 2.1, -1.5)