You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
How to draw a volxe size?
15 views (last 30 days)
Show older comments
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.
Accepted Answer
Walter Roberson
on 6 Sep 2021
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)
1 Comment
stefano chiappini
on 8 Sep 2021
Thank you for your reply, but this does not help me.
More Answers (2)
stefano chiappini
on 6 Sep 2021
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.
20 Comments
Walter Roberson
on 8 Sep 2021
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
on 8 Sep 2021
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
on 8 Sep 2021
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
on 8 Sep 2021
Edited: stefano chiappini
on 8 Sep 2021
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
on 8 Sep 2021
Which should be green and which should be red?
stefano chiappini
on 8 Sep 2021
It doesn't care, the important thing is that the result is similar to that shown in the figure.
Walter Roberson
on 8 Sep 2021
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
on 8 Sep 2021
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
on 9 Sep 2021
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
on 9 Sep 2021
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
on 9 Sep 2021
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
on 9 Sep 2021
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
on 9 Sep 2021
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
on 9 Sep 2021
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
on 9 Sep 2021
Edited: Walter Roberson
on 10 Sep 2021
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
on 9 Sep 2021
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
on 10 Sep 2021
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
on 10 Sep 2021
Thank you so much. You have answered my question.
But, where have you studied all information of Matlab?
Walter Roberson
on 10 Sep 2021
I have had a lot of practice answering questions ;-)
stefano chiappini
on 10 Sep 2021
I believe you. I will start with the documentation place on the Matlab website.
Thank you so much!!!!!!! ;- )
stefano chiappini
on 9 Sep 2021
Edited: Walter Roberson
on 9 Sep 2021
Have you idea how to plot this figure in 3D enviroment, plesae? Thank you so much
3 Comments
Walter Roberson
on 9 Sep 2021
No. There is no way to make your data look like that diagram.
stefano chiappini
on 9 Sep 2021
Edited: stefano chiappini
on 9 Sep 2021
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
on 9 Sep 2021
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.
See Also
Categories
Find more on Startup and Shutdown in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)