From 3da653a0a370b99553dbfa6556e335e48cd6482c Mon Sep 17 00:00:00 2001 From: Roxeena Date: Tue, 18 Feb 2025 13:24:18 +0100 Subject: [PATCH] Avoid overflow of impact image --- data/tasks/impactcorridor.task | 23 ++ modules/neoviz/tasks/impactcorridortask.cpp | 418 ++++++++++++-------- modules/neoviz/tasks/impactcorridortask.h | 21 + 3 files changed, 294 insertions(+), 168 deletions(-) create mode 100644 data/tasks/impactcorridor.task diff --git a/data/tasks/impactcorridor.task b/data/tasks/impactcorridor.task new file mode 100644 index 0000000000..04740b5443 --- /dev/null +++ b/data/tasks/impactcorridor.task @@ -0,0 +1,23 @@ +--[[return { + { + Type = "ImpactCorridorTask", + KernelDirectory = "C:/Users/malej60/Documents/develop/b612-oribt-uncertainty/data/2023 CX1/openspace_variants/", + OutputFilename = "C:/Users/malej60/Documents/develop/b612-oribt-uncertainty/images/new/2023-CX1-test.png", + TimeIntervalStart = "2023-02-13T02:39:00.000", + TimeIntervalEnd = "2023-02-13T03:39:00.000", + ImageWidth = 5400, + ImageHeight = 2700 + } +}]]-- + +return { + { + Type = "ImpactCorridorTask", + KernelDirectory = "C:/Users/malej60/Documents/develop/b612-oribt-uncertainty/data/2004 MN4/openspace_variants_high_ip/", + OutputFilename = "C:/Users/malej60/Documents/develop/b612-oribt-uncertainty/images/new/2004-MN4-test.png", + TimeIntervalStart = "2029-02-01", + TimeIntervalEnd = "2029-08-17", + ImageWidth = 5400, + ImageHeight = 2700 + } +} diff --git a/modules/neoviz/tasks/impactcorridortask.cpp b/modules/neoviz/tasks/impactcorridortask.cpp index 0f0afb1ea5..24ccc2717c 100644 --- a/modules/neoviz/tasks/impactcorridortask.cpp +++ b/modules/neoviz/tasks/impactcorridortask.cpp @@ -38,20 +38,50 @@ #include "SpiceZpr.h" namespace { - constexpr std::string_view _loggerCat = "RenderableTube"; + constexpr std::string_view _loggerCat = "ImpactCorridorTask"; constexpr int MaxEarthRadius = 6378; - constexpr unsigned int WinSiz = 200; + constexpr unsigned int WinSiz = 200; constexpr std::string_view TimeStringFormat = "YYYY MON DD HR:MN:SC.###"; constexpr unsigned int TimeStringLength = 41; constexpr double ToDegrees = 180.0 / std::numbers::pi; constexpr int DebugMaxVariants = 100; - constexpr bool DebugMode = true; + constexpr bool DebugMode = false; // Offset spice Id to not collide with any existing NAIF id. This is the first NAIF ID // for the variants constexpr int IdOffset = 1000000; - int pixelIndex(int pixelW, int pixelH, int nChannels, int imageWidth, int imageHeight); + std::vector> gaussianFilter(int kernelSize, double filterStrength) { + ghoul_assert(kernelSize % 2 == 1, "Kernel size must be odd"); + + double sigma = filterStrength; + double s = 2.0 * sigma * sigma; + double sum = 0.0; + std::vector> gaussianKernel = std::vector>( + kernelSize, + std::vector(kernelSize, 0.0) + ); + + // Generating a kernelSize x kernelSize kernel + int halfSize = kernelSize / 2; + double r = 0.0; + for (int x = -halfSize; x <= halfSize; x++) { + for (int y = -halfSize; y <= halfSize; y++) { + r = sqrt(x * x + y * y); + gaussianKernel[x + halfSize][y + halfSize] = (exp(-(r * r) / s)) / (std::numbers::pi * s); + sum += gaussianKernel[x + halfSize][y + halfSize]; + } + } + + // normalising the Kernel + for (int i = 0; i < 5; ++i) { + for (int j = 0; j < 5; ++j) { + gaussianKernel[i][j] /= sum; + } + } + + return gaussianKernel; + } struct [[codegen::Dictionary(ImpactCorridorTask)]] Parameters { // Path to directory with kernels for the variants of the asteroid to test. The @@ -77,6 +107,10 @@ namespace { // The distance from Earth center to consider as an impact. If not given, then the // maximum radius of Earth is used. std::optional impactDistance; + + // The step size in number of seconds used to search for impacts. If not given, + // then the number of seconds in one day is used. + std::optional stepSize; }; #include "impactcorridortask_codegen.cpp" } // namespace @@ -87,10 +121,6 @@ documentation::Documentation ImpactCorridorTask::documentation() { return codegen::doc("neoviz_documentation_task"); } -std::string ImpactCorridorTask::description() { - return "Return the NEOviz impact corridor map image for the given input"; -} - ImpactCorridorTask::ImpactCorridorTask(const ghoul::Dictionary& dictionary) { const Parameters p = codegen::bake(dictionary); @@ -101,171 +131,27 @@ ImpactCorridorTask::ImpactCorridorTask(const ghoul::Dictionary& dictionary) { _imageWidth = p.imageWidth; _imageHeight = p.imageHeight; _impactDistance = p.impactDistance.value_or(MaxEarthRadius); + _stepSize = p.stepSize.value_or(spd_c()); } -int pixelIndex(int pixelW, int pixelH, int nChannels, int imageWidth, int imageHeight) { - if (pixelH < 0) { - pixelH = std::abs(pixelH); - - // Mirror in the 0 longitude line (Greenwich) - int greenwich = static_cast(std::round(imageWidth / 2.0)); - pixelW = greenwich + (greenwich - pixelW); - - } - else if (pixelH >= imageHeight) { - pixelH = imageHeight - (pixelH - imageHeight); - - // Mirror in the 0 longitude line (Greenwich) - int greenwich = static_cast(std::round(imageWidth / 2.0)); - pixelW = greenwich + (greenwich - pixelW); - } - - if (pixelW < 0) { - // Go over to the left side of the international daytime line - pixelW = imageWidth - std::abs(pixelW); - } - else if (pixelW >= imageWidth) { - // Go over to the right side of the international daytime line - pixelW = pixelW - imageWidth; - } - - int pixelIndex = nChannels * pixelW + nChannels * imageWidth * pixelH; - if (pixelIndex < 0 || pixelIndex >= (imageWidth * imageHeight * nChannels) ) { - throw ghoul::lua::LuaError("Pixel index out of bounds"); - } - return pixelIndex; +std::string ImpactCorridorTask::description() { + return "Return the NEOviz impact corridor map image for the given input"; } -void ImpactCorridorTask::perform(const Task::ProgressCallback & progressCallback) { - int nVarieants = std::distance( +void ImpactCorridorTask::perform(const Task::ProgressCallback& progressCallback) { + int nVariants = static_cast(std::distance( std::filesystem::directory_iterator(_kernelDirectory), std::filesystem::directory_iterator() - ); - LINFO(std::format("Number of variants: {}", nVarieants)); + )); + LINFO(std::format("Number of variants: {}", nVariants)); LINFO("Finding impacts..."); // Load the general kernels for the major bodies in the solar system SpiceManager::ref().loadKernel(absPath("${SYNC}/http/general_spk/2/de430.bsp")); SpiceManager::ref().loadKernel(absPath("${SYNC}/http/general_pck/1/pck00011.tpc")); - // Use a step size of 1 day (in seconds). - SpiceDouble step = spd_c(); - - // Specify the time range to search in - SPICEDOUBLE_CELL(timerange, WinSiz); - SpiceDouble startTime; - SpiceDouble endTime; - str2et_c(_timeIntervalStart.c_str(), &startTime); - str2et_c(_timeIntervalEnd.c_str(), &endTime); - wninsd_c(startTime, endTime, &timerange); - SpiceDouble impactStart; - SpiceDouble impactEnd; - SpiceChar impactTime[TimeStringLength]; - - // Variables to calculate and store the impact position - glm::dvec3 impactPosition = glm::dvec3(0.0); - SpiceDouble lightTime; - double impactLatitude = 0.0; - double impactLongitude = 0.0; - double impactAltitude = 0.0; - double impactLatitudeDegree = 0.0; - double impactLongitudeDegree = 0.0; - - // Loop over all variants of the asteroid - int counterVariants = 0; - int variantId = IdOffset; - std::string variantNAIFName; - int managerKernelId = 0; - SPICEDOUBLE_CELL(result, WinSiz); - for (const auto& variantKernel : - std::filesystem::directory_iterator(_kernelDirectory)) - { - if (DebugMode && counterVariants > DebugMaxVariants) { - break; - } - - variantNAIFName = std::to_string(variantId); - - // Load the kernel file for this variant - managerKernelId = SpiceManager::ref().loadKernel(variantKernel); - - // Get times within the searching timerange when the variant is closer than the - // impact distance to Earth. This is then considered as an impact. - gfdist_c( - variantNAIFName.c_str(), // Name of the target body - "NONE", // Aberration correction flag - "EARTH", // Name of the observing body - "<", // Relational operator (example <, = or >) - _impactDistance, // Reference value in km - 0.0, // Adjustment value for absolute extrema searches, not used - step, // Step size used for locating extrema and roots(The number of seconds in a day) - 100, // Workspace window interval count(from example) - &timerange, // SPICE window to which the search is confined(from example)); - &result // SPICE window containing results - ); - - // Check if the variant impacts or misses Earth - if (wncard_c(&result) == 0) { - // No impact - // LDEBUG(std::format("Variant {} does not impact Earth", variantNAIFName)); - ++variantId; - ++counterVariants; - SpiceManager::ref().unloadKernel(managerKernelId); - progressCallback(counterVariants / static_cast(nVarieants)); - continue; - } - - // Get the time range for when the variant is located within the impact distance - wnfetd_c(&result, 0, &impactStart, &impactEnd); - - // Then take the start of that timerange and convert it to a string, this gives - // the time of impact - timout_c( - impactStart, // Time in J2000 seconds - TimeStringFormat.data(), // Format for the output string - TimeStringLength, // Length of the output string plus 1 - impactTime // Result - ); - - // Get the cartesian X, Y, Z position of the variant at the impact time - spkpos_c( - variantNAIFName.c_str(), // Target body name to check position for - impactStart, // Time in J2000 seconds - "IAU_EARTH", // Reference frame of output position vector - "NONE", // Aberration correction flag - "EARTH", // Observing body name, reference frame of output position - glm::value_ptr(impactPosition), // Output position vector - &lightTime - ); - - // Convert the position to a lat, long, alt coordinate - reclat_c( - glm::value_ptr(impactPosition), - &impactAltitude, - &impactLongitude, - &impactLatitude - ); - impactLatitudeDegree = impactLatitude * ToDegrees; - impactLongitudeDegree = impactLongitude * ToDegrees; - - // Flip the y axis for the image, north is up - impactLatitudeDegree *= -1; - - // Add the result to the list of impact coordinates, this will be used to create - // the impact corridor image map later - _impactCoordinates.push_back({ - .id = variantId, - .time = impactTime, - .latitude = impactLatitudeDegree, - .longitude = impactLongitudeDegree - }); - - // Reset - SpiceManager::ref().unloadKernel(managerKernelId); - ++variantId; - ++counterVariants; - progressCallback(counterVariants / static_cast(nVarieants)); - } + // Find the impacts if there are any + findImpacts(progressCallback, nVariants); progressCallback(0.0); LINFO("Creating image..."); @@ -273,11 +159,11 @@ void ImpactCorridorTask::perform(const Task::ProgressCallback & progressCallback // Create a transparent image to then plot the impact coordinates on const unsigned int Size = _imageWidth * _imageHeight; const unsigned int NChannels = 4; - std::vector image = std::vector(Size * NChannels, 0.0); + std::vector image = std::vector(Size * NChannels, 0); // Tools to plot the impact coordinates on the image // @TODO: Make these parameters configurable - const int brushSize = 16; // Somewhere between 15 and 20 is good + const int brushSize = 15; // Somewhere between 15 and 20 is good const int brushSaturation = 100; // Tweak this to get a good result const int halfBrush = static_cast(std::round(brushSize / 2.0)); @@ -309,11 +195,15 @@ void ImpactCorridorTask::perform(const Task::ProgressCallback & progressCallback )); for (int c = 0; c < NChannels; ++c) { + // Check for overflow + if (static_cast(image[index + c]) + saturation > 255) { + // Clamp the color in case it will get too saturated + image[index + c] = 255; + continue; + } + // Add the saturation to the pixel image[index + c] += saturation; - - // Clamp the color in case it gets too saturated - image[index + c] = std::min(static_cast(image[index + c]), 255); } } } @@ -322,7 +212,14 @@ void ImpactCorridorTask::perform(const Task::ProgressCallback & progressCallback progressCallback(impactCounter / static_cast(_impactCoordinates.size())); } + // @TODO: Apply gaussian filter to the image + const int kernelSize = (halfBrush % 2 == 1) ? halfBrush : halfBrush + 1; + const double filterStrength = 1.0; + const std::vector> filter = + gaussianFilter(kernelSize, filterStrength); + // @TODO: Apply a transfer function to bring color to the image + // @TODO: Use the night layer image to estimate impact risk due to population density // Create and save the resulting impact map image @@ -337,9 +234,194 @@ void ImpactCorridorTask::perform(const Task::ProgressCallback & progressCallback // Estimate the impact probability, to compare with the expected value LINFO(std::format( - "Impact probability: {}", - _impactCoordinates.size() / static_cast(nVarieants) + "Impact probability: {} ({}/{})", + _impactCoordinates.size() / static_cast(nVariants), + _impactCoordinates.size(), nVariants )); } +void ImpactCorridorTask::findImpacts(const Task::ProgressCallback& progressCallback, + int nVariants) +{ + // Specify the time range to search in + SPICEDOUBLE_CELL(timerange, WinSiz); + SpiceDouble startTime; + SpiceDouble endTime; + str2et_c(_timeIntervalStart.c_str(), &startTime); + str2et_c(_timeIntervalEnd.c_str(), &endTime); + wninsd_c(startTime, endTime, &timerange); + SpiceDouble impactStart; + SpiceDouble impactEnd; + SpiceChar impactTime[TimeStringLength]; + + // Variables to calculate and store the impact position + glm::dvec3 impactPosition = glm::dvec3(0.0); + SpiceDouble lightTime; + double impactLatitude = 0.0; + double impactLongitude = 0.0; + double impactAltitude = 0.0; + double impactLatitudeDegree = 0.0; + double impactLongitudeDegree = 0.0; + + // Loop over all variants of the asteroid + int counterVariants = 0; + int variantId = IdOffset; + std::string variantNAIFName; + int managerKernelId = 0; + SPICEDOUBLE_CELL(result, WinSiz); + for (const auto& variantKernel : + std::filesystem::directory_iterator(_kernelDirectory)) + { + if (DebugMode && counterVariants >= DebugMaxVariants) { + break; + } + + variantNAIFName = std::to_string(variantId); + + // Make sure it is a kernel file and not anything else + if (!variantKernel.exists() || variantKernel.is_directory()) { + LWARNING(std::format( + "{} is not a file or could not be found", variantKernel.path() + )); + + ++variantId; + ++counterVariants; + SpiceManager::ref().unloadKernel(managerKernelId); + progressCallback(counterVariants / static_cast(nVariants)); + continue; + } + + // Load the kernel file for this variant + managerKernelId = SpiceManager::ref().loadKernel(variantKernel); + + // Get times within the searching timerange when the variant is closer than the + // impact distance to Earth. This is then considered as an impact. + gfdist_c( + variantNAIFName.c_str(), // Name of the target body + "NONE", // Aberration correction flag + "EARTH", // Name of the observing body + "<", // Relational operator (example <, =, or >) + _impactDistance, // Reference value in km + 0.0, // Adjustment value for absolute extrema searches + _stepSize, // Step size used for locating extrema and roots + 100, // Workspace window interval count + &timerange, // Time range for which the search is confined + &result // SPICE window containing results + ); + + // Check if the current kernel file and the current variant ID match + if (failed_c()) { + // If not, then there is probably a missing kernel file and we need to skip to + // the next valid variant + LWARNING(std::format("Missing kernel file {:06}", variantId - IdOffset + 1)); + variantId += 2; + counterVariants += 2; + + SpiceManager::ref().unloadKernel(managerKernelId); + progressCallback(counterVariants / static_cast(nVariants)); + reset_c(); + continue; + } + + // Check if the variant impacts or misses Earth + if (wncard_c(&result) == 0) { + // No impact + // LDEBUG(std::format("Variant {} does not impact Earth", variantNAIFName)); + ++variantId; + ++counterVariants; + SpiceManager::ref().unloadKernel(managerKernelId); + progressCallback(counterVariants / static_cast(nVariants)); + continue; + } + + // Get the time range for when the variant is located within the impact distance + wnfetd_c(&result, 0, &impactStart, &impactEnd); + + // Then take the start of that timerange and convert it to a string, this gives + // the time of impact + timout_c( + impactStart, // Time in J2000 seconds + TimeStringFormat.data(), // Format for the output string + TimeStringLength, // Length of the output string plus 1 + impactTime // Result + ); + + // Get the cartesian X, Y, Z position of the variant at the impact time + spkpos_c( + variantNAIFName.c_str(), // Target body name to check position for + impactStart, // Time in J2000 seconds + "IAU_EARTH", // Reference frame of output position vector + "NONE", // Aberration correction flag + "EARTH", // Observing body name, frame of output + glm::value_ptr(impactPosition), // Output position vector + &lightTime // Resulting light time between the positions + ); + + // Convert the position to a lat, long, alt coordinate + reclat_c( + glm::value_ptr(impactPosition), // Cartesian position + &impactAltitude, // Altitude + &impactLongitude, // Longitude + &impactLatitude // Latitude + ); + impactLatitudeDegree = impactLatitude * ToDegrees; + impactLongitudeDegree = impactLongitude * ToDegrees; + + // Flip the y axis for the image, north is up + impactLatitudeDegree *= -1; + + // Add the result to the list of impact coordinates, this will be used to create + // the impact corridor image map later + _impactCoordinates.push_back({ + .id = variantId, + .time = impactTime, + .latitude = impactLatitudeDegree, + .longitude = impactLongitudeDegree + }); + + // Reset + SpiceManager::ref().unloadKernel(managerKernelId); + ++variantId; + ++counterVariants; + progressCallback(counterVariants / static_cast(nVariants)); + } +} + +int ImpactCorridorTask::pixelIndex(int pixelW, int pixelH, int nChannels, int imageWidth, + int imageHeight) +{ + if (pixelH < 0) { + // Mirror in the noth pole line + pixelH = std::abs(pixelH); + + // Mirror in the 0 longitude line (Greenwich) + int greenwich = static_cast(std::round(imageWidth / 2.0)); + pixelW = greenwich + (greenwich - pixelW); + + } + else if (pixelH >= imageHeight) { + // Mirror in the south pole line + pixelH = imageHeight - (pixelH - imageHeight); + + // Mirror in the 0 longitude line (Greenwich) + int greenwich = static_cast(std::round(imageWidth / 2.0)); + pixelW = greenwich + (greenwich - pixelW); + } + + if (pixelW < 0) { + // Go over to the left side of the international daytime line + pixelW = imageWidth - std::abs(pixelW); + } + else if (pixelW >= imageWidth) { + // Go over to the right side of the international daytime line + pixelW = pixelW - imageWidth; + } + + int pixelIndex = nChannels * pixelW + nChannels * imageWidth * pixelH; + if (pixelIndex < 0 || pixelIndex >= (imageWidth * imageHeight * nChannels)) { + throw ghoul::RuntimeError("Pixel index out of bounds"); + } + return pixelIndex; +} + } // namespace openspace::neoviz diff --git a/modules/neoviz/tasks/impactcorridortask.h b/modules/neoviz/tasks/impactcorridortask.h index 0f1a7ccd6a..d4fe43b431 100644 --- a/modules/neoviz/tasks/impactcorridortask.h +++ b/modules/neoviz/tasks/impactcorridortask.h @@ -49,9 +49,30 @@ private: double longitude = 0.0; }; + void findImpacts(const Task::ProgressCallback& progressCallback, int nVariants); + + /** + * Retrieve the index in the flat data list of pixels that cooresponds to a given + * pixel coordinate. If the given pixel coordinate is outside the bounds of the image + * then it will be properly wrapped/mirrored around the image considering that it + * represents a map. If the pixel croses the international day time line then it will + * be wrapped around to the other side. If instead the pixel crosses the pole, the it + * will be mirrored around Greenwich. + * + * \param pixelW The horizontal coordinate for the pixel in the image + * \param pixelH The vertical coordinate for the pixel in the image + * \param nChannels The number of channels in the image + * \param imageWidth The total width of the image in pixels + * \param imageHeight The total height of the image in pixels + * \return The index of the first channel of the given pixel in the flat data list + */ + int pixelIndex(int pixelW, int pixelH, int nChannels, int imageWidth, + int imageHeight); + std::filesystem::path _kernelDirectory; std::filesystem::path _outputFilename; int _impactDistance; + double _stepSize; std::string _timeIntervalStart; std::string _timeIntervalEnd; std::vector _impactCoordinates;