How to use R and the igraph package to segment images.

How to use R and the igraph package to segment images.

This is not intended to be an efficient or accurate method. If you have any facility in this sort of thing your first thought might be “why would anyone do it like this?” Sorry, can’t help you! There might also be some discussion of whether walktrap communities are an appropriate way to estimate lambda connectedness. Who cares? It’s fun.

Step 1. Load your images.

library(rgdal)
filename <- file.choose()
img <- readGDAL(filename)

If this doesn’t work for you for some reason, I’m not interested in hearing about it. It worked for me. Here’s the proof:

Step 2. The goal is to cover the image with a regular hexagonal tiling. Each hexagon will be assigned a color equal to the mean of the colors it covers in the original full resolution image.

CELLSIZE= 19
heximage <- spsample(img,type=”hexagonal”,cellsize=CELLSIZE)
HH <- over(HexPoints2SpatialPolygons(heximage),img,fn=mean)
HH[is.na(HH)] <- 0
rgbcolors <- apply(HH, 1, function(c)
     rgb(red=c[[1]],
     green=c[[2]],
     blue=c[[3]],
     maxColorValue=255))
plot(HexPoints2SpatialPolygons(heximage),col=rgbcolors)

Step 3. I’m using the function delvor from the alphahull package to find the edges of the triangulation of the hexagonal tiling. That doesn’t quite make sense, but my meaning should be clear.

library(alphahull)
del.hexpoints <- delvor(jitter( heximage@coords,factor=0.00001))
plot(del.hexpoints,wlines=”del”,wpoints=F,add=T)
mesh.distances <- apply(del.hexpoints$mesh,1,
      function(aa) sqrt( (aa[5]-aa[3])^2 + (aa[6] — aa[4])^2))

Step 4. I’m not going to comment on the extremely important step 4. You can figure it out.

library(colorspace)
library(igraph)
source("de2000.R")
rows <- which(round(mesh.distances) < 2*CELLSIZE — 1)
hexpoints.graph <- graph.edgelist(
     del.hexpoints$mesh[rows,1:2],
     directed=F)
labcolors <- LAB(apply(HH,2,function(c) c/255))
V(hexpoints.graph)$color <- rownames(coords(labcolors))
E(hexpoints.graph)$weight <- apply(
 get.edges(hexpoints.graph,1:ecount(hexpoints.graph)),
 1,
 function(edg) {
 c1 = labcolors[edg[1]]
 c2 = labcolors[edg[2]]
 return(deltaE2000(c1,c2,kl=1,kc=2,kh=1))
 })
E(hexpoints.graph)$weight <- 1/E(hexpoints.graph)$weight

Well perhaps some of that was not completely transparent, but my posts are seldom read anyway. If and when something happens and someone complains perhaps I can say more.

Step 5. Use walktrap communities to partition the hexagons, and find the polygons surrounding each partition.

library(rgeos)
WT <- walktrap.community(hexpoints.graph,steps=40)
memWT <- membership(WT)
classes <- unique(memWT)
polylist <- sapply(classes, function(cl) {
 newpoly <- gUnaryUnion(HexPoints2SpatialPolygons(heximage[which(memWT==cl)]))
 newpoly <- gSimplify(newpoly,tol=10,topologyPreserve=T)
 return(as(newpoly,"gpc.poly"))
})
p <- polylist[[1]]
for (x in polylist[-1]) { p <- append.poly(p,x)}

Step 6. To illustrate the process, I’ll plot the original image, the graph of the image with the edges scaled to show the difference between adjacent colors, and finally the found segments of the image drawn over the haxagonal tiling.

par(mfrow=c(1,3),bg=”black”)
image(img,red=”band1",green=”band2",blue=”band3")
plot(hexpoints.graph,
 layout=jitter(heximage@coords,factor=1),
 vertex.size=2,
 vertex.label=NA,
 vertex.frame.color=rgbcolors,
 vertex.color=rgbcolors,
 edge.color=rgb(1,1,1),
 edge.width = .01*E(hexpoints.graph)$weight)
plot(HexPoints2SpatialPolygons(heximage),col=rgbcolors)
plot(p,add=T,poly.args = list(lwd = 3))

Step 7. de2000.R The CIEDE2000 color difference function comes from Gaurav Sharma. This R code is a direct translation from his matlab function. I’ve left some of his comments in the code, but removed others in the copy/paste process. Sloppy of me, I know.

deltaE2000 <- function(Labstd,Labsample,kl=1,kc=1,kh=1) {
lstd = coords(Labstd)[,1]
 astd = coords(Labstd)[,2]
 bstd = coords(Labstd)[,3]
 Cabstd = sqrt(astd^2+bstd^2)
 
 lsample = coords(Labsample)[,1]
 asample = coords(Labsample)[,2]
 bsample = coords(Labsample)[,3]
 Cabsample = sqrt(asample^2 + bsample^2)

Cabarithmean = (Cabstd + Cabsample)/2
 G = 0.5* ( 1 — sqrt( (Cabarithmean^7)/(Cabarithmean^7 + 25^7)))
 apstd = (1+G)*astd
 apsample = (1+G)*asample
 apstd = (1+G)*astd
 apsample = (1+G)*asample
 Cpsample = sqrt(apsample^2+bsample^2)
 Cpstd = sqrt(apstd^2+bstd^2)
 # Compute product of chromas and locations at which it is zero for use later
 Cpprod = (Cpsample*Cpstd);
 zcidx = which(Cpprod == 0);

# Ensure hue is between 0 and 2pi
 # NOTE: MATLAB already defines atan2(0,0) as zero but explicitly set it
 # just in case future definitions change
 hpstd = atan2(bstd,apstd);
 hpstd = hpstd+2*pi*(hpstd < 0)
 hpstd[which( (abs(apstd)+abs(bstd))== 0) ] = 0
 hpsample = atan2(bsample,apsample)
 hpsample = hpsample+2*pi*(hpsample < 0)
 hpsample[which( (abs(apsample)+abs(bsample))==0) ] = 0
 dL = lsample-lstd
 dC = Cpsample-Cpstd
 dhp = hpsample-hpstd
 dhp = dhp — 2*pi* (dhp > pi )
 dhp = dhp + 2*pi* (dhp < (-pi) )
 dhp[zcidx ] = 0
 dH = 2*sqrt(Cpprod)*sin(dhp/2)
 Lp = (lsample+lstd)/2
 Cp = (Cpstd+Cpsample)/2
 hp = (hpstd+hpsample)/2
 hp = hp — ( abs(hpstd-hpsample) > pi ) *pi
 hp = hp+ (hp < 0) *2*pi
 hp[zcidx] = hpsample[zcidx]+hpstd[zcidx]
 Lpm502 = (Lp-50)^2
 Sl = 1 + 0.015*Lpm502/sqrt(20+Lpm502)
 Sc = 1+0.045*Cp
 T = 1–0.17*cos(hp — pi/6 ) + 0.24*cos(2*hp) + 0.32*cos(3*hp+pi/30) -0.20*cos(4*hp-63*pi/180)
 Sh = 1 + 0.015*Cp*T
 delthetarad = (30*pi/180)*exp(- ( (180/pi*hp-275)/25)^2)
 Rc = 2*sqrt((Cp^7)/(Cp^7 + 25^7))
 RT = — sin(2*delthetarad)*Rc
 klSl = kl*Sl
 kcSc = kc*Sc
 khSh = kh*Sh
 de00 = sqrt( (dL/klSl)^2 + (dC/kcSc)^2 + (dH/khSh)^2 + RT*(dC/kcSc)*(dH/khSh) )
 return(as.numeric(de00))
}