rm(list=ls(all=TRUE))  # clear all variables
graphics.off()  # clear all graphics
# Convolution
# Greg Francis
# PSY 628
# 16 August 2024

library(imager)
library(gsignal)

# Open image file
im <- load.image('Mountains.jpeg')

# convert to grayscale

grayMountains <- grayscale(im)

plot(grayMountains)

# Define convolution kernel (receptive field)

# On-center, off-surround
ksize=9 # Nicer to use an odd number
Kernal1 <-matrix(0, nrow= ksize, ncol= ksize) 
middle=floor(ksize/2)+1
sigma_c = 3
sigma_s = 8
for (x in 1: ksize){
	for(y in 1: ksize){
		
		val = 1/(sigma_c * sqrt(2*pi)) * exp(-((x-middle)^2 +(y-middle)^2 ) /(2*sigma_c*sigma_c)) - 1/(sigma_s * sqrt(2*pi)) * exp(-((x-middle)^2 +(y-middle)^2 ) /(2*sigma_s*sigma_s))
		Kernal1[x,y] = val
		
	}	
}


# Oriented filters (based on Difference of Offset Gaussians)
ksize=25 # Nicer to use an odd number
Kernal1 <-matrix(0, nrow= ksize, ncol= ksize)  
middle=floor(ksize/2)+1 # middle of kernal is (middle, middle)
theta = pi/2 # Orientation (in radians) 0 -> vertical
sigma = 4 # standard deviation of Gaussian envelope
offset= 0.5 # shifts middle in each direction along the orientaton direction

for (x in 1: ksize){
	for(y in 1: ksize){ 	
		
		xprime1 = x-(middle-offset*cos(theta)) 
		yprime1 = y-(middle-offset*sin(theta))
		xprime2 = x-(middle+ offset*cos(theta))
		yprime2 =  y-(middle+ offset*sin(theta))
			
		val = 1/(sigma * sqrt(2*pi)) * exp(-(xprime1^2 +yprime1^2 ) /(2*sigma*sigma)) - 1/(sigma * sqrt(2*pi)) * exp(-(xprime2^2 +yprime2^2 ) /(2*sigma*sigma))
		
		Kernal1[x,y] = val	
	}	
}

# Balance excitation (center) and inhibition (surround)
sumEx = sum(Kernal1[Kernal1>0])
sumIn = -sum(Kernal1[Kernal1<0])

Kernal1[Kernal1>0] = Kernal1[Kernal1>0]/sumEx

Kernal1[Kernal1<0] = Kernal1[Kernal1<0]/sumIn


dev.new()
image(Kernal1)
contour(Kernal1, add = TRUE)

convolvedImage <- conv2(grayMountains, Kernal1, "same")

# Set negative values to 0
convolvedImage[convolvedImage<0] = 0


# Convolved image seems to be flipped/rotated relative to original, this corrects
rotate <- function(x) t(apply(x, 2, rev))
# Flip across horizontal middle
temp <- apply(convolvedImage, 2, rev)
# rotate
convolvedImage<- rotate(rotate(temp))

dev.new()
image(convolvedImage, col=gray.colors(50))