demo_exercise_12_14_circumscribing_tetrahedron

% Demo Barycentric coordinates of point cloud within minimal spanning
 tetrahedron

 see Kaseorg, A. (2014). 
 How do I find the side of the largest cube completely contained inside a 
 regular % tetrahedron of side s? 
 https://www.quora.com/How-do-I-find..., last visited 10/2015. 538

 see PCV exercise 12.14

 Wolfgang Förstner 6/2016
 wfoerstn@uni-bonn.de

 last changes: Susanne Wenzel 04/18
 wenzel@igg.uni-bonn.de

0001 %% Demo Barycentric coordinates of point cloud within minimal spanning
0002 % tetrahedron
0003 %
0004 % see Kaseorg, A. (2014).
0005 % How do I find the side of the largest cube completely contained inside a
0006 % regular % tetrahedron of side s?
0007 % https://www.quora.com/How-do-I-find..., last visited 10/2015. 538
0008 %
0009 % see PCV exercise 12.14
0010 %
0011 % Wolfgang Förstner 6/2016
0012 % wfoerstn@uni-bonn.de
0013 %
0014 % last changes: Susanne Wenzel 04/18
0015 % wenzel@igg.uni-bonn.de
0016 
0017 clc
0018 close all
0019 clearvars
0020 
0021 disp('--------------------------------------------')
0022 disp('   Barycentric coordinates of point cloud   ')
0023 disp('--------------------------------------------')
0024 
0025 % set range of coordinates
0026 % ddx=5;ddy=20;ddz=0.1;
0027 ddx=1;ddy=1;ddz=1;
0028 
0029 cube=0; % Cube given by params dy, dy, dz or unitcube
0030 
0031 if cube == 1
0032     disp('Cube [-1,+1]^3')
0033     X=[[-1,1,1,-1,-1,1,1,-1]*ddx;
0034        [1,1,-1,-1,1,1,-1,-1]*ddy;
0035        [-1,-1,-1,-1,1,1,1,1]*ddz];
0036 else
0037     disp(['Cube ',num2str(ddx),' x ',num2str(ddy),' x ',num2str(ddz)])
0038     N = 1000;
0039     X = [rand(1,N)*ddx;rand(1,N)*ddy;rand(1,N)*ddz];
0040 end
0041 
0042 % determine minimal bounding box
0043 xmin = min(X(1,:));
0044 ymin = min(X(2,:));
0045 zmin = min(X(3,:));
0046 xmax = max(X(1,:));
0047 ymax = max(X(2,:));
0048 zmax = max(X(3,:));
0049 mx = ( xmax + xmin )/2;
0050 my = ( ymax + ymin )/2;
0051 mz = ( zmax + zmin )/2;
0052 
0053 % four points of minimal spanning tetrahedron
0054 dx = xmax - xmin;
0055 dy = ymax - ymin;
0056 dz = zmax - zmin;
0057 u = ( 2+sqrt(2) )/6;
0058 v = u/sqrt(2);
0059 T = [mx-u*dx, my, mz-v*dz;...
0060      mx+u*dx, my, mz-v*dz;...
0061      mx, my-u*dy, mz+v*dz;...
0062      mx, my+u*dy, mz+v*dz]';
0063 disp('Corners of Tetrahedron')
0064 disp(num2str(T));
0065 
0066 Th = [T;ones(1,4)];
0067 
0068 % check barycentric coordinates
0069 alpha = inv(Th)*[X;ones(1,N)];                                             %#ok<MINV>
0070 
0071 % max alpha should be < 1
0072 disp(['maximum absolute Barycentric coordinate =',num2str(max((alpha(:))))])
0073 
0074 figure
0075 hold on
0076 plot3(X(1,:),X(2,:),X(3,:),'.b')
0077 plot3(  [ T(1,1), T(1,2), T(1,3), T(1,4), T(1,1), T(1,3), T(1,2), T(1,4)],...
0078         [ T(2,1), T(2,2), T(2,3), T(2,4), T(2,1), T(2,3), T(2,2), T(2,4)],...
0079         [ T(3,1), T(3,2), T(3,3), T(3,4), T(3,1), T(3,3), T(3,2), T(3,4)],...
0080         '-k','LineWidth',4' );
0081