function [cHStanArray,cMaxFilt, cMinFilt,cMinSupport,cMaxSupport] = HStanMultiOrientation_Test(img, options)
%function [cHStanArray,cMaxFilt, cMinFilt] = HStanMultiOrientation(img, options)
%
%img is a grayscale image (one layer)
%cHStanArray is a cell array of hough-line-like features at increasing scales.

%Default Options
D.Verbose = 0;
D.ShowFigures = 0;
D.MaxIterations = 5;
D.ConfidenceCellSupportSize = 7;
D.MaxCellSupportSize = 7;
D.MinCellSupportSize = 7;
D.NOrientations = 18;
D.FilterResolution = 5;
D.FilterType = 'bar';
D.MaxThickness = 3;
D.DecimationSpan = 2;
D.UseLinearVersionOfContinuity = 0;
D.bBilinearDecimation = 0;
D.bIncludeThreshold = 0;
D.bIncludeOrientationSmoothing = 0;
D.bIncludeEdgeConfidence = 1;
D.bConfidenceMultiply = 1;
D.bConfidencePoolOneDir = 0;
if(nargin < 2), options = [];, end;
options = ResolveMissingOptions(options,D);
cfilters = BuildDirectionalDerivatives(options);

cConfidenceSupport = BuildConfidenceSupport(options);
cMinSupport = BuildMinimumSupport(options);
cMaxSupport = BuildMaximumSupport(options);
cHStanArray{1} = ApplyDirectionalFilters(img,cfilters);
if(options.ShowFigures)
  if(options.ShowFigures < 0);
    h = figure(-options.ShowFigures);
  else
    h = figure;
  end
end
for i = 1:options.MaxIterations
  if(options.bIncludeEdgeConfidence)
    for iOri = 1:size(cHStanArray{i},3)
      cConfidenceArray{i}(:,:,iOri) = rankfilter_arbitrarySupport(cHStanArray{i}(:,:,iOri),cConfidenceSupport{iOri});
    end
    % cConfidenceArray{i} = cConfidenceArray{i}.^2;
  end
  if(options.bConfidencePoolOneDir)
    for iOri = 1:size(cHStanArray{i},3)
      cConfidenceArray{i}(:,:,iOri) = 2 * sumfilter_arbitrarySupport(cHStanArray{i}(:,:,iOri),cMinSupport{iOri},@mean);
      cConfidenceArray{i}(:,:,iOri) = min(cConfidenceArray{i}(:,:,iOri),1);
    end
  end
  if(options.bConfidenceMultiply)
    cHStanArray{i} = cHStanArray{i} .* cConfidenceArray{i};
  end
  if(options.bIncludeThreshold)
    cHStanArray{i} = ThresholdArray(cHStanArray{i},options);
  end
  if(options.bIncludeOrientationSmoothing)
    cHStanArray{i} = SmoothOverOrientations(cHStanArray{i},options);
  end
  if(options.Verbose), fprintf('Iteration %d ',i), end;
%  cMinFilt{i} = -ApplyMaximumFilter(-cHStanArray{i}, cMinSupport);
%  cMaxFilt{i} = ApplyMaximumFilter(cMinFilt{i}, cMaxSupport);
  if(options.UseLinearVersionOfContinuity)
    for j = 1:options.NOrientations
      f = -cMaxSupport{j} + 2 * cMinSupport{j};
      cMaxFilt{i}(:,:,j) = imfilter(-cHStanArray{i}(:,:,j), f,'same');
    end
  else
    for j = 1:options.NOrientations
      if(options.Verbose), fprintf('o'), end;
      cMinFilt{i}(:,:,j) = -maxfilter_arbitrarySupport(-cHStanArray{i}(:,:,j), cMinSupport{j});
    end
    for j = 1:options.NOrientations
      if(options.Verbose), fprintf('M'), end;
      cMaxFilt{i}(:,:,j) = maxfilter_arbitrarySupport(cMinFilt{i}(:,:,j), cMaxSupport{j});
    end
  end
  if(options.Verbose), fprintf('D\n'), end;
  cHStanArray{i+1} = ApplyDecimation(cMaxFilt{i},options);
  if(options.ShowFigures)
    subplot(ceil(sqrt(options.MaxIterations)),ceil(sqrt(options.MaxIterations)),i)
    imshow(2 * FourLayerColor(cHStanArray{i},[0]));
  end
end

function o = ApplyDecimation(image,options);
sz = size(image);
sz = sz(1:2) / options.DecimationSpan;
if(options.bBilinearDecimation)
  o = imresize(image,sz,'bilinear');
else
  o = imresize(image,sz,'nearest');
end
function x = NormalizeInput(x);
x = x - min(x(:)) + 10^-6;
x = x/max(x(:));

function cConfidenceSupport = BuildConfidenceSupport(options)
%the support of the minimum are those pixels wich are within .5 of the appropriate line (L1)
for i = 1:options.NOrientations
   cConfidenceSupport{i} = zeros(options.ConfidenceCellSupportSize);
   o(i) = (pi * (i-1)) / options.NOrientations - 0 * pi/2;
   u = [sin(-o(i)),cos(-o(i))];
   u = u / norm(u);
   o(i) = mod(o(i),pi);
   if not((o(i) < (pi/4)) || (o(i) > (3*pi/4)))  % for all x find one best y
     for x = 1:options.ConfidenceCellSupportSize 
        x2 = x - options.ConfidenceCellSupportSize / 2 - .5;  
        for y = 1:options.ConfidenceCellSupportSize 
           y2 = y - options.ConfidenceCellSupportSize / 2 - .5;  
   	   r = [y2, x2];
           d(y) = abs(r * u');
	end
	[m_dy,mi_dy] = min(d);
        cConfidenceSupport{i}(mi_dy,x) = 1;
     end   
   else % for all y, find one best x
     for y = 1:options.ConfidenceCellSupportSize 
        y2 = y - options.ConfidenceCellSupportSize / 2 - .5;  
        for x = 1:options.ConfidenceCellSupportSize 
           x2 = x - options.ConfidenceCellSupportSize / 2 - .5;  
   	   r = [y2, x2];
           d(x) = abs(r * u');
	end
	[m_dx,mi_dx] = min(d);
        cConfidenceSupport{i}(y,mi_dx) = 1;
     end   
   end
end 
  
function cMinSupport = BuildMinimumSupport(options);
%the support of the minimum are those pixels wich are within .5 of the appropriate line (L1)
for i = 1:options.NOrientations
   cMinSupport{i} = zeros(options.MinCellSupportSize);
   o(i) = (pi * (i-1)) / options.NOrientations - 1 * pi/2;
   u = [sin(-o(i)),cos(-o(i))];
   u = u / norm(u);
   o(i) = mod(o(i),pi);
   if not((o(i) < (pi/4)) || (o(i) > (3*pi/4)))  % for all x find one best y
     for x = 1:options.MinCellSupportSize 
        x2 = x - options.MinCellSupportSize / 2 - .5;  
        for y = 1:options.MinCellSupportSize 
           y2 = y - options.MinCellSupportSize / 2 - .5;  
   	   r = [y2, x2];
           d(y) = abs(r * u');
	end
	[m_dy,mi_dy] = min(d);
        cMinSupport{i}(mi_dy,x) = 1;
     end   
   else % for all y, find one best x
     for y = 1:options.MinCellSupportSize 
        y2 = y - options.MinCellSupportSize / 2 - .5;  
        for x = 1:options.MinCellSupportSize 
           x2 = x - options.MinCellSupportSize / 2 - .5;  
   	   r = [y2, x2];
           d(x) = abs(r * u');
	end
	[m_dx,mi_dx] = min(d);
        cMinSupport{i}(y,mi_dx) = 1;
     end   
   end
end 
     
function cMaxSupport = BuildMaximumSupport(options);
%the support of the minimum are those pixels wich are within .5 of the appropriate line (L1)
for i = 1:options.NOrientations
   cMaxSupport{i} = zeros(options.MaxCellSupportSize);
   o(i) = (pi * (i-1)) / options.NOrientations - 1 * pi/2;
   u = [sin(-o(i)),cos(-o(i))];
   u = u / norm(u);
   for x = 1:options.MaxCellSupportSize 
      for y = 1:options.MaxCellSupportSize
        x2 = x - options.MaxCellSupportSize / 2 - .5;  
        y2 = y - options.MaxCellSupportSize / 2 - .5;  
        r = [y2, x2];
 	d = abs(r * u');
 	if(d < options.MaxThickness)
 	  cMaxSupport{i}(y,x) = 1;
 	end
      end
   end
end 
     
function ImgStack = ApplyDirectionalFilters(img,cfilters);
%function ImgStack = ApplyDirectionalFilters(img,cfilters);
nFilts = size(cfilters,2);
ImgStack = zeros(size(img,1), size(img,2),nFilts);
for i = 1:nFilts
  pimg = padimage(img,max(ceil(size(cfilters{i})/2)),'symmetric');
  pimg = abs(imfilter(pimg, cfilters{i}));
  ImgStack(:,:,i) = unpadimage(pimg,max(ceil(size(cfilters{i})/2))); 
end  
   
function M = ThresholdArray(M,options)
d.CutoffFraction = .75;
options = ResolveMissingOptions(options,d);
%% for i = 1:size(M,3)
%%   Monelayer = M(:,:,i);
%%   cutoff = percentile(Monelayer(:),options.CutoffFraction,1);
%%   Monelayer = Monelayer - cutoff;
%%   Monelayer = max(Monelayer,0);
%%   M(:,:,i) = Monelayer + cutoff;
%% end
cutoff = percentile(M(:),options.CutoffFraction,1);
M(find(M < cutoff)) = 0;
% M = M - cutoff;
% M = max(M,0);
% M = M + cutoff;

function M = SmoothOverOrientations(M,options)
d.OrientationSmoothingKernel = [.05 .25 .4 .25 .05];
options = ResolveMissingOptions(options,d);
k = options.OrientationSmoothingKernel;
M = imfilter(M,permute(k(:),[3,2,1]),'circular');

   
   
