/////////////////////////////////////////////////////////////////////////////////// // Copyright (C) 2022 Jon Beniston, M7RCE // // // // 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 . // /////////////////////////////////////////////////////////////////////////////////// #include #include "czml.h" #include "mapsettings.h" #include "mapmodel.h" #include "util/units.h" CZML::CZML(const MapSettings *settings) : m_settings(settings) { } QJsonObject CZML::init() { QString start = QDateTime::currentDateTimeUtc().toString(Qt::ISODate); QString stop = QDateTime::currentDateTimeUtc().addSecs(60*60).toString(Qt::ISODate); QString interval = QString("%1/%2").arg(start).arg(stop); QJsonObject spec { {"interval", interval}, {"currentTime", start}, {"range", "UNBOUNDED"} }; QJsonObject doc { {"id", "document"}, {"version", "1.0"}, {"clock", spec} }; return doc; } // 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 CZML::cartesian3FromDegrees(double longitude, double latitude, double height) const { return cartesianFromRadians(Units::degreesToRadians(longitude), Units::degreesToRadians(latitude), height); } // FIXME: QVector3D is only float! // See Cesium.Cartesian3.fromRadians QVector3D CZML::cartesianFromRadians(double longitude, double latitude, double height) const { 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 CZML::fromHeadingPitchRoll(double heading, double pitch, double roll) const { 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 CZML::eastNorthUpToFixedFrame(QVector3D origin) const { // 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 CZML::fromRotation(QMatrix3x3 mat) const { 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 CZML::orientation(double longitude, double latitude, double altitude, double heading, double pitch, double roll) const { // 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 = cartesian3FromDegrees(longitude, latitude, altitude); // Calculate matrix to transform from East, North, Up (ENU) frame to ECEF frame QMatrix4x4 enuToECEFTransform = eastNorthUpToFixedFrame(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; } QJsonObject CZML::update(MapItem *mapItem, bool isTarget, bool isSelected) { (void) isTarget; // Don't currently use CLIP_TO_GROUND in Cesium due to Jitter bug // https://github.com/CesiumGS/cesium/issues/4049 // Instead we implement our own clipping code in map3d.html const QStringList heightReferences = {"NONE", "CLAMP_TO_GROUND", "RELATIVE_TO_GROUND", "NONE"}; QString dt; if (mapItem->m_takenTrackDateTimes.size() > 0) { dt = mapItem->m_takenTrackDateTimes.last()->toString(Qt::ISODateWithMs); } else { dt = QDateTime::currentDateTimeUtc().toString(Qt::ISODateWithMs); } QString id = mapItem->m_name; // Keep a hash of the time we first saw each item bool existingId = m_ids.contains(id); if (!existingId) { m_ids.insert(id, dt); } bool removeObj = false; bool fixedPosition = mapItem->m_fixedPosition; float displayDistanceMax = std::numeric_limits::max(); QString image = mapItem->m_image; if ((image == "antenna.png") || (image == "antennaam.png") || (image == "antennadab.png") || (image == "antennafm.png") || (image == "antennatime.png")) { displayDistanceMax = 1000000; } if (image == "") { // Need to remove this from the map removeObj = true; } QJsonArray coords; if (!removeObj) { if (!fixedPosition && (mapItem->m_predictedTrackCoords.size() > 0)) { QListIterator i(mapItem->m_takenTrackCoords); QListIterator j(mapItem->m_takenTrackDateTimes); while (i.hasNext()) { QGeoCoordinate *c = i.next(); coords.append(j.next()->toString(Qt::ISODateWithMs)); coords.append(c->longitude()); coords.append(c->latitude()); coords.append(c->altitude()); } if (mapItem->m_predictedTrackCoords.size() > 0) { QListIterator k(mapItem->m_predictedTrackCoords); QListIterator l(mapItem->m_predictedTrackDateTimes); k.toBack(); l.toBack(); while (k.hasPrevious()) { QGeoCoordinate *c = k.previous(); coords.append(l.previous()->toString(Qt::ISODateWithMs)); coords.append(c->longitude()); coords.append(c->latitude()); coords.append(c->altitude()); } } } else { // Only send latest position, to reduce processing if (!fixedPosition && mapItem->m_positionDateTime.isValid()) { coords.push_back(mapItem->m_positionDateTime.toString(Qt::ISODateWithMs)); } coords.push_back(mapItem->m_longitude); coords.push_back(mapItem->m_latitude); coords.push_back(mapItem->m_altitude); } } else { coords = m_lastPosition.value(id); } QJsonObject position { {"cartographicDegrees", coords}, }; if (!fixedPosition) { // Don't use forward extrapolation for satellites (with predicted tracks), as // it seems to jump about. We use it for AIS and ADS-B that don't have predicted tracks if (mapItem->m_predictedTrackCoords.size() == 0) { // Need 2 different positions to enable extrapolation, otherwise entity may not appear bool hasMoved = m_hasMoved.contains(id); if (!hasMoved && m_lastPosition.contains(id) && (m_lastPosition.value(id) != coords)) { hasMoved = true; m_hasMoved.insert(id, true); } if (hasMoved) { position.insert("forwardExtrapolationType", "EXTRAPOLATE"); position.insert("forwardExtrapolationDuration", 60); // Use linear interpolation for now - other two can go crazy with aircraft on the ground //position.insert("interpolationAlgorithm", "HERMITE"); //position.insert("interpolationDegree", "2"); //position.insert("interpolationAlgorithm", "LAGRANGE"); //position.insert("interpolationDegree", "5"); } else { position.insert("forwardExtrapolationType", "HOLD"); } } else { // Interpolation goes wrong at end points //position.insert("interpolationAlgorithm", "LAGRANGE"); //position.insert("interpolationDegree", "5"); //position.insert("interpolationAlgorithm", "HERMITE"); //position.insert("interpolationDegree", "2"); } } QQuaternion q = orientation(mapItem->m_longitude, mapItem->m_latitude, mapItem->m_altitude, mapItem->m_heading, mapItem->m_pitch, mapItem->m_roll); QJsonArray quaternion; if (!fixedPosition && mapItem->m_orientationDateTime.isValid()) { quaternion.push_back(mapItem->m_orientationDateTime.toString(Qt::ISODateWithMs)); } quaternion.push_back(q.x()); quaternion.push_back(q.y()); quaternion.push_back(q.z()); quaternion.push_back(q.scalar()); QJsonObject orientation { {"unitQuaternion", quaternion}, {"forwardExtrapolationType", "HOLD"}, // If we extrapolate, aircraft tend to spin around {"forwardExtrapolationDuration", 60}, // {"interpolationAlgorithm", "LAGRANGE"} }; QJsonObject orientationPosition { {"velocityReference", "#position"}, }; QJsonObject noPosition { {"cartographicDegrees", coords}, {"forwardExtrapolationType", "NONE"} }; // Point QColor pointColor = QColor::fromRgba(mapItem->m_itemSettings->m_3DPointColor); QJsonArray pointRGBA { pointColor.red(), pointColor.green(), pointColor.blue(), pointColor.alpha() }; QJsonObject pointColorObj { {"rgba", pointRGBA} }; QJsonObject point { {"pixelSize", 8}, {"color", pointColorObj}, {"heightReference", heightReferences[mapItem->m_altitudeReference]}, {"show", mapItem->m_itemSettings->m_enabled && mapItem->m_itemSettings->m_display3DPoint} }; // If clamping to ground, we need to disable depth test, so part of the point isn't clipped // However, when the point isn't clamped to ground, we shouldn't use this, otherwise // the point will become visible through the globe if (mapItem->m_altitudeReference == 1) { point.insert("disableDepthTestDistance", 100000000); } // Model QJsonArray node0Cartesian { {0.0, mapItem->m_modelAltitudeOffset, 0.0} }; QJsonObject node0Translation { {"cartesian", node0Cartesian} }; QJsonObject node0Transform { {"translation", node0Translation} }; QJsonObject nodeTransforms { {"node0", node0Transform}, }; QJsonObject model { {"gltf", m_settings->m_modelURL + mapItem->m_model}, {"incrementallyLoadTextures", false}, // Aircraft will flash as they appear without textures if this is the default of true {"heightReference", heightReferences[mapItem->m_altitudeReference]}, {"runAnimations", false}, {"show", mapItem->m_itemSettings->m_enabled && mapItem->m_itemSettings->m_display3DModel}, {"minimumPixelSize", mapItem->m_itemSettings->m_3DModelMinPixelSize}, {"maximumScale", 20000} // Stop it getting too big when zoomed really far out }; if (mapItem->m_modelAltitudeOffset != 0.0) { model.insert("nodeTransformations", nodeTransforms); } // Path QColor pathColor = QColor::fromRgba(mapItem->m_itemSettings->m_3DTrackColor); QJsonArray pathColorRGBA { pathColor.red(), pathColor.green(), pathColor.blue(), pathColor.alpha() }; QJsonObject pathColorObj { {"rgba", pathColorRGBA} }; // Paths can't be clamped to ground, so AIS paths can be underground if terrain is used // See: https://github.com/CesiumGS/cesium/issues/7133 QJsonObject pathSolidColorMaterial { {"color", pathColorObj} }; QJsonObject pathMaterial { {"solidColor", pathSolidColorMaterial} }; bool showPath = mapItem->m_itemSettings->m_enabled && mapItem->m_itemSettings->m_display3DTrack && ( m_settings->m_displayAllGroundTracks || (m_settings->m_displaySelectedGroundTracks && isSelected)); QJsonObject path { // We want full paths for sat tracker, so leadTime and trailTime should be 0 // Should be configurable.. 6000=100mins ~> 1 orbit for LEO //{"leadTime", "6000"}, //{"trailTime", "6000"}, {"width", "3"}, {"material", pathMaterial}, {"show", showPath} }; // Label QJsonArray labelPixelOffsetArray { 20, 0 }; QJsonObject labelPixelOffset { {"cartesian2", labelPixelOffsetArray} }; QJsonArray labelEyeOffsetArray { 0, mapItem->m_labelAltitudeOffset, 0 // Position above the object, dependent on the height of the model }; QJsonObject labelEyeOffset { {"cartesian", labelEyeOffsetArray} }; QJsonObject labelHorizontalOrigin { {"horizontalOrigin", "LEFT"} }; QJsonArray labelDisplayDistance { 0, displayDistanceMax }; QJsonObject labelDistanceDisplayCondition { {"distanceDisplayCondition", labelDisplayDistance} }; QJsonObject label { {"text", mapItem->m_label}, {"show", m_settings->m_displayNames && mapItem->m_itemSettings->m_enabled && mapItem->m_itemSettings->m_display3DLabel}, {"scale", mapItem->m_itemSettings->m_3DLabelScale}, {"pixelOffset", labelPixelOffset}, {"eyeOffset", labelEyeOffset}, {"verticalOrigin", "BASELINE"}, {"horizontalOrigin", "LEFT"}, {"heightReference", heightReferences[mapItem->m_altitudeReference]}, }; if (displayDistanceMax != std::numeric_limits::max()) { label.insert("disableDepthTestDistance", 100000000.0); label.insert("distanceDisplayCondition", labelDistanceDisplayCondition); } // Use billboard for APRS as we don't currently have 3D objects QString imageURL = mapItem->m_image; if (imageURL.startsWith("qrc://")) { imageURL = imageURL.mid(6); // Redirect to our embedded webserver, which will check resources } QJsonObject billboard { {"image", imageURL}, {"heightReference", heightReferences[mapItem->m_altitudeReference]}, {"verticalOrigin", "BOTTOM"} // To stop it being cut in half when zoomed out }; if (mapItem->m_altitudeReference == 1) { billboard.insert("disableDepthTestDistance", 100000000); } QJsonObject obj { {"id", id} // id must be unique }; if (!removeObj) { obj.insert("position", position); if (!fixedPosition) { if (mapItem->m_useHeadingPitchRoll) { obj.insert("orientation", orientation); } else { obj.insert("orientation", orientationPosition); } } obj.insert("point", point); if (!mapItem->m_model.isEmpty()) { obj.insert("model", model); } else { obj.insert("billboard", billboard); } obj.insert("label", label); obj.insert("description", mapItem->m_text); if (!fixedPosition) { obj.insert("path", path); } if (!fixedPosition) { if (mapItem->m_takenTrackDateTimes.size() > 0 && mapItem->m_predictedTrackDateTimes.size() > 0) { QString availability = QString("%1/%2") .arg(mapItem->m_takenTrackDateTimes.last()->toString(Qt::ISODateWithMs)) .arg(mapItem->m_predictedTrackDateTimes.last()->toString(Qt::ISODateWithMs)); obj.insert("availability", availability); } else { QString oneMin = QDateTime::currentDateTimeUtc().addSecs(60).toString(Qt::ISODateWithMs); QString createdToNow = QString("%1/%2").arg(m_ids[id]).arg(oneMin); // From when object was created to now obj.insert("availability", createdToNow); } } m_lastPosition.insert(id, coords); } else { // Disable forward extrapolation obj.insert("position", noPosition); } // Use our own clipping routine, due to // https://github.com/CesiumGS/cesium/issues/4049 if (mapItem->m_altitudeReference == 3) { obj.insert("altitudeReference", "CLIP_TO_GROUND"); } //qDebug() << obj; return obj; }