Getting the points belonging to the convex hull

Question

I have a binary image of one granule in Matlab. I can find the convex hull of the granule with the following function:

[K, V] = convhull(granule);

How can I find all pixels which belong to the convex hull, but don't belong to the granule in the original image? I mean I'd like to do something like that:

granule2 = zeros(size(granule));
granule2(K == 1 & granule == 0) = 2;

It doesn't work, because K is of size x by 3, where x is the number of triangles in the convex hull.

Edit: according to the documentation, the convex hull should be an array with the indexes of points making up the convex hull in each row. So how can I find all points which are inside of the volume determined by these points.

Edit2: Let me put it in another words: I have an image which is a 3D array of points. It's not a sphere and it has some indents (so the convex hull doesn't lay on the surface of my image).

I want to find the convex hull and after that find all points which are inside the convex hull, but are outside the granule. Here is how it would look like in 2D (I want to find the red pixels):

Edit3: NicolaSysnet, Your algorithm should return all the pixels (their indexes) which are red in my picture (the picture is in 2D,because it was easier to draw).

Answer 1

[K, V] = convhull(granule);
granule2 = zeros(size(granule));
tmp = granule(K,:)==0; %// all points on the convex hull but not in the granule
tmp2 = tmp(:,1)==tmp(:,2); %// both indices of your granule are zero
granule2(tmp(tmp2)) = 2;

K are the row numbers of your points corresponding to the points on the convex hull, V is just the volume spanned by that convex hull. Thus, you can use this row index to find your zero indices in granule.

Using the below example:

granule = rand(1e3,2);
[K, V] = convhull(granule);
granule2 = zeros(size(granule));
tmp = granule(K,:)<0.1; %// all points on the convex hull but not in the granule
tmp2 = tmp(:,1)==tmp(:,2); %// both indices of your granule are below 0.1
granule2(tmp(tmp2)) = 2;

results in sum(tmp2)=11, thus there are 11 points in this case which are on the convex hull and have both indices below 0.1 (I could not use ==0 since there are no zeros in my array).

You might want to switch up the condition for tmp2 based on what you actually want for it.

Unsurprisingly more people have struggled with this and actually written code for this, see the MathWorks Central, code by John D'Errico.

Answer 2

This should work:

ix=find(granule);
[x,y,z]=ind2sub(size(granule),ix);
K=convhulln([x,y,z]);
% take the index of the points on the frontier
front=unique([K(:,1);K(:,2);K(:,3)]);

% calculate the minimum box containing the whole granule
minx=min(x(front));
maxx=max(x(front));
miny=min(y(front));
maxy=max(y(front));
minz=min(z(front));
maxz=max(z(front));
% Make a mesh for the points of this box
[X, Y, Z]=meshgrid(minx:maxx,miny:maxy,minz:maxz);
X=X(:);Y=Y(:);Z=Z(:);
% Remove the points of the granule from the box
ind=find(granule(sub2ind(size(granule),X,Y,Z))==0);
X=X(ind);Y=Y(ind);Z=Z(ind);

% For every face
for face=1:length(K),
    % take the coordinates of the vertices of the face
    P1=[x(K(face,1)), y(K(face,1)), z(K(face,1))];
    P2=[x(K(face,2)), y(K(face,2)), z(K(face,2))];
    P3=[x(K(face,3)), y(K(face,3)), z(K(face,3))];

    % calculate the normal to the face (croos product of two vectors on the face, 
    % we take the two sides of the face using the first point as pivot)
    normal=cross(P2-P1,P3-P1);
    % take another point on the convex hull
    inner=1;
    direction=0;
    while direction==0,
        P4=[x(inner), y(inner), z(inner)];

        % calculate the projection of the other point on the normal (always using P1 as pivot) 
        % to verify which is the sign for the projection of the points inside the convex hull
        direction=dot(P4-P1,normal);
        inner=inner+1;
    end
    % do the same for all the remaining points of the mesh
    project=sum(([X,Y,Z]-repmat(P1,size(X,1),1)).*repmat(normal,size(X,1),1),2);
    % if the sign of the projection is different to the one of the test point, 
    % the point is outside the convex hull so remove it
    ind=find(sign(project)*sign(direction)>=1);
    X=X(ind);Y=Y(ind);Z=Z(ind);
    plot3(X,Y,Z,'.')
end

What this code does is verifying, for every face, that every point in a rectangular mesh lays on the same side of this face of the convex hull.

To speed-up calculations the rectangular mesh is made only for the coordinates that touch the granule projection on the three coordinate planes, in this way, if the granule occupy the 5% of the area you can cut 90-95% of the points to check, without calculations.

Answer 3

What you need to do is use convhull to generate the watertight surface of your hull and then I see two options:

• For each point in your 1024^3 grid check if it's inside or outside the hull - maybe using a 'in polyherdron test' (such as this)

• Alternatively rather than testing on every voxel (3d pixel), you could potentially voxelise (ie rasterise) your convex hull - that is transforming this watertight surface into a set of voxel. An option for that is available here.

Once you know which points are within the convex hull you can easily remove the points from the granule.