function [imest,psfest]=bdecon2(blurry,t,d)

%	[imest,psfest]=bdecon1(blurry,t,d)
%
%		BDECON1 uses the NAS-RIF type algorithm to determine an
%	estimate to the true image from the blurred one. An initial PSF of
%	an impulse is used, which is altered through numerous iterations
%	to optimize a cost function, which is in turn related to known
%	qualities of the true image. Goes through IT iterations, changing
%	u by T*grad(J). d is the rate of change at which the loop stops
%	iterating, ie for d=2, the loop stop when the change of 
%	u/u-u < .02




u0=[0 0 1 0 0]; u0=u0'*u0;
bsize=size(blurry);
maxbsize=max(bsize);
bgcol=150;
t=t;

%initializations and constants ^^
%bgcol- background color, 0 is black, 255 is full white.


%da loop
u=u0;
usize=size(u);

%initialize dd & uk
dd=1;
uk=u0;

while (abs(dd)>=d)

FK=fft2(u,bsize(1),bsize(2)).*fft2(blurry,bsize(1),bsize(2));
fk=abs(ifft2(FK,bsize(1),bsize(2)));

shine=(1+sign(fk))/2;


% Jk loop: this could also be point mult., but way slower

bl=bgcol*ones((bsize+usize-1));

bx1=(usize(1)-1)/2+1;
bx2=bx1+bsize(1)-1;
by1=(usize(2)-1)/2+1;
by2=by1+bsize(2)-1;

bx2-bx1;
by2-by1;

bl(bx1:bx2,by1:by2)=blurry;

x1=(usize(1)-1)/2;
y1=(usize(2)-1)/2;

Jk=zeros(usize);
	for(i=-x1:x1)
		for(j=-y1:y1)

Jk(i+bx1,j+by1)=2*sum(sum((fk.*shine.*bl((bx1+i):(bx2+i),...
		((by1+j):(by2+j))) + ...
	2*((fk-bgcol).*bl((bx1+i):(bx2+i), ...
                ((by1+j):(by2+j)))))));

		end
	end


%new filter
u=u - t*(Jk/(sum(sum(Jk))));

%use to compare for autostop program
ut=u-uk;
dd=sum(sum(ut))/sum(sum(u))

uk=u;

end


%exit stuff:

psfest=u;

IMEST=fft2(u,bsize(1),bsize(2)).*fft2(blurry,bsize(1),bsize(2));
imest=abs(ifft2(IMEST));

why
