Calculating volume of a convex hull

I want to calculate convex hull volumes of various sets of points in different dimensions.

For example, if I'm working with q dimensional space, I want to calculate the volume of convex hull which is a (q-1) dimensional object. However, I receive error since the object that I want to calculate the volume of, is not full dimensional. Is there any way that I can calculate those volumes?

I inserted an example set of points.

x=[1 2 1; 3 2 2 ;1 1 1 ; 3 1 2]
[k,vol]=convhulln(x)

The error I receive is this:

Error using builtin

QH6154 qhull precision error: initial facet 3 is coplanar with the interior point

...

Answer 1

[Ayush]

Hi Sitek,

You are getting that error because all your points are coplanar, and you are trying to find the convex hull in 3D of a 2D figure.

As the error message displays. the underlying library (qhull) gives us some suggestions in case your data is lower-dimensional. As suggested in the error message, you may use the option 'QJ' to joggle the input by using convhulln(x,{'QJ'}).

Read more about it here: http://www.qhull.org/html/qh-optq.htm#QJn

There is an alternative workaround as well. Since the points are co-planar, you can essentially treat the problem as a 2D convex hull problem by projecting the points onto a 2D plane and then using a 2D convex hull algorithm.

% Defining the points matrix
x = [1 2 1; 3 2 2; 1 1 1; 3 1 2];
% STEP 1
% a. Verify co-planarity by calculating the volume of the tetrahedron
A = x(1,:);
B = x(2,:);
C = x(3,:);
D = x(4,:);
AB = B - A;
AC = C - A;
AD = D - A;
% b. Calculate the volume of the parallelepiped using scalar triple product
volume = dot(AB, cross(AC, AD));
% c. Use tolerance to account for floating-point errors, if co-planar goto STEP 2
if abs(volume) < 1e-7
    disp('The points are co-planar.');
    
    % STEP 2: Project the points onto a 2D plane and find convex hull
    % d. Choose the first point A as the reference point and AB, AC as basis vectors
    U = AB / norm(AB);
    V = AC - dot(AC, U) * U;
    V = V / norm(V);
    
    % e. Initialize the projected 2D points array
    projectedPoints = zeros(size(x, 1), 2);
    
    for i = 1:size(x, 1)
        W = x(i,:) - A;
        % f. Project point onto the new basis
        projectedPoints(i,1) = dot(W, U);
        projectedPoints(i,2) = dot(W, V);
    end
    
    % f. Compute the convex hull in 2D
    k = convhull(projectedPoints(:,1), projectedPoints(:,2));
    
    % g. Map the convex hull indices back to the original 3D points
    convexHull3D = x(k, :);
    convexHull3D = unique(convexHull3D, 'rows');
    
    % h. Display the results
    disp('The convex hull in 2D (indices of the original points):');
    disp(k);
    
    disp('The convex hull in 3D space:');
    disp(convexHull3D);
    
else
    disp('The points are not co-planar.');
end

Hope this was helpul!