graphics_toolkit("fltk");

nGC = 10;
GCs = load ("GC.pointset.dat" );
PFs = load ("GC.ParallelFibers.0.section.dat" );
AAs = load ("GC.AscendingAxons.0.section.dat" );

pfpoints = (size (PFs)(2)-1)/3
aapoints = (size (AAs)(2)-1)/3
myGCinds = randi (size(GCs)(1),nGC,1)

plot3 (GCs(myGCinds,1),GCs(myGCinds,2),GCs(myGCinds,3),'*');
hold on;
for i=1:nGC
  id = myGCinds(i)-1
  myPFinds = find (PFs(:,1) == id)
  myPFs = PFs(myPFinds,:);
  for i=1:pfpoints
      x = myPFs(:,(i-1)*3+2);
      y = myPFs(:,(i-1)*3+3);
      z = myPFs(:,(i-1)*3+4);
      plot3(x,y,z,'g.');
  endfor       
  myAAinds = find (AAs(:,1) == id)
  myAAs = AAs(myAAinds,:);
  for i=1:aapoints
      x = myAAs(:,(i-1)*3+2);
      y = myAAs(:,(i-1)*3+3);
      z = myAAs(:,(i-1)*3+4);
      plot3(x,y,z,'g.');
  endfor       
endfor