Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found

Target

Select target project
  • orekit/rugged
  • sdinot/rugged
  • yzokras/rugged
  • youngcle/rugged-mod
4 results
Show changes
Hide whitespace changes
Inline Side-by-side
Showing
with 2123 additions and 1583 deletions
src/test/resources/replay/replay-inverse-loc-01.txt
View file @ acf902ea
......@@ -8,7 +8,7 @@ DEM cell: t0 latIndex 122 lonIndex 99 elevation 1.778606495649060e+03
DEM cell: t0 latIndex 122 lonIndex 100 elevation 1.809481571226391e+03
DEM cell: t0 latIndex 123 lonIndex 99 elevation 1.772149770526270e+03
DEM cell: t0 latIndex 123 lonIndex 100 elevation 1.804520334312959e+03
inverse location: sensorName s0 latitude -3.930524233859007e-01 longitude 2.502596047891660e+00 elevation 1.785429914217195e+03 minLine 0 maxLine 2000 lightTime true aberration true
inverse location: sensorName s0 latitude -3.930524233859007e-01 longitude 2.502596047891660e+00 elevation 1.785429914217195e+03 minLine 0 maxLine 2000 lightTime true aberration true refraction false
sensor: sensorName s0 nbPixels 2000 position 1.500000000000000e+00 0.000000000000000e+00 -2.000000000000000e-01
sensor mean plane: sensorName s0 minLine 0 maxLine 2000 maxEval 50 accuracy 1.000000000000000e-02 normal 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 cachedResults 0
sensor datation: sensorName s0 lineNumber 1.000000000000000e+03 date 2012-01-01T12:30:00.00000000000000Z
......
src/test/resources/replay/replay-inverse-loc-02.txt
View file @ acf902ea
......@@ -9,7 +9,7 @@ DEM cell: t0 latIndex 123 lonIndex 102 elevation 1.870339491485550e+03
DEM cell: t0 latIndex 124 lonIndex 101 elevation 1.832853313250435e+03
DEM cell: t0 latIndex 124 lonIndex 102 elevation 1.866819199623890e+03
sensor: sensorName s0 nbPixels 2000 position 1.500000000000000e+00 0.000000000000000e+00 -2.000000000000000e-01
inverse location: sensorName s0 latitude -3.930305027145593e-01 longitude 2.502739578717502e+00 elevation 1.853045504720662e+03 minLine 0 maxLine 2000 lightTime true aberration true
inverse location: sensorName s0 latitude -3.930305027145593e-01 longitude 2.502739578717502e+00 elevation 1.853045504720662e+03 minLine 0 maxLine 2000 lightTime true aberration true refraction false
sensor mean plane: sensorName s0 minLine 0 maxLine 2000 maxEval 50 accuracy 1.000000000000000e-02 normal 1.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00 cachedResults 5 lineNumber 1.753079999999941e+03 date 2012-01-01T12:30:01.12961999999991Z target -4.733584770718903e+06 3.515842173828561e+06 -2.428300713188566e+06 targetDirection 4.440892098500626e-16 -7.496646388847980e-01 6.618178973144536e-01 -1.357198500050715e-08 -1.198160098465607e-08 -1.357198500050715e-08 lineNumber 1.753079999999936e+03 date 2012-01-01T12:30:01.12961999999990Z target -4.733634454810624e+06 3.515791377291927e+06 -2.428292161157289e+06 targetDirection 1.221245327087672e-15 -7.496479456645485e-01 6.618368058373009e-01 -1.357209409884309e-08 -1.198230643978911e-08 -1.357209409884309e-08 lineNumber 1.753080000000098e+03 date 2012-01-01T12:30:01.12962000000015Z target -4.733684171891571e+06 3.515740571824192e+06 -2.428283613116900e+06 targetDirection 2.664535259100376e-15 -7.496312519673904e-01 6.618557139390744e-01 -1.357220357852602e-08 -1.198301226854942e-08 -1.357220357852602e-08 lineNumber 1.753080000000016e+03 date 2012-01-01T12:30:01.12962000000002Z target -4.733734047015606e+06 3.515689711809538e+06 -2.428275079609979e+06 targetDirection 1.110223024625157e-15 -7.496145577933011e-01 6.618746216197999e-01 -1.357231495356307e-08 -1.198371980774984e-08 -1.357231495356307e-08 lineNumber 1.753080000000124e+03 date 2012-01-01T12:30:01.12962000000019Z target -4.733783920649958e+06 3.515638855451895e+06 -2.428266547191616e+06 targetDirection 1.665334536937735e-15 -7.495978631423109e-01 6.618935288794415e-01 -1.357242629275803e-08 -1.198442735245515e-08 -1.357242629275803e-08
sensor datation: sensorName s0 lineNumber 1.000000000000000e+03 date 2012-01-01T12:30:00.00000000000000Z
span: minDate 2012-01-01T12:29:57.50000000000000Z maxDate 2012-01-01T12:30:02.50000000000000Z tStep 1.000000000000000e-03 tolerance 5.000000000000000e+00 inertialFrame EME2000
......@@ -28,14 +28,14 @@ sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12961999999997Z pixelNumber 5
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12961999999997Z pixelNumber 6 los 0.000000000000000e+00 -7.496980238944815e-01 6.617800790056386e-01
sensor datation: sensorName s0 lineNumber 1.753079999999893e+03 date 2012-01-01T12:30:01.12961999999984Z
inverse location result: lineNumber 1.753079999999893e+03 pixelNumber 5.000000000428289e+00
inverse location: sensorName s0 latitude -3.930323254721546e-01 longitude 2.502727645691323e+00 elevation 1.847490708944211e+03 minLine 0 maxLine 2000 lightTime true aberration true
inverse location: sensorName s0 latitude -3.930323254721546e-01 longitude 2.502727645691323e+00 elevation 1.847490708944211e+03 minLine 0 maxLine 2000 lightTime true aberration true refraction false
sensor datation: sensorName s0 lineNumber 1.753080000000180e+03 date 2012-01-01T12:30:01.12962000000027Z
sensor rate: sensorName s0 lineNumber 1.753080000000180e+03 rate 6.666666666666666e+02
sensor datation: sensorName s0 lineNumber 1.753080000000150e+03 date 2012-01-01T12:30:01.12962000000022Z
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12962000000022Z pixelNumber 7 los 0.000000000000000e+00 -7.497147156838981e-01 6.617611692196904e-01
sensor datation: sensorName s0 lineNumber 1.753080000000120e+03 date 2012-01-01T12:30:01.12962000000018Z
inverse location result: lineNumber 1.753080000000120e+03 pixelNumber 6.000000000570598e+00
inverse location: sensorName s0 latitude -3.930341476421113e-01 longitude 2.502715716159622e+00 elevation 1.841953904196975e+03 minLine 0 maxLine 2000 lightTime true aberration true
inverse location: sensorName s0 latitude -3.930341476421113e-01 longitude 2.502715716159622e+00 elevation 1.841953904196975e+03 minLine 0 maxLine 2000 lightTime true aberration true refraction false
sensor datation: sensorName s0 lineNumber 1.753080000000450e+03 date 2012-01-01T12:30:01.12962000000067Z
sensor rate: sensorName s0 lineNumber 1.753080000000450e+03 rate 6.666666666666666e+02
sensor LOS: sensorName s0 date 2012-01-01T12:30:01.12962000000027Z pixelNumber 8 los 0.000000000000000e+00 -7.497314069963494e-01 6.617422590127325e-01
......
src/test/resources/replay/replay-inverse-loc-03.txt
View file @ acf902ea
# Rugged library dump file, created on 2015-07-23T14:25:51Z
# all units are SI units (m, m/s, rad ...)
sensor: sensorName s0 nbPixels 2552 position 0.000000000000000e+00 0.000000000000000e+00 0.000000000000000e+00
inverse location: sensorName s0 latitude 1.426113542156674e+00 longitude -8.868185262195961e-01 elevation 3.102209319249422e+01 minLine -23040 maxLine 39851 lightTime false aberration false
inverse location: sensorName s0 latitude 1.426113542156674e+00 longitude -8.868185262195961e-01 elevation 3.102209319249422e+01 minLine -23040 maxLine 39851 lightTime false aberration false refraction false
sensor mean plane: sensorName s0 minLine -23040 maxLine 39851 maxEval 50 accuracy 1.000000000000000e-02 normal 9.994827482027009e-01 -2.652758464258975e-02 1.818029972819145e-02 cachedResults 6 lineNumber 2.496503931691230e+04 date 2015-07-07T18:40:12.48786757028467Z target 5.829885691015603e+05 -7.150938152185429e+05 6.289925119642745e+06 targetDirection -1.574993670955663e-02 8.911384689865436e-02 9.958969132317702e-01 -2.069121568391233e-07 2.118875364268554e-08 -2.069121568391233e-07 lineNumber 2.498990311831761e+04 date 2015-07-07T18:40:12.52680428328538Z target 5.827831035018093e+05 -7.152819110623585e+05 6.290040934506435e+06 targetDirection -1.574653262138034e-02 8.923613956787135e-02 9.958860166231007e-01 -2.106025594148216e-07 6.293169679808365e-08 -2.106025594148216e-07 lineNumber 2.501446735659665e+04 date 2015-07-07T18:40:12.56527188043036Z target 5.825914009565986e+05 -7.154869255677856e+05 6.290306437975926e+06 targetDirection -1.574243698611968e-02 8.937839743168323e-02 9.958733241483490e-01 -2.139580511938570e-07 1.003165269577007e-07 -2.139580511938570e-07 lineNumber 2.503906484840334e+04 date 2015-07-07T18:40:12.60379155259964Z target 5.823982602071301e+05 -7.156902096509174e+05 6.290556499572599e+06 targetDirection -1.573844749711850e-02 8.951962023030334e-02 9.958607025403716e-01 -2.164310151145369e-07 1.276529299636964e-07 -2.164310151145369e-07 lineNumber 2.506393183649445e+04 date 2015-07-07T18:40:12.64273325595032Z target 5.821928312311379e+05 -7.158784186376551e+05 6.290673076115117e+06 targetDirection -1.573492577246402e-02 8.964463942750402e-02 9.958495120914844e-01 -2.189036239175299e-07 1.552935576937779e-07 -2.189036239175299e-07 lineNumber 2.508894303441162e+04 date 2015-07-07T18:40:12.68190079188860Z target 5.819807540563018e+05 -7.160584640312950e+05 6.290717393595084e+06 targetDirection -1.573177124337613e-02 8.976098395370098e-02 9.958390819470831e-01 -2.207684005657831e-07 1.758718590383612e-07 -2.207684005657831e-07
sensor datation: sensorName s0 lineNumber 8.405500000000000e+03 date 2015-07-07T18:39:46.55562900000000Z
span: minDate 2015-07-07T18:38:55.00000000000000Z maxDate 2015-07-07T18:40:35.80000000000000Z tStep 1.000000000000000e-01 tolerance 1.000000000000000e+01 inertialFrame EME2000
......
src/tutorials/java/AffinagePleiades/AffinageRugged.java deleted 100644 → 0
View file @ 950e0db2
/* Copyright 2013-2016 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package AffinagePleiades;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum;
import java.io.File;
import java.util.Locale;
import java.util.Collections;
import java.io.FileWriter;
import java.io.StringWriter;
import java.io.PrintWriter;
import java.io.IOException;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitExceptionWrapper;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.orbits.Orbit;
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.PVCoordinates;
/**
* Parameter refining context
* @author Jonathan Guinet
* @author Lucie LabatAllee
*/
public class AffinageRugged {
public static void main(String[] args) {
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, "workspace/data/orekit-data");
//File orekitData = new File(home, "COTS/orekit-data");
DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(orekitData));
//CreateOrbit
OrbitModel orbitmodel = new OrbitModel();
BodyShape earth = orbitmodel.createEarth();
NormalizedSphericalHarmonicsProvider gravityField = orbitmodel.createGravityField();
//create Pleiades Viewing Model
PleiadesViewingModel pleiadesViewingModel = new PleiadesViewingModel();
AbsoluteDate minDate = pleiadesViewingModel.getMinDate();
AbsoluteDate maxDate = pleiadesViewingModel.getMaxDate();
AbsoluteDate refDate = pleiadesViewingModel.getDatationReference();
Orbit orbit = orbitmodel.createOrbit(gravityField.getMu(), refDate);
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);
PVCoordinates PV = orbit.getPVCoordinates(earth.getBodyFrame());
GeodeticPoint gp = earth.transform(PV.getPosition(), earth.getBodyFrame(), orbit.getDate());
System.out.format(" **** Orbit Definition **** %n");
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()));
System.out.format(" **** Build Viewing Model **** %n");
// Build Rugged ...
RuggedBuilder ruggedBuilder = new RuggedBuilder();
LineSensor lineSensor = pleiadesViewingModel.getLineSensor();
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,
orbitmodel.orbitToPV(orbit, earth, minDate.shiftedBy(-0.0), maxDate.shiftedBy(+0.0), 0.25),
8, CartesianDerivativesFilter.USE_PV,
orbitmodel.orbitToQ(orbit, earth, minDate.shiftedBy(-0.0), maxDate.shiftedBy(+0.0), 0.25),
2, AngularDerivativesFilter.USE_R);
Rugged rugged = ruggedBuilder.build();
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.
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);
//GeodeticPoint upRightPoint = rugged.directLocation(firstLineDate, position, los);
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());
/*try {
FileWriter out = new FileWriter("geodetic_position.txt");
out.write("LATITUDE, LONGITUDE, LINE, PIX \n");
for (double line = 0; line < pleiadesViewingModel.dimension; line += 1000) {
AbsoluteDate date = lineSensor.getDate(line);
for (int pixel = 0; pixel < lineSensor.getNbPixels(); pixel += 1000) {
GeodeticPoint gp2 = rugged.directLocation(date, lineSensor.getPosition(),
lineSensor.getLOS(date, pixel));
String result = String.format("%8.10f, %8.10f, %2.1f, %2.1f \n ",FastMath.toDegrees(gp2.getLatitude()),
FastMath.toDegrees(gp2.getLongitude()), line,(double) pixel);
//System.out.println(result);
out.write(result);
}
}
out.close();
}
catch (IOException e) {
System.err.println(e.getLocalizedMessage());
System.exit(1);
}*/
int pixelPosition = pleiadesViewingModel.dimension-1;
los = lineSensor.getLOS(firstLineDate,pixelPosition);
GeodeticPoint upperRight = rugged.directLocation(firstLineDate, position, los);
AbsoluteDate lineDate_y = lineSensor.getDate(pleiadesViewingModel.dimension-1);
los = lineSensor.getLOS(lineDate_y,0);
GeodeticPoint lowerLeft = rugged.directLocation(lineDate_y, position, los);
double GSD_X = DistanceTools.computeDistanceRad(upLeftPoint.getLongitude(), upLeftPoint.getLatitude(),upperRight.getLongitude() , upperRight.getLatitude())/pleiadesViewingModel.dimension;
double GSD_Y = DistanceTools.computeDistanceRad(upLeftPoint.getLongitude(), upLeftPoint.getLatitude(),lowerLeft.getLongitude() , lowerLeft.getLatitude())/pleiadesViewingModel.dimension;
System.out.format(" GSD X %2.2f Y %2.2f **** %n", GSD_X, GSD_Y);
System.out.format(" **** Add roll and pitch values **** %n");
double rollValue = FastMath.toRadians(-0.01);
double pitchValue = FastMath.toRadians(0.02);
double factorValue = 1.05;
System.out.format("roll : %3.5f pitch : %3.5f factor : %3.5f \n",rollValue,pitchValue,factorValue);
rugged.
getLineSensor("line").
getParametersDrivers().
filter(driver -> driver.getName().equals("roll")).
findFirst().get().setValue(rollValue);
rugged.
getLineSensor("line").
getParametersDrivers().
filter(driver -> driver.getName().equals("pitch")).
findFirst().get().setValue(pitchValue);
rugged.
getLineSensor("line").
getParametersDrivers().
filter(driver -> driver.getName().equals("factor")).
findFirst().get().setValue(factorValue);
System.out.format(" **** Generate Measures **** %n");
MeasureGenerator measure = new MeasureGenerator(pleiadesViewingModel,rugged);
int lineSampling = 1000;
int pixelSampling = 1000;
//measure.CreateMeasure(lineSampling, pixelSampling);
//double pixErr = 0.0;
//double altErr = 0.0;
final double[] pixErr = new double[4];
pixErr[0] = 0; // lat mean
pixErr[1] = 0.0; // lat std
pixErr[2] = 0; // lon mean
pixErr[3] = 0.0; // lon std
final double[] altErr = new double[2];
altErr[0] = 0.0; //mean
altErr[1] = 0.0; //std
measure.CreateNoisyMeasure(lineSampling, pixelSampling,pixErr,altErr); // Test with noisy measures
System.out.format("nb TiePoints %d %n", measure.getMeasureCount());
System.out.format(" **** Initial Residuals **** %n");
LocalisationMetrics gtResiduals = new LocalisationMetrics(measure.getMapping(),rugged, false);
System.out.format("gt residuals max : %3.4f mean %3.4f meters %n",gtResiduals.getMaxResidual(),gtResiduals.getMeanResidual());
System.out.format(" **** Reset Roll/Pitch/Factor **** %n");
rugged.
getLineSensor(pleiadesViewingModel.getSensorName()).
getParametersDrivers().
filter(driver -> driver.getName().equals("roll") || driver.getName().equals("pitch")).
forEach(driver -> {
try {
driver.setSelected(true);
driver.setValue(0.0);
} catch (OrekitException e) {
throw new OrekitExceptionWrapper(e);
}
});
rugged.
getLineSensor(pleiadesViewingModel.getSensorName()).
getParametersDrivers().
filter(driver -> driver.getName().equals("factor")).
forEach(driver -> {
try {
driver.setSelected(true);
driver.setValue(1.0); // default value: no Z scale factor applied
} catch (OrekitException e) {
throw new OrekitExceptionWrapper(e);
}
});
System.out.format(" **** Initial Residuals **** %n");
LocalisationMetrics initlialResiduals = new LocalisationMetrics(measure.getMapping(),rugged, false);
System.out.format("residuals max : %3.4f mean %3.4f meters %n",initlialResiduals.getMaxResidual(),initlialResiduals.getMeanResidual());
System.out.format(" **** Start optimization **** %n");
// perform parameters estimation
int maxIterations = 100;
double convergenceThreshold = 1e-14;
System.out.format("iterations max %d convergence threshold %3.6e \n",maxIterations, convergenceThreshold);
Optimum optimum = rugged.estimateFreeParameters(Collections.singletonList(measure.getMapping()), maxIterations,convergenceThreshold);
System.out.format("max value %3.6e %n",optimum.getResiduals().getMaxValue());
System.out.format(" Optimization performed in %d iterations \n",optimum.getEvaluations());
// check estimated values
double estimatedRoll = rugged.getLineSensor(pleiadesViewingModel.getSensorName()).
getParametersDrivers().
filter(driver -> driver.getName().equals("roll")).
findFirst().get().getValue();
double rollError = (estimatedRoll - rollValue);
System.out.format("Estimated roll %3.5f roll error %3.6e %n", estimatedRoll, rollError);
double estimatedPitch = rugged.getLineSensor(pleiadesViewingModel.getSensorName()).
getParametersDrivers().
filter(driver -> driver.getName().equals("pitch")).
findFirst().get().getValue();
double pitchError = (estimatedPitch - pitchValue);
System.out.format("Estimated pitch %3.5f pitch error %3.6e %n ", estimatedPitch, pitchError);
double estimatedFactor = rugged.getLineSensor(pleiadesViewingModel.getSensorName()).
getParametersDrivers().
filter(driver -> driver.getName().equals("factor")).
findFirst().get().getValue();
double factorError = (estimatedFactor - factorValue);
System.out.format("Estimated factor %3.5f factor error %3.6e %n ", estimatedFactor, factorError);
System.out.format(" **** Compute Statistics **** %n");
LocalisationMetrics localisationResiduals = new LocalisationMetrics(measure.getMapping(),rugged, false);
System.out.format("residuals max : %3.4e mean %3.4e meters %n",localisationResiduals.getMaxResidual(),localisationResiduals.getMeanResidual());
LocalisationMetrics localisationResidualsDeg = new LocalisationMetrics(measure.getMapping(),rugged, true);
System.out.format("residuals max : %3.4e deg mean %3.4e deg %n",localisationResidualsDeg.getMaxResidual(),localisationResidualsDeg.getMeanResidual());
} catch (OrekitException oe) {
System.err.println(oe.getLocalizedMessage());
System.exit(1);
} catch (RuggedException re) {
System.err.println(re.getLocalizedMessage());
System.exit(1);
}
}
}
src/tutorials/java/AffinagePleiades/AffinageRuggedLiaison.java deleted 100644 → 0
View file @ 950e0db2
/* Copyright 2013-2016 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package AffinagePleiades;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum;
import java.io.File;
import java.util.Locale;
import java.util.Collections;
import java.io.FileWriter;
import java.io.StringWriter;
import java.io.PrintWriter;
import java.io.IOException;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.errors.OrekitExceptionWrapper;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.orbits.Orbit;
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.PVCoordinates;
/**
* Parameter refining context
* @author Jonathan Guinet
* @author Lucie LabatAllee
*/
public class AffinageRuggedLiaison {
/**
* @param args
*/
public static void main(String[] args) {
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, "workspace/data/orekit-data");
//File orekitData = new File(home, "COTS/orekit-data");
DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(orekitData));
//create 2 Pleiades Viewing Model
final String SensorNameA = "SensorA";
final double incidenceAngleA = -0.0;
final String dateA = "2016-01-01T11:59:50.0";
PleiadesViewingModel pleiadesViewingModelA = new PleiadesViewingModel(SensorNameA,incidenceAngleA,dateA);
final AbsoluteDate minDateA = pleiadesViewingModelA.getMinDate();
final AbsoluteDate maxDateA = pleiadesViewingModelA.getMaxDate();
final AbsoluteDate refDateA = pleiadesViewingModelA.getDatationReference();
//CreateOrbitA
OrbitModel orbitmodelA = new OrbitModel();
BodyShape earthA = orbitmodelA.createEarth();
NormalizedSphericalHarmonicsProvider gravityFieldA = orbitmodelA.createGravityField();
Orbit orbitA = orbitmodelA.createOrbit(gravityFieldA.getMu(), refDateA);
final double [] rollPoly = {0.0,0.0,0.0};
double[] pitchPoly = {0.0,0.0};
double[] yawPoly = {0.0,0.0,0.0};
orbitmodelA.setLOFTransform(rollPoly, pitchPoly, yawPoly, minDateA);
PVCoordinates PVA = orbitA.getPVCoordinates(earthA.getBodyFrame());
GeodeticPoint gpA = earthA.transform(PVA.getPosition(), earthA.getBodyFrame(), orbitA.getDate());
System.out.format(Locale.US, "Geodetic Point at date %s : φ = %8.10f °, λ = %8.10f %n",orbitA.getDate().toString(),
FastMath.toDegrees(gpA.getLatitude()),
FastMath.toDegrees(gpA.getLongitude()));
//create 2 Pleiades Viewing Model
final String SensorNameB = "SensorB";
final double incidenceAngleB = +0.0;
final String dateB = "2016-01-01T12:00:10.0";
PleiadesViewingModel pleiadesViewingModelB = new PleiadesViewingModel(SensorNameB,incidenceAngleB,dateB);
final AbsoluteDate minDateB = pleiadesViewingModelB.getMinDate();
final AbsoluteDate maxDateB = pleiadesViewingModelB.getMaxDate();
final AbsoluteDate refDateB = pleiadesViewingModelB.getDatationReference();
//CreateOrbitA
OrbitModel orbitmodelB = new OrbitModel();
BodyShape earthB = orbitmodelB.createEarth();
NormalizedSphericalHarmonicsProvider gravityFieldB = orbitmodelB.createGravityField();
Orbit orbitB = orbitmodelB.createOrbit(gravityFieldB.getMu(), refDateB);
orbitmodelB.setLOFTransform(rollPoly, pitchPoly, yawPoly, minDateB);
PVCoordinates PVB = orbitB.getPVCoordinates(earthB.getBodyFrame());
GeodeticPoint gpB = earthB.transform(PVB.getPosition(), earthB.getBodyFrame(), orbitB.getDate());
System.out.format(Locale.US, "Geodetic Point at date %s : φ = %8.10f °, λ = %8.10f %n",orbitA.getDate().toString(),
FastMath.toDegrees(gpB.getLatitude()),
FastMath.toDegrees(gpB.getLongitude()));
// Build Rugged A and B
RuggedBuilder ruggedBuilderA = new RuggedBuilder();
LineSensor lineSensorA = pleiadesViewingModelA.getLineSensor();
ruggedBuilderA.addLineSensor(lineSensorA);
ruggedBuilderA.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID);
ruggedBuilderA.setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF);
ruggedBuilderA.setTimeSpan(minDateA,maxDateA, 0.001, 5.0).
setTrajectory(InertialFrameId.EME2000,
orbitmodelA.orbitToPV(orbitA, earthA, minDateA.shiftedBy(-0.0), maxDateA.shiftedBy(+0.0), 0.25),
8, CartesianDerivativesFilter.USE_PV,
orbitmodelA.orbitToQ(orbitA, earthA, minDateA.shiftedBy(-0.0), maxDateA.shiftedBy(+0.0), 0.25),
2, AngularDerivativesFilter.USE_R);
Rugged ruggedA = ruggedBuilderA.build();
RuggedBuilder ruggedBuilderB = new RuggedBuilder();
LineSensor lineSensorB = pleiadesViewingModelB.getLineSensor();
ruggedBuilderB.addLineSensor(lineSensorB);
ruggedBuilderB.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID);
ruggedBuilderB.setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF);
ruggedBuilderB.setTimeSpan(minDateB,maxDateB, 0.001, 5.0).
setTrajectory(InertialFrameId.EME2000,
orbitmodelB.orbitToPV(orbitB, earthB, minDateB.shiftedBy(-0.0), maxDateB.shiftedBy(+0.0), 0.25),
8, CartesianDerivativesFilter.USE_PV,
orbitmodelB.orbitToQ(orbitB, earthB, minDateB.shiftedBy(-0.0), maxDateB.shiftedBy(+0.0), 0.25),
2, AngularDerivativesFilter.USE_R);
Rugged ruggedB = ruggedBuilderB.build();
Vector3D positionA = lineSensorA.getPosition(); // This returns a zero vector since we set the relative position of the sensor w.r.T the satellite to 0.
AbsoluteDate lineDateA = lineSensorA.getDate(pleiadesViewingModelA.dimension/2);
Vector3D losA = lineSensorA.getLOS(lineDateA,pleiadesViewingModelA.dimension/2);
GeodeticPoint centerPointA = ruggedA.directLocation(lineDateA, positionA, losA);
System.out.format(Locale.US, "center geodetic position : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(centerPointA.getLatitude()),
FastMath.toDegrees(centerPointA.getLongitude()),centerPointA.getAltitude());
Vector3D positionB = lineSensorB.getPosition(); // This returns a zero vector since we set the relative position of the sensor w.r.T the satellite to 0.
AbsoluteDate lineDateB = lineSensorB.getDate(pleiadesViewingModelB.dimension/2);
Vector3D losB = lineSensorB.getLOS(lineDateB,pleiadesViewingModelB.dimension/2);
GeodeticPoint centerPointB = ruggedB.directLocation(lineDateB, positionB, losB);
System.out.format(Locale.US, "center geodetic position : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(centerPointB.getLatitude()),
FastMath.toDegrees(centerPointB.getLongitude()),centerPointB.getAltitude());
lineDateB = lineSensorB.getDate(0);
losB = lineSensorB.getLOS(lineDateB,0);
GeodeticPoint firstPointB = ruggedB.directLocation(lineDateB, positionB, losB);
System.out.format(Locale.US, "first geodetic position : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(firstPointB.getLatitude()),
FastMath.toDegrees(firstPointB.getLongitude()),firstPointB.getAltitude());
lineDateB = lineSensorB.getDate(pleiadesViewingModelB.dimension-1);
losB = lineSensorB.getLOS(lineDateB,pleiadesViewingModelB.dimension-1);
GeodeticPoint lastPointB = ruggedB.directLocation(lineDateB, positionB, losB);
System.out.format(Locale.US, "last geodetic position : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(lastPointB.getLatitude()),
FastMath.toDegrees(lastPointB.getLongitude()),lastPointB.getAltitude());
double distance = DistanceTools.computeDistanceRad(centerPointA.getLongitude(), centerPointA.getLatitude(),
centerPointB.getLongitude(), centerPointB.getLatitude());
System.out.format("distance %f meters %n",distance);
double rollValueA = FastMath.toRadians(0.6);
double pitchValueA = FastMath.toRadians(0.0);
double factorValue = 1.00;
System.out.format("roll : %3.5f pitch : %3.5f factor : %3.5f \n",rollValueA,pitchValueA,factorValue);
ruggedA.
getLineSensor(SensorNameA).
getParametersDrivers().
filter(driver -> driver.getName().equals("roll")).
findFirst().get().setValue(rollValueA);
ruggedA.
getLineSensor(SensorNameA).
getParametersDrivers().
filter(driver -> driver.getName().equals("pitch")).
findFirst().get().setValue(pitchValueA);
ruggedA.
getLineSensor(SensorNameA).
getParametersDrivers().
filter(driver -> driver.getName().equals("factor")).
findFirst().get().setValue(factorValue);
lineDateA = lineSensorA.getDate(pleiadesViewingModelA.dimension/2);
losA = lineSensorA.getLOS(lineDateA,pleiadesViewingModelA.dimension/2);
centerPointA = ruggedA.directLocation(lineDateA, positionA, losA);
System.out.format(Locale.US, "center geodetic position : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(centerPointA.getLatitude()),
FastMath.toDegrees(centerPointA.getLongitude()),centerPointA.getAltitude());
distance = DistanceTools.computeDistanceRad(centerPointA.getLongitude(), centerPointA.getLatitude(),
centerPointB.getLongitude(), centerPointB.getLatitude());
System.out.format("distance %f meters %n",distance);
// Search the sensor pixel seeing point
int minLine = 0;
int maxLine = pleiadesViewingModelB.dimension-1;
SensorPixel sensorPixelB = ruggedB.inverseLocation(SensorNameB, centerPointA, minLine, maxLine);
// we need to test if the sensor pixel is found in the prescribed lines otherwise the sensor pixel is null
if (sensorPixelB != null){
System.out.format(Locale.US, "Sensor Pixel found : line = %5.3f, pixel = %5.3f %n", sensorPixelB.getLineNumber(), sensorPixelB.getPixelNumber());
} else {
System.out.println("Sensor Pixel is null: point cannot be seen between the prescribed line numbers\n");
}
//add mapping
SensorToSensorMeasureGenerator measure = new SensorToSensorMeasureGenerator(pleiadesViewingModelA,ruggedA,pleiadesViewingModelB,ruggedB);
int lineSampling = 1000;
int pixelSampling = 1000;
measure.CreateMeasure(lineSampling, pixelSampling);
System.out.format("nb TiePoints %d %n", measure.getMeasureCount());
} catch (OrekitException oe) {
System.err.println(oe.getLocalizedMessage());
System.exit(1);
} catch (RuggedException re) {
System.err.println(re.getLocalizedMessage());
System.exit(1);
}
}
}
src/tutorials/java/AffinagePleiades/DistanceTools.java deleted 100644 → 0
View file @ 950e0db2
package AffinagePleiades;
import org.hipparchus.util.FastMath;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
public class DistanceTools {
/**
* Earth radius in cms
*/
public static final double EARTH_RADIUS = 637100000d;
public DistanceTools() {
}
private static boolean checkDiff(double lon1, double lat1, double lon2, double lat2)
{
if(FastMath.abs(lon1 - lon2) < 1e-15 && FastMath.abs(lat1 - lat2) < 1e-15 )
return false;
return true;
}
/**
* Convert the given degrees value in meters
*
* @param distanceDeg
* @return the given degrees value in meters
*/
public static double computeDistanceDeg(double lonDeg, double latDeg, double distanceDeg) {
double lonDeg1 = lonDeg;
double latDeg1 = latDeg;
double lonDeg2 = lonDeg + distanceDeg;
double latDeg2 = latDeg;
return computeDistance(lonDeg1, latDeg1, lonDeg2, latDeg2);
}
/**
* Compute distance between point (lonDeg1, latDeg1) and point (lonDeg2,
* latDeg2) in meters
*
* @param lonDeg1
* @param latDeg1
* @param lonDeg2
* @param latDeg2
* @return distance between point (lonDeg1, latDeg1) and point (lonDeg2,
* latDeg2)
*/
public static double computeDistance(double lonDeg1, double latDeg1, double lonDeg2, double latDeg2) {
double returned;
if(!checkDiff(lonDeg1,latDeg1,lonDeg2,latDeg2) )
return 0.0;
//System.out.format("******** %n");
//System.out.format("%8.14f %8.14f %n",lonDeg1,latDeg1);
//System.out.format("%3.14e %n",lonDeg2-lonDeg1);
//System.out.format("%8.14f %8.14f %n",lonDeg2,latDeg2);
//System.out.format("%3.14e %n",latDeg2-latDeg1);
double xRad1 = FastMath.toRadians(lonDeg1);
double xRad2 = FastMath.toRadians(lonDeg2);
double yRad1 = FastMath.toRadians(latDeg1);
double yRad2 = FastMath.toRadians(latDeg2);
returned = computeDistanceRad(xRad1, yRad1, xRad2, yRad2);
return returned;
}
/**
* distance in meters between point (xRad1, yRad1) and point (xRad2, yRad2)
*
* @param xRad1
* @param xRad2
* @param yRad1
* @param yRad2
* @return
*/
public static double computeDistanceRad(double xRad1, double yRad1, double xRad2, double yRad2) {
// get vectors on unit sphere from angular coordinates
Vector3D p1 = new Vector3D(yRad1, xRad1); //
Vector3D p2 = new Vector3D(yRad2, xRad2);
return EARTH_RADIUS /100* Vector3D.angle(p1, p2);
}
}
src/tutorials/java/AffinagePleiades/LocalisationMetrics.java deleted 100644 → 0
View file @ 950e0db2
/* Copyright 2013-2016 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package AffinagePleiades;
import org.hipparchus.linear.ArrayRealVector;
import org.hipparchus.linear.RealVector;
import org.orekit.rugged.api.SensorToGroundMapping;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.time.AbsoluteDate;
import java.util.Collections;
import java.util.Collection;
import java.util.Map;
import java.util.Set;
import java.util.Iterator;
import java.util.Locale;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
/**
* class for measure generation
*
* @author Jonathan Guinet
*/
public class LocalisationMetrics {
/** mapping */
private Set<Map.Entry<SensorPixel, GeodeticPoint>> groundTruthMappings;
private Set<Map.Entry<SensorPixel, GeodeticPoint>> estimationMappings;
private Rugged rugged;
private LineSensor sensor;
private PleiadesViewingModel viewingModel;
private int measureCount;
private boolean computeInDeg;
/* max residual distance */
private double resMax;
/* mean residual distance */
private double resMean;
/**
* Simple constructor.
* <p>
*
* </p>
*/
public LocalisationMetrics(SensorToGroundMapping groundTruthMapping, Rugged rugged, boolean computeInDeg)
throws RuggedException {
groundTruthMappings = groundTruthMapping.getMappings();
this.rugged = rugged;
this.sensor = rugged.getLineSensor(groundTruthMapping.getSensorName());
this.computeInDeg = computeInDeg;
this.computeMetrics();
}
/**
* Get the maximum residual;
*
* @return max residual
*/
public double getMaxResidual() {
return resMax;
}
/**
* Get the mean residual;
*
* @return mean residual
*/
public double getMeanResidual() {
return resMean;
}
public void computeMetrics() throws RuggedException {
// final RealVector longDiffVector;
// final RealVector latDiffVector;
// final RealVector altDiffVector;
double count = 0;
resMax = 0;
int k = groundTruthMappings.size();
Iterator<Map.Entry<SensorPixel, GeodeticPoint>> gtIt = groundTruthMappings.iterator();
while (gtIt.hasNext()) {
Map.Entry<SensorPixel, GeodeticPoint> gtMapping = gtIt.next();
final SensorPixel gtSP = gtMapping.getKey();
final GeodeticPoint gtGP = gtMapping.getValue();
AbsoluteDate date = sensor.getDate(gtSP.getLineNumber());
GeodeticPoint esGP = rugged.directLocation(date, sensor.getPosition(),
sensor.getLOS(date, (int) gtSP.getPixelNumber()));
double distance = 0;
if (this.computeInDeg == true) {
double lonDiff = esGP.getLongitude() - gtGP.getLongitude();
double latDiff = esGP.getLatitude() - gtGP.getLatitude();
double altDiff = esGP.getAltitude() - gtGP.getAltitude();
distance = Math.sqrt(lonDiff * lonDiff + latDiff * latDiff);
} else {
distance = DistanceTools.computeDistanceRad(esGP.getLongitude(), esGP.getLatitude(),
gtGP.getLongitude(), gtGP.getLatitude());
}
count += distance;
if (distance > resMax) {
resMax = distance;
}
// distanceVector.append(distance);
}
// resMax = distanceVector.getMaxValue();
// System.out.format(Locale.US, "max: %3.6e %n
// ",distanceVector.getMaxValue() )
resMean = count / k;
// System.out.format("number of points %d %n", k);
}
}
src/tutorials/java/AffinagePleiades/MeasureGenerator.java deleted 100644 → 0
View file @ 950e0db2
/* Copyright 2013-2016 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package AffinagePleiades;
import org.orekit.rugged.api.SensorToGroundMapping;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.time.AbsoluteDate;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.random.UncorrelatedRandomVectorGenerator;
import org.hipparchus.random.GaussianRandomGenerator;
import org.hipparchus.random.Well19937a;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
/** class for measure generation
* @author Jonathan Guinet
*/
public class MeasureGenerator {
/** mapping */
private SensorToGroundMapping mapping;
private Rugged rugged;
private LineSensor sensor;
private PleiadesViewingModel viewingModel;
private int measureCount;
/** Simple constructor.
* <p>
*
* </p>
*/
public MeasureGenerator(PleiadesViewingModel viewingModel, Rugged rugged) throws RuggedException
{
// generate reference mapping
String sensorName = viewingModel.getSensorName();
mapping = new SensorToGroundMapping(sensorName);
this.rugged = rugged;
this.viewingModel = viewingModel;
sensor = rugged.getLineSensor(mapping.getSensorName());
measureCount = 0;
}
public SensorToGroundMapping getMapping() {
return mapping;
}
public int getMeasureCount() {
return measureCount;
}
public void CreateMeasure(final int lineSampling,final int pixelSampling) throws RuggedException
{
for (double line = 0; line < viewingModel.dimension; line += lineSampling) {
AbsoluteDate date = sensor.getDate(line);
for (int pixel = 0; pixel < sensor.getNbPixels(); pixel += pixelSampling) {
GeodeticPoint gp2 = rugged.directLocation(date, sensor.getPosition(),
sensor.getLOS(date, pixel));
mapping.addMapping(new SensorPixel(line, pixel), gp2);
measureCount++;
}
}
}
public void CreateNoisyMeasure(final int lineSampling,final int pixelSampling, double [] pixErr,double [] altErr) throws RuggedException
{
Vector3D latLongError = estimateLatLongError();
double latErrorMean = pixErr[0]*latLongError.getX(); // in line: -0.000002 deg
double lonErrorMean = pixErr[2]*latLongError.getY(); // in line: 0.000012 deg
double latErrorStd = pixErr[1]*latLongError.getX(); // in line: -0.000002 deg
double lonErrorStd = pixErr[3]*latLongError.getY(); // in line: 0.000012 deg
System.out.format("Corresponding error estimation on ground: {lat mean: %1.10f rad , lat stddev: %1.10f rad} %n",latErrorMean,latErrorStd);
System.out.format("Corresponding error estimation on ground: {lon mean: %1.10f rad, lat stddev: %1.10f rad } %n",lonErrorMean,lonErrorStd);
System.out.format("Corresponding error estimation on ground: {alt mean: %1.3f m, alt stddev: %1.3f m } %n",altErr[0], altErr[0]);
// Gaussian random generator
// Build a null mean random uncorrelated vector generator with standard deviation corresponding to the estimated error on ground
double mean[] = {latErrorMean,lonErrorMean,altErr[0]};
double std[] = {latErrorStd,lonErrorStd,altErr[1]};
System.out.format(" mean : %1.10f %1.10f %1.10f %n",mean[0],mean[1],mean[2]);
System.out.format(" std : %1.10f %1.10f %1.10f %n",std[0],std[1],std[2]);
GaussianRandomGenerator rng = new GaussianRandomGenerator(new Well19937a(0xefac03d9be4d24b9l));
UncorrelatedRandomVectorGenerator rvg = new UncorrelatedRandomVectorGenerator(mean, std, rng);
System.out.format("Add a gaussian noise to measures without biais (null mean) and standard deviation corresponding to the estimated error on ground.%n");
for (double line = 0; line < viewingModel.dimension; line += lineSampling) {
AbsoluteDate date = sensor.getDate(line);
for (int pixel = 0; pixel < sensor.getNbPixels(); pixel += pixelSampling) {
// Components of generated vector follow (independent) Gaussian distribution
Vector3D vecRandom = new Vector3D(rvg.nextVector());
GeodeticPoint gp2 = rugged.directLocation(date, sensor.getPosition(),
sensor.getLOS(date, pixel));
GeodeticPoint gpNoisy = new GeodeticPoint(gp2.getLatitude()+vecRandom.getX(),
gp2.getLongitude()+vecRandom.getY(),
gp2.getAltitude()+vecRandom.getZ());
//if(line == 0) {
// System.out.format("Init gp: (%f,%d): %s %n",line,pixel,gp2.toString());
// System.out.format("Random: (%f,%d): %s %n",line,pixel,vecRandom.toString());
// System.out.format("Final gp: (%f,%d): %s %n",line,pixel,gpNoisy.toString());
//}
mapping.addMapping(new SensorPixel(line, pixel), gpNoisy);
measureCount++;
}
}
}
private Vector3D estimateLatLongError() throws RuggedException {
System.out.format("Uncertainty in pixel (in line) for a real geometric refining: 1 pixel (assumption)%n");
final int pix =sensor.getNbPixels()/2;
final int line= (int) FastMath.floor(pix); // assumption : same number of line and pixels;
System.out.format(" pixel size estimated at position pix : %d line : %d %n", pix, line);
final AbsoluteDate date = sensor.getDate(line);
GeodeticPoint gp_pix0 = rugged.directLocation(date, sensor.getPosition(), sensor.getLOS(date, pix));
final AbsoluteDate date1 = sensor.getDate(line+1);
GeodeticPoint gp_pix1 = rugged.directLocation(date1, sensor.getPosition(), sensor.getLOS(date1, pix+1));
double latErr=FastMath.abs(gp_pix0.getLatitude()-gp_pix1.getLatitude());
double lonErr=FastMath.abs(gp_pix0.getLongitude()-gp_pix1.getLongitude());
double dist = FastMath.sqrt(lonErr*lonErr + latErr*latErr)/FastMath.sqrt(2);
final double distanceX = DistanceTools.computeDistanceRad(gp_pix0.getLongitude(), gp_pix0.getLatitude(),gp_pix1.getLongitude(), gp_pix0.getLatitude());
final double distanceY = DistanceTools.computeDistanceRad(gp_pix0.getLongitude(), gp_pix0.getLatitude(),gp_pix0.getLongitude(), gp_pix1.getLatitude());
System.out.format(" estimated distance %3.3f %3.3f %n",distanceX, distanceY);
//System.out.format(" lat : %1.10f %1.10f %n", latErr, lonErr);
return new Vector3D(latErr,lonErr,0.0);
}
}
src/tutorials/java/AffinagePleiades/OrbitModel.java deleted 100644 → 0
View file @ 950e0db2
/*
* Copyright 2013-2016 CS Systèmes d'Information Licensed to CS Systèmes
* d'Information (CS) under one or more contributor license agreements. See the
* NOTICE file distributed with this work for additional information regarding
* copyright ownership. CS licenses this file to You under the Apache License,
* Version 2.0 (the "License"); you may not use this file except in compliance
* with the License. You may obtain a copy of the License at
* http://www.apache.org/licenses/LICENSE-2.0 Unless required by applicable law
* or agreed to in writing, software distributed under the License is
* distributed on an "AS IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY
* KIND, either express or implied. See the License for the specific language
* governing permissions and limitations under the License.
*/
package AffinagePleiades;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.RotationOrder;
import org.hipparchus.geometry.euclidean.threed.RotationConvention;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.ode.nonstiff.DormandPrince853Integrator;
import org.hipparchus.util.FastMath;
import java.util.ArrayList;
import java.util.List;
import java.util.Arrays;
import org.orekit.attitudes.AttitudeProvider;
import org.orekit.attitudes.NadirPointing;
import org.orekit.attitudes.YawCompensation;
import org.orekit.attitudes.LofOffset;
import org.orekit.attitudes.TabulatedLofOffset;
import org.orekit.bodies.BodyShape;
import org.orekit.bodies.CelestialBodyFactory;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.errors.OrekitException;
import org.orekit.forces.gravity.HolmesFeatherstoneAttractionModel;
import org.orekit.forces.gravity.ThirdBodyAttraction;
import org.orekit.forces.gravity.potential.GravityFieldFactory;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.Transform;
import org.orekit.frames.LOFType;
import org.orekit.orbits.CircularOrbit;
import org.orekit.orbits.Orbit;
import org.orekit.orbits.OrbitType;
import org.orekit.orbits.PositionAngle;
import org.orekit.propagation.Propagator;
import org.orekit.propagation.SpacecraftState;
import org.orekit.propagation.analytical.KeplerianPropagator;
import org.orekit.propagation.numerical.NumericalPropagator;
import org.orekit.propagation.sampling.OrekitFixedStepHandler;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.utils.Constants;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
import org.orekit.utils.AngularDerivativesFilter;
/**
* Orbit Model class to generate PV and Q
*
* @author jguinet
*/
public class OrbitModel {
private double[] LOFTransformRollPoly;
private double[] LOFTransformPitchPoly;
private double[] LOFTransformYawPoly;
private AbsoluteDate refDate;
public static void addSatellitePV(TimeScale gps, Frame eme2000, Frame itrf,
ArrayList<TimeStampedPVCoordinates> satellitePVList,
String absDate, double px, double py,
double pz, double vx, double vy,
double vz)
throws OrekitException {
AbsoluteDate ephemerisDate = new AbsoluteDate(absDate, gps);
Vector3D position = new Vector3D(px, py, pz);
Vector3D velocity = new Vector3D(vx, vy, vz);
PVCoordinates pvITRF = new PVCoordinates(position, velocity);
Transform transform = itrf.getTransformTo(eme2000, ephemerisDate);
Vector3D pEME2000 = transform.transformPosition(pvITRF.getPosition());
Vector3D vEME2000 = transform.transformVector(pvITRF.getVelocity());
satellitePVList
.add(new TimeStampedPVCoordinates(ephemerisDate, pEME2000, vEME2000,
Vector3D.ZERO));
}
public void addSatelliteQ(TimeScale gps,
ArrayList<TimeStampedAngularCoordinates> satelliteQList,
String absDate, double q0, double q1, double q2,
double q3) {
AbsoluteDate attitudeDate = new AbsoluteDate(absDate, gps);
Rotation rotation = new Rotation(q0, q1, q2, q3, true);
TimeStampedAngularCoordinates pair = new TimeStampedAngularCoordinates(attitudeDate,
rotation,
Vector3D.ZERO,
Vector3D.ZERO);
satelliteQList.add(pair);
}
public BodyShape createEarth()
throws OrekitException {
return new OneAxisEllipsoid(Constants.WGS84_EARTH_EQUATORIAL_RADIUS,
Constants.WGS84_EARTH_FLATTENING,
FramesFactory.getITRF(
IERSConventions.IERS_2010, true));
}
public NormalizedSphericalHarmonicsProvider createGravityField()
throws OrekitException {
return GravityFieldFactory.getNormalizedProvider(12, 12);
}
public Orbit createOrbit(double mu, AbsoluteDate date)
throws OrekitException {
// the following orbital parameters have been computed using
// Orekit tutorial about phasing, using the following configuration:
//
// orbit.date = 2012-01-01T00:00:00.000
// phasing.orbits.number = 143
// phasing.days.number = 10
// sun.synchronous.reference.latitude = 0
// sun.synchronous.reference.ascending = false
// sun.synchronous.mean.solar.time = 10:30:00
// gravity.field.degree = 12
// gravity.field.order = 12
Frame eme2000 = FramesFactory.getEME2000();
// return new CircularOrbit(7173352.811913891,
return new CircularOrbit(694000.0 +
Constants.WGS84_EARTH_EQUATORIAL_RADIUS, // ajouter
// Rayon
-4.029194321683225E-4, 0.0013530362644647786,
FastMath.toRadians(98.2), // pleiades
// inclinaison 98.2
// FastMath.toRadians(-86.47),
FastMath.toRadians(-86.47 + 180),
FastMath.toRadians(135.9 + 0.3),
PositionAngle.TRUE, eme2000, date, mu);
}
public void setLOFTransform(double[] rollPoly, double[] pitchPoly,
double[] yawPoly, AbsoluteDate refDate) {
LOFTransformRollPoly = rollPoly.clone();
LOFTransformPitchPoly = pitchPoly.clone();
LOFTransformYawPoly = yawPoly.clone();
this.refDate = refDate;
}
private double getPoly(double[] poly, double shift) {
double val = 0.0;
for (int coef = 0; coef < poly.length; coef++) {
val = val + poly[coef] * Math.pow(shift, coef);
}
return val;
}
private Rotation getOffset(BodyShape earth, Orbit orbit, double shift)
throws OrekitException {
LOFType type = LOFType.VVLH;
double roll = getPoly(LOFTransformRollPoly, shift);
double pitch = getPoly(LOFTransformPitchPoly, shift);
double yaw = getPoly(LOFTransformYawPoly, shift);
/**System.out.format("at shift %f roll %f pitch %f yaw %f \n ", shift,
roll, pitch, yaw);*/
LofOffset law = new LofOffset(orbit.getFrame(), type, RotationOrder.XYZ,
roll, pitch, yaw);
Rotation offsetAtt = law
.getAttitude(orbit, orbit.getDate().shiftedBy(shift),
orbit.getFrame())
.getRotation();
NadirPointing nadirPointing = new NadirPointing(orbit.getFrame(),
earth);
Rotation nadirAtt = nadirPointing
.getAttitude(orbit, orbit.getDate().getDate().shiftedBy(shift),
orbit.getFrame())
.getRotation();
Rotation offsetProper = offsetAtt
.compose(nadirAtt.revert(), RotationConvention.VECTOR_OPERATOR);
return offsetProper;
}
public AttitudeProvider createAttitudeProvider(BodyShape earth, Orbit orbit)
throws OrekitException {
LOFType type = LOFType.VVLH;
ArrayList<TimeStampedAngularCoordinates> list = new ArrayList<TimeStampedAngularCoordinates>();
for (double shift = -10.0; shift < +10.0; shift += 1e-2) {
list.add(new TimeStampedAngularCoordinates(refDate.shiftedBy(shift),
getOffset(earth, orbit,
shift),
Vector3D.ZERO,
Vector3D.ZERO));
}
TabulatedLofOffset tabulated = new TabulatedLofOffset(orbit
.getFrame(), type, list, 2, AngularDerivativesFilter.USE_R);
return tabulated;
}
public Propagator createPropagator(BodyShape earth,
NormalizedSphericalHarmonicsProvider gravityField,
Orbit orbit)
throws OrekitException {
AttitudeProvider attidudeProvider = createAttitudeProvider(earth,
orbit);
SpacecraftState state = new SpacecraftState(orbit, attidudeProvider
.getAttitude(orbit, orbit.getDate(), orbit.getFrame()), 1180.0);
// numerical model for improving orbit
OrbitType type = OrbitType.CIRCULAR;
double[][] tolerances = NumericalPropagator.tolerances(0.1, orbit,
type);
DormandPrince853Integrator integrator = new DormandPrince853Integrator(1.0e-4 *
orbit
.getKeplerianPeriod(),
1.0e-1 * orbit
.getKeplerianPeriod(),
tolerances[0],
tolerances[1]);
integrator.setInitialStepSize(1.0e-2 * orbit.getKeplerianPeriod());
NumericalPropagator numericalPropagator = new NumericalPropagator(integrator);
numericalPropagator
.addForceModel(new HolmesFeatherstoneAttractionModel(earth
.getBodyFrame(), gravityField));
numericalPropagator
.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory
.getSun()));
numericalPropagator
.addForceModel(new ThirdBodyAttraction(CelestialBodyFactory
.getMoon()));
numericalPropagator.setOrbitType(type);
numericalPropagator.setInitialState(state);
numericalPropagator.setAttitudeProvider(attidudeProvider);
return numericalPropagator;
}
public List<TimeStampedPVCoordinates>
orbitToPV(Orbit orbit, BodyShape earth, AbsoluteDate minDate,
AbsoluteDate maxDate, double step)
throws OrekitException {
Propagator propagator = new KeplerianPropagator(orbit);
propagator.setAttitudeProvider(createAttitudeProvider(earth, orbit));
propagator.propagate(minDate);
final List<TimeStampedPVCoordinates> list = new ArrayList<TimeStampedPVCoordinates>();
propagator.setMasterMode(step, new OrekitFixedStepHandler() {
public void init(SpacecraftState s0, AbsoluteDate t) {
}
public void handleStep(SpacecraftState currentState,
boolean isLast) {
list.add(new TimeStampedPVCoordinates(currentState
.getDate(), currentState.getPVCoordinates().getPosition(),
currentState
.getPVCoordinates()
.getVelocity(),
Vector3D.ZERO));
}
});
propagator.propagate(maxDate);
return list;
}
public List<TimeStampedAngularCoordinates>
orbitToQ(Orbit orbit, BodyShape earth, AbsoluteDate minDate,
AbsoluteDate maxDate, double step)
throws OrekitException {
Propagator propagator = new KeplerianPropagator(orbit);
propagator.setAttitudeProvider(createAttitudeProvider(earth, orbit));
propagator.propagate(minDate);
final List<TimeStampedAngularCoordinates> list = new ArrayList<TimeStampedAngularCoordinates>();
propagator.setMasterMode(step, new OrekitFixedStepHandler() {
public void init(SpacecraftState s0, AbsoluteDate t) {
}
public void handleStep(SpacecraftState currentState,
boolean isLast) {
list.add(new TimeStampedAngularCoordinates(currentState
.getDate(), currentState.getAttitude().getRotation(),
Vector3D.ZERO,
Vector3D.ZERO));
}
});
propagator.propagate(maxDate);
return list;
}
}
src/tutorials/java/AffinagePleiades/SensorToSensorMeasureGenerator.java deleted 100644 → 0
View file @ 950e0db2
/* Copyright 2013-2016 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package AffinagePleiades;
import org.orekit.rugged.api.SensorToSensorMapping;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.time.AbsoluteDate;
import org.hipparchus.analysis.differentiation.DerivativeStructure;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.random.UncorrelatedRandomVectorGenerator;
import org.hipparchus.random.GaussianRandomGenerator;
import org.hipparchus.random.Well19937a;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.utils.DSGenerator;
import java.util.Locale;
/** class for measure generation
* @author Jonathan Guinet
*/
public class SensorToSensorMeasureGenerator {
/** mapping */
private SensorToSensorMapping mapping;
private Rugged ruggedA;
private Rugged ruggedB;
private LineSensor sensorA;
private LineSensor sensorB;
private PleiadesViewingModel viewingModelA;
private PleiadesViewingModel viewingModelB;
private int measureCount;
private String sensorNameA;
private String sensorNameB;
/** SensorToSensorMeasureGenerator
* @param viewingModelA
* @param ruggedA
* @param viewingModelB
* @param ruggedB
* @throws RuggedException
*/
public SensorToSensorMeasureGenerator(PleiadesViewingModel viewingModelA, Rugged ruggedA,PleiadesViewingModel viewingModelB, Rugged ruggedB) throws RuggedException
{
// generate reference mapping
sensorNameA = viewingModelA.getSensorName();
sensorNameB = viewingModelB.getSensorName();
mapping = new SensorToSensorMapping(sensorNameA, sensorNameB);
this.ruggedA = ruggedA;
this.ruggedB = ruggedB;
this.viewingModelA = viewingModelA;
this.viewingModelB = viewingModelB;
sensorA = ruggedA.getLineSensor(mapping.getSensorNameA());
sensorB = ruggedB.getLineSensor(mapping.getSensorNameB());
measureCount = 0;
}
public SensorToSensorMapping getMapping() {
return mapping;
}
public int getMeasureCount() {
return measureCount;
}
/** tie point creation from directLocation with RuugedA and inverse location with RuggedB
* @param lineSampling sampling along lines
* @param pixelSampling sampling along columns
* @throws RuggedException
*/
public void CreateMeasure(final int lineSampling,final int pixelSampling) throws RuggedException
{
// Search the sensor pixel seeing point
final int minLine = 0;
final int maxLine = viewingModelB.dimension - 1;
for (double line = 0; line < viewingModelA.dimension; line += lineSampling) {
AbsoluteDate dateA = sensorA.getDate(line);
for (int pixelA = 0; pixelA < sensorA
.getNbPixels(); pixelA += pixelSampling) {
GeodeticPoint gpA = ruggedA
.directLocation(dateA, sensorA.getPosition(),
sensorA.getLOS(dateA, pixelA));
SensorPixel sensorPixelB = ruggedB
.inverseLocation(sensorNameB, gpA, minLine, maxLine);
// we need to test if the sensor pixel is found in the
// prescribed lines otherwise the sensor pixel is null
if (sensorPixelB != null) {
System.out.format(Locale.US,
"Sensor Pixel found : line = %5.3f, pixel = %5.3f %n",
sensorPixelB.getLineNumber(),
sensorPixelB.getPixelNumber());
// TODO test if distance is 0.0 with crossing LOS
final AbsoluteDate dateB = sensorB.getDate(sensorPixelB.getLineNumber());
double pixelB = sensorPixelB.getPixelNumber();
//TODO work with double pixel values
double distance = ruggedA.distanceBetweenLOS(sensorA,dateA, pixelA, sensorB ,dateB,(int) pixelB);
System.out.format(Locale.US,"distance %f %n",distance);
mapping.addMapping(new SensorPixel(line, pixelA),
sensorPixelB, 0.0);
measureCount++;
} else {
System.out
.println("Sensor Pixel is null: point cannot be seen between the prescribed line numbers\n");
}
}
}
}
}
src/tutorials/java/fr/cs/examples/AtmosphericRefractionExamples.java 0 → 100644
View file @ acf902ea
package fr.cs.examples;
import java.io.File;
import java.net.URISyntaxException;
import java.text.DecimalFormat;
import java.text.DecimalFormatSymbols;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.stat.descriptive.DescriptiveStatistics;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.frames.Frame;
import org.orekit.frames.FramesFactory;
import org.orekit.frames.Transform;
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.LinearLineDatation;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.rugged.los.FixedRotation;
import org.orekit.rugged.los.LOSBuilder;
import org.orekit.rugged.los.TimeDependentLOS;
import org.orekit.rugged.refraction.AtmosphericRefraction;
import org.orekit.rugged.refraction.MultiLayerModel;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.IERSConventions;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
/** A simple example of how to use the atmospheric refraction correction.
*/
public class AtmosphericRefractionExamples {
public static void main(String[] args) throws URISyntaxException {
// 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, "orekit-data");
DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(orekitData));
// Sensor's definition
// ===================
// Line of sight
// -------------
// The raw viewing direction of pixel i with respect to the instrument is defined by the vector:
List<Vector3D> rawDirs = new ArrayList<Vector3D>();
for (int i = 0; i < 2000; i++) {
// 20° field of view, 2000 pixels
rawDirs.add(new Vector3D(0., i*FastMath.toRadians(20.)/2000., 1.));
}
// The instrument is oriented 10° off nadir around the X-axis, we need to rotate the viewing
// direction to obtain the line of sight in the satellite frame
LOSBuilder losBuilder = new LOSBuilder(rawDirs);
losBuilder.addTransform(new FixedRotation("10-degrees-rotation", Vector3D.PLUS_I, FastMath.toRadians(10.)));
TimeDependentLOS lineOfSight = losBuilder.build();
// Datation model
// --------------
// We use Orekit for handling time and dates, and Rugged for defining the datation model:
TimeScale gps = TimeScalesFactory.getGPS();
AbsoluteDate absDate = new AbsoluteDate("2009-12-11T16:59:30.0", gps);
LinearLineDatation lineDatation = new LinearLineDatation(absDate, 1., 20.);
// Line sensor
// -----------
// With the LOS and the datation now define, we can initialize a line sensor object in Rugged:
LineSensor lineSensor = new LineSensor("mySensor", lineDatation, Vector3D.ZERO, lineOfSight);
// Satellite position, velocity and attitude
// =========================================
// Reference frames
// ----------------
// In our application, we simply need to know the name of the frames we are working with. Positions and
// velocities are given in the ITRF terrestrial frame, while the quaternions are given in EME2000
// inertial frame.
Frame eme2000 = FramesFactory.getEME2000();
boolean simpleEOP = true; // we don't want to compute tiny tidal effects at millimeter level
Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, simpleEOP);
// Satellite attitude
// ------------------
ArrayList<TimeStampedAngularCoordinates> satelliteQList = new ArrayList<TimeStampedAngularCoordinates>();
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:58:42.592937", -0.340236, 0.333952, -0.844012, -0.245684);
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:06.592937", -0.354773, 0.329336, -0.837871, -0.252281);
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:30.592937", -0.369237, 0.324612, -0.831445, -0.258824);
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:54.592937", -0.3836, 0.319792, -0.824743, -0.265299);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:00:18.592937", -0.397834, 0.314883, -0.817777, -0.271695);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:00:42.592937", -0.411912, 0.309895, -0.810561, -0.278001);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:06.592937", -0.42581, 0.304838, -0.803111, -0.284206);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:30.592937", -0.439505, 0.299722, -0.795442, -0.290301);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:54.592937", -0.452976, 0.294556, -0.787571, -0.296279);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:02:18.592937", -0.466207, 0.28935, -0.779516, -0.302131);
// Positions and velocities
// ------------------------
ArrayList<TimeStampedPVCoordinates> satellitePVList = new ArrayList<TimeStampedPVCoordinates>();
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:58:42.592937", -726361.466, -5411878.485, 4637549.599, -2463.635, -4447.634, -5576.736);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:04.192937", -779538.267, -5506500.533, 4515934.894, -2459.848, -4312.676, -5683.906);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:25.792937", -832615.368, -5598184.195, 4392036.13, -2454.395, -4175.564, -5788.201);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:47.392937", -885556.748, -5686883.696, 4265915.971, -2447.273, -4036.368, -5889.568);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:00:08.992937", -938326.32, -5772554.875, 4137638.207, -2438.478, -3895.166, -5987.957);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:00:30.592937", -990887.942, -5855155.21, 4007267.717, -2428.011, -3752.034, -6083.317);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:00:52.192937", -1043205.448, -5934643.836, 3874870.441, -2415.868, -3607.05, -6175.600);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:13.792937", -1095242.669, -6010981.571, 3740513.34, -2402.051, -3460.291, -6264.760);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:35.392937", -1146963.457, -6084130.93, 3604264.372, -2386.561, -3311.835, -6350.751);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:56.992937", -1198331.706, -6154056.146, 3466192.446, -2369.401, -3161.764, -6433.531);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:02:18.592937", -1249311.381, -6220723.191, 3326367.397, -2350.574, -3010.159, -6513.056);
// Rugged initialization
// ---------------------
RuggedBuilder builder = new RuggedBuilder().
setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID).
setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF).
setTimeSpan(absDate, absDate.shiftedBy(60.0), 0.01, 5. / lineSensor.getRate(0.)).
setTrajectory(InertialFrameId.EME2000,
satellitePVList, 4, CartesianDerivativesFilter.USE_P,
satelliteQList, 4, AngularDerivativesFilter.USE_R).
addLineSensor(lineSensor);
// Build Rugged without atmospheric refraction model
Rugged ruggedWithout = builder.build();
// Defines atmospheric refraction model (with the default multi layers model)
AtmosphericRefraction atmosphericRefraction = new MultiLayerModel(ruggedWithout.getEllipsoid());
// One can use its own atmospheric model that must extend the abstract refraction.AtmosphericRefraction class.
// Set the optional grid steps, for the inverse location computation. Useless for direct location computation.
// This setting is optional as default values exist. But can be useful for instance to improve results by given smaller steps
// or to speed up the computation bigger steps.
int pixelStep = 100;
int lineStep = 100;
atmosphericRefraction.setGridSteps(pixelStep, lineStep);
// To compute the inverse location grid with atmospheric refraction, a default margin is available:
// atmosphericRefraction.getComputationParameters().getDefaultInverseLocMargin();
// but can be modified with:
// atmosphericRefraction.setInverseLocMargin(margin);
// Build Rugged with atmospheric refraction model
builder.setRefractionCorrection(atmosphericRefraction);
Rugged ruggedWith = builder.build();
// Direct location on a line WITHOUT and WITH atmospheric correction
// -----------------------------------------------------------------
double chosenLine = 200.;
// Without atmosphere
GeodeticPoint[] gpWithoutAtmosphericRefractionCorrection = ruggedWithout.directLocation("mySensor", chosenLine);
// With atmosphere
GeodeticPoint[] gpWithAtmosphericRefractionCorrection = ruggedWith.directLocation("mySensor", chosenLine);
double earthRadius = ruggedWithout.getEllipsoid().getEquatorialRadius();
DescriptiveStatistics statDistance = new DescriptiveStatistics();
for (int i = 0; i < gpWithAtmosphericRefractionCorrection.length; i++) {
double currentRadius = earthRadius + (gpWithAtmosphericRefractionCorrection[i].getAltitude()+ gpWithoutAtmosphericRefractionCorrection[i].getAltitude())/2.;
double distance = computeDistanceInMeter(currentRadius, gpWithAtmosphericRefractionCorrection[i], gpWithoutAtmosphericRefractionCorrection[i]);
statDistance.addValue(distance);
}
System.out.format("Distance must be > 0 and < 2m:" +
" Max = " + DFS.format(statDistance.getMax()) +
" Min = " + DFS.format(statDistance.getMin()) +
" Median = " + DFS.format(statDistance.getPercentile(50.)) +
" Mean = " + DFS.format(statDistance.getMean()) +
" Std deviation= " + DFS.format(statDistance.getStandardDeviation()) +
"\n");
// Inverse loc WITH atmospheric correction
// ==========================================================================
final double epsilonPixel = 1.e-3;
final double epsilonLine = 1.1e-2;
int minLine = (int) FastMath.floor(lineSensor.getLine(ruggedWithout.getMinDate()));
int maxLine = (int) FastMath.ceil(lineSensor.getLine(ruggedWithout.getMaxDate()));
for (int i = 0; i < gpWithAtmosphericRefractionCorrection.length; i++) {
// to check if we go back to the initial point when taking the geodetic point with atmospheric correction
GeodeticPoint gpWith = gpWithAtmosphericRefractionCorrection[i];
SensorPixel sensorPixelReverseWith = ruggedWith.inverseLocation("mySensor", gpWith, minLine, maxLine);
if (sensorPixelReverseWith != null) {
// Compare the computed pixel vs expected
if (FastMath.abs(i - sensorPixelReverseWith.getPixelNumber()) > epsilonPixel) {
System.out.format("Problem with pixel " + i + " . Delta pixel= " + DFS.format((i - sensorPixelReverseWith.getPixelNumber())) + "\n");
}
// Compare the computed line vs the expected
if (FastMath.abs(chosenLine - sensorPixelReverseWith.getLineNumber())> epsilonLine) {
System.out.format("Probem with line, for pixel " + i + ". Delta line= " + DFS.format((chosenLine - sensorPixelReverseWith.getLineNumber())) + "\n");
}
} else {
System.out.println("Inverse location failed for pixel " + i + " with atmospheric refraction correction for geodetic point computed with" );
}
} // end loop on pixel i
System.out.format("Inverse location gave same pixel and line as direct location input for" +
" epsilon pixel= " + DFS.format(epsilonPixel) +
" and epsilon line= " + DFS.format(epsilonLine) +
"\n");
}
private static void addSatellitePV(TimeScale gps, Frame eme2000, Frame itrf,
ArrayList<TimeStampedPVCoordinates> satellitePVList,
String absDate,
double px, double py, double pz, double vx, double vy, double vz) {
AbsoluteDate ephemerisDate = new AbsoluteDate(absDate, gps);
Vector3D position = new Vector3D(px, py, pz); // in ITRF, unit: m
Vector3D velocity = new Vector3D(vx, vy, vz); // in ITRF, unit: m/s
PVCoordinates pvITRF = new PVCoordinates(position, velocity);
Transform transform = itrf.getTransformTo(eme2000, ephemerisDate);
PVCoordinates pvEME2000 = transform.transformPVCoordinates(pvITRF);
satellitePVList.add(new TimeStampedPVCoordinates(ephemerisDate, pvEME2000.getPosition(), pvEME2000.getVelocity(), Vector3D.ZERO));
}
private static void addSatelliteQ(TimeScale gps, ArrayList<TimeStampedAngularCoordinates> satelliteQList, String absDate,
double q0, double q1, double q2, double q3) {
AbsoluteDate attitudeDate = new AbsoluteDate(absDate, gps);
Rotation rotation = new Rotation(q0, q1, q2, q3, true); // q0 is the scalar term
TimeStampedAngularCoordinates pair =
new TimeStampedAngularCoordinates(attitudeDate, rotation, Vector3D.ZERO, Vector3D.ZERO);
satelliteQList.add(pair);
}
private static double computeDistanceInMeter(double earthRadius, final GeodeticPoint gp1, final GeodeticPoint gp2) {
// get vectors on unit sphere from angular coordinates
final Vector3D p1 = new Vector3D(gp1.getLongitude(), gp1.getLatitude());
final Vector3D p2 = new Vector3D(gp2.getLongitude(), gp2.getLatitude());
return earthRadius * Vector3D.angle(p1, p2);
}
/** Definition of print format. */
private static final DecimalFormatSymbols SYMBOLS = new DecimalFormatSymbols(Locale.US);
/** Defined of number of digit after the point. */
private static final DecimalFormat DFS = new DecimalFormat("#.#####", SYMBOLS);
}
src/tutorials/java/fr/cs/examples/DirectLocation.java
View file @ acf902ea
/* Copyright 2013-2016 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
......@@ -16,16 +16,16 @@
*/
package fr.cs.examples;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import java.io.File;
import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.frames.Frame;
......@@ -40,8 +40,8 @@ import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.LinearLineDatation;
import org.orekit.rugged.los.LOSBuilder;
import org.orekit.rugged.los.FixedRotation;
import org.orekit.rugged.los.LOSBuilder;
import org.orekit.rugged.los.TimeDependentLOS;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
......@@ -62,12 +62,16 @@ public class DirectLocation {
// 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, "orekit-data");
DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(orekitData));
DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(orekitData));
// Sensor's definition
// ===================
// Line of sight
// -------------
// The raw viewing direction of pixel i with respect to the instrument is defined by the vector:
List<Vector3D> rawDirs = new ArrayList<Vector3D>();
for (int i = 0; i < 2000; i++) {
//20° field of view, 2000 pixels
// 20° field of view, 2000 pixels
rawDirs.add(new Vector3D(0d, i*FastMath.toRadians(20)/2000d, 1d));
}
......@@ -78,23 +82,33 @@ public class DirectLocation {
TimeDependentLOS lineOfSight = losBuilder.build();
// Datation model
// --------------
// We use Orekit for handling time and dates, and Rugged for defining the datation model:
TimeScale gps = TimeScalesFactory.getGPS();
AbsoluteDate absDate = new AbsoluteDate("2009-12-11T16:59:30.0", gps);
LinearLineDatation lineDatation = new LinearLineDatation(absDate, 1d, 20);
LinearLineDatation lineDatation = new LinearLineDatation(absDate, 1d, 20);
// With the LOS and the datation now defined , we can initialize a line sensor object in Rugged:
// Line sensor
// -----------
// With the LOS and the datation now defined, we can initialize a line sensor object in Rugged:
LineSensor lineSensor = new LineSensor("mySensor", lineDatation, Vector3D.ZERO, lineOfSight);
// Satellite position, velocity and attitude
// =========================================
// Reference frames
// ----------------
// In our application, we simply need to know the name of the frames we are working with. Positions and
// velocities are given in the ITRF terrestrial frame, while the quaternions are given in EME2000
// inertial frame.
// inertial frame.
Frame eme2000 = FramesFactory.getEME2000();
boolean simpleEOP = true; // we don't want to compute tiny tidal effects at millimeter level
Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, simpleEOP);
// Satellite attitude
// ------------------
ArrayList<TimeStampedAngularCoordinates> satelliteQList = new ArrayList<TimeStampedAngularCoordinates>();
ArrayList<TimeStampedPVCoordinates> satellitePVList = new ArrayList<TimeStampedPVCoordinates>();
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:58:42.592937", -0.340236d, 0.333952d, -0.844012d, -0.245684d);
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:06.592937", -0.354773d, 0.329336d, -0.837871d, -0.252281d);
......@@ -107,6 +121,10 @@ public class DirectLocation {
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:54.592937", -0.452976d, 0.294556d, -0.787571d, -0.296279d);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:02:18.592937", -0.466207d, 0.28935d, -0.779516d, -0.302131d);
// Positions and velocities
// ------------------------
ArrayList<TimeStampedPVCoordinates> satellitePVList = new ArrayList<TimeStampedPVCoordinates>();
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:58:42.592937", -726361.466d, -5411878.485d, 4637549.599d, -2463.635d, -4447.634d, -5576.736d);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:04.192937", -779538.267d, -5506500.533d, 4515934.894d, -2459.848d, -4312.676d, -5683.906d);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:25.792937", -832615.368d, -5598184.195d, 4392036.13d, -2454.395d, -4175.564d, -5788.201d);
......@@ -119,16 +137,20 @@ public class DirectLocation {
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:56.992937", -1198331.706d, -6154056.146d, 3466192.446d, -2369.401d, -3161.764d, -6433.531d);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:02:18.592937", -1249311.381d, -6220723.191d, 3326367.397d, -2350.574d, -3010.159d, -6513.056d);
// Rugged initialization
// ---------------------
Rugged rugged = new RuggedBuilder().
setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID).
setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID).
setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF).
setTimeSpan(absDate, absDate.shiftedBy(60.0), 0.01, 5 / lineSensor.getRate(0)).
setTimeSpan(absDate, absDate.shiftedBy(60.0), 0.01, 5 / lineSensor.getRate(0)).
setTrajectory(InertialFrameId.EME2000,
satellitePVList, 4, CartesianDerivativesFilter.USE_P,
satelliteQList, 4, AngularDerivativesFilter.USE_R).
addLineSensor(lineSensor).
build();
// Direct location of a single sensor pixel (first line, first pixel)
// ------------------------------------------------------------------
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.
AbsoluteDate firstLineDate = lineSensor.getDate(0);
Vector3D los = lineSensor.getLOS(firstLineDate, 0);
......@@ -151,21 +173,21 @@ public class DirectLocation {
private static void addSatellitePV(TimeScale gps, Frame eme2000, Frame itrf,
ArrayList<TimeStampedPVCoordinates> satellitePVList,
String absDate,
double px, double py, double pz, double vx, double vy, double vz)
throws OrekitException {
double px, double py, double pz, double vx, double vy, double vz) {
AbsoluteDate ephemerisDate = new AbsoluteDate(absDate, gps);
Vector3D position = new Vector3D(px, py, pz);
Vector3D velocity = new Vector3D(vx, vy, vz);
Vector3D position = new Vector3D(px, py, pz); // in ITRF, unit: m
Vector3D velocity = new Vector3D(vx, vy, vz); // in ITRF, unit: m/s
PVCoordinates pvITRF = new PVCoordinates(position, velocity);
Transform transform = itrf.getTransformTo(eme2000, ephemerisDate);
PVCoordinates pvEME2000 = transform.transformPVCoordinates(pvITRF);
PVCoordinates pvEME2000 = transform.transformPVCoordinates(pvITRF);
satellitePVList.add(new TimeStampedPVCoordinates(ephemerisDate, pvEME2000.getPosition(), pvEME2000.getVelocity(), Vector3D.ZERO));
}
private static void addSatelliteQ(TimeScale gps, ArrayList<TimeStampedAngularCoordinates> satelliteQList, String absDate,
double q0, double q1, double q2, double q3) {
AbsoluteDate attitudeDate = new AbsoluteDate(absDate, gps);
Rotation rotation = new Rotation(q0, q1, q2, q3, true);
Rotation rotation = new Rotation(q0, q1, q2, q3, true); // q0 is the scalar term
TimeStampedAngularCoordinates pair =
new TimeStampedAngularCoordinates(attitudeDate, rotation, Vector3D.ZERO, Vector3D.ZERO);
satelliteQList.add(pair);
......
src/tutorials/java/fr/cs/examples/DirectLocationWithDEM.java
View file @ acf902ea
/* Copyright 2013-2016 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
......@@ -25,7 +25,7 @@ import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataProvidersManager;
import org.orekit.data.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.frames.Frame;
......@@ -65,12 +65,16 @@ public class DirectLocationWithDEM {
// 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, "orekit-data");
DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(orekitData));
DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(orekitData));
// Sensor's definition
// ===================
// Line of sight
// -------------
// The raw viewing direction of pixel i with respect to the instrument is defined by the vector:
List<Vector3D> rawDirs = new ArrayList<Vector3D>();
for (int i = 0; i < 2000; i++) {
//20° field of view, 2000 pixels
// 20° field of view, 2000 pixels
rawDirs.add(new Vector3D(0d, i*FastMath.toRadians(20)/2000d, 1d));
}
......@@ -81,23 +85,33 @@ public class DirectLocationWithDEM {
TimeDependentLOS lineOfSight = losBuilder.build();
// Datation model
// --------------
// We use Orekit for handling time and dates, and Rugged for defining the datation model:
TimeScale gps = TimeScalesFactory.getGPS();
AbsoluteDate absDate = new AbsoluteDate("2009-12-11T16:59:30.0", gps);
LinearLineDatation lineDatation = new LinearLineDatation(absDate, 1d, 20);
LinearLineDatation lineDatation = new LinearLineDatation(absDate, 1d, 20);
// With the LOS and the datation now defined , we can initialize a line sensor object in Rugged:
// Line sensor
// -----------
// With the LOS and the datation now defined, we can initialize a line sensor object in Rugged:
LineSensor lineSensor = new LineSensor("mySensor", lineDatation, Vector3D.ZERO, lineOfSight);
// Satellite position, velocity and attitude
// =========================================
// Reference frames
// ----------------
// In our application, we simply need to know the name of the frames we are working with. Positions and
// velocities are given in the ITRF terrestrial frame, while the quaternions are given in EME2000
// inertial frame.
// inertial frame.
Frame eme2000 = FramesFactory.getEME2000();
boolean simpleEOP = true; // we don't want to compute tiny tidal effects at millimeter level
Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, simpleEOP);
// Satellite attitude
// ------------------
ArrayList<TimeStampedAngularCoordinates> satelliteQList = new ArrayList<TimeStampedAngularCoordinates>();
ArrayList<TimeStampedPVCoordinates> satellitePVList = new ArrayList<TimeStampedPVCoordinates>();
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:58:42.592937", -0.340236d, 0.333952d, -0.844012d, -0.245684d);
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:06.592937", -0.354773d, 0.329336d, -0.837871d, -0.252281d);
......@@ -109,7 +123,10 @@ public class DirectLocationWithDEM {
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:30.592937", -0.439505d, 0.299722d, -0.795442d, -0.290301d);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:54.592937", -0.452976d, 0.294556d, -0.787571d, -0.296279d);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:02:18.592937", -0.466207d, 0.28935d, -0.779516d, -0.302131d);
// Positions and velocities
// ------------------------
ArrayList<TimeStampedPVCoordinates> satellitePVList = new ArrayList<TimeStampedPVCoordinates>();
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:58:42.592937", -726361.466d, -5411878.485d, 4637549.599d, -2463.635d, -4447.634d, -5576.736d);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:04.192937", -779538.267d, -5506500.533d, 4515934.894d, -2459.848d, -4312.676d, -5683.906d);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:25.792937", -832615.368d, -5598184.195d, 4392036.13d, -2454.395d, -4175.564d, -5788.201d);
......@@ -122,23 +139,27 @@ public class DirectLocationWithDEM {
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:56.992937", -1198331.706d, -6154056.146d, 3466192.446d, -2369.401d, -3161.764d, -6433.531d);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:02:18.592937", -1249311.381d, -6220723.191d, 3326367.397d, -2350.574d, -3010.159d, -6513.056d);
// Rugged initialization
// ---------------------
Rugged rugged = new RuggedBuilder().
setAlgorithm(AlgorithmId.DUVENHAGE).
setDigitalElevationModel(new VolcanicConeElevationUpdater(new GeodeticPoint(FastMath.toRadians(37.58),
FastMath.toRadians(-96.95),
2463.0),
FastMath.toRadians(30.0), 16.0,
FastMath.toRadians(1.0), 1201),
8).
setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF).
setTimeSpan(absDate, absDate.shiftedBy(60.0), 0.01, 5 / lineSensor.getRate(0)).
setTrajectory(InertialFrameId.EME2000,
satellitePVList, 4, CartesianDerivativesFilter.USE_P,
satelliteQList, 4, AngularDerivativesFilter.USE_R).
addLineSensor(lineSensor).
build();
setAlgorithm(AlgorithmId.DUVENHAGE).
setDigitalElevationModel(new VolcanicConeElevationUpdater(new GeodeticPoint(FastMath.toRadians(37.58),
FastMath.toRadians(-96.95),
2463.0),
FastMath.toRadians(30.0), 16.0,
FastMath.toRadians(1.0), 1201),
8).
setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF).
setTimeSpan(absDate, absDate.shiftedBy(60.0), 0.01, 5 / lineSensor.getRate(0)).
setTrajectory(InertialFrameId.EME2000,
satellitePVList, 4, CartesianDerivativesFilter.USE_P,
satelliteQList, 4, AngularDerivativesFilter.USE_R).
addLineSensor(lineSensor).
build();
// Direct location of a single sensor pixel (first line, first pixel)
// ------------------------------------------------------------------
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.
AbsoluteDate firstLineDate = lineSensor.getDate(0);
Vector3D los = lineSensor.getLOS(firstLineDate, 0);
......@@ -159,16 +180,16 @@ public class DirectLocationWithDEM {
}
private static void addSatellitePV(TimeScale gps, Frame eme2000, Frame itrf,
ArrayList<TimeStampedPVCoordinates> satellitePVList,
String absDate,
double px, double py, double pz, double vx, double vy, double vz)
throws OrekitException {
ArrayList<TimeStampedPVCoordinates> satellitePVList,
String absDate,
double px, double py, double pz, double vx, double vy, double vz) {
AbsoluteDate ephemerisDate = new AbsoluteDate(absDate, gps);
Vector3D position = new Vector3D(px, py, pz);
Vector3D velocity = new Vector3D(vx, vy, vz);
Vector3D position = new Vector3D(px, py, pz); // in ITRF, unit: m
Vector3D velocity = new Vector3D(vx, vy, vz); // in ITRF, unit: m/s
PVCoordinates pvITRF = new PVCoordinates(position, velocity);
Transform transform = itrf.getTransformTo(eme2000, ephemerisDate);
PVCoordinates pvEME2000 = transform.transformPVCoordinates(pvITRF);
Transform transform = itrf.getTransformTo(eme2000, ephemerisDate);
PVCoordinates pvEME2000 = transform.transformPVCoordinates(pvITRF);
satellitePVList.add(new TimeStampedPVCoordinates(ephemerisDate, pvEME2000.getPosition(), pvEME2000.getVelocity(), Vector3D.ZERO));
}
......@@ -176,9 +197,9 @@ public class DirectLocationWithDEM {
private static void addSatelliteQ(TimeScale gps, ArrayList<TimeStampedAngularCoordinates> satelliteQList, String absDate,
double q0, double q1, double q2, double q3) {
AbsoluteDate attitudeDate = new AbsoluteDate(absDate, gps);
Rotation rotation = new Rotation(q0, q1, q2, q3, true);
Rotation rotation = new Rotation(q0, q1, q2, q3, true); // q0 is the scalar term
TimeStampedAngularCoordinates pair =
new TimeStampedAngularCoordinates(attitudeDate, rotation, Vector3D.ZERO, Vector3D.ZERO);
new TimeStampedAngularCoordinates(attitudeDate, rotation, Vector3D.ZERO, Vector3D.ZERO);
satelliteQList.add(pair);
}
......@@ -199,8 +220,9 @@ public class DirectLocationWithDEM {
this.n = n;
}
public void updateTile(double latitude, double longitude, UpdatableTile tile)
throws RuggedException {
@Override
public void updateTile(double latitude, double longitude, UpdatableTile tile) {
double step = size / (n - 1);
double minLatitude = size * FastMath.floor(latitude / size);
double minLongitude = size * FastMath.floor(longitude / size);
......@@ -211,8 +233,8 @@ public class DirectLocationWithDEM {
for (int j = 0; j < n; ++j) {
double cellLongitude = minLongitude + j * step;
double distance = Constants.WGS84_EARTH_EQUATORIAL_RADIUS *
FastMath.hypot(cellLatitude - summit.getLatitude(),
cellLongitude - summit.getLongitude());
FastMath.hypot(cellLatitude - summit.getLatitude(),
cellLongitude - summit.getLongitude());
double altitude = FastMath.max(summit.getAltitude() - distance * sinSlope,
base);
tile.setElevation(i, j, altitude);
......
src/tutorials/java/fr/cs/examples/InverseLocation.java
View file @ acf902ea
/* Copyright 2013-2016 CS Systèmes d'Information
* Licensed to CS Systèmes d'Information (CS) under one or more
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
......@@ -25,7 +25,8 @@ import org.hipparchus.geometry.euclidean.threed.Rotation;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.data.DataProvidersManager;
import org.orekit.bodies.OneAxisEllipsoid;
import org.orekit.data.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.frames.Frame;
......@@ -44,6 +45,7 @@ import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.rugged.los.FixedRotation;
import org.orekit.rugged.los.LOSBuilder;
import org.orekit.rugged.los.TimeDependentLOS;
import org.orekit.rugged.utils.RoughVisibilityEstimator;
import org.orekit.time.AbsoluteDate;
import org.orekit.time.TimeScale;
import org.orekit.time.TimeScalesFactory;
......@@ -63,12 +65,16 @@ public class InverseLocation {
// 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, "orekit-data");
DataProvidersManager.getInstance().addProvider(new DirectoryCrawler(orekitData));
DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(orekitData));
// Sensor's definition
// ===================
// Line of sight
// -------------
// The raw viewing direction of pixel i with respect to the instrument is defined by the vector:
List<Vector3D> rawDirs = new ArrayList<Vector3D>();
for (int i = 0; i < 2000; i++) {
//20° field of view, 2000 pixels
// 20° field of view, 2000 pixels
rawDirs.add(new Vector3D(0d, i*FastMath.toRadians(20)/2000d, 1d));
}
......@@ -79,26 +85,34 @@ public class InverseLocation {
TimeDependentLOS lineOfSight = losBuilder.build();
// Datation model
// --------------
// We use Orekit for handling time and dates, and Rugged for defining the datation model:
TimeScale gps = TimeScalesFactory.getGPS();
AbsoluteDate absDate = new AbsoluteDate("2009-12-11T16:59:30.0", gps);
LinearLineDatation lineDatation = new LinearLineDatation(absDate, 1d, 20);
LinearLineDatation lineDatation = new LinearLineDatation(absDate, 1d, 20);
// Line sensor
// -----------
// With the LOS and the datation now defined , we can initialize a line sensor object in Rugged:
String sensorName = "mySensor";
LineSensor lineSensor = new LineSensor(sensorName, lineDatation, Vector3D.ZERO, lineOfSight);
// Satellite position, velocity and attitude
// =========================================
// Reference frames
// ----------------
// In our application, we simply need to know the name of the frames we are working with. Positions and
// velocities are given in the ITRF terrestrial frame, while the quaternions are given in EME2000
// inertial frame.
// inertial frame.
Frame eme2000 = FramesFactory.getEME2000();
boolean simpleEOP = true; // we don't want to compute tiny tidal effects at millimeter level
Frame itrf = FramesFactory.getITRF(IERSConventions.IERS_2010, simpleEOP);
// Satellite attitude
// ------------------
ArrayList<TimeStampedAngularCoordinates> satelliteQList = new ArrayList<TimeStampedAngularCoordinates>();
ArrayList<TimeStampedPVCoordinates> satellitePVList = new ArrayList<TimeStampedPVCoordinates>();
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:58:42.592937", -0.340236d, 0.333952d, -0.844012d, -0.245684d);
addSatelliteQ(gps, satelliteQList, "2009-12-11T16:59:06.592937", -0.354773d, 0.329336d, -0.837871d, -0.252281d);
......@@ -111,6 +125,10 @@ public class InverseLocation {
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:01:54.592937", -0.452976d, 0.294556d, -0.787571d, -0.296279d);
addSatelliteQ(gps, satelliteQList, "2009-12-11T17:02:18.592937", -0.466207d, 0.28935d, -0.779516d, -0.302131d);
// Positions and velocities
// ------------------------
ArrayList<TimeStampedPVCoordinates> satellitePVList = new ArrayList<TimeStampedPVCoordinates>();
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:58:42.592937", -726361.466d, -5411878.485d, 4637549.599d, -2463.635d, -4447.634d, -5576.736d);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:04.192937", -779538.267d, -5506500.533d, 4515934.894d, -2459.848d, -4312.676d, -5683.906d);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T16:59:25.792937", -832615.368d, -5598184.195d, 4392036.13d, -2454.395d, -4175.564d, -5788.201d);
......@@ -123,37 +141,77 @@ public class InverseLocation {
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:01:56.992937", -1198331.706d, -6154056.146d, 3466192.446d, -2369.401d, -3161.764d, -6433.531d);
addSatellitePV(gps, eme2000, itrf, satellitePVList, "2009-12-11T17:02:18.592937", -1249311.381d, -6220723.191d, 3326367.397d, -2350.574d, -3010.159d, -6513.056d);
Rugged rugged = new RuggedBuilder().
setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID).
setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF).
setTimeSpan(absDate, absDate.shiftedBy(60.0), 0.01, 5 / lineSensor.getRate(0)).
setTrajectory(InertialFrameId.EME2000,
satellitePVList, 4, CartesianDerivativesFilter.USE_P,
satelliteQList, 4, AngularDerivativesFilter.USE_R).
addLineSensor(lineSensor).
build();
// Rugged initialization
// ---------------------
RuggedBuilder ruggedBuilder = new RuggedBuilder();
ruggedBuilder.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID);
ruggedBuilder.setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF);
ruggedBuilder.setTimeSpan(absDate, absDate.shiftedBy(60.0), 0.01, 5 / lineSensor.getRate(0));
ruggedBuilder.setTrajectory(InertialFrameId.EME2000,
satellitePVList, 4, CartesianDerivativesFilter.USE_P,
satelliteQList, 4, AngularDerivativesFilter.USE_R).
addLineSensor(lineSensor);
Rugged rugged = ruggedBuilder.build();
// Inverse location of a Geodetic Point
// ------------------------------------
// Point defined by its latitude, longitude and altitude
double latitude = FastMath.toRadians(37.585);
double longitude = FastMath.toRadians(-96.949);
double altitude = 0.0d;
// For a GeodeticPoint : angles are given in radians and altitude in meters
GeodeticPoint gp = new GeodeticPoint(latitude, longitude, altitude);
// Search the sensor pixel seeing point
// ....................................
// interval of lines where to search the point
int minLine = 0;
int maxLine = 100;
SensorPixel sensorPixel = rugged.inverseLocation(sensorName, gp, minLine, maxLine);
// we need to test if the sensor pixel is found in the prescribed lines otherwise the sensor pixel is null
// Be careful: no exception is thrown if the pixel is not found
if (sensorPixel != null){
// The pixel was found
System.out.format(Locale.US, "Sensor Pixel found : line = %5.3f, pixel = %5.3f %n", sensorPixel.getLineNumber(), sensorPixel.getPixelNumber());
} else {
// The pixel was not found: set to null
System.out.println("Sensor Pixel is null: point cannot be seen between the prescribed line numbers\n");
}
// Find the date at which the sensor sees the ground point
// .......................................................
AbsoluteDate dateLine = rugged.dateLocation(sensorName, gp, minLine, maxLine);
System.out.println("Date at which the sensor sees the ground point : " + dateLine);
// How to find min and max lines ? with the RoughVisibilityEstimator
// -------------------------------
// Create a RoughVisibilityEstimator for inverse location
OneAxisEllipsoid oneAxisEllipsoid = ruggedBuilder.getEllipsoid();
Frame pvFrame = ruggedBuilder.getInertialFrame();
// Initialize the RoughVisibilityEstimator
RoughVisibilityEstimator roughVisibilityEstimator = new RoughVisibilityEstimator(oneAxisEllipsoid, pvFrame, satellitePVList);
// Compute the approximated line with a rough estimator
AbsoluteDate roughLineDate = roughVisibilityEstimator.estimateVisibility(gp);
double roughLine = lineSensor.getLine(roughLineDate);
// Compute the min / max lines interval using a margin around the roughLine
int sensorMinLine= 0;
int sensorMaxLine = 1000;
int margin = 100;
int minLineRough = (int) FastMath.max(roughLine - margin, sensorMinLine);
int maxLineRough = (int) FastMath.min(roughLine + margin, sensorMaxLine);
SensorPixel sensorPixelRoughLine = rugged.inverseLocation(sensorName, gp, minLineRough, maxLineRough);
System.out.format(Locale.US, "Rough line found = %5.1f; InverseLocation gives (margin of %d around rough line): line = %5.3f, pixel = %5.3f %n", roughLine, margin, sensorPixelRoughLine.getLineNumber(), sensorPixel.getPixelNumber());
} catch (OrekitException oe) {
System.err.println(oe.getLocalizedMessage());
System.exit(1);
......@@ -165,25 +223,25 @@ public class InverseLocation {
}
private static void addSatellitePV(TimeScale gps, Frame eme2000, Frame itrf,
ArrayList<TimeStampedPVCoordinates> satellitePVList,
String absDate,
double px, double py, double pz, double vx, double vy, double vz)
throws OrekitException {
ArrayList<TimeStampedPVCoordinates> satellitePVList,
String absDate,
double px, double py, double pz, double vx, double vy, double vz) {
AbsoluteDate ephemerisDate = new AbsoluteDate(absDate, gps);
Vector3D position = new Vector3D(px, py, pz);
Vector3D velocity = new Vector3D(vx, vy, vz);
Vector3D position = new Vector3D(px, py, pz); // in ITRF, unit: m
Vector3D velocity = new Vector3D(vx, vy, vz); // in ITRF, unit: m/s
PVCoordinates pvITRF = new PVCoordinates(position, velocity);
Transform transform = itrf.getTransformTo(eme2000, ephemerisDate);
PVCoordinates pvEME2000 = transform.transformPVCoordinates(pvITRF);
PVCoordinates pvEME2000 = transform.transformPVCoordinates(pvITRF);
satellitePVList.add(new TimeStampedPVCoordinates(ephemerisDate, pvEME2000.getPosition(), pvEME2000.getVelocity(), Vector3D.ZERO));
}
private static void addSatelliteQ(TimeScale gps, ArrayList<TimeStampedAngularCoordinates> satelliteQList, String absDate,
double q0, double q1, double q2, double q3) {
AbsoluteDate attitudeDate = new AbsoluteDate(absDate, gps);
Rotation rotation = new Rotation(q0, q1, q2, q3, true);
Rotation rotation = new Rotation(q0, q1, q2, q3, true); // q0 is the scalar term
TimeStampedAngularCoordinates pair =
new TimeStampedAngularCoordinates(attitudeDate, rotation, Vector3D.ZERO, Vector3D.ZERO);
new TimeStampedAngularCoordinates(attitudeDate, rotation, Vector3D.ZERO, Vector3D.ZERO);
satelliteQList.add(pair);
}
......
src/tutorials/java/fr/cs/examples/refiningPleiades/GroundRefining.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package fr.cs.examples.refiningPleiades;
import java.io.File;
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.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.forces.gravity.potential.NormalizedSphericalHarmonicsProvider;
import org.orekit.orbits.Orbit;
import org.orekit.rugged.adjustment.measurements.SensorToGroundMapping;
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
import fr.cs.examples.refiningPleiades.generators.GroundMeasurementGenerator;
import fr.cs.examples.refiningPleiades.generators.Noise;
import fr.cs.examples.refiningPleiades.metrics.DistanceTools;
import fr.cs.examples.refiningPleiades.models.OrbitModel;
import fr.cs.examples.refiningPleiades.models.PleiadesViewingModel;
/**
* Class for testing refining (Ground Control Points GCP study)
* with or without noisy measurements
* @author Jonathan Guinet
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @see SensorToGroundMapping
* @see GroundMeasurementGenerator
* @since 2.0
*/
public class GroundRefining extends Refining {
/** Pleiades viewing model */
private static PleiadesViewingModel pleiadesViewingModel;
/** Sensor name */
private static String sensorName;
/** Main function
*/
public static void main(String[] args) {
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, "orekit-data");
DataContext.getDefault().getDataProvidersManager().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");
AbsoluteDate minDate = pleiadesViewingModel.getMinDate();
AbsoluteDate maxDate = pleiadesViewingModel.getMaxDate();
AbsoluteDate refDate = pleiadesViewingModel.getDatationReference();
LineSensor lineSensor = pleiadesViewingModel.getLineSensor();
// Satellite position, velocity and attitude: create an orbit model
// ----------------------------------------------------------------
OrbitModel orbitmodel = new OrbitModel();
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);
ruggedBuilder.setTrajectory(InertialFrameId.EME2000, satellitePVList,nbPVPoints,
CartesianDerivativesFilter.USE_PV, satelliteQList,
nbQPoints, AngularDerivativesFilter.USE_R);
ruggedBuilder.setName("Rugged_refining");
Rugged rugged = ruggedBuilder.build();
// Compute ground sample distance (GSD)
// ------------------------------------
double [] gsd = refining.computeGSD(rugged, lineSensor);
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");
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(rugged, sensorName, rollValue, pitchValue, factorValue);
// Generate measurements (observations) from physical model disrupted
// ------------------------------------------------------------------
int lineSampling = 1000;
int pixelSampling = 1000;
// Noise definition
// distribution: gaussian(0), vector dimension:3
final Noise noise = new Noise(0,3);
// lat mean, long mean, alt mean
final double[] mean = {0.0,0.0,0.0};
// lat std, long std, alt std}
final double[] standardDeviation = {0.0,0.0,0.0};
noise.setMean(mean);
noise.setStandardDeviation(standardDeviation);
// Ground measurements
GroundMeasurementGenerator measurements = refining.generateNoisyPoints(lineSampling, pixelSampling,
rugged, sensorName,
pleiadesViewingModel.getDimension(),
noise);
// Compute ground truth residuals
// ------------------------------
System.out.format("\n**** Ground truth residuals **** %n");
refining.computeMetrics(measurements.getGroundMapping(), rugged, false);
// Initialize physical model without disruptions
// ---------------------------------------------
System.out.format("\n**** Initialize physical model without disruptions: reset Roll/Pitch/Factor **** %n");
refining.resetModel(rugged, sensorName, true);
// Compute initial residuals
// -------------------------
System.out.format("\n**** Initial Residuals **** %n");
refining.computeMetrics(measurements.getGroundMapping(), rugged, false);
// Start optimization
// ------------------
System.out.format("\n**** Start optimization **** %n");
int maxIterations = 100;
double convergenceThreshold = 1.e-11;
refining.optimization(maxIterations, convergenceThreshold, measurements.getObservables(), rugged);
// Check estimated values
// ----------------------
System.out.format("\n**** Check parameters ajustement **** %n");
refining.paramsEstimation(rugged, sensorName, rollValue, pitchValue, factorValue);
// Compute statistics
// ------------------
System.out.format("\n**** Compute Statistics **** %n");
// Residuals computed in meters
refining.computeMetrics(measurements.getGroundMapping(), rugged, false);
// Residuals computed in degrees
refining.computeMetrics(measurements.getGroundMapping(), rugged, true);
} catch (OrekitException oe) {
System.err.println(oe.getLocalizedMessage());
System.exit(1);
} catch (RuggedException re) {
System.err.println(re.getLocalizedMessage());
System.exit(1);
}
}
/** Constructor */
private GroundRefining() {
sensorName = "line";
pleiadesViewingModel = new PleiadesViewingModel(sensorName);
}
/** Estimate ground sample distance (GSD)
* @param LineSensor line sensor
* @return the GSD
*/
private double[] computeGSD(final Rugged rugged, final LineSensor lineSensor) {
// Get number of line
int dimension = pleiadesViewingModel.getDimension();
// 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,dimension-1);
// Get center geodetic point
AbsoluteDate lineDate = lineSensor.getDate(dimension/2);
los = lineSensor.getLOS(lineDate,dimension/2);
// Get upper right geodetic point
int pixelPosition = dimension-1;
los = lineSensor.getLOS(firstLineDate,pixelPosition);
GeodeticPoint upperRight = rugged.directLocation(firstLineDate, position, los);
// Get lower left geodetic point
AbsoluteDate lineDate_y = lineSensor.getDate(dimension-1);
los = lineSensor.getLOS(lineDate_y,0);
GeodeticPoint lowerLeft = rugged.directLocation(lineDate_y, position, los);
double gsdX = DistanceTools.computeDistanceInMeter(upLeftPoint, upperRight)/dimension;
double gsdY = DistanceTools.computeDistanceInMeter(upLeftPoint, lowerLeft)/dimension;
double [] gsd = {gsdX, gsdY};
return gsd;
}
// /**
// * Get the Pleiades viewing model
// * @return the Pleiades viewing model
// */
// public PleiadesViewingModel getPleiadesViewingModel() {
// return pleiadesViewingModel;
// }
//
// /**
// * Set the Pleiades viewing model
// * @param pleiadesViewingModel Pleiades viewing model to set
// */
// public void setPleiadesViewingModel(final PleiadesViewingModel pleiadesViewingModel) {
// this.pleiadesViewingModel = pleiadesViewingModel;
// }
//
// /**
// * Get the orbit model
// * @return the orbit model
// */
// public OrbitModel getOrbitmodel() {
// return orbitmodel;
// }
//
// /**
// * Set the orbit model
// * @param orbitmodel the orbit model to set
// */
// public void setOrbitmodel(final OrbitModel orbitmodel) {
// this.orbitmodel = orbitmodel;
// }
//
// /**
// * Get the sensor name
// * @return the sensor name
// */
// public String getSensorName() {
// return sensorName;
// }
//
// /**
// * Get the Rugged instance
// * @return the rugged instance
// */
// public Rugged getRugged() {
// return rugged;
// }
//
// /**
// * Set the Rugged instance
// * @param rugged the Rugged instance to set
// */
// public void setRugged(final Rugged rugged) {
// this.rugged = rugged;
// }
}
src/tutorials/java/fr/cs/examples/refiningPleiades/InterRefining.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package fr.cs.examples.refiningPleiades;
import java.io.File;
import java.util.Arrays;
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.DataContext;
import org.orekit.data.DirectoryCrawler;
import org.orekit.errors.OrekitException;
import org.orekit.orbits.Orbit;
import org.orekit.rugged.adjustment.measurements.SensorToGroundMapping;
import org.orekit.rugged.api.AlgorithmId;
import org.orekit.rugged.api.BodyRotatingFrameId;
import org.orekit.rugged.api.EllipsoidId;
import org.orekit.rugged.api.InertialFrameId;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.api.RuggedBuilder;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.time.AbsoluteDate;
import org.orekit.utils.AngularDerivativesFilter;
import org.orekit.utils.CartesianDerivativesFilter;
import org.orekit.utils.Constants;
import org.orekit.utils.PVCoordinates;
import org.orekit.utils.TimeStampedAngularCoordinates;
import org.orekit.utils.TimeStampedPVCoordinates;
import fr.cs.examples.refiningPleiades.generators.GroundMeasurementGenerator;
import fr.cs.examples.refiningPleiades.generators.InterMeasurementGenerator;
import fr.cs.examples.refiningPleiades.generators.Noise;
import fr.cs.examples.refiningPleiades.metrics.DistanceTools;
import fr.cs.examples.refiningPleiades.models.OrbitModel;
import fr.cs.examples.refiningPleiades.models.PleiadesViewingModel;
/**
* Class for testing refining (tie points study)
* with or without noisy measurements
* @author Jonathan Guinet
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @see SensorToGroundMapping
* @see GroundMeasurementGenerator
* @since 2.0
*/
public class InterRefining extends Refining {
/** Pleiades viewing model A */
private static PleiadesViewingModel pleiadesViewingModelA;
/** Pleiades viewing model B */
private static PleiadesViewingModel pleiadesViewingModelB;
/** Sensor name A corresponding to viewing model A */
private static String sensorNameA;
/** Sensor name A corresponding to viewing model B */
private static String sensorNameB;
/** Main function
*/
public static void main(String[] args) {
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, "orekit-data");
DataContext.getDefault().getDataProvidersManager().addProvider(new DirectoryCrawler(orekitData));
// Initialize refining context
// ---------------------------
InterRefining refining = new InterRefining();
// Sensors's definition: creation of 2 Pleiades viewing models
// -----------------------------------------------------------
System.out.format("**** Build Pleiades viewing model A and its orbit definition **** %n");
// 1/- Create First Pleiades Viewing Model
final AbsoluteDate minDateA = pleiadesViewingModelA.getMinDate();
final AbsoluteDate maxDateA = pleiadesViewingModelA.getMaxDate();
final AbsoluteDate refDateA = pleiadesViewingModelA.getDatationReference();
LineSensor lineSensorA = pleiadesViewingModelA.getLineSensor();
System.out.format(Locale.US, "Viewing model A - date min: %s max: %s ref: %s %n", minDateA, maxDateA, refDateA);
// ----Satellite position, velocity and attitude: create orbit model A
OrbitModel orbitmodelA = new OrbitModel();
BodyShape earthA = orbitmodelA.createEarth();
Orbit orbitA = orbitmodelA.createOrbit(Constants.EIGEN5C_EARTH_MU, refDateA);
// ----If no LOF Transform Attitude Provider is Nadir Pointing Yaw Compensation
final double [] rollPoly = {0.0,0.0,0.0};
double[] pitchPoly = {0.025,0.0};
double[] yawPoly = {0.0,0.0,0.0};
orbitmodelA.setLOFTransform(rollPoly, pitchPoly, yawPoly, minDateA);
// ----Satellite attitude
List<TimeStampedAngularCoordinates> satelliteQListA = orbitmodelA.orbitToQ(orbitA, earthA, minDateA.shiftedBy(-0.0), maxDateA.shiftedBy(+0.0), 0.25);
int nbQPoints = 2;
// ----Position and velocities
PVCoordinates PVA = orbitA.getPVCoordinates(earthA.getBodyFrame());
List<TimeStampedPVCoordinates> satellitePVListA = orbitmodelA.orbitToPV(orbitA, earthA, minDateA.shiftedBy(-0.0), maxDateA.shiftedBy(+0.0), 0.25);
int nbPVPoints = 8;
// ----Convert frame and coordinates type in one call
GeodeticPoint gpA = earthA.transform(PVA.getPosition(), earthA.getBodyFrame(), orbitA.getDate());
System.out.format(Locale.US, "(A) Geodetic Point at date %s : φ = %8.10f °, λ = %8.10f %n",
orbitA.getDate().toString(),
FastMath.toDegrees(gpA.getLatitude()),
FastMath.toDegrees(gpA.getLongitude()));
// Rugged A initialization
// ---------------------
System.out.format("**** Rugged A initialization **** %n");
RuggedBuilder ruggedBuilderA = new RuggedBuilder();
ruggedBuilderA.addLineSensor(lineSensorA);
ruggedBuilderA.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID);
ruggedBuilderA.setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF);
ruggedBuilderA.setTimeSpan(minDateA,maxDateA, 0.001, 5.0);
ruggedBuilderA.setTrajectory(InertialFrameId.EME2000, satellitePVListA, nbPVPoints,
CartesianDerivativesFilter.USE_PV, satelliteQListA,
nbQPoints, AngularDerivativesFilter.USE_R);
ruggedBuilderA.setLightTimeCorrection(false);
ruggedBuilderA.setAberrationOfLightCorrection(false);
ruggedBuilderA.setName("RuggedA");
Rugged ruggedA = ruggedBuilderA.build();
System.out.format("\n**** Build Pleiades viewing model B and its orbit definition **** %n");
// 2/- Create Second Pleiades Viewing Model
final AbsoluteDate minDateB = pleiadesViewingModelB.getMinDate();
final AbsoluteDate maxDateB = pleiadesViewingModelB.getMaxDate();
final AbsoluteDate refDateB = pleiadesViewingModelB.getDatationReference();
LineSensor lineSensorB = pleiadesViewingModelB.getLineSensor();
System.out.format(Locale.US, "Viewing model B - date min: %s max: %s ref: %s %n", minDateB, maxDateB, refDateB);
// ----Satellite position, velocity and attitude: create orbit model B
OrbitModel orbitmodelB = new OrbitModel();
BodyShape earthB = orbitmodelB.createEarth();
Orbit orbitB = orbitmodelB.createOrbit(Constants.EIGEN5C_EARTH_MU, refDateB);
// ----Satellite attitude
List<TimeStampedAngularCoordinates> satelliteQListB = orbitmodelB.orbitToQ(orbitB, earthB, minDateB.shiftedBy(-0.0), maxDateB.shiftedBy(+0.0), 0.25);
// ----Position and velocities
PVCoordinates PVB = orbitB.getPVCoordinates(earthB.getBodyFrame());
List<TimeStampedPVCoordinates> satellitePVListB = orbitmodelB.orbitToPV(orbitB, earthB, minDateB.shiftedBy(-0.0), maxDateB.shiftedBy(+0.0), 0.25);
GeodeticPoint gpB = earthB.transform(PVB.getPosition(), earthB.getBodyFrame(), orbitB.getDate());
System.out.format(Locale.US, "(B) Geodetic Point at date %s : φ = %8.10f °, λ = %8.10f %n",orbitA.getDate().toString(),
FastMath.toDegrees(gpB.getLatitude()),
FastMath.toDegrees(gpB.getLongitude()));
// Rugged B initialization
// ---------------------
System.out.format("**** Rugged B initialization **** %n");
RuggedBuilder ruggedBuilderB = new RuggedBuilder();
ruggedBuilderB.addLineSensor(lineSensorB);
ruggedBuilderB.setAlgorithm(AlgorithmId.IGNORE_DEM_USE_ELLIPSOID);
ruggedBuilderB.setEllipsoid(EllipsoidId.WGS84, BodyRotatingFrameId.ITRF);
ruggedBuilderB.setTimeSpan(minDateB,maxDateB, 0.001, 5.0);
ruggedBuilderB.setTrajectory(InertialFrameId.EME2000, satellitePVListB, nbPVPoints,
CartesianDerivativesFilter.USE_PV, satelliteQListB,
nbQPoints, AngularDerivativesFilter.USE_R);
ruggedBuilderB.setLightTimeCorrection(false);
ruggedBuilderB.setAberrationOfLightCorrection(false);
ruggedBuilderB.setName("RuggedB");
Rugged ruggedB = ruggedBuilderB.build();
// Compute distance between LOS
// -----------------------------
double distance = refining.computeDistanceBetweenLOS(ruggedA, lineSensorA, ruggedB, lineSensorB);
System.out.format("distance %f meters %n",distance);
// Initialize disruptions:
// -----------------------
// Introduce rotations around instrument axes (roll and pitch translations, scale factor)
System.out.format("\n**** Add disruptions on both acquisitions (A,B): roll and pitch rotations, scale factor **** %n");
double rollValueA = FastMath.toRadians(0.004);
double pitchValueA = FastMath.toRadians(0.0008);
double rollValueB = FastMath.toRadians(-0.004);
double pitchValueB = FastMath.toRadians(0.0008);
double factorValue = 1.0;
System.out.format("Acquisition A roll: %3.5f \tpitch: %3.5f \tfactor: %3.5f \n",rollValueA,pitchValueA,factorValue);
System.out.format("Acquisition B roll: %3.5f \tpitch: %3.5f \tfactor: %3.5f \n",rollValueB,pitchValueB,factorValue);
// Apply disruptions on physical model (acquisition A)
refining.applyDisruptions(ruggedA, sensorNameA, rollValueA, pitchValueA, factorValue);
// Apply disruptions on physical model (acquisition B)
refining.applyDisruptions(ruggedB, sensorNameB, rollValueB, pitchValueB, factorValue);
// Generate measurements (observations) from physical model disrupted
// ------------------------------------------------------------------
int lineSampling = 1000;
int pixelSampling = 1000;
// Noise definition
// distribution: gaussian(0), vector dimension: 2
final Noise noise = new Noise(0,2);
// {pixelA mean, pixelB mean}
final double[] mean = {5.0,0.0};
// {pixelB std, pixelB std}
final double[] standardDeviation = {0.1,0.1};
noise.setMean(mean);
noise.setStandardDeviation(standardDeviation);
System.out.format("\tSensor A mean: %f std %f %n",mean[0],standardDeviation[0]);
System.out.format("\tSensor B mean: %f std %f %n",mean[1],standardDeviation[1]);
// Inter-measurements (between both viewing models)
InterMeasurementGenerator measurements = refining.generateNoisyPoints(lineSampling, pixelSampling,
ruggedA, sensorNameA,
pleiadesViewingModelA.getDimension(),
ruggedB, sensorNameB,
pleiadesViewingModelB.getDimension(),
noise);
// Compute ground truth residuals
// ------------------------------
System.out.format("\n**** Ground truth residuals **** %n");
refining.computeMetrics(measurements.getInterMapping(), ruggedA, ruggedB, false);
// Initialize physical models without disruptions
// -----------------------------------------------
System.out.format("\n**** Initialize physical models (A,B) without disruptions: reset Roll/Pitch/Factor **** %n");
refining.resetModel(ruggedA, sensorNameA, false);
refining.resetModel(ruggedB, sensorNameB, false);
// Compute initial residuals
// -------------------------
System.out.format("\n**** Initial Residuals **** %n");
refining.computeMetrics(measurements.getInterMapping(), ruggedA, ruggedB, false);
// Start optimization
// ------------------
System.out.format("\n**** Start optimization **** %n");
final int maxIterations = 100;
final double convergenceThreshold = 1.e-7;
List<Rugged> ruggedList = Arrays.asList(ruggedA, ruggedB);
refining.optimization(maxIterations, convergenceThreshold,
measurements.getObservables(), ruggedList);
// Check estimated values
// ----------------------
System.out.format("\n**** Check parameters ajustement **** %n");
System.out.format("Acquisition A:%n");
refining.paramsEstimation(ruggedA, sensorNameA, rollValueA, pitchValueA, factorValue);
System.out.format("Acquisition B:%n");
refining.paramsEstimation(ruggedB, sensorNameB, rollValueB, pitchValueB, factorValue);
// Compute statistics
// ------------------
System.out.format("\n**** Compute Statistics **** %n");
refining.computeMetrics(measurements.getInterMapping(), ruggedA, ruggedB, false);
} catch (OrekitException oe) {
System.err.println(oe.getLocalizedMessage());
System.exit(1);
} catch (RuggedException re) {
System.err.println(re.getLocalizedMessage());
System.exit(1);
}
}
/** Constructor */
private InterRefining() {
sensorNameA = "SensorA";
final double incidenceAngleA = -5.0;
final String dateA = "2016-01-01T11:59:50.0";
pleiadesViewingModelA = new PleiadesViewingModel(sensorNameA, incidenceAngleA, dateA);
sensorNameB = "SensorB";
final double incidenceAngleB = 0.0;
final String dateB = "2016-01-01T12:02:50.0";
pleiadesViewingModelB = new PleiadesViewingModel(sensorNameB, incidenceAngleB, dateB);
}
/** Estimate distance between LOS
* @param lineSensorA line sensor A
* @param lineSensorB line sensor B
* @return GSD
*/
private double computeDistanceBetweenLOS(final Rugged ruggedA, final LineSensor lineSensorA,
final Rugged ruggedB, final LineSensor lineSensorB) {
// Get number of line of sensors
int dimensionA = pleiadesViewingModelA.getDimension();
int dimensionB = pleiadesViewingModelB.getDimension();
Vector3D positionA = lineSensorA.getPosition();
// This returns a zero vector since we set the relative position of the sensor w.r.T the satellite to 0.
AbsoluteDate lineDateA = lineSensorA.getDate(dimensionA/2);
Vector3D losA = lineSensorA.getLOS(lineDateA,dimensionA/2);
GeodeticPoint centerPointA = ruggedA.directLocation(lineDateA, positionA, losA);
System.out.format(Locale.US, "\ncenter geodetic position A : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(centerPointA.getLatitude()),
FastMath.toDegrees(centerPointA.getLongitude()),centerPointA.getAltitude());
Vector3D positionB = lineSensorB.getPosition();
// This returns a zero vector since we set the relative position of the sensor w.r.T the satellite to 0.
AbsoluteDate lineDateB = lineSensorB.getDate(dimensionB/2);
Vector3D losB = lineSensorB.getLOS(lineDateB,dimensionB/2);
GeodeticPoint centerPointB = ruggedB.directLocation(lineDateB, positionB, losB);
System.out.format(Locale.US, "center geodetic position B : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(centerPointB.getLatitude()),
FastMath.toDegrees(centerPointB.getLongitude()),centerPointB.getAltitude());
lineDateB = lineSensorB.getDate(0);
losB = lineSensorB.getLOS(lineDateB,0);
GeodeticPoint firstPointB = ruggedB.directLocation(lineDateB, positionB, losB);
System.out.format(Locale.US, "first geodetic position B : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(firstPointB.getLatitude()),
FastMath.toDegrees(firstPointB.getLongitude()),firstPointB.getAltitude());
lineDateB = lineSensorB.getDate(dimensionB-1);
losB = lineSensorB.getLOS(lineDateB,dimensionB-1);
GeodeticPoint lastPointB = ruggedB.directLocation(lineDateB, positionB, losB);
System.out.format(Locale.US, "last geodetic position B : φ = %8.10f °, λ = %8.10f °, h = %8.3f m%n",
FastMath.toDegrees(lastPointB.getLatitude()),
FastMath.toDegrees(lastPointB.getLongitude()),lastPointB.getAltitude());
double distance = DistanceTools.computeDistanceInMeter(centerPointA, centerPointB);
return distance;
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/Refining.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package fr.cs.examples.refiningPleiades;
import java.util.ArrayList;
import java.util.Collection;
import java.util.Collections;
import java.util.List;
import org.hipparchus.optim.nonlinear.vector.leastsquares.LeastSquaresOptimizer.Optimum;
import org.orekit.rugged.adjustment.AdjustmentContext;
import org.orekit.rugged.adjustment.measurements.Observables;
import org.orekit.rugged.adjustment.measurements.SensorToGroundMapping;
import org.orekit.rugged.adjustment.measurements.SensorToSensorMapping;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.errors.RuggedMessages;
import fr.cs.examples.refiningPleiades.generators.GroundMeasurementGenerator;
import fr.cs.examples.refiningPleiades.generators.InterMeasurementGenerator;
import fr.cs.examples.refiningPleiades.generators.Noise;
import fr.cs.examples.refiningPleiades.metrics.LocalisationMetrics;
/**
* Class for refining problems common methods
* @author Jonathan Guinet
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @see SensorToGroundMapping
* @see SensorToSensorMapping
* @see GroundMeasurementGenerator
* @see InterMeasurementGenerator
* @since 2.0
*/
public class Refining {
/**
* Constructor
*/
public Refining() {
}
/** Apply disruptions on acquisition for roll, pitch and scale factor
* @param rugged Rugged instance
* @param sensorName line sensor name
* @param rollValue rotation on roll value
* @param pitchValue rotation on pitch value
* @param factorValue scale factor
*/
public void applyDisruptions(final Rugged rugged, final String sensorName,
final double rollValue, final double pitchValue, final double factorValue) {
rugged.
getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(sensorName + rollSuffix)).
findFirst().get().setValue(rollValue);
rugged.
getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(sensorName + pitchSuffix)).
findFirst().get().setValue(pitchValue);
rugged.
getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(factorName)).
findFirst().get().setValue(factorValue);
}
/** Generate measurements without noise (sensor to ground mapping)
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
* @param rugged Rugged instance
* @param sensorName line sensor name
* @param dimension number of line of the sensor
* @return ground measurements generator (sensor to ground mapping)
*/
public GroundMeasurementGenerator generatePoints(final int lineSampling, final int pixelSampling,
final Rugged rugged, final String sensorName,
final int dimension) {
GroundMeasurementGenerator measurements = new GroundMeasurementGenerator(rugged, sensorName, dimension);
System.out.format("\n**** Generate measurements (without noise; sensor to ground mapping) **** %n");
// Generation measurements without noise
measurements.createMeasurement(lineSampling, pixelSampling);
System.out.format("Number of tie points generated: %d %n", measurements.getMeasurementCount());
return measurements;
}
/** Generate measurements without noise (sensor to sensor mapping)
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
* @param ruggedA Rugged instance of acquisition A
* @param sensorNameA line sensor name A
* @param dimensionA dimension for acquisition A
* @param ruggedB Rugged instance of acquisition B
* @param sensorNameB line sensor name B
* @param dimensionB dimension for acquisition B
* @return inter measurements generator (sensor to sensor mapping)
*/
public InterMeasurementGenerator generatePoints(final int lineSampling, final int pixelSampling,
final Rugged ruggedA, final String sensorNameA, final int dimensionA,
final Rugged ruggedB, final String sensorNameB, final int dimensionB) {
// Outliers control
final double outlierValue = 1e+2;
// Earth constraint weight
final double earthConstraintWeight = 0.1;
// Generate measurements with constraints on outliers control and Earth distance
InterMeasurementGenerator measurements = new InterMeasurementGenerator(ruggedA, sensorNameA, dimensionA,
ruggedB, sensorNameB, dimensionB,
outlierValue,
earthConstraintWeight);
System.out.format("\n**** Generate measurements (without noise; sensor to sensor mapping) **** %n");
// Generation measurements without noise
measurements.createMeasurement(lineSampling, pixelSampling);
System.out.format("Number of tie points generated: %d %n", measurements.getMeasurementCount());
return measurements;
}
/** Generate noisy measurements (sensor to ground mapping)
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
* @param rugged Rugged instance
* @param sensorName line sensor name
* @param dimension dimension
* @param noise Noise structure to generate noisy measurements
* @return ground measurements generator (sensor to ground mapping)
*/
public GroundMeasurementGenerator generateNoisyPoints(final int lineSampling, final int pixelSampling,
final Rugged rugged, final String sensorName, final int dimension,
final Noise noise) {
// Generate ground measurements
GroundMeasurementGenerator measurements = new GroundMeasurementGenerator(rugged, sensorName, dimension);
System.out.format("\n**** Generate noisy measurements (sensor to ground mapping) **** %n");
// Generate noisy measurements
measurements.createNoisyMeasurement(lineSampling, pixelSampling, noise);
System.out.format("Number of tie points generated: %d %n", measurements.getMeasurementCount());
return measurements;
}
/** Generate noisy measurements (sensor to sensor mapping)
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
* @param ruggedA Rugged instance of acquisition A
* @param sensorNameA line sensor name A
* @param dimensionA dimension for acquisition A
* @param ruggedB Rugged instance of acquisition B
* @param sensorNameB line sensor name B
* @param dimensionB dimension for acquisition B
* @param noise noise structure to generate noisy measurements
* @return inter-measurements generator (sensor to sensor mapping)
*/
public InterMeasurementGenerator generateNoisyPoints(final int lineSampling, final int pixelSampling,
final Rugged ruggedA, final String sensorNameA, final int dimensionA,
final Rugged ruggedB, final String sensorNameB, final int dimensionB,
final Noise noise) {
// Outliers control
final double outlierValue = 1.e+2;
// Earth constraint weight
final double earthConstraintWeight = 0.1;
// Generate measurements with constraints on Earth distance and outliers control
InterMeasurementGenerator measurements = new InterMeasurementGenerator(ruggedA, sensorNameA, dimensionA,
ruggedB, sensorNameB, dimensionB,
outlierValue,
earthConstraintWeight);
System.out.format("\n**** Generate noisy measurements (sensor to sensor mapping) **** %n");
// Generation noisy measurements
measurements.createNoisyMeasurement(lineSampling, pixelSampling, noise);
System.out.format("Number of tie points generated: %d %n", measurements.getMeasurementCount());
return measurements;
}
/** Compute metrics to evaluate geometric performances in location,
* for Ground Control Points GCP study.
* @param groundMapping sensor to ground mapping
* @param rugged Rugged instance
* @param unit flag to know if distance is computed in meters (false) or with angular (true)
*/
public void computeMetrics(final SensorToGroundMapping groundMapping,
final Rugged rugged, final boolean unit) {
String stUnit = null;
if(unit) {
stUnit="degrees";
} else {
stUnit="meters";
}
LocalisationMetrics residuals = new LocalisationMetrics(groundMapping, rugged, unit);
System.out.format("Max: %3.4e Mean: %3.4e %s%n",residuals.getMaxResidual(),residuals.getMeanResidual(), stUnit);
}
/** Compute metrics to evaluate geometric performances in location,
* for tie points study.
* @param interMapping sensor to sensor mapping
* @param ruggedA Rugged instance A
* @param ruggedB Rugged instance B
* @param unit flag to know if distance is computed in meters (false) or with angular (true)
*/
public void computeMetrics(final SensorToSensorMapping interMapping,
final Rugged ruggedA, final Rugged ruggedB,
final boolean unit) {
String stUnit = null;
if(unit) stUnit="degrees";
else stUnit="meters";
LocalisationMetrics residuals = new LocalisationMetrics(interMapping, ruggedA, ruggedB, unit);
System.out.format("Max: %1.4e Mean: %1.4e %s%n",residuals.getMaxResidual(),residuals.getMeanResidual(), stUnit);
System.out.format("LOS distance Max: %1.4e Mean: %1.4e %s%n",residuals.getLosMaxDistance(),residuals.getLosMeanDistance(), stUnit);
System.out.format("Earth distance Max: %1.4e Mean: %1.4e %s%n",residuals.getEarthMaxDistance(),residuals.getEarthMeanDistance(), stUnit);
}
/** Reset a model
* @param rugged Rugged instance
* @param sensorName line sensor name
* @param isSelected flag to known if factor parameter is selected or not
*/
public void resetModel(final Rugged rugged, final String sensorName, final boolean isSelected) {
rugged.
getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(sensorName + rollSuffix)
|| driver.getName().equals(sensorName + pitchSuffix)).
forEach(driver -> {
driver.setSelected(true);
driver.setValue(0.0);
});
rugged.
getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(factorName)).
forEach(driver -> {
driver.setSelected(isSelected);
// default value: no Z scale factor applied
driver.setValue(1.0);
});
}
/** Start optimization to adjust parameters (Ground Control Points GCP study).
* @param maxIterations iterations max
* @param convergenceThreshold threshold of convergence
* @param measurements ground measurements
* @param rugged Rugged instance
*/
public void optimization(final int maxIterations, final double convergenceThreshold,
final Observables measurements, final Rugged rugged) {
System.out.format("Iterations max: %d\tconvergence threshold: %3.6e \n", maxIterations, convergenceThreshold);
AdjustmentContext adjustmentContext = new AdjustmentContext(Collections.singletonList(rugged), measurements);
Optimum optimum = adjustmentContext.estimateFreeParameters(Collections.singletonList(rugged.getName()), maxIterations, convergenceThreshold);
// Print statistics
System.out.format("Max value: %3.6e %n", optimum.getResiduals().getMaxValue());
System.out.format("Optimization performed in %d iterations \n", optimum.getIterations());
System.out.format("Optimization performed with %d evaluations of model (objective) function \n", optimum.getEvaluations());
System.out.format("RMSE: %f \n", optimum.getRMS());
}
/** Start optimization to adjust parameters (tie points study).
* @param maxIterations iterations max
* @param convergenceThreshold threshold of convergence
* @param measurements measurements
* @param ruggeds Rugged instances A and B
*/
public void optimization(final int maxIterations, final double convergenceThreshold,
final Observables measurements,
final Collection<Rugged> ruggeds) {
System.out.format("Iterations max: %d\tconvergence threshold: %3.6e \n", maxIterations, convergenceThreshold);
if(ruggeds.size()!= 2 ) {
throw new RuggedException(RuggedMessages.UNSUPPORTED_REFINING_CONTEXT,ruggeds.size());
}
AdjustmentContext adjustmentContext = new AdjustmentContext(ruggeds, measurements);
List<String> ruggedNameList = new ArrayList<String>();
for(Rugged rugged : ruggeds) {
ruggedNameList.add(rugged.getName());
}
Optimum optimum = adjustmentContext.estimateFreeParameters(ruggedNameList, maxIterations, convergenceThreshold);
// Print statistics
System.out.format("Max value: %3.6e %n", optimum.getResiduals().getMaxValue());
System.out.format("Optimization performed in %d iterations \n", optimum.getIterations());
System.out.format("Optimization performed with %d evaluations of model (objective) function \n", optimum.getEvaluations());
System.out.format("RMSE: %f \n", optimum.getRMS());
}
/** Check adjusted parameters of an acquisition
* @param rugged Rugged instance
* @param sensorName line sensor name
* @param rollValue rotation on roll value
* @param pitchValue rotation on pitch value
* @param factorValue scale factor
*/
public void paramsEstimation(final Rugged rugged, final String sensorName,
final double rollValue, final double pitchValue, final double factorValue) {
// Estimate Roll
double estimatedRoll = rugged.getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(sensorName + rollSuffix)).
findFirst().get().getValue();
double rollError = (estimatedRoll - rollValue);
System.out.format("Estimated roll: %3.5f\troll error: %3.6e %n", estimatedRoll, rollError);
// Estimate pitch
double estimatedPitch = rugged.getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(sensorName + pitchSuffix)).
findFirst().get().getValue();
double pitchError = (estimatedPitch - pitchValue);
System.out.format("Estimated pitch: %3.5f\tpitch error: %3.6e %n", estimatedPitch, pitchError);
// Estimate scale factor
double estimatedFactor = rugged.getLineSensor(sensorName).
getParametersDrivers().
filter(driver -> driver.getName().equals(factorName)).
findFirst().get().getValue();
double factorError = (estimatedFactor - factorValue);
System.out.format("Estimated factor: %3.5f\tfactor error: %3.6e %n", estimatedFactor, factorError);
}
/**
* Part of the name of parameter drivers
*/
static final String rollSuffix = "_roll";
static final String pitchSuffix = "_pitch";
static final String factorName = "factor";
public static String getRollsuffix() {
return rollSuffix;
}
public static String getPitchsuffix() {
return pitchSuffix;
}
public static String getFactorname() {
return factorName;
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/generators/GroundMeasurementGenerator.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package fr.cs.examples.refiningPleiades.generators;
import org.hipparchus.geometry.euclidean.threed.Vector3D;
import org.hipparchus.random.GaussianRandomGenerator;
import org.hipparchus.random.UncorrelatedRandomVectorGenerator;
import org.hipparchus.random.Well19937a;
import org.hipparchus.util.FastMath;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.adjustment.measurements.Observables;
import org.orekit.rugged.adjustment.measurements.SensorToGroundMapping;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.time.AbsoluteDate;
/** Ground measurements generator (sensor to ground mapping).
* @author Jonathan Guinet
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @since 2.0
*/
public class GroundMeasurementGenerator implements Measurable {
/** Sensor to ground mapping. */
private SensorToGroundMapping groundMapping;
/** Observables which contains ground mapping. */
private Observables observables;
/** Rugged instance */
private Rugged rugged;
/** Line sensor */
private LineSensor sensor;
/** Number of lines of the sensor */
private int dimension;
/** Number of measurements */
private int measurementCount;
/** Simple constructor.
* @param rugged Rugged instance
* @param sensorName sensor name
* @param dimension number of line of the sensor
*/
public GroundMeasurementGenerator(final Rugged rugged, final String sensorName, final int dimension) {
// Generate reference mapping
this.groundMapping = new SensorToGroundMapping(rugged.getName(), sensorName);
// Create observables for one model
this.observables = new Observables(1);
this.rugged = rugged;
this.sensor = rugged.getLineSensor(groundMapping.getSensorName());
this.dimension = dimension;
this.measurementCount = 0;
}
/** Get the sensor to ground mapping
* @return the sensor to ground mapping
*/
public SensorToGroundMapping getGroundMapping() {
return groundMapping;
}
/** Get the observables which contains the ground mapping
* @return the observables which contains the ground mapping
*/
public Observables getObservables() {
return observables;
}
@Override
public int getMeasurementCount() {
return measurementCount;
}
@Override
public void createMeasurement(final int lineSampling, final int pixelSampling) {
for (double line = 0; line < dimension; line += lineSampling) {
final AbsoluteDate date = sensor.getDate(line);
for (int pixel = 0; pixel < sensor.getNbPixels(); pixel += pixelSampling) {
final GeodeticPoint gp2 = rugged.directLocation(date, sensor.getPosition(),
sensor.getLOS(date, pixel));
groundMapping.addMapping(new SensorPixel(line, pixel), gp2);
// increment the number of measurements
measurementCount++;
}
}
observables.addGroundMapping(groundMapping);
}
@Override
public void createNoisyMeasurement(final int lineSampling, final int pixelSampling, final Noise noise) {
// Estimate latitude and longitude errors (rad)
final Vector3D latLongError = estimateLatLongError();
// Get noise features
final double[] mean = noise.getMean(); // [latitude, longitude, altitude] mean
final double[] std = noise.getStandardDeviation(); // [latitude, longitude, altitude] standard deviation
final double latErrorMean = mean[0] * latLongError.getX();
final double lonErrorMean = mean[1] * latLongError.getY();
final double latErrorStd = std[0] * latLongError.getX();
final double lonErrorStd = std[1] * latLongError.getY();
// Gaussian random generator
// Build a null mean random uncorrelated vector generator with standard deviation corresponding to the estimated error on ground
final double meanGenerator[] = {latErrorMean, lonErrorMean, mean[2]};
final double stdGenerator[] = {latErrorStd, lonErrorStd, std[2]};
// seed has been fixed for tests purpose
final GaussianRandomGenerator rng = new GaussianRandomGenerator(new Well19937a(0xefac03d9be4d24b9l));
final UncorrelatedRandomVectorGenerator rvg = new UncorrelatedRandomVectorGenerator(meanGenerator, stdGenerator, rng);
for (double line = 0; line < dimension; line += lineSampling) {
final AbsoluteDate date = sensor.getDate(line);
for (int pixel = 0; pixel < sensor.getNbPixels(); pixel += pixelSampling) {
// Components of generated vector follow (independent) Gaussian distribution
final Vector3D vecRandom = new Vector3D(rvg.nextVector());
final GeodeticPoint gp2 = rugged.directLocation(date, sensor.getPosition(),
sensor.getLOS(date, pixel));
final GeodeticPoint gpNoisy = new GeodeticPoint(gp2.getLatitude() + vecRandom.getX(),
gp2.getLongitude() + vecRandom.getY(),
gp2.getAltitude() + vecRandom.getZ());
groundMapping.addMapping(new SensorPixel(line, pixel), gpNoisy);
// increment the number of measurements
measurementCount++;
}
}
this.observables.addGroundMapping(groundMapping);
}
/** Compute latitude and longitude errors
* @return the latitude and longitude errors (rad)
*/
private Vector3D estimateLatLongError() {
final int pix = sensor.getNbPixels() / 2;
final int line = (int) FastMath.floor(pix); // assumption : same number of line and pixels;
final AbsoluteDate date = sensor.getDate(line);
final GeodeticPoint gp_pix0 = rugged.directLocation(date, sensor.getPosition(), sensor.getLOS(date, pix));
final AbsoluteDate date1 = sensor.getDate(line + 1);
final GeodeticPoint gp_pix1 = rugged.directLocation(date1, sensor.getPosition(), sensor.getLOS(date1, pix + 1));
final double latErr = FastMath.abs(gp_pix0.getLatitude() - gp_pix1.getLatitude());
final double lonErr = FastMath.abs(gp_pix0.getLongitude() - gp_pix1.getLongitude());
return new Vector3D(latErr, lonErr, 0.0);
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/generators/InterMeasurementGenerator.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package fr.cs.examples.refiningPleiades.generators;
import org.hipparchus.random.GaussianRandomGenerator;
import org.hipparchus.random.UncorrelatedRandomVectorGenerator;
import org.hipparchus.random.Well19937a;
import org.orekit.bodies.GeodeticPoint;
import org.orekit.rugged.api.Rugged;
import org.orekit.rugged.errors.RuggedException;
import org.orekit.rugged.errors.RuggedMessages;
import org.orekit.rugged.linesensor.LineSensor;
import org.orekit.rugged.linesensor.SensorPixel;
import org.orekit.rugged.adjustment.measurements.Observables;
import org.orekit.rugged.adjustment.measurements.SensorToSensorMapping;
import org.orekit.rugged.utils.SpacecraftToObservedBody;
import org.orekit.time.AbsoluteDate;
import fr.cs.examples.refiningPleiades.metrics.DistanceTools;
/**
* Inter-measurements generator (sensor to sensor mapping).
* @author Jonathan Guinet
* @author Lucie Labatallee
* @author Guylaine Prat
* @since 2.0
*/
public class InterMeasurementGenerator implements Measurable {
/** Mapping from sensor A to sensor B. */
private SensorToSensorMapping interMapping;
/** Observables which contains sensor to sensor mapping.*/
private Observables observables;
/** Rugged instance corresponding to the viewing model A */
private Rugged ruggedA;
/** Rugged instance corresponding to the viewing model B */
private Rugged ruggedB;
/** Sensor A */
private LineSensor sensorA;
/** Sensor B */
private LineSensor sensorB;
/** Number of measurements */
private int measurementCount;
/** Sensor name B */
private String sensorNameB;
/** Number of line for acquisition A */
private int dimensionA;
/** Number of line for acquisition B */
private int dimensionB;
/** Limit value for outlier points */
private double outlier;
/** Default constructor: measurements generation without outlier points control
* and without Earth distance constraint.
* @param ruggedA Rugged instance corresponding to the viewing model A
* @param sensorNameA sensor name A
* @param dimensionA number of line for acquisition A
* @param ruggedB Rugged instance corresponding to the viewing model B
* @param sensorNameB sensor name B
* @param dimensionB number of line for acquisition B
*/
public InterMeasurementGenerator(final Rugged ruggedA, final String sensorNameA, final int dimensionA,
final Rugged ruggedB, final String sensorNameB, final int dimensionB) {
// Initialize parameters
initParams(ruggedA, sensorNameA, dimensionA, ruggedB, sensorNameB, dimensionB);
// Generate reference mapping, without Earth distance constraint
interMapping = new SensorToSensorMapping(sensorNameA, sensorNameB);
// Create observables for two models
observables = new Observables(2);
}
/** Constructor for measurements generation taking into account outlier points control,
* without Earth distance constraint.
* @param ruggedA Rugged instance corresponding to the viewing model A
* @param sensorNameA sensor name A
* @param dimensionA dimension for acquisition A
* @param ruggedB Rugged instance corresponding to the viewing model B
* @param sensorNameB sensor name B
* @param dimensionB dimension for acquisition B
* @param outlier limit value for outlier points
*/
public InterMeasurementGenerator(final Rugged ruggedA, final String sensorNameA, final int dimensionA,
final Rugged ruggedB, final String sensorNameB, final int dimensionB,
final double outlier) {
this(ruggedA, sensorNameA, dimensionA, ruggedB, sensorNameB, dimensionB);
this.outlier = outlier;
}
/** Constructor for measurements generation taking into account outlier points control,
* and Earth distance constraint.
* @param ruggedA Rugged instance corresponding to the viewing model A
* @param sensorNameA sensor name A
* @param dimensionA dimension for acquisition A
* @param ruggedB Rugged instance corresponding to the viewing model B
* @param sensorNameB sensor name B
* @param dimensionB dimension for acquisition B
* @param outlier limit value for outlier points
* @param earthConstraintWeight weight given to the Earth distance constraint
* with respect to the LOS distance (between 0 and 1).
*/
public InterMeasurementGenerator(final Rugged ruggedA, final String sensorNameA, final int dimensionA,
final Rugged ruggedB, final String sensorNameB, final int dimensionB,
final double outlier, final double earthConstraintWeight) {
// Initialize parameters
initParams(ruggedA, sensorNameA, dimensionA, ruggedB, sensorNameB, dimensionB);
// Generate reference mapping, with Earth distance constraints.
// Weighting will be applied as follow :
// (1-earthConstraintWeight) for losDistance weighting
// earthConstraintWeight for earthDistance weighting
interMapping = new SensorToSensorMapping(sensorNameA, ruggedA.getName(), sensorNameB, ruggedB.getName(), earthConstraintWeight);
// Outlier points control
this.outlier = outlier;
// Observables which contains sensor to sensor mapping
this.observables = new Observables(2);
}
/** Get the mapping from sensor A to sensor B
* @return the mapping from sensor A to sensor B
*/
public SensorToSensorMapping getInterMapping() {
return interMapping;
}
/** Get the observables which contains sensor to sensor mapping
* @return the observables which contains sensor to sensor mapping
*/
public Observables getObservables() {
return observables;
}
@Override
public int getMeasurementCount() {
return measurementCount;
}
@Override
public void createMeasurement(final int lineSampling, final int pixelSampling) {
// Search the sensor pixel seeing point
final int minLine = 0;
final int maxLine = dimensionB - 1;
for (double line = 0; line < dimensionA; line += lineSampling) {
final AbsoluteDate dateA = sensorA.getDate(line);
for (double pixelA = 0; pixelA < sensorA.getNbPixels(); pixelA += pixelSampling) {
final GeodeticPoint gpA = ruggedA.directLocation(dateA, sensorA.getPosition(),
sensorA.getLOS(dateA, pixelA));
final SensorPixel sensorPixelB = ruggedB.inverseLocation(sensorNameB, gpA, minLine, maxLine);
// We need to test if the sensor pixel is found in the prescribed lines
// otherwise the sensor pixel is null
if (sensorPixelB != null) {
final AbsoluteDate dateB = sensorB.getDate(sensorPixelB.getLineNumber());
final double pixelB = sensorPixelB.getPixelNumber();
final SpacecraftToObservedBody scToBodyA = ruggedA.getScToBody();
final GeodeticPoint gpB = ruggedB.directLocation(dateB, sensorB.getPosition(),
sensorB.getLOS(dateB, pixelB));
final double geoDistance = DistanceTools.computeDistanceInMeter(gpA, gpB);
if (geoDistance < this.outlier) {
final SensorPixel realPixelA = new SensorPixel(line, pixelA);
final SensorPixel realPixelB = new SensorPixel(sensorPixelB.getLineNumber(), sensorPixelB.getPixelNumber());
final AbsoluteDate realDateA = sensorA.getDate(realPixelA.getLineNumber());
final AbsoluteDate realDateB = sensorB.getDate(realPixelB.getLineNumber());
final double[] distanceLOSB = ruggedB.distanceBetweenLOS(sensorA, realDateA, realPixelA.getPixelNumber(), scToBodyA,
sensorB, realDateB, realPixelB.getPixelNumber());
final double losDistance = 0.0;
final double earthDistance = distanceLOSB[1];
interMapping.addMapping(realPixelA, realPixelB, losDistance, earthDistance);
// Increment the number of measurements
this.measurementCount++;
}
}
}
}
this.observables.addInterMapping(interMapping);
}
/** Tie point creation from direct 1ocation with Rugged A and inverse location with Rugged B
* @param lineSampling sampling along lines
* @param pixelSampling sampling along columns
* @param noise errors to apply to measure for pixel A and pixel B
*/
public void createNoisyMeasurement(final int lineSampling, final int pixelSampling, final Noise noise) {
// Get noise features (errors)
// [pixelA, pixelB] mean
final double[] mean = noise.getMean();
// [pixelA, pixelB] standard deviation
final double[] std = noise.getStandardDeviation();
// Search the sensor pixel seeing point
final int minLine = 0;
final int maxLine = dimensionB - 1;
final double meanA[] = { mean[0], mean[0] };
final double stdA[] = { std[0], std[0] };
final double meanB[] = { mean[1], mean[1] };
final double stdB[] = { std[1], std[1] };
// Seed has been fixed for tests purpose
final GaussianRandomGenerator rngA = new GaussianRandomGenerator(new Well19937a(0xefac03d9be4d24b9l));
final UncorrelatedRandomVectorGenerator rvgA = new UncorrelatedRandomVectorGenerator(meanA, stdA, rngA);
// Seed has been fixed for tests purpose
final GaussianRandomGenerator rngB = new GaussianRandomGenerator(new Well19937a(0xdf1c03d9be0b34b9l));
final UncorrelatedRandomVectorGenerator rvgB = new UncorrelatedRandomVectorGenerator(meanB, stdB, rngB);
for (double line = 0; line < dimensionA; line += lineSampling) {
final AbsoluteDate dateA = sensorA.getDate(line);
for (double pixelA = 0; pixelA < sensorA.getNbPixels(); pixelA += pixelSampling) {
final GeodeticPoint gpA = ruggedA.directLocation(dateA, sensorA.getPosition(),
sensorA.getLOS(dateA, pixelA));
final SensorPixel sensorPixelB = ruggedB.inverseLocation(sensorNameB, gpA, minLine, maxLine);
// We need to test if the sensor pixel is found in the prescribed lines
// otherwise the sensor pixel is null
if (sensorPixelB != null) {
final AbsoluteDate dateB = sensorB.getDate(sensorPixelB.getLineNumber());
final double pixelB = sensorPixelB.getPixelNumber();
// Get spacecraft to body transform of Rugged instance A
final SpacecraftToObservedBody scToBodyA = ruggedA.getScToBody();
final GeodeticPoint gpB = ruggedB.directLocation(dateB, sensorB.getPosition(),
sensorB.getLOS(dateB, pixelB));
final double geoDistance = DistanceTools.computeDistanceInMeter(gpA, gpB);
// Create the inter mapping if distance is below outlier value
if (geoDistance < outlier) {
final double[] vecRandomA = rvgA.nextVector();
final double[] vecRandomB = rvgB.nextVector();
final SensorPixel realPixelA = new SensorPixel(line + vecRandomA[0], pixelA + vecRandomA[1]);
final SensorPixel realPixelB = new SensorPixel(sensorPixelB.getLineNumber() + vecRandomB[0],
sensorPixelB.getPixelNumber() + vecRandomB[1]);
final AbsoluteDate realDateA = sensorA.getDate(realPixelA.getLineNumber());
final AbsoluteDate realDateB = sensorB.getDate(realPixelB.getLineNumber());
final double[] distanceLOSB = ruggedB.distanceBetweenLOS(sensorA, realDateA, realPixelA.getPixelNumber(), scToBodyA,
sensorB, realDateB, realPixelB.getPixelNumber());
final Double losDistance = 0.0;
final Double earthDistance = distanceLOSB[1];
interMapping.addMapping(realPixelA, realPixelB, losDistance, earthDistance);
// Increment the number of measurements
this.measurementCount++;
} // end test if geoDistance < outlier
} // end test if sensorPixelB != null
} // end loop on pixel of sensorA
} // end loop on line of sensorA
this.observables.addInterMapping(interMapping);
}
/** Default constructor: measurements generation without outlier points control
* and Earth distances constraint.
* @param rA Rugged instance A
* @param sNameA sensor name A
* @param dimA dimension for acquisition A
* @param rB Rugged instance B
* @param sNameB sensor name B
* @param dimB dimension for acquisition B
*/
private void initParams(final Rugged rA, final String sNameA, final int dimA,
final Rugged rB, final String sNameB, final int dimB) {
this.sensorNameB = sNameB;
// Check that sensors's name is different
if (sNameA.contains(sNameB)) {
throw new RuggedException(RuggedMessages.DUPLICATED_PARAMETER_NAME, sNameA);
}
this.ruggedA = rA;
this.ruggedB = rB;
this.sensorA = rA.getLineSensor(sNameA);
this.sensorB = rB.getLineSensor(sNameB);
this.dimensionA = dimA;
this.dimensionB = dimB;
this.measurementCount = 0;
}
}
src/tutorials/java/fr/cs/examples/refiningPleiades/generators/Measurable.java 0 → 100644
View file @ acf902ea
/* Copyright 2013-2022 CS GROUP
* Licensed to CS GROUP (CS) under one or more
* contributor license agreements. See the NOTICE file distributed with
* this work for additional information regarding copyright ownership.
* CS licenses this file to You under the Apache License, Version 2.0
* (the "License"); you may not use this file except in compliance with
* the License. You may obtain a copy of the License at
*
* http://www.apache.org/licenses/LICENSE-2.0
*
* Unless required by applicable law or agreed to in writing, software
* distributed under the License is distributed on an "AS IS" BASIS,
* WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
* See the License for the specific language governing permissions and
* limitations under the License.
*/
package fr.cs.examples.refiningPleiades.generators;
/** For measurements generator.
* @author Lucie Labat-Allee
* @author Guylaine Prat
* @since 2.0
*/
public interface Measurable {
/** Get the number of measurements
* @return the number of measurements
*/
int getMeasurementCount();
/** Create measurements (without noise)
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
*/
void createMeasurement(int lineSampling, int pixelSampling);
/** Create noisy measurements
* @param lineSampling line sampling
* @param pixelSampling pixel sampling
* @param noise the noise to add to the measurements
*/
void createNoisyMeasurement(int lineSampling, int pixelSampling, Noise noise);
}