diff --git a/CMakeLists.txt b/CMakeLists.txt index 592a3e3..c2ebe22 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -3,7 +3,7 @@ cmake_minimum_required(VERSION 3.5) include(ExternalProject) project(meteordemod - VERSION 2.6.1 + VERSION 2.6.2 LANGUAGES CXX ) @@ -52,6 +52,7 @@ ExternalProject_Add(libcorrect ) find_package(OpenCV) +find_package(OpenCL) link_directories( ${CMAKE_BINARY_DIR}/sgp4-build/libsgp4 @@ -62,8 +63,10 @@ add_definitions(-D_USE_MATH_DEFINES -D_SILENCE_EXPERIMENTAL_FILESYSTEM_DEPRECATI add_executable(meteordemod main.cpp - imageproc/spreadimage.cpp imageproc/threatimage.cpp + imageproc/projectimage.cpp + imageproc/blendimages.cpp + imageproc/tps.cpp decoder/correlation.cpp decoder/reedsolomon.cpp decoder/viterbi.cpp @@ -83,6 +86,7 @@ add_executable(meteordemod tools/databuffer.cpp tools/iniparser.cpp tools/threadpool.cpp + tools/opencl.cpp GIS/shapereader.cpp GIS/shaperenderer.cpp GIS/dbfilereader.cpp @@ -119,6 +123,12 @@ target_include_directories(meteordemod PUBLIC 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} 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/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/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..3d98be0 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.2" #define MyAppPublisher "Digitelektro" #define MyAppURL "https://github.com/Digitelektro/MeteorDemod" #define MyAppExeName "meteordemod.exe" diff --git a/main.cpp b/main.cpp index 96a6b48..f4f3fa2 100644 --- a/main.cpp +++ b/main.cpp @@ -13,8 +13,10 @@ #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" @@ -32,7 +34,18 @@ struct ImageForSpread { std::string fileNameBase; }; -void searchForImages(std::list& imagesOut, std::list& geolocationCalculatorsOut, const std::string& channelName); +struct ImageSearchResult { + 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); @@ -148,7 +161,7 @@ int main(int argc, char* argv[]) { caduFileStream.close(); } } - } catch(std::exception ex) { + } catch(const std::exception& ex) { std::cout << ex.what() << std::endl; } @@ -175,6 +188,13 @@ 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 << std::to_string(passStart.Ticks()) << std::endl; + datFileStream << std::to_string(passLength.Ticks()) << std::endl; + datFileStream.close(); + } + std::list imagesToSpread; if(mLrptDecoder.isChannel64Available() && mLrptDecoder.isChannel65Available() && mLrptDecoder.isChannel68Available()) { @@ -219,10 +239,10 @@ 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()) { @@ -244,10 +264,10 @@ int main(int argc, char* argv[]) { imagesToSpread.push_back(ImageForSpread(thermalImage, "thermal_")); irImage = ThreatImage::invertIR(irImage); - imagesToSpread.push_back(ImageForSpread(irImage, "IR_")); + imagesToSpread.push_back(ImageForSpread(irImage, "67_")); if(mSettings.addRainOverlay()) { - imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(irImage, rainOverlay), "rain_IR_")); + 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()); @@ -320,12 +340,12 @@ int main(int argc, char* argv[]) { 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()) { 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_IR_")); + imagesToSpread.push_back(ImageForSpread(ThreatImage::addRainOverlay(ch68, rainOverlay), "rain_68_")); } } else if(mLrptDecoder.isChannel64Available() && mLrptDecoder.isChannel65Available()) { @@ -352,14 +372,14 @@ int main(int argc, char* argv[]) { 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_IR_")); + 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_")); cv::Mat thermalRef = cv::imread(mSettings.getResourcesPath() + "thermal_ref.bmp"); cv::Mat thermalImage = ThreatImage::irToTemperature(ch68, thermalRef); @@ -378,283 +398,458 @@ int main(int argc, char* argv[]) { 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, mSettings.getScanAngle(), mSettings.getRoll(), mSettings.getPitch(), mSettings.getYaw(), 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(); + + 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); + 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); + 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); + 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.generateComposite321()) { - std::list images321; - std::list geolocationCalculators321; - searchForImages(images321, geolocationCalculators321, "321"); - - if(images321.size() > 1 && images321.size() == geolocationCalculators321.size()) { + if(mSettings.generateComposite221()) { + if(images.images221.size() > 1) { if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { - for(auto& img : images321) { + for(auto& img : images.images221) { img = ThreatImage::sharpen(img); } } - SpreadImage spreadImage; if(mSettings.compositeEquadistantProjection()) { - cv::Mat composite = spreadImage.equidistantProjection(images321, geolocationCalculators321, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate equidistant channel 321 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_321_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.jpg"; + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } + if(mSettings.compositeMercatorProjection()) { - cv::Mat composite = spreadImage.mercatorProjection(images321, geolocationCalculators321, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate mercator channel 321 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_321_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.jpg"; + 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.jpg"; + 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.jpg"; + 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.jpg"; + 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.jpg"; + 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.jpg"; + 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.jpg"; + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } } } - if(mSettings.generateComposite224()) { - std::list images224; - std::list geolocationCalculators224; - searchForImages(images224, geolocationCalculators224, "224"); - if(images224.size() > 1 && images224.size() == geolocationCalculators224.size()) { - SpreadImage spreadImage; + if(mSettings.generateComposite68()) { + if(images.images68.size() > 1) { + std::list irImages; if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { - for(auto& img : images224) { - img = ThreatImage::sharpen(img); + 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(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.jpg"; + 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.jpg"; + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); + } + } + + if(images.images67.size() > 1) { + std::list irImages; + if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { + for(const auto& img : images.images67) { + auto irImage = ThreatImage::equalize(img); + irImage = ThreatImage::invertIR(irImage); + irImages.emplace_back(irImage); } } if(mSettings.compositeEquadistantProjection()) { - cv::Mat composite = spreadImage.equidistantProjection(images224, geolocationCalculators224, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate equidistant channel 224 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "equidistant_" + compositeFileNameDateSS.str() + "_224_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.jpg"; + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } + if(mSettings.compositeMercatorProjection()) { - cv::Mat composite = spreadImage.mercatorProjection(images224, geolocationCalculators224, mSettings.getCompositeProjectionScale(), [](float progress) { - std::cout << "Generate mercator channel 224 composite image " << (int)progress << "% \t\t\r" << std::flush; - }); - std::cout << std::endl; - saveImage(mSettings.getOutputPath() + "mercator_" + compositeFileNameDateSS.str() + "_224_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.jpg"; + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } } } - if(mSettings.generateComposite68()) { - 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.jpg"; + 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.jpg"; + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); + } + } + if(images.images67.size() > 1) { + std::list thermalImages; 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) { + 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 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 : 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.jpg"; + 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 : 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.jpg"; + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } } } if(mSettings.generateComposite68Rain()) { - std::list images68; - std::list geolocationCalculators68; - searchForImages(images68, geolocationCalculators68, "68"); - - if(images68.size() > 1 && images68.size() == geolocationCalculators68.size()) { + if(images.images68.size() > 1) { + std::list irImages; + cv::Mat rainRef = cv::imread(mSettings.getResourcesPath() + "rain.bmp"); if(mSettings.compositeEquadistantProjection() || mSettings.compositeMercatorProjection()) { - cv::Mat rainRef = cv::imread(mSettings.getResourcesPath() + "rain.bmp"); - for(auto& img : images68) { + for(const auto& img : images.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); + 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); } } - 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 : 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.jpg"; + 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 : 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.jpg"; + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } } - } - if(mSettings.generateCompositeThermal()) { - std::list images68; - std::list geolocationCalculators68; - searchForImages(images68, geolocationCalculators68, "68"); - - if(images68.size() > 1 && images68.size() == geolocationCalculators68.size()) { + 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.jpg"; + 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.jpg"; + std::cout << "Saving composite: " << filePath << std::endl; + saveImage(filePath, composite); } } } @@ -665,47 +860,83 @@ 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()); + TleReader::TLE tle; + reader.processFile(); + if(!reader.getTLE(mSettings.getSatNameInTLE(), tle)) { + std::cout << "TLE data not found in TLE file, unable to create composite images..." << std::endl; + return result; + } + 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() == ".gcp" && fileCreatedSec < (mSettings.getCompositeMaxAgeHours() * 3600)) { + 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); - - do { - fs::path fileJPG(folder + "/" + fileNameBase + "_" + channelName + ".jpg"); - - if(fs::exists(fileJPG)) { - map[cftime] = std::make_tuple(entry.path().generic_string(), fileJPG.generic_string()); - break; - } - - fs::path fileBMP(folder + "/" + fileNameBase + "_" + channelName + ".bmp"); - - if(fs::exists(fileBMP)) { - map[cftime] = std::make_tuple(entry.path().generic_string(), fileBMP.generic_string()); - + std::string pathBase = folder + "/" + fileNameBase; + cv::Mat img; + + 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(!img.empty()) { + result.imageSizes.emplace_back(img.size()); + std::ifstream datFileStream(entry.path().generic_string()); + if(!datFileStream) { break; } - } 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))); + std::string line1, line2; + datFileStream >> line1; + datFileStream >> line2; + datFileStream.close(); + int64_t ticks = std::stoll(line1); + DateTime passStart(ticks); + ticks = std::stoll(line2); + TimeSpan passLength(ticks); + PixelGeolocationCalculator calc(tle, passStart, passLength, mSettings.getScanAngle(), mSettings.getRoll(), mSettings.getPitch(), mSettings.getYaw(), 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..8dbec7e --- /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 distance(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] * distance(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 cb42002..4c478cf 100644 --- a/resources/settings.ini +++ b/resources/settings.ini @@ -29,8 +29,8 @@ TimeOffset=0 [METEOR-M-2-3] SatNameInTLE=METEOR-M2 3 -ScanAngle=110.3 -Roll=0 +ScanAngle=109.9 +Roll=0.2 Pitch=0 Yaw=-0.4 # Onboard time offset in seconds @@ -74,12 +74,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;