1 /*
2 * Licensed to the Apache Software Foundation (ASF) under one or more
3 * contributor license agreements. See the NOTICE file distributed with
4 * this work for additional information regarding copyright ownership.
5 * The ASF licenses this file to You under the Apache License, Version 2.0
6 * (the "License"); you may not use this file except in compliance with
7 * the License. You may obtain a copy of the License at
8 *
9 * https://www.apache.org/licenses/LICENSE-2.0
10 *
11 * Unless required by applicable law or agreed to in writing, software
12 * distributed under the License is distributed on an "AS IS" BASIS,
13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14 * See the License for the specific language governing permissions and
15 * limitations under the License.
16 */
17
18 /*
19 * This is not the original file distributed by the Apache Software Foundation
20 * It has been modified by the Hipparchus project
21 */
22 package org.hipparchus.geometry.euclidean.threed;
23
24 import org.hipparchus.exception.LocalizedCoreFormats;
25 import org.hipparchus.exception.MathRuntimeException;
26 import org.hipparchus.geometry.euclidean.oned.Vector1D;
27 import org.hipparchus.geometry.euclidean.twod.Euclidean2D;
28 import org.hipparchus.geometry.euclidean.twod.PolygonsSet;
29 import org.hipparchus.geometry.euclidean.twod.SubLine;
30 import org.hipparchus.geometry.euclidean.twod.Vector2D;
31 import org.hipparchus.geometry.partitioning.Embedding;
32 import org.hipparchus.geometry.partitioning.Hyperplane;
33 import org.hipparchus.geometry.partitioning.RegionFactory;
34 import org.hipparchus.util.FastMath;
35
36 /** The class represent planes in a three dimensional space.
37 */
38 public class Plane
39 implements Hyperplane<Euclidean3D, Vector3D, Plane, SubPlane>,
40 Embedding<Euclidean3D, Vector3D, Euclidean2D, Vector2D> {
41
42 /** Offset of the origin with respect to the plane. */
43 private double originOffset;
44
45 /** Origin of the plane frame. */
46 private Vector3D origin;
47
48 /** First vector of the plane frame (in plane). */
49 private Vector3D u;
50
51 /** Second vector of the plane frame (in plane). */
52 private Vector3D v;
53
54 /** Third vector of the plane frame (plane normal). */
55 private Vector3D w;
56
57 /** Tolerance below which points are considered identical. */
58 private final double tolerance;
59
60 /** Build a plane normal to a given direction and containing the origin.
61 * @param normal normal direction to the plane
62 * @param tolerance tolerance below which points are considered identical
63 * @exception MathRuntimeException if the normal norm is too small
64 */
65 public Plane(final Vector3D normal, final double tolerance)
66 throws MathRuntimeException {
67 setNormal(normal);
68 this.tolerance = tolerance;
69 originOffset = 0;
70 setFrame();
71 }
72
73 /** Build a plane from a point and a normal.
74 * @param p point belonging to the plane
75 * @param normal normal direction to the plane
76 * @param tolerance tolerance below which points are considered identical
77 * @exception MathRuntimeException if the normal norm is too small
78 */
79 public Plane(final Vector3D p, final Vector3D normal, final double tolerance)
80 throws MathRuntimeException {
81 setNormal(normal);
82 this.tolerance = tolerance;
83 originOffset = -p.dotProduct(w);
84 setFrame();
85 }
86
87 /** Build a plane from three points.
88 * <p>The plane is oriented in the direction of
89 * {@code (p2-p1) ^ (p3-p1)}</p>
90 * @param p1 first point belonging to the plane
91 * @param p2 second point belonging to the plane
92 * @param p3 third point belonging to the plane
93 * @param tolerance tolerance below which points are considered identical
94 * @exception MathRuntimeException if the points do not constitute a plane
95 */
96 public Plane(final Vector3D p1, final Vector3D p2, final Vector3D p3, final double tolerance)
97 throws MathRuntimeException {
98 this(p1, p2.subtract(p1).crossProduct(p3.subtract(p1)), tolerance);
99 }
100
101 /** Copy constructor.
102 * <p>The instance created is completely independent of the original
103 * one. A deep copy is used, none of the underlying object are
104 * shared.</p>
105 * @param plane plane to copy
106 */
107 public Plane(final Plane plane) {
108 originOffset = plane.originOffset;
109 origin = plane.origin;
110 u = plane.u;
111 v = plane.v;
112 w = plane.w;
113 tolerance = plane.tolerance;
114 }
115
116 /** Copy the instance.
117 * <p>The instance created is completely independant of the original
118 * one. A deep copy is used, none of the underlying objects are
119 * shared (except for immutable objects).</p>
120 * @return a new hyperplane, copy of the instance
121 */
122 @Override
123 public Plane copySelf() {
124 return new Plane(this);
125 }
126
127 /** Reset the instance as if built from a point and a normal.
128 * @param p point belonging to the plane
129 * @param normal normal direction to the plane
130 * @exception MathRuntimeException if the normal norm is too small
131 */
132 public void reset(final Vector3D p, final Vector3D normal) throws MathRuntimeException {
133 setNormal(normal);
134 originOffset = -p.dotProduct(w);
135 setFrame();
136 }
137
138 /** Reset the instance from another one.
139 * <p>The updated instance is completely independant of the original
140 * one. A deep reset is used none of the underlying object is
141 * shared.</p>
142 * @param original plane to reset from
143 */
144 public void reset(final Plane original) {
145 originOffset = original.originOffset;
146 origin = original.origin;
147 u = original.u;
148 v = original.v;
149 w = original.w;
150 }
151
152 /** Set the normal vactor.
153 * @param normal normal direction to the plane (will be copied)
154 * @exception MathRuntimeException if the normal norm is too small
155 */
156 private void setNormal(final Vector3D normal) throws MathRuntimeException {
157 final double norm = normal.getNorm();
158 if (norm < 1.0e-10) {
159 throw new MathRuntimeException(LocalizedCoreFormats.ZERO_NORM);
160 }
161 w = new Vector3D(1.0 / norm, normal);
162 }
163
164 /** Reset the plane frame.
165 */
166 private void setFrame() {
167 origin = new Vector3D(-originOffset, w);
168 u = w.orthogonal();
169 v = Vector3D.crossProduct(w, u);
170 }
171
172 /** Get the origin point of the plane frame.
173 * <p>The point returned is the orthogonal projection of the
174 * 3D-space origin in the plane.</p>
175 * @return the origin point of the plane frame (point closest to the
176 * 3D-space origin)
177 */
178 public Vector3D getOrigin() {
179 return origin;
180 }
181
182 /** Get the normalized normal vector.
183 * <p>The frame defined by ({@link #getU() getU()}, {@link #getV() getV()},
184 * {@code getNormal()}) is a right-handed orthonormalized
185 * frame).</p>
186 * @return normalized normal vector
187 * @see #getU
188 * @see #getV
189 */
190 public Vector3D getNormal() {
191 return w;
192 }
193
194 /** Get the plane first canonical vector.
195 * <p>The frame defined by ({@code getU()}, {@link #getV() getV()},
196 * {@link #getNormal() getNormal()}) is a right-handed orthonormalized
197 * frame).</p>
198 * @return normalized first canonical vector
199 * @see #getV
200 * @see #getNormal
201 */
202 public Vector3D getU() {
203 return u;
204 }
205
206 /** Get the plane second canonical vector.
207 * <p>The frame defined by ({@link #getU() getU()}, {@code getV()},
208 * {@link #getNormal() getNormal()}) is a right-handed orthonormalized
209 * frame).</p>
210 * @return normalized second canonical vector
211 * @see #getU
212 * @see #getNormal
213 */
214 public Vector3D getV() {
215 return v;
216 }
217
218 /** {@inheritDoc}
219 */
220 @Override
221 public Vector3D project(Vector3D point) {
222 return toSpace(toSubSpace(point));
223 }
224
225 /** {@inheritDoc}
226 */
227 @Override
228 public double getTolerance() {
229 return tolerance;
230 }
231
232 /** Revert the plane.
233 * <p>Replace the instance by a similar plane with opposite orientation.</p>
234 * <p>The new plane frame is chosen in such a way that a 3D point that had
235 * {@code (x, y)} in-plane coordinates and {@code z} offset with
236 * respect to the plane and is unaffected by the change will have
237 * {@code (y, x)} in-plane coordinates and {@code -z} offset with
238 * respect to the new plane. This means that the {@code u} and {@code v}
239 * vectors returned by the {@link #getU} and {@link #getV} methods are exchanged,
240 * and the {@code w} vector returned by the {@link #getNormal} method is
241 * reversed.</p>
242 */
243 public void revertSelf() {
244 final Vector3D tmp = u;
245 u = v;
246 v = tmp;
247 w = w.negate();
248 originOffset = -originOffset;
249 }
250
251 /** Transform a 3D space point into an in-plane point.
252 * @param point point of the space (must be a {@link Vector3D
253 * Vector3D} instance)
254 * @return in-plane point (really a {@link
255 * org.hipparchus.geometry.euclidean.twod.Vector2D Vector2D} instance)
256 * @see #toSpace
257 */
258 @Override
259 public Vector2D toSubSpace(final Vector3D point) {
260 return new Vector2D(point.dotProduct(u), point.dotProduct(v));
261 }
262
263 /** Transform an in-plane point into a 3D space point.
264 * @param point in-plane point (must be a {@link
265 * org.hipparchus.geometry.euclidean.twod.Vector2D Vector2D} instance)
266 * @return 3D space point (really a {@link Vector3D Vector3D} instance)
267 * @see #toSubSpace
268 */
269 @Override
270 public Vector3D toSpace(final Vector2D point) {
271 return new Vector3D(point.getX(), u, point.getY(), v, -originOffset, w);
272 }
273
274 /** Get one point from the 3D-space.
275 * @param inPlane desired in-plane coordinates for the point in the
276 * plane
277 * @param offset desired offset for the point
278 * @return one point in the 3D-space, with given coordinates and offset
279 * relative to the plane
280 */
281 public Vector3D getPointAt(final Vector2D inPlane, final double offset) {
282 return new Vector3D(inPlane.getX(), u, inPlane.getY(), v, offset - originOffset, w);
283 }
284
285 /** Check if the instance is similar to another plane.
286 * <p>Planes are considered similar if they contain the same
287 * points. This does not mean they are equal since they can have
288 * opposite normals.</p>
289 * @param plane plane to which the instance is compared
290 * @return true if the planes are similar
291 */
292 public boolean isSimilarTo(final Plane plane) {
293 final double angle = Vector3D.angle(w, plane.w);
294 return ((angle < 1.0e-10) && (FastMath.abs(originOffset - plane.originOffset) < tolerance)) ||
295 ((angle > (FastMath.PI - 1.0e-10)) && (FastMath.abs(originOffset + plane.originOffset) < tolerance));
296 }
297
298 /** Rotate the plane around the specified point.
299 * <p>The instance is not modified, a new instance is created.</p>
300 * @param center rotation center
301 * @param rotation vectorial rotation operator
302 * @return a new plane
303 */
304 public Plane rotate(final Vector3D center, final Rotation rotation) {
305
306 final Vector3D delta = origin.subtract(center);
307 final Plane plane = new Plane(center.add(rotation.applyTo(delta)),
308 rotation.applyTo(w), tolerance);
309
310 // make sure the frame is transformed as desired
311 plane.u = rotation.applyTo(u);
312 plane.v = rotation.applyTo(v);
313
314 return plane;
315
316 }
317
318 /** Translate the plane by the specified amount.
319 * <p>The instance is not modified, a new instance is created.</p>
320 * @param translation translation to apply
321 * @return a new plane
322 */
323 public Plane translate(final Vector3D translation) {
324
325 final Plane plane = new Plane(origin.add(translation), w, tolerance);
326
327 // make sure the frame is transformed as desired
328 plane.u = u;
329 plane.v = v;
330
331 return plane;
332
333 }
334
335 /** Get the intersection of a line with the instance.
336 * @param line line intersecting the instance
337 * @return intersection point between between the line and the
338 * instance (null if the line is parallel to the instance)
339 */
340 public Vector3D intersection(final Line line) {
341 final Vector3D direction = line.getDirection();
342 final double dot = w.dotProduct(direction);
343 if (FastMath.abs(dot) < 1.0e-10) {
344 return null;
345 }
346 final Vector3D point = line.toSpace(Vector1D.ZERO);
347 final double k = -(originOffset + w.dotProduct(point)) / dot;
348 return new Vector3D(1.0, point, k, direction);
349 }
350
351 /** Build the line shared by the instance and another plane.
352 * @param other other plane
353 * @return line at the intersection of the instance and the
354 * other plane (really a {@link Line Line} instance)
355 */
356 public Line intersection(final Plane other) {
357 final Vector3D direction = Vector3D.crossProduct(w, other.w);
358 if (direction.getNorm() < tolerance) {
359 return null;
360 }
361 final Vector3D point = intersection(this, other, new Plane(direction, tolerance));
362 return new Line(point, point.add(direction), tolerance);
363 }
364
365 /** Get the intersection point of three planes.
366 * @param plane1 first plane1
367 * @param plane2 second plane2
368 * @param plane3 third plane2
369 * @return intersection point of three planes, null if some planes are parallel
370 */
371 public static Vector3D intersection(final Plane plane1, final Plane plane2, final Plane plane3) {
372
373 // coefficients of the three planes linear equations
374 final double a1 = plane1.w.getX();
375 final double b1 = plane1.w.getY();
376 final double c1 = plane1.w.getZ();
377 final double d1 = plane1.originOffset;
378
379 final double a2 = plane2.w.getX();
380 final double b2 = plane2.w.getY();
381 final double c2 = plane2.w.getZ();
382 final double d2 = plane2.originOffset;
383
384 final double a3 = plane3.w.getX();
385 final double b3 = plane3.w.getY();
386 final double c3 = plane3.w.getZ();
387 final double d3 = plane3.originOffset;
388
389 // direct Cramer resolution of the linear system
390 // (this is still feasible for a 3x3 system)
391 final double a23 = b2 * c3 - b3 * c2;
392 final double b23 = c2 * a3 - c3 * a2;
393 final double c23 = a2 * b3 - a3 * b2;
394 final double determinant = a1 * a23 + b1 * b23 + c1 * c23;
395 if (FastMath.abs(determinant) < 1.0e-10) {
396 return null;
397 }
398
399 final double r = 1.0 / determinant;
400 return new Vector3D(
401 (-a23 * d1 - (c1 * b3 - c3 * b1) * d2 - (c2 * b1 - c1 * b2) * d3) * r,
402 (-b23 * d1 - (c3 * a1 - c1 * a3) * d2 - (c1 * a2 - c2 * a1) * d3) * r,
403 (-c23 * d1 - (b1 * a3 - b3 * a1) * d2 - (b2 * a1 - b1 * a2) * d3) * r);
404
405 }
406
407 /** Build a region covering the whole hyperplane.
408 * @return a region covering the whole hyperplane
409 */
410 @Override
411 public SubPlane wholeHyperplane() {
412 return new SubPlane(this, new PolygonsSet(tolerance));
413 }
414
415 /** {@inheritDoc} */
416 @Override
417 public SubPlane emptyHyperplane() {
418 final RegionFactory<Euclidean2D, Vector2D, org.hipparchus.geometry.euclidean.twod.Line, SubLine> factory = new RegionFactory<>();
419 return new SubPlane(this, factory.getComplement(new PolygonsSet(tolerance)));
420 }
421
422 /** Build a region covering the whole space.
423 * @return a region containing the instance (really a {@link
424 * PolyhedronsSet PolyhedronsSet} instance)
425 */
426 @Override
427 public PolyhedronsSet wholeSpace() {
428 return new PolyhedronsSet(tolerance);
429 }
430
431 /** Check if the instance contains a point.
432 * @param p point to check
433 * @return true if p belongs to the plane
434 */
435 public boolean contains(final Vector3D p) {
436 return FastMath.abs(getOffset(p)) < tolerance;
437 }
438
439 /** Get the offset (oriented distance) of a parallel plane.
440 * <p>This method should be called only for parallel planes otherwise
441 * the result is not meaningful.</p>
442 * <p>The offset is 0 if both planes are the same, it is
443 * positive if the plane is on the plus side of the instance and
444 * negative if it is on the minus side, according to its natural
445 * orientation.</p>
446 * @param plane plane to check
447 * @return offset of the plane
448 */
449 public double getOffset(final Plane plane) {
450 return originOffset + (sameOrientationAs(plane) ? -plane.originOffset : plane.originOffset);
451 }
452
453 /** Get the offset (oriented distance) of a point.
454 * <p>The offset is 0 if the point is on the underlying hyperplane,
455 * it is positive if the point is on one particular side of the
456 * hyperplane, and it is negative if the point is on the other side,
457 * according to the hyperplane natural orientation.</p>
458 * @param point point to check
459 * @return offset of the point
460 */
461 @Override
462 public double getOffset(final Vector3D point) {
463 return point.dotProduct(w) + originOffset;
464 }
465
466 /** {@inheritDoc} */
467 @Override
468 public Vector3D moveToOffset(final Vector3D point, final double offset) {
469 final double delta = offset - getOffset(point);
470 return new Vector3D(point.getX() + delta * w.getX(),
471 point.getY() + delta * w.getY(),
472 point.getZ() + delta * w.getZ());
473 }
474
475 /** {@inheritDoc} */
476 @Override
477 public Vector3D arbitraryPoint() {
478 return origin;
479 }
480
481 /** Check if the instance has the same orientation as another hyperplane.
482 * @param other other hyperplane to check against the instance
483 * @return true if the instance and the other hyperplane have
484 * the same orientation
485 */
486 @Override
487 public boolean sameOrientationAs(final Plane other) {
488 return other.w.dotProduct(w) > 0.0;
489 }
490
491 }