**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

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

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

Walter Roberson
on 9 Sep 2021

Edited: Walter Roberson
on 10 Sep 2021

format long g

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

### See Also

### Categories

### 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 (한국어)