function violaK = violaKernel(X1, X2, kernelParams)
%Compute the a linear kernel using the viola features on a set of images 
%This is going to be a bit slow 

if exist(kernelParams.kernelFunctionName, 'file') ~= 2
    error('Kernel function %s does not exist\n', kernelParams.kernelFunctionName); 
end

numX1Examples = size(X1, 1); 
numX2Examples = size(X2, 1); 

%We want X1 to be smaller to improve efficiency 
if numX1Examples > numX2Examples 
    tempX = X1; 
    X1 = X2; 
    X2 = tempX;
    numX1Examples = size(X1, 1);
    numX2Examples = size(X2, 1);
    clear tempX;
    pack;
    swap = 1;
else
    swap = 0;
end

kernelFunctionName = kernelParams.kernelFunctionName;
width = kernelParams.width;
height = kernelParams.height; 
X1cacheSize = min(kernelParams.kernelCacheSize, numX1Examples);
X2cacheSize = min(kernelParams.kernelCacheSize, numX2Examples);

numX1Chunks = ceil(numX1Examples/X1cacheSize);
numX2Chunks = ceil(numX2Examples/X2cacheSize);

%Free up some memory 
pack; 
violaK = zeros(numX1Examples, numX2Examples); 
ex1 = violaFeatures(reshape(X1(1, :), width, height));
%ex1 = X1(1, :); 
numFeatures = length(ex1);

for i=1:numX1Chunks
    fprintf('Calculating kernel chunk %d of %d\n', i, numX1Chunks);
    
    startX1Index = (i-1)*X1cacheSize + 1;
    endX1Index = min(i*X1cacheSize, numX1Examples);

    X1Cache = zeros(endX1Index-startX1Index, numFeatures); 
    
    %First fill up X1Cache
    for k=1:(endX1Index-startX1Index)+1
        %X1Cache(k, :) = X1(startX1Index+k-1, :);
        X1Cache(k, :) = violaFeatures(reshape(X1(startX1Index+k-1, :), width, height));
        X1Cache(k, :) = X1Cache(k, :)/norm(X1Cache(k, :)); 
    end 
   
    
    for j=1:numX2Chunks
        startX2Index = (j-1)*X2cacheSize + 1;
        endX2Index = min(j*X2cacheSize, numX2Examples);

        X2Cache = zeros(endX2Index-startX2Index, numFeatures); 
        
        %Fill up the X2Cache 
        for k=1:(endX2Index-startX2Index)+1
            %X2Cache(k, :) = X2(startX2Index+k-1, :);
            X2Cache(k, :) = violaFeatures(reshape(X2(startX2Index+k-1, :), width, height));
            X2Cache(k, :) = X2Cache(k, :)/norm(X2Cache(k, :)); 
        end
        
        violaK(startX1Index:endX1Index, startX2Index:endX2Index) = feval(kernelFunctionName, X1Cache, X2Cache, kernelParams);
    end
end

clear X1Cache X2Cache; 
pack; 

if swap == 1 
    violaK = violaK'; 
end