You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by lu...@apache.org on 2014/12/02 22:20:27 UTC

[1/8] [math] Added splitters in Boundary attributes.

Repository: commons-math
Updated Branches:
  refs/heads/master e11c00085 -> d21d3a94a


Added splitters in Boundary attributes.


Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/6525cc40
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/6525cc40
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/6525cc40

Branch: refs/heads/master
Commit: 6525cc4035fe8df3308c03f0ca9f450728e122ac
Parents: 26ee481
Author: Luc Maisonobe <lu...@apache.org>
Authored: Sun Nov 30 18:09:06 2014 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 2 15:24:31 2014 +0100

----------------------------------------------------------------------
 src/changes/changes.xml                         |  4 ++
 .../geometry/partitioning/AbstractRegion.java   | 67 +++++++++++++-----
 .../partitioning/AbstractSubHyperplane.java     | 67 +++++++++++++-----
 .../partitioning/BoundaryAttribute.java         | 32 +++++++++
 .../geometry/partitioning/BoundaryBuilder.java  | 23 ++++--
 .../geometry/partitioning/Characterization.java | 67 ++++++++++++++----
 .../math3/geometry/partitioning/NodesSet.java   | 72 +++++++++++++++++++
 .../geometry/partitioning/RegionFactory.java    | 74 ++++++++++++++++----
 8 files changed, 338 insertions(+), 68 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/6525cc40/src/changes/changes.xml
----------------------------------------------------------------------
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index c8f0afd..46ad704 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -73,6 +73,10 @@ Users are encouraged to upgrade to this version as this release not
   2. A few methods in the FastMath class are in fact slower that their
   counterpart in either Math or StrictMath (cf. MATH-740 and MATH-901).
 ">
+      <action dev="luc" type="add" >
+        Boundary attributes in regions now provides the BSP tree nodes that
+        were used to split the sub-hyperplane froming the boundary part of the facet. 
+      </action>
       <action dev="luc" type="fix" issue="MATH-1162" >
         Fixed a problem with vanishing cut sub-hyperplanes during BSP tree merging.
       </action>

http://git-wip-us.apache.org/repos/asf/commons-math/blob/6525cc40/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java
index 3fca476..fd3dcc7 100644
--- a/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java
@@ -19,7 +19,9 @@ package org.apache.commons.math3.geometry.partitioning;
 import java.util.ArrayList;
 import java.util.Collection;
 import java.util.Comparator;
+import java.util.HashMap;
 import java.util.Iterator;
+import java.util.Map;
 import java.util.TreeSet;
 
 import org.apache.commons.math3.geometry.Point;
@@ -349,7 +351,7 @@ public abstract class AbstractRegion<S extends Space, T extends Space> implement
     /** {@inheritDoc} */
     public BSPTree<S> getTree(final boolean includeBoundaryAttributes) {
         if (includeBoundaryAttributes && (tree.getCut() != null) && (tree.getAttribute() == null)) {
-            // we need to compute the boundary attributes
+            // compute the boundary attributes
             tree.visit(new BoundaryBuilder<S>());
         }
         return tree;
@@ -465,36 +467,65 @@ public abstract class AbstractRegion<S extends Space, T extends Space> implement
      * transform to the instance
      */
     public AbstractRegion<S, T> applyTransform(final Transform<S, T> transform) {
-        return buildNew(recurseTransform(getTree(false), transform));
+
+        // transform the tree, except for boundary attribute splitters
+        final Map<BSPTree<S>, BSPTree<S>> map = new HashMap<BSPTree<S>, BSPTree<S>>();
+        final BSPTree<S> transformedTree = recurseTransform(getTree(false), transform, map);
+
+        // set up the boundary attributes splitters
+        for (final Map.Entry<BSPTree<S>, BSPTree<S>> entry : map.entrySet()) {
+            if (entry.getKey().getCut() != null) {
+                @SuppressWarnings("unchecked")
+                BoundaryAttribute<S> original = (BoundaryAttribute<S>) entry.getKey().getAttribute();
+                if (original != null) {
+                    @SuppressWarnings("unchecked")
+                    BoundaryAttribute<S> transformed = (BoundaryAttribute<S>) entry.getValue().getAttribute();
+                    for (final BSPTree<S> splitter : original.getSplitters()) {
+                        transformed.getSplitters().add(map.get(splitter));
+                    }
+                }
+            }
+        }
+
+        return buildNew(transformedTree);
+
     }
 
     /** Recursively transform an inside/outside BSP-tree.
      * @param node current BSP tree node
      * @param transform transform to apply
+     * @param map transformed nodes map
      * @return a new tree
      */
     @SuppressWarnings("unchecked")
-    private BSPTree<S> recurseTransform(final BSPTree<S> node, final Transform<S, T> transform) {
+    private BSPTree<S> recurseTransform(final BSPTree<S> node, final Transform<S, T> transform,
+                                        final Map<BSPTree<S>, BSPTree<S>> map) {
 
+        final BSPTree<S> transformedNode;
         if (node.getCut() == null) {
-            return new BSPTree<S>(node.getAttribute());
-        }
+            transformedNode = new BSPTree<S>(node.getAttribute());
+        } else {
+
+            final SubHyperplane<S>  sub = node.getCut();
+            final SubHyperplane<S> tSub = ((AbstractSubHyperplane<S, T>) sub).applyTransform(transform);
+            BoundaryAttribute<S> attribute = (BoundaryAttribute<S>) node.getAttribute();
+            if (attribute != null) {
+                final SubHyperplane<S> tPO = (attribute.getPlusOutside() == null) ?
+                    null : ((AbstractSubHyperplane<S, T>) attribute.getPlusOutside()).applyTransform(transform);
+                final SubHyperplane<S> tPI = (attribute.getPlusInside()  == null) ?
+                    null  : ((AbstractSubHyperplane<S, T>) attribute.getPlusInside()).applyTransform(transform);
+                // we start with an empty list of splitters, it will be filled in out of recursion
+                attribute = new BoundaryAttribute<S>(tPO, tPI, new NodesSet<S>());
+            }
 
-        final SubHyperplane<S>  sub = node.getCut();
-        final SubHyperplane<S> tSub = ((AbstractSubHyperplane<S, T>) sub).applyTransform(transform);
-        BoundaryAttribute<S> attribute = (BoundaryAttribute<S>) node.getAttribute();
-        if (attribute != null) {
-            final SubHyperplane<S> tPO = (attribute.getPlusOutside() == null) ?
-                null : ((AbstractSubHyperplane<S, T>) attribute.getPlusOutside()).applyTransform(transform);
-            final SubHyperplane<S> tPI = (attribute.getPlusInside()  == null) ?
-                null  : ((AbstractSubHyperplane<S, T>) attribute.getPlusInside()).applyTransform(transform);
-            attribute = new BoundaryAttribute<S>(tPO, tPI);
+            transformedNode = new BSPTree<S>(tSub,
+                                             recurseTransform(node.getPlus(),  transform, map),
+                                             recurseTransform(node.getMinus(), transform, map),
+                                             attribute);
         }
 
-        return new BSPTree<S>(tSub,
-                                    recurseTransform(node.getPlus(),  transform),
-                                    recurseTransform(node.getMinus(), transform),
-                                    attribute);
+        map.put(node, transformedNode);
+        return transformedNode;
 
     }
 

http://git-wip-us.apache.org/repos/asf/commons-math/blob/6525cc40/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractSubHyperplane.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractSubHyperplane.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractSubHyperplane.java
index 8a4931f..3fd6b54 100644
--- a/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractSubHyperplane.java
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractSubHyperplane.java
@@ -16,6 +16,9 @@
  */
 package org.apache.commons.math3.geometry.partitioning;
 
+import java.util.HashMap;
+import java.util.Map;
+
 import org.apache.commons.math3.geometry.Space;
 
 /** This class implements the dimension-independent parts of {@link SubHyperplane}.
@@ -107,39 +110,67 @@ public abstract class AbstractSubHyperplane<S extends Space, T extends Space>
      */
     public AbstractSubHyperplane<S, T> applyTransform(final Transform<S, T> transform) {
         final Hyperplane<S> tHyperplane = transform.apply(hyperplane);
+
+        // transform the tree, except for boundary attribute splitters
+        final Map<BSPTree<T>, BSPTree<T>> map = new HashMap<BSPTree<T>, BSPTree<T>>();
         final BSPTree<T> tTree =
-            recurseTransform(remainingRegion.getTree(false), tHyperplane, transform);
+            recurseTransform(remainingRegion.getTree(false), tHyperplane, transform, map);
+
+        // set up the boundary attributes splitters
+        for (final Map.Entry<BSPTree<T>, BSPTree<T>> entry : map.entrySet()) {
+            if (entry.getKey().getCut() != null) {
+                @SuppressWarnings("unchecked")
+                BoundaryAttribute<T> original = (BoundaryAttribute<T>) entry.getKey().getAttribute();
+                if (original != null) {
+                    @SuppressWarnings("unchecked")
+                    BoundaryAttribute<T> transformed = (BoundaryAttribute<T>) entry.getValue().getAttribute();
+                    for (final BSPTree<T> splitter : original.getSplitters()) {
+                        transformed.getSplitters().add(map.get(splitter));
+                    }
+                }
+            }
+        }
+
         return buildNew(tHyperplane, remainingRegion.buildNew(tTree));
+
     }
 
     /** Recursively transform a BSP-tree from a sub-hyperplane.
      * @param node current BSP tree node
      * @param transformed image of the instance hyperplane by the transform
      * @param transform transform to apply
+     * @param map transformed nodes map
      * @return a new tree
      */
     private BSPTree<T> recurseTransform(final BSPTree<T> node,
                                         final Hyperplane<S> transformed,
-                                        final Transform<S, T> transform) {
-        if (node.getCut() == null) {
-            return new BSPTree<T>(node.getAttribute());
-        }
+                                        final Transform<S, T> transform,
+                                        final Map<BSPTree<T>, BSPTree<T>> map) {
 
-        @SuppressWarnings("unchecked")
-        BoundaryAttribute<T> attribute =
-            (BoundaryAttribute<T>) node.getAttribute();
-        if (attribute != null) {
-            final SubHyperplane<T> tPO = (attribute.getPlusOutside() == null) ?
-                null : transform.apply(attribute.getPlusOutside(), hyperplane, transformed);
-            final SubHyperplane<T> tPI = (attribute.getPlusInside() == null) ?
-                null : transform.apply(attribute.getPlusInside(), hyperplane, transformed);
-            attribute = new BoundaryAttribute<T>(tPO, tPI);
+        final BSPTree<T> transformedNode;
+        if (node.getCut() == null) {
+            transformedNode = new BSPTree<T>(node.getAttribute());
+        } else {
+
+            @SuppressWarnings("unchecked")
+            BoundaryAttribute<T> attribute = (BoundaryAttribute<T>) node.getAttribute();
+            if (attribute != null) {
+                final SubHyperplane<T> tPO = (attribute.getPlusOutside() == null) ?
+                    null : transform.apply(attribute.getPlusOutside(), hyperplane, transformed);
+                final SubHyperplane<T> tPI = (attribute.getPlusInside() == null) ?
+                    null : transform.apply(attribute.getPlusInside(), hyperplane, transformed);
+                // we start with an empty list of splitters, it will be filled in out of recursion
+                attribute = new BoundaryAttribute<T>(tPO, tPI, new NodesSet<T>());
+            }
+
+            transformedNode = new BSPTree<T>(transform.apply(node.getCut(), hyperplane, transformed),
+                    recurseTransform(node.getPlus(),  transformed, transform, map),
+                    recurseTransform(node.getMinus(), transformed, transform, map),
+                    attribute);
         }
 
-        return new BSPTree<T>(transform.apply(node.getCut(), hyperplane, transformed),
-                              recurseTransform(node.getPlus(), transformed, transform),
-                              recurseTransform(node.getMinus(), transformed, transform),
-                              attribute);
+        map.put(node, transformedNode);
+        return transformedNode;
 
     }
 

http://git-wip-us.apache.org/repos/asf/commons-math/blob/6525cc40/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryAttribute.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryAttribute.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryAttribute.java
index ea7d9a7..dad884c 100644
--- a/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryAttribute.java
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryAttribute.java
@@ -45,6 +45,9 @@ public class BoundaryAttribute<S extends Space> {
      */
     private final SubHyperplane<S> plusInside;
 
+    /** Sub-hyperplanes that were used to split the boundary part. */
+    private final NodesSet<S> splitters;
+
     /** Simple constructor.
      * @param plusOutside part of the node cut sub-hyperplane that
      * belongs to the boundary and has the outside of the region on
@@ -52,11 +55,33 @@ public class BoundaryAttribute<S extends Space> {
      * @param plusInside part of the node cut sub-hyperplane that
      * belongs to the boundary and has the inside of the region on the
      * plus side of its underlying hyperplane (may be null)
+     * @deprecated as of 3.4, the constructor has been replaced by a new one
+     * which is not public anymore, as it is intended to be used only by
+     * {@link BoundaryBuilder}
      */
+    @Deprecated
     public BoundaryAttribute(final SubHyperplane<S> plusOutside,
                              final SubHyperplane<S> plusInside) {
+        this(plusOutside, plusInside, null);
+    }
+
+    /** Simple constructor.
+     * @param plusOutside part of the node cut sub-hyperplane that
+     * belongs to the boundary and has the outside of the region on
+     * the plus side of its underlying hyperplane (may be null)
+     * @param plusInside part of the node cut sub-hyperplane that
+     * belongs to the boundary and has the inside of the region on the
+     * plus side of its underlying hyperplane (may be null)
+     * @param splitters sub-hyperplanes that were used to
+     * split the boundary part (may be null)
+     * @since 3.4
+     */
+    BoundaryAttribute(final SubHyperplane<S> plusOutside,
+                      final SubHyperplane<S> plusInside,
+                      final NodesSet<S> splitters) {
         this.plusOutside = plusOutside;
         this.plusInside  = plusInside;
+        this.splitters   = splitters;
     }
 
     /** Get the part of the node cut sub-hyperplane that belongs to the
@@ -81,4 +106,11 @@ public class BoundaryAttribute<S extends Space> {
         return plusInside;
     }
 
+    /** Get the sub-hyperplanes that were used to split the boundary part.
+     * @return sub-hyperplanes that were used to split the boundary part
+     */
+    public NodesSet<S> getSplitters() {
+        return splitters;
+    }
+
 }

http://git-wip-us.apache.org/repos/asf/commons-math/blob/6525cc40/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java
index 80038f1..cea4de3 100644
--- a/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java
@@ -28,11 +28,6 @@ import org.apache.commons.math3.geometry.Space;
  */
 class BoundaryBuilder<S extends Space> implements BSPTreeVisitor<S> {
 
-    /** Simple constructor.
-     */
-    public BoundaryBuilder() {
-    }
-
     /** {@inheritDoc} */
     public Order visitOrder(BSPTree<S> node) {
         return Order.PLUS_MINUS_SUB;
@@ -43,6 +38,7 @@ class BoundaryBuilder<S extends Space> implements BSPTreeVisitor<S> {
 
         SubHyperplane<S> plusOutside = null;
         SubHyperplane<S> plusInside  = null;
+        NodesSet<S>      splitters   = null;
 
         // characterize the cut sub-hyperplane,
         // first with respect to the plus sub-tree
@@ -57,6 +53,9 @@ class BoundaryBuilder<S extends Space> implements BSPTreeVisitor<S> {
                 // this part belongs to the boundary,
                 // it has the outside on its plus side and the inside on its minus side
                 plusOutside = minusChar.insideTouching();
+                splitters = new NodesSet<S>();
+                splitters.addAll(minusChar.getInsideSplitters());
+                splitters.addAll(plusChar.getOutsideSplitters());
             }
         }
 
@@ -69,11 +68,23 @@ class BoundaryBuilder<S extends Space> implements BSPTreeVisitor<S> {
                 // this part belongs to the boundary,
                 // it has the inside on its plus side and the outside on its minus side
                 plusInside = minusChar.outsideTouching();
+                if (splitters == null) {
+                    splitters = new NodesSet<S>();
+                }
+                splitters.addAll(minusChar.getOutsideSplitters());
+                splitters.addAll(plusChar.getInsideSplitters());
+            }
+        }
+
+        if (splitters != null) {
+            // the parent nodes are natural splitters for boundary sub-hyperplanes
+            for (BSPTree<S> up = node.getParent(); up != null; up = up.getParent()) {
+                splitters.add(up);
             }
         }
 
         // set the boundary attribute at non-leaf nodes
-        node.setAttribute(new BoundaryAttribute<S>(plusOutside, plusInside));
+        node.setAttribute(new BoundaryAttribute<S>(plusOutside, plusInside, splitters));
 
     }
 

http://git-wip-us.apache.org/repos/asf/commons-math/blob/6525cc40/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.java
index a67f873..58efe0f 100644
--- a/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.java
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.java
@@ -16,6 +16,9 @@
  */
 package org.apache.commons.math3.geometry.partitioning;
 
+import java.util.ArrayList;
+import java.util.List;
+
 import org.apache.commons.math3.exception.MathInternalError;
 import org.apache.commons.math3.geometry.Space;
 
@@ -32,6 +35,12 @@ class Characterization<S extends Space> {
     /** Part of the cut sub-hyperplane that touch inside cells. */
     private SubHyperplane<S> insideTouching;
 
+    /** Nodes that were used to split the outside touching part. */
+    private final NodesSet<S> outsideSplitters;
+
+    /** Nodes that were used to split the outside touching part. */
+    private final NodesSet<S> insideSplitters;
+
     /** Simple constructor.
      * <p>Characterization consists in splitting the specified
      * sub-hyperplane into several parts lying in inside and outside
@@ -45,9 +54,11 @@ class Characterization<S extends Space> {
      * @param sub sub-hyperplane to characterize
      */
     public Characterization(final BSPTree<S> node, final SubHyperplane<S> sub) {
-        outsideTouching = null;
-        insideTouching  = null;
-        characterize(node, sub);
+        outsideTouching  = null;
+        insideTouching   = null;
+        outsideSplitters = new NodesSet<S>();
+        insideSplitters  = new NodesSet<S>();
+        characterize(node, sub, new ArrayList<BSPTree<S>>());
     }
 
     /** Filter the parts of an hyperplane belonging to the boundary.
@@ -61,29 +72,33 @@ class Characterization<S extends Space> {
      * outside/inside) do belong to the boundary.</p>
      * @param node current BSP tree node
      * @param sub sub-hyperplane to characterize
+     * @param splitters nodes that did split the current one
      */
-    private void characterize(final BSPTree<S> node, final SubHyperplane<S> sub) {
+    private void characterize(final BSPTree<S> node, final SubHyperplane<S> sub,
+                              final List<BSPTree<S>> splitters) {
         if (node.getCut() == null) {
             // we have reached a leaf node
             final boolean inside = (Boolean) node.getAttribute();
             if (inside) {
-                addInsideTouching(sub);
+                addInsideTouching(sub, splitters);
             } else {
-                addOutsideTouching(sub);
+                addOutsideTouching(sub, splitters);
             }
         } else {
             final Hyperplane<S> hyperplane = node.getCut().getHyperplane();
             switch (sub.side(hyperplane)) {
             case PLUS:
-                characterize(node.getPlus(), sub);
+                characterize(node.getPlus(),  sub, splitters);
                 break;
             case MINUS:
-                characterize(node.getMinus(), sub);
+                characterize(node.getMinus(), sub, splitters);
                 break;
             case BOTH:
                 final SubHyperplane.SplitSubHyperplane<S> split = sub.split(hyperplane);
-                characterize(node.getPlus(),  split.getPlus());
-                characterize(node.getMinus(), split.getMinus());
+                splitters.add(node);
+                characterize(node.getPlus(),  split.getPlus(),  splitters);
+                characterize(node.getMinus(), split.getMinus(), splitters);
+                splitters.remove(splitters.size() - 1);
                 break;
             default:
                 // this should not happen
@@ -94,24 +109,30 @@ class Characterization<S extends Space> {
 
     /** Add a part of the cut sub-hyperplane known to touch an outside cell.
      * @param sub part of the cut sub-hyperplane known to touch an outside cell
+     * @param splitters sub-hyperplanes that did split the current one
      */
-    private void addOutsideTouching(final SubHyperplane<S> sub) {
+    private void addOutsideTouching(final SubHyperplane<S> sub,
+                                    final List<BSPTree<S>> splitters) {
         if (outsideTouching == null) {
             outsideTouching = sub;
         } else {
             outsideTouching = outsideTouching.reunite(sub);
         }
+        outsideSplitters.addAll(splitters);
     }
 
     /** Add a part of the cut sub-hyperplane known to touch an inside cell.
      * @param sub part of the cut sub-hyperplane known to touch an inside cell
+     * @param splitters sub-hyperplanes that did split the current one
      */
-    private void addInsideTouching(final SubHyperplane<S> sub) {
+    private void addInsideTouching(final SubHyperplane<S> sub,
+                                   final List<BSPTree<S>> splitters) {
         if (insideTouching == null) {
             insideTouching = sub;
         } else {
             insideTouching = insideTouching.reunite(sub);
         }
+        insideSplitters.addAll(splitters);
     }
 
     /** Check if the cut sub-hyperplane touches outside cells.
@@ -129,6 +150,17 @@ class Characterization<S extends Space> {
         return outsideTouching;
     }
 
+    /** Get the nodes that were used to split the outside touching part.
+     * <p>
+     * Splitting nodes are internal nodes (i.e. they have a non-null
+     * cut sub-hyperplane).
+     * </p>
+     * @return nodes that were used to split the outside touching part
+     */
+    public NodesSet<S> getOutsideSplitters() {
+        return outsideSplitters;
+    }
+
     /** Check if the cut sub-hyperplane touches inside cells.
      * @return true if the cut sub-hyperplane touches inside cells
      */
@@ -144,4 +176,15 @@ class Characterization<S extends Space> {
         return insideTouching;
     }
 
+    /** Get the nodes that were used to split the inside touching part.
+     * <p>
+     * Splitting nodes are internal nodes (i.e. they have a non-null
+     * cut sub-hyperplane).
+     * </p>
+     * @return nodes that were used to split the inside touching part
+     */
+    public NodesSet<S> getInsideSplitters() {
+        return insideSplitters;
+    }
+
 }

http://git-wip-us.apache.org/repos/asf/commons-math/blob/6525cc40/src/main/java/org/apache/commons/math3/geometry/partitioning/NodesSet.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/NodesSet.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/NodesSet.java
new file mode 100644
index 0000000..688279a
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/NodesSet.java
@@ -0,0 +1,72 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math3.geometry.partitioning;
+
+import java.util.ArrayList;
+import java.util.Iterator;
+import java.util.List;
+
+import org.apache.commons.math3.geometry.Space;
+
+/** Set of {@link BSPTree BSP tree} nodes.
+ * @see BoundaryAttribute
+ * @param <S> Type of the space.
+ * @since 3.4
+ */
+public class NodesSet<S extends Space> implements Iterable<BSPTree<S>> {
+
+    /** List of sub-hyperplanes. */
+    private List<BSPTree<S>> list;
+
+    /** Simple constructor.
+     */
+    public NodesSet() {
+        list = new ArrayList<BSPTree<S>>();
+    }
+
+    /** Add a node if not already known.
+     * @param node node to add
+     */
+    public void add(final BSPTree<S> node) {
+
+        for (final BSPTree<S> existing : list) {
+            if (node == existing) {
+                // the node is already known, don't add it
+                return;
+            }
+        }
+
+        // the node was not known, add it
+        list.add(node);
+
+    }
+
+    /** Add nodes if they are not already known.
+     * @param iterator nodes iterator
+     */
+    public void addAll(final Iterable<BSPTree<S>> iterator) {
+        for (final BSPTree<S> node : iterator) {
+            add(node);
+        }
+    }
+
+    /** {@inheritDoc} */
+    public Iterator<BSPTree<S>> iterator() {
+        return list.iterator();
+    }
+
+}

http://git-wip-us.apache.org/repos/asf/commons-math/blob/6525cc40/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
index 47b28bd..8cf2603 100644
--- a/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
@@ -16,6 +16,9 @@
  */
 package org.apache.commons.math3.geometry.partitioning;
 
+import java.util.HashMap;
+import java.util.Map;
+
 import org.apache.commons.math3.geometry.Point;
 import org.apache.commons.math3.geometry.Space;
 import org.apache.commons.math3.geometry.euclidean.twod.PolygonsSet;
@@ -129,6 +132,11 @@ public class RegionFactory<S extends Space> {
      * region independent region will be built
      * @return a new region, complement of the specified one
      */
+    /** Get the complement of the region (exchanged interior/exterior).
+     * @param region region to complement, it will not modified, a new
+     * region independent region will be built
+     * @return a new region, complement of the specified one
+     */
     public Region<S> getComplement(final Region<S> region) {
         return region.buildNew(recurseComplement(region.getTree(false)));
     }
@@ -138,24 +146,62 @@ public class RegionFactory<S extends Space> {
      * @return new tree, complement of the node
      */
     private BSPTree<S> recurseComplement(final BSPTree<S> node) {
-        if (node.getCut() == null) {
-            return new BSPTree<S>(((Boolean) node.getAttribute()) ? Boolean.FALSE : Boolean.TRUE);
+
+        // transform the tree, except for boundary attribute splitters
+        final Map<BSPTree<S>, BSPTree<S>> map = new HashMap<BSPTree<S>, BSPTree<S>>();
+        final BSPTree<S> transformedTree = recurseComplement(node, map);
+
+        // set up the boundary attributes splitters
+        for (final Map.Entry<BSPTree<S>, BSPTree<S>> entry : map.entrySet()) {
+            if (entry.getKey().getCut() != null) {
+                @SuppressWarnings("unchecked")
+                BoundaryAttribute<S> original = (BoundaryAttribute<S>) entry.getKey().getAttribute();
+                if (original != null) {
+                    @SuppressWarnings("unchecked")
+                    BoundaryAttribute<S> transformed = (BoundaryAttribute<S>) entry.getValue().getAttribute();
+                    for (final BSPTree<S> splitter : original.getSplitters()) {
+                        transformed.getSplitters().add(map.get(splitter));
+                    }
+                }
+            }
         }
 
-        @SuppressWarnings("unchecked")
-        BoundaryAttribute<S> attribute = (BoundaryAttribute<S>) node.getAttribute();
-        if (attribute != null) {
-            final SubHyperplane<S> plusOutside =
-                (attribute.getPlusInside() == null) ? null : attribute.getPlusInside().copySelf();
-            final SubHyperplane<S> plusInside  =
-                (attribute.getPlusOutside() == null) ? null : attribute.getPlusOutside().copySelf();
-            attribute = new BoundaryAttribute<S>(plusOutside, plusInside);
+        return transformedTree;
+
+    }
+
+    /** Recursively build the complement of a BSP tree.
+     * @param node current node of the original tree
+     * @param map transformed nodes map
+     * @return new tree, complement of the node
+     */
+    private BSPTree<S> recurseComplement(final BSPTree<S> node,
+                                         final Map<BSPTree<S>, BSPTree<S>> map) {
+
+        final BSPTree<S> transformedNode;
+        if (node.getCut() == null) {
+            transformedNode = new BSPTree<S>(((Boolean) node.getAttribute()) ? Boolean.FALSE : Boolean.TRUE);
+        } else {
+
+            @SuppressWarnings("unchecked")
+            BoundaryAttribute<S> attribute = (BoundaryAttribute<S>) node.getAttribute();
+            if (attribute != null) {
+                final SubHyperplane<S> plusOutside =
+                        (attribute.getPlusInside() == null) ? null : attribute.getPlusInside().copySelf();
+                final SubHyperplane<S> plusInside  =
+                        (attribute.getPlusOutside() == null) ? null : attribute.getPlusOutside().copySelf();
+                // we start with an empty list of splitters, it will be filled in out of recursion
+                attribute = new BoundaryAttribute<S>(plusOutside, plusInside, new NodesSet<S>());
+            }
+
+            transformedNode = new BSPTree<S>(node.getCut().copySelf(),
+                                             recurseComplement(node.getPlus(),  map),
+                                             recurseComplement(node.getMinus(), map),
+                                             attribute);
         }
 
-        return new BSPTree<S>(node.getCut().copySelf(),
-                              recurseComplement(node.getPlus()),
-                              recurseComplement(node.getMinus()),
-                              attribute);
+        map.put(node, transformedNode);
+        return transformedNode;
 
     }
 


[5/8] [math] set up a weak link between a line and its reverse to avoid duplications.

Posted by lu...@apache.org.
set up a weak link between a line and its reverse to avoid duplications.


Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/84cf29ec
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/84cf29ec
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/84cf29ec

Branch: refs/heads/master
Commit: 84cf29ecf31c67c7640f8821d4fdda625a173aa6
Parents: 6525cc4
Author: Luc Maisonobe <lu...@apache.org>
Authored: Tue Dec 2 15:37:15 2014 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 2 15:37:15 2014 +0100

----------------------------------------------------------------------
 .../math3/geometry/euclidean/twod/Line.java     | 39 ++++++++++++++++++--
 1 file changed, 36 insertions(+), 3 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/84cf29ec/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java
index b40f0fa..d6fd487 100644
--- a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java
+++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java
@@ -80,6 +80,9 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
     /** Tolerance below which points are considered identical. */
     private final double tolerance;
 
+    /** Reverse line. */
+    private Line reverse;
+
     /** Build a line from two points.
      * <p>The line is oriented from p1 to p2</p>
      * @param p1 first point
@@ -118,6 +121,7 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
         this.sin          = sin;
         this.originOffset = originOffset;
         this.tolerance    = tolerance;
+        this.reverse      = null;
     }
 
     /** Build a line from two points.
@@ -152,6 +156,7 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
         sin          = line.sin;
         originOffset = line.originOffset;
         tolerance    = line.tolerance;
+        reverse      = null;
     }
 
     /** {@inheritDoc} */
@@ -165,6 +170,7 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
      * @param p2 second point
      */
     public void reset(final Vector2D p1, final Vector2D p2) {
+        unlinkReverse();
         final double dx = p2.getX() - p1.getX();
         final double dy = p2.getY() - p1.getY();
         final double d = FastMath.hypot(dx, dy);
@@ -186,6 +192,7 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
      * @param alpha angle of the line with respect to abscissa axis
      */
     public void reset(final Vector2D p, final double alpha) {
+        unlinkReverse();
         this.angle   = MathUtils.normalizeAngle(alpha, FastMath.PI);
         cos          = FastMath.cos(this.angle);
         sin          = FastMath.sin(this.angle);
@@ -195,6 +202,7 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
     /** Revert the instance.
      */
     public void revertSelf() {
+        unlinkReverse();
         if (angle < FastMath.PI) {
             angle += FastMath.PI;
         } else {
@@ -205,14 +213,37 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
         originOffset = -originOffset;
     }
 
+    /** Unset the link between an instance and its reverse.
+     */
+    private void unlinkReverse() {
+        if (reverse != null) {
+            reverse.reverse = null;
+        }
+        reverse = null;
+    }
+
     /** Get the reverse of the instance.
      * <p>Get a line with reversed orientation with respect to the
-     * instance. A new object is built, the instance is untouched.</p>
+     * instance.</p>
+     * <p>
+     * As long as neither the instance nor its reverse are modified
+     * (i.e. as long as none of the {@link #reset(Vector2D, Vector2D)},
+     * {@link #reset(Vector2D, double)}, {@link #revertSelf()},
+     * {@link #setAngle(double)} or {@link #setOriginOffset(double)}
+     * methods are called), then the line and its reverse remain linked
+     * together so that {@code line.getReverse().getReverse() == line}.
+     * When one of the line is modified, the link is deleted as both
+     * instance becomes independent.
+     * </p>
      * @return a new line, with orientation opposite to the instance orientation
      */
     public Line getReverse() {
-        return new Line((angle < FastMath.PI) ? (angle + FastMath.PI) : (angle - FastMath.PI),
-                        -cos, -sin, -originOffset, tolerance);
+        if (reverse == null) {
+            reverse = new Line((angle < FastMath.PI) ? (angle + FastMath.PI) : (angle - FastMath.PI),
+                               -cos, -sin, -originOffset, tolerance);
+            reverse.reverse = this;
+        }
+        return reverse;
     }
 
     /** Transform a space point into a sub-space point.
@@ -383,6 +414,7 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
      * @param angle new angle of the line with respect to the abscissa axis
      */
     public void setAngle(final double angle) {
+        unlinkReverse();
         this.angle = MathUtils.normalizeAngle(angle, FastMath.PI);
         cos        = FastMath.cos(this.angle);
         sin        = FastMath.sin(this.angle);
@@ -399,6 +431,7 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
      * @param offset offset of the origin
      */
     public void setOriginOffset(final double offset) {
+        unlinkReverse();
         originOffset = offset;
     }
 


[6/8] [math] Removed development traces.

Posted by lu...@apache.org.
Removed development traces.

Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/893ef53c
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/893ef53c
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/893ef53c

Branch: refs/heads/master
Commit: 893ef53c6cb9b243ba8f28c8cad6e42c8019af0c
Parents: 84cf29e
Author: Luc Maisonobe <lu...@apache.org>
Authored: Tue Dec 2 22:06:33 2014 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 2 22:06:33 2014 +0100

----------------------------------------------------------------------
 .../commons/math3/geometry/partitioning/RegionFactory.java   | 8 --------
 1 file changed, 8 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/893ef53c/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
index 8cf2603..16d4472 100644
--- a/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
@@ -21,8 +21,6 @@ import java.util.Map;
 
 import org.apache.commons.math3.geometry.Point;
 import org.apache.commons.math3.geometry.Space;
-import org.apache.commons.math3.geometry.euclidean.twod.PolygonsSet;
-import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
 import org.apache.commons.math3.geometry.partitioning.BSPTree.VanishingCutHandler;
 import org.apache.commons.math3.geometry.partitioning.Region.Location;
 
@@ -296,12 +294,6 @@ public class RegionFactory<S extends Space> {
             // get a representative point in the degenerate cell
             final BSPTree<S> cell = node.pruneAroundConvexCell(Boolean.TRUE, Boolean.FALSE, null);
             final Region<S> r = region1.buildNew(cell);
-            for (Vector2D[] loop : ((PolygonsSet) r).getVertices()) {
-                System.out.format(java.util.Locale.US, "%n");
-                for (Vector2D v : loop) {
-                    System.out.format(java.util.Locale.US, "%14.10f %14.10f%n", v.getX(), v.getY());
-                }
-            }
             final Point<S> p = r.getBarycenter();
             return new BSPTree<S>(region1.checkPoint(p) == Location.INSIDE &&
                                   region2.checkPoint(p) == Location.OUTSIDE);


[3/8] [math] Improved accuracy of internal computation.

Posted by lu...@apache.org.
Improved accuracy of internal computation.

The improvement is obtained by using high accuracy linear combinations
and by preserving cos/sin as much as possible. This last part is
important in the frequent case of lines along y axis, i.e. when angle is
π/2. We want the cosine to remain 0 and not get about 10^-17. 

Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/5f667c03
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/5f667c03
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/5f667c03

Branch: refs/heads/master
Commit: 5f667c031cdf02adbee9ef83cb71990333c58ba3
Parents: c29e3d6
Author: Luc Maisonobe <lu...@apache.org>
Authored: Sun Nov 30 11:44:07 2014 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 2 15:24:31 2014 +0100

----------------------------------------------------------------------
 .../math3/geometry/euclidean/twod/Line.java     | 54 ++++++++++----------
 1 file changed, 28 insertions(+), 26 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/5f667c03/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java
index 70ab1ac..b40f0fa 100644
--- a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java
+++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/Line.java
@@ -31,6 +31,7 @@ import org.apache.commons.math3.geometry.partitioning.Hyperplane;
 import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
 import org.apache.commons.math3.geometry.partitioning.Transform;
 import org.apache.commons.math3.util.FastMath;
+import org.apache.commons.math3.util.MathArrays;
 import org.apache.commons.math3.util.MathUtils;
 
 /** This class represents an oriented line in the 2D plane.
@@ -147,8 +148,8 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
      */
     public Line(final Line line) {
         angle        = MathUtils.normalizeAngle(line.angle, FastMath.PI);
-        cos          = FastMath.cos(angle);
-        sin          = FastMath.sin(angle);
+        cos          = line.cos;
+        sin          = line.sin;
         originOffset = line.originOffset;
         tolerance    = line.tolerance;
     }
@@ -174,9 +175,9 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
             originOffset = p1.getY();
         } else {
             angle        = FastMath.PI + FastMath.atan2(-dy, -dx);
-            cos          = FastMath.cos(angle);
-            sin          = FastMath.sin(angle);
-            originOffset = (p2.getX() * p1.getY() - p1.getX() * p2.getY()) / d;
+            cos          = dx / d;
+            sin          = dy / d;
+            originOffset = MathArrays.linearCombination(p2.getX(), p1.getY(), -p1.getX(), p2.getY()) / d;
         }
     }
 
@@ -188,7 +189,7 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
         this.angle   = MathUtils.normalizeAngle(alpha, FastMath.PI);
         cos          = FastMath.cos(this.angle);
         sin          = FastMath.sin(this.angle);
-        originOffset = cos * p.getY() - sin * p.getX();
+        originOffset = MathArrays.linearCombination(cos, p.getY(), -sin, p.getX());
     }
 
     /** Revert the instance.
@@ -235,14 +236,14 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
     /** {@inheritDoc} */
     public Vector1D toSubSpace(final Point<Euclidean2D> point) {
         Vector2D p2 = (Vector2D) point;
-        return new Vector1D(cos * p2.getX() + sin * p2.getY());
+        return new Vector1D(MathArrays.linearCombination(cos, p2.getX(), sin, p2.getY()));
     }
 
     /** {@inheritDoc} */
     public Vector2D toSpace(final Point<Euclidean1D> point) {
         final double abscissa = ((Vector1D) point).getX();
-        return new Vector2D(abscissa * cos - originOffset * sin,
-                            abscissa * sin + originOffset * cos);
+        return new Vector2D(MathArrays.linearCombination(abscissa, cos, -originOffset, sin),
+                            MathArrays.linearCombination(abscissa, sin,  originOffset, cos));
     }
 
     /** Get the intersection point of the instance and another line.
@@ -251,12 +252,12 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
      * or null if there are no intersection points
      */
     public Vector2D intersection(final Line other) {
-        final double d = sin * other.cos - other.sin * cos;
+        final double d = MathArrays.linearCombination(sin, other.cos, -other.sin, cos);
         if (FastMath.abs(d) < tolerance) {
             return null;
         }
-        return new Vector2D((cos * other.originOffset - other.cos * originOffset) / d,
-                            (sin * other.originOffset - other.sin * originOffset) / d);
+        return new Vector2D(MathArrays.linearCombination(cos, other.originOffset, -other.cos, originOffset) / d,
+                            MathArrays.linearCombination(sin, other.originOffset, -other.sin, originOffset) / d);
     }
 
     /** {@inheritDoc}
@@ -298,7 +299,7 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
      */
     public double getOffset(final Line line) {
         return originOffset +
-               ((cos * line.cos + sin * line.sin > 0) ? -line.originOffset : line.originOffset);
+               (MathArrays.linearCombination(cos, line.cos, sin, line.sin) > 0 ? -line.originOffset : line.originOffset);
     }
 
     /** Get the offset (oriented distance) of a vector.
@@ -312,13 +313,13 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
     /** {@inheritDoc} */
     public double getOffset(final Point<Euclidean2D> point) {
         Vector2D p2 = (Vector2D) point;
-        return sin * p2.getX() - cos * p2.getY() + originOffset;
+        return MathArrays.linearCombination(sin, p2.getX(), -cos, p2.getY(), 1.0, originOffset);
     }
 
     /** {@inheritDoc} */
     public boolean sameOrientationAs(final Hyperplane<Euclidean2D> other) {
         final Line otherL = (Line) other;
-        return (sin * otherL.sin + cos * otherL.cos) >= 0.0;
+        return MathArrays.linearCombination(sin, otherL.sin, cos, otherL.cos) >= 0.0;
     }
 
     /** Get one point from the plane.
@@ -330,7 +331,8 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
     public Vector2D getPointAt(final Vector1D abscissa, final double offset) {
         final double x       = abscissa.getX();
         final double dOffset = offset - originOffset;
-        return new Vector2D(x * cos + dOffset * sin, x * sin - dOffset * cos);
+        return new Vector2D(MathArrays.linearCombination(x, cos,  dOffset, sin),
+                            MathArrays.linearCombination(x, sin, -dOffset, cos));
     }
 
     /** Check if the line contains a point.
@@ -360,14 +362,14 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
      * (they can have either the same or opposite orientations)
      */
     public boolean isParallelTo(final Line line) {
-        return FastMath.abs(sin * line.cos - cos * line.sin) < tolerance;
+        return FastMath.abs(MathArrays.linearCombination(sin, line.cos, -cos, line.sin)) < tolerance;
     }
 
     /** Translate the line to force it passing by a point.
      * @param p point by which the line should pass
      */
     public void translateToPoint(final Vector2D p) {
-        originOffset = cos * p.getY() - sin * p.getX();
+        originOffset = MathArrays.linearCombination(cos, p.getY(), -sin, p.getX());
     }
 
     /** Get the angle of the line.
@@ -457,9 +459,9 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
             cYY = m[3];
             cY1 = m[5];
 
-            c1Y = cXY * cY1 - cYY * cX1;
-            c1X = cXX * cY1 - cYX * cX1;
-            c11 = cXX * cYY - cYX * cXY;
+            c1Y = MathArrays.linearCombination(cXY, cY1, -cYY, cX1);
+            c1X = MathArrays.linearCombination(cXX, cY1, -cYX, cX1);
+            c11 = MathArrays.linearCombination(cXX, cYY, -cYX, cXY);
 
             if (FastMath.abs(c11) < 1.0e-20) {
                 throw new MathIllegalArgumentException(LocalizedFormats.NON_INVERTIBLE_TRANSFORM);
@@ -472,16 +474,16 @@ public class Line implements Hyperplane<Euclidean2D>, Embedding<Euclidean2D, Euc
             final Vector2D p2D = (Vector2D) point;
             final double  x   = p2D.getX();
             final double  y   = p2D.getY();
-            return new Vector2D(cXX * x + cXY * y + cX1,
-                               cYX * x + cYY * y + cY1);
+            return new Vector2D(MathArrays.linearCombination(cXX, x, cXY, y, cX1, 1),
+                                MathArrays.linearCombination(cYX, x, cYY, y, cY1, 1));
         }
 
         /** {@inheritDoc} */
         public Line apply(final Hyperplane<Euclidean2D> hyperplane) {
             final Line   line    = (Line) hyperplane;
-            final double rOffset = c1X * line.cos + c1Y * line.sin + c11 * line.originOffset;
-            final double rCos    = cXX * line.cos + cXY * line.sin;
-            final double rSin    = cYX * line.cos + cYY * line.sin;
+            final double rOffset = MathArrays.linearCombination(c1X, line.cos, c1Y, line.sin, c11, line.originOffset);
+            final double rCos    = MathArrays.linearCombination(cXX, line.cos, cXY, line.sin);
+            final double rSin    = MathArrays.linearCombination(cYX, line.cos, cYY, line.sin);
             final double inv     = 1.0 / FastMath.sqrt(rSin * rSin + rCos * rCos);
             return new Line(FastMath.PI + FastMath.atan2(-rSin, -rCos),
                             inv * rCos, inv * rSin,


[8/8] [math] Deprecated AVLTree and OrderedTuple.

Posted by lu...@apache.org.
Deprecated AVLTree and OrderedTuple.

These two classes were dirty hacks used only for the extraction of
boundary from 2D Euclidean BSP trees. They were not really efficient and
produces some problems as topology was not preserved.

A completely new implementation of boundary extraction has been set up
as of 3.4, that preserves more topological information and therefore
should be more robust. This new implementation does not use these
classes anymore.

As AVLTree is not really a math concept, this class is not considered to
belong to [math] scope. OrderedTuple is more a grey zone, but it is
really something inefficient and no need for it is foreseen.

Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/d21d3a94
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/d21d3a94
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/d21d3a94

Branch: refs/heads/master
Commit: d21d3a94abdbcf5a890fb1b5de28f1060e9e584d
Parents: e6aae3a
Author: Luc Maisonobe <lu...@apache.org>
Authored: Tue Dec 2 22:18:57 2014 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 2 22:18:57 2014 +0100

----------------------------------------------------------------------
 .../commons/math3/geometry/partitioning/utilities/AVLTree.java    | 3 +++
 .../math3/geometry/partitioning/utilities/OrderedTuple.java       | 3 +++
 .../math3/geometry/partitioning/utilities/AVLTreeTest.java        | 1 +
 3 files changed, 7 insertions(+)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/d21d3a94/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTree.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTree.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTree.java
index 36826c6..9412421 100644
--- a/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTree.java
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTree.java
@@ -42,7 +42,10 @@ package org.apache.commons.math3.geometry.partitioning.utilities;
  * @param <T> the type of the elements
  *
  * @since 3.0
+ * @deprecated as of 3.4, this class is not used anymore and considered
+ * to be out of scope of Apache Commons Math
  */
+@Deprecated
 public class AVLTree<T extends Comparable<T>> {
 
     /** Top level node. */

http://git-wip-us.apache.org/repos/asf/commons-math/blob/d21d3a94/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/OrderedTuple.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/OrderedTuple.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/OrderedTuple.java
index 6c46be7..2dad2d7 100644
--- a/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/OrderedTuple.java
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/utilities/OrderedTuple.java
@@ -89,7 +89,10 @@ import org.apache.commons.math3.util.FastMath;
  * </ul>
  *
  * @since 3.0
+ * @deprecated as of 3.4, this class is not used anymore and considered
+ * to be out of scope of Apache Commons Math
  */
+@Deprecated
 public class OrderedTuple implements Comparable<OrderedTuple> {
 
     /** Sign bit mask. */

http://git-wip-us.apache.org/repos/asf/commons-math/blob/d21d3a94/src/test/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTreeTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTreeTest.java b/src/test/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTreeTest.java
index 9d7f424..49cfef6 100644
--- a/src/test/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTreeTest.java
+++ b/src/test/java/org/apache/commons/math3/geometry/partitioning/utilities/AVLTreeTest.java
@@ -20,6 +20,7 @@ import org.apache.commons.math3.geometry.partitioning.utilities.AVLTree;
 import org.junit.Assert;
 import org.junit.Test;
 
+@Deprecated
 public class AVLTreeTest {
 
     @Test


[4/8] [math] Identified a new failure mode for polygons.

Posted by lu...@apache.org.
Identified a new failure mode for polygons.

Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/c29e3d68
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/c29e3d68
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/c29e3d68

Branch: refs/heads/master
Commit: c29e3d681f2a46cc3288dc4f38925efbcf30c233
Parents: e11c000
Author: Luc Maisonobe <lu...@apache.org>
Authored: Sun Nov 30 11:40:27 2014 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 2 15:24:31 2014 +0100

----------------------------------------------------------------------
 .../euclidean/twod/PolygonsSetTest.java         | 37 ++++++++++++++++++++
 1 file changed, 37 insertions(+)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/c29e3d68/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java b/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java
index 348ace8..c6c68a0 100644
--- a/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java
+++ b/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java
@@ -24,6 +24,7 @@ import org.apache.commons.math3.geometry.euclidean.oned.IntervalsSet;
 import org.apache.commons.math3.geometry.euclidean.oned.Vector1D;
 import org.apache.commons.math3.geometry.partitioning.BSPTree;
 import org.apache.commons.math3.geometry.partitioning.BoundaryProjection;
+import org.apache.commons.math3.geometry.partitioning.Hyperplane;
 import org.apache.commons.math3.geometry.partitioning.Region;
 import org.apache.commons.math3.geometry.partitioning.Region.Location;
 import org.apache.commons.math3.geometry.partitioning.RegionFactory;
@@ -1107,6 +1108,42 @@ public class PolygonsSetTest {
 
     }
 
+    @Test
+    public void testThinRectangle() {
+
+        RegionFactory<Euclidean2D> factory = new RegionFactory<Euclidean2D>();
+        Vector2D pA = new Vector2D(0.0,        1.0);
+        Vector2D pB = new Vector2D(0.0,        0.0);
+        Vector2D pC = new Vector2D(1.0 / 64.0, 0.0);
+        Vector2D pD = new Vector2D(1.0 / 64.0, 1.0);
+
+        // if tolerance is smaller than rectangle width, the rectangle is computed accurately
+        Hyperplane<Euclidean2D>[] h1 = new Line[] {
+            new Line(pA, pB, 1.0 / 256),
+            new Line(pB, pC, 1.0 / 256),
+            new Line(pC, pD, 1.0 / 256),
+            new Line(pD, pA, 1.0 / 256)
+        };
+        Region<Euclidean2D> accuratePolygon = factory.buildConvex(h1);
+        Assert.assertEquals(1.0 / 64.0, accuratePolygon.getSize(), 1.0e-10);
+        Assert.assertTrue(Double.isInfinite(new RegionFactory<Euclidean2D>().getComplement(accuratePolygon).getSize()));
+        Assert.assertEquals(2 * (1.0 + 1.0 / 64.0), accuratePolygon.getBoundarySize(), 1.0e-10);
+
+        // if tolerance is larger than rectangle width, the rectangle degenerates
+        // as of 3.3, its two long edges cannot be distinguished anymore and this part of the test fails
+        Hyperplane<Euclidean2D>[] h2 = new Line[] {
+            new Line(pA, pB, 1.0 / 16),
+            new Line(pB, pC, 1.0 / 16),
+            new Line(pC, pD, 1.0 / 16),
+            new Line(pD, pA, 1.0 / 16)
+        };
+        Region<Euclidean2D> degeneratedPolygon = factory.buildConvex(h2);
+        Assert.assertEquals(1.0 / 64.0, degeneratedPolygon.getSize(), 1.0e-10);
+        Assert.assertTrue(Double.isInfinite(new RegionFactory<Euclidean2D>().getComplement(degeneratedPolygon).getSize()));
+        Assert.assertEquals(2 * (1.0 + 1.0 / 64.0), degeneratedPolygon.getBoundarySize(), 1.0e-10);
+
+    }
+
     private PolygonsSet buildSet(Vector2D[][] vertices) {
         ArrayList<SubHyperplane<Euclidean2D>> edges = new ArrayList<SubHyperplane<Euclidean2D>>();
         for (int i = 0; i < vertices.length; ++i) {


[7/8] [math] Fixed a problem with too thin polygons considered to have infinite size.

Posted by lu...@apache.org.
Fixed a problem with too thin polygons considered to have infinite size.


Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/e6aae3a8
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/e6aae3a8
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/e6aae3a8

Branch: refs/heads/master
Commit: e6aae3a8bffb981b8dce19642cc1ab831a428cae
Parents: 893ef53
Author: Luc Maisonobe <lu...@apache.org>
Authored: Tue Dec 2 22:09:06 2014 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 2 22:09:06 2014 +0100

----------------------------------------------------------------------
 src/changes/changes.xml                         |   3 +
 .../geometry/euclidean/twod/PolygonsSet.java    | 414 +++++++++++++------
 .../euclidean/twod/PolygonsSetTest.java         |  60 +--
 3 files changed, 330 insertions(+), 147 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/e6aae3a8/src/changes/changes.xml
----------------------------------------------------------------------
diff --git a/src/changes/changes.xml b/src/changes/changes.xml
index 46ad704..a4b462c 100644
--- a/src/changes/changes.xml
+++ b/src/changes/changes.xml
@@ -73,6 +73,9 @@ Users are encouraged to upgrade to this version as this release not
   2. A few methods in the FastMath class are in fact slower that their
   counterpart in either Math or StrictMath (cf. MATH-740 and MATH-901).
 ">
+      <action dev="luc" type="fix" issue="MATH-1174" >
+        Fixed a problem with too thin polygons considered to have infinite size.
+      </action>
       <action dev="luc" type="add" >
         Boundary attributes in regions now provides the BSP tree nodes that
         were used to split the sub-hyperplane froming the boundary part of the facet. 

http://git-wip-us.apache.org/repos/asf/commons-math/blob/e6aae3a8/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java
index 532d6a0..8361831 100644
--- a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java
+++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java
@@ -20,7 +20,6 @@ import java.util.ArrayList;
 import java.util.Collection;
 import java.util.List;
 
-import org.apache.commons.math3.exception.MathInternalError;
 import org.apache.commons.math3.geometry.Point;
 import org.apache.commons.math3.geometry.euclidean.oned.Euclidean1D;
 import org.apache.commons.math3.geometry.euclidean.oned.Interval;
@@ -31,10 +30,9 @@ import org.apache.commons.math3.geometry.partitioning.AbstractSubHyperplane;
 import org.apache.commons.math3.geometry.partitioning.BSPTree;
 import org.apache.commons.math3.geometry.partitioning.BSPTreeVisitor;
 import org.apache.commons.math3.geometry.partitioning.BoundaryAttribute;
+import org.apache.commons.math3.geometry.partitioning.Hyperplane;
 import org.apache.commons.math3.geometry.partitioning.Side;
 import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
-import org.apache.commons.math3.geometry.partitioning.utilities.AVLTree;
-import org.apache.commons.math3.geometry.partitioning.utilities.OrderedTuple;
 import org.apache.commons.math3.util.FastMath;
 
 /** This class represents a 2D region: a set of polygons.
@@ -689,19 +687,34 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
                 vertices = new Vector2D[0][];
             } else {
 
-                // sort the segments according to their start point
-                final SegmentsBuilder visitor = new SegmentsBuilder();
+                // build the unconnected segments
+                final SegmentsBuilder visitor = new SegmentsBuilder(getTolerance());
                 getTree(true).visit(visitor);
-                final AVLTree<ComparableSegment> sorted = visitor.getSorted();
-
-                // identify the loops, starting from the open ones
-                // (their start segments are naturally at the sorted set beginning)
-                final ArrayList<List<ComparableSegment>> loops = new ArrayList<List<ComparableSegment>>();
-                while (!sorted.isEmpty()) {
-                    final AVLTree<ComparableSegment>.Node node = sorted.getSmallest();
-                    final List<ComparableSegment> loop = followLoop(node, sorted);
+                final List<ConnectableSegment> segments = visitor.getSegments();
+
+                // connect all segments, using topological criteria first
+                // and using Euclidean distance only as a last resort
+                int pending = segments.size();
+                pending -= naturalFollowerConnections(segments);
+                if (pending > 0) {
+                    pending -= splitEdgeConnections(segments);
+                }
+                if (pending > 0) {
+                    pending -= closeVerticesConnections(segments);
+                }
+
+                // create the segment loops
+                final ArrayList<List<ConnectableSegment>> loops = new ArrayList<List<ConnectableSegment>>();
+                for (ConnectableSegment s = getUnprocessed(segments); s != null; s = getUnprocessed(segments)) {
+                    final List<ConnectableSegment> loop = followLoop(s);
                     if (loop != null) {
-                        loops.add(loop);
+                        if (loop.get(0).getStart() == null) {
+                            // this is an open loop, we put it on the front
+                            loops.add(0, loop);
+                        } else {
+                            // this is a closed loop, we put it on the back
+                            loops.add(loop);
+                        }
                     }
                 }
 
@@ -709,7 +722,7 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
                 vertices = new Vector2D[loops.size()][];
                 int i = 0;
 
-                for (final List<ComparableSegment> loop : loops) {
+                for (final List<ConnectableSegment> loop : loops) {
                     if (loop.size() < 2 ||
                         (loop.size() == 2 && loop.get(0).getStart() == null && loop.get(1).getEnd() == null)) {
                         // single infinite line
@@ -764,126 +777,245 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
 
     }
 
-    /** Follow a boundary loop.
-     * @param node node containing the segment starting the loop
-     * @param sorted set of segments belonging to the boundary, sorted by
-     * start points (contains {@code node})
-     * @return a list of connected sub-hyperplanes starting at
-     * {@code node}
+    /** Connect the segments using only natural follower information.
+     * @param segments segments complete segments list
+     * @return number of connections performed
      */
-    private List<ComparableSegment> followLoop(final AVLTree<ComparableSegment>.Node node,
-                                               final AVLTree<ComparableSegment> sorted) {
-
-        final ArrayList<ComparableSegment> loop = new ArrayList<ComparableSegment>();
-        ComparableSegment segment = node.getElement();
-        loop.add(segment);
-        final Vector2D globalStart = segment.getStart();
-        Vector2D end = segment.getEnd();
-        node.delete();
-
-        // is this an open or a closed loop ?
-        final boolean open = segment.getStart() == null;
-
-        while ((end != null) && (open || (globalStart.distance((Point<Euclidean2D>) end) > getTolerance()))) {
-
-            // search the sub-hyperplane starting where the previous one ended
-            AVLTree<ComparableSegment>.Node selectedNode = null;
-            ComparableSegment       selectedSegment  = null;
-            double                  selectedDistance = Double.POSITIVE_INFINITY;
-            final ComparableSegment lowerLeft        = new ComparableSegment(end, -1.0e-10, -1.0e-10);
-            final ComparableSegment upperRight       = new ComparableSegment(end, +1.0e-10, +1.0e-10);
-            for (AVLTree<ComparableSegment>.Node n = sorted.getNotSmaller(lowerLeft);
-                 (n != null) && (n.getElement().compareTo(upperRight) <= 0);
-                 n = n.getNext()) {
-                segment = n.getElement();
-                final double distance = end.distance((Point<Euclidean2D>) segment.getStart());
-                if (distance < selectedDistance) {
-                    selectedNode     = n;
-                    selectedSegment  = segment;
-                    selectedDistance = distance;
+    private int naturalFollowerConnections(final List<ConnectableSegment> segments) {
+        int connected = 0;
+        for (final ConnectableSegment segment : segments) {
+            if (segment.getNext() == null) {
+                final BSPTree<Euclidean2D> node = segment.getNode();
+                final BSPTree<Euclidean2D> end  = segment.getEndNode();
+                for (final ConnectableSegment candidateNext : segments) {
+                    if (candidateNext.getPrevious()  == null &&
+                        candidateNext.getNode()      == end &&
+                        candidateNext.getStartNode() == node) {
+                        // connect the two segments
+                        segment.setNext(candidateNext);
+                        candidateNext.setPrevious(segment);
+                        ++connected;
+                        break;
+                    }
                 }
             }
+        }
+        return connected;
+    }
 
-            if (selectedDistance > 1.0e-10) {
-                // this is a degenerated loop, it probably comes from a very
-                // tiny region with some segments smaller than the threshold, we
-                // simply ignore it
-                return null;
+    /** Connect the segments resulting from a line splitting a straight edge.
+     * @param segments segments complete segments list
+     * @return number of connections performed
+     */
+    private int splitEdgeConnections(final List<ConnectableSegment> segments) {
+        int connected = 0;
+        for (final ConnectableSegment segment : segments) {
+            if (segment.getNext() == null) {
+                final Hyperplane<Euclidean2D> hyperplane = segment.getNode().getCut().getHyperplane();
+                final BSPTree<Euclidean2D> end  = segment.getEndNode();
+                for (final ConnectableSegment candidateNext : segments) {
+                    if (candidateNext.getPrevious()                      == null &&
+                        candidateNext.getNode().getCut().getHyperplane() == hyperplane &&
+                        candidateNext.getStartNode()                     == end) {
+                        // connect the two segments
+                        segment.setNext(candidateNext);
+                        candidateNext.setPrevious(segment);
+                        ++connected;
+                        break;
+                    }
+                }
             }
+        }
+        return connected;
+    }
 
-            end = selectedSegment.getEnd();
-            loop.add(selectedSegment);
-            selectedNode.delete();
+    /** Connect the segments using Euclidean distance.
+     * <p>
+     * This connection heuristic should be used last, as it relies
+     * only on a fuzzy distance criterion.
+     * </p>
+     * @param segments segments complete segments list
+     * @return number of connections performed
+     */
+    private int closeVerticesConnections(final List<ConnectableSegment> segments) {
+        int connected = 0;
+        for (final ConnectableSegment segment : segments) {
+            if (segment.getNext() == null) {
+                final Vector2D end = segment.getEnd();
+                for (final ConnectableSegment candidateNext : segments) {
+                    if (candidateNext.getPrevious() == null &&
+                        Vector2D.distance(end, candidateNext.getStart()) <= getTolerance()) {
+                        // connect the two segments
+                        segment.setNext(candidateNext);
+                        candidateNext.setPrevious(segment);
+                        ++connected;
+                        break;
+                    }
+                }
+            }
+        }
+        return connected;
+    }
 
+    /** Get first unprocessed segment from a list.
+     * @param segments segments list
+     * @return first segment that has not been processed yet
+     * or null if all segments have been processed
+     */
+    private ConnectableSegment getUnprocessed(final List<ConnectableSegment> segments) {
+        for (final ConnectableSegment segment : segments) {
+            if (!segment.isProcessed()) {
+                return segment;
+            }
         }
+        return null;
+    }
 
-        if ((loop.size() == 2) && !open) {
-            // this is a degenerated infinitely thin loop, we simply ignore it
-            return null;
+    /** Build the loop containing a segment.
+     * <p>
+     * The segment put in the loop will be marked as processed.
+     * </p>
+     * @param defining segment used to define the loop
+     * @return loop containing the segment (may be null if the loop is a
+     * degenerated infinitely thin 2 points loop
+     */
+    private List<ConnectableSegment> followLoop(final ConnectableSegment defining) {
+
+        final List<ConnectableSegment> loop = new ArrayList<ConnectableSegment>();
+        loop.add(defining);
+        defining.setProcessed(true);
+
+        // add segments in connection order
+        ConnectableSegment next = defining.getNext();
+        while (next != defining && next != null) {
+            loop.add(next);
+            next.setProcessed(true);
+            next = next.getNext();
         }
 
-        if ((end == null) && !open) {
-            throw new MathInternalError();
+        if (next == null) {
+            // the loop is open and we have found its end,
+            // we need to find its start too
+            ConnectableSegment previous = defining.getPrevious();
+            while (previous != null) {
+                loop.add(0, previous);
+                previous.setProcessed(true);
+                previous = previous.getPrevious();
+            }
+        } else {
+            if (loop.size() == 2) {
+                // this is a degenerated infinitely thin loop, we simply ignore it
+                return null;
+            }
         }
 
         return loop;
 
     }
 
-    /** Private extension of Segment allowing comparison. */
-    private static class ComparableSegment extends Segment implements Comparable<ComparableSegment> {
+    /** Private extension of Segment allowing connection. */
+    private static class ConnectableSegment extends Segment {
+
+        /** Node containing segment. */
+        private final BSPTree<Euclidean2D> node;
 
-        /** Sorting key. */
-        private OrderedTuple sortingKey;
+        /** Node whose intersection with current node defines start point. */
+        private final BSPTree<Euclidean2D> startNode;
+
+        /** Node whose intersection with current node defines end point. */
+        private final BSPTree<Euclidean2D> endNode;
+
+        /** Previous segment. */
+        private ConnectableSegment previous;
+
+        /** Next segment. */
+        private ConnectableSegment next;
+
+        /** Indicator for completely processed segments. */
+        private boolean processed;
 
         /** Build a segment.
          * @param start start point of the segment
          * @param end end point of the segment
          * @param line line containing the segment
+         * @param node node containing the segment
+         * @param startNode node whose intersection with current node defines start point
+         * @param endNode node whose intersection with current node defines end point
          */
-        public ComparableSegment(final Vector2D start, final Vector2D end, final Line line) {
+        public ConnectableSegment(final Vector2D start, final Vector2D end, final Line line,
+                                  final BSPTree<Euclidean2D> node,
+                                  final BSPTree<Euclidean2D> startNode,
+                                  final BSPTree<Euclidean2D> endNode) {
             super(start, end, line);
-            sortingKey = (start == null) ?
-                         new OrderedTuple(Double.NEGATIVE_INFINITY, Double.NEGATIVE_INFINITY) :
-                         new OrderedTuple(start.getX(), start.getY());
+            this.node      = node;
+            this.startNode = startNode;
+            this.endNode   = endNode;
+            this.previous  = null;
+            this.next      = null;
+            this.processed = false;
         }
 
-        /** Build a dummy segment.
-         * <p>
-         * The object built is not a real segment, only the sorting key is used to
-         * allow searching in the neighborhood of a point. This is an horrible hack ...
-         * </p>
-         * @param start start point of the segment
-         * @param dx abscissa offset from the start point
-         * @param dy ordinate offset from the start point
+        /** Get the node containing segment.
+         * @return node containing segment
          */
-        public ComparableSegment(final Vector2D start, final double dx, final double dy) {
-            super(null, null, null);
-            sortingKey = new OrderedTuple(start.getX() + dx, start.getY() + dy);
+        public BSPTree<Euclidean2D> getNode() {
+            return node;
         }
 
-        /** {@inheritDoc} */
-        public int compareTo(final ComparableSegment o) {
-            return sortingKey.compareTo(o.sortingKey);
+        /** Get the node whose intersection with current node defines start point.
+         * @return node whose intersection with current node defines start point
+         */
+        public BSPTree<Euclidean2D> getStartNode() {
+            return startNode;
         }
 
-        /** {@inheritDoc} */
-        @Override
-        public boolean equals(final Object other) {
-            if (this == other) {
-                return true;
-            } else if (other instanceof ComparableSegment) {
-                return compareTo((ComparableSegment) other) == 0;
-            } else {
-                return false;
-            }
+        /** Get the node whose intersection with current node defines end point.
+         * @return node whose intersection with current node defines end point
+         */
+        public BSPTree<Euclidean2D> getEndNode() {
+            return endNode;
         }
 
-        /** {@inheritDoc} */
-        @Override
-        public int hashCode() {
-            return getStart().hashCode() ^ getEnd().hashCode() ^
-                   getLine().hashCode() ^ sortingKey.hashCode();
+        /** Get the previous segment.
+         * @return previous segment
+         */
+        public ConnectableSegment getPrevious() {
+            return previous;
+        }
+
+        /** Set the previous segment.
+         * @param previous previous segment
+         */
+        public void setPrevious(final ConnectableSegment previous) {
+            this.previous = previous;
+        }
+
+        /** Get the next segment.
+         * @return next segment
+         */
+        public ConnectableSegment getNext() {
+            return next;
+        }
+
+        /** Set the next segment.
+         * @param next previous segment
+         */
+        public void setNext(final ConnectableSegment next) {
+            this.next = next;
+        }
+
+        /** Set the processed flag.
+         * @param processed processed flag to set
+         */
+        public void setProcessed(final boolean processed) {
+            this.processed = processed;
+        }
+
+        /** Check if the segment has been processed.
+         * @return true if the segment has been processed
+         */
+        public boolean isProcessed() {
+            return processed;
         }
 
     }
@@ -891,12 +1023,18 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
     /** Visitor building segments. */
     private static class SegmentsBuilder implements BSPTreeVisitor<Euclidean2D> {
 
-        /** Sorted segments. */
-        private AVLTree<ComparableSegment> sorted;
+        /** Tolerance for close nodes connection. */
+        private final double tolerance;
 
-        /** Simple constructor. */
-        public SegmentsBuilder() {
-            sorted = new AVLTree<ComparableSegment>();
+        /** Built segments. */
+        private final List<ConnectableSegment> segments;
+
+        /** Simple constructor.
+         * @param tolerance tolerance for close nodes connection
+         */
+        public SegmentsBuilder(final double tolerance) {
+            this.tolerance = tolerance;
+            this.segments  = new ArrayList<ConnectableSegment>();
         }
 
         /** {@inheritDoc} */
@@ -908,11 +1046,12 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
         public void visitInternalNode(final BSPTree<Euclidean2D> node) {
             @SuppressWarnings("unchecked")
             final BoundaryAttribute<Euclidean2D> attribute = (BoundaryAttribute<Euclidean2D>) node.getAttribute();
+            final Iterable<BSPTree<Euclidean2D>> splitters = attribute.getSplitters();
             if (attribute.getPlusOutside() != null) {
-                addContribution(attribute.getPlusOutside(), false);
+                addContribution(attribute.getPlusOutside(), node, splitters, false);
             }
             if (attribute.getPlusInside() != null) {
-                addContribution(attribute.getPlusInside(), true);
+                addContribution(attribute.getPlusInside(), node, splitters, true);
             }
         }
 
@@ -922,32 +1061,69 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
 
         /** Add the contribution of a boundary facet.
          * @param sub boundary facet
+         * @param node node containing segment
+         * @param splitters splitters for the boundary facet
          * @param reversed if true, the facet has the inside on its plus side
          */
-        private void addContribution(final SubHyperplane<Euclidean2D> sub, final boolean reversed) {
+        private void addContribution(final SubHyperplane<Euclidean2D> sub,
+                                     final BSPTree<Euclidean2D> node,
+                                     final Iterable<BSPTree<Euclidean2D>> splitters,
+                                     final boolean reversed) {
             @SuppressWarnings("unchecked")
             final AbstractSubHyperplane<Euclidean2D, Euclidean1D> absSub =
                 (AbstractSubHyperplane<Euclidean2D, Euclidean1D>) sub;
             final Line line      = (Line) sub.getHyperplane();
             final List<Interval> intervals = ((IntervalsSet) absSub.getRemainingRegion()).asList();
             for (final Interval i : intervals) {
-                final Vector2D start = Double.isInfinite(i.getInf()) ?
-                                      null : (Vector2D) line.toSpace((Point<Euclidean1D>) new Vector1D(i.getInf()));
-                final Vector2D end   = Double.isInfinite(i.getSup()) ?
-                                      null : (Vector2D) line.toSpace((Point<Euclidean1D>) new Vector1D(i.getSup()));
+
+                // find the 2D points
+                final Vector2D startV = Double.isInfinite(i.getInf()) ?
+                                        null : (Vector2D) line.toSpace((Point<Euclidean1D>) new Vector1D(i.getInf()));
+                final Vector2D endV   = Double.isInfinite(i.getSup()) ?
+                                        null : (Vector2D) line.toSpace((Point<Euclidean1D>) new Vector1D(i.getSup()));
+
+                // recover the connectivity information
+                final BSPTree<Euclidean2D> startN = selectClosest(startV, splitters);
+                final BSPTree<Euclidean2D> endN   = selectClosest(endV, splitters);
+
                 if (reversed) {
-                    sorted.insert(new ComparableSegment(end, start, line.getReverse()));
+                    segments.add(new ConnectableSegment(endV, startV, line.getReverse(),
+                                                        node, endN, startN));
                 } else {
-                    sorted.insert(new ComparableSegment(start, end, line));
+                    segments.add(new ConnectableSegment(startV, endV, line,
+                                                        node, startN, endN));
+                }
+
+            }
+        }
+
+        /** Select the node whose cut sub-hyperplane is closest to specified point.
+         * @param point reference point
+         * @param candidates candidate nodes
+         * @return node closest to point, or null if no node is closer than tolerance
+         */
+        private BSPTree<Euclidean2D> selectClosest(final Vector2D point, final Iterable<BSPTree<Euclidean2D>> candidates) {
+
+            BSPTree<Euclidean2D> selected = null;
+            double min = Double.POSITIVE_INFINITY;
+
+            for (final BSPTree<Euclidean2D> node : candidates) {
+                final double distance = FastMath.abs(node.getCut().getHyperplane().getOffset(point));
+                if (distance < min) {
+                    selected = node;
+                    min      = distance;
                 }
             }
+
+            return min <= tolerance ? selected : null;
+
         }
 
-        /** Get the sorted segments.
-         * @return sorted segments
+        /** Get the segments.
+         * @return built segments
          */
-        public AVLTree<ComparableSegment> getSorted() {
-            return sorted;
+        public List<ConnectableSegment> getSegments() {
+            return segments;
         }
 
     }

http://git-wip-us.apache.org/repos/asf/commons-math/blob/e6aae3a8/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java
----------------------------------------------------------------------
diff --git a/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java b/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java
index c6c68a0..52bbbca 100644
--- a/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java
+++ b/src/test/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSetTest.java
@@ -850,34 +850,42 @@ public class PolygonsSetTest {
         Line[] l = {
             new Line(new Vector2D(0.0, 0.625000007541172),
                      new Vector2D(1.0, 0.625000007541172), 1.0e-10),
-                     new Line(new Vector2D(-0.19204433621902645, 0.0),
-                              new Vector2D(-0.19204433621902645, 1.0), 1.0e-10),
-                              new Line(new Vector2D(-0.40303524786887,  0.4248364535319128),
-                                       new Vector2D(-1.12851149797877, -0.2634107480798909), 1.0e-10),
-                                       new Line(new Vector2D(0.0, 2.0),
-                                                new Vector2D(1.0, 2.0), 1.0e-10)
+            new Line(new Vector2D(-0.19204433621902645, 0.0),
+                     new Vector2D(-0.19204433621902645, 1.0), 1.0e-10),
+            new Line(new Vector2D(-0.40303524786887,  0.4248364535319128),
+                     new Vector2D(-1.12851149797877, -0.2634107480798909), 1.0e-10),
+            new Line(new Vector2D(0.0, 2.0),
+                     new Vector2D(1.0, 2.0), 1.0e-10)
         };
 
         BSPTree<Euclidean2D> node1 =
             new BSPTree<Euclidean2D>(new SubLine(l[0],
-                                          new IntervalsSet(intersectionAbscissa(l[0], l[1]),
-                                                           intersectionAbscissa(l[0], l[2]),
-                                                           1.0e-10)),
-                                                           new BSPTree<Euclidean2D>(Boolean.TRUE), new BSPTree<Euclidean2D>(Boolean.FALSE),
-                                                           null);
+                                                 new IntervalsSet(intersectionAbscissa(l[0], l[1]),
+                                                                  intersectionAbscissa(l[0], l[2]),
+                                                                  1.0e-10)),
+                                     new BSPTree<Euclidean2D>(Boolean.TRUE),
+                                     new BSPTree<Euclidean2D>(Boolean.FALSE),
+                                     null);
         BSPTree<Euclidean2D> node2 =
             new BSPTree<Euclidean2D>(new SubLine(l[1],
-                                          new IntervalsSet(intersectionAbscissa(l[1], l[2]),
-                                                           intersectionAbscissa(l[1], l[3]),
-                                                           1.0e-10)),
-                                                           node1, new BSPTree<Euclidean2D>(Boolean.FALSE), null);
+                                                 new IntervalsSet(intersectionAbscissa(l[1], l[2]),
+                                                                  intersectionAbscissa(l[1], l[3]),
+                                                                  1.0e-10)),
+                                     node1,
+                                     new BSPTree<Euclidean2D>(Boolean.FALSE),
+                                     null);
         BSPTree<Euclidean2D> node3 =
             new BSPTree<Euclidean2D>(new SubLine(l[2],
-                                          new IntervalsSet(intersectionAbscissa(l[2], l[3]),
-                                                           Double.POSITIVE_INFINITY, 1.0e-10)),
-                                                           node2, new BSPTree<Euclidean2D>(Boolean.FALSE), null);
+                                                 new IntervalsSet(intersectionAbscissa(l[2], l[3]),
+                                                 Double.POSITIVE_INFINITY, 1.0e-10)),
+                                     node2,
+                                     new BSPTree<Euclidean2D>(Boolean.FALSE),
+                                     null);
         BSPTree<Euclidean2D> node4 =
-            new BSPTree<Euclidean2D>(l[3].wholeHyperplane(), node3, new BSPTree<Euclidean2D>(Boolean.FALSE), null);
+            new BSPTree<Euclidean2D>(l[3].wholeHyperplane(),
+                                     node3,
+                                     new BSPTree<Euclidean2D>(Boolean.FALSE),
+                                     null);
 
         PolygonsSet set = new PolygonsSet(node4, 1.0e-10);
         Assert.assertEquals(0, set.getVertices().length);
@@ -1062,9 +1070,9 @@ public class PolygonsSetTest {
                 RegionFactory<Euclidean2D>().difference(set1.copySelf(),
                                                         set2.copySelf());
 
-        Vector2D[][] verticies = set.getVertices();
-        Assert.assertTrue(verticies[0][0] != null);
-        Assert.assertEquals(1, verticies.length);
+        Vector2D[][] vertices = set.getVertices();
+        Assert.assertTrue(vertices[0][0] != null);
+        Assert.assertEquals(1, vertices.length);
     }
 
     @Test
@@ -1167,12 +1175,8 @@ public class PolygonsSetTest {
     private SubHyperplane<Euclidean2D> buildHalfLine(Vector2D start, Vector2D end,
                                                      boolean startIsVirtual) {
         Line   line  = new Line(start, end, 1.0e-10);
-        double lower = startIsVirtual
-        ? Double.NEGATIVE_INFINITY
-        : (line.toSubSpace(start)).getX();
-        double upper = startIsVirtual
-        ? (line.toSubSpace(end)).getX()
-        : Double.POSITIVE_INFINITY;
+        double lower = startIsVirtual ? Double.NEGATIVE_INFINITY : (line.toSubSpace(start)).getX();
+        double upper = startIsVirtual ? (line.toSubSpace(end)).getX() : Double.POSITIVE_INFINITY;
         return new SubLine(line, new IntervalsSet(lower, upper, 1.0e-10));
     }
 


[2/8] [math] Reorganized code, avoiding too many internal classes.

Posted by lu...@apache.org.
Reorganized code, avoiding too many internal classes.

Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/26ee4819
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/26ee4819
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/26ee4819

Branch: refs/heads/master
Commit: 26ee48193928f525fe0885cba6b28fdc8eff6af9
Parents: 5f667c0
Author: Luc Maisonobe <lu...@apache.org>
Authored: Sun Nov 30 13:37:43 2014 +0100
Committer: Luc Maisonobe <lu...@apache.org>
Committed: Tue Dec 2 15:24:31 2014 +0100

----------------------------------------------------------------------
 .../geometry/euclidean/twod/PolygonsSet.java    |   2 +-
 .../geometry/partitioning/AbstractRegion.java   | 266 +------------------
 .../geometry/partitioning/BoundaryBuilder.java  |  84 ++++++
 .../geometry/partitioning/Characterization.java | 147 ++++++++++
 .../geometry/partitioning/InsideFinder.java     | 150 +++++++++++
 .../geometry/partitioning/RegionFactory.java    |  47 +++-
 6 files changed, 430 insertions(+), 266 deletions(-)
----------------------------------------------------------------------


http://git-wip-us.apache.org/repos/asf/commons-math/blob/26ee4819/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java
index fe7c463..532d6a0 100644
--- a/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java
+++ b/src/main/java/org/apache/commons/math3/geometry/euclidean/twod/PolygonsSet.java
@@ -784,7 +784,7 @@ public class PolygonsSet extends AbstractRegion<Euclidean2D, Euclidean1D> {
         // is this an open or a closed loop ?
         final boolean open = segment.getStart() == null;
 
-        while ((end != null) && (open || (globalStart.distance((Point<Euclidean2D>) end) > 1.0e-10))) {
+        while ((end != null) && (open || (globalStart.distance((Point<Euclidean2D>) end) > getTolerance()))) {
 
             // search the sub-hyperplane starting where the previous one ended
             AVLTree<ComparableSegment>.Node selectedNode = null;

http://git-wip-us.apache.org/repos/asf/commons-math/blob/26ee4819/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java
index abddb78..3fca476 100644
--- a/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/AbstractRegion.java
@@ -16,18 +16,15 @@
  */
 package org.apache.commons.math3.geometry.partitioning;
 
-import java.lang.reflect.Array;
 import java.util.ArrayList;
 import java.util.Collection;
 import java.util.Comparator;
 import java.util.Iterator;
 import java.util.TreeSet;
 
-import org.apache.commons.math3.exception.MathInternalError;
-import org.apache.commons.math3.geometry.Space;
 import org.apache.commons.math3.geometry.Point;
+import org.apache.commons.math3.geometry.Space;
 import org.apache.commons.math3.geometry.Vector;
-import org.apache.commons.math3.geometry.partitioning.Region.Location;
 
 /** Abstract class for all regions, independently of geometry type or dimension.
 
@@ -358,122 +355,6 @@ public abstract class AbstractRegion<S extends Space, T extends Space> implement
         return tree;
     }
 
-    /** Visitor building boundary shell tree.
-     * <p>
-     * The boundary shell is represented as {@link BoundaryAttribute boundary attributes}
-     * at each internal node.
-     * </p>
-     */
-    private static class BoundaryBuilder<S extends Space> implements BSPTreeVisitor<S> {
-
-        /** {@inheritDoc} */
-        public Order visitOrder(BSPTree<S> node) {
-            return Order.PLUS_MINUS_SUB;
-        }
-
-        /** {@inheritDoc} */
-        public void visitInternalNode(BSPTree<S> node) {
-
-            SubHyperplane<S> plusOutside = null;
-            SubHyperplane<S> plusInside  = null;
-
-            // characterize the cut sub-hyperplane,
-            // first with respect to the plus sub-tree
-            @SuppressWarnings("unchecked")
-            final SubHyperplane<S>[] plusChar = (SubHyperplane<S>[]) Array.newInstance(SubHyperplane.class, 2);
-            characterize(node.getPlus(), node.getCut().copySelf(), plusChar);
-
-            if (plusChar[0] != null && !plusChar[0].isEmpty()) {
-                // plusChar[0] corresponds to a subset of the cut sub-hyperplane known to have
-                // outside cells on its plus side, we want to check if parts of this subset
-                // do have inside cells on their minus side
-                @SuppressWarnings("unchecked")
-                final SubHyperplane<S>[] minusChar = (SubHyperplane<S>[]) Array.newInstance(SubHyperplane.class, 2);
-                characterize(node.getMinus(), plusChar[0], minusChar);
-                if (minusChar[1] != null && !minusChar[1].isEmpty()) {
-                    // this part belongs to the boundary,
-                    // it has the outside on its plus side and the inside on its minus side
-                    plusOutside = minusChar[1];
-                }
-            }
-
-            if (plusChar[1] != null && !plusChar[1].isEmpty()) {
-                // plusChar[1] corresponds to a subset of the cut sub-hyperplane known to have
-                // inside cells on its plus side, we want to check if parts of this subset
-                // do have outside cells on their minus side
-                @SuppressWarnings("unchecked")
-                final SubHyperplane<S>[] minusChar = (SubHyperplane<S>[]) Array.newInstance(SubHyperplane.class, 2);
-                characterize(node.getMinus(), plusChar[1], minusChar);
-                if (minusChar[0] != null && !minusChar[0].isEmpty()) {
-                    // this part belongs to the boundary,
-                    // it has the inside on its plus side and the outside on its minus side
-                    plusInside = minusChar[0];
-                }
-            }
-
-            // set the boundary attribute at non-leaf nodes
-            node.setAttribute(new BoundaryAttribute<S>(plusOutside, plusInside));
-
-        }
-
-        /** {@inheritDoc} */
-        public void visitLeafNode(BSPTree<S> node) {
-        }
-
-        /** Filter the parts of an hyperplane belonging to the boundary.
-         * <p>The filtering consist in splitting the specified
-         * sub-hyperplane into several parts lying in inside and outside
-         * cells of the tree. The principle is to call this method twice for
-         * each cut sub-hyperplane in the tree, once on the plus node and
-         * once on the minus node. The parts that have the same flag
-         * (inside/inside or outside/outside) do not belong to the boundary
-         * while parts that have different flags (inside/outside or
-         * outside/inside) do belong to the boundary.</p>
-         * @param node current BSP tree node
-         * @param sub sub-hyperplane to characterize
-         * @param characterization placeholder where to put the characterized parts
-         */
-        private void characterize(final BSPTree<S> node, final SubHyperplane<S> sub,
-                                  final SubHyperplane<S>[] characterization) {
-            if (node.getCut() == null) {
-                // we have reached a leaf node
-                final boolean inside = (Boolean) node.getAttribute();
-                if (inside) {
-                    if (characterization[1] == null) {
-                        characterization[1] = sub;
-                    } else {
-                        characterization[1] = characterization[1].reunite(sub);
-                    }
-                } else {
-                    if (characterization[0] == null) {
-                        characterization[0] = sub;
-                    } else {
-                        characterization[0] = characterization[0].reunite(sub);
-                    }
-                }
-            } else {
-                final Hyperplane<S> hyperplane = node.getCut().getHyperplane();
-                switch (sub.side(hyperplane)) {
-                case PLUS:
-                    characterize(node.getPlus(), sub, characterization);
-                    break;
-                case MINUS:
-                    characterize(node.getMinus(), sub, characterization);
-                    break;
-                case BOTH:
-                    final SubHyperplane.SplitSubHyperplane<S> split = sub.split(hyperplane);
-                    characterize(node.getPlus(),  split.getPlus(),  characterization);
-                    characterize(node.getMinus(), split.getMinus(), characterization);
-                    break;
-                default:
-                    // this should not happen
-                    throw new MathInternalError();
-                }
-            }
-        }
-
-    }
-
     /** {@inheritDoc} */
     public double getBoundarySize() {
         final BoundarySizeVisitor<S> visitor = new BoundarySizeVisitor<S>();
@@ -525,146 +406,11 @@ public abstract class AbstractRegion<S extends Space, T extends Space> implement
 
     /** {@inheritDoc} */
     public Side side(final Hyperplane<S> hyperplane) {
-        final Sides sides = new Sides();
-        recurseSides(tree, hyperplane.wholeHyperplane(), sides);
-        return sides.plusFound() ?
-              (sides.minusFound() ? Side.BOTH  : Side.PLUS) :
-              (sides.minusFound() ? Side.MINUS : Side.HYPER);
-    }
-
-    /** Search recursively for inside leaf nodes on each side of the given hyperplane.
-
-     * <p>The algorithm used here is directly derived from the one
-     * described in section III (<i>Binary Partitioning of a BSP
-     * Tree</i>) of the Bruce Naylor, John Amanatides and William
-     * Thibault paper <a
-     * href="http://www.cs.yorku.ca/~amana/research/bsptSetOp.pdf">Merging
-     * BSP Trees Yields Polyhedral Set Operations</a> Proc. Siggraph
-     * '90, Computer Graphics 24(4), August 1990, pp 115-124, published
-     * by the Association for Computing Machinery (ACM)..</p>
-
-     * @param node current BSP tree node
-     * @param sub sub-hyperplane
-     * @param sides object holding the sides found
-     */
-    private void recurseSides(final BSPTree<S> node, final SubHyperplane<S> sub, final Sides sides) {
-
-        if (node.getCut() == null) {
-            if ((Boolean) node.getAttribute()) {
-                // this is an inside cell expanding across the hyperplane
-                sides.rememberPlusFound();
-                sides.rememberMinusFound();
-            }
-            return;
-        }
-
-        final Hyperplane<S> hyperplane = node.getCut().getHyperplane();
-        switch (sub.side(hyperplane)) {
-        case PLUS :
-            // the sub-hyperplane is entirely in the plus sub-tree
-            if (node.getCut().side(sub.getHyperplane()) == Side.PLUS) {
-                if (!isEmpty(node.getMinus())) {
-                    sides.rememberPlusFound();
-                }
-            } else {
-                if (!isEmpty(node.getMinus())) {
-                    sides.rememberMinusFound();
-                }
-            }
-            if (!(sides.plusFound() && sides.minusFound())) {
-                recurseSides(node.getPlus(), sub, sides);
-            }
-            break;
-        case MINUS :
-            // the sub-hyperplane is entirely in the minus sub-tree
-            if (node.getCut().side(sub.getHyperplane()) == Side.PLUS) {
-                if (!isEmpty(node.getPlus())) {
-                    sides.rememberPlusFound();
-                }
-            } else {
-                if (!isEmpty(node.getPlus())) {
-                    sides.rememberMinusFound();
-                }
-            }
-            if (!(sides.plusFound() && sides.minusFound())) {
-                recurseSides(node.getMinus(), sub, sides);
-            }
-            break;
-        case BOTH :
-            // the sub-hyperplane extends in both sub-trees
-            final SubHyperplane.SplitSubHyperplane<S> split = sub.split(hyperplane);
-
-            // explore first the plus sub-tree
-            recurseSides(node.getPlus(), split.getPlus(), sides);
-
-            // if needed, explore the minus sub-tree
-            if (!(sides.plusFound() && sides.minusFound())) {
-                recurseSides(node.getMinus(), split.getMinus(), sides);
-            }
-            break;
-        default :
-            // the sub-hyperplane and the cut sub-hyperplane share the same hyperplane
-            if (node.getCut().getHyperplane().sameOrientationAs(sub.getHyperplane())) {
-                if ((node.getPlus().getCut() != null) || ((Boolean) node.getPlus().getAttribute())) {
-                    sides.rememberPlusFound();
-                }
-                if ((node.getMinus().getCut() != null) || ((Boolean) node.getMinus().getAttribute())) {
-                    sides.rememberMinusFound();
-                }
-            } else {
-                if ((node.getPlus().getCut() != null) || ((Boolean) node.getPlus().getAttribute())) {
-                    sides.rememberMinusFound();
-                }
-                if ((node.getMinus().getCut() != null) || ((Boolean) node.getMinus().getAttribute())) {
-                    sides.rememberPlusFound();
-                }
-            }
-        }
-
-    }
-
-    /** Utility class holding the already found sides. */
-    private static final class Sides {
-
-        /** Indicator of inside leaf nodes found on the plus side. */
-        private boolean plusFound;
-
-        /** Indicator of inside leaf nodes found on the plus side. */
-        private boolean minusFound;
-
-        /** Simple constructor.
-         */
-        public Sides() {
-            plusFound  = false;
-            minusFound = false;
-        }
-
-        /** Remember the fact that inside leaf nodes have been found on the plus side.
-         */
-        public void rememberPlusFound() {
-            plusFound = true;
-        }
-
-        /** Check if inside leaf nodes have been found on the plus side.
-         * @return true if inside leaf nodes have been found on the plus side
-         */
-        public boolean plusFound() {
-            return plusFound;
-        }
-
-        /** Remember the fact that inside leaf nodes have been found on the minus side.
-         */
-        public void rememberMinusFound() {
-            minusFound = true;
-        }
-
-        /** Check if inside leaf nodes have been found on the minus side.
-         * @return true if inside leaf nodes have been found on the minus side
-         */
-        public boolean minusFound() {
-            return minusFound;
-        }
-
+        final InsideFinder<S> finder = new InsideFinder<S>(this);
+        finder.recurseSides(tree, hyperplane.wholeHyperplane());
+        return finder.plusFound() ?
+              (finder.minusFound() ? Side.BOTH  : Side.PLUS) :
+              (finder.minusFound() ? Side.MINUS : Side.HYPER);
     }
 
     /** {@inheritDoc} */

http://git-wip-us.apache.org/repos/asf/commons-math/blob/26ee4819/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java
new file mode 100644
index 0000000..80038f1
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/BoundaryBuilder.java
@@ -0,0 +1,84 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math3.geometry.partitioning;
+
+import org.apache.commons.math3.geometry.Space;
+
+/** Visitor building boundary shell tree.
+ * <p>
+ * The boundary shell is represented as {@link BoundaryAttribute boundary attributes}
+ * at each internal node.
+ * </p>
+ * @param <S> Type of the space.
+ * @since 3.4
+ */
+class BoundaryBuilder<S extends Space> implements BSPTreeVisitor<S> {
+
+    /** Simple constructor.
+     */
+    public BoundaryBuilder() {
+    }
+
+    /** {@inheritDoc} */
+    public Order visitOrder(BSPTree<S> node) {
+        return Order.PLUS_MINUS_SUB;
+    }
+
+    /** {@inheritDoc} */
+    public void visitInternalNode(BSPTree<S> node) {
+
+        SubHyperplane<S> plusOutside = null;
+        SubHyperplane<S> plusInside  = null;
+
+        // characterize the cut sub-hyperplane,
+        // first with respect to the plus sub-tree
+        final Characterization<S> plusChar = new Characterization<S>(node.getPlus(), node.getCut().copySelf());
+
+        if (plusChar.touchOutside()) {
+            // plusChar.outsideTouching() corresponds to a subset of the cut sub-hyperplane
+            // known to have outside cells on its plus side, we want to check if parts
+            // of this subset do have inside cells on their minus side
+            final Characterization<S> minusChar = new Characterization<S>(node.getMinus(), plusChar.outsideTouching());
+            if (minusChar.touchInside()) {
+                // this part belongs to the boundary,
+                // it has the outside on its plus side and the inside on its minus side
+                plusOutside = minusChar.insideTouching();
+            }
+        }
+
+        if (plusChar.touchInside()) {
+            // plusChar.insideTouching() corresponds to a subset of the cut sub-hyperplane
+            // known to have inside cells on its plus side, we want to check if parts
+            // of this subset do have outside cells on their minus side
+            final Characterization<S> minusChar = new Characterization<S>(node.getMinus(), plusChar.insideTouching());
+            if (minusChar.touchOutside()) {
+                // this part belongs to the boundary,
+                // it has the inside on its plus side and the outside on its minus side
+                plusInside = minusChar.outsideTouching();
+            }
+        }
+
+        // set the boundary attribute at non-leaf nodes
+        node.setAttribute(new BoundaryAttribute<S>(plusOutside, plusInside));
+
+    }
+
+    /** {@inheritDoc} */
+    public void visitLeafNode(BSPTree<S> node) {
+    }
+
+}

http://git-wip-us.apache.org/repos/asf/commons-math/blob/26ee4819/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.java
new file mode 100644
index 0000000..a67f873
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/Characterization.java
@@ -0,0 +1,147 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math3.geometry.partitioning;
+
+import org.apache.commons.math3.exception.MathInternalError;
+import org.apache.commons.math3.geometry.Space;
+
+/** Cut sub-hyperplanes characterization with respect to inside/outside cells.
+ * @see BoundaryBuilder
+ * @param <S> Type of the space.
+ * @since 3.4
+ */
+class Characterization<S extends Space> {
+
+    /** Part of the cut sub-hyperplane that touch outside cells. */
+    private SubHyperplane<S> outsideTouching;
+
+    /** Part of the cut sub-hyperplane that touch inside cells. */
+    private SubHyperplane<S> insideTouching;
+
+    /** Simple constructor.
+     * <p>Characterization consists in splitting the specified
+     * sub-hyperplane into several parts lying in inside and outside
+     * cells of the tree. The principle is to compute characterization
+     * twice for each cut sub-hyperplane in the tree, once on the plus
+     * node and once on the minus node. The parts that have the same flag
+     * (inside/inside or outside/outside) do not belong to the boundary
+     * while parts that have different flags (inside/outside or
+     * outside/inside) do belong to the boundary.</p>
+     * @param node current BSP tree node
+     * @param sub sub-hyperplane to characterize
+     */
+    public Characterization(final BSPTree<S> node, final SubHyperplane<S> sub) {
+        outsideTouching = null;
+        insideTouching  = null;
+        characterize(node, sub);
+    }
+
+    /** Filter the parts of an hyperplane belonging to the boundary.
+     * <p>The filtering consist in splitting the specified
+     * sub-hyperplane into several parts lying in inside and outside
+     * cells of the tree. The principle is to call this method twice for
+     * each cut sub-hyperplane in the tree, once on the plus node and
+     * once on the minus node. The parts that have the same flag
+     * (inside/inside or outside/outside) do not belong to the boundary
+     * while parts that have different flags (inside/outside or
+     * outside/inside) do belong to the boundary.</p>
+     * @param node current BSP tree node
+     * @param sub sub-hyperplane to characterize
+     */
+    private void characterize(final BSPTree<S> node, final SubHyperplane<S> sub) {
+        if (node.getCut() == null) {
+            // we have reached a leaf node
+            final boolean inside = (Boolean) node.getAttribute();
+            if (inside) {
+                addInsideTouching(sub);
+            } else {
+                addOutsideTouching(sub);
+            }
+        } else {
+            final Hyperplane<S> hyperplane = node.getCut().getHyperplane();
+            switch (sub.side(hyperplane)) {
+            case PLUS:
+                characterize(node.getPlus(), sub);
+                break;
+            case MINUS:
+                characterize(node.getMinus(), sub);
+                break;
+            case BOTH:
+                final SubHyperplane.SplitSubHyperplane<S> split = sub.split(hyperplane);
+                characterize(node.getPlus(),  split.getPlus());
+                characterize(node.getMinus(), split.getMinus());
+                break;
+            default:
+                // this should not happen
+                throw new MathInternalError();
+            }
+        }
+    }
+
+    /** Add a part of the cut sub-hyperplane known to touch an outside cell.
+     * @param sub part of the cut sub-hyperplane known to touch an outside cell
+     */
+    private void addOutsideTouching(final SubHyperplane<S> sub) {
+        if (outsideTouching == null) {
+            outsideTouching = sub;
+        } else {
+            outsideTouching = outsideTouching.reunite(sub);
+        }
+    }
+
+    /** Add a part of the cut sub-hyperplane known to touch an inside cell.
+     * @param sub part of the cut sub-hyperplane known to touch an inside cell
+     */
+    private void addInsideTouching(final SubHyperplane<S> sub) {
+        if (insideTouching == null) {
+            insideTouching = sub;
+        } else {
+            insideTouching = insideTouching.reunite(sub);
+        }
+    }
+
+    /** Check if the cut sub-hyperplane touches outside cells.
+     * @return true if the cut sub-hyperplane touches outside cells
+     */
+    public boolean touchOutside() {
+        return outsideTouching != null && !outsideTouching.isEmpty();
+    }
+
+    /** Get all the parts of the cut sub-hyperplane known to touch outside cells.
+     * @return parts of the cut sub-hyperplane known to touch outside cells
+     * (may be null or empty)
+     */
+    public SubHyperplane<S> outsideTouching() {
+        return outsideTouching;
+    }
+
+    /** Check if the cut sub-hyperplane touches inside cells.
+     * @return true if the cut sub-hyperplane touches inside cells
+     */
+    public boolean touchInside() {
+        return insideTouching != null && !insideTouching.isEmpty();
+    }
+
+    /** Get all the parts of the cut sub-hyperplane known to touch inside cells.
+     * @return parts of the cut sub-hyperplane known to touch inside cells
+     * (may be null or empty)
+     */
+    public SubHyperplane<S> insideTouching() {
+        return insideTouching;
+    }
+
+}

http://git-wip-us.apache.org/repos/asf/commons-math/blob/26ee4819/src/main/java/org/apache/commons/math3/geometry/partitioning/InsideFinder.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/InsideFinder.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/InsideFinder.java
new file mode 100644
index 0000000..2b0b405
--- /dev/null
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/InsideFinder.java
@@ -0,0 +1,150 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements.  See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License.  You may obtain a copy of the License at
+ *
+ *      http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.commons.math3.geometry.partitioning;
+
+import org.apache.commons.math3.geometry.Space;
+
+/** Utility class checking if inside nodes can be found
+ * on the plus and minus sides of an hyperplane.
+ * @param <S> Type of the space.
+ * @since 3.4
+ */
+class InsideFinder<S extends Space> {
+
+    /** Region on which to operate. */
+    private final Region<S> region;
+
+    /** Indicator of inside leaf nodes found on the plus side. */
+    private boolean plusFound;
+
+    /** Indicator of inside leaf nodes found on the plus side. */
+    private boolean minusFound;
+
+    /** Simple constructor.
+     * @param region region on which to operate
+     */
+    public InsideFinder(final Region<S> region) {
+        this.region = region;
+        plusFound  = false;
+        minusFound = false;
+    }
+
+    /** Search recursively for inside leaf nodes on each side of the given hyperplane.
+
+     * <p>The algorithm used here is directly derived from the one
+     * described in section III (<i>Binary Partitioning of a BSP
+     * Tree</i>) of the Bruce Naylor, John Amanatides and William
+     * Thibault paper <a
+     * href="http://www.cs.yorku.ca/~amana/research/bsptSetOp.pdf">Merging
+     * BSP Trees Yields Polyhedral Set Operations</a> Proc. Siggraph
+     * '90, Computer Graphics 24(4), August 1990, pp 115-124, published
+     * by the Association for Computing Machinery (ACM)..</p>
+
+     * @param node current BSP tree node
+     * @param sub sub-hyperplane
+     */
+    public void recurseSides(final BSPTree<S> node, final SubHyperplane<S> sub) {
+
+        if (node.getCut() == null) {
+            if ((Boolean) node.getAttribute()) {
+                // this is an inside cell expanding across the hyperplane
+                plusFound  = true;
+                minusFound = true;
+            }
+            return;
+        }
+
+        final Hyperplane<S> hyperplane = node.getCut().getHyperplane();
+        switch (sub.side(hyperplane)) {
+        case PLUS :
+            // the sub-hyperplane is entirely in the plus sub-tree
+            if (node.getCut().side(sub.getHyperplane()) == Side.PLUS) {
+                if (!region.isEmpty(node.getMinus())) {
+                    plusFound  = true;
+                }
+            } else {
+                if (!region.isEmpty(node.getMinus())) {
+                    minusFound = true;
+                }
+            }
+            if (!(plusFound && minusFound)) {
+                recurseSides(node.getPlus(), sub);
+            }
+            break;
+        case MINUS :
+            // the sub-hyperplane is entirely in the minus sub-tree
+            if (node.getCut().side(sub.getHyperplane()) == Side.PLUS) {
+                if (!region.isEmpty(node.getPlus())) {
+                    plusFound  = true;
+                }
+            } else {
+                if (!region.isEmpty(node.getPlus())) {
+                    minusFound = true;
+                }
+            }
+            if (!(plusFound && minusFound)) {
+                recurseSides(node.getMinus(), sub);
+            }
+            break;
+        case BOTH :
+            // the sub-hyperplane extends in both sub-trees
+            final SubHyperplane.SplitSubHyperplane<S> split = sub.split(hyperplane);
+
+            // explore first the plus sub-tree
+            recurseSides(node.getPlus(), split.getPlus());
+
+            // if needed, explore the minus sub-tree
+            if (!(plusFound && minusFound)) {
+                recurseSides(node.getMinus(), split.getMinus());
+            }
+            break;
+        default :
+            // the sub-hyperplane and the cut sub-hyperplane share the same hyperplane
+            if (node.getCut().getHyperplane().sameOrientationAs(sub.getHyperplane())) {
+                if ((node.getPlus().getCut() != null) || ((Boolean) node.getPlus().getAttribute())) {
+                    plusFound  = true;
+                }
+                if ((node.getMinus().getCut() != null) || ((Boolean) node.getMinus().getAttribute())) {
+                    minusFound = true;
+                }
+            } else {
+                if ((node.getPlus().getCut() != null) || ((Boolean) node.getPlus().getAttribute())) {
+                    minusFound = true;
+                }
+                if ((node.getMinus().getCut() != null) || ((Boolean) node.getMinus().getAttribute())) {
+                    plusFound  = true;
+                }
+            }
+        }
+
+    }
+
+    /** Check if inside leaf nodes have been found on the plus side.
+     * @return true if inside leaf nodes have been found on the plus side
+     */
+    public boolean plusFound() {
+        return plusFound;
+    }
+
+    /** Check if inside leaf nodes have been found on the minus side.
+     * @return true if inside leaf nodes have been found on the minus side
+     */
+    public boolean minusFound() {
+        return minusFound;
+    }
+
+}

http://git-wip-us.apache.org/repos/asf/commons-math/blob/26ee4819/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java b/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
index ae9c3dd..47b28bd 100644
--- a/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
+++ b/src/main/java/org/apache/commons/math3/geometry/partitioning/RegionFactory.java
@@ -16,8 +16,12 @@
  */
 package org.apache.commons.math3.geometry.partitioning;
 
+import org.apache.commons.math3.geometry.Point;
 import org.apache.commons.math3.geometry.Space;
+import org.apache.commons.math3.geometry.euclidean.twod.PolygonsSet;
+import org.apache.commons.math3.geometry.euclidean.twod.Vector2D;
 import org.apache.commons.math3.geometry.partitioning.BSPTree.VanishingCutHandler;
+import org.apache.commons.math3.geometry.partitioning.Region.Location;
 
 /** This class is a factory for {@link Region}.
 
@@ -115,7 +119,7 @@ public class RegionFactory<S extends Space> {
      */
     public Region<S> difference(final Region<S> region1, final Region<S> region2) {
         final BSPTree<S> tree =
-            region1.getTree(false).merge(region2.getTree(false), new DifferenceMerger());
+            region1.getTree(false).merge(region2.getTree(false), new DifferenceMerger(region1, region2));
         tree.visit(nodeCleaner);
         return region1.buildNew(tree);
     }
@@ -206,7 +210,23 @@ public class RegionFactory<S extends Space> {
     }
 
     /** BSP tree leaf merger computing difference of two regions. */
-    private class DifferenceMerger implements BSPTree.LeafMerger<S> {
+    private class DifferenceMerger implements BSPTree.LeafMerger<S>, VanishingCutHandler<S> {
+
+        /** Region to subtract from. */
+        private final Region<S> region1;
+
+        /** Region to subtract. */
+        private final Region<S> region2;
+
+        /** Simple constructor.
+         * @param region1 region to subtract from
+         * @param region2 region to subtract
+         */
+        public DifferenceMerger(final Region<S> region1, final Region<S> region2) {
+            this.region1 = region1.copySelf();
+            this.region2 = region2.copySelf();
+        }
+
         /** {@inheritDoc} */
         public BSPTree<S> merge(final BSPTree<S> leaf, final BSPTree<S> tree,
                                 final BSPTree<S> parentTree, final boolean isPlusChild,
@@ -215,15 +235,32 @@ public class RegionFactory<S extends Space> {
                 // the leaf node represents an inside cell
                 final BSPTree<S> argTree =
                     recurseComplement(leafFromInstance ? tree : leaf);
-                argTree.insertInTree(parentTree, isPlusChild, new VanishingToLeaf(true));
+                argTree.insertInTree(parentTree, isPlusChild, this);
                 return argTree;
             }
             // the leaf node represents an outside cell
             final BSPTree<S> instanceTree =
                 leafFromInstance ? leaf : tree;
-            instanceTree.insertInTree(parentTree, isPlusChild, new VanishingToLeaf(false));
+            instanceTree.insertInTree(parentTree, isPlusChild, this);
             return instanceTree;
         }
+
+        /** {@inheritDoc} */
+        public BSPTree<S> fixNode(final BSPTree<S> node) {
+            // get a representative point in the degenerate cell
+            final BSPTree<S> cell = node.pruneAroundConvexCell(Boolean.TRUE, Boolean.FALSE, null);
+            final Region<S> r = region1.buildNew(cell);
+            for (Vector2D[] loop : ((PolygonsSet) r).getVertices()) {
+                System.out.format(java.util.Locale.US, "%n");
+                for (Vector2D v : loop) {
+                    System.out.format(java.util.Locale.US, "%14.10f %14.10f%n", v.getX(), v.getY());
+                }
+            }
+            final Point<S> p = r.getBarycenter();
+            return new BSPTree<S>(region1.checkPoint(p) == Location.INSIDE &&
+                                  region2.checkPoint(p) == Location.OUTSIDE);
+        }
+
     }
 
     /** Visitor removing internal nodes attributes. */
@@ -252,7 +289,7 @@ public class RegionFactory<S extends Space> {
         private final boolean inside;
 
         /** Simple constructor.
-         * @param inside inside/outside indocator to use for ambiguous nodes
+         * @param inside inside/outside indicator to use for ambiguous nodes
          */
         public VanishingToLeaf(final boolean inside) {
             this.inside = inside;