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.math3.geometry.spherical.twod;
018
019import org.apache.commons.math3.geometry.Point;
020import org.apache.commons.math3.geometry.euclidean.threed.Rotation;
021import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;
022import org.apache.commons.math3.geometry.partitioning.Embedding;
023import org.apache.commons.math3.geometry.partitioning.Hyperplane;
024import org.apache.commons.math3.geometry.partitioning.SubHyperplane;
025import org.apache.commons.math3.geometry.partitioning.Transform;
026import org.apache.commons.math3.geometry.spherical.oned.Arc;
027import org.apache.commons.math3.geometry.spherical.oned.ArcsSet;
028import org.apache.commons.math3.geometry.spherical.oned.S1Point;
029import org.apache.commons.math3.geometry.spherical.oned.Sphere1D;
030import org.apache.commons.math3.util.FastMath;
031
032/** This class represents an oriented great circle on the 2-sphere.
033
034 * <p>An oriented circle can be defined by a center point. The circle
035 * is the the set of points that are in the normal plan the center.</p>
036
037 * <p>Since it is oriented the two spherical caps at its two sides are
038 * unambiguously identified as a left cap and a right cap. This can be
039 * used to identify the interior and the exterior in a simple way by
040 * local properties only when part of a line is used to define part of
041 * a spherical polygon boundary.</p>
042
043 * @since 3.3
044 */
045public class Circle implements Hyperplane<Sphere2D>, Embedding<Sphere2D, Sphere1D> {
046
047    /** Pole or circle center. */
048    private Vector3D pole;
049
050    /** First axis in the equator plane, origin of the phase angles. */
051    private Vector3D x;
052
053    /** Second axis in the equator plane, in quadrature with respect to x. */
054    private Vector3D y;
055
056    /** Tolerance below which close sub-arcs are merged together. */
057    private final double tolerance;
058
059    /** Build a great circle from its pole.
060     * <p>The circle is oriented in the trigonometric direction around pole.</p>
061     * @param pole circle pole
062     * @param tolerance tolerance below which close sub-arcs are merged together
063     */
064    public Circle(final Vector3D pole, final double tolerance) {
065        reset(pole);
066        this.tolerance = tolerance;
067    }
068
069    /** Build a great circle from two non-aligned points.
070     * <p>The circle is oriented from first to second point using the path smaller than \( \pi \).</p>
071     * @param first first point contained in the great circle
072     * @param second second point contained in the great circle
073     * @param tolerance tolerance below which close sub-arcs are merged together
074     */
075    public Circle(final S2Point first, final S2Point second, final double tolerance) {
076        reset(first.getVector().crossProduct(second.getVector()));
077        this.tolerance = tolerance;
078    }
079
080    /** Build a circle from its internal components.
081     * <p>The circle is oriented in the trigonometric direction around center.</p>
082     * @param pole circle pole
083     * @param x first axis in the equator plane
084     * @param y second axis in the equator plane
085     * @param tolerance tolerance below which close sub-arcs are merged together
086     */
087    private Circle(final Vector3D pole, final Vector3D x, final Vector3D y,
088                   final double tolerance) {
089        this.pole      = pole;
090        this.x         = x;
091        this.y         = y;
092        this.tolerance = tolerance;
093    }
094
095    /** Copy constructor.
096     * <p>The created instance is completely independent from the
097     * original instance, it is a deep copy.</p>
098     * @param circle circle to copy
099     */
100    public Circle(final Circle circle) {
101        this(circle.pole, circle.x, circle.y, circle.tolerance);
102    }
103
104    /** {@inheritDoc} */
105    public Circle copySelf() {
106        return new Circle(this);
107    }
108
109    /** Reset the instance as if built from a pole.
110     * <p>The circle is oriented in the trigonometric direction around pole.</p>
111     * @param newPole circle pole
112     */
113    public void reset(final Vector3D newPole) {
114        this.pole = newPole.normalize();
115        this.x    = newPole.orthogonal();
116        this.y    = Vector3D.crossProduct(newPole, x).normalize();
117    }
118
119    /** Revert the instance.
120     */
121    public void revertSelf() {
122        // x remains the same
123        y    = y.negate();
124        pole = pole.negate();
125    }
126
127    /** Get the reverse of the instance.
128     * <p>Get a circle with reversed orientation with respect to the
129     * instance. A new object is built, the instance is untouched.</p>
130     * @return a new circle, with orientation opposite to the instance orientation
131     */
132    public Circle getReverse() {
133        return new Circle(pole.negate(), x, y.negate(), tolerance);
134    }
135
136    /** {@inheritDoc} */
137    public Point<Sphere2D> project(Point<Sphere2D> point) {
138        return toSpace(toSubSpace(point));
139    }
140
141    /** {@inheritDoc} */
142    public double getTolerance() {
143        return tolerance;
144    }
145
146    /** {@inheritDoc}
147     * @see #getPhase(Vector3D)
148     */
149    public S1Point toSubSpace(final Point<Sphere2D> point) {
150        return new S1Point(getPhase(((S2Point) point).getVector()));
151    }
152
153    /** Get the phase angle of a direction.
154     * <p>
155     * The direction may not belong to the circle as the
156     * phase is computed for the meridian plane between the circle
157     * pole and the direction.
158     * </p>
159     * @param direction direction for which phase is requested
160     * @return phase angle of the direction around the circle
161     * @see #toSubSpace(Point)
162     */
163    public double getPhase(final Vector3D direction) {
164        return FastMath.PI + FastMath.atan2(-direction.dotProduct(y), -direction.dotProduct(x));
165    }
166
167    /** {@inheritDoc}
168     * @see #getPointAt(double)
169     */
170    public S2Point toSpace(final Point<Sphere1D> point) {
171        return new S2Point(getPointAt(((S1Point) point).getAlpha()));
172    }
173
174    /** Get a circle point from its phase around the circle.
175     * @param alpha phase around the circle
176     * @return circle point on the sphere
177     * @see #toSpace(Point)
178     * @see #getXAxis()
179     * @see #getYAxis()
180     */
181    public Vector3D getPointAt(final double alpha) {
182        return new Vector3D(FastMath.cos(alpha), x, FastMath.sin(alpha), y);
183    }
184
185    /** Get the X axis of the circle.
186     * <p>
187     * This method returns the same value as {@link #getPointAt(double)
188     * getPointAt(0.0)} but it does not do any computation and always
189     * return the same instance.
190     * </p>
191     * @return an arbitrary x axis on the circle
192     * @see #getPointAt(double)
193     * @see #getYAxis()
194     * @see #getPole()
195     */
196    public Vector3D getXAxis() {
197        return x;
198    }
199
200    /** Get the Y axis of the circle.
201     * <p>
202     * This method returns the same value as {@link #getPointAt(double)
203     * getPointAt(0.5 * FastMath.PI)} but it does not do any computation and always
204     * return the same instance.
205     * </p>
206     * @return an arbitrary y axis point on the circle
207     * @see #getPointAt(double)
208     * @see #getXAxis()
209     * @see #getPole()
210     */
211    public Vector3D getYAxis() {
212        return y;
213    }
214
215    /** Get the pole of the circle.
216     * <p>
217     * As the circle is a great circle, the pole does <em>not</em>
218     * belong to it.
219     * </p>
220     * @return pole of the circle
221     * @see #getXAxis()
222     * @see #getYAxis()
223     */
224    public Vector3D getPole() {
225        return pole;
226    }
227
228    /** Get the arc of the instance that lies inside the other circle.
229     * @param other other circle
230     * @return arc of the instance that lies inside the other circle
231     */
232    public Arc getInsideArc(final Circle other) {
233        final double alpha  = getPhase(other.pole);
234        final double halfPi = 0.5 * FastMath.PI;
235        return new Arc(alpha - halfPi, alpha + halfPi, tolerance);
236    }
237
238    /** {@inheritDoc} */
239    public SubCircle wholeHyperplane() {
240        return new SubCircle(this, new ArcsSet(tolerance));
241    }
242
243    /** Build a region covering the whole space.
244     * @return a region containing the instance (really a {@link
245     * SphericalPolygonsSet SphericalPolygonsSet} instance)
246     */
247    public SphericalPolygonsSet wholeSpace() {
248        return new SphericalPolygonsSet(tolerance);
249    }
250
251    /** {@inheritDoc}
252     * @see #getOffset(Vector3D)
253     */
254    public double getOffset(final Point<Sphere2D> point) {
255        return getOffset(((S2Point) point).getVector());
256    }
257
258    /** Get the offset (oriented distance) of a direction.
259     * <p>The offset is defined as the angular distance between the
260     * circle center and the direction minus the circle radius. It
261     * is therefore 0 on the circle, positive for directions outside of
262     * the cone delimited by the circle, and negative inside the cone.</p>
263     * @param direction direction to check
264     * @return offset of the direction
265     * @see #getOffset(Point)
266     */
267    public double getOffset(final Vector3D direction) {
268        return Vector3D.angle(pole, direction) - 0.5 * FastMath.PI;
269    }
270
271    /** {@inheritDoc} */
272    public boolean sameOrientationAs(final Hyperplane<Sphere2D> other) {
273        final Circle otherC = (Circle) other;
274        return Vector3D.dotProduct(pole, otherC.pole) >= 0.0;
275    }
276
277    /** Get a {@link org.apache.commons.math3.geometry.partitioning.Transform
278     * Transform} embedding a 3D rotation.
279     * @param rotation rotation to use
280     * @return a new transform that can be applied to either {@link
281     * Point Point}, {@link Circle Line} or {@link
282     * org.apache.commons.math3.geometry.partitioning.SubHyperplane
283     * SubHyperplane} instances
284     */
285    public static Transform<Sphere2D, Sphere1D> getTransform(final Rotation rotation) {
286        return new CircleTransform(rotation);
287    }
288
289    /** Class embedding a 3D rotation. */
290    private static class CircleTransform implements Transform<Sphere2D, Sphere1D> {
291
292        /** Underlying rotation. */
293        private final Rotation rotation;
294
295        /** Build a transform from a {@code Rotation}.
296         * @param rotation rotation to use
297         */
298        CircleTransform(final Rotation rotation) {
299            this.rotation = rotation;
300        }
301
302        /** {@inheritDoc} */
303        public S2Point apply(final Point<Sphere2D> point) {
304            return new S2Point(rotation.applyTo(((S2Point) point).getVector()));
305        }
306
307        /** {@inheritDoc} */
308        public Circle apply(final Hyperplane<Sphere2D> hyperplane) {
309            final Circle circle = (Circle) hyperplane;
310            return new Circle(rotation.applyTo(circle.pole),
311                              rotation.applyTo(circle.x),
312                              rotation.applyTo(circle.y),
313                              circle.tolerance);
314        }
315
316        /** {@inheritDoc} */
317        public SubHyperplane<Sphere1D> apply(final SubHyperplane<Sphere1D> sub,
318                                             final Hyperplane<Sphere2D> original,
319                                             final Hyperplane<Sphere2D> transformed) {
320            // as the circle is rotated, the limit angles are rotated too
321            return sub;
322        }
323
324    }
325
326}