File Exchange

image thumbnail

octree - partitioning 3D points into spatial subvolumes

version 1.2 (16.3 KB) by

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

4.93333
16 Ratings

41 Downloads

Updated

View License

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.

Comments and Ratings (36)

Tank you, it was useful for my work.

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?

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

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

Very well done.

jia dzh

thank you very much

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

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

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

Thank you for the file.

elmer magsino

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

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

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

Zhouxin Xi

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?

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

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

Siyi Deng (view profile)

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

Anton Semechko (view profile)

superb job

Updates

1.2

Thanks to Jing for a BinDepths bugfix

1.1

Added ability to shrink bins to minimally cover their child points/bins

MATLAB Release
MATLAB 8.1 (R2013a)

Download apps, toolboxes, and other File Exchange content using Add-On Explorer in MATLAB.

» Watch video