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)
How to draw a volxe size?
8 次查看（过去 30 天）
显示 更早的评论
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.
采纳的回答
Walter Roberson
2021-9-6
1 个评论
更多回答（2 个）
stefano chiappini
2021-9-6
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 个评论
Walter Roberson
2021-9-8
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
2021-9-8
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
2021-9-8
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
2021-9-8
编辑：stefano chiappini
2021-9-8
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
stefano chiappini
2021-9-8
It doesn't care, the important thing is that the result is similar to that shown in the figure.
Walter Roberson
2021-9-8
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
2021-9-8
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
2021-9-9
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
2021-9-9
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
2021-9-9
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
2021-9-9
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
2021-9-9
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
2021-9-9
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
2021-9-9
编辑：Walter Roberson
2021-9-10
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
2021-9-9
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
2021-9-10
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
2021-9-10
Thank you so much. You have answered my question.
But, where have you studied all information of Matlab?
stefano chiappini
2021-9-10
I believe you. I will start with the documentation place on the Matlab website.
Thank you so much!!!!!!! ;- )
stefano chiappini
2021-9-9
编辑：Walter Roberson
2021-9-9
Have you idea how to plot this figure in 3D enviroment, plesae? Thank you so much
3 个评论
stefano chiappini
2021-9-9
编辑：stefano chiappini
2021-9-9
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
2021-9-9
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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Startup and Shutdown 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改，无法完成操作。请重新加载页面以查看其更新后的状态。
Translated by 
您也可以从以下列表中选择网站：
如何获得最佳网站性能
选择中国网站（中文或英文）以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)