function slide=qhulldemo
% QHULLDEMO Tessellation and interpolation of scattered data.
%   QHULLDEMO describes convex hulls, Delaunay tessellations,
%   and Voronoi diagrams in 3 dimensions. It also shows how
%   to interpolate 3 dimensional scattered data.

% Copyright 1984-2000 The MathWorks, Inc.
% $Revision: 1.2 $  $Date: 2000/09/11 16:19:08 $

if nargout<1,
  playshow qhulldemo
else
  %========== Slide 1 ==========

  slide(1).code={
   'n = sprintf(''\n'');',
   's =      ''if isempty(getappdata(gcf,''''lp'''')),'';',
   's = [s n ''  if isempty(get(hittest,''''buttondownfcn'''')),''];',
   's = [s n ''    set(gcf,''''pointer'''',''''arrow'''');''];',
   's = [s n ''  else,''];',
   's = [s n ''    set(gcf,''''pointer'''',''''fleur'''');''];',
   's = [s n ''  end''];',
   's = [s n ''else''];',
   's = [s n ''  set(gcf,''''pointer'''',''''fleur'''');''];',
   's = [s n ''  lp=getappdata(gcf,''''lp'''');''];',
   's = [s n ''  cp=get(gcf,''''currentpoint'''');''];',
   's = [s n ''  diff=lp-cp;camorbit(diff(1),diff(2));''];',
   's = [s n ''  setappdata(gcf,''''lp'''',cp);''];',
   's = [s n ''end''];',
   'set(gcf,''windowbuttonmotionfcn'',s);',
   'set(gcf, ''WindowButtonUpFcn'', ''setappdata(gcf,''''lp'''',[]);'')',
   'colormap(spring)',
   'cla reset;',
   'hold on;',
   'n = 50; X = randn(n,3); C = convhulln(X);',
   'h = [];',
   'for i = 1:size(C,1)',
   '  j= C(i,[1 2 3 1]);',
   '  h(i)=patch(X(j,1),X(j,2),X(j,3),rand,''FaceAlpha'',0.6);',
   'end',
   'h = [h plot3(X(:,1),X(:,2),X(:,3),''ko'',''markerfacecolor'',''k'')];',
   'view(3), axis equal off tight vis3d; camzoom(1.5)',
   'set(h, ''buttondownfcn'',''setappdata(gcf,''''lp'''',get(gcf,''''currentpoint''''));'');',
   'set(findobj(''tag'',''comments''),''backgroundcolor'',''k''); set(gcf,''visible'',''on''); drawnow; pause(0);set(findobj(''tag'',''comments''),''backgroundcolor'',''w'')'};
    %  The line above is to work around an OpenGL bug (G88236).  Remove when fixed.  Without the pause(0), it doesn't work on Unix.
  slide(1).text={
   'Many applications in science, engineering, statistics, and mathematics use structures like convex hulls, Delaunay tessellations and Voronoi diagrams for analyzing data. MATLAB enables you to geometrically analyze data sets in any dimension.',
   '',
   'Here, in 3 dimensions, we show a set of 50 points with its convex hull. Click and drag this 3-D shape, or any other shown in this demo, to rotate it.',
   ''};

  %========== Slide 2 ==========

  slide(2).code={
   'cla reset; hold on',
   'x=zeros(8,3); x([5:8,11,12,15,16,18,20,22,24])=1;',
   'd = [1 2 4 3 1 5 6 8 7 5 6 2 4 8 7 3];',
   'h = plot3(x(d,1),x(d,2),x(d,3),''b:'');',
   'x(9,:) = [0.5 0.5 0.5];     % Add center point',
   'h = [h plot3(x(:,1),x(:,2),x(:,3),''b.'',''markersize'',20)];',
   't = text(x(:,1),x(:,2),x(:,3), num2str((1:9)''));',
   ' set(t,''verticalalignment'',''bottom'',''fontweight'',''bold'', ''fontsize'',12);',
   'view(3); axis equal tight off vis3d',
   'set(h, ''buttondownfcn'',''setappdata(gcf,''''lp'''',get(gcf,''''currentpoint''''));'');',
   'camorbit(10,0);',
   'set(gcf,''DoubleBuffer'',''off'');drawnow;set(gcf,''DoubleBuffer'',''on'')',
   ''};
   % without the drawnow above, the doublebuffering doesn't seem to kick in properly.
  slide(2).text={
   'We can create a data set  X of the 8 vertices of a cube plus its center.',
   '',
   '>> X = [ 0 0 0 ; 0 0 1 ; 0 1 0 ; 0  1  1 ;  1 0 0 ; 1 0 1 ; 1 1 0 ; 1 1 1 ; ',
    '       0.5 0.5 0.5];',
   '',
   'X is a 9-by-3 matrix where each row is the 3-D coordinates of one point.'};

  %========== Slide 3 ==========

  slide(3).code={
   'cla reset; hold on',
   'tri = convhulln(x);',
   'h = [];',
   'for i = 1:size(tri,1)',
   '  c = tri(i,[1 2 3 1]);',
   '  h(i) = patch(x(c,1),x(c,2),x(c,3),i,''FaceAlpha'', 0.9);',
   'end',
   'h = [h plot3(x(:,1),x(:,2),x(:,3),''ko'',''markerfacecolor'',''k'')];',
   'view(3); axis equal tight off vis3d',
   'set(h, ''buttondownfcn'',''setappdata(gcf,''''lp'''',get(gcf,''''currentpoint''''));'');',
   '' };
  slide(3).text={
   'The convex hull of a data set is the smallest convex region that contains the data set.',
   '',
   'The convex hull of the cube data set X can be computed by',
   '>> K = convhulln(X);',
   '',
   'For this data set X, the convex hull has 12 facets, each corresponding to a row in K and plotted above with X. The cube is transparent so that you can see all the facets and data points.'};

  %========== Slide 4 ==========

  slide(4).code={
   '% Slide 5',
   'cla reset',
   'hold on',
   'tri=delaunayn(x);',
   'h=[];',
   'for i = 1:size(tri,1)',
   '  y = tri(i,[1 1 1 2; 2 2 3 3; 3 4 4 4]);',
   '  x1 = reshape(x(y,1),3,4);',
   '  x2 = reshape(x(y,2),3,4);',
   '  x3 = reshape(x(y,3),3,4);',
   '  h(i)=patch(x1,x2,x3,(1:4)*i,''FaceAlpha'',0.8);',
   'end',
   'h = [h plot3(x(:,1),x(:,2),x(:,3),''ko'',''markerfacecolor'',''k'')];',
   'view(3); axis equal tight off vis3d',
   'set(h, ''buttondownfcn'',''setappdata(gcf,''''lp'''',get(gcf,''''currentpoint''''));'');',
   'camorbit(10,0)' };
  slide(4).text={
   'A Delaunay tessellation in 3 dimensions is a set of tetrahedrons such that no data points are contained in any tetrahedron''s circumsphere.',
   '',
   'The Delaunay tessellation of the data set X can be computed by',
   '>> T = delaunayn(X);',
   '',
   'The 12 rows of T represent the 12 tetrahedrons that partition the data set X.'};

  %========== Slide 5 ==========

  slide(5).code={
   '% Slide 5',
   'cla reset; hold on',
   '[c,v] = voronoin(x);',
   'nx = c(v{9},:);',
   'tri = convhulln(nx);',
   'h=[];',
   'for i = 1:size(tri,1)',
   'h(i)=patch(nx(tri(i,:),1),nx(tri(i,:),2),nx(tri(i,:),3),rand,''FaceAlpha'',0.8);',
   'end',
   'h = [h plot3(x(d,1),x(d,2),x(d,3),''b:.'',''markersize'',20)];',
   'h = [h plot3(x(9,1),x(9,2),x(9,3),''k.'',''markersize'',20)];',
   'view(3); axis equal tight off vis3d; camzoom(1.5)',
   'set(h,''buttondownfcn'',''setappdata(gcf,''''lp'''',get(gcf,''''currentpoint''''));'');',
   'camorbit(20,0)',
   '' };
  slide(5).text={
   'A Voronoi diagram partitions the data space into polyhedral regions, with one region for each data point.  Anywhere within a region is closer to its data point than any other in the set.',
   'The Voronoi diagram of the cube data set X can be computed by',
   '>> [V,C] = voronoin(X);',
   'V is the set of Voronoi vertices. C represents the set of Voronoi regions.  For our data set X, C has 9 Voronoi regions.  Here we show one Voronoi region, the region for the center point of the cube.'};

  %========== Slide 6 ==========

  slide(6).code={
   'cla reset; hold on;',
   'X = 2*rand(500,3)-1;',
   'v = sum(X.^2,2);',
   'h = plot3([0.02, 0.47],[0.02,0.57],[0.02,0.57],''k-'',''linewidth'',2);',
   'h = [h plot3(0,0,0,''bo'',''markerfacecolor'',''b'')];',
   'h = [h plot3(x(d,1),x(d,2),x(d,3),''r:'')];',
   'h = [ h plot3(0.5,0.6,0.6,''go'',''markerfacecolor'',''g'')];',
   'text(0.02,-0.2,0,''the origin'',''fontsize'',12,''fontweight'',''bold'');',
   'text(0.55,0.6,0.6,''a point in X'',''fontsize'',12,''fontweight'',''bold'');',
   'text(0.28,0.3,0.35,''sqrt(v)'',''fontsize'',12,''fontweight'',''bold'');',
   'view(3); axis equal tight off vis3d;',
   'set(h,''buttondownfcn'',''setappdata(gcf,''''lp'''',get(gcf,''''currentpoint''''));'');',
   'camorbit(20,-10);',
   'set(gcf,''DoubleBuffer'',''off'');drawnow;set(gcf,''DoubleBuffer'',''on'')' };
   % again, double buffering doesn't seem to "kick in" without the above workaround.
  slide(6).text={
   'GRIDDATAN interpolates multidimensional scattered data. It uses DELAUNAYN to tessellate the data, and then interpolates based on the tessellation.',
   'We start with a data set of 500 random points in 3 dimensions and compute the values of a function, the squared distance from the origin, at each of these points.',
   '>> n = 500;',
   '>> X = 2*rand(n,3)-1;',
   '>> v = sum(X.^2,2);'};

  %========== Slide 7 ==========

  slide(7).code={
   'cla reset; hold on;',
   'setptr(gcf,''watch'');',
   'd  = -0.8:0.2:0.8;',
   '[x0,y0,z0] = meshgrid(d,d,d);',
   'X0 = [x0(:) y0(:) z0(:)];',
   'v0 = reshape(griddatan(X,v,X0),size(x0));',
   '',
   'h = plot3(X(:,1),X(:,2),X(:,3),''k+'',''markerfacecolor'',''k'');',
   'h = [h plot3(X0(:,1),X0(:,2),X0(:,3),''r.'',''markerfacecolor'',''r'')];',
   '',
   'view(3); axis equal tight off vis3d; camzoom(1.6);',
   'set(h, ''buttondownfcn'',''setappdata(gcf,''''lp'''',get(gcf,''''currentpoint''''));'');',
   'set(h,''erasemode'',''background'');',
   'set(gcf,''DoubleBuffer'',''off''); drawnow; set(gcf,''DoubleBuffer'',''on'')',
   'setptr(gcf,''arrow'')' };
   % again, double buffering doesn't seem to "kick in" without the above workaround.  Needed when going backwards.
  slide(7).text={
   'With GRIDDATAN we can interpolate X and the values v over a grid X0 to get the function values v0 over this grid.',
   '>> d  = -0.8:0.2:0.8;',
   '>> [x0,y0,z0] = meshgrid(d,d,d);',
   '>> X0 = [x0(:) y0(:) z0(:)];',
   '>> v0 = griddatan(X, v, X0);',
   '>> v0 = reshape(v0, size(x0));',
   '',
   'The black points are X and the red points are the grid X0.',
   ''};

  %========== Slide 8 ==========

  slide(8).code={
   'cla reset; hold on;',
   'h = patch(isosurface(x0,y0,z0,v0,0.6),''FaceColor'',''red'',''EdgeColor'',''none'');',
   'isonormals(x0,y0,z0,v0,h);',
   'view(3); axis equal tight off vis3d; camzoom(1.6)',
   'camlight; lighting phong',
   'set(h, ''buttondownfcn'',''setappdata(gcf,''''lp'''',get(gcf,''''currentpoint''''));'');',
   '' };
  slide(8).text={
   'To visualize the surface at all points where the function takes on a constant value, we can use ISOSURFACE and ISONORMALS. Since the function is the squared distance from the origin, the surface at a constant value is a sphere.',
   '>> c=0.6;  % constant value',
   '>> h=patch(isosurface(x0,y0,z0,v0,c),''FaceColor'',''r'',''EdgeColor'',''n'');',
   '>> hold on; isonormals(x0,y0,z0,v0,h);',
   '>> view(3); axis equal tight off vis3d; camzoom(1.6)',
   '>> camlight; lighting phong'};

  %========== Slide 9 ==========

  slide(9).code={
   'cla reset',
   'load qhulldemo',
   'd = -0.8:0.05:0.8;',
   '[x0,y0,z0] = meshgrid(d,d,d);',
   'h = patch(isosurface(x0,y0,z0,v0,0.6));',
   'isonormals(x0,y0,z0,v0,h);',
   'set(h,''FaceColor'',''red'',''EdgeColor'',''none'');',
   'view(3); axis equal tight off vis3d; camzoom(1.6)',
   'camlight; lighting phong',
   'set(h, ''buttondownfcn'',''setappdata(gcf,''''lp'''',get(gcf,''''currentpoint''''));'');',
   '' };
  slide(9).text={
   'With more data points in X and a denser grid X0, the sphere is smoother but takes longer to compute.',
   '',
   'Here is a precomputed sphere generated using 5000 data points in X and a distance between gridpoints of 0.05.'};
end
