mirror of
https://github.com/f4exb/sdrangel.git
synced 2024-08-28 08:46:33 -04:00
1015 lines
35 KiB
C++
1015 lines
35 KiB
C++
///////////////////////////////////////////////////////////////////////////////////
|
|
// Copyright (C) 2015-2018 Edouard Griffiths, F4EXB. //
|
|
// Copyright (C) 2021 Jon Beniston, M7RCE //
|
|
// //
|
|
// This program is free software; you can redistribute it and/or modify //
|
|
// it under the terms of the GNU General Public License as published by //
|
|
// the Free Software Foundation as version 3 of the License, or //
|
|
// (at your option) any later version. //
|
|
// //
|
|
// This program is distributed in the hope that it will be useful, //
|
|
// but WITHOUT ANY WARRANTY; without even the implied warranty of //
|
|
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the //
|
|
// GNU General Public License V3 for more details. //
|
|
// //
|
|
// You should have received a copy of the GNU General Public License //
|
|
// along with this program. If not, see <http://www.gnu.org/licenses/>. //
|
|
///////////////////////////////////////////////////////////////////////////////////
|
|
|
|
#include <algorithm>
|
|
|
|
#include <QTime>
|
|
#include <QBuffer>
|
|
#include <QDebug>
|
|
|
|
#include "maincore.h"
|
|
#include "util/units.h"
|
|
|
|
#include "aptdemod.h"
|
|
#include "aptdemodimageworker.h"
|
|
|
|
#include "SWGMapItem.h"
|
|
|
|
MESSAGE_CLASS_DEFINITION(APTDemodImageWorker::MsgConfigureAPTDemodImageWorker, Message)
|
|
MESSAGE_CLASS_DEFINITION(APTDemodImageWorker::MsgSaveImageToDisk, Message)
|
|
MESSAGE_CLASS_DEFINITION(APTDemodImageWorker::MsgSetSatelliteName, Message)
|
|
|
|
APTDemodImageWorker::APTDemodImageWorker(APTDemod *aptDemod) :
|
|
m_messageQueueToGUI(nullptr),
|
|
m_aptDemod(aptDemod),
|
|
m_sgp4(nullptr),
|
|
m_running(false),
|
|
m_mutex(QMutex::Recursive)
|
|
{
|
|
for (int y = 0; y < APT_MAX_HEIGHT; y++)
|
|
{
|
|
m_image.prow[y] = new float[APT_PROW_WIDTH];
|
|
m_tempImage.prow[y] = new float[APT_PROW_WIDTH];
|
|
}
|
|
|
|
resetDecoder();
|
|
}
|
|
|
|
APTDemodImageWorker::~APTDemodImageWorker()
|
|
{
|
|
m_inputMessageQueue.clear();
|
|
|
|
for (int y = 0; y < APT_MAX_HEIGHT; y++)
|
|
{
|
|
delete[] m_image.prow[y];
|
|
delete[] m_tempImage.prow[y];
|
|
}
|
|
|
|
delete m_sgp4;
|
|
}
|
|
|
|
void APTDemodImageWorker::reset()
|
|
{
|
|
QMutexLocker mutexLocker(&m_mutex);
|
|
m_inputMessageQueue.clear();
|
|
}
|
|
|
|
void APTDemodImageWorker::startWork()
|
|
{
|
|
QMutexLocker mutexLocker(&m_mutex);
|
|
connect(&m_inputMessageQueue, SIGNAL(messageEnqueued()), this, SLOT(handleInputMessages()));
|
|
m_running = true;
|
|
}
|
|
|
|
void APTDemodImageWorker::stopWork()
|
|
{
|
|
QMutexLocker mutexLocker(&m_mutex);
|
|
disconnect(&m_inputMessageQueue, SIGNAL(messageEnqueued()), this, SLOT(handleInputMessages()));
|
|
m_running = false;
|
|
}
|
|
|
|
void APTDemodImageWorker::handleInputMessages()
|
|
{
|
|
Message* message;
|
|
|
|
while ((message = m_inputMessageQueue.pop()) != nullptr)
|
|
{
|
|
if (handleMessage(*message)) {
|
|
delete message;
|
|
}
|
|
}
|
|
}
|
|
|
|
bool APTDemodImageWorker::handleMessage(const Message& cmd)
|
|
{
|
|
if (MsgConfigureAPTDemodImageWorker::match(cmd))
|
|
{
|
|
QMutexLocker mutexLocker(&m_mutex);
|
|
MsgConfigureAPTDemodImageWorker& cfg = (MsgConfigureAPTDemodImageWorker&) cmd;
|
|
qDebug("APTDemodImageWorker::handleMessage: MsgConfigureAPTDemodImageWorker");
|
|
applySettings(cfg.getSettings(), cfg.getForce());
|
|
return true;
|
|
}
|
|
else if (MsgSaveImageToDisk::match(cmd))
|
|
{
|
|
saveImageToDisk();
|
|
return true;
|
|
}
|
|
else if (MsgSetSatelliteName::match(cmd))
|
|
{
|
|
MsgSetSatelliteName& msg = (MsgSetSatelliteName&) cmd;
|
|
m_satelliteName = msg.getSatelliteName();
|
|
return true;
|
|
}
|
|
else if (APTDemod::MsgPixels::match(cmd))
|
|
{
|
|
QMutexLocker mutexLocker(&m_mutex);
|
|
const APTDemod::MsgPixels& pixelsMsg = (APTDemod::MsgPixels&) cmd;
|
|
const float *pixels = pixelsMsg.getPixels();
|
|
processPixels(pixels);
|
|
delete[] pixels;
|
|
return true;
|
|
}
|
|
else if (APTDemod::MsgResetDecoder::match(cmd))
|
|
{
|
|
resetDecoder();
|
|
return true;
|
|
}
|
|
else
|
|
{
|
|
return false;
|
|
}
|
|
}
|
|
|
|
void APTDemodImageWorker::applySettings(const APTDemodSettings& settings, bool force)
|
|
{
|
|
(void) force;
|
|
|
|
bool callRecalcCoords = false;
|
|
bool callProcessImage = false;
|
|
|
|
if ((settings.m_cropNoise != m_settings.m_cropNoise) ||
|
|
(settings.m_denoise != m_settings.m_denoise) ||
|
|
(settings.m_linearEqualise != m_settings.m_linearEqualise) ||
|
|
(settings.m_histogramEqualise != m_settings.m_histogramEqualise) ||
|
|
(settings.m_precipitationOverlay != m_settings.m_precipitationOverlay) ||
|
|
(settings.m_flip != m_settings.m_flip) ||
|
|
(settings.m_channels != m_settings.m_channels) ||
|
|
(settings.m_transparencyThreshold != m_settings.m_transparencyThreshold) ||
|
|
(settings.m_opacityThreshold != m_settings.m_opacityThreshold) ||
|
|
(settings.m_palettes != m_settings.m_palettes) ||
|
|
(settings.m_palette != m_settings.m_palette) ||
|
|
(settings.m_horizontalPixelsPerDegree != m_settings.m_horizontalPixelsPerDegree) ||
|
|
(settings.m_verticalPixelsPerDegree != m_settings.m_verticalPixelsPerDegree))
|
|
{
|
|
// Call after settings have been applied
|
|
callProcessImage = true;
|
|
}
|
|
|
|
if ((settings.m_satTimeOffset != m_settings.m_satTimeOffset) ||
|
|
(settings.m_satYaw != m_settings.m_satYaw))
|
|
{
|
|
callRecalcCoords = true;
|
|
callProcessImage = true;
|
|
}
|
|
|
|
if (!settings.m_decodeEnabled && m_settings.m_decodeEnabled)
|
|
{
|
|
// Decode complete - make sure we do a full image update
|
|
// so we aren't left we unprocessed lines
|
|
callProcessImage = true;
|
|
}
|
|
|
|
if (settings.m_palettes != m_settings.m_palettes)
|
|
{
|
|
// Load colour palettes
|
|
m_palettes.clear();
|
|
for (auto palette : settings.m_palettes)
|
|
{
|
|
QImage img;
|
|
img.load(palette);
|
|
if ((img.width() != 256) || (img.height() != 256)) {
|
|
qWarning() << "APT colour palette " << palette << " is not 256x256 pixels - " << img.width() << "x" << img.height();
|
|
}
|
|
m_palettes.append(img);
|
|
}
|
|
}
|
|
|
|
m_settings = settings;
|
|
|
|
if (callRecalcCoords) {
|
|
recalcCoords();
|
|
}
|
|
if (callProcessImage) {
|
|
sendImageToGUI();
|
|
}
|
|
}
|
|
|
|
void APTDemodImageWorker::resetDecoder()
|
|
{
|
|
m_image.nrow = 0;
|
|
m_image.zenith = 0;
|
|
m_tempImage.nrow = 0;
|
|
m_tempImage.zenith = 0;
|
|
m_greyImage = QImage(APT_IMG_WIDTH, APT_MAX_HEIGHT, QImage::Format_Grayscale8);
|
|
m_greyImage.fill(0);
|
|
m_colourImage = QImage(APT_IMG_WIDTH, APT_MAX_HEIGHT, QImage::Format_RGB888);
|
|
m_colourImage.fill(0);
|
|
m_satelliteName = "";
|
|
m_satCoords.clear();
|
|
m_pixelCoords.clear();
|
|
delete m_sgp4;
|
|
m_sgp4 = nullptr;
|
|
}
|
|
|
|
// Convert Qt QDataTime to QGP4 DateTime
|
|
static DateTime qDateTimeToDateTime(QDateTime qdt)
|
|
{
|
|
QDateTime utc = qdt.toUTC();
|
|
QDate date = utc.date();
|
|
QTime time = utc.time();
|
|
DateTime dt;
|
|
dt.Initialise(date.year(), date.month(), date.day(), time.hour(), time.minute(), time.second(), time.msec() * 1000);
|
|
return dt;
|
|
}
|
|
|
|
// Get heading in range [0,360)
|
|
static double normaliseHeading(double heading)
|
|
{
|
|
return fmod(heading + 360.0, 360.0);
|
|
}
|
|
|
|
// Get longitude in range -180,180
|
|
static double normaliseLongitude(double lon)
|
|
{
|
|
return fmod(lon + 540.0, 360.0) - 180.0;
|
|
}
|
|
|
|
// Calculate heading (azimuth) in degrees
|
|
double APTDemodImageWorker::calcHeading(CoordGeodetic from, CoordGeodetic to) const
|
|
{
|
|
// From https://en.wikipedia.org/wiki/Azimuth Section In Geodesy
|
|
double flattening = 1.0 / 298.257223563; // For WGS84 ellipsoid
|
|
double eSq = flattening * (2.0 - flattening);
|
|
double oneMinusESq = (1.0 - flattening) * (1.0 - flattening);
|
|
|
|
double tl1 = tan(from.latitude);
|
|
double tl2 = tan(to.latitude);
|
|
double n1 = 1.0 + oneMinusESq * tl2 * tl2;
|
|
double d1 = 1.0 + oneMinusESq * tl1 * tl1;
|
|
|
|
double l = to.longitude - from.longitude;
|
|
|
|
double alpha;
|
|
|
|
if (from.latitude == 0.0)
|
|
{
|
|
alpha = atan2(sin(l), (oneMinusESq * tan(to.latitude)));
|
|
}
|
|
else
|
|
{
|
|
double lambda = oneMinusESq * tan(to.latitude) / tan(from.latitude) + eSq * sqrt(n1/d1);
|
|
|
|
alpha = atan2(sin(l), ((lambda - cos(l)) * sin(from.latitude)));
|
|
}
|
|
|
|
double deg = Units::radiansToDegrees(alpha);
|
|
if (!m_settings.m_northToSouth) {
|
|
deg += 180.0;
|
|
}
|
|
deg = normaliseHeading(deg);
|
|
return deg;
|
|
}
|
|
|
|
// CoordGeodetic are in radians. Distance in metres. Bearing in radians.
|
|
// https://www.movable-type.co.uk/scripts/latlong.html
|
|
// This approximates Earth as spherical. If we need more accurate algorithm, see:
|
|
// https://www.movable-type.co.uk/scripts/latlong-vincenty.html
|
|
static void calcRadialEndPoint(CoordGeodetic start, double distance, double bearing, CoordGeodetic &end)
|
|
{
|
|
double earthRadius = 6378137.0; // At equator
|
|
double delta = distance/earthRadius;
|
|
end.latitude = std::asin(sin(start.latitude)*cos(delta) + cos(start.latitude)*sin(delta)*cos(bearing));
|
|
end.longitude = start.longitude + std::atan2(sin(bearing)*sin(delta)*cos(start.latitude), cos(delta) - sin(start.latitude)*sin(end.latitude));
|
|
end.longitude = normaliseLongitude(end.longitude);
|
|
}
|
|
|
|
void APTDemodImageWorker::calcPixelCoords(CoordGeodetic centreCoord, double heading)
|
|
{
|
|
// Calculate coordinates of each pixel in a row (swath)
|
|
// Assume satellite is at centre pixel, and project +-90 degrees from satellite heading
|
|
// https://www.ncei.noaa.gov/pub/data/satellite/publications/podguides/N-15%20thru%20N-19/pdf/APPENDIX%20J%20Instrument%20Scan%20Properties.pdf
|
|
// Swath for AVHRR/3 of 2926.6km at 833km altitude over spherical Earth
|
|
// Some docs say resolution is 4.0km, but it varies as per fig 4.2.3-1 in:
|
|
// https://www.ncei.noaa.gov/pub/data/satellite/publications/podguides/N-15%20thru%20N-19/pdf/2.1%20Section%204.0%20Real%20Time%20Data%20Systems%20for%20Local%20Users%20.pdf
|
|
// TODO: Could try to adjust for altitude
|
|
|
|
QVector<CoordGeodetic> pixelCoords(APT_CH_WIDTH);
|
|
pixelCoords[APT_CH_WIDTH/2] = centreCoord;
|
|
double heading1 = Units::degreesToRadians(heading + m_settings.m_satYaw + 90.0);
|
|
double heading2 = Units::degreesToRadians(heading + m_settings.m_satYaw - 90.0);
|
|
for (int i = 1; i <= APT_CH_WIDTH/2; i++)
|
|
{
|
|
double distance = i * 2926600.0/APT_CH_WIDTH;
|
|
calcRadialEndPoint(centreCoord, distance, heading1, pixelCoords[APT_CH_WIDTH/2-i]);
|
|
calcRadialEndPoint(centreCoord, distance, heading2, pixelCoords[APT_CH_WIDTH/2+i]);
|
|
}
|
|
|
|
if (m_settings.m_northToSouth) {
|
|
m_pixelCoords.append(pixelCoords);
|
|
} else {
|
|
m_pixelCoords.prepend(pixelCoords);
|
|
}
|
|
}
|
|
|
|
// Recalculate all pixel coordiantes as satTimeOffset or satYaw has changed
|
|
void APTDemodImageWorker::recalcCoords()
|
|
{
|
|
if (m_sgp4)
|
|
{
|
|
m_satCoords.clear();
|
|
m_pixelCoords.clear();
|
|
for (int row = 0; row < m_image.nrow; row++)
|
|
{
|
|
QDateTime qdt = m_settings.m_aosDateTime.addMSecs(m_settings.m_satTimeOffset * 1000.0f + row * 500);
|
|
calcCoords(qdt, row);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Calculate pixel coordinates for a single row at the given date and time
|
|
void APTDemodImageWorker::calcCoords(QDateTime qdt, int row)
|
|
{
|
|
try
|
|
{
|
|
DateTime dt = qDateTimeToDateTime(qdt);
|
|
|
|
// Calculate satellite position
|
|
Eci eci = m_sgp4->FindPosition(dt);
|
|
|
|
// Convert satellite position to geodetic coordinates (lat and long)
|
|
CoordGeodetic geo = eci.ToGeodetic();
|
|
|
|
m_satCoords.append(geo);
|
|
|
|
// Calculate satellite heading (Could convert eci.Velocity() instead)
|
|
double heading;
|
|
if (m_satCoords.size() == 2)
|
|
{
|
|
heading = calcHeading(m_satCoords[0], m_satCoords[1]);
|
|
calcPixelCoords(m_satCoords[0], heading);
|
|
calcPixelCoords(m_satCoords[1], heading);
|
|
}
|
|
else if (m_satCoords.size() > 2)
|
|
{
|
|
heading = calcHeading(m_satCoords[row-1], m_satCoords[row]);
|
|
calcPixelCoords(geo, heading);
|
|
}
|
|
}
|
|
catch (SatelliteException& se)
|
|
{
|
|
qDebug() << "APTDemodImageWorker::calcCoord: " << se.what();
|
|
}
|
|
catch (DecayedException& de)
|
|
{
|
|
qDebug() << "APTDemodImageWorker::calcCoord: " << de.what();
|
|
}
|
|
catch (TleException& tlee)
|
|
{
|
|
qDebug() << "APTDemodImageWorker::calcCoord: " << tlee.what();
|
|
}
|
|
}
|
|
|
|
|
|
// Calculate satellite's geodetic coordinates and heading
|
|
void APTDemodImageWorker::calcCoord(int row)
|
|
{
|
|
if (row == 0)
|
|
{
|
|
QStringList elements = m_settings.m_tle.trimmed().split("\n");
|
|
if (elements.size() == 3)
|
|
{
|
|
// Initalise SGP4
|
|
Tle tle(elements[0].toStdString(), elements[1].toStdString(), elements[2].toStdString());
|
|
m_sgp4 = new SGP4(tle);
|
|
|
|
// Output time so we can check time offset from when AOS is signalled
|
|
qDebug() << "APTDemod: Processing row 0 at " << QDateTime::currentDateTime();
|
|
|
|
calcCoords(m_settings.m_aosDateTime, row);
|
|
}
|
|
else
|
|
{
|
|
qDebug() << "APTDemodImageWorker::calcCoord: No TLE for satellite. Is Satellite Tracker running?";
|
|
return;
|
|
}
|
|
}
|
|
else if (m_sgp4 == nullptr)
|
|
{
|
|
return;
|
|
}
|
|
else
|
|
{
|
|
// Calculate time at which
|
|
// Don't try to use QDateTime::currentDateTime() as processing & scheduling delays mean
|
|
// it's not constant and can sometimes even be 0
|
|
// Lines should be transmitted at 2 per second, so just use number of rows since AOS
|
|
// We add a user-defined delay to account for delays in transferring SDR data and demodulation
|
|
QDateTime qdt = m_settings.m_aosDateTime.addMSecs(m_settings.m_satTimeOffset * 1000.0f + row * 500);
|
|
calcCoords(qdt, row);
|
|
}
|
|
}
|
|
|
|
void APTDemodImageWorker::processPixels(const float *pixels)
|
|
{
|
|
if (m_image.nrow < APT_MAX_HEIGHT)
|
|
{
|
|
// Calculate lat and lon of centre of row
|
|
calcCoord(m_image.nrow);
|
|
|
|
std::copy(pixels, pixels + APT_PROW_WIDTH, m_image.prow[m_image.nrow]);
|
|
m_image.nrow++;
|
|
|
|
if (m_image.nrow % m_settings.m_scanlinesPerImageUpdate == 0) { // send full image only every N lines
|
|
sendImageToGUI();
|
|
} else { // else send unprocessed line just to show stg is moving
|
|
sendLineToGUI();
|
|
}
|
|
}
|
|
}
|
|
|
|
void APTDemodImageWorker::sendImageToGUI()
|
|
{
|
|
// Send image to GUI
|
|
if (m_messageQueueToGUI)
|
|
{
|
|
QStringList imageTypes;
|
|
QImage image = processImage(imageTypes, m_settings.m_channels);
|
|
m_messageQueueToGUI->push(APTDemod::MsgImage::create(image, imageTypes, m_satelliteName));
|
|
if (m_sgp4) {
|
|
sendImageToMap(image);
|
|
}
|
|
}
|
|
}
|
|
|
|
// Find the value of the pixel closest to the given coordinates
|
|
// If we have previously found a pixel, we constrain the search to be nearby, in order to speed up the search
|
|
QRgb APTDemodImageWorker::findNearest(const QImage &image, double latitude, double longitude, int xPrevious, int yPrevious, int &xNearest, int &yNearest) const
|
|
{
|
|
double dmin = 360.0 * 360.0 + 90.0 * 90.0;
|
|
xNearest = -1;
|
|
yNearest = -1;
|
|
QRgb p = qRgba(0, 0, 0, 0); // Transparent
|
|
|
|
int xMin, xMax;
|
|
int yMin, yMax;
|
|
|
|
int yStartPostCrop;
|
|
int yEndPostCrop;
|
|
if (m_settings.m_northToSouth)
|
|
{
|
|
yStartPostCrop = abs(m_tempImage.zenith);
|
|
yEndPostCrop = yStartPostCrop + image.height();
|
|
}
|
|
else
|
|
{
|
|
yStartPostCrop = m_image.nrow - m_tempImage.nrow - abs(m_tempImage.zenith);
|
|
yEndPostCrop = yStartPostCrop + image.height();
|
|
}
|
|
|
|
if (xPrevious == -1)
|
|
{
|
|
yMin = yStartPostCrop;
|
|
yMax = yEndPostCrop;
|
|
xMin = 0;
|
|
xMax = m_pixelCoords[0].size();
|
|
}
|
|
else
|
|
{
|
|
int searchRadius = 4;
|
|
yMin = yPrevious - searchRadius;
|
|
yMax = yPrevious + searchRadius + 1;
|
|
xMin = xPrevious - searchRadius;
|
|
xMax = xPrevious + searchRadius + 1;
|
|
yMin = std::max(yMin, yStartPostCrop);
|
|
yMax = std::min(yMax, yEndPostCrop);
|
|
xMin = std::max(xMin, 0);
|
|
xMax = std::min(xMax, m_pixelCoords[0].size());
|
|
}
|
|
|
|
const int ySize = yEndPostCrop-1;
|
|
const int xSize = m_pixelCoords[0].size()-1;
|
|
for (int y = yMin; y < yMax; y++)
|
|
{
|
|
for (int x = xMin; x < xMax; x++)
|
|
{
|
|
CoordGeodetic coord = m_pixelCoords[y][x];
|
|
double dlat = coord.latitude - latitude;
|
|
double dlon = coord.longitude - longitude;
|
|
double d = dlat * dlat + dlon * dlon;
|
|
if (d < dmin)
|
|
{
|
|
dmin = d;
|
|
xNearest = x;
|
|
yNearest = y;
|
|
// Only use color of pixel if we're inside the source image
|
|
if ( ((y != yStartPostCrop) || ((y == yStartPostCrop) && (latitude <= coord.latitude)))
|
|
&& ((y != ySize) || ((y == ySize) && (latitude >= coord.latitude)))
|
|
&& ((x != 0) || ((x == 0) && (longitude >= coord.longitude)))
|
|
&& ((x != xSize) || ((x == xSize) && (longitude <= coord.longitude)))
|
|
)
|
|
{
|
|
p = image.pixel(x, y - yStartPostCrop);
|
|
}
|
|
else
|
|
{
|
|
p = qRgba(0, 0, 0, 0); // Transparent
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
return p;
|
|
}
|
|
|
|
// Calculate bounding box for projected image in terms of latitude and longitude
|
|
// TODO: Handle crossing of anti-meridian
|
|
void APTDemodImageWorker::calcBoundingBox(double &east, double &south, double &west, double &north, const QImage &image)
|
|
{
|
|
int start;
|
|
if (m_settings.m_northToSouth) {
|
|
start = abs(m_tempImage.zenith);
|
|
} else {
|
|
start = m_image.nrow - m_tempImage.nrow - abs(m_tempImage.zenith);
|
|
}
|
|
int stop = start + image.height();
|
|
|
|
east = -M_PI;
|
|
west = M_PI;
|
|
north = -M_PI/2.0;
|
|
south = M_PI/2.0;
|
|
|
|
//FILE *f = fopen("coords.txt", "w");
|
|
for (int y = start; y < stop; y++)
|
|
{
|
|
for (int x = 0; x < m_pixelCoords[y].size(); x++)
|
|
{
|
|
double latitude = m_pixelCoords[y][x].latitude;
|
|
double longitude = m_pixelCoords[y][x].longitude;
|
|
//fprintf(f, "%f,%f ", Units::radiansToDegrees(m_pixelCoords[y][x].latitude), Units::radiansToDegrees(m_pixelCoords[y][x].longitude));
|
|
south = std::min(latitude, south);
|
|
north = std::max(latitude, north);
|
|
east = std::max(longitude, east);
|
|
west = std::min(longitude, west);
|
|
}
|
|
//fprintf(f, "\n");
|
|
}
|
|
//fclose(f);
|
|
}
|
|
|
|
// Project satellite image to equidistant cyclindrical projection (Plate Carree) for use on 3D Map
|
|
// We've previously computed lat and lon for each pixel in satellite image
|
|
// so we just work through coords in projected image, trying to find closest pixel in satellite image
|
|
// FIXME: How do we handle sat going over the poles?
|
|
QImage APTDemodImageWorker::projectImage(const QImage &image)
|
|
{
|
|
double east, south, west, north;
|
|
|
|
// Calculate bounding box for image tile
|
|
calcBoundingBox(east, south, west, north, image);
|
|
m_tileEast = ceil(Units::radiansToDegrees(east));
|
|
m_tileWest = floor(Units::radiansToDegrees(west));
|
|
m_tileNorth = ceil(Units::radiansToDegrees(north));
|
|
m_tileSouth = floor(Units::radiansToDegrees(south));
|
|
|
|
double widthDeg = m_tileEast - m_tileWest;
|
|
double heightDeg = m_tileNorth - m_tileSouth;
|
|
|
|
int width = widthDeg * m_settings.m_horizontalPixelsPerDegree;
|
|
int height = heightDeg * m_settings.m_verticalPixelsPerDegree;
|
|
|
|
//image.save("source.png");
|
|
//FILE *f = fopen("mapping.txt", "w");
|
|
QImage projection(width, height, QImage::Format_ARGB32);
|
|
int xNearest, yNearest, xPrevious, yPrevious;
|
|
xPrevious = -1;
|
|
yPrevious = -1;
|
|
for (int y = 0; y < height; y++)
|
|
{
|
|
// Calculate geodetic coords of pixel in projected image
|
|
double lat = m_tileNorth - (y / (double)m_settings.m_verticalPixelsPerDegree);
|
|
// Reverse search direction in alternate rows, so we are always seaching
|
|
// close to previously found pixel
|
|
if ((y & 1) == 0)
|
|
{
|
|
for (int x = 0; x < width; x++)
|
|
{
|
|
double lon = m_tileWest + (x / (double)m_settings.m_horizontalPixelsPerDegree);
|
|
// Find closest pixel in source image
|
|
QRgb pixel = findNearest(image, Units::degreesToRadians(lat), Units::degreesToRadians(lon), xPrevious, yPrevious, xNearest, yNearest);
|
|
xPrevious = xNearest;
|
|
yPrevious = yNearest;
|
|
projection.setPixel(x, y, pixel);
|
|
//fprintf(f, "%f,%f,%d,%d,%d ", lat, lon, xNearest, yNearest, pixel==0);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
for (int x = width - 1; x >= 0; x--)
|
|
{
|
|
double lon = m_tileWest + (x / (double)m_settings.m_horizontalPixelsPerDegree);
|
|
// Find closest pixel in source image
|
|
QRgb pixel = findNearest(image, Units::degreesToRadians(lat), Units::degreesToRadians(lon), xPrevious, yPrevious, xNearest, yNearest);
|
|
xPrevious = xNearest;
|
|
yPrevious = yNearest;
|
|
projection.setPixel(x, y, pixel);
|
|
//fprintf(f, "%f,%f,%d,%d,%d ", lat, lon, xNearest, yNearest, pixel==0);
|
|
}
|
|
}
|
|
//fprintf(f, "\n");
|
|
}
|
|
//fclose(f);
|
|
|
|
return projection;
|
|
}
|
|
|
|
// Make an image transparent, so when overlaid on 3D map, we can see the underlying terrain
|
|
// Image is full transparent below m_transparencyThreshold and fully opaque above m_opacityThreshold
|
|
void APTDemodImageWorker::makeTransparent(QImage &image)
|
|
{
|
|
for (int y = 0; y < image.height(); y++)
|
|
{
|
|
for (int x = 0; x < image.width(); x++)
|
|
{
|
|
QRgb pixel = image.pixel(x, y);
|
|
int grey = qGray(pixel);
|
|
if (grey < m_settings.m_transparencyThreshold)
|
|
{
|
|
// Make fully transparent
|
|
pixel = qRgba(qRed(pixel), qGreen(pixel), qBlue(pixel), 0);
|
|
image.setPixel(x, y, pixel);
|
|
}
|
|
else if (grey < m_settings.m_opacityThreshold)
|
|
{
|
|
// Make slightly transparent
|
|
float opacity = 1.0f - ((m_settings.m_opacityThreshold - grey) / (float)(m_settings.m_opacityThreshold - m_settings.m_transparencyThreshold));
|
|
opacity = opacity * 255.0f;
|
|
opacity = std::min(255.0f, opacity);
|
|
opacity = std::max(0.0f, opacity);
|
|
pixel = qRgba(qRed(pixel), qGreen(pixel), qBlue(pixel), (int)std::round(opacity));
|
|
image.setPixel(x, y, pixel);
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void APTDemodImageWorker::sendImageToMap(QImage image)
|
|
{
|
|
// Send to Map feature
|
|
MessagePipesLegacy& messagePipes = MainCore::instance()->getMessagePipes();
|
|
QList<MessageQueue*> *mapMessageQueues = messagePipes.getMessageQueues(m_aptDemod, "mapitems");
|
|
if (mapMessageQueues)
|
|
{
|
|
// Only display one channel on map
|
|
QImage selectedChannel;
|
|
if (m_settings.m_channels == APTDemodSettings::BOTH_CHANNELS) {
|
|
selectedChannel = extractImage(image, APTDemodSettings::CHANNEL_B);
|
|
} else {
|
|
selectedChannel = image;
|
|
}
|
|
|
|
// Project image to geodetic coords (lat & lon)
|
|
selectedChannel = projectImage(selectedChannel);
|
|
//selectedChannel.save("projected.png");
|
|
|
|
// Use alpha channel to remove land & sea
|
|
makeTransparent(selectedChannel);
|
|
|
|
// Encode image as base64 PNG
|
|
QByteArray ba;
|
|
QBuffer buffer(&ba);
|
|
buffer.open(QIODevice::WriteOnly);
|
|
selectedChannel.save(&buffer, "PNG");
|
|
QByteArray data = ba.toBase64();
|
|
|
|
// Create name for the image
|
|
QString satName = m_satelliteName;
|
|
satName.replace(" ", "_");
|
|
QString name = QString("apt_%1_%2").arg(satName).arg(m_settings.m_aosDateTime.toString("yyyyMMdd_hhmmss"));
|
|
|
|
// Send name to GUI
|
|
m_messageQueueToGUI->push(APTDemod::MsgMapImageName::create(name));
|
|
|
|
QList<MessageQueue*>::iterator it = mapMessageQueues->begin();
|
|
for (; it != mapMessageQueues->end(); ++it)
|
|
{
|
|
SWGSDRangel::SWGMapItem *swgMapItem = new SWGSDRangel::SWGMapItem();
|
|
swgMapItem->setName(new QString(name));
|
|
swgMapItem->setImage(new QString(data));
|
|
swgMapItem->setAltitude(3000.0); // Typical cloud height - So it appears above objects on the ground
|
|
swgMapItem->setType(1);
|
|
swgMapItem->setImageTileEast(m_tileEast);
|
|
swgMapItem->setImageTileWest(m_tileWest);
|
|
swgMapItem->setImageTileNorth(m_tileNorth);
|
|
swgMapItem->setImageTileSouth(m_tileSouth);
|
|
|
|
MainCore::MsgMapItem *msg = MainCore::MsgMapItem::create(m_aptDemod, swgMapItem);
|
|
(*it)->push(msg);
|
|
}
|
|
}
|
|
}
|
|
|
|
void APTDemodImageWorker::sendLineToGUI()
|
|
{
|
|
if (m_messageQueueToGUI)
|
|
{
|
|
float *pixels = m_image.prow[m_image.nrow-1];
|
|
uchar *line;
|
|
APTDemod::MsgLine *msg = APTDemod::MsgLine::create(&line);
|
|
|
|
if (m_settings.m_channels == APTDemodSettings::BOTH_CHANNELS)
|
|
{
|
|
for (int i = 0; i < APT_IMG_WIDTH; i++) {
|
|
line[i] = roundAndClip(pixels[i]);
|
|
}
|
|
msg->setSize(APT_IMG_WIDTH);
|
|
}
|
|
else if (m_settings.m_channels == APTDemodSettings::CHANNEL_A)
|
|
{
|
|
for (int i = 0; i < APT_CH_WIDTH; i++) {
|
|
line[i] = roundAndClip(pixels[i + APT_CHA_OFFSET]);
|
|
}
|
|
msg->setSize(APT_CH_WIDTH);
|
|
}
|
|
else
|
|
{
|
|
for (int i = 0; i < APT_CH_WIDTH; i++) {
|
|
line[i] = roundAndClip(pixels[i + APT_CHB_OFFSET]);
|
|
}
|
|
msg->setSize(APT_CH_WIDTH);
|
|
}
|
|
|
|
m_messageQueueToGUI->push(msg);
|
|
}
|
|
}
|
|
|
|
QImage APTDemodImageWorker::processImage(QStringList& imageTypes, APTDemodSettings::ChannelSelection channels)
|
|
{
|
|
copyImage(&m_tempImage, &m_image);
|
|
|
|
// Calibrate channels according to wavelength (1.7x to stop flickering)
|
|
if (m_tempImage.nrow >= 1.7 * APT_CALIBRATION_ROWS)
|
|
{
|
|
m_tempImage.chA = apt_calibrate(m_tempImage.prow, m_tempImage.nrow, APT_CHA_OFFSET, APT_CH_WIDTH);
|
|
m_tempImage.chB = apt_calibrate(m_tempImage.prow, m_tempImage.nrow, APT_CHB_OFFSET, APT_CH_WIDTH);
|
|
QStringList channelTypes({
|
|
"", // Unknown
|
|
"Visible (0.58-0.68 um)",
|
|
"Near-IR (0.725-1.0 um)",
|
|
"Near-IR (1.58-1.64 um)",
|
|
"Thermal-infrared (10.3-11.3 um)",
|
|
"Thermal-infrared (11.5-12.5 um)",
|
|
"Mid-infrared (3.55-3.93 um)"
|
|
});
|
|
|
|
imageTypes.append(channelTypes[m_tempImage.chA]);
|
|
imageTypes.append(channelTypes[m_tempImage.chB]);
|
|
}
|
|
|
|
// Crop noise due to low elevation at top and bottom of image
|
|
if (m_settings.m_cropNoise) {
|
|
m_tempImage.zenith -= apt_cropNoise(&m_tempImage);
|
|
}
|
|
|
|
// Denoise filter
|
|
if (m_settings.m_denoise)
|
|
{
|
|
apt_denoise(m_tempImage.prow, m_tempImage.nrow, APT_CHA_OFFSET, APT_CH_WIDTH);
|
|
apt_denoise(m_tempImage.prow, m_tempImage.nrow, APT_CHB_OFFSET, APT_CH_WIDTH);
|
|
}
|
|
|
|
// Flip image if satellite pass is North to South
|
|
if (m_settings.m_flip)
|
|
{
|
|
apt_flipImage(&m_tempImage, APT_CH_WIDTH, APT_CHA_OFFSET);
|
|
apt_flipImage(&m_tempImage, APT_CH_WIDTH, APT_CHB_OFFSET);
|
|
}
|
|
|
|
// Linear equalise to improve contrast
|
|
if (m_settings.m_linearEqualise)
|
|
{
|
|
apt_linearEnhance(m_tempImage.prow, m_tempImage.nrow, APT_CHA_OFFSET, APT_CH_WIDTH);
|
|
apt_linearEnhance(m_tempImage.prow, m_tempImage.nrow, APT_CHB_OFFSET, APT_CH_WIDTH);
|
|
}
|
|
|
|
// Histogram equalise to improve contrast
|
|
if (m_settings.m_histogramEqualise)
|
|
{
|
|
apt_histogramEqualise(m_tempImage.prow, m_tempImage.nrow, APT_CHA_OFFSET, APT_CH_WIDTH);
|
|
apt_histogramEqualise(m_tempImage.prow, m_tempImage.nrow, APT_CHB_OFFSET, APT_CH_WIDTH);
|
|
}
|
|
|
|
if (m_settings.m_precipitationOverlay)
|
|
{
|
|
// Overlay precipitation
|
|
for (int r = 0; r < m_tempImage.nrow; r++)
|
|
{
|
|
uchar *l = m_colourImage.scanLine(r);
|
|
for (int i = 0; i < APT_IMG_WIDTH; i++)
|
|
{
|
|
float p = m_tempImage.prow[r][i];
|
|
|
|
if ((i >= APT_CHB_OFFSET) && (i < APT_CHB_OFFSET + APT_CH_WIDTH) && (p >= 198))
|
|
{
|
|
apt_rgb_t rgb = apt_applyPalette(apt_PrecipPalette, p - 198);
|
|
// Negative float values get converted to positive uchars here
|
|
l[i*3] = (uchar)rgb.r;
|
|
l[i*3+1] = (uchar)rgb.g;
|
|
l[i*3+2] = (uchar)rgb.b;
|
|
int a = i - APT_CHB_OFFSET + APT_CHA_OFFSET;
|
|
l[a*3] = (uchar)rgb.r;
|
|
l[a*3+1] = (uchar)rgb.g;
|
|
l[a*3+2] = (uchar)rgb.b;
|
|
}
|
|
else
|
|
{
|
|
uchar q = roundAndClip(p);
|
|
l[i*3] = q;
|
|
l[i*3+1] = q;
|
|
l[i*3+2] = q;
|
|
}
|
|
}
|
|
}
|
|
return extractImage(m_colourImage, channels);
|
|
}
|
|
else if (channels == APTDemodSettings::TEMPERATURE)
|
|
{
|
|
// Temperature calibration
|
|
int satnum = 15;
|
|
if (m_satelliteName == "NOAA 18") {
|
|
satnum = 18;
|
|
} else if (m_satelliteName == "NOAA 19") {
|
|
satnum = 19;
|
|
}
|
|
apt_temperature(satnum, &m_tempImage, APT_CHB_OFFSET, APT_CH_WIDTH);
|
|
|
|
// Apply colour palette
|
|
for (int r = 0; r < m_tempImage.nrow; r++)
|
|
{
|
|
uchar *l = m_colourImage.scanLine(r);
|
|
for (int i = 0; i < APT_CH_WIDTH; i++)
|
|
{
|
|
float p = m_tempImage.prow[r][i+APT_CHB_OFFSET];
|
|
uchar q = roundAndClip(p);
|
|
l[i*3] = apt_TempPalette[q*3];
|
|
l[i*3+1] = apt_TempPalette[q*3+1];
|
|
l[i*3+2] = apt_TempPalette[q*3+2];
|
|
}
|
|
}
|
|
return m_colourImage.copy(0, 0, APT_CH_WIDTH, m_tempImage.nrow);
|
|
}
|
|
else if (channels == APTDemodSettings::PALETTE)
|
|
{
|
|
if ((m_settings.m_palette >= 0) && (m_settings.m_palette < m_palettes.size()))
|
|
{
|
|
// Apply colour palette
|
|
for (int r = 0; r < m_tempImage.nrow; r++)
|
|
{
|
|
uchar *l = m_colourImage.scanLine(r);
|
|
for (int i = 0; i < APT_CH_WIDTH; i++)
|
|
{
|
|
float pA = m_tempImage.prow[r][i+APT_CHA_OFFSET];
|
|
float pB = m_tempImage.prow[r][i+APT_CHB_OFFSET];
|
|
uchar qA = roundAndClip(pA);
|
|
uchar qB = roundAndClip(pB);
|
|
QRgb rgb = m_palettes[m_settings.m_palette].pixel(qA, qB);
|
|
l[i*3] = qRed(rgb);
|
|
l[i*3+1] = qGreen(rgb);
|
|
l[i*3+2] = qBlue(rgb);
|
|
}
|
|
}
|
|
return m_colourImage.copy(0, 0, APT_CH_WIDTH, m_tempImage.nrow);
|
|
}
|
|
else
|
|
{
|
|
qDebug() << "APTDemodImageWorker::processImage - Illegal palette number: " << m_settings.m_palette;
|
|
return QImage();
|
|
}
|
|
}
|
|
else
|
|
{
|
|
// Extract grey-scale image
|
|
for (int r = 0; r < m_tempImage.nrow; r++)
|
|
{
|
|
uchar *l = m_greyImage.scanLine(r);
|
|
|
|
for (int i = 0; i < APT_IMG_WIDTH; i++)
|
|
{
|
|
float p = m_tempImage.prow[r][i];
|
|
l[i] = roundAndClip(p);
|
|
}
|
|
}
|
|
return extractImage(m_greyImage, channels);
|
|
}
|
|
}
|
|
|
|
QImage APTDemodImageWorker::extractImage(QImage image, APTDemodSettings::ChannelSelection channels)
|
|
{
|
|
if (channels == APTDemodSettings::BOTH_CHANNELS) {
|
|
return image.copy(0, 0, APT_IMG_WIDTH, m_tempImage.nrow);
|
|
} else if (channels == APTDemodSettings::CHANNEL_A) {
|
|
return image.copy(APT_CHA_OFFSET, 0, APT_CH_WIDTH, m_tempImage.nrow);
|
|
} else {
|
|
return image.copy(APT_CHB_OFFSET, 0, APT_CH_WIDTH, m_tempImage.nrow);
|
|
}
|
|
}
|
|
|
|
void APTDemodImageWorker::prependPath(QString &filename)
|
|
{
|
|
if (!m_settings.m_autoSavePath.isEmpty())
|
|
{
|
|
if (m_settings.m_autoSavePath.endsWith('/')) {
|
|
filename = m_settings.m_autoSavePath + filename;
|
|
} else {
|
|
filename = m_settings.m_autoSavePath + '/' + filename;
|
|
}
|
|
}
|
|
}
|
|
|
|
void APTDemodImageWorker::saveImageToDisk()
|
|
{
|
|
QStringList imageTypes;
|
|
QImage image = processImage(imageTypes, APTDemodSettings::BOTH_CHANNELS);
|
|
|
|
if (image.height() >= m_settings.m_autoSaveMinScanLines)
|
|
{
|
|
QString filename;
|
|
QDateTime dateTime;
|
|
QString dt;
|
|
if (m_settings.m_aosDateTime.isValid()) {
|
|
dateTime = m_settings.m_aosDateTime;
|
|
} else {
|
|
dateTime = QDateTime::currentDateTime();
|
|
}
|
|
dt = dateTime.toString("yyyyMMdd_hhmm");
|
|
QString sat = m_satelliteName;
|
|
sat.replace(" ", "_");
|
|
|
|
if (m_settings.m_saveCombined)
|
|
{
|
|
filename = QString("apt_%1_%2.png").arg(sat).arg(dt);
|
|
prependPath(filename);
|
|
if (!image.save(filename)) {
|
|
qCritical() << "Failed to save APT image to: " << filename;
|
|
}
|
|
}
|
|
|
|
QImage chA = extractImage(image, APTDemodSettings::CHANNEL_A);
|
|
QImage chB = extractImage(image, APTDemodSettings::CHANNEL_B);
|
|
|
|
if (m_settings.m_saveSeparate)
|
|
{
|
|
filename = QString("apt_%1_%2_cha.png").arg(sat).arg(dt);
|
|
prependPath(filename);
|
|
if (!chA.save(filename)) {
|
|
qCritical() << "Failed to save APT image to: " << filename;
|
|
}
|
|
filename = QString("apt_%1_%2_chb.png").arg(sat).arg(dt);
|
|
prependPath(filename);
|
|
if (!chB.save(filename)) {
|
|
qCritical() << "Failed to save APT image to: " << filename;
|
|
}
|
|
}
|
|
|
|
if (m_settings.m_saveProjection)
|
|
{
|
|
filename = QString("apt_%1_%2_cha_eqi_cylindrical.png").arg(sat).arg(dt);
|
|
prependPath(filename);
|
|
QImage chAProj = projectImage(chA);
|
|
if (!chAProj.save(filename)) {
|
|
qCritical() << "Failed to save APT image to: " << filename;
|
|
}
|
|
filename = QString("apt_%1_%2_chb_eqi_cylindrical.png").arg(sat).arg(dt);
|
|
prependPath(filename);
|
|
QImage chBProj = projectImage(chB);
|
|
if (!chBProj.save(filename)) {
|
|
qCritical() << "Failed to save APT image to: " << filename;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
void APTDemodImageWorker::copyImage(apt_image_t *dst, apt_image_t *src)
|
|
{
|
|
dst->nrow = src->nrow;
|
|
dst->zenith = src->zenith;
|
|
dst->chA = src->chA;
|
|
dst->chB = src->chB;
|
|
|
|
for (int i = 0; i < src->nrow; i++) {
|
|
std::copy(src->prow[i], src->prow[i] + APT_PROW_WIDTH, dst->prow[i]);
|
|
}
|
|
}
|
|
|
|
uchar APTDemodImageWorker::roundAndClip(float p)
|
|
{
|
|
int q = (int) round(p);
|
|
q = q > 255 ? 255 : q < 0 ? 0 : q;
|
|
return q;
|
|
}
|