theonlineoasis

Image forensics – CFA interpolation detection

Light reaching digital camera sensors (CCDs) is filtered by a colour filter array (CFA), or Bayer array, which allows each sensor element to measure the intensity of either red, green or blue light. The monochrome image captured by the CCD is converted into a full-colour image by interpolating the missing pair of RGB values at each sample position.

My slides on CFA interpolation detection include a more thorough introduction to the topic.

Code samples on this website come without any guarantees – use them at your own risk! I am grateful for comments or patches sent to my email address: first name.last name@cl.cam.ac.uk.

#!/usr/bin/python

import os
import sys

# dcraw command line
#   Bilinear interpolation
#     dcraw -q 0
#   Variable numbers of gradients
#     dcraw -q 1
#   Patterned pixel grouping
#     dcraw -q 2
#   Adaptive homogeneity-directed
#     dcraw -q 3
#   Bilinear with 1-pass median filter
#     dcraw -q 0 -m 1
#

# http://www.pythonware.com/library/pil/handbook/index.htm
from PIL import Image

import matplotlib.pyplot
import numpy
import scipy.signal

if len(sys.argv) != 2:
  print 'Usage:\n\t%s image_file' % (sys.argv[0],)
  sys.exit(0)
image_file_name = sys.argv[1]

# Load the green channel of the image
green_channel = Image.open(image_file_name).split()[1]
(width, height) = green_channel.size
samples_g = numpy.reshape(numpy.fromstring(green_channel.tostring(),
                                           dtype=numpy.uint8),
                          (height, width))

# High pass filter
laplace_filter = \
  numpy.matrix([[ 0,  1,  0],
                [ 1, -4,  1],
                [ 0,  1,  0]])
samples_g = scipy.signal.convolve(samples_g, laplace_filter, 'valid', False)
(width, height) = samples_g.shape

# Find the sum of all the diagonals in the filtered green channel
diagonal_sums = []
for y in range(-width + 1, height):
  diagonal_sum = 0.0
  for x in range(max(0, -y), min(width, height - y)):
    diagonal_sum += abs(samples_g[(x, y + x)])
  diagonal_sums.append(float(diagonal_sum) / 
                       float(min(width, height - y) - max(0, -y)))

# Show the DFT of the signal, with a log scale on the Y axis
matplotlib.pyplot.gca().set_yscale('log')
matplotlib.pyplot.plot(abs(numpy.fft.fft(diagonal_sums)))
matplotlib.pyplot.show()