function [a,g,e0,e,region]=gaussnewton(i1,i2,a0,g0,region)

% function [a,g,e0,e,region]=gaussnewton(i1,i2,a0,g0,region)

  [auxy,auxx]=find(region);
  centx=mean(auxx);
  centy=mean(auxy);
  a0c=[a0(1:4); a0(5:6)+[a0(1) a0(2); a0(3) a0(4)]*[centx;centy]];
[limy,limx]=size(i1);
%  i2w=aff(aff(i2,[1 0 0 1 -centx -centy]),a0c);
i2w=g0*aff(i2,a0);
[ix,iy]=grad(i2w);

ix=ix.*region;iy=iy.*region;
ix2=ix.^2;iy2=iy.^2;ixiy=ix.*iy;
it=i1-i2w;

m1=sum(sum(kron(ones(1,limy)',(1-centx:limx-centx).^2).*ix2));
m2=sum(sum(kron((1-centy:limy-centy)',(1-centx:limx-centx)).*ix2));
m3=sum(sum(kron(ones(1,limy)',(1-centx:limx-centx).^2).*ixiy));
m4=sum(sum(kron((1-centy:limy-centy)',(1-centx:limx-centx)).*ixiy));
m5=sum(sum(kron(ones(1,limy)',(1-centx:limx-centx)).*ix2));
m6=sum(sum(kron(ones(1,limy)',(1-centx:limx-centx)).*ixiy));
m7=sum(sum(kron((1-centy:limy-centy).^2',ones(1,limx)).*ix2));
m8=sum(sum(kron((1-centy:limy-centy)',ones(1,limx)).*ixiy));
m9=sum(sum(kron((1-centy:limy-centy).^2',ones(1,limx)).*ixiy));
m10=sum(sum(kron((1-centy:limy-centy)',ones(1,limx)).*ix2));
m11=sum(sum(kron(ones(1,limy)',(1-centx:limx-centx).^2).*iy2));
m12=sum(sum(kron((1-centy:limy-centy)',(1-centx:limx-centx)).*iy2));
m13=sum(sum(kron(ones(1,limy)',(1-centx:limx-centx)).*iy2));
m14=sum(sum(kron((1-centy:limy-centy).^2',ones(1,limx)).*iy2));
m15=sum(sum(kron((1-centy:limy-centy)',ones(1,limx)).*iy2));
m16=sum(sum(ix2));
m17=sum(sum(ixiy));
m18=sum(sum(iy2));
m=[m1  m2  m3  m4  m5  m6 sum(sum(kron(ones(1,limy)',(1-centx:limx-centx)).*ix.*i2w));
   m2  m7  m8  m9  m10 m8 sum(sum(kron((1-centy:limy-centy)',ones(1,limx)).*ix.*i2w));
   m3  m4  m11 m12 m6  m13 sum(sum(kron(ones(1,limy)',(1-centx:limx-centx)).*iy.*i2w));
   m4  m9  m12 m14 m8  m15 sum(sum(kron((1-centy:limy-centy)',ones(1,limx)).*iy.*i2w));
   m5  m10 m6  m8  m16 m17 sum(sum(ix.*i2w));
   m6  m8  m13 m15 m17 m18 sum(sum(iy.*i2w));
   sum(sum(kron(ones(1,limy)',(1-centx:limx-centx)).*ix.*i2w))   sum(sum(kron((1-centy:limy-centy)',ones(1,limx)).*ix.*i2w))   sum(sum(kron(ones(1,limy)',(1-centx:limx-centx)).*iy.*i2w))   sum(sum(kron((1-centy:limy-centy)',ones(1,limx)).*iy.*i2w))  sum(sum(ix.*i2w)) sum(sum(iy.*i2w)) sum(sum(i2w.^2))];
b=[sum(sum(kron(ones(1,limy)',(1-centx:limx-centx)).*ix.*it));
   sum(sum(kron((1-centy:limy-centy)',ones(1,limx)).*ix.*it));
   sum(sum(kron(ones(1,limy)',(1-centx:limx-centx)).*iy.*it));
   sum(sum(kron((1-centy:limy-centy)',ones(1,limx)).*iy.*it));
   sum(sum(ix.*it));
   sum(sum(iy.*it));
   sum(sum(i2w.*it))];

aux=m\b;
da=aux(1:6);dg=aux(7);
%dg=sum(sum(i2w.*it))/sum(sum(i2w.^2))

  ac=a0c+da;
  a=[ac(1:4);ac(5:6)-[ac(1) ac(2); ac(3) ac(4)]*[centx;centy]];
g=g0+dg
e0=sqrt(sum(sum(((i1-i2w).^2).*region))/sum(sum(region)));
e=sqrt(sum(sum(((i1-g*aff(i2,a)).^2).*region))/sum(sum(region)));

region=ones(size(i1));
im=i2;
R=[a(1) a(2); a(3) a(4)];
[m,n]=size(im);
[y,x]=meshgrid(1:m,1:n);
xy=[x(:) y(:)]*R'; % Rotate points
xy(:,1)=xy(:,1)+a(5);
xy(:,2)=xy(:,2)+a(6);
xx=xy(:,1); yy=xy(:,2);
%   xx=min(max(xx,1),n);
%   yy=min(max(yy,1),m);
aux = interp2(im,xx,yy,'*linear');
imaff=reshape(aux,[n,m])';
region(find(isnan(imaff)))=zeros(size(find(isnan(imaff))));