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:
- 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!
