hold off
%imshow(Image,256)
sizeR = size
sizeC = size

% OBTAIN FILTER PARAMETERS
% Rebuild the filters only if needed

f=input('Do you need a new filter? 0 = no, else = yes   ')
if f~=0
	alpha=input('Enter alpha in -1.5 to 1.5, the OZCO shape parameter.   ')
	beta=input('Enter beta in 0.2 to 2.0, the OZCO scale parameter.   ')


% 	BUILD the FILTERS
% 	px (py) is the smoothing function (second integral of OZCO).
% 	dx (dy) is the smoothed differentiator (first integral of OZCO).

% 	First, some constants we'll need...

	tau = 1;
	alpha0 = (alpha + 2)/beta^2;
	alpha1 = tau*(alpha+2)/beta + 0.5*tau^2*(alpha+1) - 2*(alpha+2)/beta^2;
	alpha1 = alpha1*exp(-beta*tau);
	alpha2 = -tau*(alpha+2)/beta + 0.5*tau^2*(alpha+1) + (alpha+2)/beta^2;
	alpha2 = alpha2*exp(-2*beta*tau);
	beta1 = -3*exp(-beta*tau);
	beta2 = 3*exp(-2*beta*tau);
	beta3 = -exp(-3*beta*tau);
	gamma1 = (0.5*beta*tau^2*(alpha+1) + tau)*exp(-beta*tau);
	gamma2 = (0.5*beta*tau^2*(alpha+1) - tau)*exp(-2*beta*tau);

% Build the smoothing function feedforward coefficient vector

	pb = [alpha0 alpha1 alpha2];

% Build the smoothing function feedback coefficient vector

	pa = [1 beta1 beta2 beta3];

% Build the differentiator feedforward coefficient vector

	db = [0 gamma1 gamma2];

% The differentiator feedback coefficients are the same as the smoother's

	da = pa;

% Run the filters over a unit pulse and display the response

	hold off
	test = [zeros(1,100) 1 zeros(1,100)];
	pxr = filter(pb,pa,test);
	dxr = filter(db,da,test);
	pxl(201:-1:1) = filter(pb,pa,test(201:-1:1));
	dxl(201:-1:1) = filter(db,da,test(201:-1:1));
	px = pxl + pxr - [zeros(1,100) pxr(101) zeros(1,100)];
	dx = dxr - dxl;
	plot(px,'g')
	hold on
	plot(10*dx,'b')
	title(['OZCO Infinite Impulse Responses, alpha = ',num2str(alpha),', beta = ',num2str(beta)])
	hold off
end

% OBTAIN ZC DETECTION and VALIDATION PARAMETERS

z=input('Do you need new zc parameters? 0=no, other number = yes   ')
	if z~=0
	Imult=input('Enter an image brightness multiplier for ZC in ROI   ')
	ZCThr=input('Enter a threshold for ZC validation   ')
	ZCThr = -ZCThr;
	end
%	ImThr=Imult*max(max(Image));

% APPLY FILTERS
% Only do the convolutions if either the image or the filter has changed

%if n~=0 | f~=0
if f~=0

	fprintf('Convolutions in progress\n')
% 	Make space for intermediate images

	Ximagef = zeros(sizeR,sizeC);
	Yimagef = zeros(sizeR,sizeC);
	Ximageb = zeros(sizeR,sizeC);
	Yimageb = zeros(sizeR,sizeC);
	Ximage = zeros(sizeR,sizeC);
	Yimage = zeros(sizeR,sizeC);
	GradSq = zeros(sizeR,sizeC);
	GradX = zeros(sizeR,sizeC);
	GradY = zeros(sizeR,sizeC);
	DI = zeros(sizeR,sizeC);

% 	CONVOLVE with Differentiator in X, Smoothing Function in Y
% 	This is the X component of the smoothed gradient.
% 	CONVOLVE with Differentiator in Y, Smoothing Function in X
% 	This is the Y component of the smoothed gradient.

	for r=1:sizeR
		Ximagef(r,:)=filter(db,da,Image(r,:));
		Yimagef(r,:)=filter(pb,pa,Image(r,:));
		Ximageb(r,sizeC:-1:1)=filter(db,da,Image(r,sizeC:-1:1));
		Yimageb(r,sizeC:-1:1)=filter(pb,pa,Image(r,sizeC:-1:1));
	end

	Ximage = Ximagef - Ximageb;
	Yimage = Yimagef + Yimageb - alpha0*Image;

	for c=1:sizeC
		Ximagef(:,c)=filter(pb,pa,Ximage(:,c));
		Yimagef(:,c)=filter(db,da,Yimage(:,c));
		Ximageb(sizeR:-1:1,c)=filter(pb,pa,Ximage(sizeR:-1:1,c));
		Yimageb(sizeR:-1:1,c)=filter(db,da,Yimage(sizeR:-1:1,c));
	end

	Ximage = Ximagef + Ximageb - alpha0*Ximage;
	Yimage = Yimagef - Yimageb;

	clear Ximagef;
	clear Ximageb;
	clear Yimagef;
	clear Yimageb;
end

% 	INNER PRODUCT of GRADIENT with ITSELF
% 	This is the gradient squared magnitude.
fprintf('Computing smoothed gradient magnitude squared\n')

	GradSq=(Ximage .* Ximage) + (Yimage .* Yimage);

% 	COMPUTE the GRADIENT of the GRADIENT SQUARED MAGNITUDE
% 	These are the first directional derivatives.
fprintf('Computing first directional derivatives\n')

	[GradX,GradY]=gradient(GradSq);

% 	COMPUTE the INNER PRODUCT of the GRADIENT of GRADIENT SQUARED
% 	with the ORIGINAL SMOOTHED GRADIENT
% 	This forms the second directional derivative in the direction
% 	of the smoothed gradient.  Zero crossings in this image indicate
% 	potential edges, subject to qualification and strength tests.
fprintf('Computing second directional derivatives\n')

	DI=(Ximage .* GradX) + (Yimage .* GradY);
end

% DETECT and QUALIFY ZCs

% ZCs equal border pixels in DI; corresponding Image pixel > ImThr

ZCmap = zeros(sizeR,sizeC);
ZC2 = zeros(sizeR,sizeC);
DxI= zeros(sizeR,sizeC);
GradX = zeros(sizeR,sizeC);
GradY = zeros(sizeR,sizeC);

fprintf('Computing initial ZC map\n')

r = 2:sizeR-1;
c = 2:sizeC-1;
ImThr=Imult*max(max(GradSq));

%ZCmap(r,c) = Image(r,c)>ImThr & DI(r,c) >= 0 & (DI(r-1,c)<0 | DI(r+1,c)<0 |DI(r,c+1)<0 | DI(r,c-1)<0);
ZCmap(r,c) = GradSq(r,c)>ImThr & DI(r,c) >= 0 & (DI(r-1,c)<0 | DI(r+1,c)<0 |DI(r,c+1)<0 | DI(r,c-1)<0);

% imshow(1-ZCmap)
% title(['Bright ZCs, Image # ',num2str(im),', alpha = ',num2str(alpha),',beta = ',num2str(beta),', Imult = ',num2str(Imult)])

% Smoothed gradient consistent with third directional derivative

fprintf('Computing normalized third directional derivatives\n')

[GradX,GradY]=gradient(DI);
DxI=(Ximage .* GradX) + (Yimage .* GradY);
DxI = DxI/max(max(abs(DxI)));

fprintf('Pruning weak ZCs from map\n')

ZCmap = ZCmap & (DxI < ZCThr);
% n=input('Enter any number to see thresholded ZCmap   ')
% imshow(1-ZCmap)
% title(['Strong ZCs, Image # ',num2str(im),', alpha = ',num2str(alpha),',beta = ',num2str(beta),', Imult = ',num2str(Imult),', ZCThr =',num2str(ZCThr)])

% fprintf('Pruning ZCs inconsistent with gradient orientation\n')

% DxI = (DxI .* GradSq);
% ZCmap = ZCmap & (DxI < 0);

% n=input('Enter any number to see qualified ZCmap   ')
% imshow(1-ZCmap)
% title(['Phantoms Deleted, Image # ',num2str(im),', alpha =',num2str(alpha),', beta = ',num2str(beta),', Imult = ',num2str(Imult),',ZCThr = ',num2str(ZCThr)])

% n=input('Enter any number to see qualified ZCmap overlay   ')
% imshow(max(ZCmap,Image))
% title(['Strong ZCs, Image # ',num2str(im),', alpha = ',num2str(alpha),',beta = ',num2str(beta),', Imult = ',num2str(Imult),', ZCThr =',num2str(ZCThr)])