PCA_Method.py

"""
ECE 520 Project - Image Forgery Detection
PCA Forgery Detection
Travis Miller
Narasimha N. Chakravarty
"""
from numpy import *
from scipy import stats
from math import floor
from math import sqrt

"""Reads in image to tested"""
image_file = open ( 'Forge_2.txt', 'r' )
image = image_file.read()
image_file.close() #test

"""Resizes data from txt file to be the dimensions of the image"""
image = image.split()
image = np.reshape(image, (512, 512))
image = image.transpose()
im_size = np.shape(image)

"""Variables used in this program"""
blocks = zeros( (255025,64) ) #Matrix U with image blocks
blocks_norm = zeros( (255025,64) ) #Image Blocks with mean subtracted
blocks_cov = [] #Covariance Matrix of Image Blocks
offsets= zeros( (2295144,4))
new_offsets= zeros( (2295144,3))
offset_save = []
max_freq = 0
a = []
k = 0

"""Creates the image for the blocks"""
for i in range(im_size[0] - 7):
for j in range(im_size[1] - 7):
blk = array(image[i:i+8, j:j+8], dtype = single) #Takes out a 8x8 block out of the image
blk = np.reshape(blk, (1, 64)) #Forms the blocks into a vector
blocks[k] = blk #Stores the block vector into the matrix
k += 1

"""Calculation of Covariance and Eigenvectors"""
blk_avg = np.mean(blocks, axis=0) #Finds the mean of the columns
blk_avg = np.reshape(blk_avg, (1, 64)) #reshapes the matrix of the mean 
blocks_norm = blocks - blk_avg #subtracts out the mean from the columns 
blocks_norm_t = blocks_norm.transpose()
blocks_cov = np.dot(blocks_norm_t, blocks_norm)/255024 #Finds the covariance of the blocks matrix
w, v = linalg.eig(blocks_cov) #Finds the eigenvalues and eignevector of the covariance matrix 
pca_v = v[: , 0:4] #Creates the Principal Component Matrix

"""Tranforms the Block Vectors into new block representation"""
new_blocks = np.dot(blocks_norm, pca_v) #Multiplies block vector with 
new_blocks_t = new_blocks.transpose() #Transpose 
new_blocks = np.around(new_blocks) #Round the values to nearest integer

"""Attaches coordinate to every block"""
position = zeros( (255025,1) ) 
x_cor = zeros( (255025,1) )
new_position = zeros( (255025,1) )
y_cor = zeros( (255025,1) )

for i in range(255025):
position[i] = i + 1 #Makes vector for of matrix positions

new_blocks = np.hstack((new_blocks, position)) #Attaches matrix posisiton to rows

"""Sorts the matrix with new representation vectors"""
ind = np.lexsort((new_blocks[:,2], new_blocks[:,1], new_blocks[:,0])) #Sorts rows based of column values
sorted_blks = new_blocks[ind] #Rearranges matrix with new indice

new_position[:,0] = sorted_blks[:,4] #Takes the block coordinates of the 


"""Finds x and y corrdinate of row"""
k=0
for i in range(255025):
point = new_position[i]
x_cor[i] = (point - 1) % 505 #Finds x coordinate of block in image
y_cor[i] = np.floor((point - 1) / 505) #Finds y coordinate of block in image
new_position = np.hstack((x_cor, y_cor,new_position)) #Attaches y and x coordinate to row

"""Calculating neighboring row offsets"""
for i in range(255025 - 9):
for j in range(1,10):
if new_position[i][0] - new_position[i+j][0] > 0: 
offsets[k][0] = new_position[i][0] - new_position[i+j][0]
offsets[k][1] = new_position[i][1] - new_position[i+j][1]
elif new_position[i][0] - new_position[i+j][0] < 0:
offsets[k][0] = new_position[i+j][0] - new_position[i][0]
offsets[k][1] = new_position[i][1] - new_position[i+j][1]
else:
offsets[k][0] = 0
offsets[k][1] = abs(new_position[i][1] - new_position[i+j][1])
offsets[k][2] = new_position[i][0]
offsets[k][3] = new_position[i][1]
k += 1

"""Searches for most occuring offset"""
ind = np.lexsort((offsets[:,1],offsets[:,0])) #Sorts for most occuring offset
offsets = offsets[ind] #Rearrenges offsets in order of most frequent
offsets = offsets.astype(int) #Changes offset type to integer
off_group = offsets[:,0] * 1000 + offsets[:,1] #Creates Offset to one number to count frequency
freq_count = np.bincount(off_group)

"""Checks if the most occuring offset mets the threshold and is far enough apart"""
for i in range(len(freq_count)):
x = i/1000
y = i%1000
x_2 = 0
y_2 = 0
if x > 0:
x_2 = long(x * x) 
y_2 = long(y*y)
if y > 500:
x += 1
x_2 = long(x)
y = 1000 - y
y_2 = long(y)
y_2 = y_2 * y_2
y = -y
if freq_count[i] > 320 and sqrt(x_2 + y_2)>25: #Threshold and magnitude check
print i 
print freq_count[i]
print ""
if max_freq < freq_count[i]:
max_freq = freq_count[i] 
max_cor = x*1000+y #Stores matrix row location of most occuring offset
copy_x = x #Stores x offset
copy_y = y #Stores y offset
new_offsets[:,0] = off_group
new_offsets[:,1] = offsets[:,2]
new_offsets[:,2] = offsets[:,3]

"""Outputs coordinates of the forged region for visual display"""
condition = np.where( new_offsets == max_cor ) 
output_cor = new_offsets[condition[0]]
output_corx = output_cor[:,1]
output_cory = output_cor[:,2]
offset_save.append(copy_x)
offset_save.append(copy_y) 
np.savetxt('x_cor.txt', output_corx) #Outputs x coordinates of forged pixels
np.savetxt('y_cor.txt', output_cory) #Outputs y coordinates of forged pixels
np.savetxt('offset.txt', offset_save) #Outputs the offset value between pixels