# Sort a cloud of points with respect to an unstructured mesh of hexahedral cells

### Question

How would you sort a cloud of points with respect to an unstructured mesh of hexahedral cells?

Each cell has a centre and a unique label to represent it. There are two cloud points basically (original point cloud, and a point cloud of the cell centres), but the cell geometry information (bounding box) may be of use, I am not sure.

### Results

I've done some asking around, and searching through the literature:

if the mesh is hexahedral and unstructured, the problem is reduced to an orthogonal range search. For this purpose, k-d trees are most often used. If the mesh is refined based on an octree data structure, the range search algorithm can be built around it. The goal is to avoid dealing with the direct mesh geometry and concentrate on point cloud A - point cloud relationship B. Point cloud A: query points, point cloud B: mesh cell centres.

• Can you clarify what you mean when saying "sort with respect to (any kind of) mesh"? Are you looking for a binning algorithm, (how many points are the in each cell)? – Szabolcs Feb 29 '12 at 10:06
• I do not understand your question quite clearly, what is the goal of sorting the points? Like making the mesh more regular? – Shuhao Cao Feb 29 '12 at 16:39
• There is a separate point cloud scattered accross the unstructured volume mesh. I need to communicate data from the cell centres to the point cloud and vice versa. – tmaric Feb 29 '12 at 19:02
• @tomislav-maric: Could you please write your solution as an answer, and then accept your own answer? This procedure is generally the accepted practice for answering your own question effectively, rather than adding the "[SOLVED]" tag to the question; also, it will earn you more reputation, because people can upvote your answer. – Geoff Oxberry Mar 7 '12 at 15:02

Important note: This answer does not answer the actual question, but it was left undeleted per request. Embarrassingly I confused hexahedral and hexagonal. The question is about sorting points into arbitrary hexahedral cells in 3D while this solution sorts points into regular hexagonal cells in 2D, or irregular ones that correspond to some Voronoi tesselation in any dimension. This method is applicable only if the mesh was generated as a Voronoi tesselation in the first place (which does seem to be an occasionally used approach).

I am not sure what you mean by sort here, but I assume you want to sort the point into hexagonal bins on the plane.

Mathematica is what I know, so I am going to show you how to do it in Mathematica, but the method can be ported to other systems. The idea is that a hexagonal lattice is the dual of a triangular one: it can be generated as the Voronoi diagram of a points in triangular arrangement. A point from the cloud belongs to a given hexagon if it is closer to the centre of that hexagon than to the centre of any other hexagon.

This method will work for meshes of different shapes too, as long as they can be generates as the Voronoi diagram of some point arrangement. (E.g. the hexagons don't need to be regular.)

Let's generate the mesh. This is a triangle lattice:

pts = Join @@ Table[{x, Sqrt[3] y}, {x, 0, 4}, {y, 0, 2}];

points = Join[pts, TranslationTransform[{1/2, Sqrt[3]/2}] /@ pts];

Needs["ComputationalGeometry"]
PlanarGraphPlot[points, LabelPoints -> False]


Its dual is the hexagonal one we are interested in:

DiagramPlot[points, LabelPoints -> False]


This builds a function nf which finds the index of the hexagon centre that some cloud point is closest to. It is the key to the method:

nf = Nearest[N[points] -> Range@Length[points]];


Now let's generate a cloud of 1000 random points and sort them with nf:

cloud = RandomReal[{-1/2, 5}, {1000, 2}];

indices = First /@ nf /@ cloud;


indices contains the indices of the centres that each cloud point is closest to. This is the information we needed. Now we can make a histogram out of them ...

Histogram[indices]


... or colour each of them ...

Show[
DiagramPlot[points, LabelPoints -> False],
PlotRange -> All, AspectRatio -> Automatic
]


... or do any sort of fancy visualization we want.

tally = Tally[indices];

ListDensityPlot[Join[points, List /@ Sort[tally][[All, 2]], 2],
InterpolationOrder -> 0,
Epilog -> (Text[#2, points[[#1]]] & @@@ tally),
PlotRange -> {{-.5, 5}, {-.5, 5}}, Mesh -> All,
ColorFunction -> (ColorData["BeachColors"][1 - #] &)]


The key point here was the function that finds the closest point to something (Nearest`). Mathematica has this built in, but there's a chance your system doesn't. If this is the case, please see this question on how to efficiently implement such a function (or just go with the naive linear time implementation if you don't have a huge amount of points to process).

• Thanks a lot! Basically what i need is a relationship that shows a connection between each point and a "bin" as you called it (3 dimensional hexahedral box). What you suggest seems very interesting, but I am dealing with meshes of millions of boxes and hundreds of thousands of points potentially.. The question is what costs more: dual mesh creation or working with bounding boxes of the "bins" and using a k-d tree for searching. I'm very new to this topic, so I really don't want to go in the wrong direction. – tmaric Feb 29 '12 at 10:54
• @tomislav-maric I am sorry I misunderstood your question. Unfortunately I am also new to this topic, and I wasn't even familiar with $k$-d trees a few days ago (as you can see from the SO question I linked to). Also, I am not sure that your mesh may correspond to a Voronoi tessellation at all (probably not), so my method may be useless here. I think it's best if I delete this answer as it is misleading. – Szabolcs Feb 29 '12 at 10:58
• Don't definitely delete it, someone may find it useful! :) It may turn pout to be the solution to the problem, it's just that I cannot yet accept it until I read about it. – tmaric Feb 29 '12 at 11:01
• And thanks for such a detailed answer, if I could I would give you more points! :) – tmaric Feb 29 '12 at 11:05
• @tomislav-maric Looking at the votes, I am worried that my answer will decrease the chance that you will get a useful one, or will contribute to the misunderstanding. I think it's more productive if I delete. – Szabolcs Feb 29 '12 at 15:11