function caCircs = CirclesFromOneContinuity(stim,options)
%function caCircs = CirclesFromOneContinuity(stim,options)
%
%finds circles from the continuity data

d.RadMin = 2;
d.RadMax = 32;
d.RadGrain = sqrt(2);
d.StopForDebugging = 0;
d.SigmaMultiplier = [1,1];
d.DecoupleLeftRight = 0;
d.EliminateTop10PercentOfOrientations = 1;
d.OrigSize = [size(stim,1), size(stim,2)];

if(nargin < 2),options = [];,end
options = ResolveMissingOptions(options,d);
  
stim = OnlyRetainMaximumOrientationStim(stim);
cHoughCirclesPreCombo = DetectCircles(stim,options);
if(options.StopForDebugging)
  keyboard
end
caCircs = CombineOrientations(cHoughCirclesPreCombo,options);    

function caCircles = CombineOrientations(cH,options)
nori = size(cH{1},3);
nearend = nori+1;
if(options.EliminateTop10PercentOfOrientations)
   nearend = nori - floor(.1 * nori);
end
for j = 1:length(cH)
   z = sort(cH{j},3);
   z(:,:,nearend:end) = 0;
   % cHoughCircles{j} = sum(z(:,:,1:3),3);
   caCircles{j} = sum(z,3);
end


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



function cHoughCircles = DetectCircles(stim,options);
% function cHoughCircles = DetectCircles(stim,options);
% each oriented element contributes to potential circles at different scales.
rads = options.RadMin;
orientations = 0:10:170; 
nori = length(orientations);
for i = 1:nori
  translations{i} = [cos(deg2rad(orientations(i))), sin(deg2rad(orientations(i)))];
end
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),nori*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+nori) = 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);


