1   /* Copyright 2002-2026 CS GROUP
2    * Licensed to CS GROUP (CS) under one or more
3    * contributor license agreements.  See the NOTICE file distributed with
4    * this work for additional information regarding copyright ownership.
5    * CS 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    *   http://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  package org.orekit.geometry.fov;
18  
19  import java.util.ArrayList;
20  import java.util.List;
21  
22  import org.hipparchus.geometry.enclosing.EnclosingBall;
23  import org.hipparchus.geometry.euclidean.threed.Line;
24  import org.hipparchus.geometry.euclidean.threed.Rotation;
25  import org.hipparchus.geometry.euclidean.threed.RotationConvention;
26  import org.hipparchus.geometry.euclidean.threed.Vector3D;
27  import org.hipparchus.geometry.partitioning.Region;
28  import org.hipparchus.geometry.spherical.twod.Edge;
29  import org.hipparchus.geometry.spherical.twod.S2Point;
30  import org.hipparchus.geometry.spherical.twod.Sphere2D;
31  import org.hipparchus.geometry.spherical.twod.SphericalPolygonsSet;
32  import org.hipparchus.geometry.spherical.twod.Vertex;
33  import org.hipparchus.util.FastMath;
34  import org.hipparchus.util.MathUtils;
35  import org.hipparchus.util.SinCos;
36  import org.orekit.bodies.GeodeticPoint;
37  import org.orekit.bodies.OneAxisEllipsoid;
38  import org.orekit.errors.OrekitException;
39  import org.orekit.errors.OrekitMessages;
40  import org.orekit.frames.Frame;
41  import org.orekit.frames.Transform;
42  import org.orekit.propagation.events.VisibilityTrigger;
43  
44  /** Class representing a spacecraft sensor Field Of View with polygonal shape.
45   * <p>Fields Of View are zones defined on the unit sphere centered on the
46   * spacecraft. They can have any shape, they can be split in several
47   * non-connected patches and can have holes.</p>
48   * @author Luc Maisonobe
49   * @since 10.1
50   */
51  public class PolygonalFieldOfView extends AbstractFieldOfView {
52  
53      /** Spherical zone. */
54      private final SphericalPolygonsSet zone;
55  
56      /** Spherical cap surrounding the zone. */
57      private final EnclosingBall<Sphere2D, S2Point> cap;
58  
59      /** Build a new instance.
60       * @param zone interior of the Field Of View, in spacecraft frame
61       * @param margin angular margin to apply to the zone (if positive,
62       * points outside of the raw FoV but close enough to the boundary are
63       * considered visible; if negative, points inside of the raw FoV
64       * but close enough to the boundary are considered not visible)
65       */
66      public PolygonalFieldOfView(final SphericalPolygonsSet zone, final double margin) {
67          super(margin);
68          this.zone = zone;
69          this.cap  = zone.getEnclosingCap();
70      }
71  
72      /** Build Field Of View with a regular polygon shape.
73       * @param center center of the polygon (the center is in the inside part)
74       * @param coneType type of defining cone
75       * @param meridian point defining the reference meridian for one contact
76       * point between defining cone and polygon (i.e. either a polygon edge
77       * middle point or a polygon vertex)
78       * @param radius defining cone angular radius
79       * @param n number of sides of the polygon
80       * @param margin angular margin to apply to the zone (if positive,
81       * points outside of the raw FoV but close enough to the boundary are
82       * considered visible; if negative, points inside of the raw FoV
83       * but close enough to the boundary are considered not visible)
84       * @since 10.1
85       */
86      public PolygonalFieldOfView(final Vector3D center, final DefiningConeType coneType,
87                                  final Vector3D meridian, final double radius,
88                                  final int n, final double margin) {
89  
90          super(margin);
91  
92          final double   verticesRadius = coneType.verticesRadius(radius, n);
93          final Vector3D vertex         = coneType.createVertex(center, meridian, verticesRadius, n);
94          this.zone                     = new SphericalPolygonsSet(center, vertex, verticesRadius,
95                                                                   n, 1.0e-12 * verticesRadius);
96  
97          final Rotation r = new Rotation(center, MathUtils.TWO_PI / n, RotationConvention.VECTOR_OPERATOR);
98          final S2Point[] support = new S2Point[n];
99          support[0] = new S2Point(vertex);
100         for (int i = 1; i < n; ++i) {
101             support[i] = new S2Point(r.applyTo(support[i - 1].getVector()));
102         }
103         this.cap = new EnclosingBall<>(new S2Point(center), Vector3D.angle(center, vertex), support);
104 
105     }
106 
107     /** Get the interior zone.
108      * @return the interior zone
109      */
110     public SphericalPolygonsSet getZone() {
111         return zone;
112     }
113 
114     /** {@inheritDoc} */
115     @Override
116     public double offsetFromBoundary(final Vector3D lineOfSight, final double angularRadius,
117                                      final VisibilityTrigger trigger) {
118 
119         final S2Point los             = new S2Point(lineOfSight);
120         final double  margin          = getMargin();
121         final double  correctedRadius = trigger.radiusCorrection(angularRadius);
122         final double  deadBand        = margin + angularRadius;
123 
124         // for faster computation, we start using only the surrounding cap, to filter out
125         // far away points (which correspond to most of the points if the Field Of View is small)
126         final double crudeDistance = cap.getCenter().distance(los) - cap.getRadius();
127         if (crudeDistance > deadBand + 0.01) {
128             // we know we are strictly outside of the zone,
129             // use the crude distance to compute the (positive) return value
130             return crudeDistance + correctedRadius - margin;
131         }
132 
133         // we are close, we need to compute carefully the exact offset;
134         // we project the point to the closest zone boundary
135         return zone.projectToBoundary(los).getOffset() + correctedRadius - margin;
136 
137     }
138 
139     /** {@inheritDoc} */
140     @Override
141     public Vector3D projectToBoundary(final Vector3D lineOfSight) {
142         return (zone.projectToBoundary(new S2Point(lineOfSight)).getProjected()).getVector();
143     }
144 
145     /** {@inheritDoc} */
146     @Override
147     public List<List<GeodeticPoint>> getFootprint(final Transform fovToBody,
148                                                   final OneAxisEllipsoid body,
149                                                   final double angularStep) {
150 
151         return getFootprint(fovToBody, body, angularStep, DEFAULT_EXTERNAL_FOOTPRINT, DEFAULT_MAX_DIST);
152 
153     }
154 
155     /** {@inheritDoc} */
156     @Override
157     public List<List<GeodeticPoint>> getFootprint(final Transform fovToBody,
158                                                   final OneAxisEllipsoid body,
159                                                   final double angularStep,
160                                                   final boolean extendedFootprint,
161                                                   final double maxDistance) {
162 
163         final Frame     bodyFrame = body.getBodyFrame();
164         final Vector3D  position  = fovToBody.transformPosition(Vector3D.ZERO);
165         final double    r         = position.getNorm();
166         if (body.isInside(position)) {
167             throw new OrekitException(OrekitMessages.POINT_INSIDE_ELLIPSOID);
168         }
169 
170         final List<List<GeodeticPoint>> footprint = new ArrayList<>();
171 
172         final List<Vertex> boundary = zone.getBoundaryLoops();
173         for (final Vertex loopStart : boundary) {
174             int count = 0;
175             final List<GeodeticPoint> loop  = new ArrayList<>();
176             boolean intersectionsFound      = false;
177             for (Edge edge = loopStart.getOutgoing();
178                  count == 0 || edge.getStart() != loopStart;
179                  edge = edge.getEnd().getOutgoing()) {
180                 ++count;
181                 final int    n     = (int) FastMath.ceil(edge.getLength() / angularStep);
182                 final double delta =  edge.getLength() / n;
183                 for (int i = 0; i < n; ++i) {
184                     final Vector3D awaySC      = new Vector3D(r, edge.getPointAt(i * delta));
185                     final Vector3D awayBody    = fovToBody.transformPosition(awaySC);
186                     final Line     lineOfSight = new Line(position, awayBody, 1.0e-3);
187                     GeodeticPoint  gp          = body.getIntersectionPoint(lineOfSight, position,
188                                                                            bodyFrame, null);
189                     if (gp != null &&
190                         Vector3D.dotProduct(awayBody.subtract(position),
191                                             body.transform(gp).subtract(position)) < 0) {
192                         // the intersection is in fact on the half-line pointing
193                         // towards the back side, it is a spurious intersection
194                         gp = null;
195                     }
196 
197                     if (gp != null) {
198                         // the line of sight does intersect the body
199                         intersectionsFound = true;
200                     } else {
201                         // the line of sight does not intersect body
202                         if (!extendedFootprint) {
203                             // we use a point on the limb
204                             gp = body.transform(body.pointOnLimb(position, awayBody), bodyFrame, null);
205                         } else {
206                             // we project in the proper direction (to a point in space)
207                             // such that the line's length is the requested value (maxDistance)
208                             final double temp_norm = awaySC.getNorm();
209                             Vector3D temp_vec = awaySC.scalarMultiply(1 / temp_norm).scalarMultiply(maxDistance);
210                             temp_vec = fovToBody.transformPosition(temp_vec);
211                             gp = body.transform(temp_vec, bodyFrame, null);
212                         }
213                     }
214 
215                     // add the point in front of the list
216                     // (to ensure the loop will be in trigonometric orientation)
217                     loop.addFirst(gp);
218 
219                 }
220             }
221 
222             if (!extendedFootprint) {
223                 if (intersectionsFound) {
224                     // at least some of the points did intersect the body,
225                     // this loop contributes to the footprint
226                     footprint.add(loop);
227                 }
228             } else {
229                 // we don't care about intersection details
230                 footprint.add(loop);
231             }
232 
233         }
234 
235         if (!extendedFootprint) { // only need to check if we care about non-intersection
236             if (footprint.isEmpty()) {
237                 // none of the Field Of View loops cross the body
238                 // either the body is outside of Field Of View, or it is fully contained
239                 // we check the center
240                 final Vector3D bodyCenter = fovToBody.getStaticInverse().transformPosition(Vector3D.ZERO);
241                 if (zone.checkPoint(new S2Point(bodyCenter)) != Region.Location.OUTSIDE) {
242                     // the body is fully contained in the Field Of View
243                     // we use the full limb as the footprint
244                     final Vector3D x        = bodyCenter.orthogonal();
245                     final Vector3D y        = Vector3D.crossProduct(bodyCenter, x).normalize();
246                     final double   sinEta   = body.getEquatorialRadius() / r;
247                     final double   sinEta2  = sinEta * sinEta;
248                     final double   cosAlpha = (FastMath.cos(angularStep) + sinEta2 - 1) / sinEta2;
249                     final int      n        = (int) FastMath.ceil(MathUtils.TWO_PI / FastMath.acos(cosAlpha));
250                     final double   delta    = MathUtils.TWO_PI / n;
251                     final List<GeodeticPoint> loop = new ArrayList<>(n);
252                     for (int i = 0; i < n; ++i) {
253                         final SinCos sc = FastMath.sinCos(i * delta);
254                         final Vector3D outside = new Vector3D(r * sc.cos(), x,
255                                                             r * sc.sin(), y);
256                         loop.add(body.transform(body.pointOnLimb(position, outside), bodyFrame, null));
257                     }
258                     footprint.add(loop);
259                 }
260             }
261         }
262 
263         return footprint;
264 
265     }
266 
267     /** Enumerate for cone/polygon relative position.
268      * @since 10.1
269      */
270     public enum DefiningConeType {
271 
272         /** Constant for cones inside polygons and touching it at polygon edges middle points. */
273         INSIDE_CONE_TOUCHING_POLYGON_AT_EDGES_MIDDLE() {
274 
275             /** {@inheritDoc}*/
276             @Override
277             protected double verticesRadius(final double radius, final int n) {
278                 // convert the inside (edges middle points) radius to outside (vertices) radius
279                 return FastMath.atan(FastMath.tan(radius) / FastMath.cos(FastMath.PI / n));
280             }
281 
282             /** {@inheritDoc}*/
283             @Override
284             protected Vector3D createVertex(final Vector3D center, final Vector3D meridian,
285                                             final double verticesRadius, final int n) {
286                 // convert the edge middle meridian to a vertex
287                 final SinCos scA = FastMath.sinCos(FastMath.PI / n);
288                 final SinCos scR = FastMath.sinCos(verticesRadius);
289                 final Vector3D z = center.normalize();
290                 final Vector3D y = Vector3D.crossProduct(center, meridian).normalize();
291                 final Vector3D x = Vector3D.crossProduct(y, z);
292                 return new Vector3D(scR.sin() * scA.cos(), x, scR.sin() * scA.sin(), y, scR.cos(), z);
293             }
294 
295         },
296 
297         /** Constant for cones outside polygons and touching it at polygon vertices. */
298         OUTSIDE_CONE_TOUCHING_POLYGON_AT_VERTICES() {
299 
300             /** {@inheritDoc}*/
301             @Override
302             protected double verticesRadius(final double radius, final int n) {
303                 return radius;
304             }
305 
306             /** {@inheritDoc}*/
307             @Override
308             protected Vector3D createVertex(final Vector3D center, final Vector3D meridian,
309                                             final double verticesRadius, final int n) {
310                 // convert the vertex meridian to a vertex
311                 final SinCos scR = FastMath.sinCos(verticesRadius);
312                 final Vector3D z = center.normalize();
313                 final Vector3D y = Vector3D.crossProduct(center, meridian).normalize();
314                 final Vector3D x = Vector3D.crossProduct(y, z);
315                 return new Vector3D(scR.sin(), x, scR.cos(), z);
316             }
317 
318         };
319 
320         /** Compute radius of cone going through vertices.
321          * @param radius defining cone angular radius
322          * @param n number of sides of the polygon
323          * @return radius of cone going through vertices
324          */
325         protected abstract double verticesRadius(double radius, int n);
326 
327         /** Create a vertex.
328          * @param center center of the polygon (the center is in the inside part)
329          * @param meridian point defining the reference meridian for one contact
330          * point between defining cone and polygon (i.e. either a polygon edge
331          * middle point or a polygon vertex)
332          * @param verticesRadius defining radius of cone passing through vertices
333          * @param n number of sides of the polygon
334          * @return created vertex
335          */
336         protected abstract Vector3D createVertex(Vector3D center, Vector3D meridian,
337                                                  double verticesRadius, int n);
338 
339     }
340 
341 }