File Exchange

## octree - partitioning 3D points into spatial subvolumes

version 1.2.0.0 (16.3 KB) by Sven

### Sven (view profile)

OcTree recursively splits a large set of points into smaller subvolumes. A QuadTree but in 3D.

Updated 05 Sep 2014

OcTree point decomposition in 3D
OcTree is used to create a tree data structure of bins containing 3D
points. Each bin may be recursively decomposed into 8 child bins.
http://en.wikipedia.org/wiki/Octree
OT = OcTree(PTS) creates an OcTree from an N-by-3 matrix of point
coordinates.

OT = OcTree(...,'PropertyName',VALUE,...) takes any of the following
property values:

binCapacity - Maximum number of points a bin may contain. If more
points exist, the bin will be recursively subdivided.
Defaults to ceil(numPts/10).
maxDepth - Maximum number of times a bin may be subdivided.
Defaults to INF.
maxSize - Maximum size of a bin edge. If any dimension of a bin
exceeds maxSize, it will be recursively subdivided.
Defaults to INF.
minSize - Minimum size of a bin edge. Subdivision will stop after
any dimension of a bin gets smaller than minSize.
Defaults to 1000*eps.
style - Either 'equal' (default) or 'weighted'. 'equal'
subdivision splits bins at their central coordinate
(ie, one bin subdivides into 8 equally sized bins).
'weighted' subdivision divides bins based on the mean
of all points they contain. Weighted subdivision is
slightly slower than equal subdivision for a large
number of points, but it can produce a more efficient
decomposition with fewer subdivisions.

Example 1: Decompose 200 random points into bins of 20 points or less,
then display each bin with its points in a separate colour.
pts = (rand(200,3)-0.5).^2;
OT = OcTree(pts,'binCapacity',20);
figure
boxH = OT.plot;
cols = lines(OT.BinCount);
doplot3 = @(p,varargin)plot3(p(:,1),p(:,2),p(:,3),varargin{:});
for i = 1:OT.BinCount
set(boxH(i),'Color',cols(i,:),'LineWidth', 1+OT.BinDepths(i))
doplot3(pts(OT.PointBins==i,:),'.','Color',cols(i,:))
end
axis image, view(3)

Example 2: Decompose 200 random points into bins of 10 points or less,
shrunk to minimallly encompass their points, then display.
pts = rand(200,3);
OT = OcTree(pts,'binCapacity',10,'style','weighted');
OT.shrink
figure
boxH = OT.plot;
cols = lines(OT.BinCount);
doplot3 = @(p,varargin)plot3(p(:,1),p(:,2),p(:,3),varargin{:});
for i = 1:OT.BinCount
set(boxH(i),'Color',cols(i,:),'LineWidth', 1+OT.BinDepths(i))
doplot3(pts(OT.PointBins==i,:),'.','Color',cols(i,:))
end
axis image, view(3)

OcTree methods:
shrink - Shrink each bin to tightly encompass its children
query - Ask which bins a new set of points belong to.
plot, plot3 - Plots bin bounding boxes to the current axes.

OcTree properties:
Points - The coordinate of points in the decomposition.
PointBins - Indices of the bin that each point belongs to.
BinCount - Total number of bins created.
BinBoundaries - BinCount-by-6 [MIN MAX] coordinates of bin edges.
BinDepths - The # of subdivisions to reach each bin.
BinParents - Indices of the bin that each bin belongs to.
Properties - Name/Val pairs used for creation (see help above)

Please post comments to this FEX page for this entry if you find any bugs or have any feature requests. There's plenty of room for improved efficiency and enhancement.

### Cite As

Sven (2019). octree - partitioning 3D points into spatial subvolumes (https://www.mathworks.com/matlabcentral/fileexchange/40732-octree-partitioning-3d-points-into-spatial-subvolumes), MATLAB Central File Exchange. Retrieved .

Ehsan Osloub

Jacob Wirtz

Lou

### Lou (view profile)

Can someone point me in the right direction: I would like to safe a matrix for each bin which contains all points in this bin. How would this work? I assume I can somehow use the BinCount?

Alisa Brown

### Alisa Brown (view profile)

So I organize my points into an octree and then I query a new point, but it's outside the bin boundaries, what will happen to the point? Would I get an error or would it approximate it to the closest bin?

Emily Schnorr

### Emily Schnorr (view profile)

What is the computation involved in a query? - For a point cloud with 10^5 points and a maximum depth of 5 I get rather slow computation times (~0.01 seconds). Should't this be quite fast?

Shiming Wei

Louis Lac

ali gurbuz

Francis

### Francis (view profile)

Thank you. I have a question how to use this code solving 2D problems.

ha ha

### ha ha (view profile)

please check the code for "max size" .
I tried : OT = OcTree(pts, 'maxSize',4);
%It mean that "the size each voxel edge < 4". But it run wrong.

amir keramatian

Karoly Soltesz

### Karoly Soltesz (view profile)

Tank you, it was useful for my work.

Jonathan Sjölund

Chuyen Nguyen

### Chuyen Nguyen (view profile)

Dear Sven,

Is it possible to add "parfor" in your algorithm to make it faster? Thank you very much

Julio

### Julio (view profile)

I can transform the cubes in a .stl file?

Floriane Zongo

### Floriane Zongo (view profile)

Dear Sven,

thank you so much for your code. I would like to ask if there is a way to make your code run from 3d points decomposition to 3d imagies decomposition. If not, can you redirect me to another code please?

Thank you again Sven

Yiyu Hong

### Yiyu Hong (view profile)

Hi, Sven
I want to use this octree struture to find intersections between a ray and a mesh, do you have any ideas? "Consistent Mesh Partitioning and Skeletonisation using the Shape Diameter Function" This paper said using octree structure can accelerate to calculate ray-mesh intersections.
Thank you!!!!

Thomas Czerniawski

Very well done.

jia dzh

### jia dzh (view profile)

thank you very much

Chuyen Nguyen

### Chuyen Nguyen (view profile)

Hello Sven, I try to save the "OT" into workspace. It gives me a error, it cannot save the variables of "OT". Is there any ways to save these variables?

Thank you Sven

Chuyen Nguyen

### Chuyen Nguyen (view profile)

Hello Sven, I used your software to voxelizer our terrestrial lidar data. I want to use these voxels for computation purpose. Can we access individual voxel?
Thank you very much Sven

Sven

### Sven (view profile)

@Liye Sun, the first bin will always be the bounding box of all your points. If your scene is 3x10x15 then this first bin will never be a cube. You can either:

- Add extra points OUTSIDE your scene (so that the first bin is 15x15x15. Every divided bin will then be cubic.

- Manually initialise a whole array of 1x1x1 bins inside your 3x10x15 space. You can do this by first making an OcTree object with maxDepth=0 (so that you get a single 3x10x15 bin), then calling divideBin() in a loop so that you split the OcTree 3 times, 10 times, and 15 times. Then you can set maxDepth=inf and call divide() to let the OcTree continue to make smaller bins

Liye SUN

### Liye SUN (view profile)

Dear Sven,

I am wondering how to obtain cubic bins when building octree for non-cubic world? The fact is when I run this code to build octree for a 3m*10m*15m scene, the bins are all cuboids.

Thank you!

JIA Kanghao

### JIA Kanghao (view profile)

Thank you for the file.

elmer magsino

### elmer magsino (view profile)

Hello guys,

Thank you very much for this file. This made me understand octree through visualization. Great job Sven.

From the octree, how can I form the 3D map? Thank you in advance.

Elmer

Peter Morovic

### Peter Morovic (view profile)

Nice work, thanks for sharing!

Sven

### Sven (view profile)

Hi Jonathan,

I haven't implemented specifically a nearest neighbour routine. You could use OcTree to find which bin your query location belongs to, and then find the nearest neighbour within that bin. This will be faster than querying all points in a large point cloud, however this alone won't guarantee correctness (since a point in an adjacent bin may be closer) unless the distance to the found nearest neighbour is less than the distance from your query point to any edges of its bin.

Did that help you out? If your point cloud isn't massive (or even if it is), John D'Errico's ipdm is pretty efficient, and there are some nearest neighbour kdtree implementations that may be helpful too.

Jonathan

### Jonathan (view profile)

Hello Sven,
I would like to use your Octree code in order to find the N closest point of an arbitrary position in my domain. I was using the quedtree code on the Matlab File Exchange but now that I am going to 3D, I would like to use your code for this for a fast localisation of neighbour points of a specific coordinates. Is it something already possible with your code? I found the OT.query function but this is not exactly what I want to do. Any idea of I can make it from your code? Thank you very much anyway for your great code. Best regards,Jonathan

jigsaw

Sven

### Sven (view profile)

Hi Christophe, for your second question, you can force the decomposition to go to a specific level with the optional parameters:
OT = OcTree(pts,'binCapacity',inf,'maxDepth',3,'maxSize',0);

The OT will still contain the higher levels above that, but you can easily get only the desired level bins with a mask such as:
OT.BinDepths==2

For the neighbours that's an interesting question. It will be easy to find the 3 "siblings" of a given bin (they will be the same level and have the same parent bin No. as that bin). It might need a bit more work to find all neighbours... what exactly would you define as a neighbour? Can only "leaf-node" bins be neighbours or do you include any level's bin that shares a face with your target a neighbour?

Christophe Langrenne

### Christophe Langrenne (view profile)

Hi,
An other question : Is it possible to have only one level tree (for exemple, only bins of level 3) ?

Christophe Langrenne

### Christophe Langrenne (view profile)

Hi Sven,
I have to find the neighbors of a given box (for a Fast Multipole Method application). Do you have an easy way to do this ?

Thanks and "bravo" for your great job.

Sven

### Sven (view profile)

@Jing: Yep, thanks for finding that bug. Fixed now.

Sven

### Sven (view profile)

@Wu Jun: It's strange but QuadTree decomposition in images is a different beast to QuadTree decomposition of 2d points, but they both have the same name. I'm afraid it's the same situation for OcTree decomposition (of image volumes vs 3d points). OcTree.m is targeted at points rather than image volume data.

Sven

### Sven (view profile)

@Greg: (sorry for the late reply) This is best done *after* plotting the bins. You can just find which bins don't contain a point and then delete or deactivate that bin's plot handle:

% Run example 1 and then:
binHasPts = ismember(1:OT.BinCount, OT.PointBins)
set(boxH(~binHasPts),'Visible','off')

Jing

### Jing (view profile)

A bug in the code???
should "this.BinParents(binNo)+1" be "this.BinDepths(binNo)+1" ??

% Prevent dividing beyond the maximum depth

if this.BinParents(binNo)+1 >= this.Properties.maxDepth
continue;
end

Wu Jun

### Wu Jun (view profile)

Hi Sven, Thanks for your great job. I want to use Octree method to divides the 3D images volume into little 3D cube just like Quadree method (matlab qtdecomp()function). Do you have any plan to add this feature to your code ? Regards

Greg

### Greg (view profile)

Hi Sven,

Thank you for your great job.
I am an absolute beginner in Matlab but I managed to run your code and got quite good results. I am wondering if it is possible to delete the empty bins (i.e. bins of the higher levels)and keep only the child bins of the lowest level before plotting the final OcTree.

Kind regards
Greg

Siyi Deng

Moein

### Moein (view profile)

@Sven Thank you very much.

Sven

### Sven (view profile)

@Moein: Thanks for your feedback and request, I'm looking at implementing the concept of element connectivity so that you can partition based on collections of points (elements) as well as the points themselves. Note that it seems you have an immediate need (2:1 balancing) for a particular application (FE meshing). I would suggest that if your problem is urgent then you might consider either adding to the class or inheriting and making a new class specific to your needs yourself. I will work on element connectivity first, and node balancing will be some time away, so my schedule won't necessarily fit in with yours.

Moein

### Moein (view profile)

Hi Sven
I need element connectivity in your code , How can I get element connectivity ? Thanks for your help

Moein

### Moein (view profile)

Also, in order to balance your mesh, you need information about elements connectivity. Because if you have second or higher level hanging nodes in one edge, you should subdivide all the elements which connected to that edge.you can find some algorithms for finding neighbors in octree mesh in Journals.

Moein

### Moein (view profile)

@Sven Thanks for your reply. In octree decomposition, there will arise “hanging” nodes, i.e., nodes that lie on the edge of an undivided box, rather than at one of its corners. If, during the subdivision, a second hanging node is placed on the edge of a box, that box should be itself subdivided. In this way, the final mesh does not have any “second-level” hanging nodes.
Also I think it would be better to define the binary mask to your code showing which elements of a set are to be subdivided. It gives more freedom to users. The suitable outputs for Finite element method are the positions of the all nodes , The connectivity array which shows the number of 8 points of each cubic element.(for octree case).

Sven

### Sven (view profile)

@Moein: Thanks for the feedback. Yes, it doesn't yet consider 2:1 balancing. I have a few plans for new features, starting with including the concept of elements (ie, triangles or quads or hexas... hopefully I can do N-dim elements).
I hadn't yet thought through node balancing, do you have a specific example you can provide and the type of output you'd want? For example are you simply looking for adjusting a current decomposition (swapping neigbours) based on particular criteria? If so, is that criteria based solely on points or do you need information about "element connectivity" too?

Moein

### Moein (view profile)

Hi , Thanks for your great job. I tested this code and I noticed that the code doesn't consider 2:1 rule for neighboring mesh elements. ( I mean the final tree is not balanced). Is it true? If yes, Do you have any plan to add this feature to your code ? (Its very important in FEM mesh generation). Regards

Anton Semechko

superb job