diff --git a/src/main/java/org/orekit/rugged/adjustment/AdjustmentContext.java b/src/main/java/org/orekit/rugged/adjustment/AdjustmentContext.java index bb55d5ed3ffac022a53aee90c623f8fa926f99f9..10deaa25f966206a2b363982d448e1419057d3b4 100644 --- a/src/main/java/org/orekit/rugged/adjustment/AdjustmentContext.java +++ b/src/main/java/org/orekit/rugged/adjustment/AdjustmentContext.java @@ -58,8 +58,7 @@ public class AdjustmentContext { public AdjustmentContext(final List<Rugged> viewingModel, final Observables measures) { this.viewingModel = new HashMap<String, Rugged>(); for (final Rugged r : viewingModel) { - /** TODO this.viewingModel.put(r.getName(), r); */ - this.viewingModel.put("Rugged", r); + this.viewingModel.put(r.getName(), r); } this.measures = measures; this.optimizerID = OptimizerId.GAUSS_NEWTON_QR; diff --git a/src/main/java/org/orekit/rugged/api/Rugged.java b/src/main/java/org/orekit/rugged/api/Rugged.java index 698d4e5d317cbaa0b2449c4a30c1c16a6a81822a..df4fa50609b609e693b46ea937eaa96f44b075aa 100644 --- a/src/main/java/org/orekit/rugged/api/Rugged.java +++ b/src/main/java/org/orekit/rugged/api/Rugged.java @@ -13,45 +13,24 @@ */ package org.orekit.rugged.api; -import java.util.ArrayList; import java.util.Collection; import java.util.HashMap; -import java.util.Iterator; -import java.util.List; import java.util.Map; import org.hipparchus.analysis.differentiation.DerivativeStructure; import org.hipparchus.geometry.euclidean.threed.FieldVector3D; import org.hipparchus.geometry.euclidean.threed.Vector3D; -import org.hipparchus.linear.Array2DRowRealMatrix; -import org.hipparchus.linear.ArrayRealVector; -import org.hipparchus.linear.DiagonalMatrix; -import org.hipparchus.linear.RealMatrix; -import org.hipparchus.linear.RealVector; -import org.hipparchus.optim.ConvergenceChecker; -import org.hipparchus.optim.nonlinear.vector.leastsquares.GaussNewtonOptimizer; -import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresBuilder; -import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer; -import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum; -import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresProblem; -import org.hipparchus.optim.nonlinear.vector.leastsquares.MultivariateJacobianFunction; -import org.hipparchus.optim.nonlinear.vector.leastsquares.ParameterValidator; import org.hipparchus.util.FastMath; -import org.hipparchus.util.Pair; import org.orekit.bodies.GeodeticPoint; -import org.orekit.errors.OrekitException; -import org.orekit.errors.OrekitExceptionWrapper; import org.orekit.frames.Transform; import org.orekit.rugged.errors.DumpManager; import org.orekit.rugged.errors.RuggedException; -import org.orekit.rugged.errors.RuggedExceptionWrapper; import org.orekit.rugged.errors.RuggedMessages; import org.orekit.rugged.intersection.IntersectionAlgorithm; import org.orekit.rugged.linesensor.LineSensor; import org.orekit.rugged.linesensor.SensorMeanPlaneCrossing; import org.orekit.rugged.linesensor.SensorPixel; import org.orekit.rugged.linesensor.SensorPixelCrossing; -import org.orekit.rugged.refining.measures.SensorToSensorMapping; import org.orekit.rugged.utils.DSGenerator; import org.orekit.rugged.utils.ExtendedEllipsoid; import org.orekit.rugged.utils.NormalizedGeodeticPoint; @@ -59,8 +38,6 @@ import org.orekit.rugged.utils.SpacecraftToObservedBody; import org.orekit.time.AbsoluteDate; import org.orekit.utils.Constants; import org.orekit.utils.PVCoordinates; -import org.orekit.utils.ParameterDriver; -import org.orekit.utils.ParameterDriversList; /** * Main class of Rugged library API. @@ -105,6 +82,9 @@ public class Rugged { /** Flag for aberration of light correction. */ private boolean aberrationOfLightCorrection; + /** Rugged name. */ + private final String name; + /** * Build a configured instance. * <p> @@ -127,12 +107,14 @@ public class Rugged { * for more accurate location * @param scToBody transforms interpolator * @param sensors sensors + * @param name RuggedName */ Rugged(final IntersectionAlgorithm algorithm, final ExtendedEllipsoid ellipsoid, final boolean lightTimeCorrection, final boolean aberrationOfLightCorrection, final SpacecraftToObservedBody scToBody, - final Collection<LineSensor> sensors) { + final Collection<LineSensor> sensors, + final String name) { // space reference this.ellipsoid = ellipsoid; @@ -143,6 +125,9 @@ public class Rugged { // intersection algorithm this.algorithm = algorithm; + // rugged name + this.name = name; + this.sensors = new HashMap<String, LineSensor>(); for (final LineSensor s : sensors) { this.sensors.put(s.getName(), s); @@ -154,6 +139,17 @@ public class Rugged { } + /** + * Get the ruggd name. + * + * @return rugged name + */ + public String getName() { + return name; + } + + + /** * Get the DEM intersection algorithm. * @@ -988,227 +984,6 @@ public class Rugged { return new DerivativeStructure[] {d, dEarth}; } - /** - * Estimate the free parameters from two viewing models (A and B) - * - * @param references reference mappings between two sensors pixels from two - * models and the corresponding computed distance between LOS that - * should ultimately be reached by adjusting selected viewing models - * parameters - * @param maxEvaluations maximum number of evaluations - * @param parametersConvergenceThreshold convergence threshold on normalized - * parameters (dimensionless, related to parameters scales) - * @param ruggedA rugged instance from viewing model A - * @return optimum of the least squares problem - * @exception RuggedException if several parameters with the same name - * exist, if no parameters have been selected for estimation, or - * if parameters cannot be estimated (too few measurements, - * ill-conditioned problem ...) - */ - public Optimum estimateFreeParams2Models( - final Collection<SensorToSensorMapping> references, - final int maxEvaluations, - final double parametersConvergenceThreshold, - Rugged ruggedA) - throws RuggedException { - try { - - final List<LineSensor> selectedSensors = new ArrayList<>(); - for (final SensorToSensorMapping reference : references) { - selectedSensors - .add(ruggedA.getLineSensor(reference.getSensorNameA())); // from - // ruggedA - // instance - selectedSensors - .add(this.getLineSensor(reference.getSensorNameB())); // from - // current - // ruggedB - // instance - - } - final DSGenerator generator = createGenerator(selectedSensors); - final ParameterDriversList selected = generator.getSelected(); - if (selected.getNbParams() == 0) { - throw new RuggedException(RuggedMessages.NO_PARAMETERS_SELECTED); - } - - int iStart = 0; - // get start point (as a normalized value) - final double[] start = new double[selected.getNbParams()]; - for (final ParameterDriver driver : selected.getDrivers()) { - start[iStart++] = driver.getNormalizedValue(); - } - - // set up target : distance between two LOS from both viewing models - // (A and B) - int n = 0; - for (final SensorToSensorMapping reference : references) { - n += reference.getMapping().size(); - } - - if (n == 0) { - throw new RuggedException(RuggedMessages.NO_REFERENCE_MAPPINGS); - } - - n = 2 * n; - - final double[] target = new double[n]; - final double[] weight = new double[n]; - int k = 0; - for (final SensorToSensorMapping reference : references) { - // Get Earth constraint weight - double earthConstraintWeight = reference.getEarthConstraintWeight(); - - int i=0; - for(Iterator<Map.Entry<SensorPixel, SensorPixel>> gtIt = reference.getMapping().iterator(); - gtIt.hasNext();i++){ - - if(i==reference.getMapping().size()) break; - - // Get LOS distance - Double losDistance = reference.getLosDistance(i); - weight[k] = 1.0 - earthConstraintWeight; - target[k++] = losDistance.doubleValue(); - - // Get Earth distance (constraint) - Double earthDistance = reference.getEarthDistance(i); - weight[k] = earthConstraintWeight; - target[k++] = earthDistance.doubleValue(); - } - } - - // prevent parameters to exceed their prescribed bounds - final ParameterValidator validator = params -> { - try { - int i = 0; - for (final ParameterDriver driver : selected.getDrivers()) { - // let the parameter handle min/max clipping - driver.setNormalizedValue(params.getEntry(i)); - params.setEntry(i++, driver.getNormalizedValue()); - } - return params; - } catch (OrekitException oe) { - throw new OrekitExceptionWrapper(oe); - } - }; - - // convergence checker - final ConvergenceChecker<LeastSquaresProblem.Evaluation> checker = (iteration, - previous, - current) -> current - .getPoint() - .getLInfDistance(previous - .getPoint()) <= parametersConvergenceThreshold; - - // model function - final MultivariateJacobianFunction model = point -> { - try { - - // set the current parameters values - int i = 0; - for (final ParameterDriver driver : selected.getDrivers()) { - driver.setNormalizedValue(point.getEntry(i++)); - } - - // compute distance and its partial derivatives - final RealVector value = new ArrayRealVector(target.length); - final RealMatrix jacobian = new Array2DRowRealMatrix(target.length, - selected - .getNbParams()); - int l = 0; - for (final SensorToSensorMapping reference : references) { - for (final Map.Entry<SensorPixel, SensorPixel> mapping : reference - .getMapping()) { - final SensorPixel spA = mapping.getKey(); - final SensorPixel spB = mapping.getValue(); - LineSensor lineSensorB = this - .getLineSensor(reference.getSensorNameB()); - LineSensor lineSensorA = ruggedA - .getLineSensor(reference.getSensorNameA()); - final AbsoluteDate dateA = lineSensorA - .getDate(spA.getLineNumber()); - final AbsoluteDate dateB = lineSensorB - .getDate(spB.getLineNumber()); - final double pixelA = spA.getPixelNumber(); - final double pixelB = spB.getPixelNumber(); - final SpacecraftToObservedBody scToBodyA = ruggedA - .getScToBody(); - - final DerivativeStructure[] ilResult = distanceBetweenLOSDerivatives(lineSensorA, - dateA, - pixelA, - scToBodyA, - lineSensorB, - dateB, - pixelB, - generator); - - if (ilResult == null) { - // TODO - } else { - // extract the value - value.setEntry(l, ilResult[0].getValue()); - value.setEntry(l + 1, ilResult[1].getValue()); - // extract the Jacobian - final int[] orders = new int[selected - .getNbParams()]; - int m = 0; - for (final ParameterDriver driver : selected - .getDrivers()) { - final double scale = driver.getScale(); - orders[m] = 1; - jacobian.setEntry(l, m, - ilResult[0] - .getPartialDerivative(orders) * - scale); - - jacobian.setEntry(l + 1, m, - ilResult[1] - .getPartialDerivative(orders) * - scale); - - orders[m] = 0; - m++; - } - } - l += 2; // pass to the next evaluation - - } - } - - // distance result with Jacobian for all reference points - return new Pair<RealVector, RealMatrix>(value, jacobian); - - } catch (RuggedException re) { - throw new RuggedExceptionWrapper(re); - } catch (OrekitException oe) { - throw new OrekitExceptionWrapper(oe); - } - }; - - // set up the least squares problem - final LeastSquaresProblem problem = new LeastSquaresBuilder() - .lazyEvaluation(false).maxIterations(maxEvaluations) - .maxEvaluations(maxEvaluations).weight(new DiagonalMatrix(weight)).start(start) - .target(target).parameterValidator(validator).checker(checker) - .model(model).build(); - - // set up the optimizer - //final LeastSquaresOptimizer optimizer = new - // LevenbergMarquardtOptimizer(); - final LeastSquaresOptimizer optimizer = new GaussNewtonOptimizer() - .withDecomposition(GaussNewtonOptimizer.Decomposition.QR); - - // solve the least squares problem - return optimizer.optimize(problem); - - } catch (RuggedExceptionWrapper rew) { - throw rew.getException(); - } catch (OrekitExceptionWrapper oew) { - final OrekitException oe = oew.getException(); - throw new RuggedException(oe, oe.getSpecifier(), oe.getParts()); - } - } /** * Get the mean plane crossing finder for a sensor. @@ -1273,6 +1048,7 @@ public class Rugged { * unknown * @see #inverseLocation(String, GeodeticPoint, int, int) */ + public DerivativeStructure[] inverseLocationDerivatives(final String sensorName, final GeodeticPoint point, final int minLine, @@ -1369,8 +1145,6 @@ public class Rugged { } - - /** * Get transform from spacecraft to inertial frame. * diff --git a/src/main/java/org/orekit/rugged/api/RuggedBuilder.java b/src/main/java/org/orekit/rugged/api/RuggedBuilder.java index 74578feb190b3cd2479b09de628088c93c87c284..acda0f89ca844e2d3e02e76df1131045d96fd119 100644 --- a/src/main/java/org/orekit/rugged/api/RuggedBuilder.java +++ b/src/main/java/org/orekit/rugged/api/RuggedBuilder.java @@ -16,9 +16,6 @@ */ package org.orekit.rugged.api; -import org.hipparchus.exception.LocalizedCoreFormats; -import org.hipparchus.geometry.euclidean.threed.Rotation; -import org.hipparchus.geometry.euclidean.threed.Vector3D; import java.io.IOException; import java.io.InputStream; import java.io.ObjectInputStream; @@ -28,6 +25,9 @@ import java.util.ArrayList; import java.util.Collections; import java.util.List; +import org.hipparchus.exception.LocalizedCoreFormats; +import org.hipparchus.geometry.euclidean.threed.Rotation; +import org.hipparchus.geometry.euclidean.threed.Vector3D; import org.orekit.bodies.OneAxisEllipsoid; import org.orekit.errors.OrekitException; import org.orekit.frames.Frame; @@ -157,6 +157,9 @@ public class RuggedBuilder { /** Sensors. */ private final List<LineSensor> sensors; + /** Rugged Name. */ + private String name; + /** Create a non-configured builder. * <p> * The builder <em>must</em> be configured before calling the @@ -169,6 +172,7 @@ public class RuggedBuilder { constantElevation = Double.NaN; lightTimeCorrection = true; aberrationOfLightCorrection = true; + name = "Rugged"; } /** Set the reference ellipsoid. @@ -182,7 +186,7 @@ public class RuggedBuilder { * @see #getEllipsoid() */ public RuggedBuilder setEllipsoid(final EllipsoidId ellipsoidID, final BodyRotatingFrameId bodyRotatingFrameID) - throws RuggedException { + throws RuggedException { return setEllipsoid(selectEllipsoid(ellipsoidID, selectBodyRotatingFrame(bodyRotatingFrameID))); } @@ -196,7 +200,7 @@ public class RuggedBuilder { */ // CHECKSTYLE: stop HiddenField check public RuggedBuilder setEllipsoid(final OneAxisEllipsoid ellipsoid) - throws RuggedException { + throws RuggedException { // CHECKSTYLE: resume HiddenField check this.ellipsoid = new ExtendedEllipsoid(ellipsoid.getEquatorialRadius(), ellipsoid.getFlattening(), @@ -214,6 +218,22 @@ public class RuggedBuilder { return ellipsoid; } + /** Get the rugged name. + * @return the name + */ + public String getName() { + return name; + } + + + /** Set the rugged name. + * @param name the Rugged name + */ + public void setName(final String name) { + this.name = name; + } + + /** Set the algorithm to use for Digital Elevation Model intersection. * <p> * Note that some algorithms require specific other methods to be called too: @@ -417,7 +437,7 @@ public class RuggedBuilder { final CartesianDerivativesFilter pvFilter, final List<TimeStampedAngularCoordinates> quaternions, final int aInterpolationNumber, final AngularDerivativesFilter aFilter) - throws RuggedException { + throws RuggedException { return setTrajectory(selectInertialFrame(inertialFrame), positionsVelocities, pvInterpolationNumber, pvFilter, quaternions, aInterpolationNumber, aFilter); @@ -574,7 +594,7 @@ public class RuggedBuilder { * @see #storeInterpolator(OutputStream) */ public RuggedBuilder setTrajectoryAndTimeSpan(final InputStream storageStream) - throws RuggedException { + throws RuggedException { try { this.inertial = null; this.pvSample = null; @@ -638,7 +658,7 @@ public class RuggedBuilder { */ private void checkFramesConsistency() throws RuggedException { if (ellipsoid != null && scToBody != null && - !ellipsoid.getBodyFrame().getName().equals(scToBody.getBodyFrame().getName())) { + !ellipsoid.getBodyFrame().getName().equals(scToBody.getBodyFrame().getName())) { throw new RuggedException(RuggedMessages.FRAMES_MISMATCH_WITH_INTERPOLATOR_DUMP, ellipsoid.getBodyFrame().getName(), scToBody.getBodyFrame().getName()); } @@ -698,7 +718,7 @@ public class RuggedBuilder { final List<TimeStampedAngularCoordinates> quaternions, final int aInterpolationNumber, final AngularDerivativesFilter aFilter) - throws RuggedException { + throws RuggedException { return new SpacecraftToObservedBody(inertialFrame, bodyFrame, minDate, maxDate, tStep, overshootTolerance, positionsVelocities, pvInterpolationNumber, @@ -728,14 +748,14 @@ public class RuggedBuilder { final CartesianDerivativesFilter pvFilter, final AngularDerivativesFilter aFilter, final Propagator propagator) - throws RuggedException { + throws RuggedException { try { // extract position/attitude samples from propagator final List<TimeStampedPVCoordinates> positionsVelocities = - new ArrayList<TimeStampedPVCoordinates>(); + new ArrayList<TimeStampedPVCoordinates>(); final List<TimeStampedAngularCoordinates> quaternions = - new ArrayList<TimeStampedAngularCoordinates>(); + new ArrayList<TimeStampedAngularCoordinates>(); propagator.setMasterMode(interpolationStep, new OrekitFixedStepHandler() { /** {@inheritDoc} */ @@ -746,7 +766,7 @@ public class RuggedBuilder { /** {@inheritDoc} */ @Override public void handleStep(final SpacecraftState currentState, final boolean isLast) - throws OrekitException { + throws OrekitException { final AbsoluteDate date = currentState.getDate(); final PVCoordinates pv = currentState.getPVCoordinates(inertialFrame); final Rotation q = currentState.getAttitude().getRotation(); @@ -863,24 +883,24 @@ public class RuggedBuilder { * @exception RuggedException if data needed for some frame cannot be loaded */ private static Frame selectInertialFrame(final InertialFrameId inertialFrameId) - throws RuggedException { + throws RuggedException { try { // set up the inertial frame switch (inertialFrameId) { - case GCRF : - return FramesFactory.getGCRF(); - case EME2000 : - return FramesFactory.getEME2000(); - case MOD : - return FramesFactory.getMOD(IERSConventions.IERS_1996); - case TOD : - return FramesFactory.getTOD(IERSConventions.IERS_1996, true); - case VEIS1950 : - return FramesFactory.getVeis1950(); - default : - // this should never happen - throw RuggedException.createInternalError(null); + case GCRF : + return FramesFactory.getGCRF(); + case EME2000 : + return FramesFactory.getEME2000(); + case MOD : + return FramesFactory.getMOD(IERSConventions.IERS_1996); + case TOD : + return FramesFactory.getTOD(IERSConventions.IERS_1996, true); + case VEIS1950 : + return FramesFactory.getVeis1950(); + default : + // this should never happen + throw RuggedException.createInternalError(null); } } catch (OrekitException oe) { throw new RuggedException(oe, oe.getSpecifier(), oe.getParts().clone()); @@ -894,20 +914,20 @@ public class RuggedBuilder { * @exception RuggedException if data needed for some frame cannot be loaded */ private static Frame selectBodyRotatingFrame(final BodyRotatingFrameId bodyRotatingFrame) - throws RuggedException { + throws RuggedException { try { // set up the rotating frame switch (bodyRotatingFrame) { - case ITRF : - return FramesFactory.getITRF(IERSConventions.IERS_2010, true); - case ITRF_EQUINOX : - return FramesFactory.getITRFEquinox(IERSConventions.IERS_1996, true); - case GTOD : - return FramesFactory.getGTOD(IERSConventions.IERS_1996, true); - default : - // this should never happen - throw RuggedException.createInternalError(null); + case ITRF : + return FramesFactory.getITRF(IERSConventions.IERS_2010, true); + case ITRF_EQUINOX : + return FramesFactory.getITRFEquinox(IERSConventions.IERS_1996, true); + case GTOD : + return FramesFactory.getGTOD(IERSConventions.IERS_1996, true); + default : + // this should never happen + throw RuggedException.createInternalError(null); } } catch (OrekitException oe) { throw new RuggedException(oe, oe.getSpecifier(), oe.getParts().clone()); @@ -924,19 +944,19 @@ public class RuggedBuilder { // set up the ellipsoid switch (ellipsoidID) { - case GRS80 : - return new OneAxisEllipsoid(6378137.0, 1.0 / 298.257222101, bodyFrame); - case WGS84 : - return new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, - Constants.WGS84_EARTH_FLATTENING, - bodyFrame); - case IERS96 : - return new OneAxisEllipsoid(6378136.49, 1.0 / 298.25645, bodyFrame); - case IERS2003 : - return new OneAxisEllipsoid(6378136.6, 1.0 / 298.25642, bodyFrame); - default : - // this should never happen - throw RuggedException.createInternalError(null); + case GRS80 : + return new OneAxisEllipsoid(6378137.0, 1.0 / 298.257222101, bodyFrame); + case WGS84 : + return new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS, + Constants.WGS84_EARTH_FLATTENING, + bodyFrame); + case IERS96 : + return new OneAxisEllipsoid(6378136.49, 1.0 / 298.25645, bodyFrame); + case IERS2003 : + return new OneAxisEllipsoid(6378136.6, 1.0 / 298.25642, bodyFrame); + default : + // this should never happen + throw RuggedException.createInternalError(null); } } @@ -954,19 +974,19 @@ public class RuggedBuilder { // set up the algorithm switch (algorithmID) { - case DUVENHAGE : - return new DuvenhageAlgorithm(updater, maxCachedTiles, false); - case DUVENHAGE_FLAT_BODY : - return new DuvenhageAlgorithm(updater, maxCachedTiles, true); - case BASIC_SLOW_EXHAUSTIVE_SCAN_FOR_TESTS_ONLY : - return new BasicScanAlgorithm(updater, maxCachedTiles); - case CONSTANT_ELEVATION_OVER_ELLIPSOID : - return new ConstantElevationAlgorithm(constantElevation); - case IGNORE_DEM_USE_ELLIPSOID : - return new IgnoreDEMAlgorithm(); - default : - // this should never happen - throw RuggedException.createInternalError(null); + case DUVENHAGE : + return new DuvenhageAlgorithm(updater, maxCachedTiles, false); + case DUVENHAGE_FLAT_BODY : + return new DuvenhageAlgorithm(updater, maxCachedTiles, true); + case BASIC_SLOW_EXHAUSTIVE_SCAN_FOR_TESTS_ONLY : + return new BasicScanAlgorithm(updater, maxCachedTiles); + case CONSTANT_ELEVATION_OVER_ELLIPSOID : + return new ConstantElevationAlgorithm(constantElevation); + case IGNORE_DEM_USE_ELLIPSOID : + return new IgnoreDEMAlgorithm(); + default : + // this should never happen + throw RuggedException.createInternalError(null); } } @@ -991,7 +1011,7 @@ public class RuggedBuilder { } createInterpolatorIfNeeded(); return new Rugged(createAlgorithm(algorithmID, tileUpdater, maxCachedTiles, constantElevation), ellipsoid, - lightTimeCorrection, aberrationOfLightCorrection, scToBody, sensors); + lightTimeCorrection, aberrationOfLightCorrection, scToBody, sensors, name); } } diff --git a/src/tutorials/java/AffinagePleiades/GroundRefining.java b/src/tutorials/java/AffinagePleiades/GroundRefining.java index d0fcf3ecd0b4c49c28575e2ea2b1b6da085101cf..20c561587281501bf08e71cbb01def9420a1dacd 100644 --- a/src/tutorials/java/AffinagePleiades/GroundRefining.java +++ b/src/tutorials/java/AffinagePleiades/GroundRefining.java @@ -16,13 +16,12 @@ */ package AffinagePleiades; -import org.hipparchus.geometry.euclidean.threed.Vector3D; -import org.hipparchus.util.FastMath; - import java.io.File; -import java.util.Locale; import java.util.List; +import java.util.Locale; +import org.hipparchus.geometry.euclidean.threed.Vector3D; +import org.hipparchus.util.FastMath; import org.orekit.bodies.BodyShape; import org.orekit.bodies.GeodeticPoint; import org.orekit.data.DataProvidersManager; @@ -55,40 +54,40 @@ import org.orekit.utils.TimeStampedPVCoordinates; /** * Class for testing refining (fulcrum points study) - * with or without noisy measurements + * with or without noisy measurements * @author Jonathan Guinet * @author Lucie Labat-Allee * @see SensorToGroundMapping * @see GroundMeasureGenerator */ public class GroundRefining extends Refining { - + /** Pleiades viewing model */ PleiadesViewingModel pleiadesViewingModel; - + /** Orbit model */ OrbitModel orbitmodel; - + /** Sensor name */ String sensorName; - + /** Rugged instance */ Rugged rugged; - + /** Ground measurements */ GroundMeasureGenerator measures; - - + + /** * Constructor */ public GroundRefining() throws RuggedException, OrekitException { - + sensorName = "line"; pleiadesViewingModel = new PleiadesViewingModel(sensorName); orbitmodel = new OrbitModel(); } - + /** Main function * @param args */ @@ -97,17 +96,17 @@ public class GroundRefining extends Refining { try { // Initialize Orekit, assuming an orekit-data folder is in user home directory - // --------------------------------------------------------------------------- + // --------------------------------------------------------------------------- File home = new File(System.getProperty("user.home")); File orekitData = new File(home, "COTS/orekit-data"); DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(orekitData)); - - + + // Initialize refining context // --------------------------- GroundRefining refining = new GroundRefining(); - - + + // Sensor's definition: create Pleiades viewing model // -------------------------------------------------- System.out.format("**** Build Pleiades viewing model and orbit definition **** %n"); @@ -115,83 +114,84 @@ public class GroundRefining extends Refining { AbsoluteDate minDate = pleiadesViewingModel.getMinDate(); AbsoluteDate maxDate = pleiadesViewingModel.getMaxDate(); AbsoluteDate refDate = pleiadesViewingModel.getDatationReference(); - LineSensor lineSensor = pleiadesViewingModel.getLineSensor(); + LineSensor lineSensor = pleiadesViewingModel.getLineSensor(); + - // Satellite position, velocity and attitude: create an orbit model // ---------------------------------------------------------------- OrbitModel orbitmodel = refining.getOrbitmodel(); BodyShape earth = orbitmodel.createEarth(); NormalizedSphericalHarmonicsProvider gravityField = orbitmodel.createGravityField(); Orbit orbit = orbitmodel.createOrbit(gravityField.getMu(), refDate); - + // Nadir's pointing final double [] rollPoly = {0.0,0.0,0.0}; double[] pitchPoly = {0.0,0.0}; double[] yawPoly = {0.0,0.0,0.0}; orbitmodel.setLOFTransform(rollPoly, pitchPoly, yawPoly, minDate); - + // Satellite attitude List<TimeStampedAngularCoordinates> satelliteQList = orbitmodel.orbitToQ(orbit, earth, minDate.shiftedBy(-0.0), maxDate.shiftedBy(+0.0), 0.25); int nbQPoints = 2; - + // Position and velocities PVCoordinates PV = orbit.getPVCoordinates(earth.getBodyFrame()); List<TimeStampedPVCoordinates> satellitePVList = orbitmodel.orbitToPV(orbit, earth, minDate.shiftedBy(-0.0), maxDate.shiftedBy(+0.0), 0.25); int nbPVPoints = 8; - + // Convert frame and coordinates type in one call GeodeticPoint gp = earth.transform(PV.getPosition(), earth.getBodyFrame(), orbit.getDate()); - + System.out.format(Locale.US, "Geodetic Point at date %s : φ = %8.10f °, λ = %8.10f %n", orbit.getDate().toString(), FastMath.toDegrees(gp.getLatitude()), FastMath.toDegrees(gp.getLongitude())); - - + + // Rugged initialization // --------------------- System.out.format("\n**** Rugged initialization **** %n"); RuggedBuilder ruggedBuilder = new RuggedBuilder(); - + ruggedBuilder.addLineSensor(lineSensor); ruggedBuilder.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID); ruggedBuilder.setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF); ruggedBuilder.setTimeSpan(minDate,maxDate, 0.001, 5.0). - setTrajectory(InertialFrameId.EME2000, - satellitePVList,nbPVPoints, CartesianDerivativesFilter.USE_PV, - satelliteQList,nbQPoints, AngularDerivativesFilter.USE_R); - + setTrajectory(InertialFrameId.EME2000, + satellitePVList,nbPVPoints, CartesianDerivativesFilter.USE_PV, + satelliteQList,nbQPoints, AngularDerivativesFilter.USE_R); + ruggedBuilder.setName("Rugged_refining"); + refining.setRugged(ruggedBuilder.build()); - + // Compute ground sample distance (GSD) // ------------------------------------ double [] gsd = refining.computeGSD(lineSensor); - System.out.format("GSD - X: %2.2f Y: %2.2f **** %n", gsd[0], gsd[1]); - - - // Initialize disruptions: + System.out.format("GSD - X: %2.2f Y: %2.2f **** %n", gsd[0], gsd[1]); + + + // Initialize disruptions: // ----------------------- - + // Introduce rotations around instrument axes (roll and pitch translations, scale factor) - System.out.format("\n**** Add disruptions: roll and pitch rotations, scale factor **** %n"); + System.out.format("\n**** Add disruptions: roll and pitch rotations, scale factor **** %n"); double rollValue = FastMath.toRadians(-0.01); double pitchValue = FastMath.toRadians(0.02); double factorValue = 1.05; System.out.format("roll: %3.5f \tpitch: %3.5f \tfactor: %3.5f \n",rollValue, pitchValue, factorValue); - + // Apply disruptions on physical model refining.applyDisruptions(refining.getRugged(), refining.getSensorName(), rollValue, pitchValue, factorValue); - + // Generate measures (observations) from physical model disrupted // -------------------------------------------------------------- - + int lineSampling = 1000; - int pixelSampling = 1000; - + int pixelSampling = 1000; + // Noise definition final Noise noise = new Noise(0,3); /* distribution: gaussian(0), vector dimension:3 */ final double[] mean = {0.0,0.0,0.0}; /* {lat mean, long mean, alt mean} */ @@ -199,38 +199,38 @@ public class GroundRefining extends Refining { noise.setMean(mean); noise.setStandardDeviation(standardDeviation); - - GroundMeasureGenerator measures = refining.generateNoisyPoints(lineSampling, pixelSampling, - refining.getRugged(), refining.getSensorName(), + + GroundMeasureGenerator measures = refining.generateNoisyPoints(lineSampling, pixelSampling, + refining.getRugged(), refining.getSensorName(), refining.getPleiadesViewingModel().getDimension(), noise); - -// // To test with measures without noise -// GroundMeasureGenerator measures = refining.generatePoints(lineSampling, pixelSampling, -// refining.getRugged(), refining.getSensorName(), -// refining.getPleiadesViewingModel().getDimension()); - + + // // To test with measures without noise + // GroundMeasureGenerator measures = refining.generatePoints(lineSampling, pixelSampling, + // refining.getRugged(), refining.getSensorName(), + // refining.getPleiadesViewingModel().getDimension()); + refining.setMeasures(measures); - - + + // Compute ground truth residues // ----------------------------- System.out.format("\n**** Ground truth residuals **** %n"); refining.computeMetrics(measures.getGroundMapping(), refining.getRugged(), false); - - + + // Initialize physical model without disruptions // --------------------------------------------- System.out.format("\n**** Initialize physical model without disruptions: reset Roll/Pitch/Factor **** %n"); refining.resetModel(refining.getRugged(), refining.getSensorName(), true); - - + + // Compute initial residues // ------------------------ System.out.format("\n**** Initial Residuals **** %n"); refining.computeMetrics(measures.getGroundMapping(), refining.getRugged(), false); - - + + // Start optimization // ------------------ System.out.format("\n**** Start optimization **** %n"); @@ -239,26 +239,26 @@ public class GroundRefining extends Refining { double convergenceThreshold = 1e-14; refining.optimization(maxIterations, convergenceThreshold, measures.getGroundMapping(), refining.getRugged()); - - + + // Check estimated values // ---------------------- System.out.format("\n**** Check parameters ajustement **** %n"); refining.paramsEstimation(refining.getRugged(), refining.getSensorName(), rollValue, pitchValue, factorValue); - - + + // Compute statistics // ------------------ System.out.format("\n**** Compute Statistics **** %n"); - + // Residues computed in meters refining.computeMetrics(measures.getGroundMapping(), refining.getRugged(), false); - + // Residues computed in degrees refining.computeMetrics(measures.getGroundMapping(), refining.getRugged(), true); - - + + } catch (OrekitException oe) { System.err.println(oe.getLocalizedMessage()); System.exit(1); @@ -276,7 +276,7 @@ public class GroundRefining extends Refining { return pleiadesViewingModel; } - + /** * @param pleiadesViewingModel the pleiadesViewingModel to set */ @@ -284,7 +284,7 @@ public class GroundRefining extends Refining { this.pleiadesViewingModel = pleiadesViewingModel; } - + /** * @return the orbitmodel */ @@ -292,7 +292,7 @@ public class GroundRefining extends Refining { return orbitmodel; } - + /** * @param orbitmodel the orbitmodel to set */ @@ -300,7 +300,7 @@ public class GroundRefining extends Refining { this.orbitmodel = orbitmodel; } - + /** * @return the sensorName */ @@ -308,7 +308,7 @@ public class GroundRefining extends Refining { return sensorName; } - + /** * @return the rugged */ @@ -316,15 +316,15 @@ public class GroundRefining extends Refining { return rugged; } - + /** * @param rugged the rugged to set */ public void setRugged(Rugged rugged) { this.rugged = rugged; } - - + + /** * @param measure the measure to set */ @@ -332,45 +332,45 @@ public class GroundRefining extends Refining { this.measures = measures; } - + /** Estimate ground sample distance (GSD) * @param LineSensor line sensor * @return GSD */ private double[] computeGSD(final LineSensor lineSensor) throws RuggedException { - + // Get position Vector3D position = lineSensor.getPosition(); // This returns a zero vector since we set the relative position of the sensor w.r.T the satellite to 0. - + // Get upper left geodetic point AbsoluteDate firstLineDate = lineSensor.getDate(0); Vector3D los = lineSensor.getLOS(firstLineDate,0); GeodeticPoint upLeftPoint = rugged.directLocation(firstLineDate, position, los); los = lineSensor.getLOS(firstLineDate,pleiadesViewingModel.dimension-1); - + // Get center geodetic point AbsoluteDate lineDate = lineSensor.getDate(pleiadesViewingModel.dimension/2); los = lineSensor.getLOS(lineDate,pleiadesViewingModel.dimension/2); -// GeodeticPoint centerPoint = rugged.directLocation(lineDate, position, los); -// System.out.format(Locale.US, "center geodetic position : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n", -// FastMath.toDegrees(centerPoint.getLatitude()), -// FastMath.toDegrees(centerPoint.getLongitude()),centerPoint.getAltitude()); - - // Get upper right geodetic point + // GeodeticPoint centerPoint = rugged.directLocation(lineDate, position, los); + // System.out.format(Locale.US, "center geodetic position : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n", + // FastMath.toDegrees(centerPoint.getLatitude()), + // FastMath.toDegrees(centerPoint.getLongitude()),centerPoint.getAltitude()); + + // Get upper right geodetic point int pixelPosition = pleiadesViewingModel.dimension-1; los = lineSensor.getLOS(firstLineDate,pixelPosition); GeodeticPoint upperRight = rugged.directLocation(firstLineDate, position, los); - - // Get lower left geodetic point + + // Get lower left geodetic point AbsoluteDate lineDate_y = lineSensor.getDate(pleiadesViewingModel.dimension-1); los = lineSensor.getLOS(lineDate_y,0); GeodeticPoint lowerLeft = rugged.directLocation(lineDate_y, position, los); - + double gsdX = DistanceTools.computeDistanceInMeter(upLeftPoint.getLongitude(), upLeftPoint.getLatitude(),upperRight.getLongitude() , upperRight.getLatitude())/pleiadesViewingModel.dimension; double gsdY = DistanceTools.computeDistanceInMeter(upLeftPoint.getLongitude(), upLeftPoint.getLatitude(),lowerLeft.getLongitude() , lowerLeft.getLatitude())/pleiadesViewingModel.dimension; - + double [] gsd = {gsdX, gsdY}; return gsd; } -} - +} + diff --git a/src/tutorials/java/AffinagePleiades/Refining.java b/src/tutorials/java/AffinagePleiades/Refining.java index 73c5a4921acbe734694b18ccf9c203efe6698825..7f93ae890b5f9d95fdec2d6e128862b6efa8d74a 100644 --- a/src/tutorials/java/AffinagePleiades/Refining.java +++ b/src/tutorials/java/AffinagePleiades/Refining.java @@ -328,7 +328,7 @@ public class Refining { observables.addGroundMapping(measures); AdjustmentContext adjustmentContext = new AdjustmentContext(viewingModel, observables); - Optimum optimum = adjustmentContext.estimateFreeParameters("Rugged",maxIterations,convergenceThreshold); + Optimum optimum = adjustmentContext.estimateFreeParameters("Rugged_refining",maxIterations,convergenceThreshold); // Print statistics System.out.format("Max value: %3.6e %n",optimum.getResiduals().getMaxValue());