function response = harris (image, charScale, intScale)
  
  if (size(image,3)>1)
    image=rgb2gray(image);
  end;
  image=im2double(image);
  
% response = harris (image, charScale, intScale)
  neighborhoodSigma = eye (2);
  if (nargin < 4) intScale = 2*charScale; end;

  % vytvoreni rozmazavaci diferencni masky
%  assert (neighborhoodSigma (1, 2) == neighborhoodSigma (2, 1));
  normalizedSigma2x2 = neighborhoodSigma ./ sqrt (det (neighborhoodSigma));

  blurMask = makeGauss (normalizedSigma2x2 * charScale * charScale);
  blurredImage = filter2 (blurMask, image, 'same');

  [dxMask dyMask] = createDiffMasks (normalizedSigma2x2 * charScale * charScale);

  dx = conv2 (blurredImage, dxMask, 'same');
  dy = conv2 (blurredImage, dyMask, 'same');

  integrationMask = makeGauss (normalizedSigma2x2 * intScale * intScale);

  m11 = filter2 (integrationMask, dy.^2, 'same');
  m1221 = filter2 (integrationMask, dx.*dy, 'same');
  m22 = filter2 (integrationMask, dx.^2, 'same');

  % vypocet mensiho z vl. cisel
  b = m11 + m22;
  d = abs (sqrt (b .* b - 4 * (m11 .* m22 - m1221 .* m1221)));
  response = (b - d) / 2;

return



function mask = makeGauss (sigma2x2);

  l = ceil (sqrt (sqrt (prod (eig (sigma2x2)))));
  sizex = round (8*l) + 1; sizey = round (8*l) + 1;

  [x,y] = ind2sub ([sizex sizey], 1 : sizex*sizey);
  X = [x; y];
  clear x; clear y;

  X = X - round (sizex/2);
  mask = (1/(2*pi*sqrt(det(sigma2x2))))*exp(-sum(X'*inv(sigma2x2).*X',2)/2);
  mask = reshape (mask, sizex, sizey);

  if (sigma2x2 (2, 1) ~= 0) % symetricky gausy neorezat => separabilni maska => rychlejsi konvoluce
      maxBorderVal = max ([max(mask (1, 1:end)), max(mask (1:end, 1))]);
      mask (mask <= maxBorderVal) = 0;
  end

return

function [dxMask, dyMask, laplacian] = createDiffMasks (sigma)

  % smery derivaci (pro diagonalni sigmu - parcialni derivace)
  trMatrix = chol (sigma);
  dirX = [0 1] * trMatrix;
  dirY = [1 0] * trMatrix;

  dxMask = zeros (2*(max (floor (abs (dirX))) + 1) + 1); dxCenter = ceil (size (dxMask)/2);
  dyMask = zeros (2*(max (floor (abs (dirY))) + 1) + 1); dyCenter = ceil (size (dyMask)/2);

  xy2 = dirX (1) - floor (dirX (1));
  xy1 = 1 - xy2;
  xx2 = dirX (2) - floor (dirX (2));
  xx1 = 1 - xx2;

  yy2 = dirY (1) - floor (dirY (1));
  yy1 = 1 - yy2;
  yx2 = dirY (2) - floor (dirY (2));
  yx1 = 1 - yx2;

  dxMask (dxCenter (1) + floor (dirX (1)) : dxCenter (1) + floor (dirX (1)) + 1, dxCenter (2) + floor (dirX (2)) : dxCenter (2) + floor (dirX (2)) + 1) ...
   = [xx1*xy1, xx2*xy1; xx1*xy2, xx2*xy2];
  dxMask (dxCenter (1) - floor (dirX (1)) - 1 : dxCenter (1) - floor (dirX (1)), dxCenter (2) - floor (dirX (2)) - 1 : dxCenter (2) - floor (dirX (2))) ...
   = -[xx2*xy2, xx1*xy2; xx2*xy1, xx1*xy1];

  dyMask (dyCenter (1) + floor (dirY (1)) : dyCenter (1) + floor (dirY (1)) + 1, dyCenter (2) + floor (dirY (2)) : dyCenter (2) + floor (dirY (2)) + 1) ...
   = [yx1*yy1, yx2*yy1; yx1*yy2, yx2*yy2];
  dyMask (dyCenter (1) - floor (dirY (1)) - 1 : dyCenter (1) - floor (dirY (1)), dyCenter (2) - floor (dirY (2)) - 1 : dyCenter (2) - floor (dirY (2))) ...
   = -[yx2*yy2, yx1*yy2; yx2*yy1, yx1*yy1];

  dxMask (dxCenter (1), dxCenter (2)) = 0;
  dyMask (dyCenter (1), dyCenter (2)) = 0;

  lSize = max (size (dxMask), size (dyMask));
  stub = zeros (lSize); stub (ceil (lSize (1)/2), ceil (lSize (2)/2)) = 1;
  laplacian = -filter2 (abs (dxMask), stub,  'same') - filter2 (abs (dyMask), stub, 'same');
  laplacian = laplacian / -sum (sum (laplacian));
  laplacian (ceil (lSize (1)/2), ceil (lSize (2)/2)) = 1;

return