### Weaves two patent images and a latent image.
### The latent image is revealed upon overlapping
### the patent images.


library(png)

floydSteinberg <- function(input, greys){
  output <- matrix(NA, nrow=nrow(input), ncol=ncol(input))
  input  <- rbind(input, 0)
  input  <- cbind(input, 0)
  for (x in 1:(nrow(input) - 1)){
    for (y in 1:(ncol(input) - 1)){
      errors      <- input[x,y] - greys
      idx         <- which.min(abs(errors))
      output[x,y] <- greys[idx]
      te <- errors[idx]
      input[x,y+1]   <- input[x,y+1] + 7*te/16
      input[x+1,y-1] <- input[x+1,y-1] + 3*te/16
      input[x+1,y]   <- input[x+1,y] + 5*te/16
      input[x+1,y+1] <- input[x+1,y+1] + 1*te/16
    }
  }
  return(output)
}

### Input files
mapp    <- readPNG('Mapp_semiflat.png')
lucia   <- readPNG('Lucia_semiflat.png')
message <- readPNG('message.png')

### Tables of all overlay subpixel patterns and their densities
total    <- choose(4,1) + choose(4,2) + choose(4,3) + choose(4,4) + choose(4,0)
density  <- rep(NA, total)
patterns <- list()
counter  <- 0
for (ndark in 0:4){
  pp <- combn(1:4, ndark) 
  for (p in 1:ncol(pp)){
    counter             <- counter + 1
    patterns[[counter]] <- pp[,p] 
    density[counter]    <- length(pp[,p])
  }  
}

### Table: density of pattern combinations 
combined.density <- matrix(NA, nrow=counter, ncol=counter)
for (cx in 1:counter){
  for (cy in 1:counter){
      combined.density[cx, cy] <- length(unique(c(patterns[[cx]], patterns[[cy]])))
  }
}

### Prepare images by reducing and adjusting grayscale levels ("posterizing")
mapp4    <- floydSteinberg((1-mapp)*4, 0:4) 
lucia4   <- floydSteinberg((1-lucia)*4, 0:4)
message4 <- message + 3

### Matching subpixel patterns with minimal error
mapp.output <- matrix(NA, nrow=200, ncol=200)
lucia.output <- matrix(NA, nrow=200, ncol=200)
for (x in 1:200){
  for (y in 1:200){
    
    candidates <- which(combined.density==message4[x,y], arr.ind=TRUE)  
    errors     <- (mapp4[x,y]-density[candidates[,1]])**2 +
    		       (lucia4[x,y]-density[candidates[,2]])**2
    		   
    target <- which(errors==min(errors))
    if (length(target) > 1)
      target <- sample(target, 1)
    mapp.output[x,y]  <- candidates[target,1]
    lucia.output[x,y] <- candidates[target,2]
  }
}

image(mapp.output)
image(lucia.output)

### Generate dot matrix output from subpixel pattern choices
dottify <- function(input){
  output <- matrix(1, nrow=400, ncol=400)
  for (x in 1:200){
    for (y in 1:200){
      pattern <- patterns[[input[x,y]]]
      for (dot in pattern){
        coordsx <- (x-1)*2 + 1 + (dot-1)%%2
        coordsy <- (y-1)*2 + 1 + (dot-1)%/%2
        output[coordsx, coordsy] <- 0
      }
    }
  }
  return (output)
}

mapp.dot  <- dottify(mapp.output)
lucia.dot <- dottify(lucia.output)
image(t(apply(mapp.dot, 2, rev)), col=c("black", "white"), breaks=c(0, 0.5, 1), axes=F)
image(t(apply(lucia.dot, 2, rev)), col=c("black", "white"), breaks=c(0, 0.5, 1), axes=F)

### Test the overlay
image(t(apply(mapp.dot*lucia.dot, 2, rev)), col=c("white", "black"), breaks=c(0, 0.5, 1))