From 581eed4ac314c09533449ad64de9130e9b8cd772 Mon Sep 17 00:00:00 2001 From: Aang23 Date: Mon, 11 Oct 2021 23:17:17 +0200 Subject: [PATCH] Add better geodetic system, better projection algo and retune everything --- src-cli/project.cpp | 34 +- .../common/geodetic/euler_coordinates.cpp | 50 +++ src-core/common/geodetic/euler_coordinates.h | 27 ++ src-core/common/geodetic/euler_raytrace.cpp | 232 ++++++++++ src-core/common/geodetic/euler_raytrace.h | 15 + .../common/geodetic/geodetic_coordinates.cpp | 48 ++ .../common/geodetic/geodetic_coordinates.h | 27 ++ .../geodetic/projection/geo_projection.cpp | 72 +++ .../geodetic/projection/geo_projection.h | 42 ++ src-core/common/geodetic/projection/geos.cpp | 182 ++++++++ src-core/common/geodetic/projection/geos.h | 47 ++ .../geodetic/projection/leo_projection.cpp | 75 ++++ .../geodetic/projection/leo_projection.h | 72 +++ .../common/geodetic/projection/proj_file.cpp | 207 +++++++++ .../common/geodetic/projection/proj_file.h | 113 +++++ .../projection/satellite_reprojector.cpp | 180 ++++++++ .../projection/satellite_reprojector.h | 64 +++ .../common/geodetic/projection/stereo.cpp | 284 ++++++++++++ src-core/common/geodetic/projection/stereo.h | 58 +++ src-core/common/geodetic/projection/tpers.cpp | 238 ++++++++++ src-core/common/geodetic/projection/tpers.h | 60 +++ .../geodetic/vincentys_calculations.cpp | 284 ++++++++++++ .../common/geodetic/vincentys_calculations.h | 53 +++ src-core/common/geodetic/wgs84.h | 19 + src-core/common/projection/geo_projection.cpp | 69 --- src-core/common/projection/geo_projection.h | 39 -- src-core/common/projection/geos.cpp | 178 -------- src-core/common/projection/geos.h | 44 -- src-core/common/projection/leo_projection.cpp | 144 ------ src-core/common/projection/leo_projection.h | 79 ---- src-core/common/projection/proj_file.cpp | 220 --------- src-core/common/projection/proj_file.h | 114 ----- .../projection/satellite_reprojector.cpp | 176 -------- .../common/projection/satellite_reprojector.h | 55 --- src-core/common/projection/stereo.cpp | 280 ------------ src-core/common/projection/stereo.h | 55 --- src-core/common/projection/tpers.cpp | 234 ---------- src-core/common/projection/tpers.h | 57 --- .../instruments/modis/module_eos_modis.cpp | 47 +- .../instruments/erm/module_fengyun_erm.cpp | 28 +- .../instruments/mwhs/module_fengyun_mwhs.cpp | 28 +- .../mwhs2/module_fengyun_mwhs2.cpp | 28 +- .../instruments/mwts/module_fengyun_mwts.cpp | 28 +- .../mwts2/module_fengyun_mwts2.cpp | 28 +- .../mwts3/module_fengyun_mwts3.cpp | 28 +- .../instruments/virr/module_fengyun_virr.cpp | 32 +- .../instruments/atms/module_jpss_atms.cpp | 30 +- .../instruments/msumr/module_meteor_msumr.cpp | 48 +- .../instruments/amsu/module_metop_amsu.cpp | 54 +-- .../instruments/ascat/module_metop_ascat.cpp | 65 +-- .../instruments/avhrr/module_metop_avhrr.cpp | 28 +- .../instruments/mhs/module_metop_mhs.cpp | 28 +- .../instruments/avhrr/module_noaa_avhrr.cpp | 63 +-- src-interface/projection/projection.cpp | 425 +++++++++--------- 54 files changed, 2938 insertions(+), 2277 deletions(-) create mode 100644 src-core/common/geodetic/euler_coordinates.cpp create mode 100644 src-core/common/geodetic/euler_coordinates.h create mode 100644 src-core/common/geodetic/euler_raytrace.cpp create mode 100644 src-core/common/geodetic/euler_raytrace.h create mode 100644 src-core/common/geodetic/geodetic_coordinates.cpp create mode 100644 src-core/common/geodetic/geodetic_coordinates.h create mode 100644 src-core/common/geodetic/projection/geo_projection.cpp create mode 100644 src-core/common/geodetic/projection/geo_projection.h create mode 100644 src-core/common/geodetic/projection/geos.cpp create mode 100644 src-core/common/geodetic/projection/geos.h create mode 100644 src-core/common/geodetic/projection/leo_projection.cpp create mode 100644 src-core/common/geodetic/projection/leo_projection.h create mode 100644 src-core/common/geodetic/projection/proj_file.cpp create mode 100644 src-core/common/geodetic/projection/proj_file.h create mode 100644 src-core/common/geodetic/projection/satellite_reprojector.cpp create mode 100644 src-core/common/geodetic/projection/satellite_reprojector.h create mode 100644 src-core/common/geodetic/projection/stereo.cpp create mode 100644 src-core/common/geodetic/projection/stereo.h create mode 100644 src-core/common/geodetic/projection/tpers.cpp create mode 100644 src-core/common/geodetic/projection/tpers.h create mode 100644 src-core/common/geodetic/vincentys_calculations.cpp create mode 100644 src-core/common/geodetic/vincentys_calculations.h create mode 100644 src-core/common/geodetic/wgs84.h delete mode 100644 src-core/common/projection/geo_projection.cpp delete mode 100644 src-core/common/projection/geo_projection.h delete mode 100644 src-core/common/projection/geos.cpp delete mode 100644 src-core/common/projection/geos.h delete mode 100644 src-core/common/projection/leo_projection.cpp delete mode 100644 src-core/common/projection/leo_projection.h delete mode 100644 src-core/common/projection/proj_file.cpp delete mode 100644 src-core/common/projection/proj_file.h delete mode 100644 src-core/common/projection/satellite_reprojector.cpp delete mode 100644 src-core/common/projection/satellite_reprojector.h delete mode 100644 src-core/common/projection/stereo.cpp delete mode 100644 src-core/common/projection/stereo.h delete mode 100644 src-core/common/projection/tpers.cpp delete mode 100644 src-core/common/projection/tpers.h diff --git a/src-cli/project.cpp b/src-cli/project.cpp index 9031a6a44..0811da7d9 100644 --- a/src-cli/project.cpp +++ b/src-cli/project.cpp @@ -1,12 +1,12 @@ #include "project.h" #include "logger.h" -#include "common/projection/satellite_reprojector.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" #include "common/map/map_drawer.h" #include "resources.h" -#include "common/projection/stereo.h" -#include "common/projection/geos.h" +#include "common/geodetic/projection/stereo.h" +#include "common/geodetic/projection/geos.h" #include "common/map/maidenhead.h" int project(int argc, char *argv[]) @@ -34,8 +34,8 @@ int project(int argc, char *argv[]) return 1; } - projection::StereoProjection proj_stereo; - projection::GEOSProjection proj_geos(30000000, -0); + geodetic::projection::StereoProjection proj_stereo; + geodetic::projection::GEOSProjection proj_geos(30000000, -0); std::function(float, float, int, int)> projectionFunc = [](float lat, float lon, int map_height, int map_width) -> std::pair { @@ -101,29 +101,29 @@ int project(int argc, char *argv[]) { logger->info("Projecting " + image.first + "..."); cimg_library::CImg src_image(image.first.c_str()); - std::shared_ptr geofile = projection::proj_file::readReferenceFile(image.second); + std::shared_ptr geofile = geodetic::projection::proj_file::readReferenceFile(image.second); if (geofile->file_type == 0) { logger->info("Reprojecting Equiectangular..."); - projection::projectEQUIToproj(src_image, projected_image, src_image.spectrum(), projectionFunc); + geodetic::projection::projectEQUIToproj(src_image, projected_image, src_image.spectrum(), projectionFunc); } - else if (geofile->file_type == projection::proj_file::LEO_TYPE) + else if (geofile->file_type == geodetic::projection::proj_file::LEO_TYPE) { - projection::proj_file::LEO_GeodeticReferenceFile leofile = *((projection::proj_file::LEO_GeodeticReferenceFile *)geofile.get()); - projection::LEOScanProjectorSettings settings = leoProjectionRefFile(leofile); - projection::LEOScanProjector projector(settings); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile leofile = *((geodetic::projection::proj_file::LEO_GeodeticReferenceFile *)geofile.get()); + geodetic::projection::LEOScanProjectorSettings settings = leoProjectionRefFile(leofile); + geodetic::projection::LEOScanProjector projector(settings); logger->info("Reprojecting LEO..."); - projection::reprojectLEOtoProj(src_image, projector, projected_image, src_image.spectrum(), projectionFunc); + geodetic::projection::reprojectLEOtoProj(src_image, projector, projected_image, src_image.spectrum(), projectionFunc); } - else if (geofile->file_type == projection::proj_file::GEO_TYPE) + else if (geofile->file_type == geodetic::projection::proj_file::GEO_TYPE) { src_image.normalize(0, 65535); - projection::proj_file::GEO_GeodeticReferenceFile gsofile = *((projection::proj_file::GEO_GeodeticReferenceFile *)geofile.get()); + geodetic::projection::proj_file::GEO_GeodeticReferenceFile gsofile = *((geodetic::projection::proj_file::GEO_GeodeticReferenceFile *)geofile.get()); src_image.resize(gsofile.image_width, gsofile.image_height); // Safety logger->info("Reprojecting GEO..."); - projection::GEOProjector projector = projection::proj_file::geoProjectionRefFile(gsofile); - projection::reprojectGEOtoProj(src_image, projector, projected_image, src_image.spectrum(), projectionFunc); + geodetic::projection::GEOProjector projector = geodetic::projection::proj_file::geoProjectionRefFile(gsofile); + geodetic::projection::reprojectGEOtoProj(src_image, projector, projected_image, src_image.spectrum(), projectionFunc); } } diff --git a/src-core/common/geodetic/euler_coordinates.cpp b/src-core/common/geodetic/euler_coordinates.cpp new file mode 100644 index 000000000..5532a3b54 --- /dev/null +++ b/src-core/common/geodetic/euler_coordinates.cpp @@ -0,0 +1,50 @@ +#include "euler_coordinates.h" + +namespace geodetic +{ + euler_coords_t::euler_coords_t() + { + } + + euler_coords_t::euler_coords_t(double roll, double pitch, double yaw, bool radians) + { + this->roll = roll; + this->pitch = pitch; + this->yaw = yaw; + this->is_radians = radians; + } + + euler_coords_t euler_coords_t::toRads() // Convert to Rads if not done already + { + if (!is_radians) + { + roll *= DEG_TO_RAD; + pitch *= DEG_TO_RAD; + yaw *= DEG_TO_RAD; + is_radians = true; + } + + return *this; + } + + euler_coords_t euler_coords_t::toDegs() // Convert to Degs if not done already + { + if (is_radians) + { + roll *= RAD_TO_DEG; + pitch *= RAD_TO_DEG; + yaw *= RAD_TO_DEG; + is_radians = false; + } + + return *this; + } + + std::string euler_coords_t::str() + { + if (is_radians) + return "Roll: " + std::to_string(roll * RAD_TO_DEG) + ", Pitch: " + std::to_string(pitch * RAD_TO_DEG) + ", Yaw: " + std::to_string(yaw * RAD_TO_DEG); + else + return "Roll: " + std::to_string(roll) + ", Pitch: " + std::to_string(pitch) + ", Yaw: " + std::to_string(yaw); + } +}; diff --git a/src-core/common/geodetic/euler_coordinates.h b/src-core/common/geodetic/euler_coordinates.h new file mode 100644 index 000000000..f012ef1ee --- /dev/null +++ b/src-core/common/geodetic/euler_coordinates.h @@ -0,0 +1,27 @@ +#pragma once + +#include +#include +#define DEG_TO_RAD (M_PI / 180.0) +#define RAD_TO_DEG (180.0 / M_PI) + +namespace geodetic +{ + // Simple class to hold Geodetic coordinates and easily convert between degs and rads + class euler_coords_t + { + private: + bool is_radians = false; + + public: + double roll = 0; // Roll, degrees by default + double pitch = 0; // Pitch, degrees by default + double yaw = 0; // Yaw, degrees by default + + euler_coords_t(); + euler_coords_t(double roll, double pitch, double yaw, bool radians = false); + euler_coords_t toRads(); // Convert to Rads if not done already + euler_coords_t toDegs(); // Convert to Degs if not done already + std::string str(); + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/euler_raytrace.cpp b/src-core/common/geodetic/euler_raytrace.cpp new file mode 100644 index 000000000..ca40b44da --- /dev/null +++ b/src-core/common/geodetic/euler_raytrace.cpp @@ -0,0 +1,232 @@ +#include "euler_raytrace.h" +#include "wgs84.h" + +#define EPSILON 2.2204460492503131e-016 + +namespace geodetic +{ + /* + Namespace to store some stuff used only for this function, + at least for now. + I initially wanted to avoid implementing all this, but it + would have made most of the code much harder and shouldn't + really impact performances that badly. + */ + namespace raytrace_to_earth_namespace + { + struct vector + { + double x = 0; + double y = 0; + double z = 0; + + vector() {} + vector(double x, double y, double z) : x(x), y(y), z(z) {} + + double distance_square() + { + return x * x + y * y + z * z; + } + + double magnitude() + { + return sqrt(x * x + y * y + z * z); + } + + vector &normalize() + { + double mag = magnitude(); + + if (mag <= 0) + return *this; + + x /= mag; + y /= mag; + z /= mag; + + return *this; + } + }; + + vector vectorCross(vector a, vector b) + { + vector c; + c.x = a.y * b.z - a.z * b.y; + c.y = a.z * b.x - a.x * b.z; + c.z = a.x * b.y - a.y * b.x; + return c; + } + + struct matrix + { + double m[16]; + + matrix() + { + std::fill(m, &m[16], 0); + m[0] = 1; + m[5] = 1; + m[10] = 1; + m[15] = 1; + } + + matrix(double a1, double a2, double a3, double a4, + double a5, double a6, double a7, double a8, + double a9, double a10, double a11, double a12, + double a13, double a14, double a15, double a16) + { + m[0] = a1; + m[1] = a2; + m[2] = a3; + m[3] = a4; + m[4] = a5; + m[5] = a6; + m[6] = a7; + m[7] = a8; + m[8] = a9; + m[9] = a10; + m[10] = a11; + m[11] = a12; + m[12] = a13; + m[13] = a14; + m[14] = a15; + m[15] = a16; + } + + double &operator[](const int &i) + { + return m[i]; + } + + matrix operator*(const matrix &b) + { + matrix out; + out[0] = m[0] * b.m[0] + m[1] * b.m[4] + m[2] * b.m[8] + m[3] * b.m[12]; + out[1] = m[0] * b.m[1] + m[1] * b.m[5] + m[2] * b.m[9] + m[3] * b.m[13]; + out[2] = m[0] * b.m[2] + m[1] * b.m[6] + m[2] * b.m[10] + m[3] * b.m[14]; + out[3] = m[0] * b.m[3] + m[1] * b.m[7] + m[2] * b.m[11] + m[3] * b.m[15]; + out[4] = m[4] * b.m[0] + m[5] * b.m[4] + m[6] * b.m[8] + m[7] * b.m[12]; + out[5] = m[4] * b.m[1] + m[5] * b.m[5] + m[6] * b.m[9] + m[7] * b.m[13]; + out[6] = m[4] * b.m[2] + m[5] * b.m[6] + m[6] * b.m[10] + m[7] * b.m[14]; + out[7] = m[4] * b.m[3] + m[5] * b.m[7] + m[6] * b.m[11] + m[7] * b.m[15]; + out[8] = m[8] * b.m[0] + m[9] * b.m[4] + m[10] * b.m[8] + m[11] * b.m[12]; + out[9] = m[8] * b.m[1] + m[9] * b.m[5] + m[10] * b.m[9] + m[11] * b.m[13]; + out[10] = m[8] * b.m[2] + m[9] * b.m[6] + m[10] * b.m[10] + m[11] * b.m[14]; + out[11] = m[8] * b.m[3] + m[9] * b.m[7] + m[10] * b.m[11] + m[11] * b.m[15]; + out[12] = m[12] * b.m[0] + m[13] * b.m[4] + m[14] * b.m[8] + m[15] * b.m[12]; + out[13] = m[12] * b.m[1] + m[13] * b.m[5] + m[14] * b.m[9] + m[15] * b.m[13]; + out[14] = m[12] * b.m[2] + m[13] * b.m[6] + m[14] * b.m[10] + m[15] * b.m[14]; + out[15] = m[12] * b.m[3] + m[13] * b.m[7] + m[14] * b.m[11] + m[15] * b.m[15]; + return out; + } + }; + }; + + /* + I initially wrote this function in a much, much less elegant way... + But during my research for some examples of cleaner ways, I stumbled + upon https://github.com/Digitelektro/MeteorDemod and reused some + implementations. They were also simplified / cleaned up. + */ + int raytrace_to_earth(geodetic_coords_t position_geo, euler_coords_t pointing, geodetic_coords_t &earth_point) + { + // Ensure all inputs are in radians + position_geo.toRads(); + pointing.toRads(); + + // Generate rotation matrices + raytrace_to_earth_namespace::matrix rotateX; + rotateX[5] = cos(-pointing.roll); + rotateX[6] = -sin(-pointing.roll); + rotateX[9] = sin(-pointing.roll); + rotateX[10] = cos(-pointing.roll); + + raytrace_to_earth_namespace::matrix rotateY; + rotateY[0] = cos(-pointing.pitch); + rotateY[2] = sin(-pointing.pitch); + rotateY[8] = -sin(-pointing.pitch); + rotateY[10] = cos(-pointing.pitch); + + raytrace_to_earth_namespace::matrix rotateZ; + rotateZ[0] = cos(-pointing.yaw); + rotateZ[1] = -sin(-pointing.yaw); + rotateZ[4] = sin(-pointing.yaw); + rotateZ[5] = cos(-pointing.yaw); + + // Total rotation matrice + raytrace_to_earth_namespace::matrix rotateXYZ = rotateZ * rotateY * rotateX; + + // Geodetic coordinates to vector + raytrace_to_earth_namespace::vector position; + + double N = pow(WGS84::a, 2) / sqrt(pow(WGS84::a, 2) * pow(cos(position_geo.lat), 2) + pow(WGS84::b, 2) * pow(sin(position_geo.lat), 2)); + position.x = (N + position_geo.alt) * cos(position_geo.lat) * cos(position_geo.lon); + position.y = (N + position_geo.alt) * cos(position_geo.lat) * sin(position_geo.lon); + position.z = ((pow(WGS84::b, 2) / pow(WGS84::a, 2)) * N + position_geo.alt) * sin(position_geo.lat); + + // Matrices + raytrace_to_earth_namespace::matrix look_matrix; + raytrace_to_earth_namespace::matrix final_matrix(1, 0, 0, position.x, + 0, 1, 0, position.y, + 0, 0, 1, position.z, + 0, 0, 0, 1); + + raytrace_to_earth_namespace::vector k(-position.x, -position.y, -position.z); + double m = k.distance_square(); + if (m >= EPSILON) + { + k.x *= (1.0 / sqrt(m)); + k.y *= (1.0 / sqrt(m)); + k.z *= (1.0 / sqrt(m)); + + raytrace_to_earth_namespace::vector i = raytrace_to_earth_namespace::vectorCross(raytrace_to_earth_namespace::vector(0, 0, 1), k).normalize(); + raytrace_to_earth_namespace::vector j = raytrace_to_earth_namespace::vectorCross(k, i).normalize(); + + look_matrix = raytrace_to_earth_namespace::matrix(i.x, j.x, k.x, 0.0, + i.y, j.y, k.y, 0.0, + i.z, j.z, k.z, 0.0, + 0.0, 0.0, 0.0, 1.0); + } + + final_matrix = final_matrix * look_matrix * rotateXYZ; + + // Vector + double u = final_matrix[2]; + double v = final_matrix[6]; + double w = final_matrix[10]; + + // WGS84 ellipsoid + const double &a = WGS84::a; + const double &b = WGS84::a; + const double &c = WGS84::b; + + double value = -pow(a, 2) * pow(b, 2) * w * position.z - pow(a, 2) * pow(c, 2) * v * position.y - pow(b, 2) * pow(c, 2) * u * position.x; + double radical = pow(a, 2) * pow(b, 2) * pow(w, 2) + pow(a, 2) * pow(c, 2) * pow(v, 2) - pow(a, 2) * pow(v, 2) * pow(position.z, 2) + + 2 * pow(a, 2) * v * w * position.y * position.z - pow(a, 2) * pow(w, 2) * pow(position.y, 2) + pow(b, 2) * pow(c, 2) * pow(u, 2) - + pow(b, 2) * pow(u, 2) * pow(position.z, 2) + 2 * pow(b, 2) * u * w * position.x * position.z - pow(b, 2) * pow(w, 2) * pow(position.x, 2) - + pow(c, 2) * pow(u, 2) * pow(position.y, 2) + 2 * pow(c, 2) * u * v * position.x * position.y - pow(c, 2) * pow(v, 2) * pow(position.x, 2); + double magnitude = pow(a, 2) * pow(b, 2) * pow(w, 2) + pow(a, 2) * pow(c, 2) * pow(v, 2) + pow(b, 2) * pow(c, 2) * pow(u, 2); + + if (radical < 0) + return 1; + + double d = (value - a * b * c * sqrt(radical)) / magnitude; + + if (d < 0) + return 1; + + position.x += d * u; + position.y += d * v; + position.z += d * w; + + // To geodetic + double p = sqrt(pow(position.x, 2) + pow(position.y, 2)); + double phi = atan2(position.z * WGS84::a, p * WGS84::b); + double lat = atan2(position.z + pow(WGS84::e2, 2) * WGS84::b * pow(sin(phi), 3), p - pow(WGS84::e, 2) * WGS84::a * pow(cos(phi), 3)); + double lon = atan2(position.y, position.x); + + earth_point = geodetic::geodetic_coords_t(lat, lon, 0, true); + + return 0; + } +}; diff --git a/src-core/common/geodetic/euler_raytrace.h b/src-core/common/geodetic/euler_raytrace.h new file mode 100644 index 000000000..2541c0c12 --- /dev/null +++ b/src-core/common/geodetic/euler_raytrace.h @@ -0,0 +1,15 @@ +#pragma once + +#include "euler_coordinates.h" +#include "geodetic_coordinates.h" + +namespace geodetic +{ + /* + Given a geodetic position and orientation in Euler coordinates, + "trace" a way down to the WGS84 ellipsoid's surface and return + the geodetic position. + Returns 1 on error. + */ + int raytrace_to_earth(geodetic_coords_t position, euler_coords_t pointing, geodetic_coords_t &earth_point); +}; \ No newline at end of file diff --git a/src-core/common/geodetic/geodetic_coordinates.cpp b/src-core/common/geodetic/geodetic_coordinates.cpp new file mode 100644 index 000000000..ffa1df62f --- /dev/null +++ b/src-core/common/geodetic/geodetic_coordinates.cpp @@ -0,0 +1,48 @@ +#include "geodetic_coordinates.h" + +namespace geodetic +{ + geodetic_coords_t::geodetic_coords_t() + { + } + + geodetic_coords_t::geodetic_coords_t(double lat, double lon, double alt, bool radians) + { + this->lat = lat; + this->lon = lon; + this->alt = alt; + this->is_radians = radians; + } + + geodetic_coords_t geodetic_coords_t::toRads() // Convert to Rads if not done already + { + if (!is_radians) + { + lat *= DEG_TO_RAD; + lon *= DEG_TO_RAD; + is_radians = true; + } + + return *this; + } + + geodetic_coords_t geodetic_coords_t::toDegs() // Convert to Degs if not done already + { + if (is_radians) + { + lat *= RAD_TO_DEG; + lon *= RAD_TO_DEG; + is_radians = false; + } + + return *this; + } + + std::string geodetic_coords_t::str() + { + if (is_radians) + return "Lat: " + std::to_string(lat * RAD_TO_DEG) + ", Lon: " + std::to_string(lon * RAD_TO_DEG) + ", Alt: " + std::to_string(alt); + else + return "Lat: " + std::to_string(lat) + ", Lon: " + std::to_string(lon) + ", Alt: " + std::to_string(alt); + } +}; diff --git a/src-core/common/geodetic/geodetic_coordinates.h b/src-core/common/geodetic/geodetic_coordinates.h new file mode 100644 index 000000000..5cd5762ea --- /dev/null +++ b/src-core/common/geodetic/geodetic_coordinates.h @@ -0,0 +1,27 @@ +#pragma once + +#include +#include +#define DEG_TO_RAD (M_PI / 180.0) +#define RAD_TO_DEG (180.0 / M_PI) + +namespace geodetic +{ + // Simple class to hold Geodetic coordinates and easily convert between degs and rads + class geodetic_coords_t + { + private: + bool is_radians = false; + + public: + double lat = 0; // Latitude, degrees by default + double lon = 0; // Longitude, degrees by default + double alt = 0; // Altitude, kilometers + + geodetic_coords_t(); + geodetic_coords_t(double lat, double lon, double alt, bool radians = false); + geodetic_coords_t toRads(); // Convert to Rads if not done already + geodetic_coords_t toDegs(); // Convert to Degs if not done already + std::string str(); + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/geo_projection.cpp b/src-core/common/geodetic/projection/geo_projection.cpp new file mode 100644 index 000000000..368b33f5d --- /dev/null +++ b/src-core/common/geodetic/projection/geo_projection.cpp @@ -0,0 +1,72 @@ +#include "geo_projection.h" + +#include +#include "logger.h" + +namespace geodetic +{ + namespace projection + { + GEOProjector::GEOProjector(double sat_lon, + double sat_height, + int img_width, + int img_height, + double hscale, + double vscale, + double x_offset, + double y_offset, + bool sweep_x) : hscale(hscale), + vscale(vscale), + x_offset(x_offset), + y_offset(y_offset) + { + height = img_height; + width = img_width; + pj.init(sat_height * 1000, sat_lon, sweep_x); + } + + int GEOProjector::forward(double lon, double lat, int &img_x, int &img_y) + { + if (pj.forward(lon, lat, x, y)) + { + // Error / out of the image + img_x = -1; + img_y = -1; + return 1; + } + + image_x = x * hscale * (width / 2.0); + image_y = y * vscale * (height / 2.0); + + image_x += width / 2.0 + x_offset; + image_y += height / 2.0 + y_offset; + + img_x = image_x; + img_y = (height - 1) - image_y; + + return 0; + } + + int GEOProjector::inverse(int img_x, int img_y, double &lon, double &lat) + { + image_y = (height - 1) - img_y; + image_x = img_x; + + image_y -= height / 2.0 + y_offset; + image_x -= width / 2.0 + x_offset; + + y = image_y / (vscale * (height / 2.0)); + x = image_x / (hscale * (width / 2.0)); + + if (pj.inverse(x, y, lon, lat)) + { + // Error / out of the image + img_x = -1; + img_y = -1; + return 1; + } + + return 0; + } + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/geo_projection.h b/src-core/common/geodetic/projection/geo_projection.h new file mode 100644 index 000000000..f8fc9f5c3 --- /dev/null +++ b/src-core/common/geodetic/projection/geo_projection.h @@ -0,0 +1,42 @@ +#pragma once + +#include "geos.h" + +/* +Code to reference a decoded image (or similar data) from a GEO satellite to Lat / Lon coordinates. + +The inverse function is currently broken, for an unknwon reason. +*/ +namespace geodetic +{ + namespace projection + { + class GEOProjector + { + private: + projection::GEOSProjection pj; + double height, width; + + double x, y; + double image_x, image_y; + + double hscale; + double vscale; + double x_offset; + double y_offset; + + public: + GEOProjector(double sat_lon, + double sat_height, + int img_width, + int img_height, + double hscale, + double vscale, + double x_offset, + double y_offset, + bool sweep_x); + int forward(double lon, double lat, int &img_x, int &img_y); + int inverse(int img_x, int img_y, double &lon, double &lat); + }; + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/geos.cpp b/src-core/common/geodetic/projection/geos.cpp new file mode 100644 index 000000000..c173a3942 --- /dev/null +++ b/src-core/common/geodetic/projection/geos.cpp @@ -0,0 +1,182 @@ +#include "geos.h" +#include +#include "../wgs84.h" + +/* +** This file was adapted and simplified from libproj, and the below +** notice kept as credits. +** +** libproj -- library of cartographic projections +** +** Copyright (c) 2004 Gerald I. Evenden +** Copyright (c) 2012 Martin Raspaud +** +** See also (section 4.4.3.2): +** https://www.cgms-info.org/documents/pdf_cgms_03.pdf +** +** Permission is hereby granted, free of charge, to any person obtaining +** a copy of this software and associated documentation files (the +** "Software"), to deal in the Software without restriction, including +** without limitation the rights to use, copy, modify, merge, publish, +** distribute, sublicense, and/or sell copies of the Software, and to +** permit persons to whom the Software is furnished to do so, subject to +** the following conditions: +** +** The above copyright notice and this permission notice shall be +** included in all copies or substantial portions of the Software. +** +** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +*/ + +namespace geodetic +{ + namespace projection + { + int GEOSProjection::init(double height, double longitude, bool sweep_x) + { + lon_0 = longitude; // The satellite's longitude + + // Constants, WGS84 + phi0 = 0; + a = WGS84::a * 1000; + es = WGS84::es; + one_es = WGS84::one_es; + + // Orbit Height + h = height; + + // Scan axis + flip_axis = sweep_x; + + radius_g_1 = h / a; + if (radius_g_1 <= 0 || radius_g_1 > 1e10) + { + // Illegal case. + // Kept just in case but we shouldn't end up there unless the user makes a mistake... + return 1; + } + + // Init the rest + radius_g = 1. + radius_g_1; + C = radius_g * radius_g - 1.0; + + radius_p = sqrt(one_es); + radius_p2 = one_es; + radius_p_inv2 = one_es; + + return 0; + } + + int GEOSProjection::forward(double lon, double lat, double &x, double &y) + { + x = y = 0; // Safety + + // Shift longitudes to use the sat's as a reference + lon -= lon_0; + if (lon < -180) + lon = lon + 360; + if (lon > 180) + lon = lon - 360; + + // To radians + double phi = lat * 0.01745329, lam = lon * 0.01745329; + + double r, Vx, Vy, Vz, tmp; + + // Calculation of geocentric latitude. + phi = atan(radius_p2 * tan(phi)); + + // Calculation of the three components of the vector from satellite to position on earth surface (lon,lat). + r = (radius_p) / hypot(radius_p * cos(phi), sin(phi)); + Vx = r * cos(lam) * cos(phi); + Vy = r * sin(lam) * cos(phi); + Vz = r * sin(phi); + + // Check visibility. + if (((radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * radius_p_inv2) < 0.) + { + x = y = 2e10; // Trigger error + return 1; + } + + // Calculation based on view angles from satellite. + tmp = radius_g - Vx; + + if (flip_axis) + { + x = radius_g_1 * atan(Vy / hypot(Vz, tmp)); + y = radius_g_1 * atan(Vz / tmp); + } + else + { + x = radius_g_1 * atan(Vy / tmp); + y = radius_g_1 * atan(Vz / hypot(Vy, tmp)); + } + + return 0; + } + + int GEOSProjection::inverse(double x, double y, double &lon, double &lat) + { + lon = lat = 0.0; + double phi = 0, lam = 0; + + double Vx, Vy, Vz, a, b, k; + + // Setting three components of vector from satellite to position. + Vx = -1.0; + + if (flip_axis) + { + Vz = tan(y / radius_g_1); + Vy = tan(x / radius_g_1) * hypot(1.0, Vz); + } + else + { + Vy = tan(x / radius_g_1); + Vz = tan(y / radius_g_1) * hypot(1.0, Vy); + } + + // Calculation of terms in cubic equation and determinant. + a = Vz / radius_p; + a = Vy * Vy + a * a + Vx * Vx; + b = 2 * radius_g * Vx; + const double det = (b * b) - 4 * a * C; + if (det < 0.0) + { + lon = lat = 2e10; // Trigger error + return 1; + } + + // Calculation of three components of vector from satellite to position. + k = (-b - sqrt(det)) / (2. * a); + Vx = radius_g + k * Vx; + Vy *= k; + Vz *= k; + + // Calculation of longitude and latitude. + lam = atan2(Vy, Vx); + phi = atan(Vz * cos(lam) / Vx); + phi = atan(radius_p_inv2 * tan(phi)); + + // To degs + lat = phi * 57.29578; + lon = lam * 57.29578; + + // Shift longitudes back to reference 0 + lon += lon_0; + if (lon < -180) + lon = lon + 360; + if (lon > 180) + lon = lon - 360; + + return 0; + } + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/geos.h b/src-core/common/geodetic/projection/geos.h new file mode 100644 index 000000000..54fb57692 --- /dev/null +++ b/src-core/common/geodetic/projection/geos.h @@ -0,0 +1,47 @@ +#pragma once + +/* +Implementation of a standard GEOS projection, adapted from libproj. +Some variables are hardcoded for the intended usecase, making some +degree of tuning unecessary. +Uses the WGS84 ellipsoid. +*/ +namespace geodetic +{ + namespace projection + { + class GEOSProjection + { + private: + double h; + double radius_p; + double radius_p2; + double radius_p_inv2; + double radius_g; + double radius_g_1; + double C; + int flip_axis; + + double phi0; + double a; + double es; + double one_es; + + double lon_0; + + public: + GEOSProjection() + { + init(0, 0); + } + GEOSProjection(double height, double longitude, bool sweep_x = false) + { + init(height, longitude, sweep_x); + } + + int init(double height, double longitude, bool sweep_x = false); // return value of 1 => Error + int forward(double lon, double lat, double &x, double &y); // return value of 1 => Error + int inverse(double x, double y, double &lon, double &lat); // return value of 1 => Error + }; + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/leo_projection.cpp b/src-core/common/geodetic/projection/leo_projection.cpp new file mode 100644 index 000000000..2b56c8ca3 --- /dev/null +++ b/src-core/common/geodetic/projection/leo_projection.cpp @@ -0,0 +1,75 @@ +#include "leo_projection.h" + +#include +#include "libs/predict/predict.h" +#include "logger.h" +#include "common/geodetic/euler_raytrace.h" +#include "common/geodetic/vincentys_calculations.h" + +namespace geodetic +{ + namespace projection + { + void LEOScanProjector::generateOrbit() + { + // Setup SGP4 model + predict_orbital_elements_t *satellite_object = predict_parse_tle(settings.sat_tle.line1.c_str(), + settings.sat_tle.line2.c_str()); + predict_position satellite_orbit; + predict_position satellite_pos2; + + for (int currentScan = 0; currentScan < (int)settings.utc_timestamps.size(); currentScan++) + { + double currentTimestamp = settings.utc_timestamps[currentScan] + settings.time_offset; + + // Get Julian time of the scan, with full accuracy and calculate the satellite's + // position at the time + predict_julian_date_t currentJulianTime = predict_to_julian_double(currentTimestamp); + predict_orbit(satellite_object, &satellite_orbit, currentJulianTime); + satellite_positions.push_back(geodetic::geodetic_coords_t(satellite_orbit.latitude, satellite_orbit.longitude, satellite_orbit.altitude, true)); + + // Calculate Az to use from the satellite's motion vector + predict_orbit(satellite_object, &satellite_pos2, predict_to_julian_double(currentTimestamp - 1)); + + double az_angle_vincentis = vincentys_inverse(geodetic_coords_t(satellite_pos2.latitude, satellite_pos2.longitude, satellite_pos2.altitude, true), + geodetic_coords_t(satellite_orbit.latitude, satellite_orbit.longitude, satellite_orbit.altitude, true)) + .reverse_azimuth; + az_angle_vincentis -= M_PI; + if (az_angle_vincentis >= M_2_PI) + az_angle_vincentis -= M_2_PI; + + satellite_directions.push_back(az_angle_vincentis * RAD_TO_DEG); + satellite_is_asc.push_back(satellite_pos2.latitude < satellite_orbit.latitude); + } + } + + LEOScanProjector::LEOScanProjector(LEOScanProjectorSettings settings) : settings(settings) + { + logger->info("Compute orbit..."); + generateOrbit(); + } + + int LEOScanProjector::inverse(int img_x, int img_y, geodetic_coords_t &ground_position) + { + // Check we're in bounds + if (img_y > (int)satellite_positions.size() || img_x >= settings.image_width) + return 1; + + double final_x = settings.invert_scan ? (settings.image_width - 1) - img_x : img_x; + + geodetic::geodetic_coords_t &satellite_position = satellite_positions[img_y]; + + const bool &ascending = satellite_is_asc[img_y]; // Not a fan of this, but it works OK and doesn't seem too bad so well, whatever I guess? + + geodetic::euler_coords_t satellite_pointing; + satellite_pointing.roll = -(((final_x - (settings.image_width / 2)) / settings.image_width) * settings.scan_angle) + settings.roll_offset; + satellite_pointing.pitch = settings.pitch_offset; + satellite_pointing.yaw = (90 + (ascending ? settings.yaw_offset : -settings.yaw_offset)) - satellite_directions[img_y]; + + geodetic::raytrace_to_earth(satellite_position, satellite_pointing, ground_position); + ground_position.toDegs(); + + return 0; + } + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/leo_projection.h b/src-core/common/geodetic/projection/leo_projection.h new file mode 100644 index 000000000..6f59bdee6 --- /dev/null +++ b/src-core/common/geodetic/projection/leo_projection.h @@ -0,0 +1,72 @@ +#pragma once + +#include +#include "tle.h" +#include "common/geodetic/geodetic_coordinates.h" +#include "libs/predict/predict.h" + +/* +Code to reference a decoded image (or similar data) from a LEO satellite to Lat / Lon coordinates. +It works by calculating the satellite's position at the time of each scan, and generating a projection +for each of them. That projection references any point visible from the satellite's position to a +known Lat / Lon. +Doing this for every scan, with some curvature correction to match the image to the projection allows +referencing any given point. +A look-up-table is generated in the constructor to speed up later processing. + +You may also notice there is no forward function (eg, Lat / Lon to x/y on the image). That's because +the easiest way to do it would be getting the closest point in a generated LUT, but doing this +efficiently gets complicated quickly... + +PS : I am not sure in any way that this is a good implementation, there may be very obvious mistakes in there... + +PS #2 : Currently the projection is aligned to the satellite's vector by computing previous and later positions +of the satellite. It's approximate and may cause issues later but I have not found a proper way to get velocity +vectors relative to the satellite itself from the ECI values. At least, not in a form I can use for this purpose. + +PS #3 : Currently this only handles a single timestamp per scanline, which is usually most instruments work. However, +some don't like MODIS, IASI, etc. They instead send a timestamp per "IFOV" (Group of samples), such as 64x64 at once +for IASI. Others like MERSI send a full scanline of 40 pixels at once... It would be good to handle those cases properly, +utilizing each individual timestamp and properly simulating the group of detector etc... +I guess this calls for writing variants of this code later on. +*/ +namespace geodetic +{ + namespace projection + { + struct LEOScanProjectorSettings + { + double scan_angle; // Total scan angle + double roll_offset; // Roll offset + double pitch_offset; // Pitch offset + double yaw_offset; // Yaw offset + double time_offset; // Timestamp offset relative to the provided timestamps + int image_width; // Input image width + bool invert_scan; // Invert the scan direction relative to the projection + tle::TLE sat_tle; // Satellite TLEs + std::vector utc_timestamps; // Timestamps. Must match each scanline of the image you will be working with + }; + + class LEOScanProjector + { + private: + // Settings + const LEOScanProjectorSettings settings; + + // Luts and values used for referencing each line + std::vector satellite_positions; + std::vector satellite_directions; + std::vector satellite_is_asc; + + // Internal functions + void generateOrbit(); + + public: + std::vector poss; + LEOScanProjector(LEOScanProjectorSettings settings); + + int inverse(int img_x, int img_y, geodetic_coords_t &coords); // Transform image coordinates to lat / lon. Return 1 if there was an error + //void inverse(int img_x, int img_y, double &lat, double &lon); // Transform image coordinates to lat / lon. Calls up a LUT to be faster + }; + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/proj_file.cpp b/src-core/common/geodetic/projection/proj_file.cpp new file mode 100644 index 000000000..0b54a49c9 --- /dev/null +++ b/src-core/common/geodetic/projection/proj_file.cpp @@ -0,0 +1,207 @@ +#include "proj_file.h" +#include +#include "common/utils.h" +#include "tle.h" +#include +#include +#include "logger.h" + +namespace geodetic +{ + namespace projection + { + namespace proj_file + { + void writeReferenceFile(GeodeticReferenceFile &geofile, std::string output_file) + { + std::ofstream output_stream(output_file, std::ios::binary); + + if (geofile.file_type == GEO_TYPE) + { + GEO_GeodeticReferenceFile &gsofile = (GEO_GeodeticReferenceFile &)geofile; + + // Write main header + output_stream.write((char *)&gsofile.file_type, 1); + output_stream.write((char *)&gsofile.utc_timestamp_seconds, 8); + + // Write GEO header + output_stream.write((char *)&gsofile.norad, 4); + output_stream.write((char *)&gsofile.position_longitude, 8); + output_stream.write((char *)&gsofile.position_height, 8); + output_stream.write((char *)&gsofile.projection_type, 1); + output_stream.write((char *)&gsofile.image_width, 4); + output_stream.write((char *)&gsofile.image_height, 4); + output_stream.write((char *)&gsofile.horizontal_scale, 8); + output_stream.write((char *)&gsofile.vertical_scale, 8); + output_stream.write((char *)&gsofile.horizontal_offset, 8); + output_stream.write((char *)&gsofile.vertical_offset, 8); + output_stream.write((char *)&gsofile.proj_sweep_x, 1); + } + else if (geofile.file_type == LEO_TYPE) + { + LEO_GeodeticReferenceFile &leofile = (LEO_GeodeticReferenceFile &)geofile; + + // Write main header + output_stream.write((char *)&leofile.file_type, 1); + output_stream.write((char *)&leofile.utc_timestamp_seconds, 8); + + // Compute header size and TLE sizes + leofile.tle_line1_length = leofile.tle_line1_data.size(); + leofile.tle_line2_length = leofile.tle_line2_data.size(); + + // Write LEO header + output_stream.write((char *)&leofile.norad, 4); + output_stream.write((char *)&leofile.tle_line1_length, 2); + output_stream.write((char *)leofile.tle_line1_data.data(), leofile.tle_line1_length); + output_stream.write((char *)&leofile.tle_line2_length, 2); + output_stream.write((char *)leofile.tle_line2_data.data(), leofile.tle_line2_length); + + output_stream.write((char *)&leofile.scan_angle, 8); + output_stream.write((char *)&leofile.roll_offset, 8); + output_stream.write((char *)&leofile.pitch_offset, 8); + output_stream.write((char *)&leofile.yaw_offset, 8); + output_stream.write((char *)&leofile.time_offset, 8); + output_stream.write((char *)&leofile.image_width, 4); + output_stream.write((char *)&leofile.invert_scan, 1); + + leofile.timestamp_count = leofile.utc_timestamps.size(); + output_stream.write((char *)&leofile.timestamp_count, 8); + for (int i = 0; i < (int)leofile.timestamp_count; i++) + output_stream.write((char *)&leofile.utc_timestamps[i], 8); + } + + output_stream.close(); + } + + std::shared_ptr readReferenceFile(std::string input_file) + { + std::ifstream input_stream(input_file, std::ios::binary); + + uint8_t file_type; + input_stream.read((char *)&file_type, 1); + if (file_type == GEO_TYPE) + { + GEO_GeodeticReferenceFile geofile; + + // Main header + input_stream.read((char *)&geofile.utc_timestamp_seconds, 8); + + // Read GEO header + input_stream.read((char *)&geofile.norad, 4); + input_stream.read((char *)&geofile.position_longitude, 8); + input_stream.read((char *)&geofile.position_height, 8); + input_stream.read((char *)&geofile.projection_type, 1); + input_stream.read((char *)&geofile.image_width, 4); + input_stream.read((char *)&geofile.image_height, 4); + input_stream.read((char *)&geofile.horizontal_scale, 8); + input_stream.read((char *)&geofile.vertical_scale, 8); + input_stream.read((char *)&geofile.horizontal_offset, 8); + input_stream.read((char *)&geofile.vertical_offset, 8); + input_stream.read((char *)&geofile.proj_sweep_x, 1); + input_stream.close(); + return std::make_shared(geofile); + } + else if (file_type == LEO_TYPE) + { + LEO_GeodeticReferenceFile leofile; + + // Main header + input_stream.read((char *)&leofile.utc_timestamp_seconds, 8); + + // Read LEO header + input_stream.read((char *)&leofile.norad, 4); + input_stream.read((char *)&leofile.tle_line1_length, 2); + { + char *str = new char[leofile.tle_line1_length]; + input_stream.read(str, leofile.tle_line1_length); + leofile.tle_line1_data = std::string(str); + delete[] str; + } + input_stream.read((char *)&leofile.tle_line2_length, 2); + { + char *str = new char[leofile.tle_line2_length]; + input_stream.read(str, leofile.tle_line2_length); + leofile.tle_line2_data = std::string(str); + delete[] str; + } + + input_stream.read((char *)&leofile.scan_angle, 8); + input_stream.read((char *)&leofile.roll_offset, 8); + input_stream.read((char *)&leofile.pitch_offset, 8); + input_stream.read((char *)&leofile.yaw_offset, 8); + input_stream.read((char *)&leofile.time_offset, 8); + input_stream.read((char *)&leofile.image_width, 4); + input_stream.read((char *)&leofile.invert_scan, 1); + input_stream.read((char *)&leofile.timestamp_count, 8); + for (int i = 0; i < (int)leofile.timestamp_count; i++) + { + double timestamp; + input_stream.read((char *)×tamp, 8); + leofile.utc_timestamps.push_back(timestamp); + } + + input_stream.close(); + return std::make_shared(leofile); + } + else + { + // Unknown + return std::make_shared(); + } + } + + LEO_GeodeticReferenceFile leoRefFileFromProjector(int norad, LEOScanProjectorSettings projector_settings) + { + projection::proj_file::LEO_GeodeticReferenceFile geofile; + + double rough_timestamp = average_common(projector_settings.utc_timestamps.begin(), projector_settings.utc_timestamps.end()); + + geofile.utc_timestamp_seconds = floor(rough_timestamp); + geofile.norad = norad; + geofile.tle_line1_data = projector_settings.sat_tle.line1; + geofile.tle_line2_data = projector_settings.sat_tle.line2; + + geofile.scan_angle = projector_settings.scan_angle; + geofile.roll_offset = projector_settings.roll_offset; + geofile.pitch_offset = projector_settings.pitch_offset; + geofile.yaw_offset = projector_settings.yaw_offset; + geofile.time_offset = projector_settings.time_offset; + geofile.image_width = projector_settings.image_width; + geofile.invert_scan = projector_settings.invert_scan; + geofile.utc_timestamps = projector_settings.utc_timestamps; + + return geofile; + } + + LEOScanProjectorSettings leoProjectionRefFile(LEO_GeodeticReferenceFile geofile) + { + LEOScanProjectorSettings projector_settings = { + (double)geofile.scan_angle, + (double)geofile.roll_offset, + (double)geofile.pitch_offset, + (double)geofile.yaw_offset, + (double)geofile.time_offset, + (int)geofile.image_width, + (bool)geofile.invert_scan, + {(int)geofile.norad, "UNKNOWN", geofile.tle_line1_data, geofile.tle_line2_data}, + geofile.utc_timestamps}; + + return projector_settings; + } + + GEOProjector geoProjectionRefFile(GEO_GeodeticReferenceFile geofile) + { + GEOProjector projector(geofile.position_longitude, + geofile.position_height, + geofile.image_width, + geofile.image_height, + geofile.horizontal_scale, + geofile.vertical_scale, + geofile.horizontal_offset, + geofile.horizontal_offset, + geofile.proj_sweep_x); + return projector; + } + }; + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/proj_file.h b/src-core/common/geodetic/projection/proj_file.h new file mode 100644 index 000000000..caf4722da --- /dev/null +++ b/src-core/common/geodetic/projection/proj_file.h @@ -0,0 +1,113 @@ +#pragma once + +#include +#include +#include +#include "leo_projection.h" +#include "geo_projection.h" +#include + +namespace geodetic +{ + namespace projection + { + namespace proj_file + { + enum GeoReferenceFileTypes + { + GEO_TYPE = 1, + LEO_TYPE = 2, + }; + + /* + Class used to write & read projection files used in SatDump for geo-referencing + They contain the following data : + + Header : + - File type marker : + 1 = GEO Sat Data + 2 = LEO Sat Data + Different data types (GEO, LEO, etc) require different data types, + indicated by that marker. More may come later. + + - 64-bit UTC seconds timestamps : + The timestamp is meant to give a rough time of acquisition. It + can either be an average for LEO satellites, or the time of + acquisition for GEO satellites. + */ + struct GeodeticReferenceFile + { + uint8_t file_type; // File type + uint64_t utc_timestamp_seconds; // Timestamp. Depending on the type it may be approximate or not + }; + + /* + GEO Type 1 file. + */ + struct GEO_GeodeticReferenceFile : public GeodeticReferenceFile + { + // Type 1 + uint32_t norad; // NORAD ID of the GEO Satellite + double position_longitude; // Longitude of the satellite at the time of aquisition + double position_height; // Orbit height at the time of acquisition + uint8_t projection_type; // 0 is full disk, 1 is cropped, etc + uint32_t image_width; // Width of the image + uint32_t image_height; // Height of the image + double horizontal_scale; // Horizontal scaling + double vertical_scale; // Vertical scaling + double horizontal_offset; // Horizontal offset + double vertical_offset; // Vertical offset + bool proj_sweep_x; // Scan axis + + GEO_GeodeticReferenceFile() + { + projection_type = 0; + file_type = GEO_TYPE; + } + }; + + /* + LEO Type 2 file. + */ + struct LEO_GeodeticReferenceFile : public GeodeticReferenceFile + { + // Type 2 + uint32_t norad; // NORAD ID of the LEO Satellite + uint16_t tle_line1_length; // Size of TLE Line 1 + std::string tle_line1_data; // TLE Line 1 data + uint16_t tle_line2_length; // Size of TLE Line 1 + std::string tle_line2_data; // TLE Line 1 data + + enum LEOProjectionType + { + SINGLE_SCANLINE, + }; + + uint8_t projection_type; // Projection type + + double scan_angle; + double roll_offset; + double pitch_offset; + double yaw_offset; + double time_offset; + uint32_t image_width; + bool invert_scan; + uint64_t timestamp_count; + std::vector utc_timestamps; // Timestamps. Must match each scanline of the image you will be working with + + LEO_GeodeticReferenceFile() + { + projection_type = SINGLE_SCANLINE; + file_type = LEO_TYPE; + } + }; + + // Functions + void writeReferenceFile(GeodeticReferenceFile &geofile, std::string output_file); + std::shared_ptr readReferenceFile(std::string input_file); + LEO_GeodeticReferenceFile leoRefFileFromProjector(int norad, LEOScanProjectorSettings projector_settings); + LEOScanProjectorSettings leoProjectionRefFile(LEO_GeodeticReferenceFile geofile); + GEOProjector geoProjectionRefFile(GEO_GeodeticReferenceFile geofile); + }; + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/satellite_reprojector.cpp b/src-core/common/geodetic/projection/satellite_reprojector.cpp new file mode 100644 index 000000000..a243f0761 --- /dev/null +++ b/src-core/common/geodetic/projection/satellite_reprojector.cpp @@ -0,0 +1,180 @@ +#include "satellite_reprojector.h" +#include "logger.h" +#include "resources.h" +#include "common/map/map_drawer.h" + +namespace geodetic +{ + namespace projection + { + void reprojectLEOtoProj(cimg_library::CImg image, + projection::LEOScanProjector &projector, + cimg_library::CImg &projected_image, + int channels, + std::function(float, float, int, int)> projectionFunction, + float opacity, + float *progress) + { + // Reproject + for (int currentScan = 0; currentScan < (int)image.height(); currentScan++) + { + // Now compute each pixel's lat / lon and plot it + for (double px = 0; px < image.width() - 1; px += 1) + { + geodetic_coords_t coords1, coords2; + int ret1 = projector.inverse(px, currentScan, coords1); + int ret2 = projector.inverse(px + 1, currentScan, coords2); + + if (ret1 || ret2) + continue; + + std::pair map_cc1 = projectionFunction(coords1.lat, coords1.lon, projected_image.height(), projected_image.width()); + std::pair map_cc2 = projectionFunction(coords2.lat, coords2.lon, projected_image.height(), projected_image.width()); + + unsigned char color[3] = {0, 0, 0}; + if (channels == 3) + { + color[0] = image[image.width() * image.height() * 0 + currentScan * image.width() + int(px)] >> 8; + color[1] = image[image.width() * image.height() * 1 + currentScan * image.width() + int(px)] >> 8; + color[2] = image[image.width() * image.height() * 2 + currentScan * image.width() + int(px)] >> 8; + } + else + { + color[0] = image[currentScan * image.width() + int(px)] >> 8; + color[1] = image[currentScan * image.width() + int(px)] >> 8; + color[2] = image[currentScan * image.width() + int(px)] >> 8; + } + + if (color[0] == 0 && color[1] == 0 && color[2] == 0) // Skip Black + continue; + + // This seems to glitch out sometimes... Need to check + double maxSize = projected_image.width() / 100.0; + if (abs(map_cc1.first - map_cc2.first) < maxSize && abs(map_cc1.second - map_cc2.second) < maxSize) + { + double circle_radius = sqrt(pow(int(map_cc1.first - map_cc2.first), 2) + pow(int(map_cc1.second - map_cc2.second), 2)); + projected_image.draw_circle(map_cc1.first, map_cc1.second, ceil(circle_radius), color, 0.4 * opacity); + } + + if (progress != nullptr) + *progress = float(currentScan) / float(image.height()); + } + + //logger->info(std::to_string(currentScan)); + + //logger->critical(std::to_string(currentScan)); + } + } + + void reprojectGEOtoProj(cimg_library::CImg image, + projection::GEOProjector &projector, + cimg_library::CImg &projected_image, + int channels, + std::function(float, float, int, int)> projectionFunction, + float opacity, float *progress) + { + for (double lat = -90; lat < 90; lat += 0.01) + { + for (double lon = -180; lon < 180; lon += 0.01) + { + int x, y; + if (projector.forward(lon, lat, x, y)) + continue; + + std::pair map_cc1 = projectionFunction(lat, lon, projected_image.height(), projected_image.width()); + + unsigned char color[3]; + if (channels == 3) + { + color[0] = image[image.width() * image.height() * 0 + y * image.width() + int(x)] >> 8; + color[1] = image[image.width() * image.height() * 1 + y * image.width() + int(x)] >> 8; + color[2] = image[image.width() * image.height() * 2 + y * image.width() + int(x)] >> 8; + } + else + { + color[0] = image[y * image.width() + int(x)] >> 8; + color[1] = image[y * image.width() + int(x)] >> 8; + color[2] = image[y * image.width() + int(x)] >> 8; + } + + if (color[0] == 0 && color[1] == 0 && color[2] == 0) // Skip Black + continue; + + if (color[0] >= 253 && color[1] >= 253 && color[2] >= 253) // Skip Full white, as it's usually filler on GEO (eg, xRIT) + continue; + + //logger->info(std::to_string(color[0]) + " " + std::to_string(color[1]) + " " + std::to_string(color[2])); + + projected_image.draw_point(map_cc1.first, map_cc1.second, color, opacity); + } + + if (progress != nullptr) + *progress = (lat + 90) / 180; + } + } + + void projectEQUIToproj(cimg_library::CImg image, cimg_library::CImg &projected_image, int channels, std::function(float, float, int, int)> toMapCoords, float opacity, float *progress) + { + for (double lat = -90; lat < 90; lat += 0.01) + { + for (double lon = -180; lon < 180; lon += 0.01) + { + int x = (lon / 360.0f) * image.width() + (image.width() / 2); + int y = image.height() - ((90.0f + lat) / 180.0f) * image.height(); + + if (x >= image.width() || y >= image.height()) + continue; + + std::pair map_cc1 = toMapCoords(lat, lon, projected_image.height(), projected_image.width()); + + unsigned char color[3]; + if (channels == 3) + { + color[0] = image[image.width() * image.height() * 0 + y * image.width() + int(x)] >> 8; + color[1] = image[image.width() * image.height() * 1 + y * image.width() + int(x)] >> 8; + color[2] = image[image.width() * image.height() * 2 + y * image.width() + int(x)] >> 8; + } + else + { + color[0] = image[y * image.width() + int(x)] >> 8; + color[1] = image[y * image.width() + int(x)] >> 8; + color[2] = image[y * image.width() + int(x)] >> 8; + } + + if (color[0] == 0 && color[1] == 0 && color[2] == 0) // Skip Black + continue; + + projected_image.draw_point(map_cc1.first, map_cc1.second, color, opacity); + } + + if (progress != nullptr) + *progress = (lat + 90) / 180; + } + } + + cimg_library::CImg projectLEOToEquirectangularMapped(cimg_library::CImg image, + projection::LEOScanProjector &projector, + int output_width, + int output_height, + int channels, + cimg_library::CImg projected_image, + std::function(float, float, int, int)> toMapCoords + + ) + { + // Output mapped data + if (projected_image.width() == 1 && projected_image.height() == 1) + projected_image = cimg_library::CImg(output_width, output_height, 1, 3, 0); + + reprojectLEOtoProj(image, projector, projected_image, channels, toMapCoords); + + unsigned char color[3] = {0, 255, 0}; + map::drawProjectedMapShapefile({resources::getResourcePath("maps/ne_10m_admin_0_countries.shp")}, + projected_image, + color, + toMapCoords); + + return projected_image; + } + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/satellite_reprojector.h b/src-core/common/geodetic/projection/satellite_reprojector.h new file mode 100644 index 000000000..e1727d414 --- /dev/null +++ b/src-core/common/geodetic/projection/satellite_reprojector.h @@ -0,0 +1,64 @@ +#pragma once +#define cimg_use_png +#define cimg_display 0 +#include "CImg.h" +#include +#include "leo_projection.h" +#include "geo_projection.h" + +/* +Implementation of a function capable of plotting LEO satellite image (or anything that +can be done with a LEOScanProjector) to an equirectangular projection for easier viewing. +This has some known defects right now such as leaving gaps in the image. It's something I +will have to fix at some point... Still thinking about the best approach :-) +*/ +namespace geodetic +{ + namespace projection + { + // Reproject LEO imagery + void reprojectLEOtoProj( + cimg_library::CImg image, // Input image to project + projection::LEOScanProjector &projector, // LEO Projector + cimg_library::CImg &projected_image, // Optional input image + int channels, // Channels + std::function(float, float, int, int)> projectionFunction, // Optional projection function, default is equirectangular + float opacity = 1.0f, // Optional opacity + float *progress = nullptr // Optional progress value + ); + + // Reproject GEO imagery + void reprojectGEOtoProj( + cimg_library::CImg image, // Input image to project + projection::GEOProjector &projector, // GEO Projector + cimg_library::CImg &projected_image, // Optional input image + int channels, // Channels + std::function(float, float, int, int)> projectionFunction, // Optional projection function, default is equirectangular + float opacity = 1.0f, // Optional opacity + float *progress = nullptr // Optional progress value + ); + + void projectEQUIToproj(cimg_library::CImg image, + cimg_library::CImg &projected_image, + int channels, + std::function(float, float, int, int)> toMapCoords, + float opacity = 1.0f, + float *progress = nullptr); + + // Reproject LEO imagery to an equirectangular projection + cimg_library::CImg projectLEOToEquirectangularMapped( + cimg_library::CImg image, // Input image to project + projection::LEOScanProjector &projector, // LEO Projector + int output_width, // Output map width + int output_height, // Output map height + int channels = 1, // Channel count + cimg_library::CImg projected_image = cimg_library::CImg(1, 1, 1, 1, 0), // Optional input image + std::function(float, float, int, int)> toMapCoords = [](float lat, float lon, int map_height, int map_width) -> std::pair + { + int imageLat = map_height - ((90.0f + lat) / 180.0f) * map_height; + int imageLon = (lon / 360.0f) * map_width + (map_width / 2); + return {imageLon, imageLat}; + } // Optional projection function, default is equirectangular + ); + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/stereo.cpp b/src-core/common/geodetic/projection/stereo.cpp new file mode 100644 index 000000000..770dfc32a --- /dev/null +++ b/src-core/common/geodetic/projection/stereo.cpp @@ -0,0 +1,284 @@ +#include "stereo.h" +#include +#include "../wgs84.h" +#ifndef M_PI_2 +#define M_PI_2 1.57079632679489661923 /* pi/2 */ +#endif +#ifndef M_PI_4 +#define M_PI_4 0.78539816339744830962 /* pi/4 */ +#endif + +/* +** This file was adapted and simplified from libproj, and the below +** notice kept as credits. +** +** libproj -- library of cartographic projections +** +** Copyright (c) 2004 Gerald I. Evenden +** Copyright (c) 2012 Martin Raspaud +** +** See also (section 4.4.3.2): +** https://www.cgms-info.org/documents/pdf_cgms_03.pdf +** +** Permission is hereby granted, free of charge, to any person obtaining +** a copy of this software and associated documentation files (the +** "Software"), to deal in the Software without restriction, including +** without limitation the rights to use, copy, modify, merge, publish, +** distribute, sublicense, and/or sell copies of the Software, and to +** permit persons to whom the Software is furnished to do so, subject to +** the following conditions: +** +** The above copyright notice and this permission notice shall be +** included in all copies or substantial portions of the Software. +** +** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +*/ + +#define EPS10 1.e-10 + +namespace geodetic +{ + namespace projection + { + double pj_tsfn(double phi, double sinphi, double e) + { + double cosphi = cos(phi); + return exp(e * atanh(e * sinphi)) * (sinphi > 0 ? cosphi / (1 + sinphi) : (1 - sinphi) / cosphi); + } + + static double ssfn_(double phit, double sinphi, double eccen) + { + sinphi *= eccen; + return (tan(.5 * (M_PI_2 + phit)) * pow((1. - sinphi) / (1. + sinphi), .5 * eccen)); + } + + int StereoProjection::init(double latitude, double longitude) + { + lon_0 = longitude; // The projection's longitude + + // Constants, WGS84 + e = WGS84::e; + phi0 = latitude * 0.01745329; + a = WGS84::a * 1000; + es = WGS84::es; + one_es = WGS84::one_es; + + if (es == 0.0) + { + // Illegal + return 1; + } + + k0 = .994; + //x0 = 2000000.; + //y0 = 2000000.; + phits = M_PI_2; + lam0 = 0.; + + // Setup + double t; + + if (fabs((t = fabs(phi0)) - M_PI_2) < EPS10) + mode = phi0 < 0. ? S_POLE : N_POLE; + else + mode = t > EPS10 ? OBLIQ : EQUIT; + phits = fabs(phits); + + if (es != 0.0) + { + double X; + + switch (mode) + { + case N_POLE: + case S_POLE: + if (fabs(phits - M_PI_2) < EPS10) + akm1 = 2. * k0 / + sqrt(pow(1 + e, 1 + e) * pow(1 - e, 1 - e)); + else + { + t = sin(phits); + akm1 = cos(phits) / pj_tsfn(phits, t, e); + t *= e; + akm1 /= sqrt(1. - t * t); + } + break; + case EQUIT: + case OBLIQ: + t = sin(phi0); + X = 2. * atan(ssfn_(phi0, t, e)) - M_PI_2; + t *= e; + akm1 = 2. * k0 * cos(phi0) / sqrt(1. - t * t); + sinX1 = sin(X); + cosX1 = cos(X); + break; + } + } + else + { + switch (mode) + { + case OBLIQ: + sinX1 = sin(phi0); + cosX1 = cos(phi0); + /*-fallthrough*/ + case EQUIT: + akm1 = 2. * k0; + break; + case S_POLE: + case N_POLE: + akm1 = fabs(phits - M_PI_2) >= EPS10 ? cos(phits) / tan(M_PI_4 - .5 * phits) : 2. * k0; + break; + } + } + return 0; + } + + int StereoProjection::forward(double lon, double lat, double &x, double &y) + { + x = y = 0; // Safety + + // Shift longitudes + lon -= lon_0; + if (lon < -180) + lon = lon + 360; + if (lon > 180) + lon = lon - 360; + + // To radians + double phi = lat * 0.01745329, lam = lon * 0.01745329; + + double coslam, sinlam, sinX = 0.0, cosX = 0.0, A = 0.0, sinphi; + + coslam = cos(lam); + sinlam = sin(lam); + sinphi = sin(phi); + if (mode == OBLIQ || mode == EQUIT) + { + const double X = 2. * atan(ssfn_(phi, sinphi, e)) - M_PI_2; + sinX = sin(X); + cosX = cos(X); + } + + switch (mode) + { + case OBLIQ: + { + const double denom = cosX1 * (1. + sinX1 * sinX + cosX1 * cosX * coslam); + if (denom == 0) + { + // Illegal + return 1; + } + A = akm1 / denom; + y = A * (cosX1 * sinX - sinX1 * cosX * coslam); + x = A * cosX; + break; + } + + case EQUIT: + /* avoid zero division */ + if (1. + cosX * coslam == 0.0) + { + y = HUGE_VAL; + } + else + { + A = akm1 / (1. + cosX * coslam); + y = A * sinX; + } + x = A * cosX; + break; + + case S_POLE: + phi = -phi; + coslam = -coslam; + sinphi = -sinphi; + /*-fallthrough*/ + case N_POLE: + if (fabs(phi - M_PI_2) < 1e-15) + x = 0; + else + x = akm1 * pj_tsfn(phi, sinphi, e); + y = -x * coslam; + break; + } + + x = x * sinlam; + return 0; + } + + int StereoProjection::inverse(double x, double y, double &lon, double &lat) + { + lon = lat = 0.0; + double phi = 0, lam = 0; + + double cosphi, sinphi, tp = 0.0, phi_l = 0.0, rho, halfe = 0.0, halfpi = 0.0; + + rho = hypot(x, y); + + switch (mode) + { + case OBLIQ: + case EQUIT: + tp = 2. * atan2(rho * cosX1, akm1); + cosphi = cos(tp); + sinphi = sin(tp); + if (rho == 0.0) + phi_l = asin(cosphi * sinX1); + else + phi_l = asin(cosphi * sinX1 + (y * sinphi * cosX1 / rho)); + + tp = tan(.5 * (M_PI_2 + phi_l)); + x *= sinphi; + y = rho * cosX1 * cosphi - y * sinX1 * sinphi; + halfpi = M_PI_2; + halfe = .5 * e; + break; + case N_POLE: + y = -y; + /*-fallthrough*/ + case S_POLE: + tp = -rho / akm1; + phi_l = M_PI_2 - 2. * atan(tp); + halfpi = -M_PI_2; + halfe = -.5 * e; + break; + } + + for (int i = 8; i > 0; --i) + { + sinphi = e * sin(phi_l); + phi = 2. * atan(tp * pow((1. + sinphi) / (1. - sinphi), halfe)) - halfpi; + if (fabs(phi_l - phi) < 1.e-10) + { + if (mode == S_POLE) + phi = -phi; + lam = (x == 0. && y == 0.) ? 0. : atan2(x, y); + + // To degs + lat = phi * 57.29578; + lon = lam * 57.29578; + + // Shift longitudes back to reference 0 + lon += lon_0; + if (lon < -180) + lon = lon + 360; + if (lon > 180) + lon = lon - 360; + + return 0; + } + phi_l = phi; + } + + return 1; + } + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/stereo.h b/src-core/common/geodetic/projection/stereo.h new file mode 100644 index 000000000..01459baad --- /dev/null +++ b/src-core/common/geodetic/projection/stereo.h @@ -0,0 +1,58 @@ +#pragma once + +/* +Implementation of a standard GEOS projection, adapted from libproj. +Some variables are hardcoded for the intended usecase, making some +degree of tuning unecessary. +Uses the WGS84 ellipsoid. +*/ +namespace geodetic +{ + namespace projection + { + class StereoProjection + { + private: + enum Mode + { + S_POLE = 0, + N_POLE = 1, + OBLIQ = 2, + EQUIT = 3 + }; + + double phits; + double sinX1; + double cosX1; + double akm1; + enum Mode mode; + + double e; + double phi0; + double a; + double es; + double one_es; + + double k0; + //double x0; + //double y0; + double lam0; + + double lon_0; + + public: + StereoProjection() + { + init(false, 0); + } + StereoProjection(double latitude, double longitude) + { + init(latitude, longitude); + } + + int init(double latitude, double longitude); // return value of 1 => Error + int forward(double lon, double lat, double &x, double &y); // return value of 1 => Error + int inverse(double x, double y, double &lon, double &lat); // return value of 1 => Error + }; + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/tpers.cpp b/src-core/common/geodetic/projection/tpers.cpp new file mode 100644 index 000000000..e05b0d9b3 --- /dev/null +++ b/src-core/common/geodetic/projection/tpers.cpp @@ -0,0 +1,238 @@ +#include "tpers.h" +#include +#include "../wgs84.h" +#ifndef M_PI_2 +#define M_PI_2 1.57079632679489661923 /* pi/2 */ +#endif + +/* +** This file was adapted and simplified from libproj, and the below +** notice kept as credits. +** +** libproj -- library of cartographic projections +** +** Permission is hereby granted, free of charge, to any person obtaining +** a copy of this software and associated documentation files (the +** "Software"), to deal in the Software without restriction, including +** without limitation the rights to use, copy, modify, merge, publish, +** distribute, sublicense, and/or sell copies of the Software, and to +** permit persons to whom the Software is furnished to do so, subject to +** the following conditions: +** +** The above copyright notice and this permission notice shall be +** included in all copies or substantial portions of the Software. +** +** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, +** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF +** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. +** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY +** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, +** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE +** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. +*/ + +#define EPS10 1.e-10 + +namespace geodetic +{ + namespace projection + { + int TPERSProjection::init(double altitude, double longitude, double latitude, double tilt, double azi) + { + lon_0 = longitude; + + double omega, gamma; + omega = tilt * 0.01745329; // to rads + gamma = azi * 0.01745329; // to rads + this->tilt = 1; + cg = cos(gamma); + sg = sin(gamma); + cw = cos(omega); + sw = sin(omega); + + height = altitude; + phi0 = latitude * 0.01745329; + a = WGS84::a* 1000; + + if (fabs(fabs(phi0) - M_PI_2) < EPS10) + mode = phi0 < 0. ? S_POLE : N_POLE; + else if (fabs(phi0) < EPS10) + mode = EQUIT; + else + { + mode = OBLIQ; + sinph0 = sin(phi0); + cosph0 = cos(phi0); + } + + pn1 = height / a; // normalize by radius + + if (pn1 <= 0 || pn1 > 1e10) + { + // Illegal!! + return 1; + } + + p = 1. + pn1; + rp = 1. / p; + h = 1. / pn1; + pfact = (p + 1.) * h; + es = 0.; + + return 0; + } + + int TPERSProjection::forward(double lon, double lat, double &x, double &y) + { + x = y = 0; // Safety + + // Shift longitudes to use the sat's as a reference + lon -= lon_0; + if (lon < -180) + lon = lon + 360; + if (lon > 180) + lon = lon - 360; + + // To radians + double phi = lat * 0.01745329, lam = lon * 0.01745329; + + double coslam, cosphi, sinphi; + + sinphi = sin(phi); + cosphi = cos(phi); + coslam = cos(lam); + + switch (mode) + { + case OBLIQ: + y = sinph0 * sinphi + cosph0 * cosphi * coslam; + break; + case EQUIT: + y = cosphi * coslam; + break; + case S_POLE: + y = -sinphi; + break; + case N_POLE: + y = sinphi; + break; + } + + if (y < rp) + { + x = y = 2e10; // Trigger error + return 1; + } + + y = pn1 / (p - y); + x = y * cosphi * sin(lam); + + switch (mode) + { + case OBLIQ: + y *= (cosph0 * sinphi - + sinph0 * cosphi * coslam); + break; + case EQUIT: + y *= sinphi; + break; + case N_POLE: + coslam = -coslam; + /*-fallthrough*/ + case S_POLE: + y *= cosphi * coslam; + break; + } + + if (tilt) + { + double yt, ba; + + yt = y * cg + x * sg; + ba = 1. / (yt * sw * h + cw); + x = (x * cg - y * sg) * cw * ba; + y = yt * ba; + } + + return 0; + } + + int TPERSProjection::inverse(double x, double y, double &lon, double &lat) + { + lon = lat = 0.0; + + double phi = 0, lam = 0; + double rh; + + if (tilt) + { + double bm, bq, yt; + + yt = 1. / (pn1 - y * sw); + bm = pn1 * x * yt; + bq = pn1 * y * cw * yt; + x = bm * cg + bq * sg; + y = bq * cg - bm * sg; + } + + rh = hypot(x, y); + + if (fabs(rh) <= EPS10) + { + lam = 0.; + phi = phi0; + } + else + { + double cosz, sinz; + sinz = 1. - rh * rh * pfact; + + if (sinz < 0.) + { + // Illegal + lon = lat = 2e10; // Trigger error + return 1; + } + + sinz = (p - sqrt(sinz)) / (pn1 / rh + rh / pn1); + cosz = sqrt(1. - sinz * sinz); + + switch (mode) + { + case OBLIQ: + phi = asin(cosz * sinph0 + y * sinz * cosph0 / rh); + y = (cosz - sinph0 * sin(phi)) * rh; + x *= sinz * cosph0; + break; + case EQUIT: + phi = asin(y * sinz / rh); + y = cosz * rh; + x *= sinz; + break; + case N_POLE: + phi = asin(cosz); + y = -y; + break; + case S_POLE: + phi = -asin(cosz); + break; + } + + lam = atan2(x, y); + } + + // To degs + lat = phi * 57.29578; + lon = lam * 57.29578; + + // Shift longitudes back to reference 0 + lon += lon_0; + if (lon < -180) + lon = lon + 360; + if (lon > 180) + lon = lon - 360; + + return 0; + } + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/projection/tpers.h b/src-core/common/geodetic/projection/tpers.h new file mode 100644 index 000000000..f51fb45bb --- /dev/null +++ b/src-core/common/geodetic/projection/tpers.h @@ -0,0 +1,60 @@ +#pragma once + +/* +Implementation of a titled perspective projection, for example useful to +project data from the point of a view of a LEO satellite. +This was adapted from libproj. +Uses the WGS84 ellipsoid. +*/ +namespace geodetic +{ + namespace projection + { + class TPERSProjection + { + private: + enum Mode + { + N_POLE = 0, + S_POLE = 1, + EQUIT = 2, + OBLIQ = 3 + }; + + double height; + double sinph0; + double cosph0; + double p; + double rp; + double pn1; + double pfact; + double h; + double cg; + double sg; + double sw; + double cw; + enum Mode mode; + int tilt; + + double phi0; + double a; + double es; + + double lon_0; + + public: + TPERSProjection() + { + init(0, 0, 0, 0, 0); + } + TPERSProjection(double altitude, double longitude, double latitude, double tilt, double azi) + { + init(altitude, longitude, latitude, tilt, azi); + } + + int init(double altitude, double longitude, double latitude, double tilt, double azi); // return value of 1 => Error + int forward(double lon, double lat, double &x, double &y); // return value of 1 => Error + int inverse(double x, double y, double &lon, double &lat); // return value of 1 => Error + }; + }; +}; \ No newline at end of file diff --git a/src-core/common/geodetic/vincentys_calculations.cpp b/src-core/common/geodetic/vincentys_calculations.cpp new file mode 100644 index 000000000..c034f0e3d --- /dev/null +++ b/src-core/common/geodetic/vincentys_calculations.cpp @@ -0,0 +1,284 @@ +#include "vincentys_calculations.h" +#include "wgs84.h" + +#define M_2PI (M_PI * 2) + +/* +The code in this file was ported to C++ and slightly modified +from https://github.com/airbreather/Gavaghan.Geodesy/blob/master/Source/Gavaghan.Geodesy/GeodeticCalculator.cs + +Original license : +--------------------------------------------------------------------------------------- +Gavaghan.Geodesy by Mike Gavaghan + +http://www.gavaghan.org/blog/free-source-code/geodesy-library-vincentys-formula/ + +This code may be freely used and modified on any personal or professional +project. It comes with no warranty. + +BitCoin tips graciously accepted at 1FB63FYQMy7hpC2ANVhZ5mSgAZEtY1aVLf +--------------------------------------------------------------------------------------- +*/ +namespace geodetic +{ + geodetic_coords_t vincentys_forward(geodetic_coords_t start, double initialBearing, double distance, double &finalBearing, double tolerance) + { + const double &a = WGS84::a; + const double &b = WGS84::b; + const double aSquared = a * a; + const double bSquared = b * b; + const double &f = WGS84::f; + + start.toRads(); // Ensure + + double phi1 = start.lat; + double alpha1 = initialBearing; + double cosAlpha1 = cos(alpha1); + double sinAlpha1 = sin(alpha1); + double s = distance * 1000; // To meters from Km + double tanU1 = (1.0 - f) * tan(phi1); + double cosU1 = 1.0 / sqrt(1.0 + tanU1 * tanU1); + double sinU1 = tanU1 * cosU1; + + // eq. 1 + double sigma1 = atan2(tanU1, cosAlpha1); + + // eq. 2 + double sinAlpha = cosU1 * sinAlpha1; + double sin2Alpha = sinAlpha * sinAlpha; + double cos2Alpha = 1 - sin2Alpha; + double uSquared = cos2Alpha * (aSquared - bSquared) / bSquared; + + // eq. 3 + double A = 1 + (uSquared / 16384) * (4096 + uSquared * (-768 + uSquared * (320 - 175 * uSquared))); + + // eq. 4 + double B = (uSquared / 1024) * (256 + uSquared * (-128 + uSquared * (74 - 47 * uSquared))); + + // iterate until there is a negligible change in sigma + double deltaSigma; + double sOverbA = s / (b * A); + double sigma = sOverbA; + double sinSigma; + double prevSigma = sOverbA; + double sigmaM2; + double cosSigmaM2; + double cos2SigmaM2; + + for (;;) + { + // eq. 5 + sigmaM2 = 2.0 * sigma1 + sigma; + cosSigmaM2 = cos(sigmaM2); + cos2SigmaM2 = cosSigmaM2 * cosSigmaM2; + sinSigma = sin(sigma); + double cosSignma = cos(sigma); + + // eq. 6 + deltaSigma = B * sinSigma * (cosSigmaM2 + (B / 4.0) * (cosSignma * (-1 + 2 * cos2SigmaM2) - (B / 6.0) * cosSigmaM2 * (-3 + 4 * sinSigma * sinSigma) * (-3 + 4 * cos2SigmaM2))); + + // eq. 7 + sigma = sOverbA + deltaSigma; + + // break after converging to tolerance + if (abs(sigma - prevSigma) < tolerance) + break; + + prevSigma = sigma; + } + + sigmaM2 = 2.0 * sigma1 + sigma; + cosSigmaM2 = cos(sigmaM2); + cos2SigmaM2 = cosSigmaM2 * cosSigmaM2; + + double cosSigma = cos(sigma); + sinSigma = sin(sigma); + + // eq. 8 + double sinU1sinSigma_cosU1cosSigmacosAlpha1 = sinU1 * sinSigma - cosU1 * cosSigma * cosAlpha1; + double phi2 = atan2(sinU1 * cosSigma + cosU1 * sinSigma * cosAlpha1, (1.0 - f) * sqrt(sin2Alpha + (sinU1sinSigma_cosU1cosSigmacosAlpha1 * sinU1sinSigma_cosU1cosSigmacosAlpha1))); + + // eq. 9 + // This fixes the pole crossing defect spotted by Matt Feemster. When a path + // passes a pole and essentially crosses a line of latitude twice - once in + // each direction - the longitude calculation got messed up. Using Atan2 + // instead of Atan fixes the defect. The change is in the next 3 lines. + //double tanLambda = sinSigma * sinAlpha1 / (cosU1 * cosSigma - sinU1*sinSigma*cosAlpha1); + //double lambda = atan(tanLambda); + double lambda = atan2(sinSigma * sinAlpha1, cosU1 * cosSigma - sinU1 * sinSigma * cosAlpha1); + + // eq. 10 + double C = (f / 16) * cos2Alpha * (4 + f * (4 - 3 * cos2Alpha)); + + // eq. 11 + double L = lambda - (1 - C) * f * sinAlpha * (sigma + C * sinSigma * (cosSigmaM2 + C * cosSigma * (-1 + 2 * cos2SigmaM2))); + + // eq. 12 + double alpha2 = atan2(sinAlpha, -sinU1 * sinSigma + cosU1 * cosSigma * cosAlpha1); + + // build result + finalBearing = alpha2; + + return geodetic_coords_t(phi2, start.lon + L, start.alt, true); + } + + geodetic_curve_t vincentys_inverse(geodetic_coords_t start, geodetic_coords_t end, double tolerance) + { + // + // All equation numbers refer back to Vincenty's publication: + // See http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf + // + + // get constants + const double &a = WGS84::a; + const double &b = WGS84::b; + const double &f = WGS84::f; + + start.toRads(); + end.toRads(); + + // get parameters as radians + double phi1 = start.lat; + double lambda1 = start.lon; + double phi2 = end.lat; + double lambda2 = end.lon; + + // calculations + double a2 = a * a; + double b2 = b * b; + double a2b2b2 = (a2 - b2) / b2; + + double omega = lambda2 - lambda1; + + double tanphi1 = tan(phi1); + double tanU1 = (1.0 - f) * tanphi1; + double U1 = atan(tanU1); + double sinU1 = sin(U1); + double cosU1 = cos(U1); + + double tanphi2 = tan(phi2); + double tanU2 = (1.0 - f) * tanphi2; + double U2 = atan(tanU2); + double sinU2 = sin(U2); + double cosU2 = cos(U2); + + double sinU1sinU2 = sinU1 * sinU2; + double cosU1sinU2 = cosU1 * sinU2; + double sinU1cosU2 = sinU1 * cosU2; + double cosU1cosU2 = cosU1 * cosU2; + + // eq. 13 + double lambda = omega; + + // intermediates we'll need to compute 's' + double A = 0.0; + double B = 0.0; + double sigma = 0.0; + double deltasigma = 0.0; + double lambda0; + bool converged = false; + + for (int i = 0; i < 20; i++) + { + lambda0 = lambda; + + double sinlambda = sin(lambda); + double coslambda = cos(lambda); + + // eq. 14 + double cosU1sinU2_sinU2cosU2coslambda = cosU1sinU2 - sinU1cosU2 * coslambda; + double sin2sigma = (cosU2 * sinlambda * cosU2 * sinlambda) + (cosU1sinU2_sinU2cosU2coslambda * cosU1sinU2_sinU2cosU2coslambda); + double sinsigma = sqrt(sin2sigma); + + // eq. 15 + double cossigma = sinU1sinU2 + (cosU1cosU2 * coslambda); + + // eq. 16 + sigma = atan2(sinsigma, cossigma); + + // eq. 17 Careful! sin2sigma might be almost 0! + double sinalpha = (sin2sigma == 0) ? 0.0 : cosU1cosU2 * sinlambda / sinsigma; + double alpha = asin(sinalpha); + double cosalpha = cos(alpha); + double cos2alpha = cosalpha * cosalpha; + + // eq. 18 Careful! cos2alpha might be almost 0! + double cos2sigmam = cos2alpha == 0.0 ? 0.0 : cossigma - 2 * sinU1sinU2 / cos2alpha; + double u2 = cos2alpha * a2b2b2; + + double cos2sigmam2 = cos2sigmam * cos2sigmam; + + // eq. 3 + A = 1.0 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2))); + + // eq. 4 + B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2))); + + // eq. 6 + deltasigma = B * sinsigma * (cos2sigmam + B / 4 * (cossigma * (-1 + 2 * cos2sigmam2) - B / 6 * cos2sigmam * (-3 + 4 * sin2sigma) * (-3 + 4 * cos2sigmam2))); + + // eq. 10 + double C = f / 16 * cos2alpha * (4 + f * (4 - 3 * cos2alpha)); + + // eq. 11 (modified) + lambda = omega + (1 - C) * f * sinalpha * (sigma + C * sinsigma * (cos2sigmam + C * cossigma * (-1 + 2 * cos2sigmam2))); + + if (i < 2) + continue; + + // see how much improvement we got + double change = abs((lambda - lambda0) / lambda); + + if (change < tolerance) + { + converged = true; + break; + } + } + + // eq. 19 + double s = b * A * (sigma - deltasigma); + double alpha1 = 0; + double alpha2 = 0; + + // didn't converge? must be N/S + if (!converged) + { + if (phi1 > phi2) + { + alpha1 = 180 * DEG_TO_RAD; + alpha2 = 0 * DEG_TO_RAD; + } + else if (phi1 < phi2) + { + alpha1 = 0 * DEG_TO_RAD; + alpha2 = 180 * DEG_TO_RAD; + } + else + { + //alpha1 = Angle.NaN; + //alpha2 = Angle.NaN; + //logger->error("Error"); + } + } + else + { + // eq. 20 + alpha1 = atan2(cosU2 * sin(lambda), (cosU1sinU2 - sinU1cosU2 * cos(lambda))); + if (alpha1 < 0.0) + alpha1 += M_2PI; + + // eq. 21 + alpha2 = atan2(cosU1 * sin(lambda), (-sinU1cosU2 + cosU1sinU2 * cos(lambda))) + M_PI; + if (alpha2 < 0.0) + alpha2 += M_2PI; + } + + if (alpha1 >= M_2PI) + alpha1 = alpha1 - M_2PI; + if (alpha2 >= M_2PI) + alpha2 = alpha2 - M_2PI; + + return geodetic_curve_t(s / 1000, alpha1, alpha2, true); + } +}; \ No newline at end of file diff --git a/src-core/common/geodetic/vincentys_calculations.h b/src-core/common/geodetic/vincentys_calculations.h new file mode 100644 index 000000000..763f19ef7 --- /dev/null +++ b/src-core/common/geodetic/vincentys_calculations.h @@ -0,0 +1,53 @@ +#pragma once + +#include "geodetic_coordinates.h" + +namespace geodetic +{ + struct geodetic_curve_t + { + double distance; + double azimuth; + double reverse_azimuth; + + bool is_radians = false; + + geodetic_curve_t(double distance, double azimuth, double reverse_azimuth, bool radians = false) + { + this->distance = distance; + this->azimuth = azimuth; + this->reverse_azimuth = reverse_azimuth; + this->is_radians = radians; + } + + geodetic_curve_t toRads() // Convert to Rads if not done already + { + if (!is_radians) + { + azimuth *= DEG_TO_RAD; + reverse_azimuth *= DEG_TO_RAD; + is_radians = true; + } + + return *this; + } + + geodetic_curve_t toDegs() // Convert to Degs if not done already + { + if (is_radians) + { + azimuth *= RAD_TO_DEG; + reverse_azimuth *= RAD_TO_DEG; + is_radians = false; + } + + return *this; + } + }; + + // Compute new position and bearing given start geodetic coordinates, bearing and displacement distance + geodetic_coords_t vincentys_forward(geodetic_coords_t start, double initialBearing, double distance, double &finalBearing, double tolerance = 1e-13); + + // Compute distance and bearing between 2 geodetic coordinates + geodetic_curve_t vincentys_inverse(geodetic_coords_t start, geodetic_coords_t end, double tolerance = 1e-13); +}; \ No newline at end of file diff --git a/src-core/common/geodetic/wgs84.h b/src-core/common/geodetic/wgs84.h new file mode 100644 index 000000000..2d6ff827c --- /dev/null +++ b/src-core/common/geodetic/wgs84.h @@ -0,0 +1,19 @@ +#pragma once + +namespace geodetic +{ + /* + WGS86 ellipsoid definition, used everywhere in SatDump (or should be!) + */ + namespace WGS84 + { + const double a = 6378.137; // Semimajor Axis + const double rf = 298.257223563; // Inverse flattening + const double f = 1.0 / rf; // Flattening + const double b = a * (1 - f); // Semiminor Axis + const double e = sqrt((pow(a, 2) - pow(b, 2)) / pow(a, 2)); // First eccentricity + const double e2 = sqrt((pow(a, 2) - pow(b, 2)) / pow(b, 2)); // Second eccentricity + const double es = pow(e, 2); // First eccentricity^2 + const double one_es = 1.0 - es; // First 1 - eccentricity^2 + }; +}; \ No newline at end of file diff --git a/src-core/common/projection/geo_projection.cpp b/src-core/common/projection/geo_projection.cpp deleted file mode 100644 index 1658e51b7..000000000 --- a/src-core/common/projection/geo_projection.cpp +++ /dev/null @@ -1,69 +0,0 @@ -#include "geo_projection.h" - -#include -#include "logger.h" - -namespace projection -{ - GEOProjector::GEOProjector(double sat_lon, - double sat_height, - int img_width, - int img_height, - double hscale, - double vscale, - double x_offset, - double y_offset, - bool sweep_x) : hscale(hscale), - vscale(vscale), - x_offset(x_offset), - y_offset(y_offset) - { - height = img_height; - width = img_width; - pj.init(sat_height * 1000, sat_lon, sweep_x); - } - - int GEOProjector::forward(double lon, double lat, int &img_x, int &img_y) - { - if (pj.forward(lon, lat, x, y)) - { - // Error / out of the image - img_x = -1; - img_y = -1; - return 1; - } - - image_x = x * hscale * (width / 2.0); - image_y = y * vscale * (height / 2.0); - - image_x += width / 2.0 + x_offset; - image_y += height / 2.0 + y_offset; - - img_x = image_x; - img_y = (height - 1) - image_y; - - return 0; - } - - int GEOProjector::inverse(int img_x, int img_y, double &lon, double &lat) - { - image_y = (height - 1) - img_y; - image_x = img_x; - - image_y -= height / 2.0 + y_offset; - image_x -= width / 2.0 + x_offset; - - y = image_y / (vscale * (height / 2.0)); - x = image_x / (hscale * (width / 2.0)); - - if (pj.inverse(x, y, lon, lat)) - { - // Error / out of the image - img_x = -1; - img_y = -1; - return 1; - } - - return 0; - } -}; \ No newline at end of file diff --git a/src-core/common/projection/geo_projection.h b/src-core/common/projection/geo_projection.h deleted file mode 100644 index bd74f59ca..000000000 --- a/src-core/common/projection/geo_projection.h +++ /dev/null @@ -1,39 +0,0 @@ -#pragma once - -#include "geos.h" - -/* -Code to reference a decoded image (or similar data) from a GEO satellite to Lat / Lon coordinates. - -The inverse function is currently broken, for an unknwon reason. -*/ -namespace projection -{ - class GEOProjector - { - private: - projection::GEOSProjection pj; - double height, width; - - double x, y; - double image_x, image_y; - - double hscale; - double vscale; - double x_offset; - double y_offset; - - public: - GEOProjector(double sat_lon, - double sat_height, - int img_width, - int img_height, - double hscale, - double vscale, - double x_offset, - double y_offset, - bool sweep_x); - int forward(double lon, double lat, int &img_x, int &img_y); - int inverse(int img_x, int img_y, double &lon, double &lat); - }; -}; \ No newline at end of file diff --git a/src-core/common/projection/geos.cpp b/src-core/common/projection/geos.cpp deleted file mode 100644 index b65b07cd9..000000000 --- a/src-core/common/projection/geos.cpp +++ /dev/null @@ -1,178 +0,0 @@ -#include "geos.h" -#include - -/* -** This file was adapted and simplified from libproj, and the below -** notice kept as credits. -** -** libproj -- library of cartographic projections -** -** Copyright (c) 2004 Gerald I. Evenden -** Copyright (c) 2012 Martin Raspaud -** -** See also (section 4.4.3.2): -** https://www.cgms-info.org/documents/pdf_cgms_03.pdf -** -** Permission is hereby granted, free of charge, to any person obtaining -** a copy of this software and associated documentation files (the -** "Software"), to deal in the Software without restriction, including -** without limitation the rights to use, copy, modify, merge, publish, -** distribute, sublicense, and/or sell copies of the Software, and to -** permit persons to whom the Software is furnished to do so, subject to -** the following conditions: -** -** The above copyright notice and this permission notice shall be -** included in all copies or substantial portions of the Software. -** -** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF -** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. -** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY -** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, -** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE -** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. -*/ - -namespace projection -{ - int GEOSProjection::init(double height, double longitude, bool sweep_x) - { - lon_0 = longitude; // The satellite's longitude - - // Constants, extracted from Proj - phi0 = 0; - a = 6.37814e+06; - es = 0.00669438; - one_es = 0.993306; - - // Orbit Height - h = height; - - // Scan axis - flip_axis = sweep_x; - - radius_g_1 = h / a; - if (radius_g_1 <= 0 || radius_g_1 > 1e10) - { - // Illegal case. - // Kept just in case but we shouldn't end up there unless the user makes a mistake... - return 1; - } - - // Init the rest - radius_g = 1. + radius_g_1; - C = radius_g * radius_g - 1.0; - - radius_p = sqrt(one_es); - radius_p2 = one_es; - radius_p_inv2 = one_es; - - return 0; - } - - int GEOSProjection::forward(double lon, double lat, double &x, double &y) - { - x = y = 0; // Safety - - // Shift longitudes to use the sat's as a reference - lon -= lon_0; - if (lon < -180) - lon = lon + 360; - if (lon > 180) - lon = lon - 360; - - // To radians - double phi = lat * 0.01745329, lam = lon * 0.01745329; - - double r, Vx, Vy, Vz, tmp; - - // Calculation of geocentric latitude. - phi = atan(radius_p2 * tan(phi)); - - // Calculation of the three components of the vector from satellite to position on earth surface (lon,lat). - r = (radius_p) / hypot(radius_p * cos(phi), sin(phi)); - Vx = r * cos(lam) * cos(phi); - Vy = r * sin(lam) * cos(phi); - Vz = r * sin(phi); - - // Check visibility. - if (((radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * radius_p_inv2) < 0.) - { - x = y = 2e10; // Trigger error - return 1; - } - - // Calculation based on view angles from satellite. - tmp = radius_g - Vx; - - if (flip_axis) - { - x = radius_g_1 * atan(Vy / hypot(Vz, tmp)); - y = radius_g_1 * atan(Vz / tmp); - } - else - { - x = radius_g_1 * atan(Vy / tmp); - y = radius_g_1 * atan(Vz / hypot(Vy, tmp)); - } - - return 0; - } - - int GEOSProjection::inverse(double x, double y, double &lon, double &lat) - { - lon = lat = 0.0; - double phi = 0, lam = 0; - - double Vx, Vy, Vz, a, b, k; - - // Setting three components of vector from satellite to position. - Vx = -1.0; - - if (flip_axis) - { - Vz = tan(y / radius_g_1); - Vy = tan(x / radius_g_1) * hypot(1.0, Vz); - } - else - { - Vy = tan(x / radius_g_1); - Vz = tan(y / radius_g_1) * hypot(1.0, Vy); - } - - // Calculation of terms in cubic equation and determinant. - a = Vz / radius_p; - a = Vy * Vy + a * a + Vx * Vx; - b = 2 * radius_g * Vx; - const double det = (b * b) - 4 * a * C; - if (det < 0.0) - { - lon = lat = 2e10; // Trigger error - return 1; - } - - // Calculation of three components of vector from satellite to position. - k = (-b - sqrt(det)) / (2. * a); - Vx = radius_g + k * Vx; - Vy *= k; - Vz *= k; - - // Calculation of longitude and latitude. - lam = atan2(Vy, Vx); - phi = atan(Vz * cos(lam) / Vx); - phi = atan(radius_p_inv2 * tan(phi)); - - // To degs - lat = phi * 57.29578; - lon = lam * 57.29578; - - // Shift longitudes back to reference 0 - lon += lon_0; - if (lon < -180) - lon = lon + 360; - if (lon > 180) - lon = lon - 360; - - return 0; - } -}; \ No newline at end of file diff --git a/src-core/common/projection/geos.h b/src-core/common/projection/geos.h deleted file mode 100644 index 1c845ab16..000000000 --- a/src-core/common/projection/geos.h +++ /dev/null @@ -1,44 +0,0 @@ -#pragma once - -/* -Implementation of a standard GEOS projection, adapted from libproj. -Some variables are hardcoded for the intended usecase, making some -degree of tuning unecessary. -Uses the WGS84 ellipsoid. -*/ -namespace projection -{ - class GEOSProjection - { - private: - double h; - double radius_p; - double radius_p2; - double radius_p_inv2; - double radius_g; - double radius_g_1; - double C; - int flip_axis; - - double phi0; - double a; - double es; - double one_es; - - double lon_0; - - public: - GEOSProjection() - { - init(0, 0); - } - GEOSProjection(double height, double longitude, bool sweep_x = false) - { - init(height, longitude, sweep_x); - } - - int init(double height, double longitude, bool sweep_x = false); // return value of 1 => Error - int forward(double lon, double lat, double &x, double &y); // return value of 1 => Error - int inverse(double x, double y, double &lon, double &lat); // return value of 1 => Error - }; -}; \ No newline at end of file diff --git a/src-core/common/projection/leo_projection.cpp b/src-core/common/projection/leo_projection.cpp deleted file mode 100644 index 26e2a204b..000000000 --- a/src-core/common/projection/leo_projection.cpp +++ /dev/null @@ -1,144 +0,0 @@ -#include "leo_projection.h" - -#include -#include "libs/predict/predict.h" -#include -#include "logger.h" - -namespace projection -{ - void LEOScanProjector::initCurvatureTable() - { - const float EARTH_RADIUS = 6371.0f; - float satellite_orbit_radius = EARTH_RADIUS + settings.correction_height; // Compute the satellite's orbit radius - corrected_width = round(settings.correction_swath / settings.correction_res); // Compute the output image size, or number of samples from the imager - float satellite_view_angle = settings.correction_swath / EARTH_RADIUS; // Compute the satellite's view angle - float edge_angle = -atanf(EARTH_RADIUS * sinf(satellite_view_angle / 2) / ((cosf(satellite_view_angle / 2)) * EARTH_RADIUS - satellite_orbit_radius)); // Max angle relative to the satellite - - curvature_correction_factors_inv.reserve(settings.image_width); - - // Generate them - for (int i = 0; i < corrected_width; i++) - { - float angle = ((float(i) / float(corrected_width)) - 0.5f) * satellite_view_angle; // Get the satellite's angle - float satellite_angle = -atanf(EARTH_RADIUS * sinf(angle) / ((cosf(angle)) * EARTH_RADIUS - satellite_orbit_radius)); // Convert to an angle relative to earth - float f = settings.image_width * ((satellite_angle / edge_angle + 1.0f) / 2.0f); // Convert that to a pixel from the original image - curvature_correction_factors_fwd.push_back(f); - curvature_correction_factors_inv[int(f)] = i; - } - } - - void LEOScanProjector::generateProjections() - { - // Setup SGP4 model - predict_orbital_elements_t *satellite_object = predict_parse_tle(settings.sat_tle.line1.c_str(), - settings.sat_tle.line2.c_str()); - predict_position satellite_orbit; - - // Projection we're gonna work with - projection::TPERSProjection pj; - - // Needed to compute Az - std::function(float, float, int, int)> toSatCoords = [&pj](float lat, float lon, int map_height, int map_width) -> std::pair - { - double x, y; - pj.forward(lon, lat, x, y); - - if (fabs(x) > 1e10f || fabs(y) > 1e10f) - return {-1, -1}; - - float hscale = 4.0; - float vscale = 4.0; - int image_x = x * hscale * (map_width / 2.0); - int image_y = y * vscale * (map_height / 2.0); - - image_x += map_width / 2.0; - image_y += map_height / 2.0; - - return {image_x, (map_height - 1) - image_y}; - }; - - for (int currentScan = 0; currentScan < (int)settings.utc_timestamps.size(); currentScan++) - { - double currentTimestamp = settings.utc_timestamps[currentScan] + settings.time_offset; - - // Get Julian time of the scan, with full accuracy and calculate the satellite's - // position at the time - predict_julian_date_t currentJulianTime = predict_to_julian_double(currentTimestamp); - predict_orbit(satellite_object, &satellite_orbit, currentJulianTime); - poss.push_back(satellite_orbit); - // Calculate Az to use from the satellite's motion vector - float az = 0; - { - // Get point of view from the sat at the time - pj.init(satellite_orbit.altitude * 1000, satellite_orbit.longitude * 57.29578, satellite_orbit.latitude * 57.29578, 0, 0); - - predict_position satellite_pos1; - predict_position satellite_pos2; - predict_orbit(satellite_object, &satellite_pos1, predict_to_julian_double(currentTimestamp - 200)); - predict_orbit(satellite_object, &satellite_pos2, predict_to_julian_double(currentTimestamp + 200)); - - std::pair geo_cc1 = toSatCoords(satellite_pos1.latitude * 57.29578, satellite_pos1.longitude * 57.29578, 200, 200); - std::pair geo_cc2 = toSatCoords(satellite_pos2.latitude * 57.29578, satellite_pos2.longitude * 57.29578, 200, 200); - - // This returns the angle from the vector - az = atan(((geo_cc1.second - geo_cc2.second)) / (geo_cc1.first - geo_cc2.first)) * 57.29578; - } - - bool invertOffset = az > 0; - - az -= 90; - - // If any Az offset is required - // This has to be relative to the sat vector so... - // We swap it out when required - if (invertOffset) - az -= settings.az_offset; - else - az += settings.az_offset; - - // Get real point of view aligned with the sat's vector - pj.init(satellite_orbit.altitude * 1000, satellite_orbit.longitude * 57.29578, satellite_orbit.latitude * 57.29578, settings.tilt_offset, az); - - projs.push_back(pj); // Save that projection for later use - sat_footprints.push_back(satellite_orbit.footprint); // Save footprint - } - } - - LEOScanProjector::LEOScanProjector(LEOScanProjectorSettings settings) : settings(settings) - { - logger->info("Include curvature table..."); - initCurvatureTable(); - logger->info("Generate projection..."); - generateProjections(); - } - - int LEOScanProjector::inverse(int img_x, int img_y, double &lat, double &lon, bool correct) - { - // Check we're in bounds - if (img_y > (int)projs.size() || img_x >= settings.image_width) - return 1; - - // Get what we're gonna use - projection::TPERSProjection &pj = projs[img_y]; - double &footprint = sat_footprints[img_y]; - - // Get the width and pixel to use, depending on if we have to do correction or not - double corr_x = correct ? curvature_correction_factors_inv[img_x] : img_x; - double width = correct ? corrected_width : settings.image_width; - - // Scale to the projected area - double proj_x = settings.invert_scan ? ((width - 1) - corr_x) : corr_x; - proj_x -= width / 2.0; - proj_x += settings.proj_offset; - double pjx = proj_x / (settings.proj_scale * (width / 2.0)); - - // The instrument has a fixed FOV, so its actual footprint varies with altitude / position - // We hence scale the input values, relative to a scan angle to actual match that independently - // of any variations. - // This has to be done as the projection is not to scale - pjx *= settings.instrument_swath / footprint; - - return pj.inverse(pjx, 0, lon, lat); // Forward return value - } -}; \ No newline at end of file diff --git a/src-core/common/projection/leo_projection.h b/src-core/common/projection/leo_projection.h deleted file mode 100644 index 11309a329..000000000 --- a/src-core/common/projection/leo_projection.h +++ /dev/null @@ -1,79 +0,0 @@ -#pragma once - -#include -#include "tle.h" -#include "tpers.h" -#include "libs/predict/predict.h" - -/* -Code to reference a decoded image (or similar data) from a LEO satellite to Lat / Lon coordinates. -It works by calculating the satellite's position at the time of each scan, and generating a projection -for each of them. That projection references any point visible from the satellite's position to a -known Lat / Lon. -Doing this for every scan, with some curvature correction to match the image to the projection allows -referencing any given point. -A look-up-table is generated in the constructor to speed up later processing. - -You may also notice there is no forward function (eg, Lat / Lon to x/y on the image). That's because -the easiest way to do it would be getting the closest point in a generated LUT, but doing this -efficiently gets complicated quickly... - -PS : I am not sure in any way that this is a good implementation, there may be very obvious mistakes in there... - -PS #2 : Currently the projection is aligned to the satellite's vector by computing previous and later positions -of the satellite. It's approximate and may cause issues later but I have not found a proper way to get velocity -vectors relative to the satellite itself from the ECI values. At least, not in a form I can use for this purpose. - -PS #3 : Currently this only handles a single timestamp per scanline, which is usually most instruments work. However, -some don't like MODIS, IASI, etc. They instead send a timestamp per "IFOV" (Group of samples), such as 64x64 at once -for IASI. Others like MERSI send a full scanline of 40 pixels at once... It would be good to handle those cases properly, -utilizing each individual timestamp and properly simulating the group of detector etc... -I guess this calls for writing variants of this code later on. -*/ -namespace projection -{ - struct LEOScanProjectorSettings - { - double proj_offset; // Projection pixel offset, horizontal - int correction_swath; // Curvature correction to match the titled perspective projection. Not real instrument swath in most cases - double correction_res; // Instrument resolution for curvature correction - float correction_height; // Satellite height for curvature correction. Probably what to tune in most cases - double instrument_swath; // For coverage computation. The instrument scans at an angle relative to the entire satellite FOV - double proj_scale; // Projection scalling compared to the full satellite footprint - double az_offset; // Azimuth offset, if the instrument is angled relative the sat's motion vector - double tilt_offset; // Instrument tilt if it's not pointing at NADIR. Usually 0 - double time_offset; // Timestamp offset relative to the provided timestamps - int image_width; // Input image width - bool invert_scan; // Invert the scan direction relative to the projection - tle::TLE sat_tle; // Satellite TLEs - std::vector utc_timestamps; // Timestamps. Must match each scanline of the image you will be working with - }; - - class LEOScanProjector - { - private: - // Settings - const LEOScanProjectorSettings settings; - - // Luts and values used for referencing each line - std::vector projs; - std::vector sat_footprints; - std::vector curvature_correction_factors_fwd; // Corrected input, original image output - std::vector curvature_correction_factors_inv; // Uncorrected input, corrected output - int corrected_width; - - // Global lat / lon LUT - std::vector>> latlon_lut; - - // Internal functions - void initCurvatureTable(); // Init curvature correction table - void generateProjections(); - - public: - std::vector poss; - LEOScanProjector(LEOScanProjectorSettings settings); - - int inverse(int img_x, int img_y, double &lat, double &lon, bool correct = true); // Transform image coordinates to lat / lon. Return 1 if there was an error - //void inverse(int img_x, int img_y, double &lat, double &lon); // Transform image coordinates to lat / lon. Calls up a LUT to be faster - }; -}; \ No newline at end of file diff --git a/src-core/common/projection/proj_file.cpp b/src-core/common/projection/proj_file.cpp deleted file mode 100644 index a7163fe1c..000000000 --- a/src-core/common/projection/proj_file.cpp +++ /dev/null @@ -1,220 +0,0 @@ -#include "proj_file.h" -#include -#include "common/utils.h" -#include "tle.h" -#include -#include -#include "logger.h" - -namespace projection -{ - namespace proj_file - { - void writeReferenceFile(GeodeticReferenceFile &geofile, std::string output_file) - { - std::ofstream output_stream(output_file, std::ios::binary); - - if (geofile.file_type == GEO_TYPE) - { - GEO_GeodeticReferenceFile &gsofile = (GEO_GeodeticReferenceFile &)geofile; - - // Write main header - output_stream.write((char *)&gsofile.file_type, 1); - output_stream.write((char *)&gsofile.utc_timestamp_seconds, 8); - - // Write GEO header - output_stream.write((char *)&gsofile.norad, 4); - output_stream.write((char *)&gsofile.position_longitude, 8); - output_stream.write((char *)&gsofile.position_height, 8); - output_stream.write((char *)&gsofile.projection_type, 1); - output_stream.write((char *)&gsofile.image_width, 4); - output_stream.write((char *)&gsofile.image_height, 4); - output_stream.write((char *)&gsofile.horizontal_scale, 8); - output_stream.write((char *)&gsofile.vertical_scale, 8); - output_stream.write((char *)&gsofile.horizontal_offset, 8); - output_stream.write((char *)&gsofile.vertical_offset, 8); - output_stream.write((char *)&gsofile.proj_sweep_x, 1); - } - else if (geofile.file_type == LEO_TYPE) - { - LEO_GeodeticReferenceFile &leofile = (LEO_GeodeticReferenceFile &)geofile; - - // Write main header - output_stream.write((char *)&leofile.file_type, 1); - output_stream.write((char *)&leofile.utc_timestamp_seconds, 8); - - // Compute header size and TLE sizes - leofile.tle_line1_length = leofile.tle_line1_data.size(); - leofile.tle_line2_length = leofile.tle_line2_data.size(); - - // Write LEO header - output_stream.write((char *)&leofile.norad, 4); - output_stream.write((char *)&leofile.tle_line1_length, 2); - output_stream.write((char *)leofile.tle_line1_data.data(), leofile.tle_line1_length); - output_stream.write((char *)&leofile.tle_line2_length, 2); - output_stream.write((char *)leofile.tle_line2_data.data(), leofile.tle_line2_length); - - output_stream.write((char *)&leofile.proj_offset, 8); - output_stream.write((char *)&leofile.correction_swath, 4); - output_stream.write((char *)&leofile.correction_res, 8); - output_stream.write((char *)&leofile.correction_height, 8); - output_stream.write((char *)&leofile.instrument_swath, 8); - output_stream.write((char *)&leofile.proj_scale, 8); - output_stream.write((char *)&leofile.az_offset, 8); - output_stream.write((char *)&leofile.tilt_offset, 8); - output_stream.write((char *)&leofile.time_offset, 8); - output_stream.write((char *)&leofile.image_width, 4); - output_stream.write((char *)&leofile.invert_scan, 1); - - leofile.timestamp_count = leofile.utc_timestamps.size(); - output_stream.write((char *)&leofile.timestamp_count, 8); - for (int i = 0; i < (int)leofile.timestamp_count; i++) - output_stream.write((char *)&leofile.utc_timestamps[i], 8); - } - - output_stream.close(); - } - - std::shared_ptr readReferenceFile(std::string input_file) - { - std::ifstream input_stream(input_file, std::ios::binary); - - uint8_t file_type; - input_stream.read((char *)&file_type, 1); - if (file_type == GEO_TYPE) - { - GEO_GeodeticReferenceFile geofile; - - // Main header - input_stream.read((char *)&geofile.utc_timestamp_seconds, 8); - - // Read GEO header - input_stream.read((char *)&geofile.norad, 4); - input_stream.read((char *)&geofile.position_longitude, 8); - input_stream.read((char *)&geofile.position_height, 8); - input_stream.read((char *)&geofile.projection_type, 1); - input_stream.read((char *)&geofile.image_width, 4); - input_stream.read((char *)&geofile.image_height, 4); - input_stream.read((char *)&geofile.horizontal_scale, 8); - input_stream.read((char *)&geofile.vertical_scale, 8); - input_stream.read((char *)&geofile.horizontal_offset, 8); - input_stream.read((char *)&geofile.vertical_offset, 8); - input_stream.read((char *)&geofile.proj_sweep_x, 1); - input_stream.close(); - return std::make_shared(geofile); - } - else if (file_type == LEO_TYPE) - { - LEO_GeodeticReferenceFile leofile; - - // Main header - input_stream.read((char *)&leofile.utc_timestamp_seconds, 8); - - // Read LEO header - input_stream.read((char *)&leofile.norad, 4); - input_stream.read((char *)&leofile.tle_line1_length, 2); - { - char *str = new char[leofile.tle_line1_length]; - input_stream.read(str, leofile.tle_line1_length); - leofile.tle_line1_data = std::string(str); - delete[] str; - } - input_stream.read((char *)&leofile.tle_line2_length, 2); - { - char *str = new char[leofile.tle_line2_length]; - input_stream.read(str, leofile.tle_line2_length); - leofile.tle_line2_data = std::string(str); - delete[] str; - } - - input_stream.read((char *)&leofile.proj_offset, 8); - input_stream.read((char *)&leofile.correction_swath, 4); - input_stream.read((char *)&leofile.correction_res, 8); - input_stream.read((char *)&leofile.correction_height, 8); - input_stream.read((char *)&leofile.instrument_swath, 8); - input_stream.read((char *)&leofile.proj_scale, 8); - input_stream.read((char *)&leofile.az_offset, 8); - input_stream.read((char *)&leofile.tilt_offset, 8); - input_stream.read((char *)&leofile.time_offset, 8); - input_stream.read((char *)&leofile.image_width, 4); - input_stream.read((char *)&leofile.invert_scan, 1); - input_stream.read((char *)&leofile.timestamp_count, 8); - for (int i = 0; i < (int)leofile.timestamp_count; i++) - { - double timestamp; - input_stream.read((char *)×tamp, 8); - leofile.utc_timestamps.push_back(timestamp); - } - - input_stream.close(); - return std::make_shared(leofile); - } - else - { - // Unknown - return std::make_shared(); - } - } - - LEO_GeodeticReferenceFile leoRefFileFromProjector(int norad, LEOScanProjectorSettings projector_settings) - { - projection::proj_file::LEO_GeodeticReferenceFile geofile; - - double rough_timestamp = average_common(projector_settings.utc_timestamps.begin(), projector_settings.utc_timestamps.end()); - - geofile.utc_timestamp_seconds = floor(rough_timestamp); - geofile.norad = norad; - geofile.tle_line1_data = projector_settings.sat_tle.line1; - geofile.tle_line2_data = projector_settings.sat_tle.line2; - - geofile.proj_offset = projector_settings.proj_offset; - geofile.correction_swath = projector_settings.correction_swath; - geofile.correction_res = projector_settings.correction_res; - geofile.correction_height = projector_settings.correction_height; - geofile.instrument_swath = projector_settings.instrument_swath; - geofile.proj_scale = projector_settings.proj_scale; - geofile.az_offset = projector_settings.az_offset; - geofile.tilt_offset = projector_settings.tilt_offset; - geofile.time_offset = projector_settings.time_offset; - geofile.image_width = projector_settings.image_width; - geofile.invert_scan = projector_settings.invert_scan; - geofile.utc_timestamps = projector_settings.utc_timestamps; - - return geofile; - } - - LEOScanProjectorSettings leoProjectionRefFile(LEO_GeodeticReferenceFile geofile) - { - LEOScanProjectorSettings projector_settings = { - (double)geofile.proj_offset, - (int)geofile.correction_swath, - (double)geofile.correction_res, - (float)geofile.correction_height, - (double)geofile.instrument_swath, - (double)geofile.proj_scale, - (double)geofile.az_offset, - (double)geofile.tilt_offset, - (double)geofile.time_offset, - (int)geofile.image_width, - (bool)geofile.invert_scan, - {(int)geofile.norad, "UNKNOWN", geofile.tle_line1_data, geofile.tle_line2_data}, - geofile.utc_timestamps}; - - return projector_settings; - } - - GEOProjector geoProjectionRefFile(GEO_GeodeticReferenceFile geofile) - { - GEOProjector projector(geofile.position_longitude, - geofile.position_height, - geofile.image_width, - geofile.image_height, - geofile.horizontal_scale, - geofile.vertical_scale, - geofile.horizontal_offset, - geofile.horizontal_offset, - geofile.proj_sweep_x); - return projector; - } - }; -}; \ No newline at end of file diff --git a/src-core/common/projection/proj_file.h b/src-core/common/projection/proj_file.h deleted file mode 100644 index 7fa771881..000000000 --- a/src-core/common/projection/proj_file.h +++ /dev/null @@ -1,114 +0,0 @@ -#pragma once - -#include -#include -#include -#include "leo_projection.h" -#include "geo_projection.h" -#include - -namespace projection -{ - namespace proj_file - { - enum GeoReferenceFileTypes - { - GEO_TYPE = 1, - LEO_TYPE = 2, - }; - - /* - Class used to write & read projection files used in SatDump for geo-referencing - They contain the following data : - - Header : - - File type marker : - 1 = GEO Sat Data - 2 = LEO Sat Data - Different data types (GEO, LEO, etc) require different data types, - indicated by that marker. More may come later. - - - 64-bit UTC seconds timestamps : - The timestamp is meant to give a rough time of acquisition. It - can either be an average for LEO satellites, or the time of - acquisition for GEO satellites. - */ - struct GeodeticReferenceFile - { - uint8_t file_type; // File type - uint64_t utc_timestamp_seconds; // Timestamp. Depending on the type it may be approximate or not - }; - - /* - GEO Type 1 file. - */ - struct GEO_GeodeticReferenceFile : public GeodeticReferenceFile - { - // Type 1 - uint32_t norad; // NORAD ID of the GEO Satellite - double position_longitude; // Longitude of the satellite at the time of aquisition - double position_height; // Orbit height at the time of acquisition - uint8_t projection_type; // 0 is full disk, 1 is cropped, etc - uint32_t image_width; // Width of the image - uint32_t image_height; // Height of the image - double horizontal_scale; // Horizontal scaling - double vertical_scale; // Vertical scaling - double horizontal_offset; // Horizontal offset - double vertical_offset; // Vertical offset - bool proj_sweep_x; // Scan axis - - GEO_GeodeticReferenceFile() - { - projection_type = 0; - file_type = GEO_TYPE; - } - }; - - /* - LEO Type 2 file. - */ - struct LEO_GeodeticReferenceFile : public GeodeticReferenceFile - { - // Type 2 - uint32_t norad; // NORAD ID of the LEO Satellite - uint16_t tle_line1_length; // Size of TLE Line 1 - std::string tle_line1_data; // TLE Line 1 data - uint16_t tle_line2_length; // Size of TLE Line 1 - std::string tle_line2_data; // TLE Line 1 data - - enum LEOProjectionType - { - SINGLE_SCANLINE, - }; - - uint8_t projection_type; // Projection type - - double proj_offset; - uint32_t correction_swath; - double correction_res; - double correction_height; - double instrument_swath; - double proj_scale; - double az_offset; - double tilt_offset; - double time_offset; - uint32_t image_width; - bool invert_scan; - uint64_t timestamp_count; - std::vector utc_timestamps; // Timestamps. Must match each scanline of the image you will be working with - - LEO_GeodeticReferenceFile() - { - projection_type = SINGLE_SCANLINE; - file_type = LEO_TYPE; - } - }; - - // Functions - void writeReferenceFile(GeodeticReferenceFile &geofile, std::string output_file); - std::shared_ptr readReferenceFile(std::string input_file); - LEO_GeodeticReferenceFile leoRefFileFromProjector(int norad, LEOScanProjectorSettings projector_settings); - LEOScanProjectorSettings leoProjectionRefFile(LEO_GeodeticReferenceFile geofile); - GEOProjector geoProjectionRefFile(GEO_GeodeticReferenceFile geofile); - }; -}; \ No newline at end of file diff --git a/src-core/common/projection/satellite_reprojector.cpp b/src-core/common/projection/satellite_reprojector.cpp deleted file mode 100644 index 9beb0f996..000000000 --- a/src-core/common/projection/satellite_reprojector.cpp +++ /dev/null @@ -1,176 +0,0 @@ -#include "satellite_reprojector.h" -#include "logger.h" -#include "resources.h" -#include "common/map/map_drawer.h" - -namespace projection -{ - void reprojectLEOtoProj(cimg_library::CImg image, - projection::LEOScanProjector &projector, - cimg_library::CImg &projected_image, - int channels, - std::function(float, float, int, int)> projectionFunction - - ) - { - // Reproject - for (int currentScan = 0; currentScan < (int)image.height(); currentScan++) - { - // Now compute each pixel's lat / lon and plot it - for (double px = 0; px < image.width() - 1; px++) - { - double lat1, lon1, lat2, lon2; - int ret1 = projector.inverse(px, currentScan, lat1, lon1); - int ret2 = projector.inverse(px + 1, currentScan, lat2, lon2); - - if (ret1 || ret2) - continue; - - std::pair map_cc1 = projectionFunction(lat1, lon1, projected_image.height(), projected_image.width()); - std::pair map_cc2 = projectionFunction(lat2, lon2, projected_image.height(), projected_image.width()); - - unsigned char color[3] = {0, 0, 0}; - if (channels == 3) - { - color[0] = image[image.width() * image.height() * 0 + currentScan * image.width() + int(px)] >> 8; - color[1] = image[image.width() * image.height() * 1 + currentScan * image.width() + int(px)] >> 8; - color[2] = image[image.width() * image.height() * 2 + currentScan * image.width() + int(px)] >> 8; - } - else - { - color[0] = image[currentScan * image.width() + int(px)] >> 8; - color[1] = image[currentScan * image.width() + int(px)] >> 8; - color[2] = image[currentScan * image.width() + int(px)] >> 8; - } - - if (color[0] == 0 && color[1] == 0 && color[2] == 0) // Skip Black - continue; - - // This seems to glitch out sometimes... Need to check - if (abs(map_cc1.first - map_cc2.first) < 50 && abs(map_cc1.second - map_cc2.second) < 50) - { - double circle_radius = sqrt(pow(int(map_cc1.first - map_cc2.first), 2) + pow(int(map_cc1.second - map_cc2.second), 2)); - projected_image.draw_circle(map_cc1.first, map_cc1.second, ceil(circle_radius), color, 0.4); - } - - //if (abs(map_cc1.first - map_cc2.first) < 20 && abs(map_cc1.second - map_cc2.second) < 20) - // projected_image.draw_point(map_cc1.first, map_cc1.second, color); - - //if (map_cc1.first == map_cc2.first && map_cc1.second == map_cc2.second) - // projected_image.draw_point(map_cc1.first, map_cc1.second, color); - //else if (abs(map_cc1.first - map_cc2.first) < 20 && abs(map_cc1.second - map_cc2.second) < 20) - // projected_image.draw_rectangle(map_cc1.first, map_cc1.second, map_cc2.first, map_cc2.second, color); - } - - //logger->info(std::to_string(currentScan)); - - //logger->critical(std::to_string(currentScan)); - } - } - - void reprojectGEOtoProj(cimg_library::CImg image, - projection::GEOProjector &projector, - cimg_library::CImg &projected_image, - int channels, - std::function(float, float, int, int)> projectionFunction) - { - for (double lat = -90; lat < 90; lat += 0.01) - { - for (double lon = -180; lon < 180; lon += 0.01) - { - int x, y; - if (projector.forward(lon, lat, x, y)) - continue; - - std::pair map_cc1 = projectionFunction(lat, lon, projected_image.height(), projected_image.width()); - - unsigned char color[3]; - if (channels == 3) - { - color[0] = image[image.width() * image.height() * 0 + y * image.width() + int(x)] >> 8; - color[1] = image[image.width() * image.height() * 1 + y * image.width() + int(x)] >> 8; - color[2] = image[image.width() * image.height() * 2 + y * image.width() + int(x)] >> 8; - } - else - { - color[0] = image[y * image.width() + int(x)] >> 8; - color[1] = image[y * image.width() + int(x)] >> 8; - color[2] = image[y * image.width() + int(x)] >> 8; - } - - if (color[0] == 0 && color[1] == 0 && color[2] == 0) // Skip Black - continue; - - if (color[0] >= 253 && color[1] >= 253 && color[2] >= 253) // Skip Full white, as it's usually filler on GEO (eg, xRIT) - continue; - - //logger->info(std::to_string(color[0]) + " " + std::to_string(color[1]) + " " + std::to_string(color[2])); - - projected_image.draw_point(map_cc1.first, map_cc1.second, color); - } - logger->info(lat); - } - } - - void projectEQUIToproj(cimg_library::CImg image, cimg_library::CImg &projected_image, int channels, std::function(float, float, int, int)> toMapCoords) - { - for (double lat = -90; lat < 90; lat += 0.01) - { - for (double lon = -180; lon < 180; lon += 0.01) - { - int x = (lon / 360.0f) * image.width() + (image.width() / 2); - int y = image.height() - ((90.0f + lat) / 180.0f) * image.height(); - - if (x >= image.width() || y >= image.height()) - continue; - - std::pair map_cc1 = toMapCoords(lat, lon, projected_image.height(), projected_image.width()); - - unsigned char color[3]; - if (channels == 3) - { - color[0] = image[image.width() * image.height() * 0 + y * image.width() + int(x)] >> 8; - color[1] = image[image.width() * image.height() * 1 + y * image.width() + int(x)] >> 8; - color[2] = image[image.width() * image.height() * 2 + y * image.width() + int(x)] >> 8; - } - else - { - color[0] = image[y * image.width() + int(x)] >> 8; - color[1] = image[y * image.width() + int(x)] >> 8; - color[2] = image[y * image.width() + int(x)] >> 8; - } - - if (color[0] == 0 && color[1] == 0 && color[2] == 0) // Skip Black - continue; - - projected_image.draw_point(map_cc1.first, map_cc1.second, color); - } - logger->info(lat); - } - } - - cimg_library::CImg projectLEOToEquirectangularMapped(cimg_library::CImg image, - projection::LEOScanProjector &projector, - int output_width, - int output_height, - int channels, - cimg_library::CImg projected_image, - std::function(float, float, int, int)> toMapCoords - - ) - { - // Output mapped data - if (projected_image.width() == 1 && projected_image.height() == 1) - projected_image = cimg_library::CImg(output_width, output_height, 1, 3, 0); - - reprojectLEOtoProj(image, projector, projected_image, channels, toMapCoords); - - unsigned char color[3] = {0, 255, 0}; - map::drawProjectedMapShapefile({resources::getResourcePath("maps/ne_10m_admin_0_countries.shp")}, - projected_image, - color, - toMapCoords); - - return projected_image; - } -}; \ No newline at end of file diff --git a/src-core/common/projection/satellite_reprojector.h b/src-core/common/projection/satellite_reprojector.h deleted file mode 100644 index d2183c258..000000000 --- a/src-core/common/projection/satellite_reprojector.h +++ /dev/null @@ -1,55 +0,0 @@ -#pragma once -#define cimg_use_png -#define cimg_display 0 -#include "CImg.h" -#include -#include "leo_projection.h" -#include "geo_projection.h" - -/* -Implementation of a function capable of plotting LEO satellite image (or anything that -can be done with a LEOScanProjector) to an equirectangular projection for easier viewing. -This has some known defects right now such as leaving gaps in the image. It's something I -will have to fix at some point... Still thinking about the best approach :-) -*/ -namespace projection -{ - // Reproject LEO imagery - void reprojectLEOtoProj( - cimg_library::CImg image, // Input image to project - projection::LEOScanProjector &projector, // LEO Projector - cimg_library::CImg &projected_image, // Optional input image - int channels, // Channels - std::function(float, float, int, int)> projectionFunction // Optional projection function, default is equirectangular - ); - - // Reproject GEO imagery - void reprojectGEOtoProj( - cimg_library::CImg image, // Input image to project - projection::GEOProjector &projector, // GEO Projector - cimg_library::CImg &projected_image, // Optional input image - int channels, // Channels - std::function(float, float, int, int)> projectionFunction // Optional projection function, default is equirectangular - ); - - void projectEQUIToproj(cimg_library::CImg image, - cimg_library::CImg &projected_image, - int channels, - std::function(float, float, int, int)> toMapCoords); - - // Reproject LEO imagery to an equirectangular projection - cimg_library::CImg projectLEOToEquirectangularMapped( - cimg_library::CImg image, // Input image to project - projection::LEOScanProjector &projector, // LEO Projector - int output_width, // Output map width - int output_height, // Output map height - int channels = 1, // Channel count - cimg_library::CImg projected_image = cimg_library::CImg(1, 1, 1, 1, 0), // Optional input image - std::function(float, float, int, int)> toMapCoords = [](float lat, float lon, int map_height, int map_width) -> std::pair - { - int imageLat = map_height - ((90.0f + lat) / 180.0f) * map_height; - int imageLon = (lon / 360.0f) * map_width + (map_width / 2); - return {imageLon, imageLat}; - } // Optional projection function, default is equirectangular - ); -}; \ No newline at end of file diff --git a/src-core/common/projection/stereo.cpp b/src-core/common/projection/stereo.cpp deleted file mode 100644 index daa60a432..000000000 --- a/src-core/common/projection/stereo.cpp +++ /dev/null @@ -1,280 +0,0 @@ -#include "stereo.h" -#include -#ifndef M_PI_2 -#define M_PI_2 1.57079632679489661923 /* pi/2 */ -#endif -#ifndef M_PI_4 -#define M_PI_4 0.78539816339744830962 /* pi/4 */ -#endif - -/* -** This file was adapted and simplified from libproj, and the below -** notice kept as credits. -** -** libproj -- library of cartographic projections -** -** Copyright (c) 2004 Gerald I. Evenden -** Copyright (c) 2012 Martin Raspaud -** -** See also (section 4.4.3.2): -** https://www.cgms-info.org/documents/pdf_cgms_03.pdf -** -** Permission is hereby granted, free of charge, to any person obtaining -** a copy of this software and associated documentation files (the -** "Software"), to deal in the Software without restriction, including -** without limitation the rights to use, copy, modify, merge, publish, -** distribute, sublicense, and/or sell copies of the Software, and to -** permit persons to whom the Software is furnished to do so, subject to -** the following conditions: -** -** The above copyright notice and this permission notice shall be -** included in all copies or substantial portions of the Software. -** -** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF -** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. -** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY -** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, -** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE -** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. -*/ - -#define EPS10 1.e-10 - -namespace projection -{ - double pj_tsfn(double phi, double sinphi, double e) - { - double cosphi = cos(phi); - return exp(e * atanh(e * sinphi)) * (sinphi > 0 ? cosphi / (1 + sinphi) : (1 - sinphi) / cosphi); - } - - static double ssfn_(double phit, double sinphi, double eccen) - { - sinphi *= eccen; - return (tan(.5 * (M_PI_2 + phit)) * pow((1. - sinphi) / (1. + sinphi), .5 * eccen)); - } - - int StereoProjection::init(double latitude, double longitude) - { - lon_0 = longitude; // The projection's longitude - - // Constants, extracted from Proj - e = 0.0818192; - phi0 = latitude * 0.01745329; - a = 6.37814e+06; - es = 0.00669438; - one_es = 0.993306; - - if (es == 0.0) - { - // Illegal - return 1; - } - - k0 = .994; - //x0 = 2000000.; - //y0 = 2000000.; - phits = M_PI_2; - lam0 = 0.; - - // Setup - double t; - - if (fabs((t = fabs(phi0)) - M_PI_2) < EPS10) - mode = phi0 < 0. ? S_POLE : N_POLE; - else - mode = t > EPS10 ? OBLIQ : EQUIT; - phits = fabs(phits); - - if (es != 0.0) - { - double X; - - switch (mode) - { - case N_POLE: - case S_POLE: - if (fabs(phits - M_PI_2) < EPS10) - akm1 = 2. * k0 / - sqrt(pow(1 + e, 1 + e) * pow(1 - e, 1 - e)); - else - { - t = sin(phits); - akm1 = cos(phits) / pj_tsfn(phits, t, e); - t *= e; - akm1 /= sqrt(1. - t * t); - } - break; - case EQUIT: - case OBLIQ: - t = sin(phi0); - X = 2. * atan(ssfn_(phi0, t, e)) - M_PI_2; - t *= e; - akm1 = 2. * k0 * cos(phi0) / sqrt(1. - t * t); - sinX1 = sin(X); - cosX1 = cos(X); - break; - } - } - else - { - switch (mode) - { - case OBLIQ: - sinX1 = sin(phi0); - cosX1 = cos(phi0); - /*-fallthrough*/ - case EQUIT: - akm1 = 2. * k0; - break; - case S_POLE: - case N_POLE: - akm1 = fabs(phits - M_PI_2) >= EPS10 ? cos(phits) / tan(M_PI_4 - .5 * phits) : 2. * k0; - break; - } - } - return 0; - } - - int StereoProjection::forward(double lon, double lat, double &x, double &y) - { - x = y = 0; // Safety - - // Shift longitudes - lon -= lon_0; - if (lon < -180) - lon = lon + 360; - if (lon > 180) - lon = lon - 360; - - // To radians - double phi = lat * 0.01745329, lam = lon * 0.01745329; - - double coslam, sinlam, sinX = 0.0, cosX = 0.0, A = 0.0, sinphi; - - coslam = cos(lam); - sinlam = sin(lam); - sinphi = sin(phi); - if (mode == OBLIQ || mode == EQUIT) - { - const double X = 2. * atan(ssfn_(phi, sinphi, e)) - M_PI_2; - sinX = sin(X); - cosX = cos(X); - } - - switch (mode) - { - case OBLIQ: - { - const double denom = cosX1 * (1. + sinX1 * sinX + cosX1 * cosX * coslam); - if (denom == 0) - { - // Illegal - return 1; - } - A = akm1 / denom; - y = A * (cosX1 * sinX - sinX1 * cosX * coslam); - x = A * cosX; - break; - } - - case EQUIT: - /* avoid zero division */ - if (1. + cosX * coslam == 0.0) - { - y = HUGE_VAL; - } - else - { - A = akm1 / (1. + cosX * coslam); - y = A * sinX; - } - x = A * cosX; - break; - - case S_POLE: - phi = -phi; - coslam = -coslam; - sinphi = -sinphi; - /*-fallthrough*/ - case N_POLE: - if (fabs(phi - M_PI_2) < 1e-15) - x = 0; - else - x = akm1 * pj_tsfn(phi, sinphi, e); - y = -x * coslam; - break; - } - - x = x * sinlam; - return 0; - } - - int StereoProjection::inverse(double x, double y, double &lon, double &lat) - { - lon = lat = 0.0; - double phi = 0, lam = 0; - - double cosphi, sinphi, tp = 0.0, phi_l = 0.0, rho, halfe = 0.0, halfpi = 0.0; - - rho = hypot(x, y); - - switch (mode) - { - case OBLIQ: - case EQUIT: - tp = 2. * atan2(rho * cosX1, akm1); - cosphi = cos(tp); - sinphi = sin(tp); - if (rho == 0.0) - phi_l = asin(cosphi * sinX1); - else - phi_l = asin(cosphi * sinX1 + (y * sinphi * cosX1 / rho)); - - tp = tan(.5 * (M_PI_2 + phi_l)); - x *= sinphi; - y = rho * cosX1 * cosphi - y * sinX1 * sinphi; - halfpi = M_PI_2; - halfe = .5 * e; - break; - case N_POLE: - y = -y; - /*-fallthrough*/ - case S_POLE: - tp = -rho / akm1; - phi_l = M_PI_2 - 2. * atan(tp); - halfpi = -M_PI_2; - halfe = -.5 * e; - break; - } - - for (int i = 8; i > 0; --i) - { - sinphi = e * sin(phi_l); - phi = 2. * atan(tp * pow((1. + sinphi) / (1. - sinphi), halfe)) - halfpi; - if (fabs(phi_l - phi) < 1.e-10) - { - if (mode == S_POLE) - phi = -phi; - lam = (x == 0. && y == 0.) ? 0. : atan2(x, y); - - // To degs - lat = phi * 57.29578; - lon = lam * 57.29578; - - // Shift longitudes back to reference 0 - lon += lon_0; - if (lon < -180) - lon = lon + 360; - if (lon > 180) - lon = lon - 360; - - return 0; - } - phi_l = phi; - } - - return 1; - } -}; \ No newline at end of file diff --git a/src-core/common/projection/stereo.h b/src-core/common/projection/stereo.h deleted file mode 100644 index cd905193b..000000000 --- a/src-core/common/projection/stereo.h +++ /dev/null @@ -1,55 +0,0 @@ -#pragma once - -/* -Implementation of a standard GEOS projection, adapted from libproj. -Some variables are hardcoded for the intended usecase, making some -degree of tuning unecessary. -Uses the WGS84 ellipsoid. -*/ -namespace projection -{ - class StereoProjection - { - private: - enum Mode - { - S_POLE = 0, - N_POLE = 1, - OBLIQ = 2, - EQUIT = 3 - }; - - double phits; - double sinX1; - double cosX1; - double akm1; - enum Mode mode; - - double e; - double phi0; - double a; - double es; - double one_es; - - double k0; - //double x0; - //double y0; - double lam0; - - double lon_0; - - public: - StereoProjection() - { - init(false, 0); - } - StereoProjection(double latitude, double longitude) - { - init(latitude, longitude); - } - - int init(double latitude, double longitude); // return value of 1 => Error - int forward(double lon, double lat, double &x, double &y); // return value of 1 => Error - int inverse(double x, double y, double &lon, double &lat); // return value of 1 => Error - }; -}; \ No newline at end of file diff --git a/src-core/common/projection/tpers.cpp b/src-core/common/projection/tpers.cpp deleted file mode 100644 index 33ccc0075..000000000 --- a/src-core/common/projection/tpers.cpp +++ /dev/null @@ -1,234 +0,0 @@ -#include "tpers.h" -#include -#ifndef M_PI_2 -#define M_PI_2 1.57079632679489661923 /* pi/2 */ -#endif - -/* -** This file was adapted and simplified from libproj, and the below -** notice kept as credits. -** -** libproj -- library of cartographic projections -** -** Permission is hereby granted, free of charge, to any person obtaining -** a copy of this software and associated documentation files (the -** "Software"), to deal in the Software without restriction, including -** without limitation the rights to use, copy, modify, merge, publish, -** distribute, sublicense, and/or sell copies of the Software, and to -** permit persons to whom the Software is furnished to do so, subject to -** the following conditions: -** -** The above copyright notice and this permission notice shall be -** included in all copies or substantial portions of the Software. -** -** THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, -** EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF -** MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. -** IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY -** CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, -** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE -** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. -*/ - -#define EPS10 1.e-10 - -namespace projection -{ - int TPERSProjection::init(double altitude, double longitude, double latitude, double tilt, double azi) - { - lon_0 = longitude; - - double omega, gamma; - omega = tilt * 0.01745329; // to rads - gamma = azi * 0.01745329; // to rads - this->tilt = 1; - cg = cos(gamma); - sg = sin(gamma); - cw = cos(omega); - sw = sin(omega); - - height = altitude; - phi0 = latitude * 0.01745329; - a = 6.37814e+06; - - if (fabs(fabs(phi0) - M_PI_2) < EPS10) - mode = phi0 < 0. ? S_POLE : N_POLE; - else if (fabs(phi0) < EPS10) - mode = EQUIT; - else - { - mode = OBLIQ; - sinph0 = sin(phi0); - cosph0 = cos(phi0); - } - - pn1 = height / a; // normalize by radius - - if (pn1 <= 0 || pn1 > 1e10) - { - // Illegal!! - return 1; - } - - p = 1. + pn1; - rp = 1. / p; - h = 1. / pn1; - pfact = (p + 1.) * h; - es = 0.; - - return 0; - } - - int TPERSProjection::forward(double lon, double lat, double &x, double &y) - { - x = y = 0; // Safety - - // Shift longitudes to use the sat's as a reference - lon -= lon_0; - if (lon < -180) - lon = lon + 360; - if (lon > 180) - lon = lon - 360; - - // To radians - double phi = lat * 0.01745329, lam = lon * 0.01745329; - - double coslam, cosphi, sinphi; - - sinphi = sin(phi); - cosphi = cos(phi); - coslam = cos(lam); - - switch (mode) - { - case OBLIQ: - y = sinph0 * sinphi + cosph0 * cosphi * coslam; - break; - case EQUIT: - y = cosphi * coslam; - break; - case S_POLE: - y = -sinphi; - break; - case N_POLE: - y = sinphi; - break; - } - - if (y < rp) - { - x = y = 2e10; // Trigger error - return 1; - } - - y = pn1 / (p - y); - x = y * cosphi * sin(lam); - - switch (mode) - { - case OBLIQ: - y *= (cosph0 * sinphi - - sinph0 * cosphi * coslam); - break; - case EQUIT: - y *= sinphi; - break; - case N_POLE: - coslam = -coslam; - /*-fallthrough*/ - case S_POLE: - y *= cosphi * coslam; - break; - } - - if (tilt) - { - double yt, ba; - - yt = y * cg + x * sg; - ba = 1. / (yt * sw * h + cw); - x = (x * cg - y * sg) * cw * ba; - y = yt * ba; - } - - return 0; - } - - int TPERSProjection::inverse(double x, double y, double &lon, double &lat) - { - lon = lat = 0.0; - - double phi = 0, lam = 0; - double rh; - - if (tilt) - { - double bm, bq, yt; - - yt = 1. / (pn1 - y * sw); - bm = pn1 * x * yt; - bq = pn1 * y * cw * yt; - x = bm * cg + bq * sg; - y = bq * cg - bm * sg; - } - - rh = hypot(x, y); - - if (fabs(rh) <= EPS10) - { - lam = 0.; - phi = phi0; - } - else - { - double cosz, sinz; - sinz = 1. - rh * rh * pfact; - - if (sinz < 0.) - { - // Illegal - lon = lat = 2e10; // Trigger error - return 1; - } - - sinz = (p - sqrt(sinz)) / (pn1 / rh + rh / pn1); - cosz = sqrt(1. - sinz * sinz); - - switch (mode) - { - case OBLIQ: - phi = asin(cosz * sinph0 + y * sinz * cosph0 / rh); - y = (cosz - sinph0 * sin(phi)) * rh; - x *= sinz * cosph0; - break; - case EQUIT: - phi = asin(y * sinz / rh); - y = cosz * rh; - x *= sinz; - break; - case N_POLE: - phi = asin(cosz); - y = -y; - break; - case S_POLE: - phi = -asin(cosz); - break; - } - - lam = atan2(x, y); - } - - // To degs - lat = phi * 57.29578; - lon = lam * 57.29578; - - // Shift longitudes back to reference 0 - lon += lon_0; - if (lon < -180) - lon = lon + 360; - if (lon > 180) - lon = lon - 360; - - return 0; - } -}; \ No newline at end of file diff --git a/src-core/common/projection/tpers.h b/src-core/common/projection/tpers.h deleted file mode 100644 index c03807922..000000000 --- a/src-core/common/projection/tpers.h +++ /dev/null @@ -1,57 +0,0 @@ -#pragma once - -/* -Implementation of a titled perspective projection, for example useful to -project data from the point of a view of a LEO satellite. -This was adapted from libproj. -Uses the WGS84 ellipsoid. -*/ -namespace projection -{ - class TPERSProjection - { - private: - enum Mode - { - N_POLE = 0, - S_POLE = 1, - EQUIT = 2, - OBLIQ = 3 - }; - - double height; - double sinph0; - double cosph0; - double p; - double rp; - double pn1; - double pfact; - double h; - double cg; - double sg; - double sw; - double cw; - enum Mode mode; - int tilt; - - double phi0; - double a; - double es; - - double lon_0; - - public: - TPERSProjection() - { - init(0, 0, 0, 0, 0); - } - TPERSProjection(double altitude, double longitude, double latitude, double tilt, double azi) - { - init(altitude, longitude, latitude, tilt, azi); - } - - int init(double altitude, double longitude, double latitude, double tilt, double azi); // return value of 1 => Error - int forward(double lon, double lat, double &x, double &y); // return value of 1 => Error - int inverse(double x, double y, double &lon, double &lat); // return value of 1 => Error - }; -}; \ No newline at end of file diff --git a/src-core/modules/eos/instruments/modis/module_eos_modis.cpp b/src-core/modules/eos/instruments/modis/module_eos_modis.cpp index b7dcf6d43..5282a3ccb 100644 --- a/src-core/modules/eos/instruments/modis/module_eos_modis.cpp +++ b/src-core/modules/eos/instruments/modis/module_eos_modis.cpp @@ -12,7 +12,7 @@ #include "common/image/earth_curvature.h" #include "modules/eos/eos.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" +#include "common/geodetic/projection/satellite_reprojector.h" #define BUFFER_SIZE 8192 @@ -365,24 +365,21 @@ namespace eos //pre_wb.resize(pre_wb.width() / 4, pre_wb.height() / 4); // Setup Projecition - projection::LEOScanProjector projector({ - 0, // Pixel offset - 1950, // Correction swath - EOS_MODIS_RES250, // Instrument res - 800, // Orbit height - 2220, // Instrument swath - 2.45, // Scale - -2.2, // Az offset - 0, // Tilt - -2.2, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 109.7, // Scan angle + -0.1, // Roll offset + 0, // Pitch offset + -2.5, // Yaw offset + -2.0, // Time offset pre_wb.width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs reader.timestamps_250 // Timestamps - }); + }; + geodetic::projection::LEOScanProjector projector(proj_settings); logger->info("Projected Channel 143 EQURAW..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(pre_wb, projector, 2048 * 4, 1024 * 4, 3); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(pre_wb, projector, 2048 * 4, 1024 * 4, 3); WRITE_IMAGE(projected_image, directory + "/MODIS-143-EQURAW-PROJ.png"); } } @@ -399,29 +396,27 @@ namespace eos EOS_ORBIT_HEIGHT, EOS_MODIS_SWATH, EOS_MODIS_RES1000); - WRITE_IMAGE(corrected23, directory + "/MODIS-29-EQU-CORRECTED.png"); + //WRITE_IMAGE(corrected23, directory + "/MODIS-29-EQU-CORRECTED.png"); // Reproject to an equirectangular proj { // Setup Projecition - projection::LEOScanProjector projector({ - 0, // Pixel offset - 1800, // Correction swath - EOS_MODIS_RES250, // Instrument res - 710, // Orbit height - 2230, // Instrument swath - 2.435, // Scale - -2.2, // Az offset - 0, // Tilt - -2.2, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 109.7, // Scan angle + -0.1, // Roll offset + 0, // Pitch offset + 2.5, // Yaw offset + -2.0, // Time offset image23.width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs reader.timestamps_1000 // Timestamps - }); + }; + geodetic::projection::LEOScanProjector projector(proj_settings); logger->info("Projected Channel 29..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image23, projector, 2048 * 4, 1024 * 4, 1); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image23, projector, 2048 * 4, 1024 * 4, 1); + //projected_image.crop(20778, 2853, 20778 + 5145, 2853 + 3573); WRITE_IMAGE(projected_image, directory + "/MODIS-29-EQU-PROJ.png"); } } diff --git a/src-core/modules/fengyun/instruments/erm/module_fengyun_erm.cpp b/src-core/modules/fengyun/instruments/erm/module_fengyun_erm.cpp index 517141e6a..49698ba36 100644 --- a/src-core/modules/fengyun/instruments/erm/module_fengyun_erm.cpp +++ b/src-core/modules/fengyun/instruments/erm/module_fengyun_erm.cpp @@ -6,8 +6,8 @@ #include "erm_reader.h" #include "common/ccsds/ccsds_1_0_1024/demuxer.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" // Return filesize size_t getFilesize(std::string filepath); @@ -100,31 +100,27 @@ namespace fengyun int norad = satData.contains("norad") > 0 ? satData["norad"].get() : 0; // Setup Projecition - projection::LEOScanProjectorSettings proj_settings = { - 6, // Pixel offset - 1500, // Correction swath - 16.0 / 4, // Instrument res - 827.0, // Orbit height - 2300, // Instrument swath - 2.4, // Scale - 1, // Az offset - 0, // Tilt - 1.5, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 102, // Scan angle + -0.0, // Roll offset + 0, // Pitch offset + -1, // Yaw offset + 2, // Time offset erm_reader.getChannel().width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs erm_reader.timestamps // Timestamps }; - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/ERM.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/ERM.georef"); } cimg_library::CImg image = erm_reader.getChannel(); logger->info("Projected Channel 1..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image, projector, 2048, 1024); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector, 2048, 1024); WRITE_IMAGE(projected_image, directory + "/ERM-1-PROJ.png"); } } diff --git a/src-core/modules/fengyun/instruments/mwhs/module_fengyun_mwhs.cpp b/src-core/modules/fengyun/instruments/mwhs/module_fengyun_mwhs.cpp index b1863e4a8..72e461c8b 100644 --- a/src-core/modules/fengyun/instruments/mwhs/module_fengyun_mwhs.cpp +++ b/src-core/modules/fengyun/instruments/mwhs/module_fengyun_mwhs.cpp @@ -7,8 +7,8 @@ #include "mwhs_reader.h" #include "common/ccsds/ccsds_1_0_1024/demuxer.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" // Return filesize size_t getFilesize(std::string filepath); @@ -120,33 +120,29 @@ namespace fengyun int norad = satData.contains("norad") > 0 ? satData["norad"].get() : 0; // Setup Projecition - projection::LEOScanProjectorSettings proj_settings = { - 8, // Pixel offset - 1800, // Correction swath - 16.0 / 4, // Instrument res - 827.0, // Orbit height - 2800, // Instrument swath - 2.55, // Scale - 0, // Az offset - 0, // Tilt - -4, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 108, // Scan angle + -3.5, // Roll offset + 0, // Pitch offset + 0, // Yaw offset + -1, // Time offset mwhs_reader.getChannel(0).width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs mwhs_reader.timestamps // Timestamps }; - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/MWHS.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/MWHS.georef"); } for (int i = 0; i < 5; i++) { cimg_library::CImg image = mwhs_reader.getChannel(i); logger->info("Projected Channel " + std::to_string(i + 1) + "..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image, projector, 2048, 1024); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector, 2048, 1024); WRITE_IMAGE(projected_image, directory + "/MWHS-" + std::to_string(i + 1) + "-PROJ.png"); } } diff --git a/src-core/modules/fengyun/instruments/mwhs2/module_fengyun_mwhs2.cpp b/src-core/modules/fengyun/instruments/mwhs2/module_fengyun_mwhs2.cpp index bb102ce0c..734c7e6dd 100644 --- a/src-core/modules/fengyun/instruments/mwhs2/module_fengyun_mwhs2.cpp +++ b/src-core/modules/fengyun/instruments/mwhs2/module_fengyun_mwhs2.cpp @@ -6,9 +6,9 @@ #include "mwhs2_reader.h" #include "common/ccsds/ccsds_1_0_1024/demuxer.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" #include "modules/fengyun/fengyun3.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" // Return filesize size_t getFilesize(std::string filepath); @@ -133,33 +133,29 @@ namespace fengyun if (mwhs_reader.lines > 0) { // Setup Projecition - projection::LEOScanProjectorSettings proj_settings = { - 4, // Pixel offset - 1550, // Correction swath - 16.0 / 4, // Instrument res - 827.0, // Orbit height - 2800, // Instrument swath - 2.65, // Scale - 0, // Az offset - 0, // Tilt - 3, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 106, // Scan angle + -0, // Roll offset + 0, // Pitch offset + -5, // Yaw offset + 2, // Time offset mwhs_reader.getChannel(0).width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs mwhs_reader.timestamps // Timestamps }; - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/MWHS-2.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/MWHS-2.georef"); } for (int i = 0; i < 15; i++) { cimg_library::CImg image = mwhs_reader.getChannel(i); logger->info("Projected Channel " + std::to_string(i + 1) + "..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image, projector, 2048, 1024); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector, 2048, 1024); WRITE_IMAGE(projected_image, directory + "/MWHS2-" + std::to_string(i + 1) + "-PROJ.png"); } } diff --git a/src-core/modules/fengyun/instruments/mwts/module_fengyun_mwts.cpp b/src-core/modules/fengyun/instruments/mwts/module_fengyun_mwts.cpp index 7f763f412..aadf421aa 100644 --- a/src-core/modules/fengyun/instruments/mwts/module_fengyun_mwts.cpp +++ b/src-core/modules/fengyun/instruments/mwts/module_fengyun_mwts.cpp @@ -6,8 +6,8 @@ #include "mwts_reader.h" #include "common/ccsds/ccsds_1_0_1024/demuxer.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" // Return filesize size_t getFilesize(std::string filepath); @@ -115,33 +115,29 @@ namespace fengyun int norad = satData.contains("norad") > 0 ? satData["norad"].get() : 0; // Setup Projecition - projection::LEOScanProjectorSettings proj_settings = { - 0, // Pixel offset - 1700, // Correction swath - 17.4 / 20, // Instrument res - 827.0, // Orbit height - 2200, // Instrument swath - 2.25, // Scale - 0, // Az offset - 0, // Tilt - 0, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 108, // Scan angle + 3, // Roll offset + 0, // Pitch offset + 0, // Yaw offset + -1, // Time offset mwts_reader.getChannel(0).width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs mwts_reader.timestamps // Timestamps }; - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/MWTS.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/MWTS.georef"); } for (int i = 0; i < 27; i++) { cimg_library::CImg image = mwts_reader.getChannel(i); logger->info("Projected Channel " + std::to_string(i + 1) + "..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image, projector, 2048 / 2, 1024 / 2); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector, 2048 / 2, 1024 / 2); WRITE_IMAGE(projected_image, directory + "/MWTS-" + std::to_string(i + 1) + "-PROJ.png"); } } diff --git a/src-core/modules/fengyun/instruments/mwts2/module_fengyun_mwts2.cpp b/src-core/modules/fengyun/instruments/mwts2/module_fengyun_mwts2.cpp index 028f207a6..e2ecb8a4b 100644 --- a/src-core/modules/fengyun/instruments/mwts2/module_fengyun_mwts2.cpp +++ b/src-core/modules/fengyun/instruments/mwts2/module_fengyun_mwts2.cpp @@ -6,8 +6,8 @@ #include "common/ccsds/ccsds_1_0_1024/demuxer.h" #include "mwts2_reader.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" // Return filesize size_t getFilesize(std::string filepath); @@ -140,33 +140,29 @@ namespace fengyun int norad = satData.contains("norad") > 0 ? satData["norad"].get() : 0; // Setup Projecition - projection::LEOScanProjectorSettings proj_settings = { - 60, // Pixel offset - 1400, // Correction swath - 17.4 / 20, // Instrument res - 827.0, // Orbit height - 2200, // Instrument swath - 2.42, // Scale - 0, // Az offset - 0, // Tilt - 0, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 98, // Scan angle + -2, // Roll offset + 0, // Pitch offset + -4, // Yaw offset + 1, // Time offset mwts_reader.getChannel(0).width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs mwts_reader.timestamps // Timestamps }; - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/MWTS-2.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/MWTS-2.georef"); } for (int i = 0; i < 16; i++) { cimg_library::CImg image = mwts_reader.getChannel(i); logger->info("Projected Channel " + std::to_string(i + 1) + "..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image, projector, 2048 / 2, 1024 / 2); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector, 2048 / 2, 1024 / 2); WRITE_IMAGE(projected_image, directory + "/MWTS2-" + std::to_string(i + 1) + "-PROJ.png"); } } diff --git a/src-core/modules/fengyun/instruments/mwts3/module_fengyun_mwts3.cpp b/src-core/modules/fengyun/instruments/mwts3/module_fengyun_mwts3.cpp index fa254ea1d..bf7b51425 100644 --- a/src-core/modules/fengyun/instruments/mwts3/module_fengyun_mwts3.cpp +++ b/src-core/modules/fengyun/instruments/mwts3/module_fengyun_mwts3.cpp @@ -6,8 +6,8 @@ #include "common/ccsds/ccsds_1_0_1024/demuxer.h" #include "mwts3_reader.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" // Return filesize size_t getFilesize(std::string filepath); @@ -144,33 +144,29 @@ namespace fengyun int norad = satData.contains("norad") > 0 ? satData["norad"].get() : 0; // Setup Projecition - projection::LEOScanProjectorSettings proj_settings = { - 48, // Pixel offset - 1600, // Correction swath - 17.4 / 20, // Instrument res - 827.0, // Orbit height - 2200, // Instrument swath - 2.10, // Scale - 0, // Az offset - 0, // Tilt - 0, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 104, // Scan angle + -1, // Roll offset + 0, // Pitch offset + -5, // Yaw offset + -1, // Time offset mwts_reader.getChannel(0).width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs mwts_reader.timestamps // Timestamps }; - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/MWTS-3.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/MWTS-3.georef"); } for (int i = 0; i < 18; i++) { cimg_library::CImg image = mwts_reader.getChannel(i); logger->info("Projected Channel " + std::to_string(i + 1) + "..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image, projector, 2048 / 2, 1024 / 2); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector, 2048 / 2, 1024 / 2); WRITE_IMAGE(projected_image, directory + "/MWTS3-" + std::to_string(i + 1) + "-PROJ.png"); } } diff --git a/src-core/modules/fengyun/instruments/virr/module_fengyun_virr.cpp b/src-core/modules/fengyun/instruments/virr/module_fengyun_virr.cpp index 4270acdc0..53e820e99 100644 --- a/src-core/modules/fengyun/instruments/virr/module_fengyun_virr.cpp +++ b/src-core/modules/fengyun/instruments/virr/module_fengyun_virr.cpp @@ -9,8 +9,8 @@ #include "common/image/earth_curvature.h" #include "modules/fengyun/fengyun3.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" #define BUFFER_SIZE 8192 @@ -384,26 +384,22 @@ namespace fengyun int norad = satData.contains("norad") > 0 ? satData["norad"].get() : 0; // Setup Projecition - projection::LEOScanProjectorSettings proj_settings = { - 3, // Pixel offset - 2100, // Correction swath - 1.1, // Instrument res - 830, // Orbit height - 2800, // Instrument swath - 2.406, // Scale - -2.7, // Az offset - 0, // Tilt - -0.1, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 110.4, // Scan angle + -0.0, // Roll offset + 0, // Pitch offset + -2.6, // Yaw offset + 0, // Time offset image1.width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs reader.timestamps // Timestamps }; - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/VIRR.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/VIRR.georef"); } cimg_library::CImg projected_image; @@ -421,7 +417,7 @@ namespace fengyun image197equ.normalize(0, std::numeric_limits::max()); } logger->info("Projected Channel 197..."); - projected_image = projection::projectLEOToEquirectangularMapped(image197equ, projector, 2048 * 4, 1024 * 4, 3); + projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image197equ, projector, 2048 * 4, 1024 * 4, 3); WRITE_IMAGE(projected_image, directory + "/VIRR-RGB-197-PROJ.png"); } @@ -429,7 +425,7 @@ namespace fengyun logger->info("Projected Channel 4..."); cimg_library::CImg tempImage4 = image4; tempImage4.equalize(1000); - projected_image = projection::projectLEOToEquirectangularMapped(tempImage4, projector, 2048 * 4, 1024 * 4, 1); + projected_image = geodetic::projection::projectLEOToEquirectangularMapped(tempImage4, projector, 2048 * 4, 1024 * 4, 1); WRITE_IMAGE(projected_image, directory + "/VIRR-4-PROJ.png"); } @@ -443,7 +439,7 @@ namespace fengyun image621.normalize(0, std::numeric_limits::max()); } logger->info("Projected channel 621..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image621, projector, 2048 * 4, 1024 * 4, 3); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image621, projector, 2048 * 4, 1024 * 4, 3); WRITE_IMAGE(projected_image, directory + "/VIRR-RGB-621-PROJ.png"); } } diff --git a/src-core/modules/jpss/instruments/atms/module_jpss_atms.cpp b/src-core/modules/jpss/instruments/atms/module_jpss_atms.cpp index a20d24d3f..3b8ea02e0 100644 --- a/src-core/modules/jpss/instruments/atms/module_jpss_atms.cpp +++ b/src-core/modules/jpss/instruments/atms/module_jpss_atms.cpp @@ -7,8 +7,8 @@ #include #include "imgui/imgui.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" // Return filesize size_t getFilesize(std::string filepath); @@ -390,26 +390,22 @@ namespace jpss int norad = satData.contains("norad") > 0 ? satData["norad"].get() : 0; // Setup Projecition - projection::LEOScanProjectorSettings proj_settings = { - 3, // Pixel offset - 1700, // Correction swath - 16.0 / 4, // Instrument res - 827.0, // Orbit height - 2200, // Instrument swath - 2.14, // Scale - -3, // Az offset - 0, // Tilt - -1.0, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 102, // Scan angle + -0.5, // Roll offset + 0, // Pitch offset + -3, // Yaw offset + 2, // Time offset image1.width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs reader.timestamps // Timestamps }; - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/ATMS.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/ATMS.georef"); } for (int i = 0; i < 22; i++) @@ -417,7 +413,7 @@ namespace jpss cimg_library::CImg image = reader.getImage(i); image.equalize(1000); logger->info("Projected channel " + std::to_string(i + 1) + "..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image, projector, 2048, 1024); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector, 2048, 1024); WRITE_IMAGE(projected_image, directory + "/ATMS-" + std::to_string(i + 1) + "-PROJ.png"); } @@ -429,7 +425,7 @@ namespace jpss } image3417.equalize(1000); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image3417, projector, 2048, 1024, 3); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image3417, projector, 2048, 1024, 3); WRITE_IMAGE(projected_image, directory + "/ATMS-RGB-3.4.17-PROJ.png"); } } diff --git a/src-core/modules/meteor/instruments/msumr/module_meteor_msumr.cpp b/src-core/modules/meteor/instruments/msumr/module_meteor_msumr.cpp index 3a439608a..396a01014 100644 --- a/src-core/modules/meteor/instruments/msumr/module_meteor_msumr.cpp +++ b/src-core/modules/meteor/instruments/msumr/module_meteor_msumr.cpp @@ -8,10 +8,10 @@ #include "common/image/earth_curvature.h" #include "modules/meteor/meteor.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" +#include "common/geodetic/projection/satellite_reprojector.h" #include "common/image/brightness_contrast.h" #include "common/image/image.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/proj_file.h" #include "common/utils.h" #define BUFFER_SIZE 8192 @@ -51,7 +51,7 @@ namespace meteor logger->info("Demultiplexing and deframing..."); time_t currentDay = time(0); - time_t dayValue = currentDay - (currentDay % 86400) - 86400 * 0; // Requires the day to be known from another source + time_t dayValue = currentDay - (currentDay % 86400) - 86400 * 1; // Requires the day to be known from another source std::vector timestamps; @@ -226,20 +226,16 @@ namespace meteor image::brightness_contrast(image321, 0.179 * 2, 0.253 * 2, 3); // Setup Projecition, tuned for 2-2 - projection::LEOScanProjectorSettings proj_settings = { - 4, // Pixel offset - 2070, // Correction swath - 1.1, // Instrument res - 830, // Orbit height - METEOR_MSUMR_SWATH, // Instrument swath - 2.46, // Scale - -2, // Az offset - 0, // Tilt - 1, // Time offset - image1.width(), // Image width - true, // Invert scan - tle::TLE(), // TLEs - timestamps // Timestamps + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 110.2, // Scan angle + -0.35, // Roll offset + 0, // Pitch offset + -2.9, // Yaw offset + 0.2, // Time offset + image1.width(), // Image width + true, // Invert scan + tle::TLE(), // TLEs + timestamps // Timestamps }; // Identify satellite, and apply per-sat settings... @@ -247,9 +243,11 @@ namespace meteor if (msumr_serial_number == 0) // METEOR-M 2, weirdly enough it has ID 0 { norad = 40069; - logger->info("Identified METEOR-M 2! Projection support is WIP!"); - proj_settings.proj_offset = -40; - proj_settings.proj_scale = 2.330; + logger->info("Identified METEOR-M 2!"); + proj_settings.scan_angle = 111.0; + proj_settings.roll_offset = 3.1; + proj_settings.yaw_offset = -3; + proj_settings.time_offset = 9; } else if (msumr_serial_number == 1) // METEOR-M 2-1... Launch failed of course... { @@ -269,19 +267,19 @@ namespace meteor // Load TLEs now proj_settings.sat_tle = tle::getTLEfromNORAD(norad); - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/MSU-MR.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/MSU-MR.georef"); } logger->info("Projected RGB 321..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image321, projector, 2048 * 4, 1024 * 4, 3); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image321, projector, 2048 * 4, 1024 * 4, 3); WRITE_IMAGE(projected_image, directory + "/MSU-MR-RGB-321-PROJ.png"); logger->info("Projected Channel 4..."); - projected_image = projection::projectLEOToEquirectangularMapped(image4, projector, 2048 * 4, 1024 * 4, 1); + projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image4, projector, 2048 * 4, 1024 * 4, 1); WRITE_IMAGE(projected_image, directory + "/MSU-MR-4-PROJ.png"); } } diff --git a/src-core/modules/metop/instruments/amsu/module_metop_amsu.cpp b/src-core/modules/metop/instruments/amsu/module_metop_amsu.cpp index fe017ada2..603e72564 100644 --- a/src-core/modules/metop/instruments/amsu/module_metop_amsu.cpp +++ b/src-core/modules/metop/instruments/amsu/module_metop_amsu.cpp @@ -8,8 +8,8 @@ #include #include "imgui/imgui.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" #define BUFFER_SIZE 8192 @@ -150,52 +150,45 @@ namespace metop // There is no "real" guarantee the A1 / A2 output will always be identical // Using the "/ 40" in instrument res slows things down but also avoids huge gaps // in the resulting image... - projection::LEOScanProjectorSettings proj_settings_a1 = { - 0, // Pixel offset - 1900, // Correction swath - 48.0 / 40, // Instrument res - 827.0, // Orbit height - 2250, // Instrument swath - 2.2, // Scale - 0, // Az offset - 0, // Tilt - 4.0, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings_a1 = { + 98, // Scan angle + 0, // Roll offset + 0, // Pitch offset + 0, // Yaw offset + 10, // Time offset a1reader.getChannel(0).width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs a1reader.timestamps // Timestamps }; - projection::LEOScanProjectorSettings proj_settings_a2 = { - 0, // Pixel offset - 1900, // Correction swath - 48.0 / 40, // Instrument res - 827.0, // Orbit height - 2250, // Instrument swath - 2.2, // Scale - 0, // Az offset - 0, // Tilt - 4.0, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings_a2 = { + 98, // Scan angle + 0, // Roll offset + 0, // Pitch offset + 0, // Yaw offset + 10, // Time offset a2reader.getChannel(0).width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs a2reader.timestamps // Timestamps }; - projection::LEOScanProjector projector_a1(proj_settings_a1); - projection::LEOScanProjector projector_a2(proj_settings_a2); + geodetic::projection::LEOScanProjector projector_a1(proj_settings_a1); + geodetic::projection::LEOScanProjector projector_a2(proj_settings_a2); { - projection::proj_file::LEO_GeodeticReferenceFile geofile_a1 = projection::proj_file::leoRefFileFromProjector(norad, proj_settings_a1); - projection::proj_file::writeReferenceFile(geofile_a1, directory + "/AMSU-A1.georef"); - projection::proj_file::LEO_GeodeticReferenceFile geofile_a2 = projection::proj_file::leoRefFileFromProjector(norad, proj_settings_a2); - projection::proj_file::writeReferenceFile(geofile_a2, directory + "/AMSU-A2.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile_a1 = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings_a1); + geodetic::projection::proj_file::writeReferenceFile(geofile_a1, directory + "/AMSU-A1.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile_a2 = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings_a2); + geodetic::projection::proj_file::writeReferenceFile(geofile_a2, directory + "/AMSU-A2.georef"); } for (int i = 0; i < 13; i++) { cimg_library::CImg image = a1reader.getChannel(i); image.equalize(1000); + image.normalize(0, 65535); logger->info("Projected channel A1 " + std::to_string(i + 3) + "..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image, projector_a1, 1024, 512); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector_a1, 1024, 512); WRITE_IMAGE(projected_image, directory + "/AMSU-A1-" + std::to_string(i + 3) + "-PROJ.png"); } @@ -203,8 +196,9 @@ namespace metop { cimg_library::CImg image = a2reader.getChannel(i); image.equalize(1000); + image.normalize(0, 65535); logger->info("Projected channel A2 " + std::to_string(i + 1) + "..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image, projector_a2, 1024, 512); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector_a2, 1024, 512); WRITE_IMAGE(projected_image, directory + "/AMSU-A2-" + std::to_string(i + 1) + "-PROJ.png"); } } diff --git a/src-core/modules/metop/instruments/ascat/module_metop_ascat.cpp b/src-core/modules/metop/instruments/ascat/module_metop_ascat.cpp index aac5b8596..8b48c2b82 100644 --- a/src-core/modules/metop/instruments/ascat/module_metop_ascat.cpp +++ b/src-core/modules/metop/instruments/ascat/module_metop_ascat.cpp @@ -6,9 +6,10 @@ #include "logger.h" #include #include "imgui/imgui.h" -//#include "common/projection/leo_to_equirect.h" -//#include "nlohmann/json_utils.h" -//#include "common/projection/ref_file.h" +#include "nlohmann/json_utils.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" +#include "common/image/earth_curvature.h" #define BUFFER_SIZE 8192 @@ -159,54 +160,36 @@ namespace metop } WRITE_IMAGE(imageAll, directory + "/ASCAT-EQU-ALL.png"); - /* // Reproject to an equirectangular proj - if (image1.height() > 0) + /*if (image2.height() > 0) { nlohmann::json satData = loadJsonFile(d_output_file_hint.substr(0, d_output_file_hint.rfind('/')) + "/sat_info.json"); int norad = satData.contains("norad") > 0 ? satData["norad"].get() : 0; //image4.equalize(1000); - int compo_size = 790; - - // Setup Projecition - projection::LEOScanProjector projector(0, // Pixel offset - 910, // Correction swath - 1, // Instrument res - 800, // Orbit height - 1800, // Instrument swath - 1.9, // Scale - -1, // Az offset - 0, // Tilt - -0, // Time offset - compo_size, // Image width - true, // Invert scan - tle::getTLEfromNORAD(norad), // TLEs - ascatreader.timestamps[1] // Timestamps - ); + // Setup Projecition + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 35, // Scan angle + 35, // Roll offset + 0, // Pitch offset + 0, // Yaw offset + -1, // Time offset + image2.width(), // Image width + true, // Invert scan + tle::getTLEfromNORAD(norad), // TLEs + ascatreader.timestamps[1] // Timestamps + }; + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::ref_file::LEO_GeodeticReferenceFile geofile = projection::ref_file::leoRefFileFromProjector(projector, - ascatreader.timestamps[1], - norad, - 2048, - reader.lines); - logger->info("GEO Reference file..."); - projection::ref_file::writeReferenceFile(geofile, directory + "/ASCAT.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/ASCAT-2.georef"); } - + logger->info("Projected channel 2..."); - cimg_library::CImg image_recomp(compo_size, ascatreader.lines[1], 1, 1, 0); - { - - image_recomp.draw_image(0, 0, 0, 0, image5); - image_recomp.draw_image(compo_size - 256, 0, 0, 0, image2); - } - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image_recomp, projector, 2048 * 8, 1024 * 8, 1); - WRITE_IMAGE(image_recomp, directory + "/ASCAT-MIDDLE.png"); - WRITE_IMAGE(projected_image, directory + "/ASCAT-MIDDLE-PROJ.png"); - } - */ + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(corr2, projector, 2048 * 4, 1024 * 4, 1); + WRITE_IMAGE(projected_image, directory + "/ASCAT-2-PROJ.png"); + }*/ } void MetOpASCATDecoderModule::drawUI(bool window) diff --git a/src-core/modules/metop/instruments/avhrr/module_metop_avhrr.cpp b/src-core/modules/metop/instruments/avhrr/module_metop_avhrr.cpp index 0a9e22dfd..ffc144277 100644 --- a/src-core/modules/metop/instruments/avhrr/module_metop_avhrr.cpp +++ b/src-core/modules/metop/instruments/avhrr/module_metop_avhrr.cpp @@ -11,10 +11,10 @@ #include "modules/metop/metop.h" #include "nlohmann/json_utils.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" +#include "common/geodetic/projection/satellite_reprojector.h" #include "common/image/brightness_contrast.h" #include "common/image/xfr.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/proj_file.h" #define BUFFER_SIZE 8192 @@ -297,30 +297,26 @@ namespace metop //image4.equalize(1000); // Setup Projecition - projection::LEOScanProjectorSettings proj_settings = { - 5, // Pixel offset - 2050, // Correction swath - 1, // Instrument res - 800, // Orbit height - METOP_AVHRR_SWATH, // Instrument swath - 2.515, // Scale - 0.4, // Az offset - 0, // Tilt + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 110.8, // Scan angle + -0.1, // Roll offset + 0, // Pitch offset + 0.2, // Yaw offset -0.3, // Time offset image4.width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs reader.timestamps // Timestamps }; - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/AVHRR.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/AVHRR.georef"); } logger->info("Projected channel 4..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image4, projector, 2048 * 4, 1024 * 4, 1); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image4, projector, 2048 * 4, 1024 * 4, 1); WRITE_IMAGE(projected_image, directory + "/AVHRR-4-PROJ.png"); cimg_library::CImg image321(2048, reader.lines, 1, 3); @@ -331,7 +327,7 @@ namespace metop } logger->info("Projected channel 321..."); - projected_image = projection::projectLEOToEquirectangularMapped(image321, projector, 2048 * 4, 1024 * 4, 3); + projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image321, projector, 2048 * 4, 1024 * 4, 3); WRITE_IMAGE(projected_image, directory + "/AVHRR-RGB-321-PROJ.png"); } } diff --git a/src-core/modules/metop/instruments/mhs/module_metop_mhs.cpp b/src-core/modules/metop/instruments/mhs/module_metop_mhs.cpp index 0ca085559..1936450c8 100644 --- a/src-core/modules/metop/instruments/mhs/module_metop_mhs.cpp +++ b/src-core/modules/metop/instruments/mhs/module_metop_mhs.cpp @@ -7,8 +7,8 @@ #include #include "imgui/imgui.h" #include "nlohmann/json_utils.h" -#include "common/projection/satellite_reprojector.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" #define BUFFER_SIZE 8192 @@ -123,33 +123,29 @@ namespace metop int norad = satData.contains("norad") > 0 ? satData["norad"].get() : 0; // Setup Projecition - projection::LEOScanProjectorSettings proj_settings = { - 2, // Pixel offset - 2070, // Correction swath - 16.0 / 4, // Instrument res - 827.0, // Orbit height - 2180, // Instrument swath - 2.18, // Scale - 0, // Az offset - 0, // Tilt - 1, // Time offset + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 100, // Scan angle + -0.7, // Roll offset + 0, // Pitch offset + 0.5, // Yaw offset + 3.5, // Time offset mhsreader.getChannel(0).width(), // Image width true, // Invert scan tle::getTLEfromNORAD(norad), // TLEs mhsreader.timestamps // Timestamps }; - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/MHS.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/MHS.georef"); } for (int i = 0; i < 5; i++) { cimg_library::CImg image = mhsreader.getChannel(i); logger->info("Projected channel " + std::to_string(i + 1) + "..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image, projector, 2048, 1024); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image, projector, 2048, 1024); WRITE_IMAGE(projected_image, directory + "/MHS-" + std::to_string(i + 1) + "-PROJ.png"); } } diff --git a/src-core/modules/noaa/instruments/avhrr/module_noaa_avhrr.cpp b/src-core/modules/noaa/instruments/avhrr/module_noaa_avhrr.cpp index 25b8b58ae..e6f2ae685 100644 --- a/src-core/modules/noaa/instruments/avhrr/module_noaa_avhrr.cpp +++ b/src-core/modules/noaa/instruments/avhrr/module_noaa_avhrr.cpp @@ -6,9 +6,9 @@ #include "imgui/imgui.h" #include "common/image/earth_curvature.h" #include "modules/noaa/noaa.h" -#include "common/projection/satellite_reprojector.h" +#include "common/geodetic/projection/satellite_reprojector.h" #include "nlohmann/json_utils.h" -#include "common/projection/proj_file.h" +#include "common/geodetic/projection/proj_file.h" #include "common/utils.h" #define BUFFER_SIZE 8192 @@ -160,23 +160,24 @@ namespace noaa int norad = 0; //28654; //satData.contains("norad") > 0 ? satData["norad"].get() : 0; image4.equalize(1000); - // Setup Projecition, based off N18 - projection::LEOScanProjectorSettings proj_settings = { - 2, // Pixel offset - 2050, // Correction swath - 1, // Instrument res - 800, // Orbit height - NOAA_AVHRR_SWATH, // Instrument swath - 2.399, // Scale - -2.9, // Az offset - 0, // Tilt - 0.3, // Time offset - image4.width(), // Image width - true, // Invert scan - tle::TLE(), // TLEs - timestamps // Timestamps + // Setup Projecition, based off N19 + geodetic::projection::LEOScanProjectorSettings proj_settings = { + 110.6, // Scan angle + -0.01, // Roll offset + 0, // Pitch offset + -3.65, // Yaw offset + 1, // Time offset + image4.width(), // Image width + true, // Invert scan + tle::TLE(), // TLEs + timestamps // Timestamps }; + { + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/AVHRR.georef"); + } + // Identify satellite, and apply per-sat settings... int scid = most_common(spacecraft_ids.begin(), spacecraft_ids.end()); if (scid == 7) // N15 @@ -186,15 +187,18 @@ namespace noaa } else if (scid == 13) // N18 { - norad = 25338; + norad = 28654; logger->info("Identified NOAA-18!"); + proj_settings.scan_angle = 110.4; + proj_settings.roll_offset = -0.1; + proj_settings.yaw_offset = -3.2; } else if (scid == 15) // N19 { norad = 33591; logger->info("Identified NOAA-19!"); - proj_settings.az_offset = -1.5; - proj_settings.proj_scale = 2.339; + //proj_settings.az_offset = -1.5; + //proj_settings.proj_scale = 2.339; } else { @@ -204,26 +208,27 @@ namespace noaa // Load TLEs now proj_settings.sat_tle = tle::getTLEfromNORAD(norad); - projection::LEOScanProjector projector(proj_settings); + geodetic::projection::LEOScanProjector projector(proj_settings); { - projection::proj_file::LEO_GeodeticReferenceFile geofile = projection::proj_file::leoRefFileFromProjector(norad, proj_settings); - projection::proj_file::writeReferenceFile(geofile, directory + "/AVHRR.georef"); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile geofile = geodetic::projection::proj_file::leoRefFileFromProjector(norad, proj_settings); + geodetic::projection::proj_file::writeReferenceFile(geofile, directory + "/AVHRR.georef"); } logger->info("Projected channel 4..."); - cimg_library::CImg projected_image = projection::projectLEOToEquirectangularMapped(image4, projector, 2048 * 4, 1024 * 4, 1); + cimg_library::CImg projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image4, projector, 2048 * 4, 1024 * 4, 1); WRITE_IMAGE(projected_image, directory + "/AVHRR-4-PROJ.png"); - cimg_library::CImg image321(2048, reader.lines, 1, 3); + cimg_library::CImg image221(2048, reader.lines, 1, 3); { - image321.draw_image(0, 0, 0, 0, image2); - image321.draw_image(0, 0, 0, 1, image2); - image321.draw_image(0, 0, 0, 2, image1); + image221.draw_image(0, 0, 0, 0, image2); + image221.draw_image(0, 0, 0, 1, image2); + image221.draw_image(0, 0, 0, 2, image1); } + image221.equalize(1000); logger->info("Projected channel 221..."); - projected_image = projection::projectLEOToEquirectangularMapped(image321, projector, 2048 * 4, 1024 * 4, 3); + projected_image = geodetic::projection::projectLEOToEquirectangularMapped(image221, projector, 2048 * 4, 1024 * 4, 3); WRITE_IMAGE(projected_image, directory + "/AVHRR-RGB-221-PROJ.png"); } } diff --git a/src-interface/projection/projection.cpp b/src-interface/projection/projection.cpp index 230eb12e7..e2d8570df 100644 --- a/src-interface/projection/projection.cpp +++ b/src-interface/projection/projection.cpp @@ -20,14 +20,14 @@ std::string getDirPath(); #include "common/map/map_drawer.h" #include "resources.h" -#include "common/projection/stereo.h" -#include "common/projection/geos.h" +#include "common/geodetic/projection/stereo.h" +#include "common/geodetic/projection/geos.h" #include "common/map/maidenhead.h" -#include "common/projection/proj_file.h" -#include "common/projection/satellite_reprojector.h" +#include "common/geodetic/projection/proj_file.h" +#include "common/geodetic/projection/satellite_reprojector.h" #include #include "global.h" -#include "common/projection/geo_projection.h" +#include "common/geodetic/projection/geo_projection.h" namespace projection { @@ -55,20 +55,25 @@ namespace projection std::string filename; std::string timestamp; cimg_library::CImg image; - std::shared_ptr georef; + std::shared_ptr georef; bool show; + float opacity; + + float progress = 0; }; std::vector filesToProject; // Utils - projection::StereoProjection proj_stereo; - projection::GEOSProjection proj_geos; + geodetic::projection::StereoProjection proj_stereo; + geodetic::projection::GEOSProjection proj_geos; std::function(float, float, int, int)> projectionFunc; bool isFirstUiRun = false; bool rendering = false; bool isRenderDone = false; + float render_progress = 0; + void initProjection() { projected_image = cimg_library::CImg(output_width, output_height, 1, 3, 0); @@ -146,22 +151,22 @@ namespace projection if (toProj.georef->file_type == 0) { logger->info("Reprojecting Equiectangular..."); - projection::projectEQUIToproj(toProj.image, projected_image, toProj.image.spectrum(), projectionFunc); + geodetic::projection::projectEQUIToproj(toProj.image, projected_image, toProj.image.spectrum(), projectionFunc, toProj.opacity, (float *)&toProj.progress); } - else if (toProj.georef->file_type == projection::proj_file::LEO_TYPE) + else if (toProj.georef->file_type == geodetic::projection::proj_file::LEO_TYPE) { - projection::proj_file::LEO_GeodeticReferenceFile leofile = *((projection::proj_file::LEO_GeodeticReferenceFile *)toProj.georef.get()); - projection::LEOScanProjectorSettings settings = leoProjectionRefFile(leofile); - projection::LEOScanProjector projector(settings); + geodetic::projection::proj_file::LEO_GeodeticReferenceFile leofile = *((geodetic::projection::proj_file::LEO_GeodeticReferenceFile *)toProj.georef.get()); + geodetic::projection::LEOScanProjectorSettings settings = leoProjectionRefFile(leofile); + geodetic::projection::LEOScanProjector projector(settings); logger->info("Reprojecting LEO..."); - projection::reprojectLEOtoProj(toProj.image, projector, projected_image, toProj.image.spectrum(), projectionFunc); + geodetic::projection::reprojectLEOtoProj(toProj.image, projector, projected_image, toProj.image.spectrum(), projectionFunc, toProj.opacity, (float *)&toProj.progress); } - else if (toProj.georef->file_type == projection::proj_file::GEO_TYPE) + else if (toProj.georef->file_type == geodetic::projection::proj_file::GEO_TYPE) { - projection::proj_file::GEO_GeodeticReferenceFile gsofile = *((projection::proj_file::GEO_GeodeticReferenceFile *)toProj.georef.get()); + geodetic::projection::proj_file::GEO_GeodeticReferenceFile gsofile = *((geodetic::projection::proj_file::GEO_GeodeticReferenceFile *)toProj.georef.get()); logger->info("Reprojecting GEO..."); - projection::GEOProjector projector = projection::proj_file::geoProjectionRefFile(gsofile); - projection::reprojectGEOtoProj(toProj.image, projector, projected_image, toProj.image.spectrum(), projectionFunc); + geodetic::projection::GEOProjector projector = geodetic::projection::proj_file::geoProjectionRefFile(gsofile); + geodetic::projection::reprojectGEOtoProj(toProj.image, projector, projected_image, toProj.image.spectrum(), projectionFunc, toProj.opacity, (float *)&toProj.progress); } } @@ -193,7 +198,9 @@ namespace projection updateImageTexture(textureID, textureBuffer, output_width, output_height); if (isRenderDone) + { isRenderDone = false; + } } ImGui::SetNextWindowPos({0, 0}); @@ -266,6 +273,22 @@ namespace projection { destroyProjection(); } + + // Calculate progress, if rendering + if (rendering) + { + render_progress = 0; + for (const FileToProject &toProj : filesToProject) + render_progress += toProj.progress; + render_progress /= filesToProject.size(); + } + else + { + render_progress = 0; + } + + ImGui::SameLine(); + ImGui::ProgressBar(render_progress); } ImGui::End(); @@ -293,9 +316,7 @@ namespace projection } ImGui::SameLine(); if (ImGui::Button("View")) - { toProj.show = true; - } } ImGui::EndTable(); @@ -355,21 +376,21 @@ namespace projection cimg_library::CImg new_image(newfile_image); new_image.normalize(0, 65535); - std::shared_ptr new_geofile; + std::shared_ptr new_geofile; if (use_equirectangular) { - new_geofile = std::make_shared(); + new_geofile = std::make_shared(); new_geofile->file_type = 0; new_geofile->utc_timestamp_seconds = 0; } else { - new_geofile = projection::proj_file::readReferenceFile(std::string(newfile_georef)); + new_geofile = geodetic::projection::proj_file::readReferenceFile(std::string(newfile_georef)); } - if (new_geofile->file_type == projection::proj_file::GEO_TYPE) + if (new_geofile->file_type == geodetic::projection::proj_file::GEO_TYPE) { - projection::proj_file::GEO_GeodeticReferenceFile gsofile = *((projection::proj_file::GEO_GeodeticReferenceFile *)new_geofile.get()); + geodetic::projection::proj_file::GEO_GeodeticReferenceFile gsofile = *((geodetic::projection::proj_file::GEO_GeodeticReferenceFile *)new_geofile.get()); new_image.resize(gsofile.image_width, gsofile.image_height); // Safety } @@ -386,7 +407,8 @@ namespace projection timestamp + " UTC", new_image, new_geofile, - false}); + false, + 1.0}); std::fill(newfile_image, &newfile_image[100], 0); std::fill(newfile_georef, &newfile_georef[100], 0); @@ -408,199 +430,198 @@ namespace projection { ImGui::Begin(std::string("View georef for " + std::string(toProj.filename)).c_str(), &toProj.show); { - if (ImGui::BeginTable("FileContentsLEO", 2, ImGuiTableFlags_Borders)) + if (ImGui::BeginTabBar(std::string("##tabbar" + toProj.timestamp).c_str())) { - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("File Type"); - ImGui::TableSetColumnIndex(1); - if (toProj.georef->file_type == 0) - ImGui::Text("0, Empty"); - else if (toProj.georef->file_type == 1) - ImGui::Text("1, GEO"); - else if (toProj.georef->file_type == 1) - ImGui::Text("2, LEO"); - else - ImGui::Text("Invalid"); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Timestamp"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%s", toProj.timestamp.c_str()); - - if (toProj.georef->file_type == projection::proj_file::GEO_TYPE) + if (ImGui::BeginTabItem("Contens")) { - projection::proj_file::GEO_GeodeticReferenceFile gsofile = *((projection::proj_file::GEO_GeodeticReferenceFile *)toProj.georef.get()); + if (ImGui::BeginTable("FileContentsLEO", 2, ImGuiTableFlags_Borders)) + { + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("File Type"); + ImGui::TableSetColumnIndex(1); + if (toProj.georef->file_type == 0) + ImGui::Text("0, Empty"); + else if (toProj.georef->file_type == 1) + ImGui::Text("1, GEO"); + else if (toProj.georef->file_type == 2) + ImGui::Text("2, LEO"); + else + ImGui::Text("Invalid"); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("NORAD"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%d", gsofile.norad); + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Timestamp"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%s", toProj.timestamp.c_str()); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Satellite Longitude"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", gsofile.position_longitude); + if (toProj.georef->file_type == geodetic::projection::proj_file::GEO_TYPE) + { + geodetic::projection::proj_file::GEO_GeodeticReferenceFile gsofile = *((geodetic::projection::proj_file::GEO_GeodeticReferenceFile *)toProj.georef.get()); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Satellite Altitude"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", gsofile.position_height); + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("NORAD"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%d", gsofile.norad); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Projection Type"); - ImGui::TableSetColumnIndex(1); - if (gsofile.projection_type == 0) - ImGui::Text("0, Full Disk"); - else - ImGui::Text("Invalid"); + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Satellite Longitude"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%f", gsofile.position_longitude); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Image Width"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%d", gsofile.image_width); + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Satellite Altitude"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%f", gsofile.position_height); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Image Height"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%d", gsofile.image_height); + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Projection Type"); + ImGui::TableSetColumnIndex(1); + if (gsofile.projection_type == 0) + ImGui::Text("0, Full Disk"); + else + ImGui::Text("Invalid"); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Horizontal Scale"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", gsofile.horizontal_scale); + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Image Width"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%d", gsofile.image_width); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Vertical Scale"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", gsofile.vertical_scale); + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Image Height"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%d", gsofile.image_height); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Horizontal Offset"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", gsofile.horizontal_offset); + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Horizontal Scale"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%f", gsofile.horizontal_scale); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Vertical Offset"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", gsofile.vertical_offset); + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Vertical Scale"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%f", gsofile.vertical_scale); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Sweep X"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%d", gsofile.proj_sweep_x); + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Horizontal Offset"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%f", gsofile.horizontal_offset); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Vertical Offset"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%f", gsofile.vertical_offset); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Sweep X"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%d", gsofile.proj_sweep_x); + } + else if (toProj.georef->file_type == geodetic::projection::proj_file::LEO_TYPE) + { + geodetic::projection::proj_file::LEO_GeodeticReferenceFile leofile = *((geodetic::projection::proj_file::LEO_GeodeticReferenceFile *)toProj.georef.get()); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("NORAD"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%d", leofile.norad); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("TLE 1"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%s", leofile.tle_line1_data.c_str()); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("TLE 2"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%s", leofile.tle_line2_data.c_str()); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Projection Type"); + ImGui::TableSetColumnIndex(1); + if (leofile.projection_type == 0) + ImGui::Text("0, Single Scanline"); + else + ImGui::Text("Invalid"); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Scan Angle"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%f", leofile.scan_angle); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Roll Offset"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%f", leofile.roll_offset); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Pitch Offset"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%f", leofile.pitch_offset); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Yaw Offset"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%f", leofile.yaw_offset); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Time offset"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%f", leofile.time_offset); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Image width"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%d", leofile.image_width); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Invert scan"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%d", leofile.invert_scan); + + ImGui::TableNextRow(); + ImGui::TableSetColumnIndex(0); + ImGui::Text("Timestamp count"); + ImGui::TableSetColumnIndex(1); + ImGui::Text("%d", leofile.timestamp_count); + } + + ImGui::EndTable(); + } + + ImGui::EndTabItem(); } - else if (toProj.georef->file_type == projection::proj_file::LEO_TYPE) + if (ImGui::BeginTabItem("Settings")) { - projection::proj_file::LEO_GeodeticReferenceFile leofile = *((projection::proj_file::LEO_GeodeticReferenceFile *)toProj.georef.get()); + ImGui::InputFloat(std::string("##opacity" + toProj.timestamp).c_str(), &toProj.opacity); + ImGui::SameLine(); + ImGui::Text("Opacity"); - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("NORAD"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%d", leofile.norad); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("TLE 1"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%s", leofile.tle_line1_data.c_str()); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("TLE 2"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%s", leofile.tle_line2_data.c_str()); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Projection Type"); - ImGui::TableSetColumnIndex(1); - if (leofile.projection_type == 0) - ImGui::Text("0, Single Scanline"); - else - ImGui::Text("Invalid"); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Correction Swath"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%d", leofile.correction_swath); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Correction Res"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", leofile.correction_res); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Correction Height"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", leofile.correction_height); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Instrument Swath"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", leofile.instrument_swath); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Projection Scale"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", leofile.proj_scale); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Azimuth Offset"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", leofile.az_offset); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Tilt offset"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", leofile.tilt_offset); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Time offset"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%f", leofile.time_offset); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Image width"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%d", leofile.image_width); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Invert scan"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%d", leofile.invert_scan); - - ImGui::TableNextRow(); - ImGui::TableSetColumnIndex(0); - ImGui::Text("Timestamp count"); - ImGui::TableSetColumnIndex(1); - ImGui::Text("%d", leofile.timestamp_count); + ImGui::EndTabItem(); } - - ImGui::EndTable(); + ImGui::EndTabBar(); } } ImGui::End();