Code

#-----------------------------------------------------------------------------
# Read a jpg file using image magick and scale it and convert to greyscale
#-----------------------------------------------------------------------------
filename <- 'data/mona.jpg'
im <- magick::image_read(filename) %>%
  magick::image_convert(type='grayscale') %>%
  magick::image_scale(geometry="75%") %>%
  magick::image_flip()

#-----------------------------------------------------------------------------
# Extract just the numeric matrix representing the image.
#-----------------------------------------------------------------------------
m <- magick::as_EBImage(im)@.Data


#-----------------------------------------------------------------------------
# Set up plot configuration
#-----------------------------------------------------------------------------
par(mar = c(0, 0, 0, 0))

xlim          <- c(0, dim(m)[1])
ylim          <- c(0, dim(m)[2])
width         <- 4
height        <- 6
dpi           <- 200


#-----------------------------------------------------------------------------
# how many objects in each frame?
#-----------------------------------------------------------------------------
nobjects_ <- c(
  seq(    2,     20 - 1,     1),
  seq(   20,    200 - 1,    10),
  seq(  200,   2000 - 1,   100),
  seq( 2000,  50000 - 1,  1000),
  seq(50000,  90000 - 1,  5000),
  seq(90000, 100000    ,   500)
)


#-----------------------------------------------------------------------------
# What alpha level for each frame?  The more objects there are,
# the lower the alpha to ensure the overlap looks good.
#-----------------------------------------------------------------------------
col_ <- nobjects_ %>% purrr::map_dbl(~1 - (log10(.x)/5 - log10(1.9)/5)^0.2) %>%
  purrr::map_chr(~gray(0, .x))


#-----------------------------------------------------------------------------
# Sample all the points necessary for the maximum nobjects and include
# a 10x factor since I know I'm going to filter out lots of the generated points
#-----------------------------------------------------------------------------
all_indices <- sample(seq(m), size = 10 * max(nobjects_ + 10), prob = (1-m)^2, replace = TRUE)
all_points  <- arrayInd(all_indices, .dim=dim(m))


#-----------------------------------------------------------------------------
# Now filter the successive points so that it's biased towards shorter 
# line segments
#-----------------------------------------------------------------------------
max_dist <- sqrt(nrow(m)^2 + ncol(m)^2)

pts <- as.data.frame(all_points) %>% 
  as.tbl() %>%
  set_names(c('x', 'y')) %>%
  mutate(
    val  = 1 - m[all_indices]
  )

for (i in 1:20) {
  pts %<>% mutate(
    dist = sqrt((x - lead(x))^2 + (y - lead(y))^2),
    sdist = dist/max_dist,
    keep  = val > 2.5 * sdist
  )
  
  pts %<>% filter(
    is.na(keep) | keep
  )
}

  
  
#-----------------------------------------------------------------------------
# Draw each frame
#-----------------------------------------------------------------------------
for (i in seq(nobjects_)) {
  nobjects <- nobjects_[i]
  col      <- col_[i]
  
  plot_points <- pts[seq(1+nobjects),c('x', 'y')]
  
  #-----------------------------------------------------------------------------
  # Plot lines and save
  #-----------------------------------------------------------------------------
  plot_filename <- sprintf("plot//mona2/plot_%06i.png", nobjects)
  withr::with_png(plot_filename, width=width*dpi, height=height*dpi, {
    plot(1, type="n", xlab="", ylab="", xlim=xlim, ylim=ylim, asp=1, ann=FALSE,
         axes = FALSE, frame.plot = FALSE)
    
    lines(plot_points, col=col)
  })
}