001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.geometry.core.partitioning;
018
019import java.util.ArrayList;
020import java.util.Collections;
021import java.util.Iterator;
022import java.util.List;
023import java.util.function.Function;
024
025import org.apache.commons.geometry.core.Point;
026import org.apache.commons.geometry.core.RegionLocation;
027import org.apache.commons.geometry.core.Transform;
028
029/** Base class for convex hyperplane-bounded regions. This class provides generic implementations of many
030 * algorithms related to convex regions.
031 * @param <P> Point implementation type
032 * @param <S> Hyperplane convex subset implementation type
033 */
034public abstract class AbstractConvexHyperplaneBoundedRegion<P extends Point<P>, S extends HyperplaneConvexSubset<P>>
035    implements HyperplaneBoundedRegion<P> {
036    /** List of boundaries for the region. */
037    private final List<S> boundaries;
038
039    /** Simple constructor. Callers are responsible for ensuring that the given list of boundaries
040     * define a convex region. No validation is performed.
041     * @param boundaries the boundaries of the convex region
042     */
043    protected AbstractConvexHyperplaneBoundedRegion(final List<S> boundaries) {
044        this.boundaries = Collections.unmodifiableList(boundaries);
045    }
046
047    /** Get the boundaries of the convex region. The exact ordering of the boundaries
048     * is not guaranteed.
049     * @return the boundaries of the convex region
050     */
051    public List<S> getBoundaries() {
052        return boundaries;
053    }
054
055    /** {@inheritDoc} */
056    @Override
057    public boolean isFull() {
058        // no boundaries => no outside
059        return boundaries.isEmpty();
060    }
061
062    /** {@inheritDoc}
063     *
064     * <p>This method always returns false.</p>
065     */
066    @Override
067    public boolean isEmpty() {
068        return false;
069    }
070
071    /** {@inheritDoc} */
072    @Override
073    public double getBoundarySize() {
074        double sum = 0.0;
075        for (final S boundary : boundaries) {
076            sum += boundary.getSize();
077        }
078
079        return sum;
080    }
081
082    /** {@inheritDoc} */
083    @Override
084    public RegionLocation classify(final P pt) {
085        boolean isOn = false;
086
087        HyperplaneLocation loc;
088        for (final S boundary : boundaries) {
089            loc = boundary.getHyperplane().classify(pt);
090
091            if (loc == HyperplaneLocation.PLUS) {
092                return RegionLocation.OUTSIDE;
093            } else if (loc == HyperplaneLocation.ON) {
094                isOn = true;
095            }
096        }
097
098        return isOn ? RegionLocation.BOUNDARY : RegionLocation.INSIDE;
099    }
100
101    /** {@inheritDoc} */
102    @Override
103    public P project(final P pt) {
104
105        P projected;
106        double dist;
107
108        P closestPt = null;
109        double closestDist = Double.POSITIVE_INFINITY;
110
111        for (final S boundary : boundaries) {
112            projected = boundary.closest(pt);
113            dist = pt.distance(projected);
114
115            if (projected != null && (closestPt == null || dist < closestDist)) {
116                closestPt = projected;
117                closestDist = dist;
118            }
119        }
120
121        return closestPt;
122    }
123
124    /** Trim the given hyperplane subset to the portion contained inside this instance.
125     * @param sub hyperplane subset to trim. Null is returned if the subset does not intersect the instance.
126     * @return portion of the argument that lies entirely inside the region represented by
127     *      this instance, or null if it does not intersect.
128     */
129    public HyperplaneConvexSubset<P> trim(final HyperplaneConvexSubset<P> sub) {
130        HyperplaneConvexSubset<P> remaining = sub;
131        for (final S boundary : boundaries) {
132            remaining = remaining.split(boundary.getHyperplane()).getMinus();
133            if (remaining == null) {
134                break;
135            }
136        }
137
138        return remaining;
139    }
140
141    /** {@inheritDoc} */
142    @Override
143    public String toString() {
144        final StringBuilder sb = new StringBuilder();
145        sb.append(this.getClass().getSimpleName())
146            .append("[boundaries= ")
147            .append(boundaries)
148            .append(']');
149
150        return sb.toString();
151    }
152
153    /** Generic, internal transform method. Subclasses should use this to implement their own transform methods.
154     * @param transform the transform to apply to the instance
155     * @param thisInstance a reference to the current instance; this is passed as
156     *      an argument in order to allow it to be a generic type
157     * @param boundaryType the type used for the boundary hyperplane subsets
158     * @param factory function used to create new convex region instances
159     * @param <R> Region implementation type
160     * @return the result of the transform operation
161     */
162    protected <R extends AbstractConvexHyperplaneBoundedRegion<P, S>> R transformInternal(
163            final Transform<P> transform, final R thisInstance, final Class<S> boundaryType,
164            final Function<? super List<S>, R> factory) {
165
166        if (isFull()) {
167            return thisInstance;
168        }
169
170        final List<S> origBoundaries = getBoundaries();
171
172        final int size = origBoundaries.size();
173        final List<S> tBoundaries = new ArrayList<>(size);
174
175        // determine if the hyperplanes should be reversed
176        final S boundary = origBoundaries.get(0);
177        HyperplaneConvexSubset<P> tBoundary = boundary.transform(transform);
178
179        final boolean reverseDirection = swapsInsideOutside(transform);
180
181        // transform all of the segments
182        if (reverseDirection) {
183            tBoundary = tBoundary.reverse();
184        }
185        tBoundaries.add(boundaryType.cast(tBoundary));
186
187        for (int i = 1; i < origBoundaries.size(); ++i) {
188            tBoundary = origBoundaries.get(i).transform(transform);
189
190            if (reverseDirection) {
191                tBoundary = tBoundary.reverse();
192            }
193
194            tBoundaries.add(boundaryType.cast(tBoundary));
195        }
196
197        return factory.apply(tBoundaries);
198    }
199
200    /** Return true if the given transform swaps the inside and outside of
201     * the region.
202     *
203     * <p>The default behavior of this method is to return true if the transform
204     * does not preserve spatial orientation (ie, {@link Transform#preservesOrientation()}
205     * is false). Subclasses may need to override this method to implement the correct
206     * behavior for their space and dimension.</p>
207     * @param transform transform to check
208     * @return true if the given transform swaps the interior and exterior of
209     *      the region
210     */
211    protected boolean swapsInsideOutside(final Transform<P> transform) {
212        return !transform.preservesOrientation();
213    }
214
215    /** Generic, internal split method. Subclasses should call this from their {@link #split(Hyperplane)} methods.
216     * @param splitter splitting hyperplane
217     * @param thisInstance a reference to the current instance; this is passed as
218     *      an argument in order to allow it to be a generic type
219     * @param boundaryType the type used for the boundary hyperplane subsets
220     * @param factory function used to create new convex region instances
221     * @param <R> Region implementation type
222     * @return the result of the split operation
223     */
224    protected <R extends AbstractConvexHyperplaneBoundedRegion<P, S>> Split<R> splitInternal(
225            final Hyperplane<P> splitter, final R thisInstance, final Class<S> boundaryType,
226            final Function<List<S>, R> factory) {
227
228        return isFull() ?
229                splitInternalFull(splitter, boundaryType, factory) :
230                splitInternalNonFull(splitter, thisInstance, boundaryType, factory);
231    }
232
233    /** Internal split method for use with full regions, i.e. regions that cover the entire space.
234     * @param splitter splitting hyperplane
235     * @param boundaryType the type used for the boundary hyperplane subsets
236     * @param factory function used to create new convex region instances
237     * @param <R> Region implementation type
238     * @return the result of the split operation
239     */
240    private <R extends AbstractConvexHyperplaneBoundedRegion<P, S>> Split<R> splitInternalFull(
241            final Hyperplane<P> splitter, final Class<S> boundaryType, final Function<? super List<S>, R> factory) {
242
243        final R minus = factory.apply(Collections.singletonList(boundaryType.cast(splitter.span())));
244        final R plus = factory.apply(Collections.singletonList(boundaryType.cast(splitter.reverse().span())));
245
246        return new Split<>(minus, plus);
247    }
248
249    /** Internal split method for use with non-full regions, i.e. regions that do not cover the entire space.
250     * @param splitter splitting hyperplane
251     * @param thisInstance a reference to the current instance; this is passed as
252     *      an argument in order to allow it to be a generic type
253     * @param boundaryType the type used for the boundary hyperplane subsets
254     * @param factory function used to create new convex region instances
255     * @param <R> Region implementation type
256     * @return the result of the split operation
257     */
258    private <R extends AbstractConvexHyperplaneBoundedRegion<P, S>> Split<R> splitInternalNonFull(
259            final Hyperplane<P> splitter, final R thisInstance, final Class<S> boundaryType,
260            final Function<? super List<S>, R> factory) {
261
262        final HyperplaneConvexSubset<P> trimmedSplitter = trim(splitter.span());
263
264        if (trimmedSplitter == null) {
265            // The splitter lies entirely outside of the region; we need
266            // to determine whether we lie on the plus or minus side of the splitter.
267
268            final SplitLocation regionLoc = determineRegionPlusMinusLocation(splitter);
269            return regionLoc == SplitLocation.MINUS ?
270                    new Split<>(thisInstance, null) :
271                    new Split<>(null, thisInstance);
272        }
273
274        // the splitter passes through the region; split the other region boundaries
275        // by the splitter
276        final ArrayList<S> minusBoundaries = new ArrayList<>();
277        final ArrayList<S> plusBoundaries = new ArrayList<>();
278
279        splitBoundaries(splitter, boundaryType, minusBoundaries, plusBoundaries);
280
281        // if the splitter was trimmed by the region boundaries, double-check that the split boundaries
282        // actually lie on both sides of the splitter; this is another case where floating point errors
283        // can cause a discrepancy between the results of splitting the splitter by the boundaries and
284        // splitting the boundaries by the splitter
285        if (!trimmedSplitter.isFull()) {
286            if (minusBoundaries.isEmpty()) {
287                if (plusBoundaries.isEmpty()) {
288                    return new Split<>(null, null);
289                }
290                return new Split<>(null, thisInstance);
291            } else if (plusBoundaries.isEmpty()) {
292                return new Split<>(thisInstance, null);
293            }
294        }
295
296        // we have a consistent region split; create the new plus and minus regions
297        minusBoundaries.add(boundaryType.cast(trimmedSplitter));
298        plusBoundaries.add(boundaryType.cast(trimmedSplitter.reverse()));
299
300        minusBoundaries.trimToSize();
301        plusBoundaries.trimToSize();
302
303        return new Split<>(factory.apply(minusBoundaries), factory.apply(plusBoundaries));
304    }
305
306    /** Determine whether the region lies on the plus or minus side of the given splitter. It is assumed
307     * that (1) the region is not full, and (2) the given splitter does not pass through the region.
308     *
309     * <p>In theory, this is a very simple operation: one need only test a single region boundary
310     * to see if it lies on the plus or minus side of the splitter. In practice, however, accumulated
311     * floating point errors can cause discrepancies between the splitting operations, causing
312     * boundaries to be classified as lying on both sides of the splitter when they should only lie on one.
313     * Therefore, this method examines as many boundaries as needed in order to determine the best response.
314     * The algorithm proceeds as follows:
315     * <ol>
316     *  <li>If any boundary lies completely on the minus or plus side of the splitter, then
317     *      {@link SplitLocation#MINUS MINUS} or {@link SplitLocation#PLUS PLUS} is returned, respectively.</li>
318     *  <li>If any boundary is coincident with the splitter ({@link SplitLocation#NEITHER NEITHER}), then
319     *      {@link SplitLocation#MINUS MINUS} is returned if the boundary hyperplane has the same orientation
320     *      as the splitter, otherwise {@link SplitLocation#PLUS PLUS}.</li>
321     *  <li>If no boundaries match the above conditions, then the sizes of the split boundaries are compared. If
322     *      the sum of the sizes of the boundaries on the minus side is greater than the sum of the sizes of
323     *      the boundaries on the plus size, then {@link SplitLocation#MINUS MINUS} is returned. Otherwise,
324     *      {@link SplitLocation#PLUS PLUS} is returned.
325     * </ol>
326     * @param splitter splitter to classify the region against; the splitter is assumed to lie
327     *      completely outside of the region
328     * @return {@link SplitLocation#MINUS} if the region lies on the minus side of the splitter and
329     *      {@link SplitLocation#PLUS} if the region lies on the plus side of the splitter
330     */
331    private SplitLocation determineRegionPlusMinusLocation(final Hyperplane<P> splitter) {
332        double minusSize = 0;
333        double plusSize = 0;
334
335        Split<? extends HyperplaneConvexSubset<P>> split;
336        SplitLocation loc;
337
338        for (final S boundary : boundaries) {
339            split = boundary.split(splitter);
340            loc = split.getLocation();
341
342            if (loc == SplitLocation.MINUS || loc == SplitLocation.PLUS) {
343                return loc;
344            } else if (loc == SplitLocation.NEITHER) {
345                return splitter.similarOrientation(boundary.getHyperplane()) ?
346                        SplitLocation.MINUS :
347                        SplitLocation.PLUS;
348            } else {
349                minusSize += split.getMinus().getSize();
350                plusSize += split.getPlus().getSize();
351            }
352        }
353
354        return minusSize > plusSize ? SplitLocation.MINUS : SplitLocation.PLUS;
355    }
356
357    /** Split the boundaries of the region by the given hyperplane, adding the split parts into the
358     * corresponding lists.
359     * @param splitter splitting hyperplane
360     * @param boundaryType the type used for the boundary hyperplane subsets
361     * @param minusBoundaries list that will contain the portions of the boundaries on the minus side
362     *      of the splitting hyperplane
363     * @param plusBoundaries list that will contain the portions of the boundaries on the plus side of
364     *      the splitting hyperplane
365     */
366    private void splitBoundaries(final Hyperplane<P> splitter, final Class<S> boundaryType,
367            final List<S> minusBoundaries, final List<S> plusBoundaries) {
368
369        Split<? extends HyperplaneConvexSubset<P>> split;
370        HyperplaneConvexSubset<P> minusBoundary;
371        HyperplaneConvexSubset<P> plusBoundary;
372
373        for (final S boundary : boundaries) {
374            split = boundary.split(splitter);
375
376            minusBoundary = split.getMinus();
377            plusBoundary = split.getPlus();
378
379            if (minusBoundary != null) {
380                minusBoundaries.add(boundaryType.cast(minusBoundary));
381            }
382
383            if (plusBoundary != null) {
384                plusBoundaries.add(boundaryType.cast(plusBoundary));
385            }
386        }
387    }
388
389    /** Internal class encapsulating the logic for building convex region boundaries from collections of hyperplanes.
390     * @param <P> Point implementation type
391     * @param <S> Hyperplane convex subset implementation type
392     */
393    protected static class ConvexRegionBoundaryBuilder<P extends Point<P>, S extends HyperplaneConvexSubset<P>> {
394
395        /** Hyperplane convex subset implementation type. */
396        private final Class<S> subsetType;
397
398        /** Construct a new instance for building convex region boundaries with the given hyperplane
399         * convex subset implementation type.
400         * @param subsetType Hyperplane convex subset implementation type
401         */
402        public ConvexRegionBoundaryBuilder(final Class<S> subsetType) {
403            this.subsetType = subsetType;
404        }
405
406        /** Compute a list of hyperplane convex subsets representing the boundaries of the convex region
407         * bounded by the given collection of hyperplanes.
408         * @param bounds hyperplanes defining the convex region
409         * @return a list of hyperplane convex subsets representing the boundaries of the convex region
410         * @throws IllegalArgumentException if the given hyperplanes do not form a convex region
411         */
412        public List<S> build(final Iterable<? extends Hyperplane<P>> bounds) {
413
414            final List<S> boundaries = new ArrayList<>();
415
416            // cut each hyperplane by every other hyperplane in order to get the region boundaries
417            int boundIdx = -1;
418            HyperplaneConvexSubset<P> boundary;
419
420            for (final Hyperplane<P> currentBound : bounds) {
421                ++boundIdx;
422
423                boundary = splitBound(currentBound, bounds, boundIdx);
424                if (boundary != null) {
425                    boundaries.add(subsetType.cast(boundary));
426                }
427            }
428
429            if (boundIdx > 0 && boundaries.isEmpty()) {
430                // nothing was added
431                throw nonConvexException(bounds);
432            }
433
434            return boundaries;
435        }
436
437        /** Split the given bounding hyperplane by all of the other hyperplanes in the given collection, returning the
438         * remaining hyperplane subset.
439         * @param currentBound the bound to split; this value is assumed to have come from {@code bounds}
440         * @param bounds collection of bounds to use to split {@code currentBound}
441         * @param currentBoundIdx the index of {@code currentBound} in {@code bounds}
442         * @return the part of {@code currentBound}'s hyperplane subset that lies on the minus side of all of the
443         *      splitting hyperplanes
444         * @throws IllegalArgumentException if the hyperplanes do not form a convex region
445         */
446        private HyperplaneConvexSubset<P> splitBound(final Hyperplane<P> currentBound,
447                final Iterable<? extends Hyperplane<P>> bounds, final int currentBoundIdx) {
448
449            HyperplaneConvexSubset<P> boundary = currentBound.span();
450
451            final Iterator<? extends Hyperplane<P>> boundsIt = bounds.iterator();
452
453            Hyperplane<P> splitter;
454            int splitterIdx = -1;
455
456            while (boundsIt.hasNext() && boundary != null) {
457                splitter = boundsIt.next();
458                ++splitterIdx;
459
460                if (currentBound == splitter) {
461                    // do not split the bound with itself
462
463                    if (currentBoundIdx > splitterIdx) {
464                        // this hyperplane is duplicated in the list; skip all but the
465                        // first insertion of its hyperplane subset
466                        return null;
467                    }
468                } else {
469                    // split the boundary
470                    final Split<? extends HyperplaneConvexSubset<P>> split = boundary.split(splitter);
471
472                    if (split.getLocation() != SplitLocation.NEITHER) {
473                        // retain the minus portion of the split
474                        boundary = split.getMinus();
475                    } else if (!currentBound.similarOrientation(splitter)) {
476                        // two or more splitters are coincident and have opposite
477                        // orientations, meaning that no area is on the minus side
478                        // of both
479                        throw nonConvexException(bounds);
480                    } else if (currentBoundIdx > splitterIdx) {
481                        // two or more hyperplanes are equivalent; only use the boundary
482                        // from the first one and return null for this one
483                        return null;
484                    }
485                }
486            }
487
488            return boundary;
489        }
490
491        /** Return an exception indicating that the given collection of hyperplanes do not produce a convex region.
492         * @param bounds collection of hyperplanes
493         * @return an exception indicating that the given collection of hyperplanes do not produce a convex region
494         */
495        private IllegalArgumentException nonConvexException(final Iterable<? extends Hyperplane<P>> bounds) {
496            return new IllegalArgumentException("Bounding hyperplanes do not produce a convex region: " + bounds);
497        }
498    }
499}