%Script used to generate a "dumbell" shape to test the Isoperimetric 
%algorithm code on a non-trivial graph with a known solution.
%
%
%Note: Script requires installation of the MEXed version of Triangle.c in
%order to perform a "quality" triangulation, in which the interior points
%are filled in.  The MEXed version of Triangle.c may be obtained from:
%http://eslab.bu.edu/software/software.php
%After compiling, the triangle.m function must be present in a directory
%contained in your MATLAB path.
%
%Note: Code involves randomness, so results may be slightly different from
%       those in thesis.
%
%
%8/8/03

% Copyright (C) 2002, 2003 Leo Grady <lgrady@cns.bu.edu>
%   Computer Vision and Computational Neuroscience Lab
%   Department of Cognitive and Neural Systems
%   Boston University
%   Boston, MA  02215
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; either version 2
% of the License, or (at your option) any later version.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
%
% Date - $Id: fig1.m,v 1.1 2006/02/19 21:06:16 jonnyreb Exp $
%========================================================================%

%Initialization
radius=20;
offset=50;
sampling=10;
tubeWidth=10;
t = 0:2*pi/sampling:(2*pi*(1-1/sampling));

%Generate dumbbell graph
%--------------------------------------------------------------------------
%Left hand circle (centered at origin)
px1 = radius*cos(t);
py1 = radius*sin(t);
segments1=[1:sampling;[2:sampling,1]];

%Right hand circle (centered at offset)
px2 = radius*cos(t)+offset(1);
py2 = radius*sin(t);
segments2=segments1+10;

%Generate joining tube
tubeX=radius:((offset(1)-2*radius)/(sampling-1)):(offset(1)-radius);
tubeY=tubeWidth*ones(1,sampling)./2;
segments3=[1:(sampling-1);2:sampling]+2*sampling;
segments4=[1:(sampling-1);2:sampling]+3*sampling;
segments5=[2*sampling+1,3*sampling+1;2,sampling];
segments6=[3*sampling,4*sampling;sampling+sampling/2,sampling+sampling/2+2];

%Concatenation
x=[px1,px2,tubeX,tubeX];
y=[py1,py2,tubeY,-tubeY];

%Holes - So that convex hull (i.e. from one ball to another) does not join
sampleTmp=floor(sampling/2);
h=[tubeX(sampleTmp),tubeX(sampleTmp);2*tubeY(1),-2*tubeY(1)];

%Triangulation
tristruct.pl=[x;y];
tristruct.sl=[segments1,segments2,segments3,segments4,segments5,segments6];
tristruct.hl=h;
out=triangle(tristruct,'Qnepcqa2');  %Function requires MEXed triangle.c code
points=out.pl';
edges=out.el';
faces=out.tl';

%Generate adjacency matrix
W=adjacency(edges);

%Display dumbbell graph
figure
gplot(W,points,'k')
axis equal
axis tight



%Apply isoperimetric algorithm
%--------------------------------------------------------------------------
%Find four ground points
N=length(points);
ground1=round((N-1)*rand(1)+1); %Random ground
ground2=round((N-1)*rand(1)+1); %Random ground
xaxisIndex=find(abs(points(:,2))<1);
[tmp,ground3]=min((points(xaxisIndex,1)-offset/2).^2); %Center of neck
xaxisIndex2=xaxisIndex;
xaxisIndex2(ground3)=[];  %Find 2nd closest center
[tmp,ground4]=min((points(xaxisIndex2,1)-offset/2).^2); %Offset from center
ground3=xaxisIndex(ground3);
ground4=xaxisIndex2(ground4);

%Generate partitions of four points
[part1_1,part2_1,constant,xFunction1]=partitiongraph(W,0,0,0,ground1,1);
[part1_2,part2_2,constant,xFunction2]=partitiongraph(W,0,0,0,ground2,1);
[part1_3,part2_3,constant,xFunction3]=partitiongraph(W,0,0,0,ground3,1);
[part1_4,part2_4,constant,xFunction4]=partitiongraph(W,0,0,0,ground4,1);


%Display panels
figure
subplot(4,2,1) %Random ground, potentials
showmesh(points,[1-xFunction1],faces);
subplot(4,2,2) %Random ground, threshholded
xThresh1=ones(N,1); xThresh1(part2_1)=0;
showmesh(points,[xThresh1],faces);

subplot(4,2,3) %Random ground, potentials
showmesh(points,(1-xFunction2),faces);
subplot(4,2,4) %Random ground, threshholded
xThresh2=ones(N,1); xThresh2(part2_2)=0;
showmesh(points,xThresh2,faces);

subplot(4,2,5) %Center ground, potentials
showmesh(points,(1-xFunction3),faces);
subplot(4,2,6) %Center ground, threshholded
xThresh3=ones(N,1); xThresh3(part2_3)=0;
showmesh(points,xThresh3,faces);

subplot(4,2,7) %Just off-center ground, potentials
showmesh(points,(1-xFunction4),faces);
subplot(4,2,8) %Just off-center ground, threshholded
xThresh4=ones(N,1); xThresh4(part2_4)=0;
showmesh(points,xThresh4,faces);