function [saMat,imgstruct] = SSDB_CircleHoughTransform(X,options);
% performs a crude hough transform of the image.

d.OrigSize = [128 128];
d.RadMin = 2;
d.RadMax = 32;
d.RadGrain = sqrt(2);
d.CentGradPerRad = 1;
d.StopForDebugging = 0;
d.LinearVersion = 0;
d.SigmaMultiplier = [1,1];
d.DecoupleLeftRight = 0;
if(nargin < 2),options = [];,end
options = ResolveMissingOptions(options,d);
nIm = size(X,2);
% saMat.SegFeature = zeros(size(X));
for i = 1:nIm;
  fprintf('SSDB_CircleHoughTransform image %d of %d\r',i,nIm);
  stim = FourWayFilter(reshape(X(:,i),options.OrigSize));
  stim = OnlyRetainMaximumOrientationStim(stim);
  % keyboard
  cHoughCirclesPreCombo = DetectCircles(stim,options);
  if(options.StopForDebugging)
    keyboard
  end
  if(options.LinearVersion)
     for j = 1:length(cHoughCirclesPreCombo)
       cHoughCircles{j} = sum(cHoughCirclesPreCombo{j},3);
       if((nargout > 1) && (i==1))
        imgstruct{j} = cHoughCircles{j};
       end
       cHoughCircles{j} = imresize(cHoughCircles{j},1/8,'bilinear');
     end
  else
    for j = 1:length(cHoughCirclesPreCombo)
      z = sort(cHoughCirclesPreCombo{j},3);
      cHoughCircles{j} = sum(z(:,:,1:3),3);
      if((nargout > 1) && (i==1))
        imgstruct{j} = cHoughCircles{j};
       end
      cHoughCircles{j} = imresize(cHoughCircles{j},1/8,'bilinear');
    end
  end
  [saMat.SegFeature(:,i)] = vectorizeCellArray(cHoughCircles);
  % if(i == 10)
  %   keyboard;
  % end
end
fprintf('SSDB_CircleHoughTransform image %d of %d\n',i,nIm);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function MLI = UniqueWindowed(IdxImg,yx_wind);
%function MLI = UniqueWindowed(IdxImg,yx_wind);

IdxImg = padimage(IdxImg,[floor(yx_wind),ceil(yx_wind)],'symmetric');
u = unique(IdxImg);
MLI = zeros(size(IdxImg,1), size(IdxImg,2), length(u));

nWindsY = size(IdxImg,1) - yx_wind(1) + 1;
nWindsX = size(IdxImg,2) - yx_wind(2) + 1;
T = zeros(size(IdxImg,1), size(IdxImg,2), length(u));
for i = u'
  T(:,:,i) = (IdxImg == u(i));
  MLI(:,:,i) = maxfilter(T(:,:,i),[floor(yx_wind/2),ceil(yx_wind/2)]);
  %  MLI(:,:,i) = imdilate(T(:,:,i),ones(floor(yx_wind(1)),ceil(yx_wind(2))));
end
MLI = unpadimage(MLI,[floor(yx_wind),ceil(yx_wind)]);

function outim = FourWayFilter(inim)
filters{1} = [ 1  1  1; 0  0  0;-1 -1 -1];
filters{2} = [-1  0  1;-1  0  1;-1  0  1];
filters{3} = [ 0  1  1;-1  0  1;-1 -1  0];
filters{4} = [ 1  1  0; 1  0 -1; 0 -1 -1];
outim = abs(ApplyFilterBank(inim,filters));

function cHoughCircles = DetectCircles(stim,options);
% function cHoughCircles = DetectCircles(stim,options);
% each oriented element contributes to potential circles at different scales.
rads = options.RadMin;
orientations = [90,0,45,135]; 
translations{1} = [1, 0];
translations{2} = [0, 1];
translations{3} = [.7, -.7];
translations{4} = [.7, .7];
while(rads(end) < options.RadMax);
  rads(end+1) = rads(end) * options.RadGrain;
end
iRad = 0;
for r = rads
   iRad = iRad + 1;
   cHoughCircles{iRad} = zeros(options.OrigSize(1),options.OrigSize(2),4*2^(options.DecoupleLeftRight));
   iori = 0;
   for ori = orientations
      iori = iori + 1;
      Translation = r * translations{iori};
      ConvK = BuildConvolutionKernel(r,ori,options);
      convIm = conv2(stim(:,:,iori),ConvK,'same');
      if(options.DecoupleLeftRight) ;
        cHoughCircles{iRad}(:,:,iori) = cHoughCircles{iRad}(:,:,iori) + TranslateIm(convIm, Translation);
        cHoughCircles{iRad}(:,:,iori+4) = cHoughCircles{iRad}(:,:,iori) + TranslateIm(convIm, -Translation);
      else
        cHoughCircles{iRad}(:,:,iori) = cHoughCircles{iRad}(:,:,iori) + TranslateIm(convIm, Translation);
        cHoughCircles{iRad}(:,:,iori) = cHoughCircles{iRad}(:,:,iori) + TranslateIm(convIm, -Translation);
      end
   end
end

function ConvK = BuildConvolutionKernel(r,ori,options)
%function ConvK = BuildConvolutionKernel(r,ori)
r = round(r);
ConvK = GausImg([2*r 2*r],[r,r],[4*r,r*1.75].*options.SigmaMultiplier);
% ConvK = GausImg([2*r 2*r],[r,r],[8*r,r*.15]);
% ConvK = ConvK > .5;
ConvK = 2 * min(ConvK, .5);
ConvK = im2double(imrotate(ConvK, ori,'bilinear','crop'));
ConvK = ConvK / (ConvK(:)' * ConvK(:) );


function t = TranslateIm(convIm, Translation);
%function t = TranslateIm(convIm, Translation);
Translation = round(Translation);
t = zeros(size(convIm));
t = circshift(convIm,[Translation]);
Translation = min(Translation,size(convIm));
Translation = max(Translation,-size(convIm) + 1);
if(Translation(1) < 0);
  t((end+Translation(1)):(end),:) = 0;
else
  t(1:Translation(1),:) = 0;
end
if(Translation(2) < 0);
  t(:,(end+Translation(2)):(end)) = 0;
else
  t(:,1:Translation(2)) = 0;
end

function stim = OnlyRetainMaximumOrientationStim(stim);
% function stim = OnlyRetainMaximumOrientationStim(stim);
M = max(stim(:));
stim = stim + rand(size(stim))*.01*M;
[s, si] = sort(stim,3);
v = repmat(s(:,:,end),[1,1,size(stim,3)]);
f = stim == v;
stim = 0 * stim;
stim(f) = v(f);


