diff --git a/CMakeLists.txt b/CMakeLists.txt index dc3a49b..d9730be 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -2,7 +2,15 @@ cmake_minimum_required(VERSION 3.5) include(ExternalProject) -project(Meteordemod LANGUAGES CXX) +project(meteordemod + VERSION 2.6.5 + LANGUAGES CXX +) + +configure_file(cmake/version.h.in version.h) + +set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -DVERSION_BUILD=${VERSION_BUILD_NUMBER}") +set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DVERSION_BUILD=${VERSION_BUILD_NUMBER}") set(CMAKE_CXX_STANDARD 17) set(CMAKE_CXX_STANDARD_REQUIRED ON) @@ -13,6 +21,10 @@ if (NOT CMAKE_BUILD_TYPE OR CMAKE_BUILD_TYPE STREQUAL "") endif() message(STATUS "Build type: ${CMAKE_BUILD_TYPE}") +if(CMAKE_BUILD_TYPE STREQUAL "Debug") + SET(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fsanitize=address") +endif() + if("${CMAKE_CXX_COMPILER_ID}" STREQUAL "GNU") if(CMAKE_CXX_COMPILER_VERSION VERSION_LESS "4.9.2") @@ -40,6 +52,7 @@ ExternalProject_Add(libcorrect ) find_package(OpenCV) +find_package(OpenCL) link_directories( ${CMAKE_BINARY_DIR}/sgp4-build/libsgp4 @@ -50,48 +63,33 @@ add_definitions(-D_USE_MATH_DEFINES -D_SILENCE_EXPERIMENTAL_FILESYSTEM_DEPRECATI add_executable(meteordemod main.cpp - imageproc/spreadimage.cpp - imageproc/spreadimage.h imageproc/threatimage.cpp - imageproc/threatimage.h - decoder/bitio.cpp - decoder/bitio.h + imageproc/projectimage.cpp + imageproc/blendimages.cpp + imageproc/tps.cpp decoder/correlation.cpp - decoder/correlation.h - decoder/meteorimage.cpp - decoder/meteorimage.h - decoder/packetparser.cpp - decoder/packetparser.h decoder/reedsolomon.cpp - decoder/reedsolomon.h decoder/viterbi.cpp - decoder/viterbi.h decoder/deinterleaver.cpp decoder/meteordecoder.cpp - decoder/meteordecoder.h + decoder/protocol/ccsds.cpp + decoder/protocol/vcdu.cpp + decoder/protocol/lrpt/decoder.cpp + decoder/protocol/lrpt/msumr/segment.cpp + decoder/protocol/lrpt/msumr/bitio.cpp + decoder/protocol/lrpt/msumr/image.cpp common/settings.cpp - common/settings.h - common/version.h tools/matrix.cpp - tools/matrix.h tools/tlereader.cpp - tools/tlereader.h tools/vector.cpp - tools/vector.h tools/pixelgeolocationcalculator.cpp - tools/pixelgeolocationcalculator.h tools/databuffer.cpp - tools/databuffer.h tools/iniparser.cpp - tools/iniparser.h tools/threadpool.cpp - tools/threadpool.h + tools/opencl.cpp GIS/shapereader.cpp - GIS/shapereader.h GIS/shaperenderer.cpp - GIS/shaperenderer.h GIS/dbfilereader.cpp - GIS/dbfilereader.h DSP/meteordemodulator.cpp DSP/agc.cpp DSP/pll.cpp @@ -118,9 +116,19 @@ include_directories( ${CMAKE_SOURCE_DIR}/external/libcorrect/include ) +target_include_directories(meteordemod PUBLIC + "${PROJECT_BINARY_DIR}" +) + add_dependencies(meteordemod sgp4) add_dependencies(meteordemod libcorrect) +if(OpenCL_FOUND) + target_compile_definitions(meteordemod PUBLIC OPENCL_FOUND=1) + include_directories(${OpenCL_INCLUDE_DIRS}) + target_link_libraries(meteordemod ${OpenCL_LIBRARIES}) +endif() + if(WIN32) target_link_libraries(meteordemod ${OpenCV_LIBS} @@ -150,9 +158,6 @@ if(UNIX AND NOT APPLE) find_file (DEBIAN_FOUND debian_version debconf.conf PATHS /etc) if (DEBIAN_FOUND) SET(CPACK_GENERATOR "DEB") - SET(CPACK_PACKAGE_VERSION_MAJOR "2") - SET(CPACK_PACKAGE_VERSION_MINOR "5") - SET(CPACK_PACKAGE_VERSION_PATCH "5") SET(CPACK_DEBIAN_PACKAGE_MAINTAINER "Digitelektro") SET(CPACK_DEBIAN_PACKAGE_HOMEPAGE "https://github.com/Digitelektro/MeteorDemod") SET(CPACK_DEBIAN_PACKAGE_DESCRIPTION "Russian Meteor M2 weather satellite data decoder") diff --git a/GIS/shapereader.h b/GIS/shapereader.h index 41fe13f..4b8e120 100644 --- a/GIS/shapereader.h +++ b/GIS/shapereader.h @@ -3,6 +3,7 @@ #include #include +#include #include #include #include @@ -170,6 +171,10 @@ class ShapeReader { pointBuffer.valueAtIndex(index, y, LittleEndian); } + bool operator==(const Point& rhs) { + return (std::abs(x - rhs.x) <= std::numeric_limits::epsilon()) && (std::abs(y - rhs.y) <= std::numeric_limits::epsilon()); + } + double x; double y; }; @@ -235,6 +240,19 @@ class ShapeReader { return *this; } } + PolyLineIterator operator++(int) { + auto temp = *this; + if(mNumberOfPoint < mPolyLineHeader.numberOfpoints) { + DataBuffer pointBuffer(16); + mInputStream.read(reinterpret_cast(pointBuffer.buffer()), pointBuffer.size()); + point = Point(pointBuffer); + mNumberOfPoint++; + } else { + mNumberOfPoint = 0; + point = Point(); + } + return temp; + } PolyLineIterator begin() { int filePosition = mRecordPosition + 12; // Recordpos + Recordheader diff --git a/GIS/shaperenderer.cpp b/GIS/shaperenderer.cpp index c6eb7c2..99125ff 100644 --- a/GIS/shaperenderer.cpp +++ b/GIS/shaperenderer.cpp @@ -36,6 +36,11 @@ void GIS::ShapeRenderer::drawShape(const cv::Mat& src, Transform_t transform) { for(polyLineIterator->begin(); *polyLineIterator != polyLineIterator->end(); ++(*polyLineIterator)) { if(transform(polyLineIterator->point.y, polyLineIterator->point.x)) { polyLines.push_back(cv::Point2d(polyLineIterator->point.y, polyLineIterator->point.x)); + } else { + if(polyLines.size() > 1) { + cv::polylines(src, polyLines, false, mColor, mThicknes); + polyLines.clear(); + } } } @@ -45,6 +50,42 @@ void GIS::ShapeRenderer::drawShape(const cv::Mat& src, Transform_t transform) { } } } + } else if(getShapeType() == ShapeReader::ShapeType::stPolygon) { + auto recordIterator = getRecordIterator(); + + if(recordIterator) { + for(recordIterator->begin(); *recordIterator != recordIterator->end(); ++(*recordIterator)) { + auto polyLineIterator = getPolyLineIterator(*recordIterator); + + if(polyLineIterator) { + bool isFirst = true; + ShapeReader::Point first; + std::vector polyLines; + for(polyLineIterator->begin(); *polyLineIterator != polyLineIterator->end(); ++(*polyLineIterator)) { + if(!isFirst && (first == polyLineIterator->point)) { + if(polyLines.size() > 1) { + cv::polylines(src, polyLines, false, mColor, mThicknes); + } + isFirst = true; + polyLines.clear(); + continue; + } + if(isFirst) { + first = polyLineIterator->point; + isFirst = false; + } + if(transform(polyLineIterator->point.y, polyLineIterator->point.x)) { + polyLines.push_back(cv::Point2d(polyLineIterator->point.y, polyLineIterator->point.x)); + } else { + if(polyLines.size() > 1) { + cv::polylines(src, polyLines, false, mColor, mThicknes); + polyLines.clear(); + } + } + } + } + } + } } else if(getShapeType() == ShapeReader::ShapeType::stPoint) { auto recordIterator = getRecordIterator(); diff --git a/cmake/version.h.in b/cmake/version.h.in new file mode 100644 index 0000000..b849750 --- /dev/null +++ b/cmake/version.h.in @@ -0,0 +1,3 @@ +#define VERSION_MAJOR @meteordemod_VERSION_MAJOR@ +#define VERSION_MINOR @meteordemod_VERSION_MINOR@ +#define VERSION_FIX @meteordemod_VERSION_PATCH@ diff --git a/common/settings.cpp b/common/settings.cpp index 67e341a..68b9f8d 100644 --- a/common/settings.cpp +++ b/common/settings.cpp @@ -79,9 +79,10 @@ void Settings::parseIni(const std::string& path) { ini::extract(mIniParser.sections["Program"]["CompositeProjectionScale"], mCompositeProjectionScale, 0.75f); ini::extract(mIniParser.sections["Program"]["CompositeAzimuthalEquidistantProjection"], mCompositeEquadistantProjection, true); ini::extract(mIniParser.sections["Program"]["CompositeMercatorProjection"], mCompositeMercatorProjection, false); - ini::extract(mIniParser.sections["Program"]["GenerateComposite123"], mGenerateComposite123, true); + ini::extract(mIniParser.sections["Program"]["GenerateComposite321"], mGenerateComposite321, true); ini::extract(mIniParser.sections["Program"]["GenerateComposite125"], mGenerateComposite125, true); ini::extract(mIniParser.sections["Program"]["GenerateComposite221"], mGenerateComposite221, true); + ini::extract(mIniParser.sections["Program"]["GenerateComposite224"], mGenerateComposite224, true); ini::extract(mIniParser.sections["Program"]["GenerateComposite68"], mGenerateComposite68, true); ini::extract(mIniParser.sections["Program"]["GenerateCompositeThermal"], mGenerateCompositeThermal, true); ini::extract(mIniParser.sections["Program"]["GenerateComposite68Rain"], mGenerateComposite68Rain, true); diff --git a/common/settings.h b/common/settings.h index 75a4b49..9f986c6 100644 --- a/common/settings.h +++ b/common/settings.h @@ -58,6 +58,15 @@ class Settings { uint8_t B; }; + struct ProjectionSetting { + std::string satelliteNameInTLE; + float scanAngle; + float yaw; + float pitch; + float roll; + float timeOffsetMs; + }; + friend std::istream& operator>>(std::istream& is, HTMLColor& color) { std::string rgb; is >> rgb; @@ -105,35 +114,15 @@ class Settings { return mJpegQuality; } - std::string getSatNameInTLE() { - std::string name; - ini::extract(mIniParser.sections[getSateliteName()]["SatNameInTLE"], name); - return name; - } - float getScanAngle() { - float angle; - ini::extract(mIniParser.sections[getSateliteName()]["ScanAngle"], angle, 110.8f); - return angle; - } - float getRoll() { - float roll; - ini::extract(mIniParser.sections[getSateliteName()]["Roll"], roll, 0.0f); - return roll; - } - float getPitch() { - float pitch; - ini::extract(mIniParser.sections[getSateliteName()]["Pitch"], pitch, 0.0f); - return pitch; - } - float getYaw() { - float yaw; - ini::extract(mIniParser.sections[getSateliteName()]["Yaw"], yaw, 0.0f); - return yaw; - } - int getTimeOffsetSec() { - int timeOffset; - ini::extract(mIniParser.sections["Program"]["TimeOffset"], timeOffset, 0); - return timeOffset; + ProjectionSetting getProjectionSetting(const std::string& satellite) { + ProjectionSetting result; + ini::extract(mIniParser.sections[satellite]["SatNameInTLE"], result.satelliteNameInTLE); + ini::extract(mIniParser.sections[satellite]["ScanAngle"], result.scanAngle, 110.8f); + ini::extract(mIniParser.sections[satellite]["Roll"], result.roll, 0.0f); + ini::extract(mIniParser.sections[satellite]["Pitch"], result.pitch, 0.0f); + ini::extract(mIniParser.sections[satellite]["Yaw"], result.yaw, 0.0f); + ini::extract(mIniParser.sections[satellite]["TimeOffset"], result.timeOffsetMs, 0.0f); + return result; } bool equadistantProjection() const { @@ -164,8 +153,8 @@ class Settings { bool compositeMercatorProjection() const { return mCompositeMercatorProjection; } - bool generateComposite123() const { - return mGenerateComposite123; + bool generateComposite321() const { + return mGenerateComposite321; } bool generateComposite125() const { return mGenerateComposite125; @@ -173,6 +162,9 @@ class Settings { bool generateComposite221() const { return mGenerateComposite221; } + bool generateComposite224() const { + return mGenerateComposite224; + } bool generateComposite68() const { return mGenerateComposite68; } @@ -310,9 +302,10 @@ class Settings { bool mCompositeEquadistantProjection; bool mCompositeMercatorProjection; - bool mGenerateComposite123; + bool mGenerateComposite321; bool mGenerateComposite125; bool mGenerateComposite221; + bool mGenerateComposite224; bool mGenerateComposite68; bool mGenerateCompositeThermal; bool mGenerateComposite68Rain; diff --git a/common/version.h b/common/version.h deleted file mode 100644 index 273615c..0000000 --- a/common/version.h +++ /dev/null @@ -1,8 +0,0 @@ -#ifndef VERSION_H -#define VERSION_H - -#define VERSION_MAJOR 2 -#define VERSION_MINOR 5 -#define VERSION_FIX 5 - -#endif // VERSION_H diff --git a/decoder/bitio.h b/decoder/bitio.h deleted file mode 100644 index 969045c..0000000 --- a/decoder/bitio.h +++ /dev/null @@ -1,22 +0,0 @@ -#ifndef BITIO_H -#define BITIO_H - -// Based on: https://github.com/artlav/meteor_decoder/blob/master/alib/bitop.pas - -#include - -class BitIOConst { - public: - BitIOConst(const uint8_t* bytes); - - uint32_t peekBits(int n); - void advanceBits(int n); - uint32_t fetchBits(int n); - - private: - const uint8_t* mpBytes; - int mPos; -}; - - -#endif // BITIO_H diff --git a/decoder/meteordecoder.cpp b/decoder/meteordecoder.cpp index e5eb77f..8374b74 100644 --- a/decoder/meteordecoder.cpp +++ b/decoder/meteordecoder.cpp @@ -1,12 +1,13 @@ #include "meteordecoder.h" +#include MeteorDecoder::MeteorDecoder(bool deInterleave, bool oqpsk, bool differentialDecode) : mDeInterleave(deInterleave) , mDifferentialDecode(differentialDecode) , mCorrelation(differentialDecode ? sSynchWordOQPSK : sSynchWordQPSK, oqpsk) {} -size_t MeteorDecoder::decode(uint8_t* softBits, size_t length) { +size_t MeteorDecoder::decode(uint8_t* softBits, size_t length, std::function callback) { size_t decodedPacketCounter = 0; size_t syncWordFound = 0; @@ -17,7 +18,7 @@ size_t MeteorDecoder::decode(uint8_t* softBits, size_t length) { length = outLen; } - mCorrelation.correlate(softBits, length, [&softBits, length, &decodedPacketCounter, &syncWordFound, this](Correlation::CorellationResult correlationResult, Correlation::PhaseShift phaseShift) { + mCorrelation.correlate(softBits, length, [&softBits, length, &decodedPacketCounter, &syncWordFound, &callback, this](Correlation::CorellationResult correlationResult, Correlation::PhaseShift phaseShift) { bool packetOk; uint32_t processedBits = 0; @@ -28,7 +29,7 @@ size_t MeteorDecoder::decode(uint8_t* softBits, size_t length) { return processedBits; } - memcpy(mDataTodecode, &softBits[correlationResult.pos + processedBits], 16384); + std::copy(&softBits[correlationResult.pos + processedBits], &softBits[correlationResult.pos + processedBits] + 16384, mDataTodecode); mCorrelation.rotateSoftIqInPlace(mDataTodecode, 16384, phaseShift); @@ -53,17 +54,18 @@ size_t MeteorDecoder::decode(uint8_t* softBits, size_t length) { for(int i = 0; i < 4; i++) { mReedSolomon.deinterleave(mViterbiResult + 4, i, 4); rsResult[i] = mReedSolomon.decode(); - mReedSolomon.interleave(mDecodedPacket, i, 4); + mReedSolomon.interleave(mDecodedPacket + 4, i, 4); } std::cout << "SyncWordFound:" << syncWordFound << " | Decoded Packets:" << decodedPacketCounter << " | Current Pos:" << (correlationResult.pos + processedBits) << " | Phase:" << phaseShift << " | synch:" << std::hex - << last_sync_ << " | RS: (" << std::dec << rsResult[0] << ", " << rsResult[1] << ", " << rsResult[2] << ", " << rsResult[3] << ")" + << last_sync_ << " | BER: " << mViterbi.getLastBER() << " | RS: (" << std::dec << rsResult[0] << ", " << rsResult[1] << ", " << rsResult[2] << ", " << rsResult[3] << ")" << "\t\t\r"; packetOk = (rsResult[0] != -1) && (rsResult[1] != -1) && (rsResult[2] != -1) && (rsResult[3] != -1); if(packetOk) { - parseFrame(mDecodedPacket, 892); + std::copy(mViterbiResult, mViterbiResult + 4, mDecodedPacket); + callback(mDecodedPacket, sizeof(mDecodedPacket)); decodedPacketCounter++; processedBits += 16384; } diff --git a/decoder/meteordecoder.h b/decoder/meteordecoder.h index 9f4e5ec..fb43264 100644 --- a/decoder/meteordecoder.h +++ b/decoder/meteordecoder.h @@ -4,14 +4,14 @@ #include #include +#include #include "correlation.h" #include "deinterleaver.h" -#include "packetparser.h" #include "reedsolomon.h" #include "viterbi.h" -class MeteorDecoder : public PacketParser { +class MeteorDecoder { private: static constexpr uint8_t PRAND[] = {0xff, 0x48, 0x0e, 0xc0, 0x9a, 0x0d, 0x70, 0xbc, 0x8e, 0x2c, 0x93, 0xad, 0xa7, 0xb7, 0x46, 0xce, 0x5a, 0x97, 0x7d, 0xcc, 0x32, 0xa2, 0xbf, 0x3e, 0x0a, 0x10, 0xf1, 0x88, 0x94, 0xcd, 0xea, 0xb1, 0xfe, 0x90, 0x1d, 0x81, 0x34, @@ -27,7 +27,7 @@ class MeteorDecoder : public PacketParser { MeteorDecoder() = delete; MeteorDecoder(bool deInterleave, bool oqpsk, bool differentialDecode); - size_t decode(uint8_t* softBits, size_t length); + size_t decode(uint8_t* softBits, size_t length, std::function callback); private: bool mDeInterleave; diff --git a/decoder/packetparser.cpp b/decoder/packetparser.cpp deleted file mode 100644 index 5f2dd8e..0000000 --- a/decoder/packetparser.cpp +++ /dev/null @@ -1,142 +0,0 @@ -#include "packetparser.h" -const int PacketParser::PACKET_FULL_MARK = 2047; - -PacketParser::PacketParser() - : MeteorImage() - , mLastFrame(0) - , mPacketOff(0) - , mPartialPacket(false) - , mFirstTimeStamp(0) - , mLastTimeStamp(0) - , mFirstTime(true) {} - -void PacketParser::parseFrame(const uint8_t* frame, int len) { - int n; - - uint16_t w = (frame[0] << 8) | frame[1]; - int ver = w >> 14; - int ssid = (w >> 6) & 0xff; - int fid = w & 0x3f; - - int frameCount = (frame[2] << 16) | (frame[3] << 8) | frame[4]; - - w = (frame[8] << 8) | frame[9]; - uint8_t hdr_mark = w >> 11; - uint16_t hdr_off = w & 0x7ff; - - if(ver == 0 || fid == 0) - return; // Empty packet - - int data_len = len - 10; - if(frameCount == mLastFrame + 1) { - if(mPartialPacket) { - if(hdr_off == PACKET_FULL_MARK) { - hdr_off = len - 10; - std::move(frame + 10, frame + 10 + hdr_off, mPacketBuffer.begin() + mPacketOff); - mPacketOff += hdr_off; - } else { - std::move(frame + 10, frame + 10 + hdr_off, mPacketBuffer.begin() + mPacketOff); - n = parsePartial(mPacketBuffer.data(), mPacketOff + hdr_off); - } - } - } else { - if(hdr_off == PACKET_FULL_MARK) { - return; - } - mPartialPacket = false; - mPacketOff = 0; - } - mLastFrame = frameCount; - - data_len -= hdr_off; - int off = hdr_off; - while(data_len > 0) { - n = parsePartial(frame + 10 + off, data_len); - if(mPartialPacket) { - mPacketOff = data_len; - std::move(frame + 10 + off, frame + 10 + off + mPacketOff, mPacketBuffer.begin()); - break; - } else { - off += n; - data_len -= n; - } - } -} - -int PacketParser::parsePartial(const uint8_t* packet, int len) { - if(len < 6) { - mPartialPacket = true; - return 0; - } - - int len_pck = (packet[4] << 8) | packet[5]; - if(len_pck >= len - 6) { - mPartialPacket = true; - return 0; - } - - parseAPD(packet, len_pck + 1); - - mPartialPacket = false; - return len_pck + 6 + 1; -} - -void PacketParser::parseAPD(const uint8_t* packet, int len) { - uint16_t w = (packet[0] << 8) | packet[1]; - int sec = (w >> 11) & 1; - int apd = w & 0x7ff; - - int pck_cnt = ((packet[2] << 8) | packet[3]) & 0x3fff; - int len_pck = (packet[4] << 8) | packet[5]; - - int ms = (packet[8] << 24) | (packet[9] << 16) | (packet[10] << 8) | packet[11]; - - if(apd == 70) { - parse70(packet + 14, len - 14); - } else { - actAPD(packet + 14, len - 14, apd, pck_cnt); - } -} - -void PacketParser::actAPD(const uint8_t* packet, int len, int apd, int pck_cnt) { - int mcu_id = packet[0]; - int scan_hdr = (packet[1] << 8) | packet[2]; - int seg_hdr = (packet[3] << 8) | packet[4]; - int q = packet[5]; - - decMCUs(packet + 6, len - 6, apd, pck_cnt, mcu_id, q); -} - -void PacketParser::parse70(const uint8_t* packet, int len) { - static int prevY; - static TimeSpan prevtimeStamp(0); - - if(len < 11) { - return; - } - - int h = packet[8]; - int m = packet[9]; - int s = packet[10]; - int ms = packet[11] * 4; - - mLastTimeStamp = TimeSpan(0, h, m, s, ms * 1000); - mLastHeightAtTimeStamp = getLastY(); - if(mFirstTime) { - mFirstTime = false; - mFirstTimeStamp = mLastTimeStamp; - mFirstHeightAtTimeStamp = getLastY(); - } else { - int lines = mLastHeightAtTimeStamp - prevY; - TimeSpan elapsedTime = mLastTimeStamp - prevtimeStamp; - - // std::cout << "Lines between timestamp : " << lines << " Time elapsed " << elapsedTime << " PixelTime:" << TimeSpan(elapsedTime.Ticks() / lines) << std::endl; - } - - prevY = getLastY(); - prevtimeStamp = mLastTimeStamp; - - // std::cout << "LastHeight: " << mLastHeightAtTimeStamp << " " << mLastTimeStamp << std::endl; - - // std::cout << "Onboard time: " << h << ":" << m << ":" << s << "." << ms; -} diff --git a/decoder/packetparser.h b/decoder/packetparser.h deleted file mode 100644 index f51bde2..0000000 --- a/decoder/packetparser.h +++ /dev/null @@ -1,56 +0,0 @@ -#ifndef PACKETPARSER_H -#define PACKETPARSER_H - -#include - -#include - -#include "TimeSpan.h" -#include "meteorimage.h" - -// Ported from https://github.com/artlav/meteor_decoder/blob/master/met_packet.pas - -class PacketParser : public MeteorImage { - public: - PacketParser(); - - void parseFrame(const uint8_t* frame, int len); - - public: - const TimeSpan getFirstTimeStamp() const { - int64_t pixelTime = (mLastTimeStamp - mFirstTimeStamp).Ticks() / (mLastHeightAtTimeStamp - mFirstHeightAtTimeStamp); - TimeSpan missingTime(pixelTime * mFirstHeightAtTimeStamp); - - return mFirstTimeStamp.Subtract(missingTime); - } - - TimeSpan getLastTimeStamp() const { - int64_t pixelTime = (mLastTimeStamp - mFirstTimeStamp).Ticks() / (mLastHeightAtTimeStamp - mFirstHeightAtTimeStamp); - int missingPixelsTime = ((getLastY() + 8) - mLastHeightAtTimeStamp); - TimeSpan missingTime(pixelTime * missingPixelsTime); - - return mLastTimeStamp.Add(missingTime); - } - - private: - int parsePartial(const uint8_t* packet, int len); - void parseAPD(const uint8_t* packet, int len); - void actAPD(const uint8_t* packet, int len, int apd, int pck_cnt); - void parse70(const uint8_t* packet, int len); - - private: - std::array mPacketBuffer; - int mLastFrame; - int mPacketOff; - bool mPartialPacket; - TimeSpan mFirstTimeStamp; - TimeSpan mLastTimeStamp; - int mFirstHeightAtTimeStamp; - int mLastHeightAtTimeStamp; - bool mFirstTime; - - private: - static const int PACKET_FULL_MARK; -}; - -#endif // PACKETPARSER_H diff --git a/decoder/protocol/ccsds.cpp b/decoder/protocol/ccsds.cpp new file mode 100644 index 0000000..040972d --- /dev/null +++ b/decoder/protocol/ccsds.cpp @@ -0,0 +1,46 @@ +#include "ccsds.h" + +namespace decoder { +namespace protocol { + +CCSDS::CCSDS(const uint8_t* data) + : mpData(data) { + mVersion = data[0] >> 5; + mType = (data[0] >> 4) & 0x01; + mSeconderyHeaderPresent = (data[0] >> 3) & 1; + mAPID = static_cast(data[0] & 0x03) << 8 | data[1]; + mSequenceFlag = data[2] >> 6; + mPacketSequenceCount = static_cast(data[2] & 0x3F) << 8 | data[3]; + mPacketLength = static_cast(data[4]) << 8 | data[5]; +} + +uint8_t CCSDS::getVersion() const { + return mVersion; +} + +bool CCSDS::getType() const { + return mType; +} + +bool CCSDS::isSecondaryHeaderPresent() const { + return mSeconderyHeaderPresent; +} + +uint16_t CCSDS::getAPID() const { + return mAPID; +} + +uint8_t CCSDS::getSequenceFlag() const { + return mSequenceFlag; +} + +uint16_t CCSDS::getgetSequenceCount() const { + return mPacketSequenceCount; +} + +uint16_t CCSDS::getPacketLength() const { + return mPacketLength + 1; // Packet length in CCSDS is length - 1 +} + +} // namespace protocol +} // namespace decoder \ No newline at end of file diff --git a/decoder/protocol/ccsds.h b/decoder/protocol/ccsds.h new file mode 100644 index 0000000..6b1c643 --- /dev/null +++ b/decoder/protocol/ccsds.h @@ -0,0 +1,48 @@ +#pragma once + +#include +#include + +namespace decoder { +namespace protocol { + +class CCSDS { + public: + CCSDS(const uint8_t* data); + + public: + uint8_t getVersion() const; + bool getType() const; + bool isSecondaryHeaderPresent() const; + uint16_t getAPID() const; + uint8_t getSequenceFlag() const; + uint16_t getgetSequenceCount() const; + uint16_t getPacketLength() const; + + uint16_t size() const { + return 6 + mPacketLength + 1; // header + length + } + uint16_t packetLength() const { + return mPacketLength + 1; + } + + const uint8_t* data() const { + return mpData; + } + const uint8_t* packetData() const { + return mpData + 6; // header + } + + private: + const uint8_t* mpData; + uint8_t mVersion; + bool mType; + bool mSeconderyHeaderPresent; + uint16_t mAPID; + uint8_t mSequenceFlag; + uint16_t mPacketSequenceCount; + uint16_t mPacketLength; +}; + +} // namespace protocol +} // namespace decoder \ No newline at end of file diff --git a/decoder/protocol/lrpt/decoder.cpp b/decoder/protocol/lrpt/decoder.cpp new file mode 100644 index 0000000..c8f31ca --- /dev/null +++ b/decoder/protocol/lrpt/decoder.cpp @@ -0,0 +1,148 @@ +#include "decoder.h" + +#include + +#include "msumr/segment.h" +#include "protocol/vcdu.h" + +namespace decoder { +namespace protocol { +namespace lrpt { + +std::string Decoder::serialNumberToSatName(uint8_t serialNumber) { + switch(serialNumber) { + case 1: + return "METEOR M2-1"; + break; + case 2: + return "METEOR M2-2"; + break; + case 3: + return "METEOR M2-3"; + break; + case 4: + return "METEOR M2-4"; + break; + case 5: + return "METEOR M2-5"; + break; + default: + return "METEOR M2"; + } +} + +void Decoder::process(const uint8_t* cadu) { + VCDU vcdu(cadu); + + if(vcdu.getVersion() == 0) { + return; + } + + if(vcdu.getVCID() == cVCIDAVHRR) { + cadu += 4; // Sync word + cadu += VCDU::size(); + + uint16_t firstHeaderPointer = static_cast(cadu[0] & 0b111) << 8 | cadu[1]; + cadu += 2; // remove MPDU Header, now ptr points to MPDU Data + + if(firstHeaderPointer < cNoHeaderMark && firstHeaderPointer >= cPDUDataSize) { + return; // invalid packet + } + + if(vcdu.getCounter() == (mLastFrame + 1)) { + if(mPartialPacket) { + if(firstHeaderPointer == cNoHeaderMark) { + // Current frame doesn't have header + mCcsdsPayload.insert(mCcsdsPayload.end(), cadu, cadu + cPDUDataSize); + } else { + mCcsdsPayload.insert(mCcsdsPayload.end(), cadu, cadu + firstHeaderPointer); + parsePartial(mCcsdsPayload.data(), mCcsdsPayload.size()); + mCcsdsPayload.clear(); + } + } + } else { + // We lost synch, reset state + mPartialPacket = false; + mCcsdsPayload.clear(); + if(firstHeaderPointer == cNoHeaderMark) { + // No header in this packet, cannot start collecting data + return; + } + } + + mLastFrame = vcdu.getCounter(); + + cadu += firstHeaderPointer; + uint16_t dataSize = cPDUDataSize - firstHeaderPointer; + while(dataSize > 0) { + int usedBytes = parsePartial(cadu, dataSize); + if(mPartialPacket) { + // store partial data for next round + mCcsdsPayload.insert(mCcsdsPayload.end(), cadu, cadu + dataSize); + dataSize = 0; + } else { + dataSize -= usedBytes; + cadu += usedBytes; + } + } + } +} + +uint16_t Decoder::parsePartial(const uint8_t* data, uint32_t length) { + if(length < 6) { + mPartialPacket = true; + return 0; // no enough bytes to decode header yet + } + + CCSDS ccsds(data); + + if((ccsds.size()) > length) { + mPartialPacket = true; + return 0; + } else { + processPacket(ccsds); + mPartialPacket = false; + } + return ccsds.size(); +} + +void Decoder::processPacket(const CCSDS& ccsds) { + if(ccsds.getAPID() == 64 || ccsds.getAPID() == 65 || ccsds.getAPID() == 66 || ccsds.getAPID() == 67 || ccsds.getAPID() == 68 || ccsds.getAPID() == 69) { + msumr::Segment segment(ccsds.packetData(), ccsds.getPacketLength()); + if(segment.isValid()) { + decode(ccsds.getAPID(), ccsds.getgetSequenceCount(), segment); + + /*DateTime onboardTime(2000, 12, 31); + onboardTime = onboardTime.AddDays(segment.getDay()); + onboardTime = onboardTime.AddMicroseconds(segment.getMsec() * 1000 + segment.getUsec()); + std::cout << "Onboard time: " << onboardTime << std::endl;*/ + } + } else if(ccsds.getAPID() == 70) { + parse70(ccsds); + } +} + +void Decoder::parse70(const CCSDS& ccsds) { + uint16_t h = ccsds.packetData()[16]; + uint16_t m = ccsds.packetData()[17]; + uint16_t s = ccsds.packetData()[18]; + uint16_t ms = static_cast(ccsds.packetData()[19]) * 4; + + mSerialNumber = ccsds.packetData()[20] >> 4; + + // std::cout << "Serial number : " << (int)mSerialNumber << std::endl; + + mLastTimeStamp = TimeSpan(0, h, m, s, ms * 1000); + mLastHeightAtTimeStamp = getLastY(); + if(mFirstTime) { + mFirstTime = false; + mFirstTimeStamp = mLastTimeStamp; + mFirstHeightAtTimeStamp = getLastY(); + } + + // std::cout << "Onboard time: " << h << ":" << m << ":" << s << "." << ms << std::endl; +} + +} // namespace lrpt +} // namespace protocol +} // namespace decoder \ No newline at end of file diff --git a/decoder/protocol/lrpt/decoder.h b/decoder/protocol/lrpt/decoder.h new file mode 100644 index 0000000..a713013 --- /dev/null +++ b/decoder/protocol/lrpt/decoder.h @@ -0,0 +1,69 @@ +#pragma once + +#include +#include +#include + +#include "DateTime.h" +#include "TimeSpan.h" +#include "msumr/image.h" +#include "protocol/ccsds.h" + +namespace decoder { +namespace protocol { +namespace lrpt { + +// http://planet.iitp.ru/retro/index.php?lang=en&page_type=spacecraft&page=meteor_m_n2_structure_2 +class Decoder : public msumr::Image { + public: + static std::string serialNumberToSatName(uint8_t serialNumber); + + public: + void process(const uint8_t* cadu); + + public: + const TimeSpan getFirstTimeStamp() const { + int64_t pixelTime = (mLastTimeStamp - mFirstTimeStamp).Ticks() / (mLastHeightAtTimeStamp - mFirstHeightAtTimeStamp); + TimeSpan missingTime(pixelTime * mFirstHeightAtTimeStamp); + + return mFirstTimeStamp.Subtract(missingTime); + } + + TimeSpan getLastTimeStamp() const { + int64_t pixelTime = (mLastTimeStamp - mFirstTimeStamp).Ticks() / (mLastHeightAtTimeStamp - mFirstHeightAtTimeStamp); + int missingPixelsTime = ((getLastY() + 8) - mLastHeightAtTimeStamp); + TimeSpan missingTime(pixelTime * missingPixelsTime); + + return mLastTimeStamp.Add(missingTime); + } + + uint8_t getSerialNumber() const { + return mSerialNumber; + } + + private: + uint16_t parsePartial(const uint8_t* data, uint32_t length); + void processPacket(const CCSDS& ccsds); + void parse70(const CCSDS& ccsds); + + private: + bool mHeaderFound = false; + bool mPartialPacket = false; + uint32_t mLastFrame = 0; + std::vector mCcsdsPayload; + + TimeSpan mFirstTimeStamp{0}; + TimeSpan mLastTimeStamp{0}; + int mFirstHeightAtTimeStamp = 0; + int mLastHeightAtTimeStamp = 0; + bool mFirstTime = true; + uint8_t mSerialNumber = 0; + + private: + static constexpr uint8_t cVCIDAVHRR = 5; + static constexpr int cNoHeaderMark = 0x7FF; + static constexpr uint16_t cPDUDataSize = 882; +}; +} // namespace lrpt +} // namespace protocol +} // namespace decoder \ No newline at end of file diff --git a/decoder/bitio.cpp b/decoder/protocol/lrpt/msumr/bitio.cpp similarity index 51% rename from decoder/bitio.cpp rename to decoder/protocol/lrpt/msumr/bitio.cpp index dd3a790..ff2e2c5 100644 --- a/decoder/bitio.cpp +++ b/decoder/protocol/lrpt/msumr/bitio.cpp @@ -1,14 +1,25 @@ #include "bitio.h" -BitIOConst::BitIOConst(const uint8_t* bytes) +namespace decoder { +namespace protocol { +namespace lrpt { +namespace msumr { + +BitIOConst::BitIOConst(const uint8_t* bytes, uint32_t size) : mpBytes(bytes) + , mSize(size) , mPos(0) {} uint32_t BitIOConst::peekBits(int n) { uint32_t result = 0; + int bit = 0; for(int i = 0; i < n; i++) { int p = mPos + i; - int bit = (mpBytes[p >> 3] >> (7 - (p & 7))) & 1; + if((p >> 3) < mSize) { + bit = (mpBytes[p >> 3] >> (7 - (p & 7))) & 1; + } else { + bit = 0; + } result = (result << 1) | bit; } return result; @@ -23,3 +34,8 @@ uint32_t BitIOConst::fetchBits(int n) { advanceBits(n); return result; } + +} // namespace msumr +} // namespace lrpt +} // namespace protocol +} // namespace decoder \ No newline at end of file diff --git a/decoder/protocol/lrpt/msumr/bitio.h b/decoder/protocol/lrpt/msumr/bitio.h new file mode 100644 index 0000000..b181698 --- /dev/null +++ b/decoder/protocol/lrpt/msumr/bitio.h @@ -0,0 +1,30 @@ +#pragma once + +// Based on: https://github.com/artlav/meteor_decoder/blob/master/alib/bitop.pas + +#include + +namespace decoder { +namespace protocol { +namespace lrpt { +namespace msumr { + +class BitIOConst { + public: + BitIOConst(const uint8_t* bytes, uint32_t size); + + uint32_t peekBits(int n); + void advanceBits(int n); + uint32_t fetchBits(int n); + + private: + const uint8_t* mpBytes; + uint32_t mSize; + uint32_t mPos; +}; + + +} // namespace msumr +} // namespace lrpt +} // namespace protocol +} // namespace decoder diff --git a/decoder/meteorimage.cpp b/decoder/protocol/lrpt/msumr/image.cpp similarity index 66% rename from decoder/meteorimage.cpp rename to decoder/protocol/lrpt/msumr/image.cpp index 829f0f7..f046d19 100644 --- a/decoder/meteorimage.cpp +++ b/decoder/protocol/lrpt/msumr/image.cpp @@ -1,4 +1,4 @@ -#include "meteorimage.h" +#include "image.h" #include #include @@ -7,8 +7,10 @@ #include "bitio.h" -const uint MCU_PER_PACKET = 14; -const uint MCU_PER_LINE = 196; +namespace decoder { +namespace protocol { +namespace lrpt { +namespace msumr { const std::array STANDARD_QUANTIZATION_TABLE{16, 11, 10, 16, 24, 40, 51, 61, 12, 12, 14, 19, 26, 58, 60, 55, 14, 13, 16, 24, 40, 57, 69, 56, 14, 17, 22, 29, 51, 87, 80, 62, 18, 22, 37, 56, 68, 109, 103, 77, 24, 35, 55, 64, 81, 104, 113, 92, 49, 64, 78, 87, 103, 121, 120, 101, 72, 92, 95, 98, 112, 100, 103, 99}; @@ -25,7 +27,7 @@ const std::array T_AC_0{0, 2, 1, 3, 3, 2, 4, 3, 5, const std::array DC_CAT_OFF{2, 3, 3, 3, 3, 3, 4, 5, 6, 7, 8, 9}; -MeteorImage::MeteorImage() +Image::Image() : mIsChannel64Available(false) , mIsChannel65Available(false) , mIsChannel66Available(false) @@ -39,23 +41,18 @@ MeteorImage::MeteorImage() initCos(); } -void MeteorImage::initCos() { +void Image::initCos() { for(int y = 0; y < 8; y++) { for(int x = 0; x < 8; x++) { mCosine[y][x] = cos(M_PI / 16 * (2 * y + 1) * x); } } for(int x = 0; x < 8; x++) { - if(x == 0) - mAlpha[x] = 1 / sqrt(2); - else - mAlpha[x] = 1; + mAlpha[x] = (x == 0) ? 1.0f / std::sqrt(2.0f) : 1.0f; } - - // for (auto x : cosine_) for (auto y: x) std::cout << std::dec << y << ' '; } -void MeteorImage::initHuffmanTable() { +void Image::initHuffmanTable() { std::array v{}; std::array min_code{}, maj_code{}; @@ -106,7 +103,7 @@ void MeteorImage::initHuffmanTable() { } } -int MeteorImage::getDcReal(uint16_t word) { +int Image::getDcReal(uint16_t word) { switch(word >> 14) { case 0: return 0; @@ -140,7 +137,7 @@ int MeteorImage::getDcReal(uint16_t word) { return -1; } -int MeteorImage::getAcReal(uint16_t word) { +int Image::getAcReal(uint16_t word) { for(int i = 0; i < 162; i++) { if(((word >> (16 - mAcTable[i].len)) & mAcTable[i].mask) == mAcTable[i].code) { return i; @@ -149,50 +146,27 @@ int MeteorImage::getAcReal(uint16_t word) { return -1; } -MeteorImage::~MeteorImage() {} +Image::~Image() {} -cv::Mat MeteorImage::getChannelImage(ChannelIDs APID, bool fillBlackLines) { - if(mFullImage.size() == 0) { +cv::Mat Image::getChannelImage(APIDs channelID, bool fillBlackLines) { + if(mChannels.size() == 0) { return cv::Mat(); } - const int width = 8 * MCU_PER_LINE; + const int width = 8 * cMCUPerLine; const int height = mCurY + 8; cv::Mat image = cv::Mat(height, width, CV_8UC3); for(int y = 0; y < image.rows; y++) { for(int x = 0; x < image.cols; x++) { - uint off = x + y * MCU_PER_LINE * 8; + uint off = x + y * cMCUPerLine * 8; cv::Vec3b& color = image.at(y, x); - image.at(cv::Point(x, y)) = color; - - switch(APID) { - case APID_68: - color[0] = mFullImage[off].r; - color[1] = mFullImage[off].r; - color[2] = mFullImage[off].r; - break; - case APID_66: - color[0] = mFullImage[off].r; - color[1] = mFullImage[off].r; - color[2] = mFullImage[off].r; - break; - case APID_65: - color[0] = mFullImage[off].g; - color[1] = mFullImage[off].g; - color[2] = mFullImage[off].g; - break; - case APID_64: - color[0] = mFullImage[off].b; - color[1] = mFullImage[off].b; - color[2] = mFullImage[off].b; - break; - } - - image.at(cv::Point(x, y)) = color; + color[0] = mChannels[off][channelID - 64]; + color[1] = mChannels[off][channelID - 64]; + color[2] = mChannels[off][channelID - 64]; } } @@ -204,24 +178,23 @@ cv::Mat MeteorImage::getChannelImage(ChannelIDs APID, bool fillBlackLines) { return image; } -cv::Mat MeteorImage::getRGBImage(ChannelIDs redAPID, ChannelIDs greenAPID, ChannelIDs blueAPID, bool fillBlackLines) { - if(mFullImage.size() == 0) { +cv::Mat Image::getRGBImage(APIDs redAPID, APIDs greenAPID, APIDs blueAPID, bool fillBlackLines, bool invertR, bool invertG, bool invertB) { + if(mChannels.size() == 0) { return cv::Mat(); } - int width = 8 * MCU_PER_LINE; + int width = 8 * cMCUPerLine; int height = mCurY + 8; cv::Mat image = cv::Mat(height, width, CV_8UC3); - cv::Vec3b color; for(int y = 0; y < image.rows; y++) { for(int x = 0; x < image.cols; x++) { - uint off = x + y * MCU_PER_LINE * 8; + uint off = x + y * cMCUPerLine * 8; + cv::Vec3b& color = image.at(cv::Point(x, y)); - color[0] = (blueAPID == APID_66 || blueAPID == APID_68) ? mFullImage[off].r : (blueAPID == APID_65) ? mFullImage[off].g : mFullImage[off].b; - color[1] = (greenAPID == APID_66 || greenAPID == APID_68) ? mFullImage[off].r : (greenAPID == APID_65) ? mFullImage[off].g : mFullImage[off].b; - color[2] = (redAPID == APID_66 || redAPID == APID_68) ? mFullImage[off].r : (redAPID == APID_65) ? mFullImage[off].g : mFullImage[off].b; - image.at(cv::Point(x, y)) = color; + color[0] = invertB ? ~mChannels[off][blueAPID - 64] : mChannels[off][blueAPID - 64]; + color[1] = invertG ? ~mChannels[off][greenAPID - 64] : mChannels[off][greenAPID - 64]; + color[2] = invertR ? ~mChannels[off][redAPID - 64] : mChannels[off][redAPID - 64]; } } @@ -233,62 +206,63 @@ cv::Mat MeteorImage::getRGBImage(ChannelIDs redAPID, ChannelIDs greenAPID, Chann return image; } -bool MeteorImage::progressImage(int apd, int mcu_id, int pck_cnt) { - if(apd == 0 || apd == 70) +bool Image::progressImage(int apd, int mcu_id, int pck_cnt) { + if(apd == 0 || apd == 70) { return false; + } if(mLastMCU == -1) { - if(mcu_id != 0) + if(mcu_id != 0) { return false; + } mPrevPacket = pck_cnt; mFirstPacket = pck_cnt; - if(apd == 65) + if(apd == 65) { mFirstPacket -= 14; - if(apd == 66) + } + if(apd == 66) { mFirstPacket -= 28; - if(apd == 68) + } + if(apd == 68) { mFirstPacket -= 28; + } mLastMCU = 0; mCurY = -1; } - if(pck_cnt < mPrevPacket) + if(pck_cnt < mPrevPacket) { mFirstPacket -= 16384; + } mPrevPacket = pck_cnt; mCurY = 8 * ((pck_cnt - mFirstPacket) / (14 + 14 + 14 + 1)); if(mCurY > mLastY) { - mFullImage.resize(MCU_PER_LINE * 8 * (mCurY + 8)); + mChannels.resize(cMCUPerLine * 8 * (mCurY + 8)); } mLastY = mCurY; return true; } -void MeteorImage::fillDqtByQ(std::array& dqt, int q) { +void Image::fillDqtByQ(std::array& dqt, int q) { float f; - if(q > 20 && q < 50) - f = 5000. / q; - else - f = 200. - 2 * q; + f = (q > 20 && q < 50) ? 5000.0 / q : 200.0 - 2.0 * q; for(int i = 0; i < 64; i++) { dqt[i] = (int)round(f / 100. * STANDARD_QUANTIZATION_TABLE[i]); - if(dqt[i] < 1) + if(dqt[i] < 1) { dqt[i] = 1; + } } } -int MeteorImage::mapRange(int cat, int vl) { +int Image::mapRange(int cat, int vl) { int maxval = (1 << cat) - 1; bool sig = (vl >> (cat - 1)) != 0; - if(sig) - return vl; - else - return vl - maxval; + return sig ? vl : vl - maxval; } -void MeteorImage::filtIdct8x8(std::array& res, std::array& inp) { +void Image::filtIdct8x8(std::array& res, std::array& inp) { for(int y = 0; y < 8; y++) { for(int x = 0; x < 8; x++) { float s = 0; @@ -303,70 +277,92 @@ void MeteorImage::filtIdct8x8(std::array& res, std::array& } } -void MeteorImage::fillPix(std::array& img_dct, int apd, int mcu_id, int m) { +void Image::fillPix(std::array& img_dct, int apd, int mcu_id, int m) { for(int i = 0; i < 64; i++) { - int t = round(img_dct[i] + 128); - if(t < 0) + int t = std::round(img_dct[i] + 128.0f); + if(t < 0) { t = 0; - if(t > 255) + } + if(t > 255) { t = 255; + } int x = (mcu_id + m) * 8 + i % 8; int y = mCurY + i / 8; - uint off = x + y * MCU_PER_LINE * 8; + uint off = x + y * cMCUPerLine * 8; - if(apd == APID_68 || apd == APID_66) { - mFullImage[off].r = t; - } - if(apd == APID_65) { - mFullImage[off].g = t; - } - if(apd == APID_64) { - mFullImage[off].b = t; + switch(apd) { + case 64: + mChannels[off][APIDs::APID64 - 64] = t; + break; + case 65: + mChannels[off][APIDs::APID65 - 64] = t; + break; + case 66: + mChannels[off][APIDs::APID66 - 64] = t; + break; + case 67: + mChannels[off][APIDs::APID67 - 64] = t; + break; + case 68: + mChannels[off][APIDs::APID68 - 64] = t; + break; + case 69: + mChannels[off][APIDs::APID69 - 64] = t; + break; + + default: + break; } } } -void MeteorImage::decMCUs(const uint8_t* packet, int len, int apd, int pck_cnt, int mcu_id, uint8_t q) { - BitIOConst b(packet); +void Image::decode(uint16_t apid, uint16_t packetCount, const Segment& segment) { + BitIOConst bitIO(segment.getPayloadData(), segment.getPayloadSize()); - if(!progressImage(apd, mcu_id, pck_cnt)) + if(!progressImage(apid, segment.getID(), packetCount)) { return; + } - if(apd == APID_64) { + if(apid == APIDs::APID64) { mIsChannel64Available = true; } - if(apd == APID_65) { + if(apid == APIDs::APID65) { mIsChannel65Available = true; } - - if(apd == APID_66) { + if(apid == APIDs::APID66) { mIsChannel66Available = true; } - if(apd == APID_68) { + if(apid == APIDs::APID67) { + mIsChannel67Available = true; + } + if(apid == APIDs::APID68) { mIsChannel68Available = true; } + if(apid == APIDs::APID69) { + mIsChannel69Available = true; + } std::array dqt{}; std::array zdct{}, dct{}, img_dct{}; - fillDqtByQ(dqt, q); + fillDqtByQ(dqt, segment.getQF()); - float prev_dc = 0; + float prevDC = 0; int m = 0; - while(m < MCU_PER_PACKET) { - int dc_cat = mDcLookup[b.peekBits(16)]; + while(m < cMCUPerPacket) { + int dc_cat = mDcLookup[bitIO.peekBits(16)]; if(dc_cat == -1) { std::cerr << "Bad DC Huffman code!" << std::endl; return; } - b.advanceBits(DC_CAT_OFF[dc_cat]); - uint32_t n = b.fetchBits(dc_cat); + bitIO.advanceBits(DC_CAT_OFF[dc_cat]); + uint32_t n = bitIO.fetchBits(dc_cat); - zdct[0] = mapRange(dc_cat, n) + prev_dc; - prev_dc = zdct[0]; + zdct[0] = mapRange(dc_cat, n) + prevDC; + prevDC = zdct[0]; int k = 1; while(k < 64) { - int ac = mAcLookup[b.peekBits(16)]; + int ac = mAcLookup[bitIO.peekBits(16)]; if(ac == -1) { std::cerr << "Bad AC Huffman code!" << std::endl; return; @@ -374,7 +370,7 @@ void MeteorImage::decMCUs(const uint8_t* packet, int len, int apd, int pck_cnt, int ac_len = mAcTable[ac].len; int ac_size = mAcTable[ac].size; int ac_run = mAcTable[ac].run; - b.advanceBits(ac_len); + bitIO.advanceBits(ac_len); if(ac_run == 0 && ac_size == 0) { for(int i = k; i < 64; i++) { @@ -389,7 +385,7 @@ void MeteorImage::decMCUs(const uint8_t* packet, int len, int apd, int pck_cnt, } if(ac_size != 0) { - uint16_t n = b.fetchBits(ac_size); + uint16_t n = bitIO.fetchBits(ac_size); zdct[k] = mapRange(ac_size, n); k++; } else if(ac_run == 15) { @@ -403,8 +399,13 @@ void MeteorImage::decMCUs(const uint8_t* packet, int len, int apd, int pck_cnt, } filtIdct8x8(img_dct, dct); - fillPix(img_dct, apd, mcu_id, m); + fillPix(img_dct, apid, segment.getID(), m); m++; } } + +} // namespace msumr +} // namespace lrpt +} // namespace protocol +} // namespace decoder \ No newline at end of file diff --git a/decoder/meteorimage.h b/decoder/protocol/lrpt/msumr/image.h similarity index 57% rename from decoder/meteorimage.h rename to decoder/protocol/lrpt/msumr/image.h index 9b337d1..2ab2847 100644 --- a/decoder/meteorimage.h +++ b/decoder/protocol/lrpt/msumr/image.h @@ -1,57 +1,37 @@ -#ifndef METEORIMAGE_H -#define METEORIMAGE_H +#pragma once #include #include #include #include +#include "segment.h" #include "threatimage.h" -// Based on: https://github.com/artlav/meteor_decoder/blob/master/met_jpg.pas - - -#ifdef __GNUC__ -#define PACK(__Declaration__) __Declaration__ __attribute__((__packed__)) -#endif - -#ifdef _MSC_VER -#define PACK(__Declaration__) __pragma(pack(push, 1)) __Declaration__ __pragma(pack(pop)) -#endif - -union Pixel { - PACK(struct { - uint8_t b; - uint8_t g; - uint8_t r; - uint8_t a; - }); - uint32_t pixel; -}; +namespace decoder { +namespace protocol { +namespace lrpt { +namespace msumr { + +class Image { + struct ac_table_rec { + int run; + int size; + int len; + uint32_t mask; + uint32_t code; + }; -struct ac_table_rec { - int run; - int size; - int len; - uint32_t mask; - uint32_t code; -}; -class MeteorImage { public: - enum ChannelIDs { - APID_68 = 68, // R - APID_66 = 66, // B - APID_65 = 65, // G - APID_64 = 64 // R - }; + enum APIDs { APID64 = 64, APID65, APID66, APID67, APID68, APID69 }; public: - MeteorImage(); - virtual ~MeteorImage(); + Image(); + virtual ~Image(); - cv::Mat getRGBImage(ChannelIDs redAPID, ChannelIDs greenAPID, ChannelIDs blueAPID, bool fillBlackLines = true); - cv::Mat getChannelImage(ChannelIDs APID, bool fillBlackLines = true); + cv::Mat getRGBImage(APIDs redAPID, APIDs greenAPID, APIDs blueAPID, bool fillBlackLines = true, bool invertR = false, bool invertG = false, bool invertB = false); + cv::Mat getChannelImage(APIDs APID, bool fillBlackLines = true); public: bool isChannel64Available() const { @@ -63,11 +43,18 @@ class MeteorImage { bool isChannel66Available() const { return mIsChannel66Available; } + bool isChannel67Available() const { + return mIsChannel67Available; + } bool isChannel68Available() const { return mIsChannel68Available; } + bool isChannel69Available() const { + return mIsChannel69Available; + } protected: + void decode(uint16_t apid, uint16_t packetCount, const Segment& segment); void decMCUs(const uint8_t* packet, int len, int apd, int pck_cnt, int mcu_id, uint8_t q); int getLastY() const { @@ -93,14 +80,24 @@ class MeteorImage { bool mIsChannel64Available; bool mIsChannel65Available; bool mIsChannel66Available; + bool mIsChannel67Available; bool mIsChannel68Available; + bool mIsChannel69Available; - std::vector mFullImage; + std::vector> mChannels; int mLastMCU, mCurY, mLastY, mFirstPacket, mPrevPacket; std::array mAcLookup{}, mDcLookup{}; std::array mAcTable{}; std::array, 8> mCosine{}; std::array mAlpha; + + private: + static constexpr uint16_t cMCUPerPacket = 14; + static constexpr uint16_t cMCUPerLine = 196; }; -#endif // METEORIMAGE_H + +} // namespace msumr +} // namespace lrpt +} // namespace protocol +} // namespace decoder diff --git a/decoder/protocol/lrpt/msumr/segment.cpp b/decoder/protocol/lrpt/msumr/segment.cpp new file mode 100644 index 0000000..568b621 --- /dev/null +++ b/decoder/protocol/lrpt/msumr/segment.cpp @@ -0,0 +1,28 @@ +#include "segment.h" + +namespace decoder { +namespace protocol { +namespace lrpt { +namespace msumr { + +Segment::Segment(const uint8_t* data, uint16_t size) + : mpData(data) + , mDataSize(size) { + + if(size >= 14) { + mDay = static_cast(data[0]) << 8 | data[1]; + mMsec = static_cast(data[2]) << 24 | static_cast(data[3]) << 16 | static_cast(data[4]) << 8 | data[5]; + mUsec = static_cast(data[6]) << 8 | data[7]; + mMCUN = data[8]; + mQT = data[9]; + mDC = data[10] & 0xF0 >> 4; + mAC = data[10] & 0x0F; + mQFM = static_cast(data[11]) << 8 | data[12]; + mQF = data[13]; + } +} + +} // namespace msumr +} // namespace lrpt +} // namespace protocol +} // namespace decoder \ No newline at end of file diff --git a/decoder/protocol/lrpt/msumr/segment.h b/decoder/protocol/lrpt/msumr/segment.h new file mode 100644 index 0000000..0c4edce --- /dev/null +++ b/decoder/protocol/lrpt/msumr/segment.h @@ -0,0 +1,62 @@ +#pragma once +#include + +namespace decoder { +namespace protocol { +namespace lrpt { +namespace msumr { + +class Segment { + public: + Segment(const uint8_t* data, uint16_t size); + + const uint8_t* getPayloadData() const { + return mpData + 14; + } + uint16_t getPayloadSize() const { + return mDataSize - 14; + } + uint8_t getID() const { + return mMCUN; + } + uint8_t getQF() const { + return mQF; + } + + uint16_t getDay() const { + return mDay; + } + + uint32_t getMsec() const { + return mMsec; + } + + uint16_t getUsec() const { + return mUsec; + } + + bool isValid() const { + // Note: These indexes are not used in meteor and QFM is always 0xFFF0 + return mQT == 0x00 && mDC == 0x00 && mDC == 0x00 && mQFM == 0xFFF0 && ((mDataSize - 14) > 0); + } + + + private: + const uint8_t* mpData; + uint16_t mDataSize; + uint16_t mDay; + uint32_t mMsec; + uint16_t mUsec; + + uint8_t mMCUN; + uint8_t mQT; + uint8_t mDC; + uint8_t mAC; + uint16_t mQFM; + uint8_t mQF; +}; + +} // namespace msumr +} // namespace lrpt +} // namespace protocol +} // namespace decoder \ No newline at end of file diff --git a/decoder/protocol/vcdu.cpp b/decoder/protocol/vcdu.cpp new file mode 100644 index 0000000..e69de29 diff --git a/decoder/protocol/vcdu.h b/decoder/protocol/vcdu.h new file mode 100644 index 0000000..a2109d5 --- /dev/null +++ b/decoder/protocol/vcdu.h @@ -0,0 +1,51 @@ +#pragma once + +#include +#include + +namespace decoder { +namespace protocol { + +class VCDU { + public: + VCDU(const uint8_t* cadu) { + mVersion = cadu[4] >> 6; + mSpaceCraftID = static_cast(cadu[4] & 0b111111) << 2 | cadu[5] >> 6; + mVCID = cadu[5] & 0x3F; + mCounter = static_cast(cadu[6]) << 16 | static_cast(cadu[7]) << 8 | cadu[8]; + mReplayFlag = cadu[9] >> 7; + mInsertZone = static_cast(cadu[10]) << 16 | cadu[11]; + } + + public: + uint8_t getVersion() const { + return mVersion; + } + uint16_t getSpaceCraftID() const { + return mSpaceCraftID; + } + uint8_t getVCID() const { + return mVCID; + } + uint32_t getCounter() const { + return mCounter; + } + bool getReplayFlag() const { + return mReplayFlag; + } + + static constexpr uint32_t size() { + return 8; + } + + private: + uint8_t mVersion; + uint16_t mSpaceCraftID; + uint8_t mVCID; + uint32_t mCounter; + bool mReplayFlag; + uint16_t mInsertZone; +}; + +} // namespace protocol +} // namespace decoder \ No newline at end of file diff --git a/decoder/reedsolomon.cpp b/decoder/reedsolomon.cpp index c79713c..01328a2 100644 --- a/decoder/reedsolomon.cpp +++ b/decoder/reedsolomon.cpp @@ -24,5 +24,15 @@ void ReedSolomon::interleave(uint8_t* output, int pos, int n) { } int ReedSolomon::decode() { - return correct_reed_solomon_decode(mpReedSolomon, mWorkBuffer, sizeof(mWorkBuffer), mResultBuffer); + int result = correct_reed_solomon_decode(mpReedSolomon, mWorkBuffer, sizeof(mWorkBuffer), mResultBuffer); + if(result != -1) { + result = 0; + for(int i = 0; i < 223; i++) { + if(mWorkBuffer[i] != mResultBuffer[i]) { + result++; + } + } + } + + return result; } diff --git a/decoder/viterbi.cpp b/decoder/viterbi.cpp index c130393..17af963 100644 --- a/decoder/viterbi.cpp +++ b/decoder/viterbi.cpp @@ -15,5 +15,38 @@ Viterbi::~Viterbi() { } size_t Viterbi::decodeSoft(const uint8_t* data, uint8_t* result, size_t blockSize) { - return correct_convolutional_decode_soft(mpConvolutional, data, blockSize, result); + if(mpConvolutional == nullptr) { + return -1; + } + + size_t reencodedSize = blockSize / 8; + if(mReencodedBuffer.size() != reencodedSize) { + mReencodedBuffer.resize(reencodedSize); + } + + size_t decodedBytes = correct_convolutional_decode_soft(mpConvolutional, data, blockSize, result); + if(decodedBytes != -1) { + correct_convolutional_encode(mpConvolutional, result, decodedBytes, mReencodedBuffer.data()); + calculateBer(data, mReencodedBuffer.data(), blockSize); + } + return decodedBytes; +} + +void Viterbi::calculateBer(const uint8_t* original, const uint8_t* reEncoded, ssize_t blockSize) { + float errors = 0, total = 0; + for(int i = 0; i < blockSize / 8; i++) { + for(int bit = 7; bit >= 0; bit--) { + if(*original != 128) { + errors += (*original > 128) != ((reEncoded[i] >> bit) & 0x01); + total++; + } + original++; + } + } + + if(total == 0) { + mLastBER = 0.0f; + } else { + mLastBER = (errors / total); + } } diff --git a/decoder/viterbi.h b/decoder/viterbi.h index d123b5c..320b7ed 100644 --- a/decoder/viterbi.h +++ b/decoder/viterbi.h @@ -2,6 +2,8 @@ #define VITERBI_H #include + +#include extern "C" { #include "correct.h" } @@ -13,9 +15,18 @@ class Viterbi { size_t decodeSoft(const uint8_t* data, uint8_t* result, size_t blockSize); + float getLastBER() { + return mLastBER; + } + + private: + void calculateBer(const uint8_t* original, const uint8_t* reEncoded, ssize_t blockSize); + private: correct_convolutional_polynomial_t mPolynomials[2]; correct_convolutional* mpConvolutional; + std::vector mReencodedBuffer; + float mLastBER = 0.0f; }; #endif // VITERBI_H diff --git a/imageproc/blendimages.cpp b/imageproc/blendimages.cpp new file mode 100644 index 0000000..0b0ff24 --- /dev/null +++ b/imageproc/blendimages.cpp @@ -0,0 +1,110 @@ +#include "blendimages.h" + +cv::Mat BlendImages::merge(std::list& images) { + std::list::iterator it = images.begin(); + + cv::Mat composite = it->clone(); + ++it; + + composite.convertTo(composite, CV_32FC3); + + for(; it != images.end(); ++it) { + it->convertTo(*it, CV_32FC3); + + cv::Mat grayScale1; + cv::Mat alpha1; + cv::medianBlur(composite, grayScale1, 5); + cv::cvtColor(grayScale1, grayScale1, cv::COLOR_BGR2GRAY); + + cv::threshold(grayScale1, alpha1, 0, 255, cv::THRESH_BINARY); + grayScale1.release(); + + cv::Mat grayScale2; + cv::Mat alpha2; + cv::medianBlur(*it, grayScale2, 5); + cv::cvtColor(grayScale2, grayScale2, cv::COLOR_BGR2GRAY); + + cv::threshold(grayScale2, alpha2, 0, 255, cv::THRESH_BINARY); + grayScale2.release(); + + cv::Mat mask; + cv::bitwise_and(alpha1, alpha2, mask); + alpha1.release(); + alpha2.release(); + + std::vector channels; + channels.push_back(mask); + channels.push_back(mask); + channels.push_back(mask); + cv::merge(channels, mask); + + mask.convertTo(mask, CV_32FC3, 1 / 255.0); + + int start0 = findImageStart(composite); + int start1 = findImageStart(*it); + bool leftToRight = start0 < start1; + + cv::Mat blendmask = blendMask(mask, leftToRight); + cv::multiply(cv::Scalar::all(1.0) - blendmask, composite, composite); + blendmask = blendMask(mask, !leftToRight); + cv::multiply(cv::Scalar::all(1.0) - blendmask, *it, *it); + + cv::add(composite, *it, composite); + } + + return composite; +} + +int BlendImages::findImageStart(const cv::Mat& img) { + int i = img.size().width; + for(int y = 0; y < img.size().height; ++y) { + for(int x = 0; x < img.size().width; x++) { + if(img.at(y, x) != cv::Vec3f(0, 0, 0)) { + if(x < i) { + i = x; + } + } + } + } + return i; +} + +cv::Mat BlendImages::blendMask(const cv::Mat& mask, bool leftToRight) { + int xStart = 0; + int xEnd = 0; + int blendWidth = 0; + float alpha; + + cv::Mat blendedMask = cv::Mat::zeros(mask.size().height, mask.size().width, mask.type()); + + for(int y = 0; y < mask.size().height; ++y) { + bool foundStart = false; + bool foundEnd = false; + for(int x = 0; x < mask.size().width; x++) { + if(!foundStart && mask.at(y, x) != cv::Vec3f(0, 0, 0)) { + foundStart = true; + xStart = x; + + for(int temp = x; temp < mask.size().width; temp++) { + if(mask.at(y, temp) != cv::Vec3f(1, 1, 1)) { + xEnd = temp; + foundEnd = true; + blendWidth = xEnd - xStart; + break; + } + } + } + + if(foundStart && foundEnd && (x >= xStart && x <= xEnd)) { + alpha = static_cast(x - xStart) / blendWidth; + alpha = leftToRight ? alpha : 1.0f - alpha; + + blendedMask.at(y, x) = mask.at(y, x) * alpha; + } else if(foundStart && foundEnd && (x > xEnd)) { + foundStart = false; + foundEnd = false; + } + } + } + return blendedMask; +} \ No newline at end of file diff --git a/imageproc/blendimages.h b/imageproc/blendimages.h new file mode 100644 index 0000000..68cbbaf --- /dev/null +++ b/imageproc/blendimages.h @@ -0,0 +1,13 @@ +#pragma once + +#include +#include + +class BlendImages { + public: + static cv::Mat merge(std::list& images); + + private: + static int findImageStart(const cv::Mat& img); + static cv::Mat blendMask(const cv::Mat& mask, bool leftToRight); +}; \ No newline at end of file diff --git a/imageproc/projectimage.cpp b/imageproc/projectimage.cpp new file mode 100644 index 0000000..0f6ccf4 --- /dev/null +++ b/imageproc/projectimage.cpp @@ -0,0 +1,377 @@ +#include "projectimage.h" + +#include +#include +#include +#include + +#include "GIS/shaperenderer.h" +#include "settings.h" +#include "threadpool.h" + +std::map ProjectImage::MarkerLookup{{"STAR", cv::MARKER_STAR}, + {"CROSS", cv::MARKER_CROSS}, + {"SQUARE", cv::MARKER_SQUARE}, + {"DIAMOND", cv::MARKER_DIAMOND}, + {"TRIANGLE_UP", cv::MARKER_TRIANGLE_UP}, + {"TRIANGLE_DOWN", cv::MARKER_TRIANGLE_DOWN}, + {"TILTED_CROSS", cv::MARKER_TILTED_CROSS}}; + +std::list ProjectImage::createCompositeProjector(Projection projection, const std::list& gcpCalclulators, float scale, int earthRadius, int altitude) { + float MinX = std::numeric_limits::max(); + float MinY = std::numeric_limits::max(); + float MaxX = std::numeric_limits::min(); + float MaxY = std::numeric_limits::min(); + float corner; + + Vector3 centerVector; + for(const auto& c : gcpCalclulators) { + centerVector += Vector3(c.getCenterCoordinate()); + } + CoordGeodetic center = centerVector.toCoordinate(); + + if(projection == Projection::Equidistant) { + for(const auto& gcpCalc : gcpCalclulators) { + auto topLeft = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(gcpCalc.getCoordinateTopLeft(), center, gcpCalc.getSatelliteHeight(), scale); + auto topRight = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(gcpCalc.getCoordinateTopRight(), center, gcpCalc.getSatelliteHeight(), scale); + auto bottomLeft = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(gcpCalc.getCoordinateBottomLeft(), center, gcpCalc.getSatelliteHeight(), scale); + auto bottomRight = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(gcpCalc.getCoordinateBottomRight(), center, gcpCalc.getSatelliteHeight(), scale); + + corner = std::min({topLeft.x, topRight.x, bottomLeft.x, bottomRight.x}); + MinX = corner < MinX ? corner : MinX; + + corner = std::min({topLeft.y, topRight.y, bottomLeft.y, bottomRight.y}); + MinY = corner < MinY ? corner : MinY; + + corner = std::max({topLeft.x, topRight.x, bottomLeft.x, bottomRight.x}); + MaxX = corner > MaxX ? corner : MaxX; + + corner = std::max({topLeft.y, topRight.y, bottomLeft.y, bottomRight.y}); + MaxY = corner > MaxY ? corner : MaxY; + } + } else if(projection == Projection::Mercator) { + for(const auto& gcpCalc : gcpCalclulators) { + PixelGeolocationCalculator::CartesianCoordinateF topLeft = PixelGeolocationCalculator::coordinateToMercatorProjection(gcpCalc.getCoordinateTopLeft(), gcpCalc.getSatelliteHeight(), scale); + PixelGeolocationCalculator::CartesianCoordinateF topRight = PixelGeolocationCalculator::coordinateToMercatorProjection(gcpCalc.getCoordinateTopRight(), gcpCalc.getSatelliteHeight(), scale); + PixelGeolocationCalculator::CartesianCoordinateF bottomLeft = PixelGeolocationCalculator::coordinateToMercatorProjection(gcpCalc.getCoordinateBottomLeft(), gcpCalc.getSatelliteHeight(), scale); + PixelGeolocationCalculator::CartesianCoordinateF bottomRight = PixelGeolocationCalculator::coordinateToMercatorProjection(gcpCalc.getCoordinateBottomRight(), gcpCalc.getSatelliteHeight(), scale); + + corner = std::min({topLeft.x, topRight.x, bottomLeft.x, bottomRight.x}); + MinX = corner < MinX ? corner : MinX; + + corner = std::min({topLeft.y, topRight.y, bottomLeft.y, bottomRight.y}); + MinY = corner < MinY ? corner : MinY; + + corner = std::max({topLeft.x, topRight.x, bottomLeft.x, bottomRight.x}); + MaxX = corner > MaxX ? corner : MaxX; + + corner = std::max({topLeft.y, topRight.y, bottomLeft.y, bottomRight.y}); + MaxY = corner > MaxY ? corner : MaxY; + } + } + + int width = static_cast(std::abs(MaxX - MinX)); + int height = static_cast(std::abs(MaxY - MinY)); + float xStart = std::min(MinX, MaxX); + float yStart = std::min(MinY, MaxY); + + std::list projectors; + for(const auto& gcpCalc : gcpCalclulators) { + ProjectImage projector(projection, gcpCalc, scale, earthRadius, altitude); + projector.mWidth = width; + projector.mHeight = height; + projector.mXStart = xStart; + projector.mYStart = yStart; + projector.mCenterCoordinate = center; + projector.mBoundariesCalcNeeded = false; + projectors.emplace_back(projector); + } + return projectors; +} + +ProjectImage::ProjectImage(Projection projection, const PixelGeolocationCalculator& geolocationCalculator, float scale, int earthRadius, int altitude) + : mProjection(projection) + , mGeolocationCalculator(geolocationCalculator) + , mScale(scale) + , mEarthRadius(earthRadius) + , mAltitude(altitude) + , mTransformer(cv::createThinPlateSplineShapeTransformer2(Settings::getInstance().getResourcesPath() + "kernels/tps.cl")) { + + mTheta = 0.5 * SWATH / earthRadius; // Maximum half-Angle subtended by Swath from Earth's centre + mHalfChord = static_cast(earthRadius * std::sin(mTheta)); // Maximum Length of chord subtended at Centre of Earth + int h = static_cast(earthRadius * std::cos(mTheta)); // Minimum Altitude from Centre of Earth to chord from std::costd::sine + mInc = earthRadius - h; // Maximum Distance from Arc to Chord + mScanAngle = std::atan(mHalfChord / static_cast(altitude + mInc)); // Maximum Angle subtended by half-chord from satellite +} + +void ProjectImage::calculateTransformation(const cv::Size& imageSize) { + if(mBoundariesCalcNeeded) { + calculateImageBoundaries(); + } + + std::vector sourcePoints, targetPoints; + for(int y = 0; y < imageSize.height; y += 100) { + for(int x = 0; x < imageSize.width; x += 100) { + if(mProjection == Projection::Equidistant || mProjection == Projection::Rectify) { + const auto p1 = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(mGeolocationCalculator.getCoordinateAt(x, y), mCenterCoordinate, mGeolocationCalculator.getSatelliteHeight(), mScale); + sourcePoints.push_back(cv::Point2f(x, y)); + targetPoints.push_back(cv::Point2f(p1.x + (-mXStart), p1.y + (-mYStart))); + } else { + // Mercator + const auto p1 = PixelGeolocationCalculator::coordinateToMercatorProjection(mGeolocationCalculator.getCoordinateAt(x, y), mGeolocationCalculator.getSatelliteHeight(), mScale); + sourcePoints.push_back(cv::Point2f(x, y)); + targetPoints.push_back(cv::Point2f(p1.x + (-mXStart), p1.y + (-mYStart))); + } + } + + int x = imageSize.width; + if(mProjection == Projection::Equidistant || mProjection == Projection::Rectify) { + const auto p1 = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(mGeolocationCalculator.getCoordinateAt(x, y), mCenterCoordinate, mGeolocationCalculator.getSatelliteHeight(), mScale); + sourcePoints.push_back(cv::Point2f(x, y)); + targetPoints.push_back(cv::Point2f(p1.x + (-mXStart), p1.y + (-mYStart))); + } else { + // Mercator + const auto p1 = PixelGeolocationCalculator::coordinateToMercatorProjection(mGeolocationCalculator.getCoordinateAt(x, y), mGeolocationCalculator.getSatelliteHeight(), mScale); + sourcePoints.push_back(cv::Point2f(x, y)); + targetPoints.push_back(cv::Point2f(p1.x + (-mXStart), p1.y + (-mYStart))); + } + } + + std::vector matches; + for(std::size_t i = 0; i < sourcePoints.size(); i++) { + matches.push_back(cv::DMatch(i, i, 0)); + } + mTransformer->estimateTransformation(targetPoints, sourcePoints, matches); + + if(mProjection == Projection::Rectify) { + rectify(imageSize); + } else { + mMapX = cv::Mat(mHeight, mWidth, CV_32FC1); + mMapY = cv::Mat(mHeight, mWidth, CV_32FC1); + + std::vector dst; + std::vector src(mHeight * mWidth); + + for(int row = 0; row < mHeight; row++) { + for(int col = 0; col < mWidth; col++) { + src[col + (row * mWidth)] = cv::Point2f(col, row); + } + } + mTransformer->applyTransformation(src, dst); + for(int row = 0; row < mHeight; row++) { + for(int col = 0; col < mWidth; col++) { + mMapX.at(row, col) = dst[col + (row * mWidth)].x; + mMapY.at(row, col) = dst[col + (row * mWidth)].y; + } + } + } +} + +cv::Mat ProjectImage::project(const cv::Mat& image) { + cv::Mat newImage; + cv::remap(image, newImage, mMapX, mMapY, cv::INTER_LINEAR); + + drawMapOverlay(newImage); + + return newImage; +} + +void ProjectImage::calculateImageBoundaries() { + float minX; + float minY; + float maxX; + float maxY; + + mCenterCoordinate = mGeolocationCalculator.getCenterCoordinate(); + + if(mProjection == Projection::Equidistant || mProjection == Projection::Rectify) { + auto topLeft = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(mGeolocationCalculator.getCoordinateTopLeft(), mCenterCoordinate, mGeolocationCalculator.getSatelliteHeight(), mScale); + auto topRight = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(mGeolocationCalculator.getCoordinateTopRight(), mCenterCoordinate, mGeolocationCalculator.getSatelliteHeight(), mScale); + auto bottomLeft = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(mGeolocationCalculator.getCoordinateBottomLeft(), mCenterCoordinate, mGeolocationCalculator.getSatelliteHeight(), mScale); + auto bottomRight = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(mGeolocationCalculator.getCoordinateBottomRight(), mCenterCoordinate, mGeolocationCalculator.getSatelliteHeight(), mScale); + + minX = std::min({topLeft.x, topRight.x, bottomLeft.x, bottomRight.x}); + minY = std::min({topLeft.y, topRight.y, bottomLeft.y, bottomRight.y}); + maxX = std::max({topLeft.x, topRight.x, bottomLeft.x, bottomRight.x}); + maxY = std::max({topLeft.y, topRight.y, bottomLeft.y, bottomRight.y}); + } else if(mProjection == Projection::Mercator) { + auto topLeft = PixelGeolocationCalculator::coordinateToMercatorProjection(mGeolocationCalculator.getCoordinateTopLeft(), mGeolocationCalculator.getSatelliteHeight(), mScale); + auto topRight = PixelGeolocationCalculator::coordinateToMercatorProjection(mGeolocationCalculator.getCoordinateTopRight(), mGeolocationCalculator.getSatelliteHeight(), mScale); + auto bottomLeft = PixelGeolocationCalculator::coordinateToMercatorProjection(mGeolocationCalculator.getCoordinateBottomLeft(), mGeolocationCalculator.getSatelliteHeight(), mScale); + auto bottomRight = PixelGeolocationCalculator::coordinateToMercatorProjection(mGeolocationCalculator.getCoordinateBottomRight(), mGeolocationCalculator.getSatelliteHeight(), mScale); + + minX = std::min({topLeft.x, topRight.x, bottomLeft.x, bottomRight.x}); + minY = std::min({topLeft.y, topRight.y, bottomLeft.y, bottomRight.y}); + maxX = std::max({topLeft.x, topRight.x, bottomLeft.x, bottomRight.x}); + maxY = std::max({topLeft.y, topRight.y, bottomLeft.y, bottomRight.y}); + } else { + } + + mWidth = static_cast(std::abs(maxX - minX)); + mHeight = static_cast(std::abs(maxY - minY)); + mXStart = static_cast(std::min(minX, maxX)); + mYStart = static_cast(std::min(minY, maxY)); +} + +void ProjectImage::rectify(const cv::Size& imageSize) { + int i, k, ni; + int DupPix; + float phi; + + mWidth = static_cast(imageSize.width * (0.902 + std::sin(mScanAngle))); + mHeight = imageSize.height; + mOriginalImageHalfWidth = (imageSize.width + 1) / 2; // Half-width of original image ('A') + mNewImageHalfWidth = (mWidth + 1) / 2; // Half-Width of stretched image ('Z') + + // Create lookup table of std::sines + std::unique_ptr lut = std::make_unique(mNewImageHalfWidth); + + for(i = 0; i < mNewImageHalfWidth - 1; i++) { + phi = std::atan((i / static_cast(mNewImageHalfWidth)) * mHalfChord / (mAltitude + mInc)); + lut[i] = static_cast(mOriginalImageHalfWidth * (std::sin(phi) / std::sin(mScanAngle))); + } + + mMapX = cv::Mat(mHeight, mWidth, CV_32FC1); + mMapY = cv::Mat(mHeight, mWidth, CV_32FC1); + + mFlip = mGeolocationCalculator.isNorthBoundPass(); + + i = 0; + ni = 0; + while(i < mNewImageHalfWidth) { + k = lut[i]; + ni = i; + while(true) { + if(lut[ni] == 0 || lut[ni] > k) + break; + ni += 1; + } + + DupPix = ni - i; // DupPix = number of repeated (spread) pixels + + if(lut[ni] == 0) { + DupPix = 1; + } + + for(int y = 0; y < mHeight; y++) { + for(int x = 0; x < DupPix; x++) { + mMapX.at(y, mFlip ? mWidth - (mNewImageHalfWidth + i - 1 + x) : mNewImageHalfWidth + i - 1 + x) = mOriginalImageHalfWidth + k - 1; + mMapX.at(y, mFlip ? mWidth - (mNewImageHalfWidth - i - x) : mNewImageHalfWidth - i - x) = mOriginalImageHalfWidth - k; + mMapY.at(y, mNewImageHalfWidth + i - 1 + x) = mFlip ? (mHeight - y) : y; + mMapY.at(y, mNewImageHalfWidth - i - x) = mFlip ? (mHeight - y) : y; + } + } + i += DupPix; + } +} + +void ProjectImage::drawMapOverlay(cv::Mat& image) { + Settings& settings = Settings::getInstance(); + + GIS::ShapeRenderer graticules(settings.getResourcesPath() + settings.getShapeGraticulesFile(), cv::Scalar(settings.getShapeGraticulesColor().B, settings.getShapeGraticulesColor().G, settings.getShapeGraticulesColor().R)); + graticules.setThickness(settings.getShapeGraticulesThickness()); + graticules.drawShape(image, std::bind(&ProjectImage::transform, this, std::placeholders::_1, std::placeholders::_2)); + + GIS::ShapeRenderer countryBorders(settings.getResourcesPath() + settings.getShapeBoundaryLinesFile(), + cv::Scalar(settings.getShapeBoundaryLinesColor().B, settings.getShapeBoundaryLinesColor().G, settings.getShapeBoundaryLinesColor().R)); + countryBorders.setThickness(settings.getShapeBoundaryLinesThickness()); + countryBorders.drawShape(image, std::bind(&ProjectImage::transform, this, std::placeholders::_1, std::placeholders::_2)); + + GIS::ShapeRenderer coastLines(settings.getResourcesPath() + settings.getShapeCoastLinesFile(), cv::Scalar(settings.getShapeCoastLinesColor().B, settings.getShapeCoastLinesColor().G, settings.getShapeCoastLinesColor().R)); + coastLines.setThickness(settings.getShapeCoastLinesThickness()); + coastLines.drawShape(image, std::bind(&ProjectImage::transform, this, std::placeholders::_1, std::placeholders::_2)); + + GIS::ShapeRenderer cities(settings.getResourcesPath() + settings.getShapePopulatedPlacesFile(), + cv::Scalar(settings.getShapePopulatedPlacesColor().B, settings.getShapePopulatedPlacesColor().G, settings.getShapePopulatedPlacesColor().R)); + cities.setFontHeight(settings.getShapePopulatedPlacesFontSize() * mScale); + cities.setFontLineWidth(settings.getShapePopulatedPlacesFontWidth()); + cities.setPointRadius(settings.getShapePopulatedPlacesPointradius() * mScale); + cities.addNumericFilter(settings.getShapePopulatedPlacesFilterColumnName(), settings.getShapePopulatedPlacesNumbericFilter()); + cities.setTextFieldName(settings.getShapePopulatedPlacesTextColumnName()); + cities.drawShape(image, std::bind(&ProjectImage::transform, this, std::placeholders::_1, std::placeholders::_2)); + + if(settings.drawReceiver()) { + double x = settings.getReceiverLatitude(); + double y = settings.getReceiverLongitude(); + bool draw = transform(x, y); + + if(draw) { + cv::drawMarker(image, cv::Point2d(x, y), cv::Scalar(0, 0, 0), stringToMarkerType(settings.getReceiverMarkType()), settings.getReceiverSize(), settings.getReceiverThickness() + 1); + cv::drawMarker(image, + cv::Point2d(x, y), + cv::Scalar(settings.getReceiverColor().B, settings.getReceiverColor().G, settings.getReceiverColor().R), + stringToMarkerType(settings.getReceiverMarkType()), + settings.getReceiverSize(), + settings.getReceiverThickness()); + } + } +} + +cv::MarkerTypes ProjectImage::stringToMarkerType(const std::string& markerType) { + auto itr = MarkerLookup.find(markerType); + if(itr != MarkerLookup.end()) { + return itr->second; + } + return cv::MARKER_TRIANGLE_UP; +} + +bool ProjectImage::transform(double& x, double& y) { + bool valid = true; + float centerLatitude = static_cast(mGeolocationCalculator.getCenterCoordinate().latitude * (180.0 / M_PI)); + float centerLongitude = static_cast(mGeolocationCalculator.getCenterCoordinate().longitude * (180.0 / M_PI)); + + if(mProjection == Projection::Equidistant) { + auto location = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection({x, y, 0}, mCenterCoordinate, mGeolocationCalculator.getSatelliteHeight(), mScale); + valid = equidistantCheck(x, y, centerLatitude, centerLongitude); + x = location.x + (-mXStart); + y = location.y + (-mYStart); + } else if(mProjection == Projection::Mercator) { + auto location = PixelGeolocationCalculator::coordinateToMercatorProjection({x, y, 0}, mGeolocationCalculator.getSatelliteHeight(), mScale); + x = location.x + (-mXStart); + y = location.y + (-mYStart); + } else if(mProjection == Projection::Rectify) { + auto location = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection({x, y, 0}, mCenterCoordinate, mGeolocationCalculator.getSatelliteHeight(), mScale); + valid = equidistantCheck(x, y, centerLatitude, centerLongitude); + x = location.x + (-mXStart); + y = location.y + (-mYStart); + + std::vector dst; + std::vector src; + src.emplace_back(cv::Point2f(x, y)); + mTransformer->applyTransformation(src, dst); + + if(dst[0].x < 0 || dst[0].x > mOriginalImageHalfWidth * 2) { + valid = false; + } + + if(dst[0].x >= mOriginalImageHalfWidth) { + int k = std::floor(dst[0].x) - mOriginalImageHalfWidth; + double phi = std::asin(static_cast(k) / mOriginalImageHalfWidth * std::sin(mScanAngle)); + double original = std::tan(phi) * (mAltitude + mInc) / mHalfChord * mNewImageHalfWidth; + x = mFlip ? mWidth - (mNewImageHalfWidth + original) : mNewImageHalfWidth + original; + } else { + int k = mOriginalImageHalfWidth - std::floor(dst[0].x); + double phi = std::asin(static_cast(k) / mOriginalImageHalfWidth * std::sin(mScanAngle)); + double original = std::tan(phi) * (mAltitude + mInc) / mHalfChord * mNewImageHalfWidth; + x = mFlip ? mWidth - (mNewImageHalfWidth - original) : mNewImageHalfWidth - original; + } + + y = mFlip ? mHeight - dst[0].y : dst[0].y; + } + return valid; +} + +bool ProjectImage::equidistantCheck(float latitude, float longitude, float centerLatitude, float centerLongitude) { + // Degree To radian + latitude = M_PI * latitude / 180.0f; + longitude = M_PI * longitude / 180.0f; + centerLatitude = M_PI * centerLatitude / 180.0f; + centerLongitude = M_PI * centerLongitude / 180.0f; + + float deltaSigma = std::sin(centerLatitude) * std::sin(latitude) + std::cos(latitude) * std::cos(longitude - centerLongitude); + if(deltaSigma < 0.0) { + return false; + } + + return true; +} diff --git a/imageproc/projectimage.h b/imageproc/projectimage.h new file mode 100644 index 0000000..b8367aa --- /dev/null +++ b/imageproc/projectimage.h @@ -0,0 +1,59 @@ +#pragma once + +#include +#include +#include +//#include +#include + +#include "pixelgeolocationcalculator.h" + +class ProjectImage { + public: + enum Projection { Equidistant, Mercator, Rectify }; + + public: + static std::list createCompositeProjector(Projection projection, const std::list& gcpCalclulators, float scale, int earthRadius = 6378, int altitude = 825); + + public: + ProjectImage(Projection projection, const PixelGeolocationCalculator& geolocationCalculator, float scale, int earthRadius = 6378, int altitude = 825); + + void calculateTransformation(const cv::Size& imageSize); + cv::Mat project(const cv::Mat& image); + + private: + void calculateImageBoundaries(); + void rectify(const cv::Size& imageSize); + void drawMapOverlay(cv::Mat& image); + cv::MarkerTypes stringToMarkerType(const std::string& markerType); + bool transform(double& x, double& y); + bool equidistantCheck(float latitude, float longitude, float centerLatitude, float centerLongitude); + + protected: + Projection mProjection; + const PixelGeolocationCalculator& mGeolocationCalculator; + float mScale; + int mEarthRadius; + int mAltitude; // Meteor MN2 orbit altitude + float mTheta; + int mInc; + int mHalfChord; + float mScanAngle; + int mNewImageHalfWidth; + int mOriginalImageHalfWidth; + CoordGeodetic mCenterCoordinate; + int mWidth; + int mHeight; + int mXStart; + int mYStart; + bool mBoundariesCalcNeeded = true; + cv::Mat mMapX; + cv::Mat mMapY; + bool mFlip = false; + + cv::Ptr mTransformer; + + static std::map MarkerLookup; + + static constexpr int SWATH = 2800; // Meteor M2M swath width +}; diff --git a/imageproc/spreadimage.cpp b/imageproc/spreadimage.cpp index cdf2de7..75d1af8 100644 --- a/imageproc/spreadimage.cpp +++ b/imageproc/spreadimage.cpp @@ -189,93 +189,6 @@ cv::Mat SpreadImage::mercatorProjection(const cv::Mat& image, const PixelGeoloca return newImage; } -// Concept for using ThinPlateSplineShapeTransform. Unfortunately it is extremly slow for big images -/*cv::Mat SpreadImage::equidistantProjection(const cv::Mat &image, const PixelGeolocationCalculator &geolocationCalculator, float scale, ProgressCallback progressCallback) -{ - cv::Point2f srcTri[3]; - cv::Point2f dstTri[3]; - - const CoordGeodetic ¢er = geolocationCalculator.getCenterCoordinate(); - - PixelGeolocationCalculator::CartesianCoordinateF topLeft = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(geolocationCalculator.getCoordinateTopLeft(), center, geolocationCalculator.getSatelliteHeight(), -scale); PixelGeolocationCalculator::CartesianCoordinateF topRight = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(geolocationCalculator.getCoordinateTopRight(), center, -geolocationCalculator.getSatelliteHeight(), scale); PixelGeolocationCalculator::CartesianCoordinateF bottomLeft = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(geolocationCalculator.getCoordinateBottomLeft(), -center, geolocationCalculator.getSatelliteHeight(), scale); PixelGeolocationCalculator::CartesianCoordinateF bottomRight = -PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(geolocationCalculator.getCoordinateBottomRight(), center, geolocationCalculator.getSatelliteHeight(), scale); - - double MinX = std::min(topLeft.x, std::min(topRight.x, std::min(bottomLeft.x, bottomRight.x))); - double MinY = std::min(topLeft.y, std::min(topRight.y, std::min(bottomLeft.y, bottomRight.y))); - double MaxX = std::max(topLeft.x, std::max(topRight.x, std::max(bottomLeft.x, bottomRight.x))); - double MaxY = std::max(topLeft.y, std::max(topRight.y, std::max(bottomLeft.y, bottomRight.y))); - - int width = static_cast(std::abs(MaxX - MinX)); - int height = static_cast(std::abs(MaxY - MinY)); - - int xStart = static_cast(std::min(MinX, MaxX)); - int yStart = static_cast(std::min(MinY, MaxY)); - - int imageWithGeorefHeight = geolocationCalculator.getGeorefMaxImageHeight() < image.size().height ? geolocationCalculator.getGeorefMaxImageHeight() : image.size().height; - - cv::Mat newImage = cv::Mat::zeros(height, width, image.type()); - //cv::Mat paddedImage = cv::Mat::zeros(height, width, image.type()); - //image.copyTo(paddedImage.rowRange(0, image.size().height).colRange(0, image.size().width)); - - auto tpsTransform = cv::createAffineTransformer(true); - std::vector sourcePoints, targetPoints; - - for (int y = 0; y < imageWithGeorefHeight; y += 50) { - for (int x = 0; x < image.size().width; x += 50) { - const PixelGeolocationCalculator::CartesianCoordinateF p1 = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(geolocationCalculator.getCoordinateAt(x, y), center, -geolocationCalculator.getSatelliteHeight(), scale); - - sourcePoints.push_back(cv::Point2f(x, y)); - targetPoints.push_back(cv::Point2f((int)(p1.x + (-xStart)), (int)(p1.y + (-yStart)))); - } - } - - std::vector matches; - for (unsigned int i = 0; i < sourcePoints.size(); i++) { - matches.push_back(cv::DMatch(i, i, 0)); - } - //tpsTransform->estimateTransformation(sourcePoints, targetPoints, matches); - tpsTransform->estimateTransformation(targetPoints, sourcePoints, matches); - std::vector transPoints; - tpsTransform->applyTransformation(sourcePoints, transPoints); - tpsTransform->warpImage(image, newImage, cv::INTER_LINEAR); - - float centerLatitude = static_cast(geolocationCalculator.getCenterCoordinate().latitude * (180.0 / M_PI)); - float centerLongitude = static_cast(geolocationCalculator.getCenterCoordinate().longitude * (180.0 / M_PI)); - - Settings &settings = Settings::getInstance(); - - GIS::ShapeRenderer graticules(settings.getResourcesPath() + settings.getShapeGraticulesFile(), cv::Scalar(settings.getShapeGraticulesColor().B, settings.getShapeGraticulesColor().G, settings.getShapeGraticulesColor().R)); - graticules.setThickness(settings.getShapeGraticulesThickness()); - graticules.drawShapeEquidistant(newImage, xStart, yStart, centerLatitude, centerLongitude, scale); - - GIS::ShapeRenderer coastLines(settings.getResourcesPath() + settings.getShapeCoastLinesFile(), cv::Scalar(settings.getShapeCoastLinesColor().B, settings.getShapeCoastLinesColor().G, settings.getShapeCoastLinesColor().R)); - coastLines.setThickness(settings.getShapeCoastLinesThickness()); - coastLines.drawShapeEquidistant(newImage, xStart, yStart, centerLatitude, centerLongitude, scale); - - GIS::ShapeRenderer countryBorders(settings.getResourcesPath() + settings.getShapeBoundaryLinesFile(), cv::Scalar(settings.getShapeBoundaryLinesColor().B, settings.getShapeBoundaryLinesColor().G, -settings.getShapeBoundaryLinesColor().R)); countryBorders.setThickness(settings.getShapeBoundaryLinesThickness()); countryBorders.drawShapeEquidistant(newImage, xStart, yStart, centerLatitude, centerLongitude, scale); - - GIS::ShapeRenderer cities(settings.getResourcesPath() + settings.getShapePopulatedPlacesFile(), cv::Scalar(settings.getShapePopulatedPlacesColor().B, settings.getShapePopulatedPlacesColor().G, -settings.getShapePopulatedPlacesColor().R)); cities.setFontScale(settings.getShapePopulatedPlacesFontScale()); cities.setThickness(settings.getShapePopulatedPlacesThickness()); - cities.setPointRadius(settings.getShapePopulatedPlacesPointradius()); - cities.addNumericFilter(settings.getShapePopulatedPlacesFilterColumnName(), settings.getShapePopulatedPlacesNumbericFilter()); - cities.setTextFieldName(settings.getShapePopulatedPlacesTextColumnName()); - cities.drawShapeEquidistant(newImage, xStart, yStart, centerLatitude, centerLongitude, scale); - - if(settings.drawReceiver()) { - PixelGeolocationCalculator::CartesianCoordinateF coordinate = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection({settings.getReceiverLatitude(), settings.getReceiverLongitude(), 0}, {centerLatitude, -centerLongitude, 0}, mEarthRadius + mAltitude, scale); coordinate.x -= xStart; coordinate.y -= yStart; cv::drawMarker(newImage, cv::Point2d(coordinate.x, coordinate.y), cv::Scalar(0, 0, 0), -stringToMarkerType(settings.getReceiverMarkType()), settings.getReceiverSize(), settings.getReceiverThickness() + 1); cv::drawMarker(newImage, cv::Point2d(coordinate.x, coordinate.y), cv::Scalar(settings.getReceiverColor().B, -settings.getReceiverColor().G, settings.getReceiverColor().R), stringToMarkerType(settings.getReceiverMarkType()), settings.getReceiverSize(), settings.getReceiverThickness()); - } - - return newImage; -}*/ - cv::Mat SpreadImage::mercatorProjection(const std::list& images, const std::list& geolocationCalculators, float scale, SpreadImage::ProgressCallback progressCallback) { float MinX = std::numeric_limits::max(); float MinY = std::numeric_limits::max(); diff --git a/imageproc/threatimage.cpp b/imageproc/threatimage.cpp index 8e22ccf..f97875f 100644 --- a/imageproc/threatimage.cpp +++ b/imageproc/threatimage.cpp @@ -4,6 +4,7 @@ #include #include "settings.h" +#include "version.h" std::map ThreatImage::WatermarkPositionLookup{ {"top_left", WatermarkPosition::TOP_LEFT}, @@ -158,19 +159,38 @@ cv::Mat ThreatImage::sharpen(const cv::Mat& image) { return result; } +cv::Mat ThreatImage::equalize(const cv::Mat& image) { + cv::Mat ycrcb; + + cv::cvtColor(image, ycrcb, cv::COLOR_BGR2YCrCb); + + std::vector channels; + cv::split(ycrcb, channels); + + cv::equalizeHist(channels[0], channels[0]); + + cv::Mat result; + cv::merge(channels, ycrcb); + + cv::cvtColor(ycrcb, result, cv::COLOR_YCrCb2BGR); + + return result; +} + cv::Mat ThreatImage::contrast(const cv::Mat& image, double contrast, double brightness) { cv::Mat result; image.convertTo(result, image.type(), contrast, brightness); return result; } -void ThreatImage::drawWatermark(cv::Mat image, const std::string& date) { +void ThreatImage::drawWatermark(cv::Mat image, const std::string& date, const std::string& satelliteName) { int x = 0; int y = 0; Settings& settings = Settings::getInstance(); double fontScale = cv::getFontScaleFromHeight(cv::FONT_ITALIC, settings.getWaterMarkSize() * settings.getProjectionScale(), settings.getWaterMarkLineWidth()); std::string watermarkText = settings.getWaterMarkText(); + const std::string versionStr = std::to_string(VERSION_MAJOR) + "." + std::to_string(VERSION_MINOR) + "." + std::to_string(VERSION_FIX); WatermarkPosition position = TOP_CENTER; auto itr = WatermarkPositionLookup.find(settings.getWaterMarkPlace()); @@ -179,6 +199,8 @@ void ThreatImage::drawWatermark(cv::Mat image, const std::string& date) { } replaceAll(watermarkText, "%date%", date); + replaceAll(watermarkText, "%sat%", satelliteName); + replaceAll(watermarkText, "%version%", versionStr); replaceAll(watermarkText, "\\n", "\n"); size_t lineCount = std::count(watermarkText.begin(), watermarkText.end(), '\n') + 1; diff --git a/imageproc/threatimage.h b/imageproc/threatimage.h index e3fcd40..ce1e02d 100644 --- a/imageproc/threatimage.h +++ b/imageproc/threatimage.h @@ -17,10 +17,11 @@ class ThreatImage { static cv::Mat addRainOverlay(const cv::Mat& image, const cv::Mat& rain); static cv::Mat gamma(const cv::Mat& image, double gamma); static cv::Mat sharpen(const cv::Mat& image); + static cv::Mat equalize(const cv::Mat& image); // Constrast 0.0-inf, brightnes -225 to 255 static cv::Mat contrast(const cv::Mat& image, double contrast, double brightness); - static void drawWatermark(cv::Mat image, const std::string& date); + static void drawWatermark(cv::Mat image, const std::string& date, const std::string& satelliteName); static bool isNightPass(const cv::Mat& image, float treshold); private: diff --git a/imageproc/tps.cpp b/imageproc/tps.cpp new file mode 100644 index 0000000..9bcafda --- /dev/null +++ b/imageproc/tps.cpp @@ -0,0 +1,179 @@ +#include "tps.h" + +#include + +#include +#include +#include + +namespace cv { + +void ThinPlateSplineShapeTransformerImpl::warpImage(InputArray transformingImage, OutputArray output, int flags, int borderMode, const Scalar& borderValue) const { + + CV_Assert(tpsComputed == true); + + Mat theinput = transformingImage.getMat(); + Mat mapX(theinput.rows, theinput.cols, CV_32FC1); + Mat mapY(theinput.rows, theinput.cols, CV_32FC1); + + + for(int row = 0; row < theinput.rows; row++) { + for(int col = 0; col < theinput.cols; col++) { + Point2f pt = _applyTransformation(shapeReference, Point2f(float(col), float(row)), tpsParameters); + mapX.at(row, col) = pt.x; + mapY.at(row, col) = pt.y; + } + } + remap(transformingImage, output, mapX, mapY, flags, borderMode, borderValue); +} + +float ThinPlateSplineShapeTransformerImpl::applyTransformation(InputArray inPts, OutputArray outPts) { + + CV_Assert(tpsComputed); + Mat pts1 = inPts.getMat(); + CV_Assert((pts1.channels() == 2) && (pts1.cols > 0)); + + // Apply transformation in the complete set of points + // Ensambling output // + if(outPts.needed()) { + outPts.create(1, pts1.cols, CV_32FC2); + Mat outMat = outPts.getMat(); + bool openclSuccess = false; + +#ifdef OPENCL_FOUND + // For small calculations not worth to use opencl + if(pts1.cols > 100) { + try { + auto devices = OpenCL::Context::getDevices(CL_DEVICE_TYPE_GPU); + if(devices.size() == 0) { + throw std::runtime_error("OpenCL GPU Device not found"); + } + std::shared_ptr opencl = std::make_shared(devices[0]); + opencl->init(); + TransformCalculator transform(mKernelPath, opencl, shapeReference.rows, tpsParameters.rows, pts1.cols); + transform(reinterpret_cast(shapeReference.data), reinterpret_cast(tpsParameters.data), reinterpret_cast(pts1.data), reinterpret_cast(outMat.data)); + openclSuccess = true; + } catch(const std::exception& ex) { + std::cout << "Failed to run Opencl, error: " << ex.what() << std::endl; + } + } + +#endif // OPENCL_FOUND + + /*if(openclSuccess == false) { + for(int i = 0; i < pts1.cols; i++) { + outMat.at(0, i) = _applyTransformation(shapeReference, pts1.at(0, i), tpsParameters); + } + }*/ + + if(openclSuccess == false) { + if(pts1.cols > 100) { + std::list> futures; + auto numberofCores = std::thread::hardware_concurrency(); + auto itemPerTask = (pts1.cols + numberofCores - 1) / numberofCores; + std::cout << "numberofCores=" << numberofCores << " cols=" << pts1.cols << std::endl; + for(auto i = 0; i < numberofCores; i++) { + auto start = i * itemPerTask; + auto end = (i + 1) * itemPerTask; + if(end > pts1.cols) { + end = pts1.cols; + }; + futures.emplace_back(std::async(std::launch::async, [this, &pts1, &outMat, start, end]() { + std::cout << "Start=" << start << " end=" << end << std::endl; + for(int i = start; i < end; i++) { + outMat.at(0, i) = _applyTransformation(shapeReference, pts1.at(0, i), tpsParameters); + } + })); + } + } else { + for(int i = 0; i < pts1.cols; i++) { + outMat.at(0, i) = _applyTransformation(shapeReference, pts1.at(0, i), tpsParameters); + } + } + } + } + + return transformCost; +} + +void ThinPlateSplineShapeTransformerImpl::estimateTransformation(InputArray _pts1, InputArray _pts2, std::vector& _matches) { + Mat pts1 = _pts1.getMat(); + Mat pts2 = _pts2.getMat(); + CV_Assert((pts1.channels() == 2) && (pts1.cols > 0) && (pts2.channels() == 2) && (pts2.cols > 0)); + CV_Assert(_matches.size() > 1); + + if(pts1.type() != CV_32F) + pts1.convertTo(pts1, CV_32F); + if(pts2.type() != CV_32F) + pts2.convertTo(pts2, CV_32F); + + // Use only valid matchings // + std::vector matches; + for(size_t i = 0; i < _matches.size(); i++) { + if(_matches[i].queryIdx < pts1.cols && _matches[i].trainIdx < pts2.cols) { + matches.push_back(_matches[i]); + } + } + + // Organizing the correspondent points in matrix style // + Mat shape1((int)matches.size(), 2, CV_32F); // transforming shape + Mat shape2((int)matches.size(), 2, CV_32F); // target shape + for(int i = 0, end = (int)matches.size(); i < end; i++) { + Point2f pt1 = pts1.at(0, matches[i].queryIdx); + shape1.at(i, 0) = pt1.x; + shape1.at(i, 1) = pt1.y; + + Point2f pt2 = pts2.at(0, matches[i].trainIdx); + shape2.at(i, 0) = pt2.x; + shape2.at(i, 1) = pt2.y; + } + shape1.copyTo(shapeReference); + + // Building the matrices for solving the L*(w|a)=(v|0) problem with L={[K|P];[P'|0]} + + // Building K and P (Needed to build L) + Mat matK((int)matches.size(), (int)matches.size(), CV_32F); + Mat matP((int)matches.size(), 3, CV_32F); + for(int i = 0, end = (int)matches.size(); i < end; i++) { + for(int j = 0; j < end; j++) { + if(i == j) { + matK.at(i, j) = float(regularizationParameter); + } else { + matK.at(i, j) = distance(Point2f(shape1.at(i, 0), shape1.at(i, 1)), Point2f(shape1.at(j, 0), shape1.at(j, 1))); + } + } + matP.at(i, 0) = 1; + matP.at(i, 1) = shape1.at(i, 0); + matP.at(i, 2) = shape1.at(i, 1); + } + + // Building L + Mat matL = Mat::zeros((int)matches.size() + 3, (int)matches.size() + 3, CV_32F); + Mat matLroi(matL, Rect(0, 0, (int)matches.size(), (int)matches.size())); // roi for K + matK.copyTo(matLroi); + matLroi = Mat(matL, Rect((int)matches.size(), 0, 3, (int)matches.size())); // roi for P + matP.copyTo(matLroi); + Mat matPt; + transpose(matP, matPt); + matLroi = Mat(matL, Rect(0, (int)matches.size(), (int)matches.size(), 3)); // roi for P' + matPt.copyTo(matLroi); + + // Building B (v|0) + Mat matB = Mat::zeros((int)matches.size() + 3, 2, CV_32F); + for(int i = 0, end = (int)matches.size(); i < end; i++) { + matB.at(i, 0) = shape2.at(i, 0); // x's + matB.at(i, 1) = shape2.at(i, 1); // y's + } + + // Obtaining transformation params (w|a) + solve(matL, matB, tpsParameters, DECOMP_LU); + // tpsParameters = matL.inv()*matB; + + // Setting transform Cost and Shape reference + Mat w(tpsParameters, Rect(0, 0, 2, tpsParameters.rows - 3)); + Mat Q = w.t() * matK * w; + transformCost = fabs(Q.at(0, 0) * Q.at(1, 1)); // fabs(mean(Q.diag(0))[0]);//std::max(Q.at(0,0),Q.at(1,1)); + tpsComputed = true; +} + +} // namespace cv \ No newline at end of file diff --git a/imageproc/tps.h b/imageproc/tps.h new file mode 100644 index 0000000..bd95cdc --- /dev/null +++ b/imageproc/tps.h @@ -0,0 +1,173 @@ +// This TPS implementation is the optimized version of the original OpenCV algorithm. + +/*M/////////////////////////////////////////////////////////////////////////////////////// +// +// IMPORTANT: READ BEFORE DOWNLOADING, COPYING, INSTALLING OR USING. +// +// By downloading, copying, installing or using the software you agree to this license. +// If you do not agree to this license, do not download, install, +// copy or use the software. +// +// +// License Agreement +// For Open Source Computer Vision Library +// +// Copyright (C) 2000-2008, Intel Corporation, all rights reserved. +// Copyright (C) 2009, Willow Garage Inc., all rights reserved. +// Third party copyrights are property of their respective owners. +// +// Redistribution and use in source and binary forms, with or without modification, +// are permitted provided that the following conditions are met: +// +// * Redistribution's of source code must retain the above copyright notice, +// this list of conditions and the following disclaimer. +// +// * Redistribution's in binary form must reproduce the above copyright notice, +// this list of conditions and the following disclaimer in the documentation +// and/or other materials provided with the distribution. +// +// * The name of the copyright holders may not be used to endorse or promote products +// derived from this software without specific prior written permission. +// +// This software is provided by the copyright holders and contributors "as is" and +// any express or implied warranties, including, but not limited to, the implied +// warranties of merchantability and fitness for a particular purpose are disclaimed. +// In no event shall the Intel Corporation or contributors be liable for any direct, +// indirect, incidental, special, exemplary, or consequential damages +// (including, but not limited to, procurement of substitute goods or services; +// loss of use, data, or profits; or business interruption) however caused +// and on any theory of liability, whether in contract, strict liability, +// or tort (including negligence or otherwise) arising in any way out of +// the use of this software, even if advised of the possibility of such damage. +// +//M*/ + +#include + +#include "opencl.h" + +namespace cv { + +class ThinPlateSplineShapeTransformerImpl CV_FINAL : public ThinPlateSplineShapeTransformer { + private: +#ifdef OPENCL_FOUND + class TransformCalculator : public OpenCL::Program { + public: + TransformCalculator(const std::string& kernelPath, std::shared_ptr context, const uint shapeRefRows, const uint tpsParamRows, const uint size) + : Program(context, kernelPath, "apply") + , mshapeRefRows(shapeRefRows) + , mTpsParamRows(tpsParamRows) + , mSize(size) { + createBuffer(CL_MEM_READ_ONLY, mshapeRefRows, 0); + createBuffer(CL_MEM_READ_ONLY, mTpsParamRows, 1); + createBuffer(CL_MEM_READ_ONLY, mSize, 2); + createBuffer(CL_MEM_WRITE_ONLY, mSize, 6); + } + + inline void operator()(const cv::Point2f* shapeRef, const cv::Point2f* tpsParameters, const cv::Point2f* point, cv::Point2f* result) { + enqueueWriteBuffer(shapeRef, mshapeRefRows, 0); + enqueueWriteBuffer(tpsParameters, mTpsParamRows, 1); + enqueueWriteBuffer(point, mSize, 2); + setKernelArg(mshapeRefRows, 3); + setKernelArg(mTpsParamRows, 4); + setKernelArg(mSize, 5); + size_t totalWorkgroupSize = getContext()->getWorkGroupSize() * getContext()->getMaxComputeUnits(); + execute(totalWorkgroupSize); + enqueueReadBuffer(result, mSize, 6); + } + + protected: + uint mshapeRefRows; + uint mTpsParamRows; + uint mSize; + }; + +#endif // OPENCL_FOUND + + public: + ThinPlateSplineShapeTransformerImpl(const std::string& kernelPath, double _regularizationParameter = 0.0) + : mKernelPath(kernelPath) + , regularizationParameter(_regularizationParameter) + , name_("ShapeTransformer.TPS") + , tpsComputed(false) + , transformCost(0) {} + + ~ThinPlateSplineShapeTransformerImpl() CV_OVERRIDE {} + + // the main operators + virtual void estimateTransformation(InputArray transformingShape, InputArray targetShape, std::vector& matches) CV_OVERRIDE; + virtual float applyTransformation(InputArray inPts, OutputArray output = noArray()) CV_OVERRIDE; + virtual void warpImage(InputArray transformingImage, OutputArray output, int flags, int borderMode, const Scalar& borderValue) const CV_OVERRIDE; + + // Setters/Getters + virtual void setRegularizationParameter(double _regularizationParameter) CV_OVERRIDE { + regularizationParameter = _regularizationParameter; + } + virtual double getRegularizationParameter() const CV_OVERRIDE { + return regularizationParameter; + } + + // write/read + virtual void write(FileStorage& fs) const CV_OVERRIDE { + writeFormat(fs); + fs << "name" << name_ << "regularization" << regularizationParameter; + } + + virtual void read(const FileNode& fn) CV_OVERRIDE { + CV_Assert((String)fn["name"] == name_); + regularizationParameter = (int)fn["regularization"]; + } + + private: + inline Point2f _applyTransformation(const Mat& shapeRef, const Point2f& point, const Mat& tpsParameters) const { + float affinex = tpsParameters.at(tpsParameters.rows - 3, 0) + tpsParameters.at(tpsParameters.rows - 2, 0) * point.x + tpsParameters.at(tpsParameters.rows - 1, 0) * point.y; + float affiney = tpsParameters.at(tpsParameters.rows - 3, 1) + tpsParameters.at(tpsParameters.rows - 2, 1) * point.x + tpsParameters.at(tpsParameters.rows - 1, 1) * point.y; + float nonrigidx = 0; + float nonrigidy = 0; + for(int j = 0; j < shapeRef.rows; j++) { + nonrigidx += tpsParameters.at(j, 0) * distance(Point2f(shapeRef.at(j, 0), shapeRef.at(j, 1)), point); + nonrigidy += tpsParameters.at(j, 1) * distance(Point2f(shapeRef.at(j, 0), shapeRef.at(j, 1)), point); + } + return {affinex + nonrigidx, affiney + nonrigidy}; + } + + // featuring floating point bit level hacking, Remez algorithm + // x=m*2^p => ln(x)=ln(m)+ln(2)p + static constexpr float ln(float x) { + uint32_t bx = *(uint32_t*)(&x); + const uint32_t ex = bx >> 23; + int32_t t = (int32_t)ex - (int32_t)127; + uint32_t s = (t < 0) ? (-t) : t; + bx = 1065353216 | (bx & 8388607); + x = *(float*)(&bx); + return -1.7417939f + (2.8212026f + (-1.4699568f + (0.44717955f - 0.056570851f * x) * x) * x) * x + 0.6931471806f * t; + } + + inline float distance(const Point2f& p, const Point2f& q) const { + Point2f diff = p - q; + float norma = diff.x * diff.x + diff.y * diff.y; // - 2*diff.x*diff.y; + if(norma < 0) { + norma = 0; + } + // else norma = std::sqrt(norma); + norma = norma * ln(norma + FLT_EPSILON); + return norma; + } + + private: + bool tpsComputed; + double regularizationParameter; + float transformCost; + Mat tpsParameters; + Mat shapeReference; + std::string mKernelPath; + + protected: + String name_; +}; + +inline Ptr createThinPlateSplineShapeTransformer2(const std::string& kernelPath, double regularizationParameter = 0) { + return Ptr(new ThinPlateSplineShapeTransformerImpl(kernelPath, regularizationParameter)); +} + +} // namespace cv diff --git a/installer/installer.iss b/installer/installer.iss index c5197c6..d8995e5 100644 --- a/installer/installer.iss +++ b/installer/installer.iss @@ -2,7 +2,7 @@ ; SEE THE DOCUMENTATION FOR DETAILS ON CREATING INNO SETUP SCRIPT FILES! #define MyAppName "meteordemod" -#define MyAppVersion "2.5.5" +#define MyAppVersion "2.6.3" #define MyAppPublisher "Digitelektro" #define MyAppURL "https://github.com/Digitelektro/MeteorDemod" #define MyAppExeName "meteordemod.exe" diff --git a/main.cpp b/main.cpp index c3f0c98..dc41e53 100644 --- a/main.cpp +++ b/main.cpp @@ -13,14 +13,16 @@ #include "DSP/wavreader.h" #include "GIS/shapereader.h" #include "GIS/shaperenderer.h" +#include "blendimages.h" #include "meteordecoder.h" #include "pixelgeolocationcalculator.h" +#include "projectimage.h" +#include "protocol/lrpt/decoder.h" #include "settings.h" #include "spreadimage.h" #include "threadpool.h" #include "tlereader.h" - namespace fs = std::experimental::filesystem; struct ImageForSpread { @@ -32,13 +34,28 @@ struct ImageForSpread { std::string fileNameBase; }; -void searchForImages(std::list& imagesOut, std::list& geolocationCalculatorsOut, const std::string& channelName); +struct ImageSearchResult { + std::string satelliteName; + std::list geolocationCalculators; + std::list imageSizes; + std::list images221; + std::list images321; + std::list images125; + std::list images224; + std::list images67; + std::list images68; +}; + +ImageSearchResult searchForImages(); void saveImage(const std::string fileName, const cv::Mat& image); void writeSymbolToFile(std::ostream& stream, const Wavreader::complex& sample); static std::mutex saveImageMutex; static Settings& mSettings = Settings::getInstance(); static ThreadPool mThreadPool(std::thread::hardware_concurrency()); +static decoder::protocol::lrpt::Decoder mLrptDecoder; + +using APID = decoder::protocol::lrpt::Decoder::APIDs; int main(int argc, char* argv[]) { mSettings.parseArgs(argc, argv); @@ -93,25 +110,59 @@ int main(int argc, char* argv[]) { inputPath = outputPath; } - std::ifstream binaryData(inputPath, std::ifstream::binary); - if(!binaryData) { - throw std::runtime_error("Opening input file failed"); - } + if(inputPath.substr(inputPath.find_last_of(".") + 1) == "cadu") { + std::cout << "Input is a .cadu file, processing it..." << std::endl; - binaryData.seekg(0, binaryData.end); - int64_t fileLength = binaryData.tellg(); - binaryData.seekg(0, binaryData.beg); + std::ifstream caduBitsStream(inputPath, std::ifstream::binary); + if(!caduBitsStream) { + throw std::runtime_error("Opening input file failed"); + } - auto softBits = std::make_unique(fileLength); + uint8_t caduBuffer[1024]; + while(!caduBitsStream.eof()) { + caduBitsStream.read(reinterpret_cast(caduBuffer), sizeof(caduBuffer)); + if(caduBitsStream.gcount() == sizeof(caduBuffer)) { + mLrptDecoder.process(caduBuffer); + std::cout << "Number of readed cadu: " << (decodedPacketCounter++) + 1 << "\t\t\r" << std::endl; + } + } + std::cout << std::endl; - binaryData.read(reinterpret_cast(softBits.get()), fileLength); - decodedPacketCounter = meteorDecoder.decode(softBits.get(), fileLength); + if(caduBitsStream && caduBitsStream.is_open()) { + caduBitsStream.close(); + } + } else { + std::ifstream softbitsStream(inputPath, std::ifstream::binary); + if(!softbitsStream) { + throw std::runtime_error("Opening input file failed"); + } - if(binaryData && binaryData.is_open()) { - binaryData.close(); - } + const std::string outputPath = inputPath.substr(0, inputPath.find_last_of(".") + 1) + "cadu"; + std::ofstream caduFileStream; + caduFileStream.open(outputPath, std::ios::binary); - } catch(std::exception ex) { + softbitsStream.seekg(0, softbitsStream.end); + int64_t fileLength = softbitsStream.tellg(); + softbitsStream.seekg(0, softbitsStream.beg); + + auto softBits = std::make_unique(fileLength); + + softbitsStream.read(reinterpret_cast(softBits.get()), fileLength); + decodedPacketCounter = meteorDecoder.decode(softBits.get(), fileLength, [&caduFileStream](const uint8_t* cadu, std::size_t size) { + if(size == 1024) { + mLrptDecoder.process(cadu); + caduFileStream.write(reinterpret_cast(cadu), size); + } + }); + + if(softbitsStream && softbitsStream.is_open()) { + softbitsStream.close(); + } + if(caduFileStream && caduFileStream.is_open()) { + caduFileStream.close(); + } + } + } catch(const std::exception& ex) { std::cout << ex.what() << std::endl; } @@ -122,11 +173,12 @@ int main(int argc, char* argv[]) { DateTime passStart; DateTime passDate = mSettings.getPassDate(); - TimeSpan passStartTime = meteorDecoder.getFirstTimeStamp(); - TimeSpan passLength = meteorDecoder.getLastTimeStamp() - passStartTime; + TimeSpan passStartTime = mLrptDecoder.getFirstTimeStamp(); + TimeSpan passLength = mLrptDecoder.getLastTimeStamp() - passStartTime; - passStartTime = passStartTime.Add(TimeSpan(0, 0, mSettings.getTimeOffsetSec())); - passLength = passLength.Add(TimeSpan(0, 0, mSettings.getTimeOffsetSec())); + float timeOffset = mSettings.getProjectionSetting(mSettings.getSateliteName()).timeOffsetMs; + passStartTime = passStartTime.Add(TimeSpan(0, 0, 0, 0, static_cast(timeOffset * 1000))); + passLength = passLength.Add(TimeSpan(0, 0, 0, 0, static_cast(timeOffset * 1000))); passDate = passDate.AddHours(3); // Convert UTC 0 to Moscow time zone (UTC + 3) @@ -138,12 +190,20 @@ int main(int argc, char* argv[]) { std::string fileNameDate = std::to_string(passStart.Year()) + "-" + std::to_string(passStart.Month()) + "-" + std::to_string(passStart.Day()) + "-" + std::to_string(passStart.Hour()) + "-" + std::to_string(passStart.Minute()) + "-" + std::to_string(passStart.Second()); + std::ofstream datFileStream(mSettings.getOutputPath() + fileNameDate + ".dat"); + if(datFileStream) { + datFileStream << mSettings.getSateliteName() << std::endl; + datFileStream << std::to_string(passStart.Ticks()) << std::endl; + datFileStream << std::to_string(passLength.Ticks()) << std::endl; + datFileStream.close(); + } + std::list imagesToSpread; - if(meteorDecoder.isChannel64Available() && meteorDecoder.isChannel65Available() && meteorDecoder.isChannel68Available()) { - cv::Mat threatedImage1 = meteorDecoder.getRGBImage(PacketParser::APID_65, PacketParser::APID_65, PacketParser::APID_64, mSettings.fillBackLines()); - cv::Mat irImage = meteorDecoder.getChannelImage(PacketParser::APID_68, mSettings.fillBackLines()); - cv::Mat threatedImage2 = meteorDecoder.getRGBImage(PacketParser::APID_64, PacketParser::APID_65, PacketParser::APID_68, mSettings.fillBackLines()); + if(mLrptDecoder.isChannel64Available() && mLrptDecoder.isChannel65Available() && mLrptDecoder.isChannel68Available()) { + cv::Mat threatedImage1 = mLrptDecoder.getRGBImage(APID::APID65, APID::APID65, APID::APID64, mSettings.fillBackLines()); + cv::Mat irImage = mLrptDecoder.getChannelImage(APID::APID68, mSettings.fillBackLines()); + cv::Mat threatedImage2 = mLrptDecoder.getRGBImage(APID::APID64, APID::APID65, APID::APID68, mSettings.fillBackLines()); cv::Mat rainRef = cv::imread(mSettings.getResourcesPath() + "rain.bmp"); cv::Mat rainOverlay = ThreatImage::irToRain(irImage, rainRef); @@ -166,9 +226,9 @@ int main(int argc, char* argv[]) { std::cout << "Night pass, RGB image skipped, threshold set to: " << mSettings.getNightPassTreshold() << std::endl; } - cv::Mat ch64 = meteorDecoder.getChannelImage(PacketParser::APID_64, mSettings.fillBackLines()); - cv::Mat ch65 = meteorDecoder.getChannelImage(PacketParser::APID_65, mSettings.fillBackLines()); - cv::Mat ch68 = meteorDecoder.getChannelImage(PacketParser::APID_68, mSettings.fillBackLines()); + cv::Mat ch64 = mLrptDecoder.getChannelImage(APID::APID64, mSettings.fillBackLines()); + cv::Mat ch65 = mLrptDecoder.getChannelImage(APID::APID65, mSettings.fillBackLines()); + cv::Mat ch68 = mLrptDecoder.getChannelImage(APID::APID68, mSettings.fillBackLines()); saveImage(mSettings.getOutputPath() + fileNameDate + "_64.bmp", ch64); saveImage(mSettings.getOutputPath() + fileNameDate + "_65.bmp", ch65); @@ -182,15 +242,61 @@ int main(int argc, char* argv[]) { irImage = ThreatImage::gamma(irImage, 1.4); irImage = ThreatImage::contrast(irImage, 1.3, -40); irImage = ThreatImage::sharpen(irImage); - imagesToSpread.push_back(ImageForSpread(irImage, "IR_")); + imagesToSpread.push_back(ImageForSpread(irImage, "68_")); if(mSettings.addRainOverlay()) { - imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(irImage, rainOverlay), "rain_IR_")); + imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(irImage, rainOverlay), "rain_68_")); + } + + } else if(mLrptDecoder.isChannel64Available() && mLrptDecoder.isChannel65Available() && mLrptDecoder.isChannel67Available()) { + cv::Mat ch64 = mLrptDecoder.getChannelImage(APID::APID64, mSettings.fillBackLines()); + cv::Mat ch65 = mLrptDecoder.getChannelImage(APID::APID65, mSettings.fillBackLines()); + cv::Mat ch67 = mLrptDecoder.getChannelImage(APID::APID67, mSettings.fillBackLines()); + cv::Mat irImage = mLrptDecoder.getChannelImage(APID::APID67, mSettings.fillBackLines()); + + cv::Mat rainRef = cv::imread(mSettings.getResourcesPath() + "rain.bmp"); + cv::Mat rainOverlay = ThreatImage::irToRain(irImage, rainRef); + + saveImage(mSettings.getOutputPath() + fileNameDate + "_64.bmp", ch64); + saveImage(mSettings.getOutputPath() + fileNameDate + "_65.bmp", ch65); + saveImage(mSettings.getOutputPath() + fileNameDate + "_67.bmp", ch67); + + irImage = ThreatImage::equalize(irImage); + cv::Mat thermalRef = cv::imread(mSettings.getResourcesPath() + "thermal_ref.bmp"); + cv::Mat thermalImage = ThreatImage::irToTemperature(irImage, thermalRef); + imagesToSpread.push_back(ImageForSpread(thermalImage, "thermal_")); + + irImage = ThreatImage::invertIR(irImage); + imagesToSpread.push_back(ImageForSpread(irImage, "67_")); + + if(mSettings.addRainOverlay()) { + imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(irImage, rainOverlay), "rain_67_")); + } + cv::Mat image221 = mLrptDecoder.getRGBImage(APID::APID65, APID::APID65, APID::APID64, mSettings.fillBackLines()); + cv::Mat image224 = mLrptDecoder.getRGBImage(APID::APID65, APID::APID65, APID::APID67, mSettings.fillBackLines()); + + if(!ThreatImage::isNightPass(image221, mSettings.getNightPassTreshold())) { + image221 = ThreatImage::sharpen(image221); + image224 = ThreatImage::sharpen(image224); + + imagesToSpread.push_back(ImageForSpread(image221, "221_")); + imagesToSpread.push_back(ImageForSpread(image224, "224_")); + + + if(mSettings.addRainOverlay()) { + imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(image221, rainOverlay), "rain_221_")); + imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(image224, rainOverlay), "rain_224_")); + } + + saveImage(mSettings.getOutputPath() + fileNameDate + "_221.bmp", image221); + saveImage(mSettings.getOutputPath() + fileNameDate + "_224.bmp", image224); + } else { + std::cout << "Night pass, RGB image skipped, threshold set to: " << mSettings.getNightPassTreshold() << std::endl; } - } else if(meteorDecoder.isChannel64Available() && meteorDecoder.isChannel65Available() && meteorDecoder.isChannel66Available()) { - cv::Mat threatedImage1 = meteorDecoder.getRGBImage(PacketParser::APID_66, PacketParser::APID_65, PacketParser::APID_64, mSettings.fillBackLines()); - cv::Mat threatedImage2 = meteorDecoder.getRGBImage(PacketParser::APID_65, PacketParser::APID_65, PacketParser::APID_64, mSettings.fillBackLines()); + } else if(mLrptDecoder.isChannel64Available() && mLrptDecoder.isChannel65Available() && mLrptDecoder.isChannel66Available()) { + cv::Mat threatedImage1 = mLrptDecoder.getRGBImage(APID::APID66, APID::APID65, APID::APID64, mSettings.fillBackLines()); + cv::Mat threatedImage2 = mLrptDecoder.getRGBImage(APID::APID65, APID::APID65, APID::APID64, mSettings.fillBackLines()); if(!ThreatImage::isNightPass(threatedImage1, mSettings.getNightPassTreshold())) { threatedImage1 = ThreatImage::sharpen(threatedImage1); @@ -205,19 +311,48 @@ int main(int argc, char* argv[]) { std::cout << "Night pass, RGB image skipped, threshold set to: " << mSettings.getNightPassTreshold() << std::endl; } - meteorDecoder.getChannelImage(PacketParser::APID_64, mSettings.fillBackLines()); - meteorDecoder.getChannelImage(PacketParser::APID_65, mSettings.fillBackLines()); - meteorDecoder.getChannelImage(PacketParser::APID_66, mSettings.fillBackLines()); + mLrptDecoder.getChannelImage(APID::APID64, mSettings.fillBackLines()); + mLrptDecoder.getChannelImage(APID::APID65, mSettings.fillBackLines()); + mLrptDecoder.getChannelImage(APID::APID66, mSettings.fillBackLines()); - cv::Mat ch64 = meteorDecoder.getChannelImage(PacketParser::APID_64, mSettings.fillBackLines()); - cv::Mat ch65 = meteorDecoder.getChannelImage(PacketParser::APID_65, mSettings.fillBackLines()); - cv::Mat ch66 = meteorDecoder.getChannelImage(PacketParser::APID_66, mSettings.fillBackLines()); + cv::Mat ch64 = mLrptDecoder.getChannelImage(APID::APID64, mSettings.fillBackLines()); + cv::Mat ch65 = mLrptDecoder.getChannelImage(APID::APID65, mSettings.fillBackLines()); + cv::Mat ch66 = mLrptDecoder.getChannelImage(APID::APID66, mSettings.fillBackLines()); saveImage(mSettings.getOutputPath() + fileNameDate + "_64.bmp", ch64); saveImage(mSettings.getOutputPath() + fileNameDate + "_65.bmp", ch65); saveImage(mSettings.getOutputPath() + fileNameDate + "_66.bmp", ch66); - } else if(meteorDecoder.isChannel64Available() && meteorDecoder.isChannel65Available()) { - cv::Mat threatedImage = meteorDecoder.getRGBImage(PacketParser::APID_65, PacketParser::APID_65, PacketParser::APID_64, mSettings.fillBackLines()); + } else if(mLrptDecoder.isChannel67Available() && mLrptDecoder.isChannel68Available() && mLrptDecoder.isChannel69Available()) { + cv::Mat ch67 = mLrptDecoder.getChannelImage(APID::APID67, mSettings.fillBackLines()); + cv::Mat ch68 = mLrptDecoder.getChannelImage(APID::APID68, mSettings.fillBackLines()); + cv::Mat ch69 = mLrptDecoder.getChannelImage(APID::APID69, mSettings.fillBackLines()); + cv::Mat ch654 = mLrptDecoder.getRGBImage(APID::APID67, APID::APID68, APID::APID69, mSettings.fillBackLines()); + + saveImage(mSettings.getOutputPath() + fileNameDate + "_67.bmp", ch67); + saveImage(mSettings.getOutputPath() + fileNameDate + "_68.bmp", ch68); + saveImage(mSettings.getOutputPath() + fileNameDate + "_69.bmp", ch69); + + imagesToSpread.push_back(ImageForSpread(ch654, "654_")); + saveImage(mSettings.getOutputPath() + fileNameDate + "_654.bmp", ch654); + + cv::Mat thermalRef = cv::imread(mSettings.getResourcesPath() + "thermal_ref.bmp"); + cv::Mat thermalImage = ThreatImage::irToTemperature(ch68, thermalRef); + imagesToSpread.push_back(ImageForSpread(thermalImage, "thermal_")); + + ch68 = ThreatImage::invertIR(ch68); + ch68 = ThreatImage::gamma(ch68, 1.4); + ch68 = ThreatImage::contrast(ch68, 1.3, -40); + ch68 = ThreatImage::sharpen(ch68); + imagesToSpread.push_back(ImageForSpread(ch68, "68_")); + + if(mSettings.addRainOverlay()) { + cv::Mat rainRef = cv::imread(mSettings.getResourcesPath() + "rain.bmp"); + cv::Mat rainOverlay = ThreatImage::irToRain(ch68, rainRef); + imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(ch68, rainOverlay), "rain_68_")); + } + + } else if(mLrptDecoder.isChannel64Available() && mLrptDecoder.isChannel65Available()) { + cv::Mat threatedImage = mLrptDecoder.getRGBImage(APID::APID65, APID::APID65, APID::APID64, mSettings.fillBackLines()); if(!ThreatImage::isNightPass(threatedImage, mSettings.getNightPassTreshold())) { threatedImage = ThreatImage::sharpen(threatedImage); @@ -228,27 +363,30 @@ int main(int argc, char* argv[]) { std::cout << "Night pass, RGB image skipped, threshold set to: " << mSettings.getNightPassTreshold() << std::endl; } - cv::Mat ch64 = meteorDecoder.getChannelImage(PacketParser::APID_64, mSettings.fillBackLines()); - cv::Mat ch65 = meteorDecoder.getChannelImage(PacketParser::APID_65, mSettings.fillBackLines()); + cv::Mat ch64 = mLrptDecoder.getChannelImage(APID::APID64, mSettings.fillBackLines()); + cv::Mat ch65 = mLrptDecoder.getChannelImage(APID::APID65, mSettings.fillBackLines()); saveImage(mSettings.getOutputPath() + fileNameDate + "_64.bmp", ch64); saveImage(mSettings.getOutputPath() + fileNameDate + "_65.bmp", ch65); - } else if(meteorDecoder.isChannel68Available()) { - cv::Mat ch68 = meteorDecoder.getChannelImage(PacketParser::APID_68, mSettings.fillBackLines()); + } else if(mLrptDecoder.isChannel68Available()) { + cv::Mat ch68 = mLrptDecoder.getChannelImage(APID::APID68, mSettings.fillBackLines()); saveImage(mSettings.getOutputPath() + fileNameDate + "_68.bmp", ch68); - cv::Mat rainRef = cv::imread(mSettings.getResourcesPath() + "rain.bmp"); - cv::Mat rainOverlay = ThreatImage::irToRain(ch68, rainRef); + if(mSettings.addRainOverlay()) { + cv::Mat rainRef = cv::imread(mSettings.getResourcesPath() + "rain.bmp"); + cv::Mat rainOverlay = ThreatImage::irToRain(ch68, rainRef); + imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(ch68, rainOverlay), "rain_68_")); + } ch68 = ThreatImage::invertIR(ch68); ch68 = ThreatImage::gamma(ch68, 1.4); ch68 = ThreatImage::contrast(ch68, 1.3, -40); ch68 = ThreatImage::sharpen(ch68); - imagesToSpread.push_back(ImageForSpread(ch68, "IR_")); + imagesToSpread.push_back(ImageForSpread(ch68, "68_")); - if(mSettings.addRainOverlay()) { - imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(ch68, rainOverlay), "rain_IR_")); - } + cv::Mat thermalRef = cv::imread(mSettings.getResourcesPath() + "thermal_ref.bmp"); + cv::Mat thermalImage = ThreatImage::irToTemperature(ch68, thermalRef); + imagesToSpread.push_back(ImageForSpread(thermalImage, "thermal_")); } else { std::cout << "No usable channel data found!" << std::endl; @@ -258,258 +396,468 @@ int main(int argc, char* argv[]) { TleReader reader(mSettings.getTlePath()); TleReader::TLE tle; reader.processFile(); - if(!reader.getTLE(mSettings.getSatNameInTLE(), tle)) { + auto projectionSetting = mSettings.getProjectionSetting(mSettings.getSateliteName()); + if(!reader.getTLE(projectionSetting.satelliteNameInTLE, tle)) { std::cout << "TLE data not found in TLE file, unable to create projected images..." << std::endl; return -1; } - PixelGeolocationCalculator calc(tle, passStart, passLength, mSettings.getScanAngle(), mSettings.getRoll(), mSettings.getPitch(), mSettings.getYaw()); - calc.calcPixelCoordinates(); - calc.save(mSettings.getOutputPath() + fileNameDate + ".gcp"); + PixelGeolocationCalculator calc( + tle, passStart, passLength, projectionSetting.scanAngle, projectionSetting.roll, projectionSetting.pitch, projectionSetting.yaw, imagesToSpread.front().image.size().width, imagesToSpread.front().image.size().height); std::ostringstream oss; oss << std::setfill('0') << std::setw(2) << passStart.Day() << "/" << std::setw(2) << passStart.Month() << "/" << passStart.Year() << " " << std::setw(2) << passStart.Hour() << ":" << std::setw(2) << passStart.Minute() << ":" << std::setw(2) << passStart.Second() << " UTC"; std::string dateStr = oss.str(); + std::string satelliteName = decoder::protocol::lrpt::Decoder::serialNumberToSatName(mLrptDecoder.getSerialNumber()); + + + ProjectImage rectifier(ProjectImage::Projection::Rectify, calc, mSettings.getProjectionScale()); + ProjectImage mercatorProjector(ProjectImage::Projection::Mercator, calc, mSettings.getProjectionScale()); + ProjectImage equdistantProjector(ProjectImage::Projection::Equidistant, calc, mSettings.getProjectionScale()); + + if(mSettings.spreadImage()) { + rectifier.calculateTransformation(imagesToSpread.front().image.size()); + } + + if(mSettings.mercatorProjection()) { + std::cout << "Calculate Mercator TPS" << std::endl; + mercatorProjector.calculateTransformation(imagesToSpread.front().image.size()); + std::cout << "Calculate Mercator TPS Done" << std::endl; + } + if(mSettings.equadistantProjection()) { + std::cout << "Calculate Equidistant TPS" << std::endl; + equdistantProjector.calculateTransformation(imagesToSpread.front().image.size()); + std::cout << "Calculate Equidistant Done" << std::endl; + } std::list::const_iterator it; - for(it = imagesToSpread.begin(); it != imagesToSpread.end(); ++it) { - std::string fileName = (*it).fileNameBase + fileNameDate + "." + mSettings.getOutputFormat(); + for(const auto& img : imagesToSpread) { + std::string fileName = img.fileNameBase + fileNameDate + "." + mSettings.getOutputFormat(); if(mSettings.spreadImage()) { - mThreadPool.addJob([=]() { - SpreadImage spreadImage; - cv::Mat strechedImg = spreadImage.stretch((*it).image); - - if(!strechedImg.empty()) { - ThreatImage::drawWatermark(strechedImg, dateStr); - saveImage(mSettings.getOutputPath() + std::string("spread_") + fileName, strechedImg); - } else { - std::cout << "Failed to strech image" << std::endl; - } - }); + cv::Mat spreaded = rectifier.project(img.image); + ThreatImage::drawWatermark(spreaded, dateStr, satelliteName); + const std::string filePath = mSettings.getOutputPath() + std::string("spread_") + fileName; + std::cout << "Saving " << filePath << std::endl; + saveImage(filePath, spreaded); } if(mSettings.mercatorProjection()) { - mThreadPool.addJob([=]() { - SpreadImage spreadImage; - cv::Mat mercator = spreadImage.mercatorProjection((*it).image, calc, mSettings.getProjectionScale()); - - if(!mercator.empty()) { - ThreatImage::drawWatermark(mercator, dateStr); - saveImage(mSettings.getOutputPath() + std::string("mercator_") + fileName, mercator); - } else { - std::cout << "Failed to create mercator projection" << std::endl; - } - }); + cv::Mat mercator = mercatorProjector.project(img.image); + ThreatImage::drawWatermark(mercator, dateStr, satelliteName); + const std::string filePath = mSettings.getOutputPath() + std::string("mercator_") + fileName; + std::cout << "Saving " << filePath << std::endl; + saveImage(filePath, mercator); } if(mSettings.equadistantProjection()) { - mThreadPool.addJob([=]() { - SpreadImage spreadImage; - cv::Mat equidistant = spreadImage.equidistantProjection((*it).image, calc, mSettings.getProjectionScale()); - - if(!equidistant.empty()) { - ThreatImage::drawWatermark(equidistant, dateStr); - saveImage(mSettings.getOutputPath() + std::string("equidistant_") + fileName, equidistant); - } else { - std::cout << "Failed to create equidistant projection" << std::endl; - } - }); + cv::Mat equidistant = equdistantProjector.project(img.image); + ThreatImage::drawWatermark(equidistant, dateStr, satelliteName); + const std::string filePath = mSettings.getOutputPath() + std::string("equidistant_") + fileName; + std::cout << "Saving " << filePath << std::endl; + saveImage(filePath, equidistant); } } - std::cout << "Generate images" << std::endl; - mThreadPool.waitForAllJobsDone(); - std::cout << "Generate images done" << std::endl; - imagesToSpread.clear(); + std::cout << "Save images done" << std::endl; } std::cout << "Generate composite images" << std::endl; std::time_t now = std::time(nullptr); std::stringstream compositeFileNameDateSS; compositeFileNameDateSS << std::put_time(std::localtime(&now), "%Y-%m-%d-%H-%M-%S"); + std::list equidistantTransform; + std::list mercatorTransform; + ImageSearchResult images = searchForImages(); + if(images.geolocationCalculators.size() > 1) { + if(mSettings.compositeEquadistantProjection()) { + equidistantTransform = ProjectImage::createCompositeProjector(ProjectImage::Projection::Equidistant, images.geolocationCalculators, mSettings.getCompositeProjectionScale()); + auto imgSizeIt = images.imageSizes.begin(); + for(auto& transform : equidistantTransform) { + std::cout << "Calculate Composite Equidistant TPS" << std::endl; + transform.calculateTransformation(*imgSizeIt++); + std::cout << "Calculate Composite Equidistant TPS done" << std::endl; + } + } + if(mSettings.compositeMercatorProjection()) { + mercatorTransform = ProjectImage::createCompositeProjector(ProjectImage::Projection::Mercator, images.geolocationCalculators, mSettings.getCompositeProjectionScale()); + auto imgSizeIt = images.imageSizes.begin(); + for(auto& transform : mercatorTransform) { + std::cout << "Calculate Composite Mercator TPS" << std::endl; + transform.calculateTransformation(*imgSizeIt++); + std::cout << "Calculate Composite Mercator TPS done" << std::endl; + } + } + } - if(mSettings.generateComposite123()) { - std::list images123; - std::list geolocationCalculators123; - searchForImages(images123, geolocationCalculators123, "123"); - - if(images123.size() > 1 && images123.size() == geolocationCalculators123.size()) { + if(mSettings.generateComposite221()) { + if(images.images221.size() > 1) { if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { - for(auto& img : images123) { + for(auto& img : images.images221) { img = ThreatImage::sharpen(img); } } - SpreadImage spreadImage; if(mSettings.compositeEquadistantProjection()) { - cv::Mat composite = spreadImage.equidistantProjection(images123, geolocationCalculators123, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate equidistant channel 123 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_123_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = equidistantTransform.begin(); + for(auto& img : images.images221) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_221_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } + if(mSettings.compositeMercatorProjection()) { - cv::Mat composite = spreadImage.mercatorProjection(images123, geolocationCalculators123, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate mercator channel 123 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_123_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = mercatorTransform.begin(); + for(auto& img : images.images221) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_221_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } } } - if(mSettings.generateComposite125()) { - std::list images125; - std::list geolocationCalculators125; - searchForImages(images125, geolocationCalculators125, "125"); + if(mSettings.generateComposite321()) { + if(images.images321.size() > 1) { + if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { + for(auto& img : images.images321) { + img = ThreatImage::sharpen(img); + } + } - if(images125.size() > 1 && images125.size() == geolocationCalculators125.size()) { - SpreadImage spreadImage; if(mSettings.compositeEquadistantProjection()) { - cv::Mat composite = spreadImage.equidistantProjection(images125, geolocationCalculators125, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate equidistant channel 125 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_125_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = equidistantTransform.begin(); + for(auto& img : images.images321) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_321_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } + if(mSettings.compositeMercatorProjection()) { - cv::Mat composite = spreadImage.mercatorProjection(images125, geolocationCalculators125, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate mercator channel 125 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_125_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = mercatorTransform.begin(); + for(auto& img : images.images321) { + + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_321_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } } } - if(mSettings.generateComposite221()) { - std::list images221; - std::list geolocationCalculators221; - searchForImages(images221, geolocationCalculators221, "221"); + if(mSettings.generateComposite125()) { + if(images.images125.size() > 1) { + if(mSettings.compositeEquadistantProjection()) { + std::list imagesToBlend; + auto transformIt = equidistantTransform.begin(); + for(auto& img : images.images125) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_125_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); + } + + if(mSettings.compositeMercatorProjection()) { + std::list imagesToBlend; + auto transformIt = mercatorTransform.begin(); + for(auto& img : images.images125) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_125_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); + } + } + } - if(images221.size() > 1 && images221.size() == geolocationCalculators221.size()) { - SpreadImage spreadImage; + if(mSettings.generateComposite224()) { + if(images.images224.size() > 1) { if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { - for(auto& img : images221) { + for(auto& img : images.images224) { img = ThreatImage::sharpen(img); } } if(mSettings.compositeEquadistantProjection()) { - cv::Mat composite = spreadImage.equidistantProjection(images221, geolocationCalculators221, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate equidistant channel 221 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_221_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = equidistantTransform.begin(); + for(const auto& img : images.images224) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_224_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } + if(mSettings.compositeMercatorProjection()) { - cv::Mat composite = spreadImage.mercatorProjection(images221, geolocationCalculators221, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate mercator channel 221 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_221_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = mercatorTransform.begin(); + for(auto& img : images.images224) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_224_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } } } + if(mSettings.generateComposite68()) { - std::list images68; - std::list geolocationCalculators68; - searchForImages(images68, geolocationCalculators68, "68"); + if(images.images68.size() > 1) { + std::list irImages; + if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { + for(const auto& img : images.images68) { + auto ir = ThreatImage::invertIR(img); + ir = ThreatImage::gamma(ir, 1.4); + ir = ThreatImage::contrast(ir, 1.3, -40); + ir = ThreatImage::sharpen(ir); + irImages.emplace_back(ir); + } + } - if(images68.size() > 1 && images68.size() == geolocationCalculators68.size()) { + if(mSettings.compositeEquadistantProjection()) { + std::list imagesToBlend; + auto transformIt = equidistantTransform.begin(); + for(const auto& img : irImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_68_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); + } + + if(mSettings.compositeMercatorProjection()) { + std::list imagesToBlend; + auto transformIt = mercatorTransform.begin(); + for(const auto& img : irImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_68_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); + } + } + + if(images.images67.size() > 1) { + std::list irImages; if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { - for(auto& img : images68) { - img = ThreatImage::invertIR(img); - img = ThreatImage::gamma(img, 1.4); - img = ThreatImage::contrast(img, 1.3, -40); - img = ThreatImage::sharpen(img); + for(const auto& img : images.images67) { + auto irImage = ThreatImage::equalize(img); + irImage = ThreatImage::invertIR(irImage); + irImages.emplace_back(irImage); } } - SpreadImage spreadImage; if(mSettings.compositeEquadistantProjection()) { - cv::Mat composite = spreadImage.equidistantProjection(images68, geolocationCalculators68, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate equidistant channel 68 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_68_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = equidistantTransform.begin(); + for(const auto& img : irImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_67_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } + if(mSettings.compositeMercatorProjection()) { - cv::Mat composite = spreadImage.mercatorProjection(images68, geolocationCalculators68, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate mercator channel 68 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_68_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = mercatorTransform.begin(); + for(const auto& img : irImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_67_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } } } - if(mSettings.generateComposite68Rain()) { - std::list images68; - std::list geolocationCalculators68; - searchForImages(images68, geolocationCalculators68, "68"); + if(mSettings.generateCompositeThermal()) { + cv::Mat thermalRef = cv::imread(mSettings.getResourcesPath() + "thermal_ref.bmp"); + if(images.images68.size() > 1) { + std::list thermalImages; + if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { + for(const auto& img : images.images68) { + thermalImages.emplace_back(ThreatImage::irToTemperature(img, thermalRef)); + } + } - if(images68.size() > 1 && images68.size() == geolocationCalculators68.size()) { + if(mSettings.compositeEquadistantProjection()) { + std::list imagesToBlend; + auto transformIt = equidistantTransform.begin(); + for(const auto& img : thermalImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_68_thermal_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); + } + + if(mSettings.compositeMercatorProjection()) { + std::list imagesToBlend; + auto transformIt = mercatorTransform.begin(); + for(const auto& img : thermalImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_68_thermal_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); + } + } + if(images.images67.size() > 1) { + std::list thermalImages; if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { - cv::Mat rainRef = cv::imread(mSettings.getResourcesPath() + "rain.bmp"); - for(auto& img : images68) { - cv::Mat rainOverlay = ThreatImage::irToRain(img, rainRef); - img = ThreatImage::invertIR(img); - img = ThreatImage::gamma(img, 1.4); - img = ThreatImage::contrast(img, 1.3, -40); - img = ThreatImage::sharpen(img); - img = ThreatImage::addRainOverlay(img, rainOverlay); + for(const auto& img : images.images67) { + thermalImages.emplace_back(ThreatImage::irToTemperature(ThreatImage::equalize(img), thermalRef)); } } - SpreadImage spreadImage; if(mSettings.compositeEquadistantProjection()) { - cv::Mat composite = spreadImage.equidistantProjection(images68, geolocationCalculators68, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate equidistant channel 68 rain composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_68_rain_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = equidistantTransform.begin(); + for(const auto& img : thermalImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_67_thermal_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } + if(mSettings.compositeMercatorProjection()) { - cv::Mat composite = spreadImage.mercatorProjection(images68, geolocationCalculators68, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate mercator channel 68 rain composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_68_rain_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = mercatorTransform.begin(); + for(const auto& img : thermalImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_67_thermal_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } } } - if(mSettings.generateCompositeThermal()) { - std::list images68; - std::list geolocationCalculators68; - searchForImages(images68, geolocationCalculators68, "68"); + if(mSettings.generateComposite68Rain()) { + if(images.images68.size() > 1) { + std::list irImages; + cv::Mat rainRef = cv::imread(mSettings.getResourcesPath() + "rain.bmp"); + if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { + for(const auto& img : images.images68) { + cv::Mat rainOverlay = ThreatImage::irToRain(img, rainRef); + cv::Mat ir = ThreatImage::invertIR(img); + ir = ThreatImage::gamma(ir, 1.4); + ir = ThreatImage::contrast(ir, 1.3, -40); + ir = ThreatImage::sharpen(ir); + ir = ThreatImage::addRainOverlay(ir, rainOverlay); + irImages.emplace_back(ir); + } + } + + if(mSettings.compositeEquadistantProjection()) { + std::list imagesToBlend; + auto transformIt = equidistantTransform.begin(); + for(const auto& img : irImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_68_rain_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); + } - if(images68.size() > 1 && images68.size() == geolocationCalculators68.size()) { + if(mSettings.compositeMercatorProjection()) { + std::list imagesToBlend; + auto transformIt = mercatorTransform.begin(); + for(const auto& img : irImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_68_rain_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); + } + } + + if(images.images67.size() > 1) { + std::list irImages; + cv::Mat rainRef = cv::imread(mSettings.getResourcesPath() + "rain.bmp"); if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { - cv::Mat thermalRef = cv::imread(mSettings.getResourcesPath() + "thermal_ref.bmp"); - for(auto& img : images68) { - img = ThreatImage::irToTemperature(img, thermalRef); + for(const auto& img : images.images67) { + cv::Mat rainOverlay = ThreatImage::irToRain(img, rainRef); + irImages.emplace_back(ThreatImage::addRainOverlay(ThreatImage::invertIR(ThreatImage::equalize(img)), rainOverlay)); } } - SpreadImage spreadImage; if(mSettings.compositeEquadistantProjection()) { - cv::Mat composite = spreadImage.equidistantProjection(images68, geolocationCalculators68, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate equidistant thermal composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_thermal_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = equidistantTransform.begin(); + for(const auto& img : irImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_67_rain_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } + if(mSettings.compositeMercatorProjection()) { - cv::Mat composite = spreadImage.mercatorProjection(images68, geolocationCalculators68, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate mercator thermal composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_thermal_composite.jpg", composite); + std::list imagesToBlend; + auto transformIt = mercatorTransform.begin(); + for(const auto& img : irImages) { + imagesToBlend.emplace_back(transformIt->project(img)); + transformIt++; + } + cv::Mat composite = BlendImages::merge(imagesToBlend); + const std::string filePath = mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_67_rain_composite." + mSettings.getOutputFormat(); + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } } } @@ -520,47 +868,95 @@ int main(int argc, char* argv[]) { return 0; } -void searchForImages(std::list& imagesOut, std::list& geolocationCalculatorsOut, const std::string& channelName) { +// 221, 321, 125, 224, 68, 67 +ImageSearchResult searchForImages() { + ImageSearchResult result; std::time_t now = std::time(nullptr); std::map> map; + TleReader reader(mSettings.getTlePath()); + reader.processFile(); + + std::map entriesSortByTime; for(const auto& entry : fs::directory_iterator(mSettings.getOutputPath())) { auto ftime = fs::last_write_time(entry); std::time_t cftime = std::chrono::system_clock::to_time_t((ftime)); - std::time_t fileCreatedSec = now - cftime; + if(entry.path().extension() == ".dat") { + entriesSortByTime[cftime] = entry; + // std::cout << "Enty" << entry << ", time=" << cftime << std::endl; + } + } - if(entry.path().extension() == ".gcp" && fileCreatedSec < (mSettings.getCompositeMaxAgeHours() * 3600)) { + for(const auto& [timeStamp, entry] : entriesSortByTime) { + std::time_t fileCreatedSec = now - timeStamp; + if(entry.path().extension() == ".dat" && fileCreatedSec < (mSettings.getCompositeMaxAgeHours() * 3600)) { std::string folder = entry.path().parent_path().generic_string(); std::string gcpFileName = entry.path().filename().generic_string(); std::string fileNameBase = gcpFileName.substr(0, gcpFileName.size() - 4); + std::string pathBase = folder + "/" + fileNameBase; + cv::Mat img; - do { - fs::path fileJPG(folder + "/" + fileNameBase + "_" + channelName + ".jpg"); + fs::path imagePath(pathBase + "_221.bmp"); + if(fs::exists(imagePath)) { + img = cv::imread(imagePath.generic_string()); + result.images221.emplace_back(img); + } + imagePath = pathBase + "_321.bmp"; + if(fs::exists(imagePath)) { + img = cv::imread(imagePath.generic_string()); + result.images321.emplace_back(img); + } + imagePath = pathBase + "_125.bmp"; + if(fs::exists(imagePath)) { + img = cv::imread(imagePath.generic_string()); + result.images125.emplace_back(img); + } + imagePath = pathBase + "_224.bmp"; + if(fs::exists(imagePath)) { + img = cv::imread(imagePath.generic_string()); + result.images224.emplace_back(img); + } + imagePath = pathBase + "_67.bmp"; + if(fs::exists(imagePath)) { + img = cv::imread(imagePath.generic_string()); + result.images67.emplace_back(img); + } + imagePath = pathBase + "_68.bmp"; + if(fs::exists(imagePath)) { + img = cv::imread(imagePath.generic_string()); + result.images68.emplace_back(img); + } - if(fs::exists(fileJPG)) { - map[cftime] = std::make_tuple(entry.path().generic_string(), fileJPG.generic_string()); - break; + if(!img.empty()) { + result.imageSizes.emplace_back(img.size()); + std::ifstream datFileStream(entry.path().generic_string()); + if(!datFileStream) { + continue; } - - fs::path fileBMP(folder + "/" + fileNameBase + "_" + channelName + ".bmp"); - - if(fs::exists(fileBMP)) { - map[cftime] = std::make_tuple(entry.path().generic_string(), fileBMP.generic_string()); - - break; + std::string linePassStart; + std::string linePassLength; + datFileStream >> result.satelliteName; + datFileStream >> linePassStart; + datFileStream >> linePassLength; + datFileStream.close(); + int64_t ticks = std::stoll(linePassStart); + DateTime passStart(ticks); + ticks = std::stoll(linePassLength); + TimeSpan passLength(ticks); + + auto projectionSetting = mSettings.getProjectionSetting(result.satelliteName); + + TleReader::TLE tle; + if(!reader.getTLE(projectionSetting.satelliteNameInTLE, tle)) { + std::cout << "TLE data not found in TLE file, unable to create composite images..." << std::endl; + return result; } - } while(false); - } - } - - if(map.size() > 1) { - for(auto const& [time, paths] : map) { - std::cout << std::get<1>(paths) << std::endl; - - geolocationCalculatorsOut.emplace_back(PixelGeolocationCalculator::load(std::get<0>(paths))); - imagesOut.emplace_back(cv::imread(std::get<1>(paths))); + PixelGeolocationCalculator calc(tle, passStart, passLength, projectionSetting.scanAngle, projectionSetting.roll, projectionSetting.pitch, projectionSetting.yaw, img.size().width, img.size().height); + result.geolocationCalculators.emplace_back(calc); + } } } + return result; } void saveImage(const std::string fileName, const cv::Mat& image) { diff --git a/resources/ShapeFiles/ne_10m_admin_0_countries.dbf b/resources/ShapeFiles/ne_10m_admin_0_countries.dbf new file mode 100644 index 0000000..dedacfc Binary files /dev/null and b/resources/ShapeFiles/ne_10m_admin_0_countries.dbf differ diff --git a/resources/ShapeFiles/ne_10m_admin_0_countries.shp b/resources/ShapeFiles/ne_10m_admin_0_countries.shp new file mode 100644 index 0000000..a10056b Binary files /dev/null and b/resources/ShapeFiles/ne_10m_admin_0_countries.shp differ diff --git a/resources/ShapeFiles/ne_10m_admin_0_countries.shx b/resources/ShapeFiles/ne_10m_admin_0_countries.shx new file mode 100644 index 0000000..9e417c2 Binary files /dev/null and b/resources/ShapeFiles/ne_10m_admin_0_countries.shx differ diff --git a/resources/ShapeFiles/ne_10m_coastline.shp b/resources/ShapeFiles/ne_10m_coastline.shp new file mode 100644 index 0000000..2066748 Binary files /dev/null and b/resources/ShapeFiles/ne_10m_coastline.shp differ diff --git a/resources/kernels/tps.cl b/resources/kernels/tps.cl new file mode 100644 index 0000000..6751685 --- /dev/null +++ b/resources/kernels/tps.cl @@ -0,0 +1,58 @@ +// featuring floating point bit level hacking, Remez algorithm +// x=m*2^p => ln(x)=ln(m)+ln(2)p +inline float ln(float x) { + uint bx = *(uint*)(&x); + const uint ex = bx >> 23; + int t = (int)ex - (int)127; + uint s = (t < 0) ? (-t) : t; + bx = 1065353216 | (bx & 8388607); + x = *(float*)(&bx); + return -1.7417939f + (2.8212026f + (-1.4699568f + (0.44717955f - 0.056570851f * x) * x) * x) * x + 0.6931471806f * t; +} + +inline float pow2(const float x) { + return x * x; +} + +inline float dist(const float x1, const float y1, const float x2, const float y2) { + float norm = pow2(x1 - x2) + pow2(y1 - y2); + if(norm < 0) { + norm = 0; + } + // return sqrt(norm); + return norm * log(norm + FLT_EPSILON); + // return norm * ln(norm + FLT_EPSILON); +} + +void kernel apply(global const float* shapeRef, global const float* tpsParameters, global const float* point, const uint shapeRefRows, const uint tpsParamRows, const uint cols, global float* out) { + int id = get_global_id(0); + int size = get_global_size(0); + + int itemPerTask = (cols + size - 1) / size; + int start = id * itemPerTask; + int end = (id + 1) * itemPerTask; + if(end > cols) { + end = cols; + }; + + for(int x = start; x < end; x++) { + for(int i = 0; i < 2; i++) { + const float a1 = tpsParameters[(tpsParamRows - 3) * 2 + i]; + const float ax = tpsParameters[(tpsParamRows - 2) * 2 + i]; + const float ay = tpsParameters[(tpsParamRows - 1) * 2 + i]; + + const float affine = a1 + ax * point[x * 2] + ay * point[x * 2 + 1]; + float nonrigid = 0; + + for(int j = 0; j < shapeRefRows; j++) { + nonrigid += tpsParameters[j * 2 + i] * dist(shapeRef[j * 2], shapeRef[j * 2 + 1], point[x * 2], point[x * 2 + 1]); + } + if(i == 0) { + out[x * 2] = affine + nonrigid; + } + if(i == 1) { + out[x * 2 + 1] = affine + nonrigid; + } + } + } +} \ No newline at end of file diff --git a/resources/settings.ini b/resources/settings.ini index 5bb7623..b7615d9 100644 --- a/resources/settings.ini +++ b/resources/settings.ini @@ -10,9 +10,10 @@ ProjectionScale=0.75 CompositeProjectionScale=0.7 CompositeAzimuthalEquidistantProjection=true CompositeMercatorProjection=false -GenerateComposite123=true +GenerateComposite321=true GenerateComposite125=true GenerateComposite221=true +GenerateComposite224=true GenerateComposite68=true GenerateCompositeThermal=true GenerateComposite68Rain=true @@ -28,12 +29,21 @@ TimeOffset=0 [METEOR-M-2-3] SatNameInTLE=METEOR-M2 3 -ScanAngle=110.3 -Roll=0 -Pitch=0 +ScanAngle=109.55 +Roll=0.3 +Pitch=0.8 Yaw=0 # Onboard time offset in seconds -TimeOffset=0 +TimeOffset=0.6 + +[METEOR-M-2-4] +SatNameInTLE=METEOR-M2 4 +ScanAngle=109.5 +Roll=0 +Pitch=0.8 +Yaw=-0.3 +# Onboard time offset in seconds +TimeOffset=0.0 [Demodulator] ;Increasing bandwidth allow the pll to lock faster on signal, while decreasing it makes the pll lock more stable @@ -54,8 +64,8 @@ Color=#AD880C Size=80 # Font line width in pixels Width=3 -;Text options: use '%date%' for date, '\n' for new line -Text=MeteorDemod by Digitelektro\nMETEOR M2-3\n%date% +;Text options: use '%date%' for date %sat% for satellite name, %version% for current version of meteordemod '\n' for new line +Text=MeteorDemod %version% by Digitelektro\n%sat%\n%date% [ReceiverLocation] Draw=true @@ -73,12 +83,12 @@ Color=#C8C8C8 Thickness=1 [ShapeFileCoastLines] -FileName=ShapeFiles/ne_50m_coastline.shp +FileName=ShapeFiles/ne_10m_coastline.shp Color=#808000 Thickness=2 [ShapeFileBoundaryLines] -FileName=ShapeFiles/ne_50m_admin_0_boundary_lines_land.shp +FileName=ShapeFiles/ne_10m_admin_0_countries.shp Color=#C8C8C8 Thickness=2 diff --git a/tools/opencl.cpp b/tools/opencl.cpp new file mode 100644 index 0000000..b35e53b --- /dev/null +++ b/tools/opencl.cpp @@ -0,0 +1,114 @@ +#include "opencl.h" + +#include + +#ifdef OPENCL_FOUND + +namespace OpenCL { + +std::vector Context::getDevices(cl_device_type deviceType) { + std::vector devices; + cl_platform_id platformsIDs[cMaximumDeviceCount]; + cl_uint platformsCount = 0; + cl_device_id deviceIDs[cMaximumDeviceCount]; + cl_uint deviceCount = 0; + char deviceName[cMaximumDeviceNameLength]; + size_t deviceNameLen = 0; + + cl_int error = clGetPlatformIDs(cMaximumDeviceCount, platformsIDs, &platformsCount); + if(error == CL_SUCCESS) { + for(cl_uint p = 0; p < platformsCount; p++) { + if(clGetDeviceIDs(platformsIDs[p], deviceType, cMaximumDeviceCount, deviceIDs, &deviceCount) == CL_SUCCESS) { + for(cl_uint d = 0; d < deviceCount; d++) { + size_t workGroupSize = 0; + size_t maxComputeUnits = 0; + cl_device_type deviceType; + + error = clGetDeviceInfo(deviceIDs[d], CL_DEVICE_TYPE, sizeof(cl_device_type), &deviceType, NULL); + if(error != CL_SUCCESS) { + continue; + } + error = clGetDeviceInfo(deviceIDs[d], CL_DEVICE_MAX_WORK_GROUP_SIZE, sizeof(size_t), &workGroupSize, NULL); + if(error != CL_SUCCESS) { + continue; + } + error = clGetDeviceInfo(deviceIDs[d], CL_DEVICE_MAX_COMPUTE_UNITS, sizeof(size_t), &maxComputeUnits, NULL); + if(error != CL_SUCCESS) { + continue; + } + error = clGetDeviceInfo(deviceIDs[d], CL_DEVICE_NAME, cMaximumDeviceNameLength, deviceName, &deviceNameLen); + if(error != CL_SUCCESS) { + continue; + } + + devices.emplace_back(Device{platformsIDs[p], deviceIDs[d], deviceType, std::string(&deviceName[0], &deviceName[deviceNameLen]), workGroupSize, maxComputeUnits}); + } + } + } + } + + return devices; +} + +Context::~Context() { + if(mContext) { + clReleaseContext(mContext); + } +} + +void Context::init() { + cl_int error = 0; + + mContext = clCreateContext(NULL, 1, &mDevice.deviceID, NULL, NULL, &error); + if(error != CL_SUCCESS) { + throw std::runtime_error("Could not init OpenCL context!"); + } +} + +cl_program Context::buildKernel(const std::string& path) { + std::ifstream ifStream(path); + std::string kernelSrc(std::istreambuf_iterator{ifStream}, {}); + + const char* srcs[1] = {kernelSrc.c_str()}; + const size_t lens[1] = {kernelSrc.length()}; + cl_int err = 0; + char errorMessage[10000]; + size_t errorMessgeLen = 0; + + cl_program program = clCreateProgramWithSource(mContext, 1, srcs, lens, &err); + err = clBuildProgram(program, 1, &mDevice.deviceID, NULL, NULL, NULL); + + if(err != CL_SUCCESS) { + if(clGetProgramBuildInfo(program, mDevice.deviceID, CL_PROGRAM_BUILD_LOG, sizeof(errorMessage), errorMessage, &errorMessgeLen) == CL_SUCCESS) { + throw std::runtime_error("OpenCV Build program error: " + std::string(&errorMessage[0], &errorMessage[errorMessgeLen])); + } else { + throw std::runtime_error("OpenCV Build program error!"); + } + } + + return program; +} + +cl_program Context::buildKernel(const char* kernelSrc, size_t srcLen) { + cl_int err = 0; + char errorMessage[10000]; + size_t errorMessgeLen = 0; + const char* srcs[1] = {kernelSrc}; + const size_t lens[1] = {srcLen}; + cl_program program = clCreateProgramWithSource(mContext, 1, srcs, lens, &err); + err = clBuildProgram(program, 1, &mDevice.deviceID, NULL, NULL, NULL); + + if(err != CL_SUCCESS) { + if(clGetProgramBuildInfo(program, mDevice.deviceID, CL_PROGRAM_BUILD_LOG, sizeof(errorMessage), errorMessage, &errorMessgeLen) == CL_SUCCESS) { + throw std::runtime_error("OpenCV Build program error: " + std::string(&errorMessage[0], &errorMessage[errorMessgeLen])); + } else { + throw std::runtime_error("OpenCV Build program error!"); + } + } + + return program; +} + +} // namespace OpenCL + +#endif // OPENCL_FOUND diff --git a/tools/opencl.h b/tools/opencl.h new file mode 100644 index 0000000..e05d51a --- /dev/null +++ b/tools/opencl.h @@ -0,0 +1,172 @@ +#pragma once + +#ifdef OPENCL_FOUND + +#define CL_TARGET_OPENCL_VERSION 300 + +#include + +#include +#include +#include +#include +#include + +namespace OpenCL { + +class Context { + friend class Program; + + public: + struct Device { + cl_platform_id platformID; + cl_device_id deviceID; + cl_device_type type; + std::string name; + size_t workGroupSize; + size_t maxComputeUnits; + }; + + public: + static std::vector getDevices(cl_device_type deviceType = CL_DEVICE_TYPE_ALL); + + public: + public: + Context() = delete; + Context(const Device& device) + : mDevice(device) {} + + ~Context(); + + void init(); + cl_program buildKernel(const std::string& path); + cl_program buildKernel(const char* kernelSrc, size_t srcLen); + + size_t getWorkGroupSize() const { + return mDevice.workGroupSize; + } + size_t getMaxComputeUnits() const { + return mDevice.workGroupSize; + } + + private: + Device mDevice; + cl_context mContext = nullptr; + + + private: + static constexpr int cMaximumDeviceCount = 100; + static constexpr int cMaximumDeviceNameLength = 300; +}; + +class Program { + public: + explicit Program(std::shared_ptr context, const std::string& path, const std::string& name) + : mContext(context) { + cl_int err = 0; + mProgram = context->buildKernel(path); + + mQueue = clCreateCommandQueueWithProperties(mContext->mContext, mContext->mDevice.deviceID, 0, &err); + if(err != CL_SUCCESS) { + throw std::runtime_error("Create queue failed!"); + } + mKernel = clCreateKernel(mProgram, name.c_str(), &err); + if(err != CL_SUCCESS) { + throw std::runtime_error("Create kernel failed!"); + } + } + + explicit Program(std::shared_ptr context, const char* kernelSrc, size_t srcLen, const std::string& name) + : mContext(context) { + cl_int err = 0; + mProgram = context->buildKernel(kernelSrc, srcLen); + + mQueue = clCreateCommandQueueWithProperties(mContext->mContext, mContext->mDevice.deviceID, 0, &err); + if(err != CL_SUCCESS) { + throw std::runtime_error("Create queue failed!"); + } + mKernel = clCreateKernel(mProgram, name.c_str(), &err); + if(err != CL_SUCCESS) { + throw std::runtime_error("Create kernel failed!"); + } + } + + ~Program() { + clFlush(mQueue); + clFinish(mQueue); + clReleaseKernel(mKernel); + clReleaseProgram(mProgram); + for(auto& [index, memory] : mBuffers) { + clReleaseMemObject(memory); + } + clReleaseCommandQueue(mQueue); + } + + std::shared_ptr getContext() { + return mContext; + } + + template + inline void createBuffer(cl_mem_flags flags, size_t length, int index) { + cl_int err = 0; + cl_mem buffer = clCreateBuffer(mContext->mContext, flags, sizeof(T) * length, NULL, &err); + if(err != CL_SUCCESS) { + throw std::runtime_error("Create OpenCL buffer failed!"); + } + mBuffers[index] = buffer; + + err = clSetKernelArg(mKernel, index, sizeof(cl_mem), &buffer); + if(err != CL_SUCCESS) { + throw std::runtime_error("Set kernel arg failed!"); + } + } + + template + inline void enqueueWriteBuffer(const T* params, size_t length, size_t paramIndex) { + cl_int err = 0; + err = clEnqueueWriteBuffer(mQueue, mBuffers[paramIndex], CL_TRUE, 0, sizeof(T) * length, params, 0, NULL, NULL); + if(err != CL_SUCCESS) { + throw std::runtime_error("Write OpenCL buffer failed!"); + } + } + + template + inline void enqueueReadBuffer(T* params, size_t length, size_t paramIndex) { + cl_int err = 0; + err = clEnqueueReadBuffer(mQueue, mBuffers[paramIndex], true, 0, sizeof(T) * length, params, 0, NULL, NULL); + if(err != CL_SUCCESS) { + throw std::runtime_error("Write OpenCL buffer failed!"); + } + } + + template + inline void setKernelArg(T arg, cl_uint index) { + cl_int err = 0; + clSetKernelArg(mKernel, index, sizeof(T), &arg); + if(err != CL_SUCCESS) { + throw std::runtime_error("Set kernel arg failed!"); + } + } + + void execute(size_t itemSize) { + size_t globalItemSize = itemSize; + cl_int err = clEnqueueNDRangeKernel(mQueue, mKernel, 1, NULL, &globalItemSize, NULL, 0, NULL, NULL); + if(err != CL_SUCCESS) { + throw std::runtime_error("Couldn't execute kernel!"); + } + } + + + private: + std::shared_ptr mContext; + cl_program mProgram = nullptr; + cl_command_queue mQueue = nullptr; + cl_kernel mKernel = nullptr; + std::map mBuffers; +}; + + +} // namespace OpenCL + + +#endif // OPENCL_FOUND diff --git a/tools/pixelgeolocationcalculator.cpp b/tools/pixelgeolocationcalculator.cpp index c0b37cb..fda40fe 100644 --- a/tools/pixelgeolocationcalculator.cpp +++ b/tools/pixelgeolocationcalculator.cpp @@ -1,103 +1,87 @@ #include "pixelgeolocationcalculator.h" -#include - #include "settings.h" -PixelGeolocationCalculator PixelGeolocationCalculator::load(const std::string& path) { - Settings& settings = Settings::getInstance(); - TleReader reader(settings.getTlePath()); - reader.processFile(); - TleReader::TLE tle; - reader.getTLE(settings.getSatNameInTLE(), tle); - PixelGeolocationCalculator calc(tle, DateTime(), TimeSpan(0), settings.getScanAngle(), settings.getRoll(), settings.getPitch(), settings.getYaw()); - std::ifstream gcpReader(path); - - if(!gcpReader) { - std::cout << "Open GCP file failed"; - return calc; - } - - calc.mCoordinates.clear(); - int i, n; - double longitude, latitude; - while(gcpReader >> i >> n >> latitude >> longitude) { - calc.mCoordinates.push_back(CoordGeodetic(latitude, longitude, 0, false)); - } - - return calc; -} - -PixelGeolocationCalculator::PixelGeolocationCalculator(const TleReader::TLE& tle, const DateTime& passStart, const TimeSpan& passLength, double scanAngle, double roll, double pitch, double yaw, int earthRadius, int satelliteAltitude) +PixelGeolocationCalculator::PixelGeolocationCalculator(const TleReader::TLE& tle, + const DateTime& passStart, + const TimeSpan& passLength, + double scanAngle, + double roll, + double pitch, + double yaw, + int imageWidth, + int imageHeight, + int earthRadius, + int satelliteAltitude) : mTle(tle.satellite, tle.line1, tle.line2) , mSgp4(mTle) , mPassStart(passStart) , mPassLength(passLength) , mScanAngle(scanAngle) - , mRoll(roll) - , mPitch(pitch) - , mYaw(yaw) + , mRollOffset(roll) + , mPitchOffset(pitch) + , mYawOffset(yaw) + , mImageWidth(imageWidth) + , mImageHeight(imageHeight) , mEarthradius(earthRadius) - , mSatelliteAltitude(satelliteAltitude) {} + , mSatelliteAltitude(satelliteAltitude) { + mCenterCoordinate = getCoordinateAt(mImageWidth / 2, mImageHeight / 2); +} -void PixelGeolocationCalculator::calcPixelCoordinates() { - double angle; - CoordGeodetic satOnGroundPrev; - CoordGeodetic satOnGround; - std::vector passList; +CoordGeodetic PixelGeolocationCalculator::getCoordinateAt(unsigned int x, unsigned int y) const { + if(x > mImageWidth) { + x = mImageWidth; + } + if(y > mImageHeight) { + y = mImageHeight; + } - DateTime passEnd = DateTime(mPassStart).Add(mPassLength); + DateTime currentTime = mPassStart; + currentTime = currentTime.AddMicroseconds(PIXELTIME_MS * 1000 * y); + Eci position = mSgp4.FindPosition(currentTime.AddMicroseconds(PIXELTIME_MS * 1000 * 10)); + Eci prevPosition = mSgp4.FindPosition(currentTime); + CoordGeodetic satOnGround = Eci(currentTime, position.Position(), position.Velocity()).ToGeodetic(); + CoordGeodetic satOnGroundPrev = Eci(currentTime, prevPosition.Position(), prevPosition.Velocity()).ToGeodetic(); - for(DateTime currentTime = mPassStart; currentTime <= passEnd; currentTime = currentTime.AddMicroseconds(PIXELTIME_MS * 10 * 1000)) { - passList.push_back(mSgp4.FindPosition(currentTime)); - } + double angle = calculateBearingAngle(satOnGround, satOnGroundPrev); + angle = degreeToRadian(90) - angle; - std::vector::const_iterator it = passList.begin(); - std::vector::const_iterator prevIt = passList.begin(); - ++it; + double currentScanAngle = (mScanAngle / mImageWidth) * ((mImageWidth / 2.0) - x); + // std::cout << "CurrentScan Angle2 " << currentScanAngle << std::endl; - DateTime currentTime = mPassStart; - for(unsigned int i = 0; it != passList.end(); prevIt = it, ++it, i++) { - satOnGroundPrev = Eci(currentTime, prevIt->Position(), prevIt->Velocity()).ToGeodetic(); - satOnGround = Eci(currentTime, it->Position(), it->Velocity()).ToGeodetic(); - currentTime = currentTime.AddMicroseconds(PIXELTIME_MS * 10 * 1000); - - angle = calculateBearingAngle(satOnGround, satOnGroundPrev); - angle = degreeToRadian(90) - angle; - - for(int n = 79; n > -79; n--) { - Vector3 posVector(los_to_earth(Vector3(satOnGround), degreeToRadian((((mScanAngle / 2.0) / 79.0)) * n), 0, angle)); - mCoordinates.push_back(posVector.toCoordinate()); - } - } + Vector3 posVector(raytraceToEarth(Vector3(satOnGround), degreeToRadian(currentScanAngle), 0, angle)); + return posVector.toCoordinate(); } -void PixelGeolocationCalculator::save(const std::string& path) { - std::ofstream file(path); +const CoordGeodetic& PixelGeolocationCalculator::getCenterCoordinate() const { + return mCenterCoordinate; +} - if(!file) { - return; - } +CoordGeodetic PixelGeolocationCalculator::getCoordinateTopLeft() const { + return getCoordinateAt(0, 0); +} - std::vector::const_iterator it; - int i = 0; - int n = 0; +CoordGeodetic PixelGeolocationCalculator::getCoordinateTopRight() const { + return getCoordinateAt(mImageWidth, 0); +} - for(it = mCoordinates.begin(); it != mCoordinates.end(); ++it, n++) { - file << (n * 10) << " " << (i * 10) << " " << std::setprecision(15) << radioanToDegree(it->latitude) << " " << radioanToDegree(it->longitude) << std::endl; +CoordGeodetic PixelGeolocationCalculator::getCoordinateBottomLeft() const { + return getCoordinateAt(0, mImageHeight); +} - if(n == 157) { - n = -1; - i++; - } - } +CoordGeodetic PixelGeolocationCalculator::getCoordinateBottomRight() const { + return getCoordinateAt(mImageWidth, mImageHeight); +} - file.close(); +bool PixelGeolocationCalculator::isNorthBoundPass() const { + auto coord1 = getCoordinateAt(0, 0); + auto coord2 = getCoordinateAt(0, 100); + return (coord1.latitude < coord2.latitude); } -Vector PixelGeolocationCalculator::los_to_earth(const Vector& position, double roll, double pitch, double yaw) { +Vector PixelGeolocationCalculator::raytraceToEarth(const Vector& position, double roll, double pitch, double yaw) const { double a = 6371.0087714; double b = 6371.0087714; double c = 6356.752314245; @@ -110,9 +94,9 @@ Vector PixelGeolocationCalculator::los_to_earth(const Vector& position, double r Vector lookVector(0, 0, 0); Matrix4x4 lookMatrix = lookAt(position, lookVector, Vector(0, 0, 1)); - Matrix4x4 rotateX = Matrix4x4::CreateRotationX(roll + degreeToRadian(mRoll)); - Matrix4x4 rotateY = Matrix4x4::CreateRotationY(pitch + degreeToRadian(mPitch)); - Matrix4x4 rotateZ = Matrix4x4::CreateRotationZ(yaw + degreeToRadian(mYaw)); + Matrix4x4 rotateX = Matrix4x4::CreateRotationX(roll + degreeToRadian(mRollOffset)); + Matrix4x4 rotateY = Matrix4x4::CreateRotationY(pitch + degreeToRadian(mPitchOffset)); + Matrix4x4 rotateZ = Matrix4x4::CreateRotationZ(yaw + degreeToRadian(mYawOffset)); matrix = matrix * lookMatrix * rotateZ * rotateY * rotateX; Vector vector3(matrix.mElements[2], matrix.mElements[6], matrix.mElements[10]); diff --git a/tools/pixelgeolocationcalculator.h b/tools/pixelgeolocationcalculator.h index 6218d2e..b09f42a 100644 --- a/tools/pixelgeolocationcalculator.h +++ b/tools/pixelgeolocationcalculator.h @@ -39,43 +39,27 @@ class PixelGeolocationCalculator { PixelGeolocationCalculator(); public: - static PixelGeolocationCalculator load(const std::string& path); - - public: - PixelGeolocationCalculator(const TleReader::TLE& tle, const DateTime& passStart, const TimeSpan& passLength, double scanAngle, double roll, double pitch, double yaw, int earthRadius = 6378, int satelliteAltitude = 825); - - void calcPixelCoordinates(); - - void save(const std::string& path); - - public: - int getGeorefMaxImageHeight() const { - return (mCoordinates.size() / 158) * 10; - } - - const CoordGeodetic& getCenterCoordinate() const { - return mCoordinates[mCoordinates.size() / 2 + 79]; - } - - inline const CoordGeodetic& getCoordinateAt(unsigned int x, unsigned int y) const { - return mCoordinates[((x / 10)) + ((y / 10) * 158)]; - } - - inline const CoordGeodetic& getCoordinateTopLeft() const { - return mCoordinates[0]; - } - - inline const CoordGeodetic& getCoordinateTopRight() const { - return mCoordinates[157]; - } - - inline const CoordGeodetic& getCoordinateBottomLeft() const { - return mCoordinates[mCoordinates.size() - 158]; - } - - inline const CoordGeodetic& getCoordinateBottomRight() const { - return mCoordinates[mCoordinates.size() - 1]; - } + PixelGeolocationCalculator(const TleReader::TLE& tle, + const DateTime& passStart, + const TimeSpan& passLength, + double scanAngle, + double roll, + double pitch, + double yaw, + int imageWidth, + int imageHeight, + int earthRadius = 6378, + int satelliteAltitude = 825); + + + const CoordGeodetic& getCenterCoordinate() const; + CoordGeodetic getCoordinateAt(unsigned int x, unsigned int y) const; + CoordGeodetic getCoordinateTopLeft() const; + CoordGeodetic getCoordinateTopRight() const; + CoordGeodetic getCoordinateBottomLeft() const; + CoordGeodetic getCoordinateBottomRight() const; + + bool isNorthBoundPass() const; inline int getEarthRadius() const { return mEarthradius; @@ -144,9 +128,9 @@ class PixelGeolocationCalculator { private: void calculateCartesionCoordinates(); - Vector los_to_earth(const Vector& position, double roll, double pitch, double yaw); - double calculateBearingAngle(const CoordGeodetic& start, const CoordGeodetic& end); - Matrix4x4 lookAt(const Vector3& position, const Vector3& target, const Vector3& up); + Vector raytraceToEarth(const Vector& position, double roll, double pitch, double yaw) const; + static double calculateBearingAngle(const CoordGeodetic& start, const CoordGeodetic& end); + static Matrix4x4 lookAt(const Vector3& position, const Vector3& target, const Vector3& up); static inline double degreeToRadian(double degree) { return (M_PI * degree / 180.0); @@ -165,11 +149,15 @@ class PixelGeolocationCalculator { SGP4 mSgp4; DateTime mPassStart; TimeSpan mPassLength; - double mScanAngle, mRoll, mPitch, mYaw; + double mScanAngle; + double mRollOffset; + double mPitchOffset; + double mYawOffset; + int mImageWidth; + int mImageHeight; int mEarthradius; int mSatelliteAltitude; - std::vector mCoordinates; - + CoordGeodetic mCenterCoordinate; static constexpr double PIXELTIME_MINUTES = 0.02564876089324618736383442265795; // Just a rough calculation for every 10 pixel in minutes static constexpr double PIXELTIME_MS = 154.0;