How to compare two 3D volumes using dice metric?

Hello everyone,

I am trying to compare two segmentations regarding their overlap.

Each Segmentation is defined by N 3D coordinates (2 examples attached)

1) I want to interpolate a volume for every segmentation.

2) I want to create a binary mask for each segmentation (everything inside the volume = 1, outside = 0).

3) I want to calculate the dice metric between volume A and volume B.

So far I have:

% Extract points for the current segmentation

points = ROIs(19).points; %ROIs contains n=30 segmentation,

% Remove duplicate points

unique_points = unique(points, 'rows');

% Calculate the volume of the convex hull

[k1, volume] = convhull(unique_points(:,1), unique_points(:,2), unique_points(:,3));

volumes(i) = volume;

% visualize segmentation volume

figure(1)

trisurf(k1,unique_points(:,1), unique_points(:,2), unique_points(:,3))

How would you create a 3D binary mask and compare two volumes using dice metric.

Many thanks in advance!

Max

Answer 1

[Jaynik]

Hi Max,

To create a 3D binary mask and compare two volumes using the Dice metric, you have already made progress by extracting points and calculating the volume of the convex hull. The next step is to create a binary mask for each volume. You can use the 'inpolyhedron' function available on File Exchange to obtain the binary mask: https://www.mathworks.com/matlabcentral/fileexchange/37856-inpolyhedron-are-points-inside-a-triangulated-volume

Based on the provided code and the sample files, you can refer the following code to calculate the "Sørensen–Dice" coefficient:

fileID_1 = fopen('Seg1.txt', 'r');
fileID_2 = fopen('Seg2.txt', 'r');
formatSpec = '%f %f %f';
data_1 = textscan(fileID_1, formatSpec);
data_2 = textscan(fileID_2, formatSpec);
fclose(fileID_1);
fclose(fileID_2);

unique_points_1 = [data_1{1}, data_1{2}, data_1{3}];
unique_points_2 = [data_2{1}, data_2{2}, data_2{3}];

[k1, volume1] = convhull(unique_points_1(:,1), unique_points_1(:,2), unique_points_1(:,3));
[k2, volume2] = convhull(unique_points_2(:,1), unique_points_2(:,2), unique_points_2(:,3));

% Plotting the convex hulls
figure;
trisurf(k1, unique_points_1(:,1), unique_points_1(:,2), unique_points_1(:,3), 'FaceColor', 'red', 'FaceAlpha', 0.5);
hold on;
trisurf(k2, unique_points_2(:,1), unique_points_2(:,2), unique_points_2(:,3), 'FaceColor', 'blue', 'FaceAlpha', 0.5);
hold off;

all_points = [unique_points_1; unique_points_2];

% Define the 3D grid bounds (you need to define gridX, gridY, gridZ)
% For example, you can use the min and max of the combined points
gridX = linspace(min(all_points(:,1)), max(all_points(:,1)), 100);
gridY = linspace(min(all_points(:,2)), max(all_points(:,2)), 100);
gridZ = linspace(min(all_points(:,3)), max(all_points(:,3)), 100);

binary_mask_1 = inpolyhedron(k1, unique_points_1, gridX, gridY, gridZ);
binary_mask_2 = inpolyhedron(k2, unique_points_2, gridX, gridY, gridZ);

intersection = binary_mask_1 & binary_mask_2;
dice_metric = 2 * nnz(intersection) / (nnz(binary_mask_1) + nnz(binary_mask_2));

disp(['Dice metric: ', num2str(dice_metric)]);

The grid resolution ("gridX", "gridY", "gridZ") is set to 100 in the sample code can be adjusted based on the input data. A higher resolution will provide a more detailed mask but will require more computation.

Hope this helps!

[Maximilian Fenski]

Hi Jaynik,

many thanks for your answer.

Works like a charm.

In the meantime i came up with a different way: I calculate the interpolated and the intersected volume directly using alphashape and inshape functions. I tested it and the DSC is quit close.

% Create volume 1 and volume 2, calculate volume 1, leave no hole: Alpha=1, plot volumes

%Volume 1
    shp1 = alphaShape(data1);
    shp1.Alpha = 1;
    % Volume 2
    shp2 = alphaShape(data2)
    shp2.Alpha = 1
    % Find shp1's and shp2's common points ("id")
    x1 = data1(:,1);
    x2 = data2(:,1);
    y1 = data1(:,2);
    y2 = data2(:,2);
    z1 = data1(:,3);
    z2 = data2(:,3);
    id1=inShape(shp2,data1);
    id2=inShape(shp1,data2);
    % create intersection volume 
    shp3=alphaShape([x1(id1); x2(id2)], [y1(id1); y2(id2)], [z1(id1); z2(id2)]);
    shp3.Alpha = 1
    V3 = volume(shp3);
    % Calculate Dice score metric 

    DSC = (2*V3)/(V1+V2)
   
    % to visualize the volumes and intersection
    figure(1) 
    plot(shp1)
    hold on 
    plot(shp2) 
    hold on 
    plot(shp3) 
    hold off 

[Jaynik]

Great to know that it worked for you.