mirror of
https://github.com/f4exb/sdrangel.git
synced 2024-11-22 08:04:49 -05:00
319 lines
12 KiB
C++
319 lines
12 KiB
C++
///////////////////////////////////////////////////////////////////////////////////
|
|
// Copyright (C) 2022 Jon Beniston, M7RCE <jon@beniston.com> //
|
|
// Copyright (C) 2011-2020 Cesium Contributors //
|
|
// //
|
|
// This program is free software; you can redistribute it and/or modify //
|
|
// it under the terms of the GNU General Public License as published by //
|
|
// the Free Software Foundation as version 3 of the License, or //
|
|
// (at your option) any later version. //
|
|
// //
|
|
// This program is distributed in the hope that it will be useful, //
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
|
|
// GNU General Public License V3 for more details. //
|
|
// //
|
|
// You should have received a copy of the GNU General Public License //
|
|
// along with this program. If not, see <http://www.gnu.org/licenses/>. //
|
|
///////////////////////////////////////////////////////////////////////////////////
|
|
|
|
#include "coordinates.h"
|
|
#include "units.h"
|
|
|
|
// Scale cartesian position on to surface of ellipsoid
|
|
QVector3D Coordinates::scaleToGeodeticSurface(QVector3D cartesian, QVector3D oneOverRadii, QVector3D oneOverRadiiSquared)
|
|
{
|
|
float centerToleranceSquared = 0.1;
|
|
|
|
double x2 = cartesian.x() * cartesian.x() * oneOverRadii.x() * oneOverRadii.x();
|
|
double y2 = cartesian.y() * cartesian.y() * oneOverRadii.y() * oneOverRadii.y();
|
|
double z2 = cartesian.z() * cartesian.z() * oneOverRadii.z() * oneOverRadii.z();
|
|
|
|
double squaredNorm = x2 + y2 + z2;
|
|
double ratio = sqrt(1.0 / squaredNorm);
|
|
|
|
QVector3D intersection = cartesian * ratio;
|
|
|
|
if (squaredNorm < centerToleranceSquared) {
|
|
return intersection;
|
|
}
|
|
|
|
QVector3D gradient(
|
|
intersection.x() * oneOverRadiiSquared.x() * 2.0,
|
|
intersection.y() * oneOverRadiiSquared.y() * 2.0,
|
|
intersection.z() * oneOverRadiiSquared.z() * 2.0
|
|
);
|
|
|
|
double lambda = ((1.0 - ratio) * cartesian.length()) / (0.5 * gradient.length());
|
|
|
|
double correction = 0.0;
|
|
double func;
|
|
double denominator;
|
|
double xMultiplier;
|
|
double yMultiplier;
|
|
double zMultiplier;
|
|
double xMultiplier2;
|
|
double yMultiplier2;
|
|
double zMultiplier2;
|
|
double xMultiplier3;
|
|
double yMultiplier3;
|
|
double zMultiplier3;
|
|
|
|
do
|
|
{
|
|
lambda -= correction;
|
|
|
|
xMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquared.x());
|
|
yMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquared.y());
|
|
zMultiplier = 1.0 / (1.0 + lambda * oneOverRadiiSquared.z());
|
|
|
|
xMultiplier2 = xMultiplier * xMultiplier;
|
|
yMultiplier2 = yMultiplier * yMultiplier;
|
|
zMultiplier2 = zMultiplier * zMultiplier;
|
|
|
|
xMultiplier3 = xMultiplier2 * xMultiplier;
|
|
yMultiplier3 = yMultiplier2 * yMultiplier;
|
|
zMultiplier3 = zMultiplier2 * zMultiplier;
|
|
|
|
func = x2 * xMultiplier2 + y2 * yMultiplier2 + z2 * zMultiplier2 - 1.0;
|
|
|
|
denominator =
|
|
x2 * xMultiplier3 * oneOverRadiiSquared.x() +
|
|
y2 * yMultiplier3 * oneOverRadiiSquared.y() +
|
|
z2 * zMultiplier3 * oneOverRadiiSquared.z();
|
|
|
|
double derivative = -2.0 * denominator;
|
|
|
|
correction = func / derivative;
|
|
}
|
|
while (abs(func) > 0.000000000001);
|
|
|
|
QVector3D result(
|
|
cartesian.x() * xMultiplier,
|
|
cartesian.y() * yMultiplier,
|
|
cartesian.z() * zMultiplier
|
|
);
|
|
return result;
|
|
}
|
|
|
|
// QVector3D.normalized doesn't work with small numbers
|
|
QVector3D Coordinates::normalized(QVector3D vec)
|
|
{
|
|
QVector3D result;
|
|
float magnitude = vec.length();
|
|
result.setX(vec.x() / magnitude);
|
|
result.setY(vec.y() / magnitude);
|
|
result.setZ(vec.z() / magnitude);
|
|
return result;
|
|
}
|
|
|
|
// Convert ECEF position to geodetic coordinates
|
|
void Coordinates::ecefToGeodetic(double x, double y, double z, double &latitude, double &longitude, double &height)
|
|
{
|
|
QVector3D wgs84OneOverRadix(1.0 / 6378137.0,
|
|
1.0 / 6378137.0,
|
|
1.0 / 6356752.3142451793);
|
|
QVector3D wgs84OneOverRadiiSquared(1.0 / (6378137.0 * 6378137.0),
|
|
1.0 / (6378137.0 * 6378137.0),
|
|
1.0 / (6356752.3142451793 * 6356752.3142451793));
|
|
|
|
QVector3D cartesian(x, y, z);
|
|
|
|
QVector3D p = scaleToGeodeticSurface(cartesian, wgs84OneOverRadix, wgs84OneOverRadiiSquared);
|
|
|
|
QVector3D n = p * wgs84OneOverRadiiSquared;
|
|
n = normalized(n);
|
|
|
|
QVector3D h = cartesian - p;
|
|
|
|
longitude = atan2(n.y(), n.x());
|
|
latitude = asin(n.z());
|
|
|
|
longitude = Units::radiansToDegrees(longitude);
|
|
latitude = Units::radiansToDegrees(latitude);
|
|
|
|
double t = QVector3D::dotProduct(h, cartesian);
|
|
double sign = t >= 0.0 ? 1.0 : 0.0;
|
|
height = sign * h.length();
|
|
}
|
|
|
|
// Convert ECEF velocity to speed and heading
|
|
void Coordinates::ecefVelToSpeedHeading(double latitude, double longitude,
|
|
double velX, double velY, double velZ,
|
|
double &speed, double &verticalRate, double &heading)
|
|
{
|
|
if ((velX == 0.0) && (velY == 0.0) && (velZ == 0.0))
|
|
{
|
|
speed = 0.0;
|
|
heading = 0.0;
|
|
verticalRate = 0.0;
|
|
return;
|
|
}
|
|
|
|
double latRad = Units::degreesToRadians(latitude);
|
|
double lonRad = Units::degreesToRadians(longitude);
|
|
|
|
double sinLat = sin(latRad);
|
|
double cosLat = cos(latRad);
|
|
double sinLon = sin(lonRad);
|
|
double cosLon = cos(lonRad);
|
|
|
|
double velEast = -velX * sinLon + velY * cosLon;
|
|
double velNorth = -velX * sinLat * cosLon - velY * sinLat * sinLon + velZ * cosLat;
|
|
double velUp = velX * cosLat * cosLon + velY * cosLat * sinLon + velZ * sinLat;
|
|
|
|
speed = sqrt(velNorth * velNorth + velEast * velEast);
|
|
verticalRate = velUp;
|
|
|
|
double headingRad = atan2(velEast, velNorth);
|
|
heading = Units::radiansToDegrees(headingRad);
|
|
if (heading < 0.0) {
|
|
heading += 360.0;
|
|
} else if (heading >= 360.0) {
|
|
heading -= 360.0;
|
|
}
|
|
}
|
|
|
|
// Convert a position specified in longitude, latitude in degrees and height in metres above WGS84 ellipsoid in to
|
|
// Earth Centered Earth Fixed frame cartesian coordinates
|
|
// See Cesium.Cartesian3.fromDegrees
|
|
QVector3D Coordinates::geodeticToECEF(double longitude, double latitude, double height)
|
|
{
|
|
return geodeticRadiansToECEF(Units::degreesToRadians(longitude), Units::degreesToRadians(latitude), height);
|
|
}
|
|
|
|
// FIXME: QVector3D is only float!
|
|
// See Cesium.Cartesian3.fromRadians
|
|
QVector3D Coordinates::geodeticRadiansToECEF(double longitude, double latitude, double height)
|
|
{
|
|
QVector3D wgs84RadiiSquared(6378137.0 * 6378137.0, 6378137.0 * 6378137.0, 6356752.3142451793 * 6356752.3142451793);
|
|
|
|
double cosLatitude = cos(latitude);
|
|
QVector3D n;
|
|
n.setX(cosLatitude * cos(longitude));
|
|
n.setY(cosLatitude * sin(longitude));
|
|
n.setZ(sin(latitude));
|
|
n.normalize();
|
|
QVector3D k;
|
|
k = wgs84RadiiSquared * n;
|
|
double gamma = sqrt(QVector3D::dotProduct(n, k));
|
|
k = k / gamma;
|
|
n = n * height;
|
|
return k + n;
|
|
}
|
|
|
|
// Convert heading, pitch and roll in degrees to a quaternoin
|
|
// See: Cesium.Quaternion.fromHeadingPitchRoll
|
|
QQuaternion Coordinates::fromHeadingPitchRoll(double heading, double pitch, double roll)
|
|
{
|
|
QVector3D xAxis(1, 0, 0);
|
|
QVector3D yAxis(0, 1, 0);
|
|
QVector3D zAxis(0, 0, 1);
|
|
|
|
QQuaternion rollQ = QQuaternion::fromAxisAndAngle(xAxis, roll);
|
|
|
|
QQuaternion pitchQ = QQuaternion::fromAxisAndAngle(yAxis, -pitch);
|
|
|
|
QQuaternion headingQ = QQuaternion::fromAxisAndAngle(zAxis, -heading);
|
|
|
|
QQuaternion temp = rollQ * pitchQ;
|
|
|
|
return headingQ * temp;
|
|
}
|
|
|
|
// Calculate a transformation matrix from a East, North, Up frame at the given position to Earth Centered Earth Fixed frame
|
|
// See: Cesium.Transforms.eastNorthUpToFixedFrame
|
|
QMatrix4x4 Coordinates::eastNorthUpToECEF(QVector3D origin)
|
|
{
|
|
// TODO: Handle special case at centre of earth and poles
|
|
QVector3D up = origin.normalized();
|
|
QVector3D east(-origin.y(), origin.x(), 0.0);
|
|
east.normalize();
|
|
QVector3D north = QVector3D::crossProduct(up, east);
|
|
QMatrix4x4 result(
|
|
east.x(), north.x(), up.x(), origin.x(),
|
|
east.y(), north.y(), up.y(), origin.y(),
|
|
east.z(), north.z(), up.z(), origin.z(),
|
|
0.0, 0.0, 0.0, 1.0
|
|
);
|
|
return result;
|
|
}
|
|
|
|
// Convert 3x3 rotation matrix to a quaternoin
|
|
// Although there is a method for this in Qt: QQuaternion::fromRotationMatrix, it seems to
|
|
// result in different signs, so the following is based on Cesium code
|
|
QQuaternion Coordinates::fromRotation(QMatrix3x3 mat)
|
|
{
|
|
QQuaternion q;
|
|
|
|
double trace = mat(0, 0) + mat(1, 1) + mat(2, 2);
|
|
|
|
if (trace > 0.0)
|
|
{
|
|
double root = sqrt(trace + 1.0);
|
|
q.setScalar(0.5 * root);
|
|
root = 0.5 / root;
|
|
|
|
q.setX((mat(2,1) - mat(1,2)) * root);
|
|
q.setY((mat(0,2) - mat(2,0)) * root);
|
|
q.setZ((mat(1,0) - mat(0,1)) * root);
|
|
}
|
|
else
|
|
{
|
|
double next[] = {1, 2, 0};
|
|
int i = 0;
|
|
if (mat(1,1) > mat(0,0)) {
|
|
i = 1;
|
|
}
|
|
if (mat(2,2) > mat(0,0) && mat(2,2) > mat(1,1)) {
|
|
i = 2;
|
|
}
|
|
int j = next[i];
|
|
int k = next[j];
|
|
|
|
double root = sqrt(mat(i,i) - mat(j,j) - mat(k,k) + 1);
|
|
double quat[] = {0.0, 0.0, 0.0};
|
|
quat[i] = 0.5 * root;
|
|
root = 0.5 / root;
|
|
|
|
q.setScalar((mat(j,k) - mat(k,j)) * root);
|
|
quat[j] = (mat(i,j) + mat(j,i)) * root;
|
|
quat[k] = (mat(i,k) + mat(k,i)) * root;
|
|
q.setX(-quat[0]);
|
|
q.setY(-quat[1]);
|
|
q.setZ(-quat[2]);
|
|
}
|
|
return q;
|
|
}
|
|
|
|
// Calculate orientation quaternion for a model (such as an aircraft) based on position and (HPR) heading, pitch and roll (in degrees)
|
|
// While Cesium supports specifying orientation as HPR, CZML doesn't currently. See https://github.com/CesiumGS/cesium/issues/5184
|
|
// CZML requires the orientation to be in the Earth Centered Earth Fixed (geocentric) reference frame (https://en.wikipedia.org/wiki/Local_tangent_plane_coordinates)
|
|
// The orientation therefore depends not only on HPR but also on position
|
|
//
|
|
// glTF uses a right-handed axis convention; that is, the cross product of right and forward yields up. glTF defines +Y as up, +Z as forward, and -X as right.
|
|
// Cesium.Quaternion.fromHeadingPitchRoll Heading is the rotation about the negative z axis. Pitch is the rotation about the negative y axis. Roll is the rotation about the positive x axis.
|
|
QQuaternion Coordinates::orientation(double longitude, double latitude, double altitude, double heading, double pitch, double roll)
|
|
{
|
|
// Forward direction for gltf models in Cesium seems to be Eastward, rather than Northward, so we adjust heading by -90 degrees
|
|
heading = -90 + heading;
|
|
|
|
// Convert position to Earth Centered Earth Fixed (ECEF) frame
|
|
QVector3D positionECEF = geodeticToECEF(longitude, latitude, altitude);
|
|
|
|
// Calculate matrix to transform from East, North, Up (ENU) frame to ECEF frame
|
|
QMatrix4x4 enuToECEFTransform = eastNorthUpToECEF(positionECEF);
|
|
|
|
// Calculate rotation based on HPR in ENU frame
|
|
QQuaternion hprENU = fromHeadingPitchRoll(heading, pitch, roll);
|
|
|
|
// Transform rotation from ENU to ECEF
|
|
QMatrix3x3 hprENU3 = hprENU.toRotationMatrix();
|
|
QMatrix4x4 hprENU4(hprENU3);
|
|
QMatrix4x4 transform = enuToECEFTransform * hprENU4;
|
|
|
|
// Convert from 4x4 matrix to 3x3 matrix then to a quaternion
|
|
QQuaternion oq = fromRotation(transform.toGenericMatrix<3,3>());
|
|
|
|
return oq;
|
|
}
|