3
\$\begingroup\$

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) {
        this.add(point);
    }

}

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
 */

public void add(Point_3 p) {
    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)) {
                    child.add(this.p);
                    break;
                }
            }

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

        // Look for the child to add the point p.
        for(OctreeNode child : children) {
            if(child.inBoundary(p)) {
                child.add(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)) {
        H.add(new Pair<OctreeNode, OctreeNode>(u, v));
        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?

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

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.

\$\endgroup\$
1
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.