# Java implementation of well-separated pair decomposition

I was trying to implement the algorithm of the 3D well-separated pair decomposition WSPD using an octree.

First, I begin by implementing the class OctreeNode as follows:

public class OctreeNode {

public static final int BSW = 0;
public static final int BSE = 1;
public static final int BNW = 2;
public static final int BNE = 3;
public static final int ASW = 4;
public static final int ASE = 5;
public static final int ANW = 6;
public static final int ANE = 7;

public int level = 0; // level set to 0.
public OctreeNode[] children = null;
public OctreeNode father = null; // father is null
public Point_3 p = null; //point stored in a leaf
public Point_3 center = null; // point to define the center of the box of the node.
public double height = 0.0, width = 0.0, depth = 0.0;


Every non-leaf Node has exactly 8 children. I implement a constructor for the class which takes a list of points in 3D.

/**
* Create the octree for storing an input point cloud.
*/
public OctreeNode(List<Point_3> points) {
// We construct the octree only for non empty point cloud.
assert(points.size() > 0);
/**
* We start by computing a box containing the point cloud.
*/
// The minimum value of x,y,z in the point cloud.
double minX = points.get(0).getX().doubleValue();
double minY = points.get(0).getY().doubleValue();
double minZ = points.get(0).getZ().doubleValue();

// The maximum value of x,y,z in the point cloud.
double maxX = points.get(0).getX().doubleValue();
double maxY = points.get(0).getY().doubleValue();
double maxZ = points.get(0).getZ().doubleValue();

for(Point_3 point : points) {

// update the minimum.
if(minX > point.getX().doubleValue()) {
minX = point.getX().doubleValue();
}

// update the minimum.
if(minY > point.getY().doubleValue()) {
minY = point.getY().doubleValue();
}

// update the minimum.
if(minZ > point.getZ().doubleValue()) {
minZ = point.getZ().doubleValue();
}

// update the maximum.
if(maxX < point.getX().doubleValue()) {
maxX = point.getX().doubleValue();
}

// update the maximum.
if(maxY < point.getY().doubleValue()) {
maxY = point.getY().doubleValue();
}

// update the maximum.
if(maxZ < point.getZ().doubleValue()) {
maxZ = point.getZ().doubleValue();
}
}

this.center = new Point_3((minX + maxX) / 2, (minY + maxY) / 2, (minZ + maxZ) / 2);
this.width = 0.75*(maxX - minX);
this.height = 0.75*(maxY - minY);
this.depth = 0.75*(maxZ - minZ);

for(Point_3 point : points) {
}

}


After, I implement the add function of my class:

/**
* @return true if the current node is a leaf.
*/
public boolean isLeaf() {
return (this.children == null);
}
/**
* @return true if the current node is empty.
*/
public boolean isEmpty() {
return (this.children == null && this.p == null);
}

/**
* @return true if the current node is single point.
*/
public boolean isSinglePoint() {
return (this.children == null && this.p != null);
}

/**
* @param o an Object.
* @return true if the current node is equals to o
*/
@Override
public boolean equals(Object o) {
OctreeNode n = (OctreeNode) o;
return (n != null) && (this.center.equals(n.center));
}

/**
* @return the diameter of the node : 0 if is empty or single point and the diameter of the box otherwise.
*/
public double diam() {
if(this.isLeaf()) {
return 0.0;
}
return 2*Math.sqrt(this.width*this.width
+ this.height*this.height
+ this.depth*this.depth);
}
/**
* Check if the point is in the boundary of the OctreeNode
* @param p a Point_3.
* @return true if p is in the cube of the OctreeNode
*/
public boolean inBoundary(Point_3 p) {
Vector_3 v = (Vector_3) this.center.minus(p);
return (Math.abs(v.getX().doubleValue()) <= this.width && Math.abs(v.getY().doubleValue()) <= this.height && Math.abs(v.getZ().doubleValue()) <= this.depth);
}
/**
* Add a node into the OctreeNode
*/

assert(this.center != null);
if(!this.inBoundary(p))
return;
// Check if the current node is a leaf.
if(this.isLeaf()) {
// Check if the current node is empty.
if(this.isEmpty()) {
this.p = p;
return;
}
else {
// Check if the node contains the same point already.
if(this.p.equals(p)) {
return;
}
// The current node contains only one point and the new point
// is different from the point of the current OctreeNode.
// We create the set of children for the current OctreeNode.
this.children = new OctreeNode[8];

// Initialize the children.
for(int i = 0; i < 8; i++) {
this.children[i] = new OctreeNode();
}

// For all children we put the current OctreeNode as father and
// we increment the level.
for(OctreeNode child : this.children) {
child.father = this;
child.level = this.level + 1;
}

// We compute then the center points for every child

Vector_3 vi = new Vector_3(this.width / 2, 0.0, 0.0);
Vector_3 vj = new Vector_3(0.0, this.height / 2, 0.0);
Vector_3 vk = new Vector_3(0.0, 0.0, this.depth / 2);

this.children[BSW].center = this.center.sum(vk.opposite()).sum(vj.opposite()).sum(vi.opposite());
this.children[BSE].center = this.center.sum(vk.opposite()).sum(vj.opposite()).sum(vi);
this.children[BNW].center = this.center.sum(vk.opposite()).sum(vj).sum(vi.opposite());
this.children[BNE].center = this.center.sum(vk.opposite()).sum(vj).sum(vi);
this.children[ASW].center = this.center.sum(vk).sum(vj.opposite()).sum(vi.opposite());
this.children[ASE].center = this.center.sum(vk).sum(vj.opposite()).sum(vi);
this.children[ANW].center = this.center.sum(vk).sum(vj).sum(vi.opposite());
this.children[ANE].center = this.center.sum(vk).sum(vj).sum(vi);

// We put half of the dimensions of the cube for every child.
for(OctreeNode child : children) {
child.width = this.width / 2;
child.depth = this.depth / 2;
child.height = this.height / 2;
}

// Look for the child to add the point of the current OctreeNode.
for(OctreeNode child : children) {
if(child.inBoundary(this.p)) {
break;
}
}

// Look for the child to add the point p.
for(OctreeNode child : children) {
if(child.inBoundary(p)) {
break;
}
}
// this.p = null;
}
}
else {

// Look for the child to add the point p.
for(OctreeNode child : children) {
if(child.inBoundary(p)) {
break;
}
}
}
}


I also implement a function distance to compute the distance between two nodes in the Octree which is the distance between the two boxes of the nodes:

/**
* @param o an OctreeNode.
* @param u an OctreeNode.
* @return the distance between the two nodes.
*/
public static double distance(OctreeNode o, OctreeNode u) {
if(o == null || o.isEmpty() || u == null || u.isEmpty()) {
return 0.0;
}
if(u.isLeaf() && o.isLeaf()) {
return u.p.distanceFrom(o.p).doubleValue();
}
double x = 0.0;
double y = 0.0;
double z = 0.0;
Vector_3 v;
if(u.isLeaf()) {
v = (Vector_3) u.p.minus(o.center);
x = Math.min(Math.abs(v.getX().doubleValue()) - o.width, 0.0);
y = Math.min(Math.abs(v.getX().doubleValue()) - o.height, 0.0);
z = Math.min(Math.abs(v.getX().doubleValue()) - o.depth, 0.0);
return Math.sqrt(x*x + y*y + z*z);
}
if(o.isLeaf()) {
return distance(u, o);
}
v = (Vector_3) u.center.minus(o.center);
x = Math.min(Math.abs(v.getX().doubleValue()) - o.width - u.width, 0.0);
y = Math.min(Math.abs(v.getX().doubleValue()) - o.height - u.height, 0.0);
z = Math.min(Math.abs(v.getX().doubleValue()) - o.depth - u.depth, 0.0);
return Math.sqrt(x*x + y*y + z*z);
}
}


Now, I implement the class Octree representing and Octree

public class Octree {
public OctreeNode root;

/**
* Create an octree from the array points.
*/
public Octree(Point_3[] points){
/**
* We use the constructor provided by the class OctreeNode to
* construct the root.
*/
this.root = new OctreeNode(Arrays.asList(points));
}
}


Here I have a question. Do you think my implementation is the best way to do it?

After this part I implement the well-separated pair decomposition class:

public class WellSeparatedPairDecomposition {
class Pair<X,Y> {
X x; Y y;
Pair(X xx, Y yy) {
x = xx;
y = yy;
}
public X getFirst() {
return x;
}
public Y getSecond() {
return y;
}

@Override
public boolean equals(Object o) {
Pair<X,Y> p = (Pair<X,Y>) o;
return (p!=null && p.getFirst() == x && p.getSecond() == y);
}
}
/**
* Compute the WSPD from an octree and a threshold s.
* @param T the Octree,
* @param s threshold.
* @return the pairs of WSPD.
*/
public OctreeNode[][] wspd (Octree T, double epsilon) {
Set<Pair<OctreeNode, OctreeNode>> H = new HashSet<Pair<OctreeNode, OctreeNode>>();

wspd(T.root, T.root, epsilon, H);

int n = H.size();
OctreeNode[][] result = new OctreeNode[n][2];
int i = 0;
for(Pair<OctreeNode, OctreeNode> pair : H) {
result[i][0] = pair.getFirst();
result[i][1] = pair.getSecond();
i++;
}
return result;
}

boolean isWS(OctreeNode u, OctreeNode v, double epsilon) {
if(u == null || v == null || u.isEmpty() || v.isEmpty()) {
return false;
}
double distance = OctreeNode.distance(u, v);
return (u.diam() < epsilon*distance && v.diam() < epsilon*distance);
}

void wspd(OctreeNode u, OctreeNode v, double epsilon, Set<Pair<OctreeNode, OctreeNode>> H) {
if(u.isEmpty() || v.isEmpty() || (v.isLeaf() && u.isLeaf() && v.equals(u)))
return;
if(isWS(u, v, epsilon)) {
return;
}
if(u.level > v.level) {
OctreeNode temp = u;
u = v;
v = temp;
}
if(!u.isLeaf()) {
for(OctreeNode uchild : u.children) {
wspd(uchild, v, epsilon, H);
}
}
}
}


My implementation works without errors. However, when I tested, it is very slow and for big size tests it gives an outOfMemoryError. How can I improve this implementation?

• The lack of indentation makes parts of the code hard to follow. – 1201ProgramAlarm Dec 21 '18 at 21:23

Overall the code is okay but can definitely be improved.

Some things I would however change.

Instead of writing short variable names and then explaining them in comments, a better approach is to have "self-documenting code". What this means is to make your variables descriptive enough to convey to the reader what they represent.

An example is in your method

/**
* @param o an OctreeNode.
* @param u an OctreeNode.
* @return the distance between the two nodes.
*/
public static double distance(OctreeNode o, OctreeNode u) {


Here you can simply use firstNode and secondNode since we are returning the distance between them.

Here is another example:

/**
* Compute the WSPD from an octree and a threshold s.
* @param T the Octree,
* @param s threshold.
* @return the pairs of WSPD.
*/
public OctreeNode[][] wspd (Octree T, double epsilon) {


Here the T would be better called a tree, or the threshold would be better called exactly that, the threshold. Also, I don't see the parameter you mentioned passed into the function (possibly a bug)?

Pair(X xx, Y yy) {
x = xx;
y = yy;
}


Here I would change xx and yy, as they could be better represented with a different variable name.

Over here, what are these nodes?

OctreeNode u, OctreeNode v


It would be preferred to describe them, so that it's easier to see what's going on when you use these variables throughout your code.

@param p a Point_3.


Why is a Point_3 here is the only object which has an underscore? You want to pick a naming style and consistently use it. Inconsistency in code leads to possible subtle bugs down the road - along with difficulty for an outside reader to follow it.

You should make an attempt to vectorize much of your code. It's an important way to re-think your code in terms that allow for computers to efficiently execute SIMD (single-instruction multiple-data) and achieve massive speedup as compared to naive implementations. There are many libraries to allow you to vectorize your code. So far the most promising-looking one (caveat: I haven't tried it) seems to be JBLAS.

Specifically, this:

double minX = points.get(0).getX().doubleValue();
double minY = points.get(0).getY().doubleValue();
double minZ = points.get(0).getZ().doubleValue();


should be represented by one vector, and this:

double maxX = points.get(0).getX().doubleValue();
double maxY = points.get(0).getY().doubleValue();
double maxZ = points.get(0).getZ().doubleValue();


should be represented by one vector. This loop:

for(Point_3 point : points) {

// update the minimum.
// ...


should not exist at all, and you should be calling a vectorized routine instead.

Similar to the above, this:

this.width = 0.75*(maxX - minX);
this.height = 0.75*(maxY - minY);
this.depth = 0.75*(maxZ - minZ);


should not be represented by three separate members, nor should you be repeating the same operation three times. Instead, represent it as one vector perhaps called this.size.

This:

public double diam() {
if(this.isLeaf()) {
return 0.0;
}
return 2*Math.sqrt(this.width*this.width +
this.height*this.height +
this.depth*this.depth);
}


is you rolling your own Euclidean norm, which you shouldn't do. BLAS has functions for this. Similar vectorization should be done elsewhere in your code, such as inBoundary. This:

for(OctreeNode child : children) {
child.width = this.width / 2;
child.depth = this.depth / 2;
child.height = this.height / 2;
}


should treat every child's size as a single vector rather than separate components. You should also compute this.size / 2 outside of the loop. It may even be possible to represent children as a matrix with each row representing a child and the columns representing x, y and z respectively, which would eliminate the need for a loop altogether.

Along the same lines, distance can be greatly tightened up by dealing with vectors rather than individual components, particularly lines like these:

x = Math.min(Math.abs(v.getX().doubleValue()) - o.width - u.width, 0.0);
y = Math.min(Math.abs(v.getX().doubleValue()) - o.height - u.height, 0.0);
z = Math.min(Math.abs(v.getX().doubleValue()) - o.depth - u.depth, 0.0);
return Math.sqrt(x*x + y*y + z*z);


This can be done in one line of vectorized code.