function initialConditionsTest3D()
% initialConditionsTest3D: test initial condition routines in 3 dimensions.
%
%   initialConditionsTest3D (no arguments)
%
% This function (basically a script file) generates a sequence of 
%   shapes built by constructive solid geometry methods
%
% In 3D, it is rather hard to visualize the implicit surface function
%   (a scalar function of three variables).
%
% Consequently, we just show the implicit surface using a isosurface plot.

% Copyright 2004 Ian M. Mitchell (mitchell@cs.ubc.ca).
% This software is used, copied and distributed under the licensing 
%   agreement contained in the file LICENSE in the top directory of 
%   the distribution.
%
% Ian Mitchell, 6/23/04

run('../addPathToKernel');

%---------------------------------------------------------------------------
% Standard output file id.
fid = 1;

%---------------------------------------------------------------------------
fprintf(fid, '\n%s %s\n%s', ...
             'A sequence of implicit surface functions in 3D', ...
             'along with their descriptions.', ...
             'Press any key after each display to move to the next');

%---------------------------------------------------------------------------
% A 3D grid.
g3.dim = 3;
g3.min = -2;
g3.max = +2;
g3.N = 101;
g3 = processGrid(g3);

%---------------------------------------------------------------------------
% Some basic shapes.
sphere1 = shapeSphere(g3);
show3D(g3, sphere1);
fprintf(fid, '\nUnit circle at the origin, generated by shapeSphere');
pause;

sphere2 = shapeSphere(g3, [ +0.5; +0.5; +0.5 ], 2);
show3D(g3, sphere2);
fprintf(fid, '\nradius 2 sphere in octant 1 at [ 0.5 0.5 0.5 ], %s', ...
             'generated by shapeSphere');
pause;

cylinder1 = shapeCylinder(g3, 1, [ 0.0; 0.0; 1.0; ], 0.5);
show3D(g3, cylinder1);
fprintf(fid, '\nHorizontal cylinder radius 0.5, offset vertical 1.0, %s', ...
             'generated by shapeCylinder');
pause;

rectangle1 = shapeRectangleByCorners(g3);
show3D(g3, rectangle1);
fprintf(fid, '\nUnit square in octant 1, %s', ...
             'generated by shapeRectangleByCorners');
pause;

rectangle2 = shapeRectangleByCorners(g3, [ -Inf; -1.0;  0.0 ], ...
                                         [ +Inf; +1.0; +1.0 ]);
show3D(g3, rectangle2);
fprintf(fid, '\nHorizontal rectangular cross-section prism %s, %s', ...
             'offset vertically up', 'generated by shapeRectangleByCorners');
pause;

rectangle3 = shapeRectangleByCenter(g3, [ 0.0; 0.0; -0.5 ], [ Inf; 2.0; 1.0 ]);
show3D(g3, rectangle3);
fprintf(fid, '\nThe same prism, offset vertically down, %s', ...
             'generated by shapeRectangleByCenter');
pause;

hyperplane1 = shapeHyperplane(g3, [ 0.0; 3.0; 4.0 ], [ 1; 1; 1 ]);
show3D(g3, hyperplane1);
fprintf(fid, '\nNormal [ 0, 3, 4 ] through [ 1, 1, 1 ], %s', ...
             'generated by shapeHyperplane');
pause;

hyperplane2 = shapeHyperplane(g3, [ -1.0; +1.0; +1.0 ]);
show3D(g3, hyperplane2);
fprintf(fid, '\nNormal [ -1, +1, +1 ] through origin, %s', ...
             'generated by shapeHyperplane');
pause;

%---------------------------------------------------------------------------
% Constructive geometry operations.
union1 = shapeUnion(sphere1, rectangle1);
show3D(g3, union1);
fprintf(fid, '\nUnion of unit sphere and unit cube, %s', ...
             'generated by shapeUnion (NOTE: not signed distance)');
pause;

intersection1 = shapeIntersection(cylinder1, hyperplane2);
show3D(g3, intersection1);
fprintf(fid, '\nIntersection of cylinder and second hyperplane, %s', ...
             'generated by shapeIntersection');
pause;

difference1 = shapeDifference(rectangle3, sphere2);
show3D(g3, difference1);
fprintf(fid, '\nRectangular prism subtract second sphere, %s', ...
             'generated by shapeDifference');
pause;

%---------------------------------------------------------------------------
% Now let's build an octohedron (convex 3D polyhedron)
%   by intersections of hyperplanes.
plusMinusOne = [ -1 +1 ];
% We're building by intersections, so start with everything.
octo = -inf * ones(g3.shape);
for n1 = plusMinusOne
  for n2 = plusMinusOne
    for n3 = plusMinusOne
      normal = [ n1; n2; n3 ];
      point = [ n1; n2; n3 ] / 2;
      octo = shapeIntersection(octo, shapeHyperplane(g3, normal, point));
    end
  end
end

show3D(g3, octo);
fprintf(fid, '\nOctohedron centered at the origin, %s', ...
             'generated by combining shapeHyperplane and shapeIntersection');


%---------------------------------------------------------------------------
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%---------------------------------------------------------------------------
function h = show3D(g, data)
% show3D: display a 3D implicit surface function.
%
%   show3D(g, data)
%
% 3D implicit surface functions are hard to visualize directly, so
%   we'll just look at their zero isosurface.
%
% parameters:
%   g   	Grid structure (see processGrid.m for details).
%   data        Array containing the implicit surface function.

%---------------------------------------------------------------------------
% Keep track of isosurfaces from previous calls, and delete them.
persistent handle
if(~isempty(handle))
  delete(handle);
  clf;
end

%---------------------------------------------------------------------------
% Isosurface plot of implicit surface.
handle = patch(isosurface(g.xs{1}, g.xs{2}, g.xs{3}, data, 0));
set(handle, 'FaceColor', 'red', 'EdgeColor', 'none');
camlight left;  camlight right;
lighting phong;
view(3);
daspect([ 1 1 1 ]);
axis(g.axis);
grid on;
xlabel('x');  ylabel('y');  zlabel('z');
rotate3d on;

h = handle;
