sigma = 0.5;
gauss = fspecial('gaussian', [5 5], sigma);
blur1 = imfilter(img, gauss, 'replicate');
dog1 = img - blur1;
%Next level
blur2 = imfilter(blur1, gauss, 'replicate');
dog2 = blur1 - blur2;

sigma = 0.5;
gauss1 = fspecial('gaussian', round([10*sigma 10*sigma]), sigma);
sigma = 1;
gauss2 = fspecial('gaussian', round([10*sigma 10*sigma]), sigma);
blur1 = imfilter(img, gauss1, 'replicate', 'same');
blur2 = imfilter(img, gauss2, 'replicate', 'same');
dog2 = blur1 - blur2;

%% Filter using DoG
stepsPerOctave = 5;
octaves = 4;
mult = nthroot(2,stepsPerOctave);

% Create blurry images
sigma = 0.5;
kernelSize = [10*sigma*2^(octaves),10*sigma*2^(octaves)]
for k = 1:octaves*stepsPerOctave+1
    disp(['Sigma is ' num2str(sigma)]);
    gauss = fspecial('gaussian', kernelSize, sigma);
    blur(:,:,k) = imfilter(I, gauss, 'replicate', 'same');
    imagesc(blur(:,:,k)); colorbar; title(['Gaussian ' num2str(k)]); pause;
    sigma = sigma * mult;
end

% Create DoG
for k = 1:octaves*stepsPerOctave
    dog(:,:,k) = blur(:,:,k+1) - blur(:,:,k);
    imagesc(dog(:,:,k)); colorbar; title(['DoG ' num2str(k)]); pause;
end