Merge pull request #29 from Digitelektro/beta

Beta

CMakeLists.txt

@@ -139,8 +139,8 @@ if(UNIX AND NOT APPLE)
     if (DEBIAN_FOUND)
         SET(CPACK_GENERATOR "DEB")
         SET(CPACK_PACKAGE_VERSION_MAJOR "2")
-        SET(CPACK_PACKAGE_VERSION_MINOR "0")
-        SET(CPACK_PACKAGE_VERSION_PATCH "1")
+        SET(CPACK_PACKAGE_VERSION_MINOR "1")
+        SET(CPACK_PACKAGE_VERSION_PATCH "0")
         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")

GIS/shaperenderer.cpp

@@ -128,7 +128,7 @@ void GIS::ShapeRenderer::drawShapeMercator(cv::Mat &src, float xStart, float ySt
     }
 }
 
-void GIS::ShapeRenderer::drawShapeEquidistant(cv::Mat &src, float xStart, float yStart, float xCenter, float yCenter)
+void GIS::ShapeRenderer::drawShapeEquidistant(cv::Mat &src, float xStart, float yStart, float centerLatitude, float centerLongitude)
 {
     if(!load()) {
         return;
@@ -146,12 +146,12 @@ void GIS::ShapeRenderer::drawShapeEquidistant(cv::Mat &src, float xStart, float
                     for(polyLineIterator->begin(); *polyLineIterator != polyLineIterator->end(); ++(*polyLineIterator)) {
                         //std::cout << polyLineIterator->point.x << " " << polyLineIterator->point.y << std::endl;
 
-                        PixelGeolocationCalculator::CartesianCoordinateF coordinate = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(polyLineIterator->point.y, polyLineIterator->point.x, xCenter, yCenter, mEarthRadius + mAltitude);
+                        PixelGeolocationCalculator::CartesianCoordinateF coordinate = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(polyLineIterator->point.y, polyLineIterator->point.x, centerLatitude, centerLongitude, mEarthRadius + mAltitude);
 
                         coordinate.x += -xStart;
                         coordinate.y += -yStart;
 
-                        if(equidistantCheck(polyLineIterator->point.y, polyLineIterator->point.x, xCenter, yCenter)) {
+                        if(equidistantCheck(polyLineIterator->point.y, polyLineIterator->point.x, centerLatitude, centerLongitude)) {
                             polyLines.push_back(cv::Point2d(coordinate.x, coordinate.y));
                         }
                     }
@@ -174,11 +174,11 @@ void GIS::ShapeRenderer::drawShapeEquidistant(cv::Mat &src, float xStart, float
                 for(recordIterator->begin(); *recordIterator != recordIterator->end(); ++(*recordIterator)) {
                     ShapeReader::Point point(*recordIterator);
 
-                    PixelGeolocationCalculator::CartesianCoordinateF coordinate = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(point.y, point.x, xCenter, yCenter, mEarthRadius + mAltitude);
+                    PixelGeolocationCalculator::CartesianCoordinateF coordinate = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(point.y, point.x, centerLatitude, centerLongitude, mEarthRadius + mAltitude);
                     coordinate.x += -xStart;
                     coordinate.y += -yStart;
 
-                    if(equidistantCheck(point.y, point.x, xCenter, yCenter) == false) {
+                    if(equidistantCheck(point.y, point.x, centerLatitude, centerLongitude) == false) {
                         continue;
                     }
 
@@ -196,11 +196,11 @@ void GIS::ShapeRenderer::drawShapeEquidistant(cv::Mat &src, float xStart, float
                     ShapeReader::Point point(*recordIterator);
                     std::vector<std::string> fieldValues = dbFilereader.getFieldValues(i);
 
-                    PixelGeolocationCalculator::CartesianCoordinateF coordinate = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(point.y, point.x, xCenter, yCenter, mEarthRadius + mAltitude);
+                    PixelGeolocationCalculator::CartesianCoordinateF coordinate = PixelGeolocationCalculator::coordinateToAzimuthalEquidistantProjection(point.y, point.x, centerLatitude, centerLongitude, mEarthRadius + mAltitude);
                     coordinate.x += -xStart;
                     coordinate.y += -yStart;
 
-                    if(equidistantCheck(point.y, point.x, xCenter, yCenter) == false) {
+                    if(equidistantCheck(point.y, point.x, centerLatitude, centerLongitude) == false) {
                         continue;
                     }
 
@@ -241,34 +241,38 @@ void GIS::ShapeRenderer::drawShapeEquidistant(cv::Mat &src, float xStart, float
     }
 }
 
-bool GIS::ShapeRenderer::equidistantCheck(float x, float y, float centerLongitude, float centerLatitude)
+bool GIS::ShapeRenderer::equidistantCheck(float latitude, float longitude, float centerLatitude, float centerLongitude)
 {
+    bool longResult = true;
+    bool latResult = true;
+
     int minLongitude = static_cast<int>(centerLongitude - 90);
     int maxLongitude = static_cast<int>(centerLongitude + 90);
     int minLatitude = static_cast<int>(centerLatitude - 45);
     int maxLatitude = static_cast<int>(centerLatitude + 45);
 
     //Normalize
-    minLongitude = (minLongitude + 540) % 360 - 180;
-    maxLongitude = (maxLongitude + 540) % 360 - 180;
+    minLongitude = ((minLongitude + 540) % 360 - 180);
+    maxLongitude = ((maxLongitude + 540) % 360 - 180);
+    minLatitude = ((minLatitude + 270) % 180 - 90);
+    maxLatitude = ((maxLatitude + 270) % 180 - 90);
 
-    if(minLatitude < -90)
+    if(maxLatitude < minLatitude)
     {
-        minLatitude = ((minLatitude + 270) % 180 - 90) * -1;
+        latResult = latitude > minLatitude || latitude < maxLatitude;
     }
-    if(maxLatitude > 90)
+    else
     {
-        maxLatitude = ((maxLatitude + 270) % 180 - 90) * -1;
+        latResult = latitude > minLatitude && latitude < maxLatitude;
     }
-
-
-
-    if (x < minLongitude || x > maxLongitude || y < minLatitude || y > maxLatitude)
+    if(maxLongitude < minLongitude)
     {
-        return false;
+        longResult = longitude < minLongitude || longitude < maxLongitude;
     }
     else
     {
-        return true;
+        longResult = longitude > minLongitude && longitude < maxLongitude;
     }
+
+    return longResult && latResult;
 }

GIS/shaperenderer.h

@@ -21,7 +21,7 @@ class ShapeRenderer : public ShapeReader
     void setTextFieldName(const std::string &name);
 
     void drawShapeMercator(cv::Mat &src, float xStart, float yStart);
-    void drawShapeEquidistant(cv::Mat &src, float xStart, float yStart, float xCenter, float yCenter);
+    void drawShapeEquidistant(cv::Mat &src, float xStart, float yStart, float centerLatitude, float centerLongitude);
 
 public: //setters
     void setThickness(int thickness) {
@@ -35,7 +35,7 @@ class ShapeRenderer : public ShapeReader
     }
 
 private:
-    bool equidistantCheck(float x, float y, float centerLongitude, float centerLatitude);
+    bool equidistantCheck(float latitude, float longitude, float centerLatitude, float centerLongitude);
 
 private:
     cv::Scalar mColor;

common/settings.cpp

@@ -54,6 +54,7 @@ void Settings::parseIni(const std::string &path)
     ini::extract(mIniParser.sections["Program"]["AzimuthalEquidistantProjection"], mEquidistantProjection, true);
     ini::extract(mIniParser.sections["Program"]["MercatorProjection"], mMercatorProjection, true);
     ini::extract(mIniParser.sections["Program"]["SpreadImage"], mSpreadImage, true);
+    ini::extract(mIniParser.sections["Program"]["AddRainOverlay"], mAddRainOverlay, true);
     ini::extract(mIniParser.sections["Program"]["JpgQuality"], mJpegQuality, 90);
     ini::extract(mIniParser.sections["Program"]["AlfaM2"], mAlfaM2, 110.8f);
     ini::extract(mIniParser.sections["Program"]["DeltaM2"], DeltaM2, -3.2f);
@@ -213,7 +214,7 @@ DateTime Settings::getPassDate() const
         if(std::regex_search(dateTimeStr, dateTimeRegex)) {
             std::replace( dateTimeStr.begin(), dateTimeStr.end(), '-', ' ');
             std::istringstream( dateTimeStr ) >> day >> month >> year;
-            dateTime.Initialise(year, month, day, today.tm_hour, today.tm_min, today.tm_sec, 0);
+            dateTime.Initialise(year, month, day, 12, 0, 0, 0);
         } else {
             std::cout << "Invalid given Date format, using today's date" << std::endl;
         }

common/settings.h

@@ -108,6 +108,7 @@ class Settings
     bool equadistantProjection() const { return mEquidistantProjection; }
     bool mercatorProjection() const { return mMercatorProjection; }
     bool spreadImage() const { return mSpreadImage; }
+    bool addRainOverlay() const { return mAddRainOverlay; }
     float getNightPassTreshold() const { return mNightPassTreshold; }
 
     int getCostasBandwidth() const { return mCostasBw; }
@@ -163,6 +164,7 @@ class Settings
     bool mEquidistantProjection;
     bool mMercatorProjection;
     bool mSpreadImage;
+    bool mAddRainOverlay;
     float mNightPassTreshold;
 
     //ini section: Demodulator

common/version.h

@@ -2,7 +2,7 @@
 #define VERSION_H
 
 #define VERSION_MAJOR 2
-#define VERSION_MINOR 0
-#define VERSION_FIX 1
+#define VERSION_MINOR 1
+#define VERSION_FIX 0
 
 #endif // VERSION_H

imageproc/threatimage.cpp

@@ -1,6 +1,7 @@
 #include "threatimage.h"
 #include <iostream>
 #include "settings.h"
+#include <opencv2/imgcodecs.hpp>
 
 std::map<std::string, ThreatImage::WatermarkPosition> ThreatImage::WatermarkPositionLookup {
     {"top_left", WatermarkPosition::TOP_LEFT},
@@ -46,7 +47,7 @@ void ThreatImage::fillBlackLines(cv::Mat &bitmap, int minimumHeight, int maximum
     }
 }
 
-cv::Mat ThreatImage::irToTemperature(const cv::Mat &irImage, const cv::Mat ref)
+cv::Mat ThreatImage::irToTemperature(const cv::Mat &irImage, const cv::Mat &ref)
 {
     if(ref.cols != 256) {
         return cv::Mat();
@@ -63,6 +64,43 @@ cv::Mat ThreatImage::irToTemperature(const cv::Mat &irImage, const cv::Mat ref)
     return thermalImage;
 }
 
+cv::Mat ThreatImage::irToRain(const cv::Mat &irImage, const cv::Mat &ref)
+{
+    if(ref.cols != 256) {
+        return cv::Mat();
+    }
+
+    cv::Mat rainImage = cv::Mat::zeros(irImage.size(), irImage.type());
+
+    for (int x = 0; x < irImage.cols; x++) {
+        for (int y = 0; y < irImage.rows; y++) {
+            uint8_t temp = irImage.at<cv::Vec3b>(y, x)[0];
+            rainImage.at<cv::Vec3b>(y, x) = ref.at<cv::Vec3b>(0,temp);
+        }
+    }
+
+    return rainImage;
+}
+
+cv::Mat ThreatImage::addRainOverlay(const cv::Mat &image, const cv::Mat &rain)
+{
+    cv::Mat rainImage = cv::Mat::zeros(image.size(), image.type());
+    cv::Mat grayScale;
+    cv::Mat alpha;
+
+    //create mask
+    cv::cvtColor(rain, grayScale, cv::COLOR_BGR2GRAY);
+    cv::threshold(grayScale, alpha, 0, 255, cv::THRESH_BINARY);
+
+    //create masked image
+    cv::bitwise_and(image, image, rainImage, cv::Scalar::all(1.0)-alpha);
+
+    //add rain overlay
+    cv::bitwise_and(rain, rain, rainImage, alpha);
+
+    return rainImage;
+}
+
 cv::Mat ThreatImage::gamma(const cv::Mat &image, double gamma)
 {
     cv::Mat newImage = image.clone();

imageproc/threatimage.h

@@ -19,7 +19,9 @@ class ThreatImage
 
 public:
     static void fillBlackLines(cv::Mat &image, int minimumHeight, int maximumHeight);
-    static cv::Mat irToTemperature(const cv::Mat &irImage, const cv::Mat ref);
+    static cv::Mat irToTemperature(const cv::Mat &irImage, const cv::Mat &ref);
+    static cv::Mat irToRain(const cv::Mat &irImage, const cv::Mat &ref);
+    static cv::Mat addRainOverlay(const cv::Mat &image, const cv::Mat &rain);
     static cv::Mat gamma(const cv::Mat &image, double gamma);
     static void drawWatermark(cv::Mat image, const std::string &date);
     static bool isNightPass(const cv::Mat &image, float treshold);

main.cpp

@@ -229,10 +229,16 @@ int main(int argc, char *argv[])
 
     DateTime passStart;
     DateTime passDate = mSettings.getPassDate();
-    TimeSpan passFirstTime = mPacketParser.getFirstTimeStamp();
-    TimeSpan passLength = mPacketParser.getLastTimeStamp() - passFirstTime;
+    TimeSpan passStartTime = mPacketParser.getFirstTimeStamp();
+    TimeSpan passLength = mPacketParser.getLastTimeStamp() - passStartTime;
+
+    passDate = passDate.AddHours(3); //Convert UTC 0 to Moscow time zone (UTC + 3)
+
+    //Satellite's date time
+    passStart.Initialise(passDate.Year(), passDate.Month(), passDate.Day(), passStartTime.Hours(), passStartTime.Minutes(), passStartTime.Seconds(), passStartTime.Microseconds());
+    //Convert satellite's Moscow time zone to UTC 0
+    passStart = passStart.AddHours(-3);
 
-    passStart.Initialise(passDate.Year(), passDate.Month(), passDate.Day(), passFirstTime.Hours()-3, passFirstTime.Minutes(), passFirstTime.Seconds(),passFirstTime.Microseconds());
     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());
 
     PixelGeolocationCalculator calc(tle, passStart, passLength, mSettings.getM2Alfa() / 2.0f, mSettings.getM2Delta());
@@ -246,9 +252,18 @@ int main(int argc, char *argv[])
 	cv::Mat irImage = mPacketParser.getChannelImage(PacketParser::APID_68, mSettings.fillBackLines());
 	cv::Mat threatedImage2 = mPacketParser.getRGBImage(PacketParser::APID_64, PacketParser::APID_65, PacketParser::APID_68, mSettings.fillBackLines());
 
+        cv::Mat rainRef = cv::imread(mSettings.getResourcesPath() + "rain.bmp");
+        cv::Mat rainOverlay = ThreatImage::irToRain(irImage, rainRef);
+
         if(!ThreatImage::isNightPass(threatedImage1, mSettings.getNightPassTreshold())) {
             imagesToSpread.push_back(ImageForSpread(threatedImage1, "221_"));
-	    imagesToSpread.push_back(ImageForSpread(threatedImage2, "125_"));
+            imagesToSpread.push_back(ImageForSpread(threatedImage2, "125_"));
+
+            if(mSettings.addRainOverlay()) {
+                imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(threatedImage1, rainOverlay), "rain_221_"));
+                imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(threatedImage2, rainOverlay), "rain_125_"));
+            }
+
 	    saveImage(mSettings.getOutputPath() + fileNameDate + "_221.bmp", threatedImage1);
 	    saveImage(mSettings.getOutputPath() + fileNameDate + "_125.bmp", threatedImage2);
         } else {
@@ -271,6 +286,10 @@ int main(int argc, char *argv[])
         irImage = ThreatImage::gamma(irImage, 1.8);
         imagesToSpread.push_back(ImageForSpread(irImage, "IR_"));
 
+        if(mSettings.addRainOverlay()) {
+            imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(irImage, rainOverlay), "rain_IR_"));
+        }
+
     } else if(mPacketParser.isChannel64Available() && mPacketParser.isChannel65Available() && mPacketParser.isChannel66Available()) {
         cv::Mat threatedImage = mPacketParser.getRGBImage(PacketParser::APID_66, PacketParser::APID_65, PacketParser::APID_64, mSettings.fillBackLines());
 
@@ -309,8 +328,18 @@ int main(int argc, char *argv[])
         saveImage(mSettings.getOutputPath() + fileNameDate + "_65.bmp", ch65);
     } else if(mPacketParser.isChannel68Available()) {
         cv::Mat ch68 = mPacketParser.getChannelImage(PacketParser::APID_68, 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);
+
+        cv::bitwise_not(ch68, ch68);
+        ch68 = ThreatImage::gamma(ch68, 1.8);
+        imagesToSpread.push_back(ImageForSpread(ch68, "IR_"));
+
+        if(mSettings.addRainOverlay()) {
+            imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(ch68, rainOverlay), "rain_IR_"));
+        }
     } else {
         std::cout << "No usable channel data found!" << std::endl;
 

resources/settings.ini

@@ -2,6 +2,7 @@
 AzimuthalEquidistantProjection=true
 MercatorProjection=true
 SpreadImage=true
+AddRainOverlay=true
 JpgQuality=100
 AlfaM2=110.8
 DeltaM2=-3.2