function [c, M] = rigid2Dreg(im1,im2)

  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % calculate partial derivatives of images
  %

  M = 5; Mmax = (M-1)/2;
  s = size(im1); sx = s(1)-M+1; sy = s(2)-M+1;
  %  x = 1:sx; y = 1:sy;
  x = -sx/2 + 1:sx/2; y = -sy/2 + 1:sy/2;

  % build filters (both filtdiff & filtsmooth smooth with a hamming
  % window, filtdiff also does derivative)
  Md = 9; n = -(Md-1)/2:(Md-1)/2;
  diff = cos(pi*n)./n - sin(pi*n)./(pi*n.^2);  diff(Md/2) = 0;
  smooth = hamming(M)'; smooth = smooth/sum(smooth);
  hdiff = conv(diff,smooth);
  hdiff = hdiff((Md-1)/2+1:(Md-1)/2+M);
  filtdiff = hdiff'*smooth;
  filtsmooth = smooth'*smooth;

  % calculate partial derivatives of images
  fxim1 = conv2(im1,filtdiff','valid');
  fxim2 = conv2(im2,filtdiff','valid');
  fxim = (fxim1+fxim2)/2;
  fyim1 = conv2(im1,filtdiff,'valid');
  fyim2 = conv2(im2,filtdiff,'valid');
  fyim = (fyim1 + fyim2)/2;
  smim1 = conv2(im1,filtsmooth,'valid');
  smim2 = conv2(im2,filtsmooth,'valid');
  diffim = smim1-smim2;
  xfyim = fyim*diag(x);
  yfxim = diag(y)*fxim;


  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  %
  % calculate pseudoinverse
  %
  
  % build matrices for overdetermined system
  df = reshape(diffim, sx*sy, 1);
  fx = reshape(fxim, sx*sy, 1);
  fy = reshape(fyim, sx*sy, 1);
  xfy = reshape(xfyim, sx*sy, 1);
  yfx = reshape(yfxim, sx*sy, 1);
  F = [xfy - yfx, fx, fy];
  c = cond(F);
  if (rank(F) ~= 3) warning('Matrix F is of low rank.'); end;

  % find optimal d = [tx,ty,w]
  Finv = pinv(F);
  d = Finv*df;

  M = [ 1 -d(1) d(2); d(1) 1 d(3); 0 0 1];

return;











