function ConvertToContinuity(oldfn,newfn)
%function ConvertToC1(oldfn,newfn)
%
% reads the data as 81x81 columns in one matrix 
% returns a new matrix in Cont format
% only retain levels 4 5 and 6

load(oldfn,'X');%-->X

options.OrigSize = [81 81];
saMat = SSDB_Continuity(X,options);
clear X;
X = saMat.HStanLevel6;
X = [X;saMat.HStanLevel5];
X = [X;saMat.HStanLevel4];
% X = [X;saMat.HStanLevel3];
% X = [X;saMat.HStanLevel2];
save(newfn,'X');



function saMat = SSDB_Continuity(X,options);
%an image representation which makes explicit the line segments

d.OrigSize = [128 128];
d.IncludeStageOne = 0;
d.MaxIterations = 6;
d.nUpperLimitFeatPerGroup = inf;  
if(nargin < 2), options = [];, end;
options = ResolveMissingOptions(options,d);

TotalSize = 18*prod(options.OrigSize);
sz = options.OrigSize;


if(options.IncludeStageOne)
  saMat.HStanLevel1 = zeros(prod(sz)*18 ,size(X,2));
end
sz = ceil(sz/2);
if(options.MaxIterations >= 2)
  saMat.HStanLevel2 = zeros(prod(sz)*18 ,size(X,2));
end
sz = ceil(sz/2);
if(options.MaxIterations >= 3)
  saMat.HStanLevel3 = zeros(prod(sz)*18 ,size(X,2));
end
sz = ceil(sz/2);
if(options.MaxIterations >= 4)
  saMat.HStanLevel4 = zeros(prod(sz)*18 ,size(X,2));
end
sz = ceil(sz/2);
if(options.MaxIterations >= 5)
  saMat.HStanLevel5 = zeros(prod(sz)*18 ,size(X,2));
end
sz = ceil(sz/2);
if(options.MaxIterations >= 6)
  saMat.HStanLevel6 = zeros(prod(sz)*18 ,size(X,2));
end
sz = ceil(sz/2);
for imgIdx = 1:size(X,2)
  fprintf('%d of %d \r',imgIdx, size(X,2));
  img = reshape(X(:,imgIdx),options.OrigSize);
  cHS = HStanMultiOrientation(img,options);
  if(options.IncludeStageOne)
    saMat.HStanLevel1(:,imgIdx) = cHS{1}(:);
  end
  if(options.MaxIterations >= 2)
    saMat.HStanLevel2(:,imgIdx) = cHS{2}(:);
  end
  if(options.MaxIterations >= 3)
    saMat.HStanLevel3(:,imgIdx) = cHS{3}(:);
  end
  if(options.MaxIterations >= 4)
    saMat.HStanLevel4(:,imgIdx) = cHS{4}(:);
  end
  if(options.MaxIterations >= 5)
    saMat.HStanLevel5(:,imgIdx) = cHS{5}(:);
  end
  if(options.MaxIterations >= 6)
    saMat.HStanLevel6(:,imgIdx) = cHS{6}(:);
  end
end
if(options.nUpperLimitFeatPerGroup)
  if(options.IncludeStageOne)
    if(size(saMat.HStanLevel1,1) > options.nUpperLimitFeatPerGroup)
      span = ceil(size(saMat.HStanLevel1,1)/options.nUpperLimitFeatPerGroup);
      saMat.HStanLevel1 = saMat.HStanLevel1(1:span:end,:);
    end
  end
  if(size(saMat.HStanLevel2,1) > options.nUpperLimitFeatPerGroup)
    span = ceil(size(saMat.HStanLevel2,1)/options.nUpperLimitFeatPerGroup);
    saMat.HStanLevel2 = saMat.HStanLevel2(1:span:end,:);
  end
  if(size(saMat.HStanLevel3,1) > options.nUpperLimitFeatPerGroup)
    span = ceil(size(saMat.HStanLevel3,1)/options.nUpperLimitFeatPerGroup);
    saMat.HStanLevel3 = saMat.HStanLevel3(1:span:end,:);
  end
  if(size(saMat.HStanLevel4,1) > options.nUpperLimitFeatPerGroup)
    span = ceil(size(saMat.HStanLevel4,1)/options.nUpperLimitFeatPerGroup);
    saMat.HStanLevel4 = saMat.HStanLevel4(1:span:end,:);
  end
  if(size(saMat.HStanLevel5,1) > options.nUpperLimitFeatPerGroup)
    span = ceil(size(saMat.HStanLevel5,1)/options.nUpperLimitFeatPerGroup);
    saMat.HStanLevel5 = saMat.HStanLevel5(1:span:end,:);
  end
  if(size(saMat.HStanLevel6,1) > options.nUpperLimitFeatPerGroup)
    span = ceil(size(saMat.HStanLevel6,1)/options.nUpperLimitFeatPerGroup);
    saMat.HStanLevel6 = saMat.HStanLevel6(1:span:end,:);
  end
end

function [cHStanArray,cMaxFilt, cMinFilt,cMinSupport,cMaxSupport] = HStanMultiOrientation(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.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;
if(nargin < 2), options = [];, end;
options = ResolveMissingOptions(options,D);
cfilters = BuildDirectionalDerivatives(options);
cMinSupport = BuildMinimumSupport(options);
cMaxSupport = BuildMaximumSupport(options);
cHStanArray{1} = ApplyDirectionalFilters(img,cfilters);
if(options.ShowFigures)
  h = figure;
end
for i = 1:options.MaxIterations
  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(FourLayerColor(cHStanArray{i},[4]));
  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 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 = 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');

function I = maxfilter_arbitrarySupport(I,Msupport);
%function I = maxfilter_arbitrarySupport(I,Msupport);
%
%support must be a binary matrix.  The centroid of the support will be at
%[floor(size(support,1)/2), floor(size(support,2)/2)]

sz = size(Msupport);
centroid = floor(sz);
fsupport = find(Msupport);
ndxsupport = NDind2sub(sz,fsupport);
nPoints = size(fsupport,1);
ndxsupport = ndxsupport - repmat(centroid, [nPoints,1]);
padamnt = max(ceil(sz/2));
ndxsupport = ndxsupport + padamnt;
ndxsupport = ndxsupport + padamnt;
padI = padimage(I,padamnt,-inf);
BigMat = zeros(nPoints, prod(size(I)));
sI = size(I);
for i = 1:nPoints
  crop = padI( ndxsupport(i,1) +  (1:sI(1))-1,ndxsupport(i,2) + (1:sI(2)) - 1 );
  BigMat(i,:) = crop(:);
end
I = reshape(max(BigMat,[],1),sI);
   
   
