What’s the problem?

There was no problem until yesterday’s supermoon, when I decided to look up at the sky sometime at around 9PM.

It was a perfectly clear night for some moon gazing!. I then pulled out my camera and started taking some pics of the moon.

I continued taking pictures every few minutes, without really noticing what was going on. It wasn’t until I started reviewing the first pictures and doing some Google searching that the eclipse became evident. I was right in the middle of it.

Half way thru the eclipse, I started to realize (and regret) that I wasn’t being very methodical about taking my pictures. All of the sudden I started to see how, if I had more pictures, I could put together a nice little time-lapse video of this astronomical phenomenon.

The problem is obvious, I’m missing data!.

Once I realized this, I started to thinker about how cool it would be to try to use some linear algebra in order to extrapolate the missing pictures.

The Approach

Let’s code.

The pictures: Basic Photo Editing

The idea is to accomplish the goal of estimating the missing pictures with minimum photo editing.

As you can see from below, the pictures were re-sized and zoomed so that the moon appears in the center of the picture with a consistent size.

All pictures were re-sized to 500x500 pixels using the sips command:

sips IMG_0* -z 500 500

These are the 19 pictures taken during the eclipse:

supermoon 2015 pics

Obtain the precise time when each picture was taken

Most digital cameras embed plenty of meta-data about the pictures as EXIF Attributes.

In order to extract these attributes, I used the ExifTool by Phil Harvey

exiftool *.jpg |grep "\(Create Date\|File Name\)"

Here is a snippet of the output:

File Name                       : IMG_0084.jpg
Create Date                     : 2015:09:27 21:24:53.00
File Name                       : IMG_0085.jpg
Create Date                     : 2015:09:27 21:25:47.10
File Name                       : IMG_0086.jpg
Create Date                     : 2015:09:27 21:26:32.60
File Name                       : IMG_0087.jpg
Create Date                     : 2015:09:27 21:27:51.60
.
.
.
.
.
File Name                       : IMG_0100.jpg
Create Date                     : 2015:09:27 22:07:08.40
File Name                       : IMG_0101.jpg
Create Date                     : 2015:09:27 22:15:00.10
File Name                       : IMG_0102.jpg
Create Date                     : 2015:09:27 22:37:17.00

With this list, some copy and paste and excel formulas, we end up with the following R data-frame:

library(knitr)
pic_no=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19)
filename = c('IMG_0084.jpg','IMG_0085.jpg','IMG_0086.jpg','IMG_0087.jpg','IMG_0088.jpg','IMG_0089.jpg','IMG_0090.jpg','IMG_0091.jpg','IMG_0092.jpg','IMG_0093.jpg','IMG_0094.jpg','IMG_0095.jpg','IMG_0096.jpg','IMG_0097.jpg','IMG_0098.jpg','IMG_0099.jpg','IMG_0100.jpg','IMG_0101.jpg','IMG_0102.jpg')
elapsed_since_first_pic_sec= c(0,54,99,178,288,329,380,570,946,1106,1439,1592,1726,1851,2015,2198,2535,3007,4344)
elapsed_since_first_pic_min =c(0,0,1,2,4,5,6,9,15,18,23,26,28,30,33,36,42,50,72)
exif_meta_data <- data.frame(pic_no,filename,elapsed_since_first_pic_sec,elapsed_since_first_pic_min, stringsAsFactors=FALSE)

kable(exif_meta_data)
pic_no filename elapsed_since_first_pic_sec elapsed_since_first_pic_min
1 IMG_0084.jpg 0 0
2 IMG_0085.jpg 54 0
3 IMG_0086.jpg 99 1
4 IMG_0087.jpg 178 2
5 IMG_0088.jpg 288 4
6 IMG_0089.jpg 329 5
7 IMG_0090.jpg 380 6
8 IMG_0091.jpg 570 9
9 IMG_0092.jpg 946 15
10 IMG_0093.jpg 1106 18
11 IMG_0094.jpg 1439 23
12 IMG_0095.jpg 1592 26
13 IMG_0096.jpg 1726 28
14 IMG_0097.jpg 1851 30
15 IMG_0098.jpg 2015 33
16 IMG_0099.jpg 2198 36
17 IMG_0100.jpg 2535 42
18 IMG_0101.jpg 3007 50
19 IMG_0102.jpg 4344 72

Denoise the pictures using SVD factorization

Let’s use SVD Factorization in order to compress the images and reduce random noise.

First, let’s pick and image and find a good compression factor that we can apply to all the images in the data-set:

library(jpeg)
library(ripa)
## Loading required package: tcltk
## Loading required package: parallel
library(RCurl)
## Loading required package: bitops
#select first image in dataset
pic_url <- paste("https://raw.githubusercontent.com/rmalarc/supermoon2015/master/", exif_meta_data[1,"filename"],sep="")

pic_data <- readJPEG(getURLContent(pic_url))

pic_matrix <- imagematrix(pic_data, type = "grey")


pic_matrix_svd <- svd(pic_matrix)
d <- diag(pic_matrix_svd$d)
u <- pic_matrix_svd$u
v <- pic_matrix_svd$v
plot(1:length(pic_matrix_svd$d),pic_matrix_svd$d)

Let’s try a factor of 100 to de-noise and compress our images.

# de-noised and compressed pic
depth <- 100
us <- as.matrix(u[, 1:depth])
vs <- as.matrix(v[, 1:depth])
ds <- as.matrix(d[1:depth, 1:depth])

pic_denoised <- us %*% ds %*% t(vs)

pic_matrix_denoised <- imagematrix(pic_denoised, type = "grey")

par(mfrow=c(1,2)) 

# original picture
plot(pic_matrix, useRaster = TRUE)

# Denoised picture
plot(pic_matrix_denoised, useRaster = TRUE)

A factor of 100 seems visually suitable for these set of pictures.

Can we estimate a missing picture using matrix linear combination?

The idea is to obtain the matrix \(A_{j}\) that corresponds to the picture that is in between \(P_{i}\) and \(P_{k}\) from the corresponding matrices \(A_{i}\) and \(A_{k}\). The strategy is to do a linear combination by multiplying each matrix for a constant factor and adding up the two as follows:

\(A_{j} = c_{i}A_{i} + c_{k}A_{k}\)

where:

\(c_{k} = \frac{j-i}{k-i}\)

and

\(c_{i} = 1-c_{k}\)

As we can see here, let’s try with image 1 and 9 and see if we can get anything that looks like image 6

kable(exif_meta_data[c(1,6,9),])
pic_no filename elapsed_since_first_pic_sec elapsed_since_first_pic_min
1 1 IMG_0084.jpg 0 0
6 6 IMG_0089.jpg 329 5
9 9 IMG_0092.jpg 946 15

Let’s code:

plotMatrix <- function(matrix){
  # It takes an image as a matrix and plots it out
  plot(imagematrix(matrix, type = "grey"), useRaster = TRUE) 
}

recomposeMatrix <- function(matrix_svd){
  # reassembles a previously decomposed image from getImageSVDProjections
  return (matrix_svd$us %*% matrix_svd$sigma_k %*% matrix_svd$vs_t)
}

getImageSVDProjections <- function(exif_meta_data,image_no){
  # This function takes the dataframe containing the list of key images and an image number
  # It then downloads the image in question, decomposes it using svd and projects it to a 
  # lower dimensional space (preset by the depth constant), returning a list with the
  # u,v and sigma_k components
  
  depth <- 100
  pic_url <- paste("https://raw.githubusercontent.com/rmalarc/supermoon2015/master/"
                   , exif_meta_data[image_no,"filename"]
                   ,sep=""
                   )
  pic_data <- readJPEG(getURLContent(pic_url))

  pic_matrix <- imagematrix(pic_data, type = "grey")

  pic_matrix_svd <- svd(pic_matrix)
  d <- diag(pic_matrix_svd$d)
  u <- pic_matrix_svd$u
  v <- pic_matrix_svd$v
  
  us <- as.matrix(u[, 1:depth])
  vs_t <- t(as.matrix(v[, 1:depth]))
  sigma_k <- as.matrix(d[1:depth, 1:depth])

  return(list(us=us,sigma_k=sigma_k,vs_t=vs_t))
}

# image i
image_1 <- getImageSVDProjections(exif_meta_data,1)
# image k
image_9 <- getImageSVDProjections(exif_meta_data,9)

# our CONTROL image j
image_6 <- getImageSVDProjections(exif_meta_data,6)


# i and k are the seconds elapsed since first pic of each pic
i <- exif_meta_data[1,"elapsed_since_first_pic_sec"]
j <- exif_meta_data[6,"elapsed_since_first_pic_sec"]
k <- exif_meta_data[9,"elapsed_since_first_pic_sec"]

c_k <- (j-i)/(k-i)
c_i <- 1-c_k

image_1_matrix <- recomposeMatrix(image_1)
image_6_matrix <- recomposeMatrix(image_6)
image_9_matrix <- recomposeMatrix(image_9)
image_j <- c_i*image_1_matrix + c_k*image_9_matrix


#estimated image
par(mfrow=c(1,3)) 
plotMatrix(image_1_matrix)
plotMatrix(image_j)
plotMatrix(image_9_matrix)

As we can see above, our estimated image results in a nice cross-fade between image 1 and 9. This raises the question, how does this image J compare to our control image #6?

par(mfrow=c(1,2)) 
plotMatrix(image_6_matrix)
plotMatrix(image_j)

Close enough I would say!

Estimate the missing frames (1 per minute at least) by linear combination of matrices

So, if we look at the table with the pictures taken, we have a span of 72 minutes with pictures. The idea is to generate one picture per minute, create an animation and save it as an animated gif.

Let’s code:

library(animation)

LoadImageMatrix <- function(exif_meta_data,pic_no){
  # This function calls the necesary functions that loads an image, denoises it and 
  # returns the corresponding image matrix

  image <- getImageSVDProjections(exif_meta_data,pic_no)
  return(recomposeMatrix(image))
}

par(mfrow=c(1,1)) 
saveGIF({
  prev_i_pic_no <- -1
  
  # Set the time range I want to cover. We have 72 min worth of data, let's shoot for
  # generating 4 frames per minute
  range <- seq(from=0, to=72, by=0.25)
  range <- range*60 # turn to seconds
  
  # initialize variables
  image_i_matrix <- NA
  image_k_matrix <- NA
  
  i <- NA
  k <- NA

  # iterate thru every desired time lapse in the range
  for (j in range) {
    # find the nearest key image i
    
    # the first time we enter the loop, image i is initialized picture # 1
    i_pic_no <- 1
    
    # once we iterate further, we need to look for the closest images i and k
    if (j > 0 ){
      i_pic_no<-max(exif_meta_data[exif_meta_data$elapsed_since_first_pic_sec<j,"pic_no"])
    }
    
    # check and see if we should be using a new key image
    if(i_pic_no != prev_i_pic_no) {
      # We should be using a new set of key images i and k, load them up
      k_pic_no<-i_pic_no + 1
      image_i_matrix <- LoadImageMatrix(exif_meta_data,i_pic_no)
      image_k_matrix <- LoadImageMatrix(exif_meta_data,k_pic_no)
      
      # i and k are the seconds elapsed since first pic of each pic
      i <- exif_meta_data[i_pic_no,"elapsed_since_first_pic_sec"]
      k <- exif_meta_data[k_pic_no,"elapsed_since_first_pic_sec"]
    }
    prev_i_pic_no <- i_pic_no
    
    # Do the linear combination of the key images i and k in order to generate image j
    c_k <- (j-i)/(k-i)
    c_i <- 1-c_k
    image_j_matrix <- c_i*image_i_matrix + c_k*image_k_matrix
    
    # Plot image J
    plotMatrix(image_j_matrix)
    }
  }, interval = 0.25, movie.name = "supermoon_2015.gif", ani.width = 500, ani.height = 500)

The resulting animated gif is here

Converting the Animated GIF to a Video file

The animated gif was then converved to MP4 video format using Cloud Convert and uploaded to YouTube.

Conclusion