Code covered by the BSD License

### Highlights from octree - partitioning 3D points into spatial subvolumes

4.90909
4.9 | 11 ratings Rate this file 53 Downloads (last 30 days) File Size: 16.3 KB File ID: #40732 Version: 1.2

# octree - partitioning 3D points into spatial subvolumes

by

### Sven (view profile)

11 Mar 2013 (Updated )

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

File Information
Description

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.

Acknowledgements
Required Products MATLAB
MATLAB release MATLAB 8.1 (R2013a)
05 Jul 2016 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

Comment only
01 Jul 2016 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

Comment only
24 Jun 2016 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

Comment only
24 Jun 2016 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!

14 Mar 2016 JIA Kanghao

### JIA Kanghao (view profile)

Thank you for the file.

17 Sep 2015 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

31 Oct 2014 Peter Morovic

### Peter Morovic (view profile)

Nice work, thanks for sharing!

30 Oct 2014 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.

Comment only
10 Oct 2014 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

04 Oct 2014 Zhouxin Xi

30 Sep 2014 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?

Comment only
29 Sep 2014 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) ?

Comment only
29 Sep 2014 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.

Comment only
05 Sep 2014 Sven

### Sven (view profile)

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

Comment only
05 Sep 2014 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.

Comment only
05 Sep 2014 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')

Comment only
04 Sep 2014 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

19 Dec 2013 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

08 Nov 2013 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

Comment only
17 Apr 2013 Siyi Deng

### Siyi Deng (view profile)

04 Apr 2013 Moein

### Moein (view profile)

@Sven Thank you very much.

Comment only
02 Apr 2013 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.

Comment only
01 Apr 2013 Moein

### Moein (view profile)

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

Comment only
29 Mar 2013 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.

Comment only
29 Mar 2013 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).

Comment only
28 Mar 2013 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?

Comment only
28 Mar 2013 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

13 Mar 2013 Anton Semechko

superb job