%kvgef 'First Derivative Operator for Symmetric Exponential Filter (K1)'
% This MatLab function was automatically generated by a converter (KhorosToMatLab) from the Khoros vgef.pane file
%
% Parameters: 
% InputFile: i 'Input Image', required: 'input image'
% Double: a0 ' a0 ', default: 0.45: 'Exponential Filter Parameter for GEF'
% Integer: t1 ' T1 ', default: 3: '1st Hysteresis Threshold'
% Integer: t2 ' T2 ', default: 4: '2nd Hysteresis Threshold'
% Integer: l 'Length', default: 5: 'Minimum Number of Pixels in a Segment'
% OutputFile: o 'Output Image', required: 'output image'
%
% Example: o = kvgef(i, {'i','';'a0',0.45;'t1',3;'t2',4;'l',5;'o',''})
%
% Khoros helpfile follows below:
%
%  PROGRAM
% vgef - First Derivative Operator for Symmetric Exponential Filter  (K1)
%
%  DESCRIPTION
% 
% This Algorithm has been written by the team of:
% 
% 
% Professor SERGE CASTAN
% 
% 
% IRIT, URA 1399
% 
% 
% 118, route de Narbonne  31062 Toulouse  FRANCE
% 
% OPTIMAL FILTER(ISEF) FOR EDGE DETECTION
% 
% PRINCIPLE :
% I. Introduction 
% Edge detection is one of the most important subjects in image processing, 
% which finds wide applications in the pattern recognition, the scene analysis and 
% the 3-D vision, because the edges correspond in general to the important changes 
% of physical or geometrical properties of objects in the scene and they are widely 
% used as primitives in the pattern recognition, the image matching etc.
% 
% The edges coincide, generally speaking, grey level transition, they can 
% be detected by maxima of gradient or the zero-crossing of the second 
% derivatives calculated by some differential operators. Because the 
% differential operators are sensitive to noise, a preprocessing such as 
% smoothing is often necessary to eliminate the noise. A well-known smoothing 
% filter is Gaussian filter and the edges can therefore be detected by a 
% Laplacian-Gaussian filter. But there is an essential difficulty of 
% the Laplacian-Gaussian filter which is the contradiction between the 
% smoothing effect and the precision of edge localization. To overcome this 
% difficulty,  We proposed the optimal linear filter based on one step model 
% (a step edge and the white noise) and the multi-edge model [9][10][11]. 
% This optimal smoothing filter is a symmetric exponential filter of an 
% infinitely large window size and can be realized by very simple recursive 
% algorithm. It  is proved that the band limited Laplacian of an input image 
% filtered by this filter can be calculated from the Difference between the input 
% and the output of this Recursive Filter (DRF). The edges detected by DRF 
% method are less noisy and with a much better precision of localization.
% 
% The maxima of gradient or zeros of the second directional derivative 
% along the gradient are a natural definition of intensity edges. Zeros of the 
% Laplacian are only extensively used for their computational convenience. 
% However, we must stress here that the zeros crossing of the Laplacian are 
% not always coincided with the maxima of gradient, for example, the zeros 
% of the Laplacian are farther apart than the maxima of gradient for circularly 
% symmetric patterns, this lack of localization by the Laplacian can also be 
% seen in the fact that zeros of Laplacian "swing wide" of corners. Therefore, 
% it had better to detect the edges by maxima of gradient or zeros of the second 
% directional derivative along the gradient .
% 
% So we propose two methods for edges detection, one uses maxima of Gradient 
% (GEF), another uses the zeros crossing of Second directional Derivative along 
% the gradient (SDEF).
% 
% II. The First and the Second Directional Derivative Operators for Symmetric 
% Exponential Filter
% A normalized symmetric exponential filter on 1-D can be written :
% fL(x)=C*a0*(1-a0)**|x|=f1(x)#f2(x)=C*(f1(x)+f2(x)-a0*d(x)) (1)
% 
% where: 
% C=1/(2-a0), # means the convolution, d(x) is dirac function.
% f1(x)={0 if x<0  and a0*(1-a0)**x   if x>= 0  (2)
% f2(x)={0 if x>0  and a0*(1-a0)**-x  if x<= 0  (2)
% we can write the first derivative operator of exponential filter:
% 
% fL'(x)=f2(x)-f1(x)         (3)
% And we can obtain the normalized second derivative operator of 
% exponential filter :
% fL"(x)=f1(x)+f2(x)-2*d(x)  (4)
% 
% Because the exponential function is separable, we can write out 2-D 
% exponential filter:
% f(x,y) = fL(x)*fL(y)                                            (5)
% 
% From the equations (18),(19) and (20),  the first and the second 
% directional derivative operators for symmetric exponential filter can be 
% written like this:
% 
% 
% fx(x,y) = fL(y)*(f2(x)-f1(x))          (6)
% 
% 
% fy(x,y) = fL(x)*(f2(y)-f1(y))          (7)
% 
% 
% fxx(x,y) = fL(y)*(f1(x)+f2(x)-2*d(x))  (8)
% 
% 
% fyy(x,y) = fL(x)*(f1(y)+f2(y)-2*d(y))  (9) 
% 
% III. The Recursive Algorithm for realizing the These Directional Derivative 
% Operators of Symmetric Exponential Filter
% 
% The exponential filter is an IIR filter corresponding to an infinite 
% window size, so we should realize the functions f1(x) and f2(x) (see formula 
% (2)) by a recursive algorithm.
% 
% Supposing I(x,y) is the input image, 
% I1(x,y)=I(x,y)#f1(x) and I2(x,y)=I(x,y)#f2(x), 
% we have the recursive algorithm :
% I1(x,y)=I1(x-1,y)+a0*(I(x,y)-I1(x-1,y))    (10)
% I2(x,y)=I2(x+1,y)+a0*(I(x,y)-I2(x+1,y))
% 
% From the equations (1),(6),(7),(8) and (9), the band-limited first 
% and second directional derivative of input image can be calculated by the 
% recursive algorithm f1 and f2 as follows 
% Ix(x,y)=I(x,y)#f1(y)#f2(y)#(f2(x)-f1(x))   (11)
% Ixx(x,y)=I(x,y)#f1(y)#f2(y)#(f2(x)+f1(x))-2*I(x,y)#f1(y)#f2(y)  (12)
% Iy(x,y) = I(x,y)#f1(x)#f2(x)#(f2(y)-f1(y))                        (13)
% Iyy(x,y) = I(x,y)#f1(x)#f2(x)#(f2(y)+f1(y))-2*I(x,y)#f1(x)#f2(x)  (14)
% 
% With this algorithm, we can calculate at the same time the band-limited 
% first and second directional derivative Ix and Ixx ( or Iy and Iyy ) of 
% input image.
% 
% IV. Edges Detection
% The band-limited first and second directional derivative of input image 
% can be obtained by the algorithms as stated above. Using them, we can then 
% realize the edge detection for an image. 
% 
% The maxima of gradient or zeros of the second directional derivative 
% along the gradient are a natural way of characterizing and localizing intensity 
% edges,  so we present here to detect the edges from the maxima of gradient or 
% zeros of the second directional derivative along the gradient by using the 
% differential operators of exponential filter.
% 
% 
% 1. Edges from the maxima of gradient
% Using the first directional derivative operator of the exponential filter, 
% the two band-limited first directional derivatives Ix and Iy can be calculated, 
% and the gradient vector can be therefore determined approximatively for every 
% point in image. The gradient magnitude image is then non maxima suppressed in 
% the gradient direction and thresholded with hysteresis, i.e. if the entire 
% segment of the contour lies above a low threshold T1, and at least one of part of 
% which is above a high threshold T2, that contour is output. The non maxima 
% suppression scheme requires three points, one of which will be the current 
% point, and the other two should be estimated of the gradient magnitude at 
% points displaced from the current point by vector normal to the edge direction.
% 
% 
% 2. Edges from the zero-crossings of the second directional derivative along the 
% gradient direction
% Because edges detected from local gradient maxima can not be a pixel 
% width (less good precision of localization), we propose an other method which 
% detect the edges from the zeros crossing of the second directional derivative 
% along the gradient direction.
% 
% We can calculate Ix,Iy,Ixx and Iyy by using the method shown in 
% paragraph III, and therefore obtain approximatively the gradient vector 
% and the second derivative in the gradient direction for every point in image. 
% We extract at first the zero crossing of second derivative along the gradient 
% direction on which the gradient magnitude must be above a low threshold, so an 
% edges image is obtained. To this image, the entire segment of the contour will 
% be kept, if the gradient magnitude on at least one part of this contour is above
% a high threshold.
% V. Comparison of Performance of the Filters
% Filtering is a problem of estimation from noisy signal, and edge 
% detection is a problem of estimating the position of maximal local signal 
% change. Up to now, many works are done for edges detection in image, and 
% different filters are proposed, for example, Gaussian filter, Canny filter, 
% exponential filter, Deriche filter etc...
% 
% We appreciate the performance of the filters as follows :
% 
% 
% (1) Precision of edge localization
% According to our analysis [10], we can calculate the average 
% localization error xe for Gaussian filter, Canny filter [15], Deriche 
% filter [16] and the exponential filter : 
% 
% 
% xeG = (4*(2*e*3.14)**0.5)/a
% 
% 
% xec = 0.81/a
% 
% 
% xeD = 4*exp(-1)/a = 1.47/ a
% 
% 
% xeE = 0
% 
% 
% i.e. xeG > xeD > xec > xeE = 0.
% 
% So,  we can see that the exponential filter localizes edge points with 
% the best precision.
% 
% 
% (2) Signal/Noise ratio on the edge point detected
% Because xe is the average estimation for the position of the edge point 
% detected, we propose to calculate Signal/Noise ratio (Eq.(7)) at the point xe.
% 
% And the signal/noise ratio for the Gaussian filter, Canny filter, 
% Deriche filter and the exponential filter is : 
% 
% 
% SNRG = 2*s*exp(-32*s)/ ((3.14)**0.5)
% 
% 
% SNRc = 0.39/ a
% 
% 
% SNRD = 0.64/ a
% 
% 
% SNRE =1/ a
% 
% 
% i.e. SNRE > SNRD > SNRc > SNRG.
% 
% Then, we see that the exponential filter has the best noise eliminating 
% effect among the above four.
% 
% 
% (3) Complexity of calculation
% For the complexity of calculation, we only tell the difference from 
% exponential filter and Deriche filter [16], because they are implemented by 
% recursive algorithms which have a simpler calculation.
% 
% Because ISEF can be realized by first order recursive filter, the 
% ISEF algorithms are much simpler than that of Deriche filter. Besides, the 
% ISEF algorithms can be implemented independently to every line and every column,
% it can be easily realized by a parallel system.
% 
% According to the analysis results above, the ISEF filter is superior 
% to the others at the 3 principal aspects of the performance of the filter.
% 
% VI. Conclusion
% The symmetric exponential filter of an infinite large window size is 
% an optimal linear filter deduced from one step edge model and the multi-edge 
% model, now we further prove that the symmetric exponential filter is the optimal
% edge detection filter in the criteria of the signal to noise ratio, the 
% localization precision and unique maximum. Obviously, the real images will be 
% still more complicated than these models, however DRF method has already provided 
% good results for different type of images. The results obtained through the 
% two new methods further show the superior performance of this filter. 
% The theoretical analysis for the performance of the filters shows also that 
% the exponential filter is superior to the other current filters.
% 
% The first and second directional derivative operators can be realized by 
% recursive algorithm and calculated at the same time. The new algorithms 
% are therefore very simple as well as DRF algorithm, and they are also easy 
% to implemented in a parallel way.
%
%  
%
%  EXAMPLES
% vgef -i cross.xv -o output.xv
% 
% will compute a edge extraction by using the first directional
% derivative operator for symmetric exponential filter on
% image cross.xv and write the result in output.xv
% 
% For the GEF edge extraction the filter parameter value is generally
% a0 between 0.40 and 0.9 (0.45  [default] is a good value)
% 
% The hysteresis Thresholds T1 and T2 can be set to 12 and 15.
% 
% The minimum Length of segment is up to the user. It is set by default to 10.
% 
%
%  "SEE ALSO"
%
%  RESTRICTIONS 
% works only with single band byte images.
% 
% Note that this routine was converted directly from Khoros 1.0.  Due to
% the mixing of unsigned versus signed char, and the difference in the rules
% that resolve the mixing of these data types between ANSI C and the original
% K&R C that the code was written in, the results of running this program
% on an image will not be the same as with vdrf in Khoros 1.0.  There is
% really no way to resolve which version is more correct, but the ANSI C version
% should at least be more consistent across machine architectures.
%
%  REFERENCES 
% 
% [1] W.K. PRATT, Digital Image Processing, New York, 1978.
% 
% [2] J. PREWITT, Object Enhancement and Extraction,Picture Processing and
% Psychopictories, Etd. by B. Lipkin and A. Rosenfeld, New York, pp.75-149, 1970.
% 
% [3] M. HUECKEL, An Operator Which Locates Edges in Digitized Pictures.
% J.A.C.M., Vol. 18, pp 113-125, 1971.
% 
% [4] R.O. DUDA and P.E. HART, Pattern Classification and Scene Analysis.
% Wiley, New York, 1973.
% 
% [5] R. HARALICK,Edge and Region Analysis for Digital Image Data. C.G.I.P.,
% Vol. 12, pp 60-73, 1980.
% 
% [6] R.HARALICK and L.WATSON, A Facet Model for Image Data. C.G.I.P., Vol. 15,
% pp 113-129, 1981.
% 
% [8] D. MARR and E.C. HILDRETH, Theory of Edge Detection. Proc. R. Soc.
% Lond. B, Vol. 207, pp 187-217, 1980.
% 
% [9] J. SHEN and S. CASTAN, Un nouvel algorithme de detection de contours,
% proceedings of 5th Conf. on P.R.&.A.I. (in French), Grenoble, 1985.
% 
% [10] J. SHEN and S. CASTAN, An Optimal Linear Operator for Edge Detection.
% Proc. CVPR'86,  Miami,1986.
% 
% [11] J. SHEN and S. CASTAN, Edge Detection Based on Multi-Edge Models.Proc.
% SPIE'87, Cannes, 1987.
% 
% [12] J. SHEN and S. CASTAN, Further Results on DRF Method for Edge Detection.
% 9th I.C.P.R., ROME, 1988.
% 
% [13] V.TORRE and T.A.POGGIO, On Edge Detection IEEE Transaction on Pattern
% Analysis and    Machine Intelligence, Vol. Pami-8, N 2, March 1986.
% 
% [14] J.S.CHEN and G.MEDIONI, Detection, Localization, and Estimation of Edges.
% IEEE  Transaction on Pattern  Analysis and Machine Intelligence, Vol. 11, N 2,
% February  1989.
% 
% [15] J.F.CANNY, Finding Edges And Lines in Images.
% MIT Technical Report N 720, 1983.
% 
% [16] R. DERICHE, Optimal Edge Detection Using Recursive Filtering.
% In proc. First International Conference on Computer Vision,
% London, June 8-12 1987.
%
%  COPYRIGHT
% Copyright (C) 1993 - 1997, Khoral Research, Inc. ("KRI")  All rights reserved.
% 


function varargout = kvgef(varargin)
if nargin ==0
  Inputs={};arglist={'',''};
elseif nargin ==1
  Inputs=varargin{1};arglist={'',''};
elseif nargin ==2
  Inputs=varargin{1}; arglist=varargin{2};
else error('Usage: [out1,..] = kvgef(Inputs,arglist).');
end
if size(arglist,2)~=2
  error('arglist must be of form {''ParameterTag1'',value1;''ParameterTag2'',value2}')
 end
narglist={'i', '__input';'a0', 0.45;'t1', 3;'t2', 4;'l', 5;'o', '__output'};
maxval={0,1,255,255,100000,0};
minval={0,0,0,0,0,0};
istoggle=[0,1,1,1,1,0];
was_set=istoggle * 0;
paramtype={'InputFile','Double','Integer','Integer','Integer','OutputFile'};
% identify the input arrays and assign them to the arguments as stated by the user
if ~iscell(Inputs)
Inputs = {Inputs};
end
NumReqOutputs=1; nextinput=1; nextoutput=1;
  for ii=1:size(arglist,1)
  wasmatched=0;
  for jj=1:size(narglist,1)
   if strcmp(arglist{ii,1},narglist{jj,1})  % a given argument was matched to the possible arguments
     wasmatched = 1;
     was_set(jj) = 1;
     if strcmp(narglist{jj,2}, '__input')
      if (nextinput > length(Inputs)) 
        error(['Input ' narglist{jj,1} ' has no corresponding input!']); 
      end
      narglist{jj,2} = 'OK_in';
      nextinput = nextinput + 1;
     elseif strcmp(narglist{jj,2}, '__output')
      if (nextoutput > nargout) 
        error(['Output nr. ' narglist{jj,1} ' is not present in the assignment list of outputs !']); 
      end
      if (isempty(arglist{ii,2}))
        narglist{jj,2} = 'OK_out';
      else
        narglist{jj,2} = arglist{ii,2};
      end

      nextoutput = nextoutput + 1;
      if (minval{jj} == 0)  
         NumReqOutputs = NumReqOutputs - 1;
      end
     elseif isstr(arglist{ii,2})
      narglist{jj,2} = arglist{ii,2};
     else
        if strcmp(paramtype{jj}, 'Integer') & (round(arglist{ii,2}) ~= arglist{ii,2})
            error(['Argument ' arglist{ii,1} ' is of integer type but non-integer number ' arglist{ii,2} ' was supplied']);
        end
        if (minval{jj} ~= 0 | maxval{jj} ~= 0)
          if (minval{jj} == 1 & maxval{jj} == 1 & arglist{ii,2} < 0)
            error(['Argument ' arglist{ii,1} ' must be bigger or equal to zero!']);
          elseif (minval{jj} == -1 & maxval{jj} == -1 & arglist{ii,2} > 0)
            error(['Argument ' arglist{ii,1} ' must be smaller or equal to zero!']);
          elseif (minval{jj} == 2 & maxval{jj} == 2 & arglist{ii,2} <= 0)
            error(['Argument ' arglist{ii,1} ' must be bigger than zero!']);
          elseif (minval{jj} == -2 & maxval{jj} == -2 & arglist{ii,2} >= 0)
            error(['Argument ' arglist{ii,1} ' must be smaller than zero!']);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} < minval{jj})
            error(['Argument ' arglist{ii,1} ' must be bigger than ' num2str(minval{jj})]);
          elseif (minval{jj} ~= maxval{jj} & arglist{ii,2} > maxval{jj})
            error(['Argument ' arglist{ii,1} ' must be smaller than ' num2str(maxval{jj})]);
          end
        end
     end
     if ~strcmp(narglist{jj,2},'OK_out') &  ~strcmp(narglist{jj,2},'OK_in') 
       narglist{jj,2} = arglist{ii,2};
     end
   end
   end
   if (wasmatched == 0 & ~strcmp(arglist{ii,1},''))
        error(['Argument ' arglist{ii,1} ' is not a valid argument for this function']);
   end
end
% match the remaining inputs/outputs to the unused arguments and test for missing required inputs
 for jj=1:size(narglist,1)
     if  strcmp(paramtype{jj}, 'Toggle')
        if (narglist{jj,2} ==0)
          narglist{jj,1} = ''; 
        end;
        narglist{jj,2} = ''; 
     end;
     if  ~strcmp(narglist{jj,2},'__input') && ~strcmp(narglist{jj,2},'__output') && istoggle(jj) && ~ was_set(jj)
          narglist{jj,1} = ''; 
          narglist{jj,2} = ''; 
     end;
     if strcmp(narglist{jj,2}, '__input')
      if (minval{jj} == 0)  % meaning this input is required
        if (nextinput > size(Inputs)) 
           error(['Required input ' narglist{jj,1} ' has no corresponding input in the list!']); 
        else
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        end
      else  % this is an optional input
        if (nextinput <= length(Inputs)) 
          narglist{jj,2} = 'OK_in';
          nextinput = nextinput + 1;
        else 
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end;
     else 
     if strcmp(narglist{jj,2}, '__output')
      if (minval{jj} == 0) % this is a required output
        if (nextoutput > nargout & nargout > 1) 
           error(['Required output ' narglist{jj,1} ' is not stated in the assignment list!']); 
        else
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
          NumReqOutputs = NumReqOutputs-1;
        end
      else % this is an optional output
        if (nargout - nextoutput >= NumReqOutputs) 
          narglist{jj,2} = 'OK_out';
          nextoutput = nextoutput + 1;
        else 
          narglist{jj,1} = '';
          narglist{jj,2} = '';
        end;
      end
     end
  end
end
if nargout
   varargout = cell(1,nargout);
else
  varargout = cell(1,1);
end
global KhorosRoot
if exist('KhorosRoot') && ~isempty(KhorosRoot)
w=['"' KhorosRoot];
else
if ispc
  w='"C:\Program Files\dip\khorosBin\';
else
[s,w] = system('which cantata');
w=['"' w(1:end-8)];
end
end
[varargout{:}]=callKhoros([w 'vgef"  '],Inputs,narglist);
