Category Archives: Image Processing

Frankenstein

Remember me, remember me, but ah! forget my fate (Dido’s Lament, Henry Purcell)

A Voronoi diagram divides a plane based on a set of original points. Each polygon, or Voronoi cell, contains an original point and all that are closer to that point than any other.

This is a nice example of a Voronoi tesselation. You can find good explanations of Voronoi diagrams and Delaunay triangulations here (in English) or here (in Spanish).

A grayscale image is simply a matrix where darkness of pixel located in coordinates (i, j) is represented by the value of its corresponding element of the matrix: a grayscale image is a dataset. This is a Voronoi diagraman of Frankenstein:

To do it I followed the next steps:

  1. Read this image
  2. Convert it to gray scale
  3. Turn it into a pure black and white image
  4. Obtain a random sample of black pixels (previous image corresponds to a sample of 6.000 points)
  5. Computes the Voronoi tesselation

Steps 1 to 3 were done with imager, a very appealing package to proccess and analice images. Step 5 was done with deldir, also a convenient package which computes Delaunay triangulation and the Dirichlet or Voronoi tessellations.

The next grid shows tesselations for sample size from 500 to 12.000 points and step equal to 500:

I gathered all previous images in this gif created with magick, another amazing package of R I discovered recently:

This is the code:

library(imager)
library(dplyr)
library(deldir)
library(ggplot2)
library(scales)

# Download the image
file="http://ereaderbackgrounds.com/movies/bw/Frankenstein.jpg"
download.file(file, destfile = "frankenstein.jpg", mode = 'wb')

# Read and convert to grayscale
load.image("frankenstein.jpg") %>% grayscale() -> x

# This is just to define frame limits
x %>% 
  as.data.frame() %>% 
  group_by() %>% 
  summarize(xmin=min(x), xmax=max(x), ymin=min(y), ymax=max(y)) %>% 
  as.vector()->rw

# Filter image to convert it to bw
x %>%
  threshold("45%") %>% 
  as.cimg() %>% 
  as.data.frame() -> df

# Function to compute and plot Voronoi tesselation depending on sample size
doPlot = function(n)
{
  #Voronoi tesselation
  df %>% 
  sample_n(n, weight=(1-value)) %>% 
  select(x,y) %>% 
  deldir(rw=rw, sort=TRUE) %>% 
  .$dirsgs -> data

  # This is just to add some alpha to lines depending on its longitude
  data %>% 
    mutate(long=sqrt((x1-x2)^2+(y1-y2)^2),
         alpha=findInterval(long, quantile(long, probs = seq(0, 1, length.out = 20)))/21)-> data

  # A little bit of ggplot to plot results
  data %>% 
    ggplot(aes(alpha=(1-alpha))) +
    geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), color="black", lwd=1) +
    scale_x_continuous(expand=c(0,0))+
    scale_y_continuous(expand=c(0,0), trans=reverse_trans())+
    theme(legend.position  = "none",
            panel.background = element_rect(fill="white"),
            axis.ticks       = element_blank(),
            panel.grid       = element_blank(),
            axis.title       = element_blank(),
            axis.text        = element_blank())->plot

return(plot)
}

# I call the previous function and store resulting plot in jpeg format
i=5000
name=paste0("frankie",i,".jpeg")
jpeg(name, width = 600, height = 800, units = "px", quality = 100)
doPlot(i)
dev.off()

# Once all images are stored I can create gif
library(magick)
frames=c()
images=list.files(pattern="jpeg")

for (i in length(images):1)
{
  x=image_read(images[i])
  x=image_scale(x, "300")
  c(x, frames) -> frames
}
animation=image_animate(frames, fps = 2)
image_write(animation, "Frankenstein.gif")

Playing With Julia (Set)

Viento, me pongo en movimiento y hago crecer las olas del mar que tienes dentro (Tercer Movimiento: Lo de Dentro, Extremoduro)

I really enjoy drawing complex numbers: it is a huge source of entertainment for me. In this experiment I play with the Julia Set, another beautiful fractal like this one. This is what I have done:

  • Choosing the function f(z)=exp(z3)-0.621
  • Generating a grid of complex numbers with both real and imaginary parts in [-2, 2]
  • Iterating f(z) over the grid a number of times so zn+1 = f(zn)
  • Drawing the resulting grid as I did here
  • Gathering all plots into a GIF with ImageMagick as I did in my previous post: each frame corresponds to a different number of iterations

This is the result:

julia

I love how easy is doing difficult things in R. You can play with the code changing f(z) as well as color palettes. Be ready to get surprised:

library(ggplot2)
library(dplyr)
library(RColorBrewer)
setwd("YOUR WORKING DIRECTORY HERE")
dir.create("output")
setwd("output")
f = function(z,c) exp(z^3)+c
# Grid of complex
z0 <- outer(seq(-2, 2, length.out = 1200),1i*seq(-2, 2, length.out = 1200),'+') %>% c()
opt <-  theme(legend.position="none",
              panel.background = element_rect(fill="white"),
              plot.margin=grid::unit(c(1,1,0,0), "mm"),
              panel.grid=element_blank(),
              axis.ticks=element_blank(),
              axis.title=element_blank(),
              axis.text=element_blank())
for (i in 1:35)
{
  z=z0
  # i iterations of f(z)
  for (k in 1:i) z <- f(z, c=-0.621) df=data.frame(x=Re(z0), y=Im(z0), z=as.vector(exp(-Mod(z)))) %>% na.omit() 
  p=ggplot(df, aes(x=x, y=y, color=z)) + 
    geom_tile() + 
    scale_x_continuous(expand=c(0,0))+
    scale_y_continuous(expand=c(0,0))+
    scale_colour_gradientn(colours=brewer.pal(8, "Paired")) + opt
  ggsave(plot=p, file=paste0("plot", stringr::str_pad(i, 4, pad = "0"),".png"), width = 1.2, height = 1.2)
}
# Place the exact path where ImageMagick is installed
system('"C:\\Program Files\\ImageMagick-6.9.3-Q16\\convert.exe" -delay 20 -loop 0 *.png julia.gif')
# cleaning up
file.remove(list.files(pattern=".png"))

Zooming

You don’t have to be beautiful to turn me on (Kiss, Prince)

I discovered recently how easy is to create GIFs with R using ImageMagick and I feel like a kid with a new toy. To begin this new era of my life as R programmer I have done this:

zooming
First of all, read this article: it explains very well how to start doing GIFs from scratch. The one I have done is inspired in this previous post where I take a set of complex numbers to transform and color it using HSV technique. In this case I use this next transformation: f(z)= -Im(z)+(Re(z)+0.5*Im(z))*1i

Modifying the range of Real and Imaginary parts of complex numbers I obtain the zooming  effect. The code is very simple. Play with it changing the transformation or the animation options. Send me your creations, I would love to see them:

library(dplyr)
library(ggplot2)
dir.create("output")
setwd("output")
id=1 # label tO name plots
for (i in seq(from=320, to=20, length.out = 38)){
z=outer(seq(from = -i, to = i, length.out = 300),1i*seq(from = -i, to = i, length.out = 500),'+') %>% c()
z0=z
for (k in 1:100) z <- -Im(z)+(Re(z)+0.5*Im(z))*1i
df=data.frame(x=Re(z0),
y=Im(z0),
h=(Arg(z)<0)*1+Arg(z)/(2*pi), s=(1+sin(2*pi*log(1+Mod(z))))/2, v=(1+cos(2*pi*log(1+Mod(z))))/2) %>% mutate(col=hsv(h,s,v))
ggplot(df, aes(x, y)) +
geom_tile(fill=df$col)+
scale_x_continuous(expand=c(0,0))+
scale_y_continuous(expand=c(0,0))+
labs(x=NULL, y=NULL)+
theme(legend.position="none",
panel.background = element_rect(fill="white"),
plot.margin=grid::unit(c(1,1,0,0), "mm"),
panel.grid=element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text=element_blank())
ggsave(file=paste0("plot",stringr::str_pad(id, 4, pad = "0"),".png"), width = 1, height = 1)
id=id+1
}
system('"C:\\Program Files\\ImageMagick-6.9.3-Q16\\convert.exe" -delay 10 -loop 0 -duplicate 1,-2-1 *.png zooming.gif')
# cleaning up
file.remove(list.files(pattern=".png"))

butteRfly

Float like a butterfly, sting like a bee (Muhammad Ali)

The Butterfly Curve was discovered by Temple H. Fay when he was in Southern University, Mississippi, and rapidly gained the attention of students and mathematicians because of its beautiful simmetry. Small dots of this plot are generated according to parametric equations of the Butterfly Curve. Big dots are randomdly distributed over the canvas:

butterfly4

This is the code to create butterflies:

library(ggplot2)
npoints=500
npointsb=1200
opt=theme(legend.position="none",
panel.background = element_blank(),
panel.grid = element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text =element_blank())
t=seq(0,10*pi,length=npointsb)
butterfly=data.frame(x=sin(t)*(exp(1)^cos(t)-2*cos(4*t)-(sin(t/12))^5), y=cos(t)*(exp(1)^cos(t)-2*cos(4*t)-(sin(t/12))^5), s=runif(npointsb, min=.1, max=10), f=factor(sample(1:10,npointsb,TRUE)), a=runif(npointsb,min=.1, max=.4))
points=data.frame(x=runif(npoints,-4,4), y=runif(npoints,-3,5), s=runif(npoints,min=30, max=50), f=factor(sample(1:10,npoints,TRUE)), a=runif(npoints,min=.05, max=.15))
data=rbind(butterfly, points)
ggplot(data, aes(x, y, colour=f))+geom_point(alpha=data$a,size=data$s)+opt

The Zebra Of Riemann

Mathematics is the art of giving the same name to different things (Henri Poincare)

Many surveys among experts point that demonstration of the Riemann Hypothesis is the most important pending mathematical issue in this world. This hypothesis is related to Riemann zeta function, which is supossed to be zero only for those complex whose real part is equal to 1/2 (this is the conjecture of Riemann itself). Confirming the conjecture would imply deep consequences in prime numbers teory and also in our knowledge of their properties. Next plot represents the argument of zeta function over complex with real part between -75 and 5 (x axis) and imaginary part between -40 and 40 (y axis). An explanation of this kind of graph can be found here. Does not it remind you of something?
zebra
Solving the hypothesis of Riemann is one of the seven Clay Mathematics Institute Millennium Prize Problems. Maybe zebras keep the secret to solve it in their skins.

This is the code to plot the graph:

require(pracma)
z=outer(seq(-75, 5, by =.1),1i*seq(-40, 40, by =.1),'+')
z=apply(z, c(1,2), function(x) Arg(zeta(x)))
opt=theme(legend.position="none",
panel.background = element_blank(),
panel.grid = element_blank(),
axis.ticks=element_blank(),
axis.title=element_blank(),
axis.text =element_blank())
z=data.frame(expand.grid(x=seq(ncol(z)), y=seq(nrow(z))), z=as.vector(z))
require(ggplot2)
ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradientn(colours=c("black","white", "gray80", "gray40")) + opt

floweR

It is the time you have wasted for your rose that makes your rose so important (Antoine de Saint-Exupéry, The Little Prince)

Yesterday I found a package called circular and I could not suppress to do this:

rose

Make your own floweRs:

require(circular)
for (i in 1:8)
{
  rose.diag(circular(runif(5000, 0, 2*pi)), 
            bins = (25-round(abs(jitter(0, amount=2*i)))), 
            axes=FALSE,
            border=rgb(255,50+204*((i-1)/8),50+204*((i-1)/8), alpha=255, max=255),
            ticks = FALSE,
            col=rgb(250,0+204*((i-1)/8),0+204*((i-1)/8), alpha=255, max=255), 
            control.circle=circle.control(lty=0), 
            shrink=0.22+(i-1)*(2-0.22)/8)
  par(new=TRUE)  
}
for (i in 1:150)
  {
  q = runif(1)*pi*2
  r = sqrt(runif(1))
  x = (0.18*r)* cos(q)
  y = (0.18*r)* sin(q)
  points(x, y, 
         col = rgb(255 ,sample(80:120, 1), 0, alpha= 95, max=255),
         bg = rgb(255 ,sample(100:200, 1), 0, alpha= 95, max=255), pch=21, cex=2)
  }

Blurry Fractals

Beauty is the first test; there is no permanent place in the world for ugly mathematics (G. H. Hardy)

Newton basin fractals are the result of iterating Newton’s method to find roots of a polynomial over the complex plane. It maybe sound a bit complicated but is actually quite simple to understand. Those who would like to read some more about Newton basin fractals can visit this page.

This fractals are very easy to generate in R and produce very nice images. Making a small number of iterations, resulting images seems to be blurred when are represented with tile geometry in ggplot. Combined with palettes provided by RColorBrewer give rise to very interesting images. Here you have some examples:

Result for f(z)=z3-1 and palette equal to Set3:Blurry1-Set3Result for f(z)=z4+z-1 and palette equal to Paired:Blurry2-PairedResult for f(z)=z5+z3+z-1 and palette equal to Dark2:Blurry3-Dark2Here you have the code. If you generate nice pictures I will be very grateful if you send them to me:

library(ggplot2)
library(numDeriv)
library(RColorBrewer)
library(gridExtra)
## Polynom: choose only one or try yourself
f  <- function (z) {z^3-1}        #Blurry 1
#f  <- function (z) {z^4+z-1}     #Blurry 2
#f  <- function (z) {z^5+z^3+z-1} #Blurry 3
z <- outer(seq(-2, 2, by = 0.01),1i*seq(-2, 2, by = 0.01),'+')
for (k in 1:5) z <- z-f(z)/matrix(grad(f, z), nrow=nrow(z))
## Supressing texts, titles, ticks, background and legend.
opt <- theme(legend.position="none",
             panel.background = element_blank(),
             axis.ticks=element_blank(), 
             axis.title=element_blank(), 
             axis.text =element_blank())
z <- data.frame(expand.grid(x=seq(ncol(z)), y=seq(nrow(z))), z=as.vector(exp(-Mod(f(z)))))
# Create plots. Choose a palette with display.brewer.all()
p1 <- ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradientn(colours=brewer.pal(8, "Paired")) + opt
p2 <- ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradientn(colours=brewer.pal(7, "Paired")) + opt
p3 <- ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradientn(colours=brewer.pal(6, "Paired")) + opt
p4 <- ggplot(z, aes(x=x, y=y, color=z)) + geom_tile() + scale_colour_gradientn(colours=brewer.pal(5, "Paired")) + opt
# Arrange four plots in a 2x2 grid
grid.arrange(p1, p2, p3, p4, ncol=2)

Warholing Grace With Clara

Do not believe anything: what artists really do is to hang around all day (Paco de Lucia)

Andy Warhol was mathematician. At least, he knew how clustering algorithms work. I am pretty sure of this after doing this experiment.  First of all, let me introduce you to the breathtaking Grace Kelly:

936full-grace-kelly

In my previous post I worked also with images showing how simple is to operate with them since they are represented by matrices. This is another example of this. Third dimension of an image matrix is an 3D array representing color of pixels in (r, g, b) format. Applying a cluster algorithm over this information generates groups of pixels with similar color. I used cluster package and because of the high size of picture I decided to use clara algorithm which is extremely fast. Apart of its high speed, another advantage of clara is that clusters are represented by real elements of the population, called medoids, instead of being by average individuals as k-means do. It fits very well with my purposes because once clusters are calculated I only have to change each pixel by its medoid and plot it. Setting clara to divide pixels into 2 groups, generates a 2 colored image. Setting it to 3 groups, generates a 3 colored one and so on. Following, you can find results from 2 to 7 groups:

Warholing2 Warholing3 Warholing4
Warholing5 Warholing6 Warholing7

Working with samples can be a handicap, maybe less important than the speed it produces. Sometimes images generated by n groups seems to be worse fitted than the one generated by n-1 groups. You can see it in this video, where results from 1 to 60 groups are presented sequentially. It only takes 42 seconds.

Here you have the code. Feel free to warholing:

library("biOps")
library("abind")
library("reshape")
library("reshape2")
library("cluster")
library("sp")
#######################################################################################################
#Initialization
#######################################################################################################
x     <- readJpeg("936full-grace-kelly.jpg")
plot(x)
#######################################################################################################
#Data
#######################################################################################################
data <- merge(merge(melt(x[,,1]), melt(x[,,2]), by=c("X1", "X2")), melt(x[,,3]), by=c("X1", "X2"))
colnames(data) <- c("X1", "X2", "r", "g", "b")
#######################################################################################################
#Clustering
#######################################################################################################
colors <- 5
clarax <- clara(data[,3:5], colors)
datacl   <- data.frame(data, clarax$cluster)
claradf2 <- as.data.frame(clarax$medoids)
claradf2$id <- as.numeric(rownames(claradf2))
colnames(claradf2) <- c("r", "g", "b", "clarax.cluster")
claradf <- merge(datacl, claradf2, by=c("clarax.cluster"))
colnames(claradf) <-c("clarax.cluster", "X1", "X2", "r.x", "g.x", "b.x", "r.y", "g.y", "b.y")
datac <- claradf[do.call("order", claradf[c("X1", "X2")]), ]
x1<-acast(datac[,c(2,3,7)], X1~X2, value.var="r.y") 
x2<-acast(datac[,c(2,3,8)], X1~X2, value.var="g.y") 
x3<-acast(datac[,c(2,3,9)], X1~X2, value.var="b.y") 
warhol <- do.call(abind, c(list(x1,x2,x3), along = 3))
plot(imagedata(warhol))
writeJpeg(paste("Warholing", as.character(colors), ".jpg", sep=""), imagedata(warhol))

Face To Face With Marilyn Monroe

Symmetry is what we see at a glance (Blaise Pascal)

Ladies and gentlement, the beautiful Marilyn Monroe:

marilyn-monroe3

There are several image processing packages in R. In this experiment I used biOps, which turns images into 3D matrices. The third dimension is a 3-array corresponding to (r, g, b) color of pixel defined by two other dimensions. Since images are defined by matrices, you can do simple operations to produce interesting images from the original one. For example, this is what happens swaping last half of columns by first one and preserving its order:

IMG01-Interchange

It is also very simple to generate two artificial symmetrical faces matching each half with itself as a mirror. By the way, I prefer the first one: maybe two sexy moles are better than just one.

IMG03-Symmetric2 IMG02-Symmetric1

Let’s introduce a bit of randomness. This is what happens when you take an uniform sample of rows and columns (first image) or an uniform sample of pixels of the image (second one):

IMG04-Uniform1 IMG05-Uniform2

And this is what happens when you divide image into blocks and mix them randomly:

IMG07-Mosaic

There is a funny and useful function called jitter which adds a small amount of noise to a numeric vector. What happens when you jitter every pixel of the image? As you can see, It becomes very vintage:

IMG06-Jitter

What if you transpose matrix? What if you change every color by another one? What if you change only a small range of them? What if you sum two images? What if you translate rgb colors into a grey scale? What if …? I answered some of these questions already and results are nice as well. After all, Marilyn can be represented as a simple matrix. Or maybe not.

Make your own experiments:

library("biOps")
library("abind")
#############################################################
#Read Original Image
#############################################################
x     <- readJpeg("marilyn-monroe3.jpg")
plot(x)
#############################################################
#1. Swap
#############################################################
plot(imagedata(abind(x[,(ncol(x)/2):ncol(x),], x[,1:(ncol(x)/2),] , along=2)))
dev.copy(jpeg,filename="IMG01-Swap.jpg");
dev.off ();
#############################################################
#2. Artificial Symmetrical faces
#############################################################
plot(imagedata(abind(x[,1:(ncol(x)/2),], x[,(ncol(x)/2):1,] , along=2)))
dev.copy(jpeg,filename="IMG02-Symmetric1.jpg");
dev.off ();
plot(imagedata(abind(x[,ncol(x):(ncol(x)/2),], x[,(ncol(x)/2):ncol(x),] , along=2)))
dev.copy(jpeg,filename="IMG03-Symmetric2.jpg");
dev.off ();
#############################################################
#3. Uniform sampling over axis points
#############################################################
x2   <- aperm(array(255, dim = c(3, ncol(x), nrow(x))))
rows <- sample(1:nrow(x), round(nrow(x)*0.80), replace = FALSE)
cols <- sample(1:ncol(x), round(ncol(x)*0.80), replace = FALSE)
for (i in 1:length(rows)) 
{
  for (j in 1: length(cols)) 
    {
    x2[rows[i], cols[j],1]<-x[rows[i], cols[j],1] 
    x2[rows[i], cols[j],2]<-x[rows[i], cols[j],2] 
    x2[rows[i], cols[j],3]<-x[rows[i], cols[j],3] 
  }
}
plot(imagedata(x2))
dev.copy(jpeg,filename="IMG04-Uniform1.jpg");
dev.off ();
#############################################################
#4. Uniform sampling over pixels
#############################################################
m2 <- matrix(rbinom(nrow(x)*ncol(x),1,0.5),nrow(x),ncol(x))
x4<- do.call(abind, c(list(x[,,1]*m2+(m2==0)*255,x[,,2]*m2+(m2==0)*255,x[,,3]*m2+(m2==0)*255), along = 3))
plot(imagedata(x4))
dev.copy(jpeg,filename="IMG05-Uniform2.jpg");
dev.off ();
#############################################################
#6. Jittering
#############################################################
x1<-mapply(as.matrix(x[,,1]), FUN=function(x) 
  {z<-round(x+jitter(0, amount=50)) 
  if(z<0|z>255) x else z})
x1 <- matrix(x1, nrow = nrow(x),ncol = ncol(x))
x2<-mapply(as.matrix(x[,,2]), FUN=function(x) 
{z<-round(x+jitter(0, amount=50)) 
 if(z<0|z>255) x else z})
x2 <- matrix(x2, nrow = nrow(x),ncol = ncol(x))
x3<-mapply(as.matrix(x[,,3]), FUN=function(x) 
{z<-round(x+jitter(0, amount=50)) 
 if(z<0|z>255) x else z})
x3 <- matrix(x3, nrow = nrow(x),ncol = ncol(x))
x4<- do.call(abind, c(list(x1,x2,x3), along = 3))
plot(imagedata(x4))
dev.copy(jpeg,filename="IMG06-Jitter.jpg");
dev.off ();
#############################################################
#7. Mosaic
#############################################################
sptr <- 6 #Row splits
rnkr <- sample(1:sptr, size = sptr, replace = FALSE)
wthr <- floor(nrow(x)/sptr) #Splits width (row)
rnkr <- as.vector(sapply(rnkr, function (x) rep(x,wthr)))
rnkr <- rnkr*10E6+seq(1, length(rnkr), by=1)
rnkr <- rank(rnkr)
sptc <- round(ncol(x)/wthr)
rnkc <- sample(1:sptc, size = sptc, replace = FALSE)
wthc <- floor(ncol(x)/sptc) #Splits width (row)
rnkc <- as.vector(sapply(rnkc, function (x) rep(x,wthc)))
rnkc <- rnkc*10E6+seq(1, length(rnkc), by=1)
rnkc <- rank(rnkc)
x2<-x[1:length(rnkr),1:length(rnkc),]
x2<-x[rank(rnkr),rank(rnkc),]
plot(imagedata(x2))
dev.copy(jpeg,filename="IMG07-Mosaic.jpg");
dev.off ();