IsmailPeters

Back to Psych 221 Projects 2013

Background

Camera Basics

A digital camera captures color images by recording data in different channels of the visible spectrum (red, green, blue). Most cameras accomplish this by using a single type of CCD or CMOS sensor at each pixel. This sensor measures the intensity of light that is focused on it, but cannot distinguish between colors. Therefore, a color filter is placed over each sensor, and the sensor only records the intensity of a specific color at that pixel. The arrangement of these color filters over an image is known as a Color Filter Array. Figure 1 shows a very common CFA known as the Bayer Array.

Figure 1: Bayer Array

In the Bayer Array, green color filters are arranged in a checkerboard pattern, whereas the red and blue color filters are arranged by alternating rows. The original reasoning for using this distribution of colors was to mimic the physiology of the human eye. During daylight vision, the luminance perception of the retina uses both L and M cones, which are more sensitive to green light. The red and blue filters control the sensitivity of the eye to chrominance.

CFA Interpolation

Since each color channel is only tallied at specific coordinates, the remaining pixels in that channel must be estimated in some way. There are many interpolation schemes - a few simple examples are the bilinear,bicubic, and smooth hue transition interpolations shown below. Other interpolation schemes, such as median filter, gradient-based, and adaptive color plane, become increasingly complex. A wide survey of interpolation methods is available at (REFERENCE), along with an analysis of the advantages and disadvantages of each.

Bilinear/Bicubic

A linear combination of the nearest N neighbors in any direction for N = 1 (bilinear) or N = 3 (bicubic). For a given matrix of known values R(x,y), the interpolated set of values r(x,y) can be computed:

$r(x,y)=\sum _{u,v=-N}^{N}h_{r}(u,v)R(x-u,y-v)$

where h_r is the 2-D kernel containing interpolation weights based on relative neighbors of a given pixel.

$h_{1}={\frac {1}{4}}{\begin{bmatrix}0&1&0\\1&4&1\\0&1&0\end{bmatrix}},h_{3}={\frac {1}{256}}{\begin{bmatrix}0&0&0&1&0&0&0\\0&0&-9&0&-9&0&0\\0&-9&0&81&0&-9&0\\1&0&81&256&81&0&1\\0&-9&0&81&0&-9&0\\0&0&-9&0&-9&0&0\\0&0&0&1&0&0&0\end{bmatrix}}$

Smooth Hue Transition

Operates on the assumption that the hue of a natural image (chrominance/luminance, or red/green and blue/green) changes smoothly in local regions of that image. In order to calculate this, the green channel is first interpolated using bilinear interpolation. A missing red point r(x,y) is calculated from known red values R(x,y) and green values G(x,y).

$r(x,y)=G(x,y)*{\frac {1}{2}}\left({\frac {R(x-1,y)}{G(x-1,y)}}+{\frac {R(x+1,y)}{G(x+1,y)}}\right)$

The sum shown above is valid for when the interpolation pixel has two adjacent pixels in the same row. The equation is altered slightly for two adjacent pixels in the same column, or when there are four adjacent pixels in the corners. The missing blue pixels are calculated in the same way.

Detecting Forgeries

Regardless of the type of interpolation used in a given image, the estimated pixels should exhibit a strong correlation, or dependence, on their surrounding pixels. If one can categorize each pixel as either correlated or independent with respect to other pixels, one should see a periodic pattern that mimics the CFA used to construct the image. Even if an altered image does not have any visual cues that point to its forgery, inspecting the correlation of pixels to one another can potentially expose which parts of an image were altered. It is important to note that, even if an image has been altered, the mentioned correlations may still be kept intact. In this way, detecting forgeries through CFA interpolation is one tool out of many used to detect forgeries.

Methods

Finding Correlations

The first step in detecting forgeries is categorizing the pixels in an image into two groups: those that are correlated with neighboring pixels and those that independent. Since the camera, CFA, an interpolation method are unknown for a given image, a robust algorithm must be used to iteratively deduce this categorization.

Expectation/Maximization

The EM algorithm has a wide variety of uses for working on a statistical set of data (REFERENCE). In the context of this problem, each step of the algorithm has a specific purpose. In the expectation step, the probability of each pixel in a given image being correlated to its neighbors is calculated based on the latest guess of interpolation coefficients. In the maximization step, new interpolation coefficients are calculated based on the probability map that reveals which pixels are indeed interpolated. Iterating through this process converges to a stable estimate of both the probability map and interpolation coefficients. If a pixel is correlated to neighboring pixels, it takes the form:

$f(x,y)=\sum _{u,v=-N}^{N}\alpha _{u,v}f(x+u,y+v)+n(x,y)$

where $\alpha _{u,v}$ are the interpolation coefficients and n(x,y) is the residual error of the linear model and ideally a Gaussian distribution with zero mean. The probability that a sample f(x,y) exists given that it is correlated to its neighbors is given by the equation:

$P(x,y)={\frac {1}{\sigma {\sqrt {2\pi }}}}exp\left[-{\frac {1}{2\sigma ^{2}}}\left(f(x,y)-\sum _{u,v=-N}^{N}\alpha _{u,v}f(x+u,y+v)\right)^{2}\right]$

where $\sigma$ is the standard deviation of this Gaussian distribution and is updated in the M-step. The new estimate of the interpolation coefficients is found by minimizing the quadratic error function:

$E(\alpha )=\sum _{x,y}w(x,y)\left(f(x,y)-\sum _{u,v=-N}^{N}\alpha _{u,v}f(x+u,y+v)\right)^{2}$

By setting the derivative of this equation equal to zero and rearranging terms, the following equation exists for each relative index (s,t) of the set interpolation coefficients $\alpha (u,v)$:

$\sum _{u,v=-N}^{N}\alpha _{u,v}\left(\sum _{x,y}w(x,y)f(x+s,y+t)f(x+u,y+v)\right)=\sum _{x,y}w(x,y)f(x+s,y+t)f(x,y)$

where w(x,y) = P(x,y). By performing this for each coefficient, a system of M equations and M unknowns, M being the number of interpolation coefficients, is set up. In this way, a new series of coefficients can be computed.
In order to make the computation optimized for MATLAB, the summations were converted into matrix multiplications. For example, for a given interpolation coefficient $\alpha _{s,t}$ , the w(x,y) , f(x+s,y+t) , and f(x,y) matrices can be written as row vectors of size 1 x N where N is the number of pixels under analysis. Thus the right hand side of the equation can be represented with matrix multiplations:

$\sum _{x,y}w(x,y)f(x+s,y+t)f(x,y)={\begin{bmatrix}w(x,y)\end{bmatrix}}.*{\begin{bmatrix}f(x+s,y+t)\end{bmatrix}}*{\begin{bmatrix}f(x,y)\end{bmatrix}}^{T}$

Similarly, the left hand side of the equation can be re-written. If all interpolation coefficients are taken into account, the entire system of equations can be represented:

${\begin{bmatrix}w(x,y)\\.\\.\\w(x,y)\end{bmatrix}}.*{\begin{bmatrix}f(x+1,y+0)\\f(x+0,y+1)\\.\\f(x+s,y+t)\end{bmatrix}}*{\begin{bmatrix}f(x+1,y+0)\\f(x+0,y+1)\\.\\f(x+u,y+v)\end{bmatrix}}^{T}*{\begin{bmatrix}\alpha _{1,0}\\\alpha _{0,1}\\.\\\alpha _{u,v}\end{bmatrix}}={\begin{bmatrix}w(x,y)\\.\\.\\w(x,y)\end{bmatrix}}.*{\begin{bmatrix}f(x+1,y+0)\\f(x+0,y+1)\\.\\f(x+s,y+t)\end{bmatrix}}*{\begin{bmatrix}f(x,y)\\f(x,y)\\.\\f(x,y)\end{bmatrix}}^{T}$

For clarity, the resulting sizes of these matrices is contained below.

${\begin{bmatrix}MxM\end{bmatrix}}*{\begin{bmatrix}\alpha _{1}\\\alpha _{2}\\.\\\alpha _{M}\end{bmatrix}}={\begin{bmatrix}Mx1\end{bmatrix}}$

In this way, a new set of coefficients can then be easily computed. This rearrangement of the data provided much faster computations, and it was very easy to add more coefficients to the system. The EM algorithm was allowed to run until the coefficients generated in successive iterations differed by an amount less than some arbitrarily small threshold.

2-D FFT

Once a probability map was generated from the EM algorithm, the next step was to determine if the periodic pattern of the CFA was present in the map. To do this, a two-dimensional FFT was performed on the probability map. The default 2-D FFT magnitude plot in MATLAB contains all of the low frequency information in the upper left corner, and the high frequency components along the left and bottom edges. In order to make the desired periodicity easier to observe, the FFT was up-sampled, high-pass filtered, and shifted so that the low frequency components are at the center of the image.

FFT BEFORE/AFTER PICTURE HERE

Similarity Measure

In order to build a classifier, a reliable quantity must be measured that is different by some threshold for real and fake images. The method chosen utilizes the similarities between the FFT of the probability map of an image and the FFT of an ideal synthetic grid. To generate a synthetic grid for a given color channel, any pixel that corresponds to that channel's color is marked with a 1 at that location; otherwise, all other locations contain a zero. In order to avoid the event that the probability map and synthetic grid are out of phase, which would yield a low similarity when compared directly, the magnitudes of the FFT results are compared. If $P(\omega _{x},\omega _{y})$ is the FFT of the probability map, and $S(\omega _{x},\omega _{y})$ is the FFT of the synthetic grid, then the similarity measurement M is:

$M=\sum _{\omega _{x},\omega _{y}}|P(\omega _{x},\omega _{y})|\cdot |S(\omega _{x},\omega _{y})|$

For windows in which the probability map is close to that expected, this similarity measure should remain constant. If the probability map distribution varies significantly, the similarity measure will deviate from the average value. This can be true for any channel.

Results

WE WON.

Retinotopic models in native space

Some text. Some analysis. Some figures.

Retinotopic models in individual subjects transformed into MNI space

Some text. Some analysis. Some figures.

Retinotopic models in group-averaged data on the MNI template brain

Some text. Some analysis. Some figures. Maybe some equations.

Equations

If you want to use equations, you can use the same formats that are use on wikipedia.
See wikimedia help on formulas for help.
This example of equation use is copied and pasted from wikipedia's article on the DFT.

The sequence of N complex numbers x₀, ..., x_N−1 is transformed into the sequence of N complex numbers X₀, ..., X_N−1 by the DFT according to the formula:

$X_{k}=\sum _{n=0}^{N-1}x_{n}e^{-{\frac {2\pi i}{N}}kn}\quad \quad k=0,\dots ,N-1$

where i is the imaginary unit and $e^{\frac {2\pi i}{N}}$ is a primitive N'th root of unity. (This expression can also be written in terms of a DFT matrix; when scaled appropriately it becomes a unitary matrix and the X_k can thus be viewed as coefficients of x in an orthonormal basis.)

The transform is sometimes denoted by the symbol ${\mathcal {F}}$, as in $\mathbf {X} ={\mathcal {F}}\left\{\mathbf {x} \right\}$ or ${\mathcal {F}}\left(\mathbf {x} \right)$ or ${\mathcal {F}}\mathbf {x}$.

The inverse discrete Fourier transform (IDFT) is given by

$x_{n}={\frac {1}{N}}\sum _{k=0}^{N-1}X_{k}e^{{\frac {2\pi i}{N}}kn}\quad \quad n=0,\dots ,N-1.$

Retinotopic models in group-averaged data projected back into native space

Some text. Some analysis. Some figures.

Conclusions

Here is where you say what your results mean.

References - Resources and related work

References

Software

Appendix I - Code and Data

Code

File:CodeFile.zip

Data

zip file with my data

Appendix II - Work partition (if a group project)

Brian and Bob gave the lectures. Jon mucked around on the wiki.