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.euclidean.threed;
018
019import java.util.Objects;
020
021import org.apache.commons.geometry.core.Transform;
022import org.apache.commons.geometry.core.partitioning.AbstractHyperplane;
023import org.apache.commons.geometry.core.partitioning.Hyperplane;
024import org.apache.commons.geometry.euclidean.threed.line.Line3D;
025import org.apache.commons.geometry.euclidean.threed.line.Lines3D;
026import org.apache.commons.geometry.euclidean.threed.rotation.QuaternionRotation;
027import org.apache.commons.geometry.euclidean.twod.ConvexArea;
028import org.apache.commons.numbers.core.Precision;
029
030/** Class representing a plane in 3 dimensional Euclidean space. Each plane is defined by a
031 * {@link #getNormal() normal} and an {@link #getOriginOffset() origin offset}. If \(\vec{n}\) is the plane normal,
032 * \(d\) is the origin offset, and \(p\) and \(q\) are any points in the plane, then the following are true:
033 * <ul>
034 *  <li>\(\lVert \vec{n} \rVert\) = 1</li>
035 *  <li>\(\vec{n} \cdot (p - q) = 0\)</li>
036 *  <li>\(d = - (\vec{n} \cdot q)\)</li>
037 *  </ul>
038 *  In other words, the normal is a unit vector such that the dot product of the normal and the difference of
039 *  any two points in the plane is always equal to \(0\). Similarly, the {@code origin offset} is equal to the
040 *  negation of the dot product of the normal and any point in the plane. The projection of the origin onto the
041 *  plane (given by {@link #getOrigin()}), is computed as \(-d \vec{n}\).
042 *
043 * <p>Instances of this class are guaranteed to be immutable.</p>
044 * @see Planes
045 */
046public class Plane extends AbstractHyperplane<Vector3D> {
047
048    /** Plane normal. */
049    private final Vector3D.Unit normal;
050
051    /** Offset of the origin with respect to the plane. */
052    private final double originOffset;
053
054    /** Construct a plane from its component parts.
055     * @param normal unit normal vector
056     * @param originOffset offset of the origin with respect to the plane
057     * @param precision precision context used to compare floating point values
058     */
059    Plane(final Vector3D.Unit normal, final double originOffset,
060          final Precision.DoubleEquivalence precision) {
061
062        super(precision);
063
064        this.normal = normal;
065        this.originOffset = originOffset;
066    }
067
068    /** Get the orthogonal projection of the 3D-space origin in the plane.
069     * @return the origin point of the plane frame (point closest to the 3D-space
070     *         origin)
071     */
072    public Vector3D getOrigin() {
073        return normal.multiply(-originOffset);
074    }
075
076    /** Get the offset of the spatial origin ({@code 0, 0, 0}) with respect to the plane.
077     * @return the offset of the origin with respect to the plane.
078     */
079    public double getOriginOffset() {
080        return originOffset;
081    }
082
083    /** Get the plane normal vector.
084     * @return plane normal vector
085     */
086    public Vector3D.Unit getNormal() {
087        return normal;
088    }
089
090    /** Return an {@link EmbeddingPlane} instance suitable for embedding 2D geometric objects
091     * into this plane. Returned instances are guaranteed to be equal between invocations.
092     * @return a plane instance suitable for embedding 2D subspaces
093     */
094    public EmbeddingPlane getEmbedding() {
095        final Vector3D.Unit u = normal.orthogonal();
096        final Vector3D.Unit v = normal.cross(u).normalize();
097
098        return new EmbeddingPlane(u, v, normal, originOffset, getPrecision());
099    }
100
101    /** {@inheritDoc} */
102    @Override
103    public double offset(final Vector3D point) {
104        return point.dot(normal) + originOffset;
105    }
106
107    /** Get the offset (oriented distance) of the given line with respect to the plane. The value
108     * closest to zero is returned, which will always be zero if the line is not parallel to the plane.
109     * @param line line to calculate the offset of
110     * @return the offset of the line with respect to the plane or 0.0 if the line
111     *      is not parallel to the plane.
112     */
113    public double offset(final Line3D line) {
114        if (!isParallel(line)) {
115            return 0.0;
116        }
117        return offset(line.getOrigin());
118    }
119
120    /** Get the offset (oriented distance) of the given plane with respect to this instance. The value
121     * closest to zero is returned, which will always be zero if the planes are not parallel.
122     * @param plane plane to calculate the offset of
123     * @return the offset of the plane with respect to this instance or 0.0 if the planes
124     *      are not parallel.
125     */
126    public double offset(final Plane plane) {
127        if (!isParallel(plane)) {
128            return 0.0;
129        }
130        return originOffset + (similarOrientation(plane) ? -plane.originOffset : plane.originOffset);
131    }
132
133    /** Check if the instance contains a point.
134     * @param p point to check
135     * @return true if p belongs to the plane
136     */
137    @Override
138    public boolean contains(final Vector3D p) {
139        return getPrecision().eqZero(offset(p));
140    }
141
142    /** Check if the instance contains a line.
143     * @param line line to check
144     * @return true if line is contained in this plane
145     */
146    public boolean contains(final Line3D line) {
147        return isParallel(line) && contains(line.getOrigin());
148    }
149
150    /** Check if the instance contains another plane. Planes are considered similar if they contain
151     * the same points. This does not mean they are equal since they can have opposite normals.
152     * @param plane plane to which the instance is compared
153     * @return true if the planes are similar
154     */
155    public boolean contains(final Plane plane) {
156        final double angle = normal.angle(plane.normal);
157        final Precision.DoubleEquivalence precision = getPrecision();
158
159        return ((precision.eqZero(angle)) && precision.eq(originOffset, plane.originOffset)) ||
160                ((precision.eq(angle, Math.PI)) && precision.eq(originOffset, -plane.originOffset));
161    }
162
163    /** {@inheritDoc} */
164    @Override
165    public Vector3D project(final Vector3D point) {
166        return getOrigin().add(point.reject(normal));
167    }
168
169    /** Project a 3D line onto the plane.
170     * @param line the line to project
171     * @return the projection of the given line onto the plane.
172     */
173    public Line3D project(final Line3D line) {
174        final Vector3D direction = line.getDirection();
175        final Vector3D projection = normal.multiply(direction.dot(normal) * (1 / normal.normSq()));
176
177        final Vector3D projectedLineDirection = direction.subtract(projection);
178        final Vector3D p1 = project(line.getOrigin());
179        final Vector3D p2 = p1.add(projectedLineDirection);
180
181        return Lines3D.fromPoints(p1, p2, getPrecision());
182    }
183
184    /** {@inheritDoc} */
185    @Override
186    public PlaneConvexSubset span() {
187        return Planes.subsetFromConvexArea(getEmbedding(), ConvexArea.full());
188    }
189
190    /** Check if the line is parallel to the instance.
191     * @param line line to check.
192     * @return true if the line is parallel to the instance, false otherwise.
193     */
194    public boolean isParallel(final Line3D line) {
195        final double dot = normal.dot(line.getDirection());
196
197        return getPrecision().eqZero(dot);
198    }
199
200    /** Check if the plane is parallel to the instance.
201     * @param plane plane to check.
202     * @return true if the plane is parallel to the instance, false otherwise.
203     */
204    public boolean isParallel(final Plane plane) {
205        return getPrecision().eqZero(normal.cross(plane.normal).norm());
206    }
207
208    /** {@inheritDoc} */
209    @Override
210    public boolean similarOrientation(final Hyperplane<Vector3D> other) {
211        return (((Plane) other).normal).dot(normal) > 0;
212    }
213
214    /** Get the intersection of a line with this plane.
215     * @param line line intersecting the instance
216     * @return intersection point between between the line and the instance (null if
217     *         the line is parallel to the instance)
218     */
219    public Vector3D intersection(final Line3D line) {
220        final Vector3D direction = line.getDirection();
221        final double dot = normal.dot(direction);
222
223        if (getPrecision().eqZero(dot)) {
224            return null;
225        }
226
227        final Vector3D point = line.pointAt(0);
228        final double k = -(originOffset + normal.dot(point)) / dot;
229
230        return Vector3D.Sum.of(point)
231                .addScaled(k, direction)
232                .get();
233    }
234
235    /** Get the line formed by the intersection of this instance with the given plane.
236     * The returned line lies in both planes and points in the direction of
237     * the cross product <code>n<sub>1</sub> x n<sub>2</sub></code>, where <code>n<sub>1</sub></code>
238     * is the normal of the current instance and <code>n<sub>2</sub></code> is the normal
239     * of the argument.
240     *
241     * <p>Null is returned if the planes are parallel.</p>
242     *
243     * @param other other plane
244     * @return line at the intersection of the instance and the other plane, or null
245     *      if no such line exists
246     */
247    public Line3D intersection(final Plane other) {
248        final Vector3D direction = normal.cross(other.normal);
249
250        if (getPrecision().eqZero(direction.norm())) {
251            return null;
252        }
253
254        final Vector3D point = intersection(this, other, Planes.fromNormal(direction, getPrecision()));
255
256        return Lines3D.fromPointAndDirection(point, direction, getPrecision());
257    }
258
259    /** Build a new reversed version of this plane, with opposite orientation.
260     * @return a new reversed plane
261     */
262    @Override
263    public Plane reverse() {
264        return new Plane(normal.negate(), -originOffset, getPrecision());
265    }
266
267    /** {@inheritDoc}
268     *
269     * <p>Instances are transformed by selecting 3 representative points from the
270     * plane, transforming them, and constructing a new plane from the transformed points.
271     * Since the normal is not transformed directly, but rather is constructed new from the
272     * transformed points, the relative orientations of points in the plane are preserved,
273     * even for transforms that do not
274     * {@link Transform#preservesOrientation() preserve orientation}. The example below shows
275     * a plane being transformed by a non-orientation-preserving transform. The normal of the
276     * transformed plane retains its counterclockwise relationship to the points in the plane,
277     * in contrast with the normal that is transformed directly by the transform.
278     * </p>
279     * <pre>
280     * // construct a plane from 3 points; the normal will be selected such that the
281     * // points are ordered counterclockwise when looking down the plane normal.
282     * Vector3D p1 = Vector3D.of(0, 0, 0);
283     * Vector3D p2 = Vector3D.of(+1, 0, 0);
284     * Vector3D p3 = Vector3D.of(0, +1, 0);
285     *
286     * Plane plane = Planes.fromPoints(p1, p2, p3, precision); // normal is (0, 0, +1)
287     *
288     * // create a transform that negates all x-values; this transform does not
289     * // preserve orientation, i.e. it will convert a right-handed system into a left-handed
290     * // system and vice versa
291     * AffineTransformMatrix3D transform = AffineTransformMatrix3D.createScale(-1, 1,  1);
292     *
293     * // transform the plane
294     * Plane transformedPlane = plane.transform(transform);
295     *
296     * // the plane normal is oriented such that transformed points are still ordered
297     * // counterclockwise when looking down the plane normal; since the point (1, 0, 0) has
298     * // now become (-1, 0, 0), the normal has flipped to (0, 0, -1)
299     * transformedPlane.getNormal();
300     *
301     * // directly transform the original plane normal; the normal is unchanged by the transform
302     * // since the target space of the transform is left-handed
303     * AffineTransformMatrix3D normalTransform = transform.normalTransform();
304     * Vector3D directlyTransformedNormal = normalTransform.apply(plane.getNormal()); // (0, 0, +1)
305     * </pre>
306     */
307    @Override
308    public Plane transform(final Transform<Vector3D> transform) {
309        // create 3 representation points lying on the plane, transform them,
310        // and use the transformed points to create a new plane
311
312        final Vector3D u = normal.orthogonal();
313        final Vector3D v = normal.cross(u);
314
315        final Vector3D p1 = getOrigin();
316        final Vector3D p2 = p1.add(u);
317        final Vector3D p3 = p1.add(v);
318
319        final Vector3D t1 = transform.apply(p1);
320        final Vector3D t2 = transform.apply(p2);
321        final Vector3D t3 = transform.apply(p3);
322
323        return Planes.fromPoints(t1, t2, t3, getPrecision());
324    }
325
326    /** Translate the plane by the specified amount.
327     * @param translation translation to apply
328     * @return a new plane
329     */
330    public Plane translate(final Vector3D translation) {
331        final Vector3D tOrigin = getOrigin().add(translation);
332
333        return Planes.fromPointAndNormal(tOrigin, normal, getPrecision());
334    }
335
336    /** Rotate the plane around the specified point.
337     * @param center rotation center
338     * @param rotation 3-dimensional rotation
339     * @return a new plane
340     */
341    public Plane rotate(final Vector3D center, final QuaternionRotation rotation) {
342        final Vector3D delta = getOrigin().subtract(center);
343        final Vector3D tOrigin = center.add(rotation.apply(delta));
344
345        // we can directly apply the rotation to the normal since it will transform
346        // it properly (there is no translation or scaling involved)
347        final Vector3D.Unit tNormal = rotation.apply(normal).normalize();
348
349        return Planes.fromPointAndNormal(tOrigin, tNormal, getPrecision());
350    }
351
352    /** Return true if this instance should be considered equivalent to the argument, using the
353     * given precision context for comparison. Instances are considered equivalent if they contain
354     * the same points, which is determined by comparing the plane {@code origins} and {@code normals}.
355     * @param other the point to compare with
356     * @param precision precision context to use for the comparison
357     * @return true if this instance should be considered equivalent to the argument
358     * @see Vector3D#eq(Vector3D, Precision.DoubleEquivalence)
359     */
360    public boolean eq(final Plane other, final Precision.DoubleEquivalence precision) {
361        return getOrigin().eq(other.getOrigin(), precision) &&
362                normal.eq(other.normal, precision);
363    }
364
365    /** {@inheritDoc} */
366    @Override
367    public int hashCode() {
368        return Objects.hash(normal, originOffset, getPrecision());
369    }
370
371    /** {@inheritDoc} */
372    @Override
373    public boolean equals(final Object obj) {
374        if (this == obj) {
375            return true;
376        } else if (obj == null || obj.getClass() != this.getClass()) {
377            return false;
378        }
379
380        final Plane other = (Plane) obj;
381
382        return Objects.equals(this.normal, other.normal) &&
383                Double.compare(this.originOffset, other.originOffset) == 0 &&
384                Objects.equals(this.getPrecision(), other.getPrecision());
385    }
386
387    /** {@inheritDoc} */
388    @Override
389    public String toString() {
390        final StringBuilder sb = new StringBuilder();
391        sb.append(getClass().getSimpleName())
392            .append("[origin= ")
393            .append(getOrigin())
394            .append(", normal= ")
395            .append(normal)
396            .append(']');
397
398        return sb.toString();
399    }
400
401    /** Get the intersection point of three planes. Returns null if no unique intersection point
402     * exists (ie, there are no intersection points or an infinite number).
403     * @param plane1 first plane1
404     * @param plane2 second plane2
405     * @param plane3 third plane2
406     * @return intersection point of the three planes or null if no unique intersection point exists
407     */
408    public static Vector3D intersection(final Plane plane1, final Plane plane2, final Plane plane3) {
409
410        // coefficients of the three planes linear equations
411        final double a1 = plane1.normal.getX();
412        final double b1 = plane1.normal.getY();
413        final double c1 = plane1.normal.getZ();
414        final double d1 = plane1.originOffset;
415
416        final double a2 = plane2.normal.getX();
417        final double b2 = plane2.normal.getY();
418        final double c2 = plane2.normal.getZ();
419        final double d2 = plane2.originOffset;
420
421        final double a3 = plane3.normal.getX();
422        final double b3 = plane3.normal.getY();
423        final double c3 = plane3.normal.getZ();
424        final double d3 = plane3.originOffset;
425
426        // direct Cramer resolution of the linear system
427        // (this is still feasible for a 3x3 system)
428        final double a23 = (b2 * c3) - (b3 * c2);
429        final double b23 = (c2 * a3) - (c3 * a2);
430        final double c23 = (a2 * b3) - (a3 * b2);
431        final double determinant = (a1 * a23) + (b1 * b23) + (c1 * c23);
432
433        // use the precision context of the first plane to determine equality
434        if (plane1.getPrecision().eqZero(determinant)) {
435            return null;
436        }
437
438        final double r = 1.0 / determinant;
439        return Vector3D.of((-a23 * d1 - (c1 * b3 - c3 * b1) * d2 - (c2 * b1 - c1 * b2) * d3) * r,
440                (-b23 * d1 - (c3 * a1 - c1 * a3) * d2 - (c1 * a2 - c2 * a1) * d3) * r,
441                (-c23 * d1 - (b1 * a3 - b3 * a1) * d2 - (b2 * a1 - b1 * a2) * d3) * r);
442    }
443}