eval('d;','d = 2;');
eval('N;','N = 21;');
clf;
eval('Csi;','Csi=1.25*d;');
eval('Cse;','Cse=d/1.25;');
eval('rng;','rng=d/4;');
eval('steps;','steps=10;');
H=fspecial('log',N,d);
H=H/sum(sum(H.^2));
subplot(2,2,1);
p=mesh(H,'LineWidth',[1.5]); 
set(get(p,'Parent'),'XLim',[0 N]);
set(get(p,'Parent'),'YLim',[0 N]);
set(get(p,'Parent'),'ZLim',[min(min(H)) max(max(H))]);

subplot(2,2,2);
winsize = floor(N/2+1)-1;
p=plot((0-winsize):winsize,H(floor(N/2+1),:));
Ylim=get(get(p,'Parent'),'YLim');
set(p,'LineWidth',[2]);
subplot(2,2,3);
%imshow(H*128/max(max(abs(H)))+128); axis('image'); 
%pause;
%colormap('default');
%We need to find the best si and se
min_err = 100000*N*N;
step = rng/steps;
for si=(Csi-rng):step:(Csi+rng),
   for se=(Cse-rng):step:(Cse+rng),
      Gsi = fspecial('gaussian',N,si);
      Gse = fspecial('gaussian',N,se);
      dog = Gsi-Gse;
      if sum(sum(dog.^2)) > 0.001,
       dog = dog/sum(sum(dog.^2));
       err = sum(sum((dog-H).^2));
       if err < min_err,
         bsi = si;
         bse = se;
         min_err = err;
      end
     end
   
   end
end
si = bsi;
se = bse;
Gsi = fspecial('gaussian',N,si);
Gse = fspecial('gaussian',N,se);
dog = Gsi-Gse;
dog = dog/sum(sum(dog.^2));
disp([sprintf('si=%12.10f se=%12.10f error =%12.10f',si,se,min_err)]);
subplot(2,2,3);
p=mesh(dog,'LineWidth',[1.5]); 
set(get(p,'Parent'),'XLim',[0 N]);
set(get(p,'Parent'),'YLim',[0 N]);
set(get(p,'Parent'),'ZLim',[min(min(H)) max(max(H))]);
subplot(2,2,4);
p=plot((0-winsize):winsize,dog(floor(N/2+1),:));
set(p,'LineWidth',[2]);
hold on;
p=plot((0-winsize):winsize,H(floor(N/2+1),:),'r--');
set(p,'LineWidth',[2]);

set(get(p,'Parent'),'YLim',Ylim);
set(get(get(p,'Parent'),'Parent'),'Color',[1 1 1]);
%pause;
figure(2);
subplot(1,2,1);
p=mesh(dog-H,'LineWidth',[1.5]);
set(get(p,'Parent'),'XLim',[0 N]);
set(get(p,'Parent'),'YLim',[0 N]);
set(get(p,'Parent'),'ZLim',[min(min(H)) max(max(H))]);
subplot(1,2,2);
p=plot((0-winsize):winsize,dog(floor(N/2+1),:)-H(floor(N/2+1),:));
set(get(p,'Parent'),'YLim',Ylim);

set(p,'LineWidth',[2]);

rng = rng/steps*2;
Cse = bse;
Csi = bsi;