How to calculate volume change of specific ROIs between two point clouds?

I need to calculate the volume change between two point clouds representing consecutive timesteps of an experiment. The workflow I have come up with is as follows:

1. import point clouds and decorrelation map
2. generate mesh of pcs
3. subtract pc1 from pc2 to generate z-change
4. create mask to isolate regions of correlation map where value < 0.7
5. apply mask to z-change to pick out ROIs
6. use area of ROI and z-change to calculate volume of each decorrelated region

So far I've managed steps 1 through 3, and have attempted 4 using some code I have from a prior segmentation project, but I'm not managing to get any further. I have included the code I have here, attached the workspace with pc data, and the correlation map is inserted below the code. Eventually this needs to be put into a loop to process sets of ca. 1200 timesteps, but for now I'm just trying to get it to work for one.

% Generate mesh

x=linspace(min([Xpre;Xpost]),max([Xpre;Xpost]),1000);
y=linspace(min([Ypre;Ypost]),max([Ypre;Ypost]),1000);
[X,Y]=meshgrid(x,y);

Fpre=griddata(Xpre,Ypre,Zpre,X,Y);
Fpost=griddata(Xpost,Ypost,Zpost,X,Y);

% Point cloud subtraction and visualisation

figure()
colormap('jet')
mesh(X,Y,zeros(size(X)),Fpost-Fpre)
view([0 90])
colorbar
exportgraphics(figure,'t0775.png')
changeMap = imread('t0775.png');

% Read correlation map and generate mask

corrMap = imread('B0775.png');
mask = im2gray(corrMap);
BW = bwlabel(mask);
label_bg = BW(1,1);
uniquelabels = sort(unique(BW(:)));

% % This is where I get stuck and think I've gone wrong

nb_pixel_in_label = [];

for i1 = 1 : length(uniquelabels)
    if uniquelabels(i1) ~= label_bg
        nb_pixel_in_label(i1) = numel(find(BW==uniquelabels(i1)));
    else
        nb_pixel_in_label(i1) = 0;
    end
end

idx_roi = find(max(nb_pixel_in_label) == nb_pixel_in_label);

new_mask = zeros(size(mask));

for i1 = 1 : length(uniquelabels)
    if uniquelabels(i1) == uniquelabels(idx_roi)
        new_mask(find(BW==uniquelabels(i1))) = 1;
    else
        new_mask(find(BW==uniquelabels(i1))) = 0;
    end
end

zoneChange = changeMap .* repmat(new_mask,[1,1,3]);

Correlation map (img B0775.png)

Answer 1

[Ayush]

Hi Amy,

I understand that you are working on workflow to calculate the volume change between two-point clouds using a decorrelation map and a mesh grid representation. I will explain the potential issues and guide you to resolve the issue.

Problem in the given approach:

The issue seems to be associated with extracting the relevant regions from the correlation map, applying a mask and calculating the volume change.

The logic in your code for creating “new_mask” is likely problematic because:

• You are checking for the maximum number of pixels in each region but then assigning a mask to only that region. This might not be selecting the correct regions.
• The mask generation logic could be incorrect as it might not effectively isolate the regions with a correlation value less than 0.7.

Here is another approach to refine your code:

1. Instead of focusing on the largest regions in your correlation map, you can consider using a simple threshold to isolate regions where correlation is less than 0.7. This will help in creating a mask that excludes well-correlated areas.

corrMap = imread('B0775.png'); 
corrMapGray = im2gray(corrMap);   
% Apply a threshold to isolate regions with correlation < 0.7 
threshold = 0.7; 
mask = corrMapGray < threshold;  % Create binary mask where correlation is less than 0.7 
imshow(mask); 

2. Now that you have the mask, you can apply it directly to change map (“changeMap”). Here is a pseudo code for your reference:

% Ensure changeMap is in the same size and format as mask 
changeMap = imread('t0775.png');   
changeMapGray = rgb2gray(changeMap);   
% Apply the mask to extract only the regions of interest (ROIs) 
maskedChangeMap = changeMapGray .* double(mask);  % Element-wise multiplication to mask the change map 
figure; 
imshow(maskedChangeMap, []); 
colorbar; 

3. After you have successfully applied the mask, the remaining task is to calculate the volume. To do this, you can multiply the “z-change” (which corresponds to the |Fpost – Fpre| difference) by the area of each pixel (based on the grid resolution) to compute the volume change in each region.

Here is the pseudo code for your reference:

volumeChange = sum(maskedChangeMap(:));  % This gives the sum of the z-changes in the ROI 

% If your grid has spacing dx and dy, you need to multiply by dx*dy for accurate volume 
dx = x(2) - x(1);   
dy = y(2) - y(1);   
pixelArea = dx * dy;   

totalVolumeChange = volumeChange * pixelArea; 
disp(['Total Volume Change: ', num2str(totalVolumeChange)]); 

4. Once you have this working for one timestep, you can easily loop over the complete set of timesteps. Assuming the point cloud and correlation maps are stored in arrays or file names that can be iterated over, you can implement a loop to process the 1200 timesteps.

Here is the pseudo code that should perform the entire process from reading the point clouds and correlation map, generating the mesh, calculating the point cloud difference, creating, applying the mask, and calculating the volume change.

% variables present in given ‘mat’ file 
% 1. Generate mesh grid 
x = linspace(min([Xpre; Xpost]), max([Xpre; Xpost]), 1000); 
y = linspace(min([Ypre; Ypost]), max([Ypre; Ypost]), 1000); 
[X, Y] = meshgrid(x, y); 

% 2. Interpolate the point clouds onto the mesh grid 
Fpre = griddata(Xpre, Ypre, Zpre, X, Y); 
Fpost = griddata(Xpost, Ypost, Zpost, X, Y); 

% 3. Calculate the z-change (difference between Fpost and Fpre) 
changeMap = Fpost - Fpre; 

% 4. Visualize the change map 
figure; 
colormap('jet'); 
mesh(X, Y, zeros(size(X)), changeMap); 
view([0, 90]); 
colorbar; 
title('Z-Change (Fpost - Fpre)'); 

% 5. Read the correlation map (replace with your actual image) 
corrMap = imread('B0775.png');  
corrMapGray = im2gray(corrMap); 

% 6. Create mask by thresholding correlation map (value < 0.7) 
threshold = 0.7; 
mask = corrMapGray < threshold;   

% 7. Apply mask to change map 
maskedChangeMap = changeMap .* mask; 

% 8. Visualize the masked change map 
figure; 
imshow(maskedChangeMap, []); 
colorbar; 
title('Masked Change Map'); 

% 9. Calculate volume change 
dx = x(2) - x(1);   
dy = y(2) - y(1);   
pixelArea = dx * dy;  

% Sum the z-changes in the masked region to get the volume change 
volumeChange = sum(maskedChangeMap(:));   

% Calculate total volume change by multiplying with pixel area 
totalVolumeChange = volumeChange * pixelArea; 

% 10. Display the total volume change 
disp(['Total Volume Change: ', num2str(totalVolumeChange)]); 

Hope it helps!

[Amy M]

Hi Ayush,

This is great, thank you so much for your help!

I'm having an issue with applying the mask to the changeMap. I've attached an image that shows the changeMap visualisation, the mask, and the error when I try to run this line:

maskedChangeMap = changeMap .* mask;

I'm not sure if maybe cropping the corrMap would be helpful... essentially the mask region in the circle corresponds to the teal and yellow region at the western edge of the changeMap.