How to draw a volxe size?

Hi, i am looking for script to calculate the crown volume of tree point cloud by the voxel size. The file input is .txt type.

I thank you in advance.

Answer 1

[Walter Roberson]

T = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/730384/test%20chioma.txt');
scatter3(T(:,1), T(:,2), T(:,3));

[~, crown_volume] = boundary(T)

crown_volume = 42.2327

[stefano chiappini]

Thank you for your reply, but this does not help me.

Answer 2

[stefano chiappini]

Hi, thank you for your answer.

I want draw a plot like in the following figure allowed in this answer. I will set up size voxel about 0.3 meter.

I hope that you will could help me.

Thank you so much.

[Walter Roberson]

I do not know what the third column of input represents. If it represents a z coordinate like what I plotted, then how do you want to calculate the values that are represented by color in your output ?

Are you wanting to calculate data density? https://www.mathworks.com/matlabcentral/fileexchange/31726-data-density-plot Are you wanting to calculate 2D histogram for a particular height (in order to do that, we need to know a thickness of the slice.)

[stefano chiappini]

My main aim is to switch by pointcloud to voxelization. I read a lot of paper but i didn't found script on Matlab. Have you know how i can this transformation?

[Walter Roberson]

vox = 0.3;
cmap = colormap(parula);

T = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/730384/test%20chioma.txt');
x = T(:,1); y = T(:,2); z = T(:,3);
xidx = floor(x/vox); xbin = xidx - min(xidx) + 1;
yidx = floor(y/vox); ybin = yidx - min(yidx) + 1;
zidx = floor(z/vox); zbin = zidx - min(zidx) + 1;
xq = xidx * vox; yq = yidx * vox; zq = zidx * vox;
scatter3(xq, yq, zq)

counts = accumarray([ybin(:), xbin(:), zbin(:)], 1);

xv = (min(xidx):max(xidx)) * vox;
yv = (min(yidx):max(yidx)) * vox;
zv = (min(zidx):max(zidx)) * vox;

LV = floor(min(zv)):max(zv);
for K = 1 : length(LV)
  L = LV(K);
  [~, zidx] = min(abs(zv - L));
  Z = counts(:,:,zidx);
  Z(Z == 0) = nan;
  figure();
  imagesc(xv, yv, Z); colormap(cmap);
  xlabel('x'); ylabel('y'); zlabel('z');
  title("Z = " + LV(K));
end

[stefano chiappini]

Thank you for you optimal answer. Your tip is very useful.

But i want draw the plot like the fiigure. I would like with the differente voxel size: 0.3 and 1.0, and set the colour voxel with green or red. Do you have idea how setting this customization? and how to calculate volume with different voxel size?

Thank you so much

[Walter Roberson]

Which should be green and which should be red?

[stefano chiappini]

It doesn't care, the important thing is that the result is similar to that shown in the figure.

[Walter Roberson]

No, I have no ideas on how to produce those diagrams. I could make guesses, but I do not know what they represent. I especially have no idea what the red one represents.

If a, b, c are intended to represent the same tree from https://www.mathworks.com/matlabcentral/answers/1447449-how-to-draw-a-volxe-size#answer_781809 then that tree has a span of about 3 metres in height. If you count the number of vertical bands in c, there are 9, suggesting a voxelization of 0.3. But then when we look at b, with its greater subdivisions, the natural inference would be that b is intended to show voxelization with 0.3 and c is intended to show voxelization with 1.0 .

[stefano chiappini]

The difference between the green and red trees is in the size of the voxels. The figure was just to show the example of the plot I wanted. I made a plot with the voxelSurf.m command but was unable to customise the figure.

Also for the volume calculation, taking into account the number and size of voxels, how can I do it?

thanks

[Walter Roberson]

counts = accumarray([ybin(:), xbin(:), zbin(:)], 1);

nnz(counts) is the number of voxels that have something in them. counts at any location tells you how many hits there were, but if you are looking at "volume" from the point of view of whether a box sized vox * vox * vox could be put into the space without hitting any branch, then nnz(counts) will tell you the number of obstructed areas without telling you anything about how much obstruction is present there. You would then multiply nnz() by vox^3 to get a volume estimate.

But I am not sure that is the meaning of "volume" for this purpose. Another plausible meaning of "volume" would be to ask, if you were to "shrink-wrap" a surface over the points, then what would be the enclosed space? For that, the volume output of boundary() times vox^3 would give you the volume estimate.

[Walter Roberson]

It is not at all clear in diagram c what the differences in color mean. It looks sort of like it is depth cued with lighting coming from the right hand side casting shadows. So the below diagrams are effectively depth maps to the closest surface. "Front occupancy" is brightest where the closest surface is nearer to you; "Back" occupancy is darkest where the furthest surface is furthest from you. So when you look at the front, it is brightnest near the centre, reflecting the fact that it is closest to you near the centre; when you look at the back, it is darkest near the centre, reflecting the fact that it is furthest from you near the centre.

T = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/730384/test%20chioma.txt');
x = T(:,1); y = T(:,2); z = T(:,3);

voxlist = [1, 0.3];
for voxidx = 1 : length(voxlist)
    vox = voxlist(voxidx);
    
    xidx = floor(x/vox); xbin = xidx - min(xidx) + 1;
    yidx = floor(y/vox); ybin = yidx - min(yidx) + 1;
    zidx = floor(z/vox); zbin = zidx - min(zidx) + 1;
    
    xb = [min(xidx), max(xidx)] * vox; 
    yb = [min(yidx), max(yidx)] * vox;
    zb = [min(zidx), max(zidx)] * vox;
    
    counts = accumarray([ybin(:), xbin(:), zbin(:)], 1);
    occupied = counts > 0;
    
    oidx = occupied .* (1:size(occupied,2));
    oidx(oidx == 0) = nan;
    froidx = permute(max(oidx,[],2), [1 3 2]);
    fralpha = double(~isnan(froidx));
    baoidx = permute(min(oidx,[],2), [1 3 2]);
    baalpha = double(~isnan(baoidx));
    maxidx = max(baoidx(:));
    
    c = linspace(0, 1, maxidx).';
    cmap = zeros(maxidx,3);
    cmap(:,voxidx) = c;
    
    figure
    imagesc(xb, zb, froidx./maxidx, 'alphadata', fralpha); colormap(cmap);
    xlabel('x'); ylabel('height'); title("front occupancy vox = " + vox);
    
    figure
    imagesc(xb, zb, baoidx./maxidx, 'alphadata', baalpha); colormap(cmap)
    xlabel('x'); ylabel('height'); title("back occupancy vox = " + vox)
end

[stefano chiappini]

Hi, thank you so much for your answer. I tthank you again for your valuable tips for this issue.

I have tried to run the following code, but the volume is not correct. The score with your script is 9.35 mc, but the exact value should be 68.33 mc (the last score i obtained by VoxR package in R programming)

T = testchioma;
x = T(:,1); y = T(:,2); z = T(:,3);
voxlist = [1, 0.05];
for voxidx = 1 : length(voxlist)
vox = voxlist(voxidx);
xidx = floor(x/vox); xbin = xidx - min(xidx) + 1;
yidx = floor(y/vox); ybin = yidx - min(yidx) + 1;
zidx = floor(z/vox); zbin = zidx - min(zidx) + 1;
xb = [min(xidx), max(xidx)] * vox;
yb = [min(yidx), max(yidx)] * vox;
zb = [min(zidx), max(zidx)] * vox;
counts = accumarray([ybin(:), xbin(:), zbin(:)], 1);
occupied = counts > 0;
oidx = occupied .* (1:size(occupied,2));
oidx(oidx == 0) = nan;
froidx = permute(max(oidx,[],2), [1 3 2]);
fralpha = double(~isnan(froidx));
baoidx = permute(min(oidx,[],2), [1 3 2]);
baalpha = double(~isnan(baoidx));
maxidx = max(baoidx(:));
c = linspace(0, 1, maxidx).';
cmap = zeros(maxidx,3);
cmap(:,voxidx) = c;
figure
imagesc(xb, zb, froidx./maxidx, 'alphadata', fralpha); colormap(cmap);
xlabel('x'); ylabel('height'); title("front occupancy vox = " + vox);
figure
imagesc(xb, zb, baoidx./maxidx, 'alphadata', baalpha); colormap(cmap)
xlabel('x'); ylabel('height'); title("back occupancy vox = " + vox)
end

d=nnz(counts);
volume=d*(vox^3);

Where do i wrong with my script?

[stefano chiappini]

An output request by my project is like the following figure.

The point cloud is trasformed to voxelization with voxe size 0.05 m. The window figure is 3D and editing by user. Here i have setted all cube of green colour and grid axis X,Y and Z.

My superadvisor request me of process this data in Matlab language.

[Walter Roberson]

T = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/730384/test%20chioma.txt');
x = T(:,1); y = T(:,2); z = T(:,3);

voxlist = [1, 0.3];
for voxidx = 1 : length(voxlist)
    vox = voxlist(voxidx);
    
    xidx = floor(x/vox); xbin = xidx - min(xidx) + 1;
    yidx = floor(y/vox); ybin = yidx - min(yidx) + 1;
    zidx = floor(z/vox); zbin = zidx - min(zidx) + 1;
    
    xb = [min(xidx), max(xidx)] * vox; 
    yb = [min(yidx), max(yidx)] * vox;
    zb = [min(zidx), max(zidx)] * vox;
    
    counts = accumarray([ybin(:), xbin(:), zbin(:)], 1);
    voxcount = nnz(counts)
    volestimate = voxcount * vox^3
    
    occupied = counts > 0;
    
    oidx = occupied .* (1:size(occupied,2));
    oidx(oidx == 0) = nan;
    froidx = permute(max(oidx,[],2), [1 3 2]);
    fralpha = double(~isnan(froidx));
    baoidx = permute(min(oidx,[],2), [1 3 2]);
    baalpha = double(~isnan(baoidx));
    maxidx = max(baoidx(:));
    
    c = linspace(0, 1, maxidx).';
    cmap = zeros(maxidx,3);
    cmap(:,voxidx) = c;
    
    figure
    imagesc(xb, zb, froidx./maxidx, 'alphadata', fralpha); colormap(cmap);
    xlabel('x'); ylabel('height'); title("front occupancy vox = " + vox);
    
    figure
    imagesc(xb, zb, baoidx./maxidx, 'alphadata', baalpha); colormap(cmap)
    xlabel('x'); ylabel('height'); title("back occupancy vox = " + vox)
end

voxcount = 127

volestimate = 127

voxcount = 2101

volestimate = 56.7270

[Walter Roberson]

You would get a slightly different estimate if you were to start your voxels at different places, such as min(x) or max(x) and use relative positions. The estimate I do works with absolute positions -- notice the Heights are absolute units, not relative to the bottom-most location in the point-cloud.

[Walter Roberson]

format long g
T = readmatrix('https://www.mathworks.com/matlabcentral/answers/uploaded_files/730384/test%20chioma.txt');
x = T(:,1); y = T(:,2); z = T(:,3);

voxlist = [1, 0.3];
for voxidx = 1 : length(voxlist)
    vox = voxlist(voxidx)
    
    xidx = floor(x/vox); xbin = xidx - min(xidx) + 1;
    yidx = floor(y/vox); ybin = yidx - min(yidx) + 1;
    zidx = floor(z/vox); zbin = zidx - min(zidx) + 1;
    xq = xidx * vox; yq = yidx * vox; zq = zidx * vox;
    
    counts = accumarray([ybin(:), xbin(:), zbin(:)], 1);
    occupied = counts > 0;
    oind = find(occupied);
    
    vol_estimate = length(oind) * vox.^3
    
    [ox, oy, oz] = ind2sub(size(occupied), oind); 
    oxq = (ox + min(xidx) - 1) * vox;
    oyq = (oy + min(yidx) - 1) * vox;
    ozq = (oz + min(zidx) - 1) * vox;
    noq = length(oxq);
    
    rF = [1 2 3 4 1; 8 7 6 5 8; 1 4 6 7 1; 2 8 5 3 2; 1 7 8 2 1; 3 5 6 4 3];
    rV = [0 0 0; 1 0 0; 1 0 1; 0 0 1; 1 1 1; 0 1 1; 0 1 0; 1 1 0] * vox;
    
    %noq = 2;
    poxq = oxq(1:noq); poyq = oyq(1:noq); pozq = ozq(1:noq);
    
    allF = repmat(rF, noq, 1) + 8*repelem((0:noq-1).', size(rF,1), 1);
    allV = repelem([poxq, poyq, pozq], size(rV,1), 1) + repmat(rV, noq, 1);
    
    cmap = zeros(1,3);
    cmap(:,voxidx) = 1;
    figure
    
    patch('Faces', allF, 'Vertices', allV, 'FaceColor', cmap, 'LineWidth', 0.1)
    title("vox = " + vox)
    view(3)
end

vox = 
     1

vol_estimate = 
   127

vox = 
                       0.3

vol_estimate = 
                    56.727

[stefano chiappini]

Great, you are fantastic. And last question: how do i calculate the numbers of voxel plotting in figure?

thank you so much for you strong help.

[Walter Roberson]

I updated to output the volume estimate.

By the way, the reason that vox 1 is done before vox 0.3, is that it made it easier to figure out which color to use for the drawing. The first vox processed writes 1.0 into column 1 of the colormap, producing red. The second vox processed writes 1.0 into the column 2 of the colormap, producing green. If you had another vox then it would write into column 3, producing blue; if you went to more colors than that, you would need a minor revision to be able to specify the colors.

[stefano chiappini]

Thank you so much. You have answered my question.

But, where have you studied all information of Matlab?

[Walter Roberson]

I have had a lot of practice answering questions ;-)

[stefano chiappini]

I believe you. I will start with the documentation place on the Matlab website.

Thank you so much!!!!!!! ;- )

Answer 3

[stefano chiappini]

https://it.mathworks.com/matlabcentral/answers/1447449-how-to-draw-a-volxe-size#comment_1728079

Have you idea how to plot this figure in 3D enviroment, plesae? Thank you so much

[Walter Roberson]

No. There is no way to make your data look like that diagram.

[stefano chiappini]

i work with abosolute units.

Infact i heave set early my script in the following way

data=fopen('test chioma.txt');
x = data(:, 1);
y = data(:, 2);
z = data(:, 3);
% Subtract means
x = x - mean(x);
y = y - mean(y);
z = z - mean(z);
%--------------------------------------------------------------------------------------------------------
% Display the data.
subplot(1, 2, 1);
plot3(x, y, z, '.', 'MarkerSize', 3);
grid on;
xlabel('Column 1');
ylabel('Column 2');
zlabel('Column 3');
title('Original data not classified yet.');
hFig = gcf;
hFig.WindowState = 'maximized'; % May not work in earlier versions of MATLAB.
drawnow;

[Walter Roberson]

The result of fopen() is a file identifier, which is a scalar integer. You then try to take the first three columns of that scalar integer.