R중 마르코프 확률장의 교체 조건 모드

PRML8.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)