Ytdyklly

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?

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.