a CRF model for denoising

R의 package CRF를 활용하여,

Bishop(2006, pp.387-389)에 소개된 denoising MRF를 구현해 보았다.

다소 큰 이미지는 adjacency matrix를 만들기에 메모리가 부족하기에,

작은 이미지 화일로 실험을 해 보았다.

node potential 결정하는 게 어려워서 그냥 0으로 했다.

edge potential를 포함해서 parameter를 적절히 정하는 법을 생각해봐야 할 듯 하다.

근데 train은 왜 안 되지?

결과. 왼쪽 위 : 원본,  왼쪽 아래 : 10% noise, 오른쪽 위 : icm으로 보정, 오른쪽 아래 : lbp로 보정

source : crf_denoising.R, image1 : alpha.jpg, image2 : R.jpg

crf_denoising.R

pr.noise=     # noising 과정에서 원본이 그대로 보존되는 확률

crf$par=       # parameter 설정, 원본의 potential, 원본의 pixel 사이의 edge potential, 원본과 noised 그림 간 edge potential

<source>

library(jpeg)

library(RgoogleMaps)

image.ideal <- readJPEG('alpha.jpg') # 0.jpg 파일에서 그림을 읽어서 image.ideal에 저장한다.

#image.ideal <- readJPEG('R.jpg')

image.ideal <- RGB2GRAY(image.ideal) # image.ideal를 흑백 사진으로 변환한다.

image.ideal[image.ideal>0.5]=1

image.ideal[image.ideal<=0.5]=0

image.ideal=1-image.ideal

str(image.ideal)

par(mfcol=c(2,2))

pr.noise=0.9

noise=runif(length(image.ideal))>pr.noise

image.noised=image.ideal+noise

image.noised[image.noised==2]=0

image(t(image.ideal), col=c(0,1), ylim=c(1,0))

image(t(image.noised), col=c(0,1), ylim=c(1,0))

library(CRF)

w=ncol(image.noised); h=nrow(image.noised)

n.nodes=2*w*h; n.states=2; n.nf=1; n.ef=1; n.edges=((w-1)*h+(h-1)*w+h*w);

mat.adj = matrix(0, nrow=n.nodes, ncol=n.nodes)

for (x in 2:(w-1)) 

  for (y in 2:(h-1)) {

    mat.adj[(x-1)*h+y,(x+1-1)*h+y]=1

    mat.adj[(x-1)*h+y,(x-1)*h+y+1]=1

    mat.adj[(x-1)*h+y,(x-1-1)*h+y]=1

    mat.adj[(x-1)*h+y,(x-1)*h+y-1]=1

  }  

for (x in 1:(w-1)) {

  y=1; mat.adj[(x-1)*h+y, (x-1+1)*h+y]=1

  y=h; mat.adj[(x-1)*h+y, (x-1+1)*h+y]=1

  }  

for (y in 1:(h-1)) {

  x=w; mat.adj[(x-1)*h+y,(x-1)*h+y+1]=1  

  x=1; mat.adj[(x-1)*h+y,(x-1)*h+y+1]=1  

}  

for (x in 1:w) 

  for (y in 1:h) {

    mat.adj[(x-1)*h+y,w*h+(x-1)*h+y]=1

  }

# testing

#library(igraph)

#g<-graph.adjacency(mat.adj, mode="undirected")

#plot(g)

#= making edges b2 i,j

#mat.adj[i,j]=1

#= making edges b2 x,y(ideal image)

#mat.adj[(x1-1)*h+y1, (x2-1)*h+y2]=1

mat.adj.t=t(mat.adj)

mat.adj=mat.adj | mat.adj.t

mat.adj=mat.adj*1

crf <- make.crf(mat.adj, n.states=n.states)

crf <- make.features(crf, n.nf, n.ef)

# node potential, edge potential among pixels of the ideal image 1,1

# edge potential between the pixels of the ideal image and the actual image 2

crf <- make.par(crf,3)

crf$node.par[,1,1]<- 1

#crf$node.par[,2,1]<- 

#1 ~ (w-1)*h+(h-1)*w,  edges between the pixels of the ideal image

#(w-1)*h+(h-1)*w+1 ~ (w-1)*h+(h-1)*w+h*w,  edges between the ideal and the actual image

edge.par.temp1 <-matrix(c(2,0,0,2), nrow=n.states, ncol=n.states)

edge.par.temp2 <-matrix(c(3,0,0,3), nrow=n.states, ncol=n.states)

edge.par1=array(edge.par.temp1, c(nrow(edge.par.temp1), ncol(edge.par.temp1),1))

edge.par2=array(edge.par.temp2, c(nrow(edge.par.temp2), ncol(edge.par.temp2),1))

for (i in 1:((w-1)*h+(h-1)*w)) crf$edge.par[[i]]=edge.par1

for (i in ((w-1)*h+(h-1)*w+1):((w-1)*h+(h-1)*w+h*w)) 

  crf$edge.par[[i]]=edge.par2

instances=matrix(c(image.ideal, image.noised)+1, nrow=1)

for (i in 1:10) {

  noise=runif(length(image.ideal))>pr.noise

  image.noised=image.ideal+noise

  image.noised[image.noised==2]=0

  instance = matrix(c(image.ideal, image.noised)+1, nrow=1)

  instances = rbind(instances,instance)

}

n.edges=crf$n.edges

node.fea <- lapply(1:11, function(i) matrix(1, crf$n.nf, crf$n.nodes))

edge.fea <- lapply(1:11, function(i) matrix(1, crf$n.ef, crf$n.edges))

#instances=matrix(c(image.ideal, image.noised)+1, nrow=11)

#train.crf(crf, instances, node.fea, edge.fea, trace=0)

node.pot = sum(image.noised==0)/sum(image.noised==1)

node.par = log(node.pot) 

node.par = 0

edge.par = log(pr.noise/(1-pr.noise))

edge2.par = ifelse(pr.noise>=0.95, 5, 

                   ifelse(pr.noise>=0.9, 4, 

                      ifelse(pr.noise>=0.8, 3, 2)))

node.par; edge.par

crf$par=c(node.par,3,edge.par); n1=node.fea[[1]]; e1=edge.fea[[1]]

crf.update(crf, node.fea=n1, edge.fea=e1)

v.clamped <- c(rep(0,length(image.ideal)),image.noised+1)

crf.clamped <- clamp.crf(crf, v.clamped)

image.ideal.recovered <- decode.icm(crf.clamped)

image.ideal.recovered <- matrix(image.ideal.recovered , ncol=w, nrow=h)

image(t(image.ideal.recovered), col=c(0,1), ylim=c(1,0))

image.ideal.recovered <- decode.lbp(crf.clamped, max.iter=100000)

image.ideal.recovered <- matrix(image.ideal.recovered , ncol=w, nrow=h)

image(t(image.ideal.recovered), col=c(0,1), ylim=c(1,0))