function test3
% Segment an image by using a transitive grouping rule based
% upon edge pair likelihoods.
%
% Copyright (C) 2001-2009 Center for Perceptual Systems

% Read histogram
fn='kodak_counts_16_8_72_72.txt';
[same,diff]=readhist(fn,16,72,72);

% Get distance dependent priors
b=beta(same,diff);

% Create a table of distance dependent likelihoods
l=ddl(same,diff);

% Apply beta parameter
for i=1:16
    l(i,:)=l(i,:)/b(i);
end

% Create a colorful map
MAPSZ=40;
map=hsv(MAPSZ);

% Randomize the order
rand('seed',10);
map=map(randperm(MAPSZ),:);

% Read some edges from a data file
fn='testdata.txt';
e=load(fn);
x=e(:,1);
y=e(:,2);
th=mod(e(:,3)-pi/2,2*pi);
[x2,y2]=pol2cart(th,2);
x2=x+x2;
y2=y+y2;
g=e(:,5);

% Get a subset
d=hypot((x-200),(y-320));
i=find(d<64);
x=x(i); y=y(i); x2=x2(i); y2=y2(i);
th=th(i); g=g(i);

% Make group numbers sequential
g=renumber(g);

N=length(x(:));
fprintf('%d total edges\n',N);

% Show the original image
figure(1)
title('Hand segmented image');
plot(x,y,x2,y2,g,map);

% Generate indices
[d,phi,thd]=ndgrid([0:15]*8+4,[0:71]*2*pi/72,[0:71]*2*pi/72);

% Get weights
w=lookup(d,phi,thd,l,8);

% Threshold
T=2.0;
ind=find(w>T);
w=w(ind); d=d(ind); phi=phi(ind); thd=thd(ind);

% Clamp likelihoods to a range
w(find(w<0.01))=0.01;
w(find(w>100))=100;

% Put likelihoods on a log scale
w=log10(w);

% Map likelihoods to (0,1)
w=(w(:)+2)/4;

% Get indices of same and opposite polarity edge pairs
si=find(thd<pi/2|thd>=3*pi/2);
oi=find(thd>=pi/2&thd<3*pi/2);

% Get indices of edge pairs on left and right
ri=find(phi<pi/2|phi>=3*pi/2);
li=find(phi>=pi/2&phi<3*pi/2);

% Get indices of edges to plots
ind=[intersect(si,ri);intersect(oi,li)];

% Plot the rule used to segment
figure(2)
title(sprintf('Transitive grouping rule (likelihood threshold=%0.1f)',T));
plotgeometry(d(ind),phi(ind),thd(ind),w(ind),3.5,true);
xlim([-50 50])
ylim([-10 10])
drawnow

% Maintain a matrix of connections
connected=sparse(N,N);

% Allocate space
d=zeros(N,1);
phi=zeros(N,1);
thd=zeros(N,1);

% Determine which pairs are connected
fprintf('Getting edge pair connections\n');
for i=1:N
    for j=1:N
        [d(j),phi(j),thd(j)]=getgeometry(x(i),y(i),th(i),x(j),y(j),th(j));
    end
    % Lookup values from the likelihood table
    w=lookup(d,phi,thd,l,8);
    connected(i,find(w>T))=1;
end

% Get group numbers
disp 'Segmenting'
g=conncomp(connected);
g=renumber(g);

% Get indices of connected pairs
[i,j]=ind2sub([N N],find(connected));

% Show them
figure(3)
title(sprintf('Transitive grouping rule segmented image (likelihood threshold=%0.1f)',T));
% Draw the connections
%h=line([x(i)';x(j)'],[y(i)';y(j)']);
%set(h,'color',[0.8 0.8 0.8]);
plot(x,y,x2,y2,g,map);

function g=renumber(g)
% Renumber groups so that they are sequential

% Get rid of 0's
g=g+1;
% Renumber them
u=unique(g);
n=1;
for i=1:length(u(:))
    g(find(g==u(i)))=n;
    n=n+1;
end

function plot(x1,y1,x2,y2,g,map)
% Plot the lines
h=line([x1';x2'],[y1';y2']);
% Color by group
for j=1:length(x1)
    c=mod(g(j),size(map,1))+1;
    set(h(j),'color',map(c,:));
    set(h(j),'linewidth',2);
end
daspect([1 1 1])
[u,gh]=unique(g);
c=[];
for i=1:length(gh);
    c=[c;sprintf('group %3d',i)];
end
legend(h(gh),cellstr(c));
drawnow
