From d2c3d9c4bed18cf65a5b6a3c6b435e915592b112 Mon Sep 17 00:00:00 2001 From: Cheng-Chin Chiang Date: Wed, 10 Apr 2024 02:50:25 +0800 Subject: [PATCH 1/5] Store the global settings in a singleton class (#1321) * Add a singleton object that stores the ProgramSettings as global settings * Apply the global setting for log protocol messages * Modify the changelog * Combine two classes GlobalValues and Global * Directly access global values from the Session without duplicating them * Remove Global class and add an GetInstance static function in ProgramSettings * Fix a build error * Rename global_values to settings * Add a static function Initialise to ProgramSettings * Restore Session functions SetExitTimeout and SetInitExitTimeout * Read top level folder and starting folder from the global settings object inside the TableController * Add a _settings, which refers to the global settings object, to the Session * Creating a local reference variable for ProgramSettings::GetInstance() * Restore _settings in the SessionManager * Add a _settings property to the TableController * Add a settings to logger functions * Fix a bug * Initialize the _settings in TableController via ProgramSettings::GetInstance() * Pass program settings via ProgramSettings::GetInstance() in logger functions * Set a static variable log_protocol_messages and copy its value from program settings * Remove _settings references from Session * Remove _settings references from TableController * Copy global setting parameters in the InitLogger function --- CHANGELOG.md | 1 + src/Logger/Logger.cc | 18 +++++++++++++----- src/Logger/Logger.h | 2 +- src/Main/Main.cc | 5 ++--- src/Main/ProgramSettings.h | 9 +++++++++ src/Session/Session.cc | 18 +++++++++++------- src/Session/Session.h | 24 +++++++++++------------- src/Session/SessionManager.cc | 3 +-- src/Table/TableController.cc | 8 ++++++-- src/Table/TableController.h | 7 ++++--- test/TestFileInfo.cc | 2 +- test/TestMain.cc | 20 +++++++++----------- 12 files changed, 69 insertions(+), 48 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index da7139c7b..e1ffbf6db 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -12,6 +12,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Remove program settings equality operators ([#1001](https://github.com/CARTAvis/carta-backend/issues/1001)). * Normalize the style of guard names in header files ([#1023](https://github.com/CARTAvis/carta-backend/issues/1023)). * Improved file IDs for generated images ([#1224](https://github.com/CARTAvis/carta-frontend/issues/1224)). +* Store global settings in a singleton class ([#1302](https://github.com/CARTAvis/carta-backend/issues/1302)). ## [4.1.0] diff --git a/src/Logger/Logger.cc b/src/Logger/Logger.cc index 9278e26d1..f2ea51392 100644 --- a/src/Logger/Logger.cc +++ b/src/Logger/Logger.cc @@ -6,6 +6,8 @@ #include "Logger.h" +#include "Main/ProgramSettings.h" + #include namespace carta { @@ -13,8 +15,14 @@ namespace logger { static bool log_protocol_messages(false); -void InitLogger(bool no_log_file, int verbosity, bool log_performance, bool log_protocol_messages_, fs::path user_directory) { - log_protocol_messages = log_protocol_messages_; +void InitLogger() { + // Copy parameters from the global settings + auto& settings = ProgramSettings::GetInstance(); + log_protocol_messages = settings.log_protocol_messages; + auto no_log = settings.no_log; + auto user_directory = settings.user_directory; + auto verbosity = settings.verbosity; + auto log_performance = settings.log_performance; // Set the stdout/stderr console auto console_sink = std::make_shared(); @@ -26,7 +34,7 @@ void InitLogger(bool no_log_file, int verbosity, bool log_performance, bool log_ // Set a log file with its full name, maximum size and the number of rotated files std::string log_fullname; - if (!no_log_file) { + if (!no_log) { log_fullname = (user_directory / "log/carta.log").string(); auto stdout_log_file_sink = std::make_shared(log_fullname, LOG_FILE_SIZE, ROTATED_LOG_FILES); stdout_log_file_sink->set_formatter( @@ -72,7 +80,7 @@ void InitLogger(bool no_log_file, int verbosity, bool log_performance, bool log_ // Set as the default logger spdlog::set_default_logger(default_logger); - if (!no_log_file) { + if (!no_log) { spdlog::info("Writing to the log file: {}", log_fullname); } @@ -87,7 +95,7 @@ void InitLogger(bool no_log_file, int verbosity, bool log_performance, bool log_ // Set a log file with its full name, maximum size and the number of rotated files std::string perf_log_fullname; - if (!no_log_file) { + if (!no_log) { perf_log_fullname = (user_directory / "log/performance.log").string(); auto perf_log_file_sink = std::make_shared(perf_log_fullname, LOG_FILE_SIZE, ROTATED_LOG_FILES); diff --git a/src/Logger/Logger.h b/src/Logger/Logger.h index 67ff23135..63125d3df 100644 --- a/src/Logger/Logger.h +++ b/src/Logger/Logger.h @@ -110,7 +110,7 @@ class carta_sink : public ansicolor_sink { namespace carta { namespace logger { -void InitLogger(bool no_log_file, int verbosity, bool log_performance, bool log_protocol_messages_, fs::path user_directory); +void InitLogger(); void LogReceivedEventType(const CARTA::EventType& event_type); void LogSentEventType(const CARTA::EventType& event_type); void FlushLogFile(); diff --git a/src/Main/Main.cc b/src/Main/Main.cc index 1c9f944c6..4cf4e43fb 100644 --- a/src/Main/Main.cc +++ b/src/Main/Main.cc @@ -46,14 +46,13 @@ int main(int argc, char* argv[]) { sigaction(SIGINT, &sig_handler, nullptr); // Main - carta::ProgramSettings settings(argc, argv); + auto settings = ProgramSettings::Initialise(argc, argv); if (settings.help || settings.version) { exit(0); } - carta::logger::InitLogger( - settings.no_log, settings.verbosity, settings.log_performance, settings.log_protocol_messages, settings.user_directory); + carta::logger::InitLogger(); settings.FlushMessages(); // flush log messages produced during Program Settings setup // Send casacore log messages (global and local) to sink. diff --git a/src/Main/ProgramSettings.h b/src/Main/ProgramSettings.h index aa0822289..313ef1c0d 100644 --- a/src/Main/ProgramSettings.h +++ b/src/Main/ProgramSettings.h @@ -31,6 +31,15 @@ namespace carta { struct ProgramSettings { + static ProgramSettings& GetInstance() { + static ProgramSettings settings; + return settings; + } + + static ProgramSettings& Initialise(int argc, char** argv) { + return GetInstance() = std::move(carta::ProgramSettings(argc, argv)); + } + bool version = false; bool help = false; std::vector port; diff --git a/src/Session/Session.cc b/src/Session/Session.cc index 423f81069..d57fa40a6 100644 --- a/src/Session/Session.cc +++ b/src/Session/Session.cc @@ -48,17 +48,12 @@ bool Session::_controller_deployment = false; std::thread* Session::_animation_thread = nullptr; Session::Session(uWS::WebSocket* ws, uWS::Loop* loop, uint32_t id, std::string address, - std::string top_level_folder, std::string starting_folder, std::shared_ptr file_list_handler, bool read_only_mode, - bool enable_scripting) + std::shared_ptr file_list_handler) : _socket(ws), _loop(loop), _id(id), _address(address), - _top_level_folder(top_level_folder), - _starting_folder(starting_folder), - _table_controller(std::make_unique(_top_level_folder, _starting_folder)), - _read_only_mode(read_only_mode), - _enable_scripting(enable_scripting), + _table_controller(std::make_unique()), _region_handler(nullptr), _file_list_handler(file_list_handler), _sync_id(0), @@ -66,6 +61,10 @@ Session::Session(uWS::WebSocket* ws, uWS::Loop* loop _animation_active(false), _cursor_settings(this), _loaders(LOADER_CACHE_SIZE) { + auto& settings = ProgramSettings::GetInstance(); + _top_level_folder = settings.top_level_folder; + _read_only_mode = settings.read_only_mode; + _enable_scripting = settings.enable_scripting; _histogram_progress = 1.0; _ref_count = 0; _animation_object = nullptr; @@ -125,6 +124,11 @@ Session::~Session() { logger::FlushLogFile(); } +void Session::SetExitTimeout(int secs) { + _exit_after_num_seconds = secs; + _exit_when_all_sessions_closed = true; +} + void Session::SetInitExitTimeout(int secs) { __exit_backend_timer = secs; struct sigaction sig_handler; diff --git a/src/Session/Session.h b/src/Session/Session.h index 697a73e7b..b6455ccee 100644 --- a/src/Session/Session.h +++ b/src/Session/Session.h @@ -29,13 +29,13 @@ #include "FileList/FileListHandler.h" #include "Frame/Frame.h" #include "ImageData/StokesFilesConnector.h" +#include "Main/ProgramSettings.h" #include "Region/RegionHandler.h" #include "SessionContext.h" +#include "Table/TableController.h" #include "ThreadingManager/Concurrency.h" #include "Util/Message.h" -#include "Table/TableController.h" - #define HISTOGRAM_CANCEL -1.0 #define UPDATE_HISTOGRAM_PROGRESS_PER_SECONDS 2.0 #define LOADER_CACHE_SIZE 25 @@ -52,9 +52,8 @@ struct PerSocketData { class Session { public: - Session(uWS::WebSocket* ws, uWS::Loop* loop, uint32_t id, std::string address, std::string top_level_folder, - std::string starting_folder, std::shared_ptr file_list_handler, bool read_only_mode = false, - bool enable_scripting = false); + Session(uWS::WebSocket* ws, uWS::Loop* loop, uint32_t id, std::string address, + std::shared_ptr file_list_handler); ~Session(); // CARTA ICD @@ -181,10 +180,8 @@ class Session { return _animation_object && !_animation_object->_stop_called; } int CalculateAnimationFlowWindow(); - static void SetExitTimeout(int secs) { - _exit_after_num_seconds = secs; - _exit_when_all_sessions_closed = true; - } + + static void SetExitTimeout(int secs); static void SetInitExitTimeout(int secs); static void SetControllerDeploymentFlag(bool controller_deployment) { @@ -278,10 +275,6 @@ class Session { uint32_t _id; std::string _address; - std::string _top_level_folder; - std::string _starting_folder; - bool _read_only_mode; - bool _enable_scripting; // File browser std::shared_ptr _file_list_handler; @@ -340,6 +333,11 @@ class Session { // Timestamp for the last protobuf message std::chrono::high_resolution_clock::time_point _last_message_timestamp; + + // Parameters which are copied from the global settings + std::string _top_level_folder; + bool _read_only_mode; + bool _enable_scripting; }; } // namespace carta diff --git a/src/Session/SessionManager.cc b/src/Session/SessionManager.cc index 01ba862c9..acf0804fa 100644 --- a/src/Session/SessionManager.cc +++ b/src/Session/SessionManager.cc @@ -87,8 +87,7 @@ void SessionManager::OnConnect(WSType* ws) { // create a Session std::unique_lock ulock(_sessions_mutex); - _sessions[session_id] = new Session(ws, loop, session_id, address, _settings.top_level_folder, _settings.starting_folder, - _file_list_handler, _settings.read_only_mode, _settings.enable_scripting); + _sessions[session_id] = new Session(ws, loop, session_id, address, _file_list_handler); _sessions[session_id]->IncreaseRefCount(); diff --git a/src/Table/TableController.cc b/src/Table/TableController.cc index 253a4371f..0268d1d6f 100644 --- a/src/Table/TableController.cc +++ b/src/Table/TableController.cc @@ -9,6 +9,7 @@ #include #include "Logger/Logger.h" +#include "Main/ProgramSettings.h" #include "Timer/ListProgressReporter.h" #include "Util/File.h" @@ -18,8 +19,11 @@ using namespace carta; -TableController::TableController(const std::string& top_level_folder, const std::string& starting_folder) - : _top_level_folder(top_level_folder), _starting_folder(starting_folder) {} +TableController::TableController() { + auto& settings = ProgramSettings::GetInstance(); + _top_level_folder = settings.top_level_folder; + _starting_folder = settings.starting_folder; +} void TableController::OnOpenFileRequest(const CARTA::OpenCatalogFile& open_file_request, CARTA::OpenCatalogFileAck& open_file_response) { int file_id = open_file_request.file_id(); diff --git a/src/Table/TableController.h b/src/Table/TableController.h index 01c56bb6d..a2685b4b9 100644 --- a/src/Table/TableController.h +++ b/src/Table/TableController.h @@ -16,6 +16,7 @@ #include #include +#include "Main/ProgramSettings.h" #include "Table.h" #include "Util/FileSystem.h" @@ -32,7 +33,7 @@ struct TableViewCache { class TableController { public: - TableController(const std::string& top_level_folder, const std::string& starting_folder); + TableController(); void OnFileListRequest(const CARTA::CatalogListRequest& file_list_request, CARTA::CatalogListResponse& file_list_response); void OnFileInfoRequest(const CARTA::CatalogFileInfoRequest& file_info_request, CARTA::CatalogFileInfoResponse& file_info_response); void OnOpenFileRequest(const CARTA::OpenCatalogFile& open_file_request, CARTA::OpenCatalogFileAck& open_file_response); @@ -53,8 +54,6 @@ class TableController { static bool FilterParamsChanged(const std::vector& filter_configs, std::string sort_column, CARTA::SortingType sorting_type, const TableViewCache& cached_config); fs::path GetPath(std::string directory, std::string name = ""); - std::string _top_level_folder; - std::string _starting_folder; std::unordered_map _tables; std::unordered_map _view_cache; @@ -62,6 +61,8 @@ class TableController { volatile bool _stop_getting_file_list; volatile bool _first_report_made; std::function _progress_callback; + std::string _top_level_folder; + std::string _starting_folder; }; } // namespace carta #endif // CARTA_SRC_TABLE_TABLECONTROLLER_H_ diff --git a/test/TestFileInfo.cc b/test/TestFileInfo.cc index 4c982cafc..d1d60fefc 100644 --- a/test/TestFileInfo.cc +++ b/test/TestFileInfo.cc @@ -179,7 +179,7 @@ class SessionFileInfoTest : public ::testing::Test { class TestSession : public Session { public: - TestSession() : Session(nullptr, nullptr, 0, "", "/", "", nullptr, -1, false) {} + TestSession() : Session(nullptr, nullptr, 0, "", nullptr) {} void TestFileInfo(const std::string& request_filename, const CARTA::FileType& request_file_type, const std::string& request_hdu = "") { auto request = Message::FileInfoRequest(SAMPLE_FILES_PATH, request_filename, request_hdu); diff --git a/test/TestMain.cc b/test/TestMain.cc index 1440d3896..668fff30a 100644 --- a/test/TestMain.cc +++ b/test/TestMain.cc @@ -12,6 +12,7 @@ #include "CommonTestUtilities.h" #include "Logger/Logger.h" +#include "Main/ProgramSettings.h" #include "ThreadingManager/ThreadingManager.h" #define TASK_THREAD_COUNT 3 @@ -22,10 +23,6 @@ int main(int argc, char** argv) { testing::AddGlobalTestEnvironment(new CartaEnvironment()); // Set global environment - int verbosity(0); - bool no_log(true); - bool log_performance(false); - bool log_protocol_messages(false); int omp_threads(omp_get_num_procs()); cxxopts::Options options("carta-icd-test", "CARTA ICD test"); @@ -34,7 +31,7 @@ int main(int argc, char** argv) { options.add_options() ("h,help", "print usage") ("verbosity", "display verbose logging from this level", - cxxopts::value()->default_value(std::to_string(verbosity)), "") + cxxopts::value()->default_value(std::to_string(0)), "") ("no_log", "do not log output to a log file", cxxopts::value()->default_value("true")) ("log_performance", "enable performance debug logs", cxxopts::value()->default_value("false")) ("log_protocol_messages", "enable protocol message debug logs", cxxopts::value()->default_value("false")) @@ -49,10 +46,11 @@ int main(int argc, char** argv) { return 0; } - verbosity = result["verbosity"].as(); - no_log = result["no_log"].as(); - log_performance = result["log_performance"].as(); - log_protocol_messages = result["log_protocol_messages"].as(); + auto& settings = ProgramSettings::GetInstance(); + settings.verbosity = result["verbosity"].as(); + settings.no_log = result["no_log"].as(); + settings.log_performance = result["log_performance"].as(); + settings.log_protocol_messages = result["log_protocol_messages"].as(); omp_threads = result["omp_threads"].as(); if (omp_threads < 0) { @@ -62,8 +60,8 @@ int main(int argc, char** argv) { carta::ThreadManager::StartEventHandlingThreads(TASK_THREAD_COUNT); carta::ThreadManager::SetThreadLimit(omp_threads); - fs::path user_directory = fs::path(getenv("HOME")) / CARTA_USER_FOLDER_PREFIX; - carta::logger::InitLogger(no_log, verbosity, log_performance, log_protocol_messages, user_directory); + ProgramSettings::GetInstance().user_directory = fs::path(getenv("HOME")) / CARTA_USER_FOLDER_PREFIX; + carta::logger::InitLogger(); int run_all_tests = RUN_ALL_TESTS(); From de4da5280f09650edec5567362a0cce6b3a54e23 Mon Sep 17 00:00:00 2001 From: pford Date: Tue, 9 Apr 2024 13:16:17 -0600 Subject: [PATCH 2/5] Refactor Region code to RegionState and RegionConverter (#1352) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Rename spatial profiles test, for cursor only * Rename region (line and point) spatial profile test * Remove unused methods, update test CMakeLists * Add region histogram and stats tests * Tweak tests and rename test suite * Add tests for Region in reference image * Add tests for matched regions and rotbox * Convert test utils CmpValues and CmpVectors to template * Add region spectral profile test * Remove float versions of test utils CmpValues and CmpVectors * Fix formatting * Fix includes after removing Message.h from Region.cc * Use FileFinder for image path in tests * Move rotbox tests into region and matched region tests * Refactor Region into RegionConverter and add tests * Remove debug output * Add performance time for line spatial profile * Remove tab from CMakeLists * Add spatial profile test for point outside image * Add region import/export test and bug fixes from refactoring * Fix formatting * Add changelog entry * Fix PvPreviewCut profile using nchan in preview image * Fix ellipse record from control points for export Enforce maj>min as in WCEllipsoid * Fix point region spectral profile and add test * Fix format * Update protobuf * Use double for ellipse angle in record for accurate unit conversion --------- Co-authored-by: Pam Harris Co-authored-by: Adrianna Pińska --- CHANGELOG.md | 1 + CMakeLists.txt | 1 + src/ImageGenerators/PvPreviewCube.cc | 5 +- src/ImageGenerators/PvPreviewCube.h | 4 +- src/ImageGenerators/PvPreviewCut.h | 2 + src/Region/CrtfImportExport.cc | 9 +- src/Region/LineBoxRegions.cc | 1 + src/Region/Region.cc | 1325 ++--------------- src/Region/Region.h | 205 +-- src/Region/RegionConverter.cc | 1012 +++++++++++++ src/Region/RegionConverter.h | 103 ++ src/Region/RegionHandler.cc | 26 +- src/Region/RegionHandler.h | 1 - src/Region/RegionImportExport.cc | 27 +- src/Region/RegionState.h | 149 ++ test/CMakeLists.txt | 10 +- test/CommonTestUtilities.cc | 23 +- test/CommonTestUtilities.h | 7 +- test/CommonTestUtilities.tcc | 17 + ...ofiles.cc => TestCursorSpatialProfiles.cc} | 104 +- test/TestExprImage.cc | 4 +- test/TestLineSpatialProfiles.cc | 222 --- test/TestPvGenerator.cc | 21 +- test/TestRegion.cc | 418 ++++++ test/TestRegionHistogram.cc | 94 ++ test/TestRegionImportExport.cc | 362 +++++ test/TestRegionMatched.cc | 447 ++++++ test/TestRegionSpatialProfiles.cc | 346 +++++ test/TestRegionSpectralProfiles.cc | 173 +++ test/TestRegionStats.cc | 108 ++ 30 files changed, 3543 insertions(+), 1684 deletions(-) create mode 100644 src/Region/RegionConverter.cc create mode 100644 src/Region/RegionConverter.h create mode 100644 src/Region/RegionState.h rename test/{TestSpatialProfiles.cc => TestCursorSpatialProfiles.cc} (88%) delete mode 100644 test/TestLineSpatialProfiles.cc create mode 100644 test/TestRegion.cc create mode 100644 test/TestRegionHistogram.cc create mode 100644 test/TestRegionImportExport.cc create mode 100644 test/TestRegionMatched.cc create mode 100644 test/TestRegionSpatialProfiles.cc create mode 100644 test/TestRegionSpectralProfiles.cc create mode 100644 test/TestRegionStats.cc diff --git a/CHANGELOG.md b/CHANGELOG.md index e1ffbf6db..3f7162879 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Normalize the style of guard names in header files ([#1023](https://github.com/CARTAvis/carta-backend/issues/1023)). * Improved file IDs for generated images ([#1224](https://github.com/CARTAvis/carta-frontend/issues/1224)). * Store global settings in a singleton class ([#1302](https://github.com/CARTAvis/carta-backend/issues/1302)). +* Move Region code to RegionConverter for maintenance and performance ([#1347](https://github.com/CARTAvis/carta-backend/issues/1347)). ## [4.1.0] diff --git a/CMakeLists.txt b/CMakeLists.txt index a2807bcd3..5b8c7d2de 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -233,6 +233,7 @@ set(SOURCE_FILES src/Region/Ds9ImportExport.cc src/Region/LineBoxRegions.cc src/Region/Region.cc + src/Region/RegionConverter.cc src/Region/RegionHandler.cc src/Region/RegionImportExport.cc src/Session/CursorSettings.cc diff --git a/src/ImageGenerators/PvPreviewCube.cc b/src/ImageGenerators/PvPreviewCube.cc index 0f66b96f3..394716d93 100644 --- a/src/ImageGenerators/PvPreviewCube.cc +++ b/src/ImageGenerators/PvPreviewCube.cc @@ -12,6 +12,7 @@ #include #include "DataStream/Smoothing.h" +#include "Logger/Logger.h" #include "Timer/Timer.h" #define LOAD_DATA_PROGRESS_INTERVAL 1000 @@ -162,8 +163,8 @@ bool PvPreviewCube::GetRegionProfile(const casacore::Slicer& region_bounding_box auto box_start = region_bounding_box.start(); auto box_length = region_bounding_box.length(); - // Initialize profile - size_t nchan = box_length(_preview_image->coordinates().spectralAxisNumber()); + // Initialize profile to channels in preview image + size_t nchan = _preview_image->shape()(_preview_image->coordinates().spectralAxisNumber()); profile.resize(nchan, NAN); std::vector npix_per_chan(nchan, 0.0); auto data_shape = _cube_data.shape(); diff --git a/src/ImageGenerators/PvPreviewCube.h b/src/ImageGenerators/PvPreviewCube.h index d39bf2702..32622317d 100644 --- a/src/ImageGenerators/PvPreviewCube.h +++ b/src/ImageGenerators/PvPreviewCube.h @@ -7,13 +7,13 @@ #ifndef CARTA_SRC_IMAGEGENERATORS_PVPREVIEWCUBE_H_ #define CARTA_SRC_IMAGEGENERATORS_PVPREVIEWCUBE_H_ +#include + #include "ImageGenerators/ImageGenerator.h" #include "Region/Region.h" #include "Util/File.h" #include "Util/Image.h" -#include - namespace carta { struct PreviewCubeParameters { diff --git a/src/ImageGenerators/PvPreviewCut.h b/src/ImageGenerators/PvPreviewCut.h index 9e2f9f9ed..6e4fc667a 100644 --- a/src/ImageGenerators/PvPreviewCut.h +++ b/src/ImageGenerators/PvPreviewCut.h @@ -7,6 +7,8 @@ #ifndef CARTA_SRC_IMAGEGENERATORS_PVPREVIEWCUT_H_ #define CARTA_SRC_IMAGEGENERATORS_PVPREVIEWCUT_H_ +#include + #include "Region/Region.h" #include "Util/File.h" diff --git a/src/Region/CrtfImportExport.cc b/src/Region/CrtfImportExport.cc index c1727857b..ba3957628 100644 --- a/src/Region/CrtfImportExport.cc +++ b/src/Region/CrtfImportExport.cc @@ -451,7 +451,6 @@ void CrtfImportExport::ProcessFileLines(std::vector& lines) { std::vector parameters; std::unordered_map properties; ParseRegionParameters(line, parameters, properties); - // Coordinate frame for world coordinates conversion RegionState region_state; CARTA::RegionStyle region_style; @@ -499,7 +498,13 @@ void CrtfImportExport::ProcessFileLines(std::vector& lines) { if (region == "text") { // Set text label for text region if (parameters.size() == 4) { // text, x, y, "text" - region_style.mutable_annotation_style()->set_text_label0(parameters[3]); + std::string label(parameters[3]); + if (!label.empty() && + ((label.front() == '"' && label.back() == '"') || (label.front() == '\'' && label.back() == '\''))) { + label.pop_back(); + label = label.substr(1); + } + region_style.mutable_annotation_style()->set_text_label0(label); } } else { // Set text position for textbox region diff --git a/src/Region/LineBoxRegions.cc b/src/Region/LineBoxRegions.cc index 7d14f38b3..f0e893483 100644 --- a/src/Region/LineBoxRegions.cc +++ b/src/Region/LineBoxRegions.cc @@ -11,6 +11,7 @@ #include #include "Logger/Logger.h" +#include "Util/Message.h" namespace carta { diff --git a/src/Region/Region.cc b/src/Region/Region.cc index 339f6ed2d..157d7334f 100644 --- a/src/Region/Region.cc +++ b/src/Region/Region.cc @@ -8,28 +8,16 @@ #include "Region.h" -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include #include -#include -#include +#include -#include "../Logger/Logger.h" -#include "Util/Image.h" +#include "RegionConverter.h" +// #include "Util/Image.h" using namespace carta; Region::Region(const RegionState& state, std::shared_ptr csys) - : _coord_sys(csys), _valid(false), _region_changed(false), _reference_region_set(false), _region_state(state) { + : _coord_sys(csys), _valid(false), _region_changed(false), _lcregion_set(false), _region_state(state) { _valid = CheckPoints(state.control_points, state.type); } @@ -60,12 +48,10 @@ bool Region::UpdateRegion(const RegionState& new_state) { void Region::ResetRegionCache() { // Invalid when region changes - _reference_region_set = false; - std::lock_guard guard(_region_mutex); - _wcs_control_points.clear(); - _reference_region.reset(); - _applied_regions.clear(); - _polygon_regions.clear(); + std::lock_guard guard(_lcregion_mutex); + _lcregion.reset(); + _lcregion_set = false; + _region_converter.reset(); } // ************************************************************************* @@ -128,6 +114,7 @@ bool Region::PointsFinite(const std::vector& points) { return points_finite; } +// ************************************************************************* // Region connection state (disconnected when region closed) bool Region::IsConnected() { @@ -139,645 +126,94 @@ void Region::WaitForTaskCancellation() { // to interrupt the running jobs in the std::unique_lock lock(GetActiveTaskMutex()); } -// ****************************************************************************************** -// Apply region to reference image in world coordinates (WCRegion) and save wcs control points - -bool Region::ReferenceRegionValid() { - return _reference_region_set && bool(_reference_region); -} - -void Region::SetReferenceRegion() { - // Create WCRegion (world coordinate region) in the reference image according to type using wcs control points - // Sets _reference_region (maybe to nullptr) or _reference_annotation - casacore::WCRegion* region(nullptr); - auto region_state = GetRegionState(); - std::vector pixel_points(region_state.control_points); - std::vector world_points; // point holder; one CARTA point is two world points (x, y) - casacore::IPosition pixel_axes(2, 0, 1); - casacore::Vector abs_rel; - auto type(region_state.type); - try { - switch (type) { - case CARTA::POINT: - case CARTA::ANNPOINT: { - if (ConvertCartaPointToWorld(pixel_points[0], _wcs_control_points)) { - // WCBox blc and trc are same point - std::lock_guard guard(_region_mutex); - casacore::Vector const wcs_points(_wcs_control_points); - region = new casacore::WCBox(wcs_points, wcs_points, pixel_axes, *_coord_sys, abs_rel); - } - break; - } - case CARTA::RECTANGLE: // 4 corners - case CARTA::POLYGON: // vertices - case CARTA::ANNRECTANGLE: // 4 corners - case CARTA::ANNPOLYGON: // vertices - case CARTA::ANNTEXT: { // 4 corners of text box - if (type == CARTA::RECTANGLE || type == CARTA::ANNRECTANGLE || type == CARTA::ANNTEXT) { - if (!RectanglePointsToWorld(pixel_points, _wcs_control_points)) { - _wcs_control_points.clear(); - } - } else { - for (auto& point : pixel_points) { - if (ConvertCartaPointToWorld(point, world_points)) { - _wcs_control_points.push_back(world_points[0]); - _wcs_control_points.push_back(world_points[1]); - } else { - _wcs_control_points.clear(); - break; - } - } - } - - if (!_wcs_control_points.empty()) { - // separate x and y in control points, convert from Vector to Quantum - std::vector x, y; - for (size_t i = 0; i < _wcs_control_points.size(); i += 2) { - x.push_back(_wcs_control_points[i].getValue()); - y.push_back(_wcs_control_points[i + 1].getValue()); - } - casacore::Vector vx(x), vy(y); - casacore::Quantum> qx, qy; - - casacore::Vector world_units = _coord_sys->worldAxisUnits(); - - qx = vx; // set values - qx.setUnit(world_units(0)); // set unit - qy = vy; // set values - qy.setUnit(world_units(1)); // set unit - - std::lock_guard guard(_region_mutex); - region = new casacore::WCPolygon(qx, qy, pixel_axes, *_coord_sys); - } - break; - } - case CARTA::ELLIPSE: - case CARTA::ANNELLIPSE: // [(cx, cy), (bmaj, bmin)] - case CARTA::ANNCOMPASS: { // [(cx, cy), (length, length)} - float ellipse_rotation; - if (EllipsePointsToWorld(pixel_points, _wcs_control_points, ellipse_rotation)) { - // control points are in order: xcenter, ycenter, major axis, minor axis - casacore::Quantity theta(ellipse_rotation, "deg"); - theta.convert("rad"); - std::lock_guard guard(_region_mutex); - region = new casacore::WCEllipsoid(_wcs_control_points[0], _wcs_control_points[1], _wcs_control_points[2], - _wcs_control_points[3], theta, 0, 1, *_coord_sys); - } - break; - } - default: - break; - } - } catch (casacore::AipsError& err) { // region failed - spdlog::error("{} region failed: {}", RegionName(type), err.getMesg()); - } - - std::shared_ptr shared_region = std::shared_ptr(region); - std::atomic_store(&_reference_region, shared_region); - - // Flag indicates that attempt was made, to avoid repeated attempts - _reference_region_set = true; -} - -bool Region::RectanglePointsToWorld(std::vector& pixel_points, std::vector& wcs_points) { - // Convert CARTA rectangle points (cx, cy), (width, height) to world coordinates - if (pixel_points.size() != 2) { - return false; - } - - // Get 4 corner points in pixel coordinates - float rotation(GetRegionState().rotation); - casacore::Vector x, y; - RectanglePointsToCorners(pixel_points, rotation, x, y); - - // Convert corners to wcs in one call for efficiency - size_t num_points(x.size()), num_axes(_coord_sys->nPixelAxes()); - casacore::Matrix pixel_coords(num_axes, num_points); - casacore::Matrix world_coords(num_axes, num_points); - pixel_coords = 0.0; - pixel_coords.row(0) = x; - pixel_coords.row(1) = y; - casacore::Vector failures; - if (!_coord_sys->toWorldMany(world_coords, pixel_coords, failures)) { - return false; - } - - // Save x and y values as Quantities - casacore::Vector world_units = _coord_sys->worldAxisUnits(); - casacore::Vector x_wcs = world_coords.row(0); - casacore::Vector y_wcs = world_coords.row(1); - // reference points: corners (x0, y0, x1, y1, x2, y2, x3, y3) in world coordinates - wcs_points.resize(num_points * 2); - for (int i = 0; i < num_points; ++i) { - wcs_points[i * 2] = casacore::Quantity(x_wcs(i), world_units(0)); - wcs_points[(i * 2) + 1] = casacore::Quantity(y_wcs(i), world_units(1)); - } - return true; -} - -void Region::RectanglePointsToCorners( - std::vector& pixel_points, float rotation, casacore::Vector& x, casacore::Vector& y) { - // Convert rectangle control points to 4 corner points - float center_x(pixel_points[0].x()), center_y(pixel_points[0].y()); - float width(pixel_points[1].x()), height(pixel_points[1].y()); - x.resize(4); - y.resize(4); - - if (rotation == 0.0) { - float x_min(center_x - width / 2.0f), x_max(center_x + width / 2.0f); - float y_min(center_y - height / 2.0f), y_max(center_y + height / 2.0f); - // Bottom left - x(0) = x_min; - y(0) = y_min; - // Bottom right - x(1) = x_max; - y(1) = y_min; - // Top right - x(2) = x_max; - y(2) = y_max; - // Top left - x(3) = x_min; - y(3) = y_max; - } else { - // Apply rotation matrix to get width and height vectors in rotated basis - float cos_x = cos(rotation * M_PI / 180.0f); - float sin_x = sin(rotation * M_PI / 180.0f); - float width_vector_x = cos_x * width; - float width_vector_y = sin_x * width; - float height_vector_x = -sin_x * height; - float height_vector_y = cos_x * height; - - // Bottom left - x(0) = center_x + (-width_vector_x - height_vector_x) / 2.0f; - y(0) = center_y + (-width_vector_y - height_vector_y) / 2.0f; - // Bottom right - x(1) = center_x + (width_vector_x - height_vector_x) / 2.0f; - y(1) = center_y + (width_vector_y - height_vector_y) / 2.0f; - // Top right - x(2) = center_x + (width_vector_x + height_vector_x) / 2.0f; - y(2) = center_y + (width_vector_y + height_vector_y) / 2.0f; - // Top left - x(3) = center_x + (-width_vector_x + height_vector_x) / 2.0f; - y(3) = center_y + (-width_vector_y + height_vector_y) / 2.0f; - } -} - -bool Region::EllipsePointsToWorld( - std::vector& pixel_points, std::vector& wcs_points, float& ellipse_rotation) { - // Convert CARTA ellipse points (cx, cy), (bmaj, bmin) to world coordinates - if (pixel_points.size() != 2) { - return false; - } - - std::vector center_world; - if (!ConvertCartaPointToWorld(pixel_points[0], center_world)) { - return false; - } - - // Store center point quantities - wcs_points = center_world; - - // Convert bmaj, bmin from pixel length to world length - float bmaj(pixel_points[1].x()), bmin(pixel_points[1].y()); - casacore::Quantity bmaj_world = _coord_sys->toWorldLength(bmaj, 0); - casacore::Quantity bmin_world = _coord_sys->toWorldLength(bmin, 1); - ellipse_rotation = GetRegionState().rotation; - - // Check if bmaj/bmin units conform (false for pV image: arcsec and Hz) - if (!bmaj_world.isConform(bmin_world.getUnit())) { - return false; - } - - // bmaj > bmin (world coords) required for WCEllipsoid; adjust rotation angle - if (bmaj_world > bmin_world) { - // carta rotation is from y-axis, ellipse rotation is from x-axis - ellipse_rotation += 90.0; - } else { - // swapping takes care of 90 deg adjustment - std::swap(bmaj_world, bmin_world); - } - - wcs_points.push_back(bmaj_world); - wcs_points.push_back(bmin_world); - return true; +std::shared_mutex& Region::GetActiveTaskMutex() { + return _active_task_mutex; } // ************************************************************************* -// Apply region to any image +// Apply region to image and return LCRegion, mask or Record -std::shared_ptr Region::GetImageRegion(int file_id, std::shared_ptr output_csys, - const casacore::IPosition& output_shape, const StokesSource& stokes_source, bool report_error) { - // Apply region to non-reference image as converted polygon vertices - // Will return nullptr if outside image or is not a closed LCRegion (line or polyline) - std::shared_ptr lc_region; - if (IsAnnotation()) { - return lc_region; - } +std::shared_ptr Region::GetImageRegion(int file_id, std::shared_ptr csys, + const casacore::IPosition& image_shape, const StokesSource& stokes_source, bool report_error) { + // Return lattice-coordinate region applied to image and/or computed stokes. + // Returns nullptr if is annotation, is not a closed region (line/polyline), or outside image. + std::shared_ptr lcregion; - std::lock_guard guard(_region_approx_mutex); - // The cache of lattice coordinate region is only for the original image (not computed stokes image). In order to avoid the ambiguity - if (stokes_source.IsOriginalImage()) { - lc_region = GetCachedLCRegion(file_id); + if (IsAnnotation() || IsLineType()) { + return lcregion; } - if (!lc_region) { - auto region_state = GetRegionState(); - if (file_id == region_state.reference_file_id) { - // Convert reference WCRegion to LCRegion and cache it - lc_region = GetConvertedLCRegion(file_id, output_csys, output_shape, stokes_source, report_error); - } else { - bool use_polygon = UseApproximatePolygon(output_csys); // check region distortion - - if (!use_polygon) { - // No distortion, do direct region conversion if possible (unless outside image or rotbox) - lc_region = GetConvertedLCRegion(file_id, output_csys, output_shape, stokes_source, report_error); - } - - if (lc_region) { - spdlog::debug("Using direct region conversion for {}", RegionName(region_state.type)); - } else { - // Use polygon approximation of reference region to translate to another image - spdlog::debug("Using polygon approximation for matched {} region", RegionName(region_state.type)); - lc_region = GetAppliedPolygonRegion(file_id, output_csys, output_shape); + // Check cache + lcregion = GetCachedLCRegion(file_id, stokes_source); - // Cache converted polygon - // Only for the original image (not computed stokes image). In order to avoid the ambiguity - if (lc_region && stokes_source.IsOriginalImage()) { - _polygon_regions[file_id] = lc_region; + if (!lcregion) { + if (IsInReferenceImage(file_id)) { + if (!_lcregion_set) { + // Create LCRegion from TableRecord + casacore::TableRecord region_record; + if (GetRegionState().IsRotbox()) { + // ControlPointsRecord is unrotated box for export, apply rotation for LCRegion + region_record = GetRotboxRecordForLCRegion(image_shape); + } else { + region_record = GetControlPointsRecord(image_shape); } - } - } - } - - return lc_region; -} - -bool Region::UseApproximatePolygon(std::shared_ptr output_csys) { - // Determine whether to convert region directly, or approximate it as a polygon in the output image. - bool use_polygon(true); - auto region_state = GetRegionState(); - CARTA::RegionType region_type = region_state.type; - - // Check ellipse and rectangle distortion; always true for other regions (polygon and point) - if ((region_type == CARTA::RegionType::ELLIPSE) || (region_type == CARTA::RegionType::RECTANGLE)) { - CARTA::Point center_point = region_state.control_points[0]; - double x_length(region_state.control_points[1].x()), y_length(region_state.control_points[1].y()); - - // Ratio of vector lengths in reference image region - double ref_length_ratio; - if (region_type == CARTA::RegionType::ELLIPSE) { - ref_length_ratio = x_length / y_length; - } else { - ref_length_ratio = y_length / x_length; - } - - std::vector points; - if (region_type == CARTA::RegionType::ELLIPSE) { - // Make polygon with 4 points - points = GetApproximateEllipsePoints(4); - } else { - // Get midpoints of 4 sides of rectangle - points = GetRectangleMidpoints(); - } - - // add center point; points = [p0, p1, p2, p3, center] - points.push_back(center_point); - - // convert points to output image pixels - casacore::Vector x, y; - if (ConvertPointsToImagePixels(points, output_csys, x, y)) { - // vector0 is (center, p0), vector1 is (center, p1) - auto v0_delta_x = x[0] - x[4]; - auto v0_delta_y = y[0] - y[4]; - auto v1_delta_x = x[1] - x[4]; - auto v1_delta_y = y[1] - y[4]; - - // Ratio of vector lengths in converted region - auto v0_length = sqrt((v0_delta_x * v0_delta_x) + (v0_delta_y * v0_delta_y)); - auto v1_length = sqrt((v1_delta_x * v1_delta_x) + (v1_delta_y * v1_delta_y)); - double converted_length_ratio = v1_length / v0_length; - // Compare reference to converted length ratio - double length_ratio_difference = fabs(ref_length_ratio - converted_length_ratio); - // spdlog::debug("{} distortion check: length ratio difference={:.3e}", RegionName(region_type), length_ratio_difference); - - if (length_ratio_difference < 1e-4) { - // Passed ratio check; check dot product of converted region - double converted_dot_product = (v0_delta_x * v1_delta_x) + (v0_delta_y * v1_delta_y); - // spdlog::debug("{} distortion check: dot product={:.3e}", RegionName(region_type), converted_dot_product); - - if (fabs(converted_dot_product) < 1e-2) { - // passed distortion tests, do not use polygon approximation - use_polygon = false; + try { + lcregion.reset(casacore::LCRegion::fromRecord(region_record, "")); + } catch (const casacore::AipsError& err) { + // Region is outside image } - } - } - } - - return use_polygon; -} - -std::vector Region::GetRectangleMidpoints() { - // Return midpoints of 4 sides of rectangle - // Find corners with rotation: blc, brc, trc, tlc - auto region_state = GetRegionState(); - casacore::Vector x, y; - RectanglePointsToCorners(region_state.control_points, region_state.rotation, x, y); - - std::vector midpoints; - // start with right side brc, trc - midpoints.push_back(Message::Point((x[1] + x[2]) / 2.0, (y[1] + y[2]) / 2.0)); - midpoints.push_back(Message::Point((x[2] + x[3]) / 2.0, (y[2] + y[3]) / 2.0)); - midpoints.push_back(Message::Point((x[3] + x[0]) / 2.0, (y[3] + y[0]) / 2.0)); - midpoints.push_back(Message::Point((x[0] + x[1]) / 2.0, (y[0] + y[1]) / 2.0)); + _lcregion_set = true; - return midpoints; -} - -std::shared_ptr Region::GetCachedPolygonRegion(int file_id) { - // Return cached polygon region applied to image with file_id - if (_polygon_regions.count(file_id)) { - std::lock_guard guard(_region_mutex); - return _polygon_regions.at(file_id); - } - - return std::shared_ptr(); -} - -std::shared_ptr Region::GetAppliedPolygonRegion( - int file_id, std::shared_ptr output_csys, const casacore::IPosition& output_shape) { - // Approximate region as polygon pixel vertices, and convert to given csys - std::shared_ptr lc_region; - - auto region_state = GetRegionState(); - bool is_point(region_state.type == CARTA::RegionType::POINT); - size_t nvertices(is_point ? 1 : DEFAULT_VERTEX_COUNT); - - // Set reference region as polygon vertices - auto polygon_points = GetReferencePolygonPoints(nvertices); - if (polygon_points.empty()) { - return lc_region; - } - - // Set x and y for matched polygon region - casacore::Vector x, y; - - if (polygon_points.size() == 1) { - // Point and ellipse have one vector for all points - if (!ConvertPointsToImagePixels(polygon_points[0], output_csys, x, y)) { - spdlog::error("Error approximating {} as polygon in matched image.", RegionName(region_state.type)); - return lc_region; - } - - if (!is_point) { - RemoveHorizontalPolygonPoints(x, y); - } - } else { - // Rectangle and polygon have one vector for each side of rectangle or segment of original polygon - for (auto& segment : polygon_points) { - casacore::Vector segment_x, segment_y; - if (!ConvertPointsToImagePixels(segment, output_csys, segment_x, segment_y)) { - spdlog::error("Error approximating {} as polygon in matched image.", RegionName(region_state.type)); - return lc_region; - } - - // For each polygon segment, if horizontal then remove points to fix LCPolygon mask - RemoveHorizontalPolygonPoints(segment_x, segment_y); - - auto old_size = x.size(); - x.resize(old_size + segment_x.size(), true); - y.resize(old_size + segment_y.size(), true); - - // Append selected segment points - for (auto i = 0; i < segment_x.size(); ++i) { - x[old_size + i] = segment_x[i]; - y[old_size + i] = segment_y[i]; - } - } - } - - try { - if (is_point) { - // Point is not a polygon (needs at least 3 points), use LCBox instead - // Form blc, trc - size_t ndim(output_shape.size()); - casacore::Vector blc(ndim, 0.0), trc(ndim); - blc(0) = x(0); - blc(1) = y(0); - trc(0) = x(0); - trc(1) = y(0); - - for (size_t i = 2; i < ndim; ++i) { - trc(i) = output_shape(i) - 1; + // Cache LCRegion + if (lcregion && stokes_source.IsOriginalImage()) { + std::lock_guard guard(_lcregion_mutex); + _lcregion = lcregion; + } } - - lc_region.reset(new casacore::LCBox(blc, trc, output_shape)); } else { - // Need 2-dim shape - casacore::IPosition keep_axes(2, 0, 1); - casacore::IPosition region_shape(output_shape.keepAxes(keep_axes)); - lc_region.reset(new casacore::LCPolygon(x, y, region_shape)); + if (!_region_converter) { + _region_converter.reset(new RegionConverter(GetRegionState(), _coord_sys)); + } + return _region_converter->GetImageRegion(file_id, csys, image_shape, stokes_source, report_error); } - } catch (const casacore::AipsError& err) { - spdlog::error("Cannot apply {} to file {}: {}", RegionName(region_state.type), file_id, err.getMesg()); } - return lc_region; + return lcregion; } -std::vector> Region::GetReferencePolygonPoints(int num_vertices) { - // Approximates reference region as polygon with input number of vertices. - // Sets _polygon_control_points in reference image pixel coordinates. - // Returns points as long as region type is supported and a closed region. - auto region_state = GetRegionState(); - switch (region_state.type) { - case CARTA::POINT: { - std::vector> points; - points.push_back(region_state.control_points); - return points; - } - case CARTA::RECTANGLE: - case CARTA::POLYGON: { - return GetApproximatePolygonPoints(num_vertices); - } - case CARTA::ELLIPSE: { - std::vector> points; - points.push_back(GetApproximateEllipsePoints(num_vertices)); - return points; - } - default: - return {}; +std::shared_ptr Region::GetCachedLCRegion(int file_id, const StokesSource& stokes_source) { + // Return cached LCRegion applied to image + std::shared_ptr lcregion; + if (!stokes_source.IsOriginalImage()) { + return lcregion; } -} - -std::vector> Region::GetApproximatePolygonPoints(int num_vertices) { - // Approximate RECTANGLE or POLYGON region as polygon with num_vertices. - // Returns vector of points for each segment - std::vector> polygon_points; - auto region_state = GetRegionState(); - std::vector region_points; - CARTA::RegionType region_type(region_state.type); - - if (region_type == CARTA::RegionType::RECTANGLE) { - // convert control points to corners to create 4-point polygon - casacore::Vector x, y; - RectanglePointsToCorners(region_state.control_points, region_state.rotation, x, y); - - for (size_t i = 0; i < x.size(); ++i) { - region_points.push_back(Message::Point(x(i), y(i))); - } - } else if (region_type == CARTA::RegionType::POLYGON) { - region_points = region_state.control_points; + if (IsInReferenceImage(file_id)) { + // Return _lcregion even if unassigned + std::lock_guard guard(_lcregion_mutex); + return _lcregion; } else { - spdlog::error("Error approximating {} as polygon: region type not supported", RegionName(region_state.type)); - return polygon_points; - } - - // Close polygon - CARTA::Point first_point(region_points[0]); - region_points.push_back(first_point); - - double total_length = GetTotalSegmentLength(region_points); - double target_segment_length = total_length / num_vertices; - - // Divide each region polygon segment into target number of segments with target length - for (size_t i = 1; i < region_points.size(); ++i) { - // Handle segment from point[i-1] to point[i] - std::vector segment_points; - - auto delta_x = region_points[i].x() - region_points[i - 1].x(); - auto delta_y = region_points[i].y() - region_points[i - 1].y(); - auto segment_length = sqrt((delta_x * delta_x) + (delta_y * delta_y)); - auto dir_x = delta_x / segment_length; - auto dir_y = delta_y / segment_length; - auto target_nsegment = round(segment_length / target_segment_length); - auto target_length = segment_length / target_nsegment; - - auto first_segment_point(region_points[i - 1]); - segment_points.push_back(first_segment_point); - - auto first_x(first_segment_point.x()); - auto first_y(first_segment_point.y()); - - for (size_t j = 1; j < target_nsegment; ++j) { - auto length_from_first = j * target_length; - auto x_offset = dir_x * length_from_first; - auto y_offset = dir_y * length_from_first; - segment_points.push_back(Message::Point(first_x + x_offset, first_y + y_offset)); - } - - polygon_points.push_back(segment_points); - } - - return polygon_points; -} - -std::vector Region::GetApproximateEllipsePoints(int num_vertices) { - // Approximate ELLIPSE region as polygon with num_vertices, return points - std::vector polygon_points; - auto region_state = GetRegionState(); - - auto cx = region_state.control_points[0].x(); - auto cy = region_state.control_points[0].y(); - auto bmaj = region_state.control_points[1].x(); - auto bmin = region_state.control_points[1].y(); - - auto delta_theta = 2.0 * M_PI / num_vertices; - auto rotation = region_state.rotation * M_PI / 180.0; - auto cos_rotation = cos(rotation); - auto sin_rotation = sin(rotation); - - for (int i = 0; i < num_vertices; ++i) { - auto theta = i * delta_theta; - auto rot_bmin = bmin * cos(theta); - auto rot_bmaj = bmaj * sin(theta); - - auto x_offset = (cos_rotation * rot_bmin) - (sin_rotation * rot_bmaj); - auto y_offset = (sin_rotation * rot_bmin) + (cos_rotation * rot_bmaj); - - polygon_points.push_back(Message::Point(cx + x_offset, cy + y_offset)); - } - - return polygon_points; -} - -double Region::GetTotalSegmentLength(std::vector& points) { - // Accumulate length of each point-to-point segment; returns total length. - double total_length(0.0); - - for (size_t i = 1; i < points.size(); ++i) { - auto delta_x = points[i].x() - points[i - 1].x(); - auto delta_y = points[i].y() - points[i - 1].y(); - total_length += sqrt((delta_x * delta_x) + (delta_y * delta_y)); - } - - return total_length; -} - -bool Region::ConvertPointsToImagePixels(const std::vector& points, std::shared_ptr output_csys, - casacore::Vector& x, casacore::Vector& y) { - // Convert pixel coords in reference image (points) to pixel coords in output image - // Coordinates returned in x and y vectors - bool converted(true); - - try { - // Convert each pixel point to output csys pixel point - size_t npoints(points.size()); - x.resize(npoints); - y.resize(npoints); - for (auto i = 0; i < npoints; ++i) { - // Convert pixel to world in reference csys - std::vector world_point; // [x, y] - if (ConvertCartaPointToWorld(points[i], world_point)) { - // Convert reference world to output csys pixel - casacore::Vector pixel_point; // [x, y] - if (ConvertWorldToPixel(world_point, output_csys, pixel_point)) { - x(i) = pixel_point(0); - y(i) = pixel_point(1); - } else { // world to pixel failed - spdlog::error("Error converting region to output image pixel coords."); - converted = false; - break; - } - } else { // pixel to world failed - spdlog::error("Error converting region to reference image world coords."); - converted = false; - break; - } + if (!_region_converter) { + _region_converter.reset(new RegionConverter(GetRegionState(), _coord_sys)); } - } catch (const casacore::AipsError& err) { - spdlog::error("Error converting region to output image: {}", err.getMesg()); - converted = false; + return _region_converter->GetCachedLCRegion(file_id); } - return converted; + // Not cached + return lcregion; } casacore::ArrayLattice Region::GetImageRegionMask(int file_id) { - // Return pixel mask for this region; requires that lcregion for this file id has been set. - // Otherwise mask is empty array. + // Return pixel mask for region applied to image. + // Requires that LCRegion for this file id has been set and cached (via GetImageRegion), + // else returns empty mask. casacore::ArrayLattice mask; - if (IsAnnotation()) { - return mask; - } - - std::shared_ptr lcregion; - if ((file_id == GetRegionState().reference_file_id) && _applied_regions.count(file_id)) { - if (_applied_regions.at(file_id)) { - std::lock_guard guard(_region_mutex); - lcregion = _applied_regions.at(file_id); - } - } else if (_polygon_regions.count(file_id)) { - if (_polygon_regions.at(file_id)) { - std::lock_guard guard(_region_mutex); - lcregion = _polygon_regions.at(file_id); - } - } + auto stokes_source = StokesSource(); + auto lcregion = GetCachedLCRegion(file_id, stokes_source); if (lcregion) { - std::lock_guard guard(_region_mutex); - // Region can either be an extension region or a fixed region, depending on whether image is matched or not + // LCRegion is an extension region or a fixed region, depending on whether image is reference or matched. auto extended_region = dynamic_cast(lcregion.get()); if (extended_region) { auto& fixed_region = static_cast(extended_region->region()); @@ -793,188 +229,75 @@ casacore::ArrayLattice Region::GetImageRegionMask(int file_id) { return mask; } -// *************************************************************** -// Apply region to any image and return LCRegion Record for export - casacore::TableRecord Region::GetImageRegionRecord( - int file_id, std::shared_ptr output_csys, const casacore::IPosition& output_shape) { - // Return Record describing Region applied to output coord sys and image_shape in pixel coordinates + int file_id, std::shared_ptr csys, const casacore::IPosition& image_shape) { + // Return Record describing any region type applied to image, in pixel coordinates. + // Can be a line or annotation region. Rotated box is a polygon describing corners *without rotation*, for export. casacore::TableRecord record; - if (IsLineType()) { - record = GetLineRecord(file_id, output_csys, output_shape); + if (IsInReferenceImage(file_id)) { + record = GetControlPointsRecord(image_shape); } else { - // Get converted LCRegion (only for enclosed regions) - // Check applied regions cache - std::shared_ptr lc_region = GetCachedLCRegion(file_id); - - if (!lc_region) { - // Convert reference region to output image - lc_region = GetConvertedLCRegion(file_id, output_csys, output_shape); - } - - if (lc_region) { - // Get LCRegion definition in Record - record = lc_region->toRecord("region"); - if (record.isDefined("region")) { - record = record.asRecord("region"); - } + if (!_region_converter) { + _region_converter.reset(new RegionConverter(GetRegionState(), _coord_sys)); } + record = _region_converter->GetImageRegionRecord(file_id, csys, image_shape); } - if (record.empty()) { - // LCRegion failed, is outside the image or a rotated rectangle. - // Manually convert control points and put in Record. - record = GetRegionPointsRecord(file_id, output_csys, output_shape); + if (!record.isDefined("isRegion")) { + // Record created from control points instead of LCRegion::toRecord + CompleteRegionRecord(record, image_shape); } - return record; } -std::shared_ptr Region::GetCachedLCRegion(int file_id) { - // Return cached region applied to image with file_id - if (_applied_regions.count(file_id)) { - std::lock_guard ulock(_region_mutex); - return _applied_regions.at(file_id); - } - - return std::shared_ptr(); -} - -std::shared_ptr Region::GetConvertedLCRegion(int file_id, std::shared_ptr output_csys, - const casacore::IPosition& output_shape, const StokesSource& stokes_source, bool report_error) { - // Convert 2D reference WCRegion to LCRegion in output coord_sys and shape - std::shared_ptr lc_region; - bool is_reference_image(file_id == GetRegionState().reference_file_id); - - if (!is_reference_image && IsRotbox()) { - // Cannot convert rotbox region, it is a polygon type. - return lc_region; - } - - try { - // Convert reference WCRegion to LCRegion in image with output csys and shape - // Create reference WCRegion if not set - if (!_reference_region_set) { - SetReferenceRegion(); - } - - // Convert to LCRegion in output csys and shape - std::lock_guard guard(_region_mutex); - if (ReferenceRegionValid()) { - std::shared_ptr reference_region = std::atomic_load(&_reference_region); - lc_region.reset(reference_region->toLCRegion(*output_csys.get(), output_shape)); - } else if (is_reference_image) { - // Create LCRegion with control points for reference image - casacore::TableRecord region_record = GetControlPointsRecord(output_shape); - lc_region.reset(casacore::LCRegion::fromRecord(region_record, "")); - } - } catch (const casacore::AipsError& err) { - if (report_error) { - spdlog::error("Error converting {} to file {}: {}", RegionName(GetRegionState().type), file_id, err.getMesg()); - } - } - - // Cache the lattice coordinate region only for the original image (not computed stokes image). In order to avoid the ambiguity - if (lc_region && stokes_source.IsOriginalImage()) { - // Make a copy and cache LCRegion in map - std::lock_guard guard(_region_mutex); - _applied_regions[file_id] = lc_region; - } - - return lc_region; -} +// *************************************************************** -casacore::TableRecord Region::GetRegionPointsRecord( - int file_id, std::shared_ptr output_csys, const casacore::IPosition& output_shape) { - // Convert control points to output coord sys if needed, and return completed record. +casacore::TableRecord Region::GetControlPointsRecord(const casacore::IPosition& image_shape) { + // Return region Record in pixel coords in format of LCRegion::toRecord() from control points. + // Rotated box is returned as unrotated LCBox, rotation retrieved from RegionState. + // For rotated box for analytics, use GetRotboxRecordForLCRegion() as polygon. casacore::TableRecord record; - auto region_state = GetRegionState(); - - if (file_id == region_state.reference_file_id) { - record = GetControlPointsRecord(output_shape); - } else { - if (_wcs_control_points.empty()) { - SetReferenceRegion(); - } - - switch (region_state.type) { - case CARTA::RegionType::POINT: - record = GetPointRecord(output_csys, output_shape); - break; - case CARTA::RegionType::RECTANGLE: - case CARTA::RegionType::POLYGON: { - // Rectangle is LCPolygon with 4 corners; - // Rotbox requires special handling to get (unrotated) rectangle corners instead of rotated polygon vertices - record = (IsRotbox() ? GetRotboxRecord(output_csys) : GetPolygonRecord(output_csys)); - break; - } - case CARTA::RegionType::ELLIPSE: - record = GetEllipseRecord(output_csys); - break; - default: - break; - } - - // Add common record fields - CompleteLCRegionRecord(record, output_shape); - } - return record; -} - -casacore::TableRecord Region::GetControlPointsRecord(const casacore::IPosition& shape) { - // Return region Record in pixel coords in format of LCRegion::toRecord() for reference image (no conversion) - // Shape needed for LCBox Record for point region - casacore::TableRecord record; auto region_state = GetRegionState(); auto region_type = region_state.type; switch (region_type) { case CARTA::RegionType::POINT: case CARTA::RegionType::ANNPOINT: { - // Box with blc=trc - auto ndim = shape.size(); + // Box with blc=trc, same ndim as image and all chan/stokes + auto ndim = image_shape.size(); casacore::Vector blc(ndim, 0.0), trc(ndim, 0.0); - blc(0) = region_state.control_points[0].x(); blc(1) = region_state.control_points[0].y(); trc(0) = region_state.control_points[0].x(); trc(1) = region_state.control_points[0].y(); + for (size_t i = 2; i < ndim; ++i) { + trc(i) = image_shape(i) - 1.0f; + } + record.define("name", "LCBox"); record.define("blc", blc); record.define("trc", trc); break; } case CARTA::RegionType::RECTANGLE: - case CARTA::RegionType::ANNRECTANGLE: { + case CARTA::RegionType::ANNRECTANGLE: + case CARTA::RegionType::ANNTEXT: { // Rectangle is LCPolygon with 4 corners; calculate from center and width/height - casacore::Vector x(5), y(5); - float center_x = region_state.control_points[0].x(); - float center_y = region_state.control_points[0].y(); - float width = region_state.control_points[1].x(); - float height = region_state.control_points[1].y(); - float x_min(center_x - width / 2.0f), x_max(center_x + width / 2.0f); - float y_min(center_y - height / 2.0f), y_max(center_y + height / 2.0f); - // Bottom left - x(0) = x_min; - y(0) = y_min; - // Bottom right - x(1) = x_max; - y(1) = y_min; - // Top right - x(2) = x_max; - y(2) = y_max; - // Top left - x(3) = x_min; - y(3) = y_max; - // LCPolygon::toRecord includes last point as first point to close region - x(4) = x(0); - y(4) = y(0); + casacore::Vector x, y; + bool apply_rotation(false); + if (GetRegionState().GetRectangleCorners(x, y, apply_rotation)) { + // LCPolygon::toRecord includes last point as first point to close region + x.resize(5, true); + y.resize(5, true); + x(4) = x(0); + y(4) = y(0); - record.define("name", "LCPolygon"); - record.define("x", x); - record.define("y", y); + record.define("name", "LCPolygon"); + record.define("x", x); + record.define("y", y); + } break; } case CARTA::RegionType::LINE: @@ -1001,316 +324,86 @@ casacore::TableRecord Region::GetControlPointsRecord(const casacore::IPosition& y(npoints) = region_state.control_points[0].y(); record.define("name", "LCPolygon"); - } else if (region_type == CARTA::RegionType::LINE || region_type == CARTA::RegionType::ANNLINE) { - // CARTA USE ONLY, name not implemented in casacore - record.define("name", "Line"); - } else if (region_type == CARTA::RegionType::ANNVECTOR) { - // CARTA USE ONLY, name not implemented in casacore - record.define("name", "Vector"); - } else if (region_type == CARTA::RegionType::ANNRULER) { - // CARTA USE ONLY, name not implemented in casacore - record.define("name", "Ruler"); } else { - // CARTA USE ONLY, name not implemented in casacore - record.define("name", "Polyline"); + // CARTA USE ONLY, line region names not in casacore because not an LCRegion + record.define("name", _region_state.GetLineRegionName()); } record.define("x", x); record.define("y", y); break; } case CARTA::RegionType::ELLIPSE: - case CARTA::RegionType::ANNELLIPSE: { + case CARTA::RegionType::ANNELLIPSE: + case CARTA::RegionType::ANNCOMPASS: { casacore::Vector center(2), radii(2); center(0) = region_state.control_points[0].x(); center(1) = region_state.control_points[0].y(); - radii(0) = region_state.control_points[1].x(); - radii(1) = region_state.control_points[1].y(); + auto major = region_state.control_points[1].x(); + auto minor = region_state.control_points[1].y(); + auto ellipse_rotation = region_state.rotation; + + // Enforce major > minor (as in WCEllipsoid) for ellipse. + bool is_compass(region_type == CARTA::RegionType::ANNCOMPASS); + if (major > minor || is_compass) { + radii(0) = major; + radii(1) = minor; + // carta rotation is from y-axis, ellipse rotation is from x-axis + ellipse_rotation += 90.0; + } else { + // swapping takes care of 90 deg adjustment + radii(0) = minor; + radii(1) = major; + } - record.define("name", "LCEllipsoid"); + if (region_type == CARTA::RegionType::ANNCOMPASS) { + record.define("name", "compass"); + } else { + record.define("name", "LCEllipsoid"); + } record.define("center", center); record.define("radii", radii); // LCEllipsoid measured from major (x) axis - casacore::Quantity theta = casacore::Quantity(region_state.rotation + 90.0, "deg"); + casacore::Quantity theta = casacore::Quantity(ellipse_rotation, "deg"); theta.convert("rad"); record.define("theta", theta.getValue()); break; } - default: + default: // Annulus not implemented break; } - - CompleteLCRegionRecord(record, shape); - return record; -} - -void Region::CompleteLCRegionRecord(casacore::TableRecord& record, const casacore::IPosition& shape) { - // Add common Record fields for record defining region - if (!record.empty()) { - record.define("isRegion", casacore::RegionType::LC); - record.define("comment", ""); - record.define("oneRel", false); - - casacore::Vector region_shape; - if (GetRegionState().type == CARTA::RegionType::POINT) { - region_shape = shape.asVector(); // LCBox uses entire image shape - } else { - region_shape.resize(2); - region_shape(0) = shape(0); - region_shape(1) = shape(1); - } - record.define("shape", region_shape); - } -} - -casacore::TableRecord Region::GetPointRecord( - std::shared_ptr output_csys, const casacore::IPosition& output_shape) { - // Return point applied to output_csys in format of LCBox::toRecord() - casacore::TableRecord record; - try { - // wcs control points is single point (x, y) - casacore::Vector pixel_point; - if (ConvertWorldToPixel(_wcs_control_points, output_csys, pixel_point)) { - auto ndim = output_shape.size(); - casacore::Vector blc(ndim), trc(ndim); - blc = 0.0; - blc(0) = pixel_point(0); - blc(1) = pixel_point(1); - trc(0) = pixel_point(0); - trc(1) = pixel_point(1); - - for (size_t i = 2; i < ndim; ++i) { - trc(i) = output_shape(i) - 1; - } - - record.define("name", "LCBox"); - record.define("blc", blc); - record.define("trc", trc); - } else { - spdlog::error("Error converting point to image."); - } - } catch (const casacore::AipsError& err) { - spdlog::error("Error converting point to image: {}", err.getMesg()); - } - + CompleteRegionRecord(record, image_shape); return record; } -casacore::TableRecord Region::GetPolygonRecord(std::shared_ptr output_csys) { - // Return region applied to output_csys in format of LCPolygon::toRecord() - // This is for POLYGON or RECTANGLE (points are four corners of box) +casacore::TableRecord Region::GetRotboxRecordForLCRegion(const casacore::IPosition& image_shape) { + // Convert rotated box corners to polygon and create LCPolygon-type record casacore::TableRecord record; - auto type = GetRegionState().type; - - try { - size_t npoints(_wcs_control_points.size() / 2); - casacore::Vector x(npoints), y(npoints); // Record fields - - // Convert each wcs control point to pixel coords in output csys - for (size_t i = 0; i < _wcs_control_points.size(); i += 2) { - std::vector world_point(2); - world_point[0] = _wcs_control_points[i]; - world_point[1] = _wcs_control_points[i + 1]; - casacore::Vector pixel_point; - if (ConvertWorldToPixel(world_point, output_csys, pixel_point)) { - // Add to x and y Vectors - int index(i / 2); - x(index) = pixel_point(0); - y(index) = pixel_point(1); - } else { - spdlog::error("Error converting {} to image pixels.", RegionName(type)); - return record; - } - } - - if (type == CARTA::RegionType::POLYGON) { - // LCPolygon::toRecord adds first point as last point to close region - x.resize(npoints + 1, true); - x(npoints) = x(0); - y.resize(npoints + 1, true); - y(npoints) = y(0); - } + casacore::Vector x, y; + if (GetRegionState().GetRectangleCorners(x, y)) { + // LCPolygon::toRecord includes last point as first point to close region + x.resize(5, true); + y.resize(5, true); + x(4) = x(0); + y(4) = y(0); - // Add fields for this region type record.define("name", "LCPolygon"); record.define("x", x); record.define("y", y); - } catch (const casacore::AipsError& err) { - spdlog::error("Error converting () to image: {}", RegionName(type), err.getMesg()); } - - return record; -} - -casacore::TableRecord Region::GetRotboxRecord(std::shared_ptr output_csys) { - // Determine corners of unrotated box (control points) applied to output_csys. - // Return region applied to output_csys in format of LCPolygon::toRecord() - casacore::TableRecord record; - - try { - // Get 4 corner points in pixel coordinates - std::vector pixel_points(GetRegionState().control_points); - float center_x(pixel_points[0].x()), center_y(pixel_points[0].y()); - float width(pixel_points[1].x()), height(pixel_points[1].y()); - - int num_points(4); - casacore::Vector x(num_points), y(num_points); - float x_min(center_x - width / 2.0f), x_max(center_x + width / 2.0f); - float y_min(center_y - height / 2.0f), y_max(center_y + height / 2.0f); - // Bottom left - x(0) = x_min; - y(0) = y_min; - // Bottom right - x(1) = x_max; - y(1) = y_min; - // Top right - x(2) = x_max; - y(2) = y_max; - // Top left - x(3) = x_min; - y(3) = y_max; - - // Convert corners to reference world coords - int num_axes(_coord_sys->nPixelAxes()); - casacore::Matrix pixel_coords(num_axes, num_points); - casacore::Matrix world_coords(num_axes, num_points); - pixel_coords = 0.0; - pixel_coords.row(0) = x; - pixel_coords.row(1) = y; - casacore::Vector failures; - if (!_coord_sys->toWorldMany(world_coords, pixel_coords, failures)) { - spdlog::error("Error converting rectangle pixel coordinates to world."); - return record; - } - casacore::Vector x_wcs = world_coords.row(0); - casacore::Vector y_wcs = world_coords.row(1); - - // Units for Quantities - casacore::Vector world_units = output_csys->worldAxisUnits(); - casacore::Vector corner_x(num_points), corner_y(num_points); - for (size_t i = 0; i < num_points; i++) { - // Convert x and y reference world coords to Quantity - std::vector world_point; - world_point.push_back(casacore::Quantity(x_wcs(i), world_units(0))); - world_point.push_back(casacore::Quantity(y_wcs(i), world_units(1))); - - // Convert reference world point to output pixel point - casacore::Vector pixel_point; - if (ConvertWorldToPixel(world_point, output_csys, pixel_point)) { - // Add to x and y Vectors - corner_x(i) = pixel_point(0); - corner_y(i) = pixel_point(1); - } else { - spdlog::error("Error converting rectangle coordinates to image."); - return record; - } - } - - // Add fields for this region type - record.define("name", "LCPolygon"); - record.define("x", corner_x); - record.define("y", corner_y); - } catch (const casacore::AipsError& err) { - spdlog::error("Error converting rectangle to image: {}", err.getMesg()); - } - - return record; -} - -casacore::TableRecord Region::GetEllipseRecord(std::shared_ptr output_csys) { - // Return region applied to output_csys in format of LCEllipsoid::toRecord() - casacore::TableRecord record; - - // Values to set in Record - casacore::Vector center(2), radii(2); - - // Center point - std::vector world_point(2); - world_point[0] = _wcs_control_points[0]; - world_point[1] = _wcs_control_points[1]; - casacore::Vector pixel_point; - try { - if (ConvertWorldToPixel(world_point, output_csys, pixel_point)) { - center(0) = pixel_point(0); - center(1) = pixel_point(1); - - // Convert radii to output world units, then to pixels - casacore::Vector increments(output_csys->increment()); - casacore::Vector world_units = output_csys->worldAxisUnits(); - casacore::Quantity bmaj = _wcs_control_points[2]; - bmaj.convert(world_units(0)); - casacore::Quantity bmin = _wcs_control_points[3]; - bmin.convert(world_units(1)); - radii(0) = fabs(bmaj.getValue() / increments(0)); - radii(1) = fabs(bmin.getValue() / increments(1)); - - // Add fields for this region type - record.define("name", "LCEllipsoid"); - record.define("center", center); - record.define("radii", radii); - - // LCEllipsoid measured from major (x) axis - // TODO: adjust angle for output csys - casacore::Quantity theta = casacore::Quantity(GetRegionState().rotation + 90.0, "deg"); - theta.convert("rad"); - record.define("theta", theta.getValue()); - } else { - spdlog::error("Incompatible coordinate systems for ellipse conversion."); - } - } catch (const casacore::AipsError& err) { - spdlog::error("Error converting ellipse to image: {}", err.getMesg()); - } - - return record; -} - -casacore::TableRecord Region::GetLineRecord( - int file_id, std::shared_ptr image_csys, const casacore::IPosition& image_shape) { - // Return line region (not closed LCRegion) in pixel coords - // for region applied to input image_csys with input image_shape - casacore::TableRecord record; - auto region_state = GetRegionState(); - - if (file_id == region_state.reference_file_id) { - record = GetControlPointsRecord(image_shape); - } else { - std::vector line_points = region_state.control_points; - auto region_type = region_state.type; - casacore::Vector x, y; - - if (ConvertPointsToImagePixels(line_points, image_csys, x, y)) { - // CARTA USE ONLY, names not implemented in casacore - if (region_type == CARTA::RegionType::LINE || region_type == CARTA::RegionType::ANNLINE) { - record.define("name", "Line"); - } else if (region_type == CARTA::RegionType::POLYLINE || region_type == CARTA::RegionType::ANNPOLYLINE) { - record.define("name", "Polyline"); - } else if (region_type == CARTA::RegionType::ANNVECTOR) { - record.define("name", "Vector"); - } else if (region_type == CARTA::RegionType::ANNRULER) { - record.define("name", "Ruler"); - } else { - return record; - } - - record.define("x", x); - record.define("y", y); - } - } - CompleteRegionRecord(record, image_shape); - return record; } void Region::CompleteRegionRecord(casacore::TableRecord& record, const casacore::IPosition& image_shape) { + // Add common Record fields for record defining region if (!record.empty()) { - // Complete Record with common fields - record.define("isRegion", 1); + record.define("isRegion", casacore::RegionType::LC); record.define("comment", ""); - record.define("oneRel", false); + record.define("oneRel", false); // control points are 0-based - casacore::Vector record_shape; - if (GetRegionState().type == CARTA::RegionType::POINT) { + casacore::Vector record_shape; + if (_region_state.IsPoint()) { // LCBox uses entire image shape record_shape = image_shape.asVector(); } else { @@ -1322,133 +415,3 @@ void Region::CompleteRegionRecord(casacore::TableRecord& record, const casacore: record.define("shape", record_shape); } } - -// *************************************************************** -// Conversion utilities - -bool Region::ConvertCartaPointToWorld(const CARTA::Point& point, std::vector& world_point) { - // Converts a CARTA point(x, y) in pixel coordinates to a Quantity vector [x, y] in world coordinates - // Returns whether conversion was successful - - // Vectors must be same number of axes as in coord system for conversion: - int naxes(_coord_sys->nPixelAxes()); - casacore::Vector pixel_values(naxes), world_values(naxes); - pixel_values = 0.0; // set "extra" axes to 0, not needed - pixel_values(0) = point.x(); - pixel_values(1) = point.y(); - - // convert pixel vector to world vector - if (!_coord_sys->toWorld(world_values, pixel_values)) { - return false; - } - - // Set Quantities from world values and units - casacore::Vector world_units = _coord_sys->worldAxisUnits(); - - world_point.resize(2); - world_point[0] = casacore::Quantity(world_values(0), world_units(0)); - world_point[1] = casacore::Quantity(world_values(1), world_units(1)); - return true; -} - -bool Region::ConvertWorldToPixel(std::vector& world_point, std::shared_ptr output_csys, - casacore::Vector& pixel_point) { - // Convert input reference world coord to output world coord, then to pixel coord - // Exception should be caught in calling function for creating error message - bool success(false); - - if (_coord_sys->hasDirectionCoordinate() && output_csys->hasDirectionCoordinate()) { - // Input and output direction reference frames - casacore::MDirection::Types reference_dir_type = _coord_sys->directionCoordinate().directionType(); - casacore::MDirection::Types output_dir_type = output_csys->directionCoordinate().directionType(); - - // Convert world point from reference to output coord sys - casacore::MDirection world_direction(world_point[0], world_point[1], reference_dir_type); - if (reference_dir_type != output_dir_type) { - world_direction = casacore::MDirection::Convert(world_direction, output_dir_type)(); - } - - // Convert output world point to pixel point - output_csys->directionCoordinate().toPixel(pixel_point, world_direction); - success = true; - } else if (_coord_sys->hasLinearCoordinate() && output_csys->hasLinearCoordinate()) { - // Get linear axes indices - auto indices = output_csys->linearAxesNumbers(); - if (indices.size() != 2) { - return false; - } - // Input and output linear frames - casacore::Vector output_units = output_csys->worldAxisUnits(); - casacore::Vector world_point_value(output_csys->nWorldAxes(), 0); - world_point_value(indices(0)) = world_point[0].get(output_units(indices(0))).getValue(); - world_point_value(indices(1)) = world_point[1].get(output_units(indices(1))).getValue(); - - // Convert world point to output pixel point - casacore::Vector tmp_pixel_point; - output_csys->toPixel(tmp_pixel_point, world_point_value); - - // Only fill the pixel coordinate results - pixel_point.resize(2); - pixel_point(0) = tmp_pixel_point(indices(0)); - pixel_point(1) = tmp_pixel_point(indices(1)); - success = true; - } - - return success; -} - -std::shared_mutex& Region::GetActiveTaskMutex() { - return _active_task_mutex; -} - -void Region::RemoveHorizontalPolygonPoints(casacore::Vector& x, casacore::Vector& y) { - // When polygon points have close y-points (horizontal segment), the x-range is masked only to the next point. - // Remove points not near integral pixel. - std::vector keep_x, keep_y; - size_t npoints(x.size()); - - for (int i = 0; i < npoints - 2; ++i) { - if (i == 0) { - // always include first point of segment - keep_x.push_back(x[i]); - keep_y.push_back(y[i]); - continue; - } - - float this_y = y[i]; - float next_y = y[i + 1]; - if (!ValuesNear(this_y, next_y)) { - // Line connecting points not ~horizontal - keep point - keep_x.push_back(x[i]); - keep_y.push_back(y[i]); - continue; - } - - // Line connecting points ~horizontal - keep point nearest integral pixel - int pixel_y = static_cast(this_y); - - if (!ValuesNear(this_y, float(pixel_y))) { - // Skip point not near pixel - continue; - } - - if ((static_cast(next_y) == pixel_y) && ((this_y - pixel_y) > (next_y - pixel_y))) { - // Skip point if next point nearer to pixel - continue; - } - - keep_x.push_back(x[i]); - keep_y.push_back(y[i]); - } - - if (keep_x.size() < npoints) { - // Set to new vector with points removed - x = casacore::Vector(keep_x); - y = casacore::Vector(keep_y); - } -} - -bool Region::ValuesNear(float val1, float val2) { - // near and nearAbs in casacore Math - return val1 == 0 || val2 == 0 ? casacore::nearAbs(val1, val2) : casacore::near(val1, val2); -} diff --git a/src/Region/Region.h b/src/Region/Region.h index bc041ff80..a0ede9615 100644 --- a/src/Region/Region.h +++ b/src/Region/Region.h @@ -4,7 +4,7 @@ SPDX-License-Identifier: GPL-3.0-or-later */ -//# Region.h: class for managing 2D region parameters +//# Region.h: class for managing 2D region parameters in reference image #ifndef CARTA_SRC_REGION_REGION_H_ #define CARTA_SRC_REGION_REGION_H_ @@ -15,121 +15,53 @@ #include #include -#include -#include #include #include #include -#include "Util/Message.h" +#include "RegionConverter.h" +#include "RegionState.h" #include "Util/Stokes.h" -#define DEFAULT_VERTEX_COUNT 1000 - namespace carta { -inline std::string RegionName(CARTA::RegionType type) { - std::unordered_map region_names = {{CARTA::POINT, "point"}, {CARTA::LINE, "line"}, - {CARTA::POLYLINE, "polyline"}, {CARTA::RECTANGLE, "rectangle"}, {CARTA::ELLIPSE, "ellipse"}, {CARTA::ANNULUS, "annulus"}, - {CARTA::POLYGON, "polygon"}, {CARTA::ANNPOINT, "ann point"}, {CARTA::ANNLINE, "ann line"}, {CARTA::ANNPOLYLINE, "ann polyline"}, - {CARTA::ANNRECTANGLE, "ann rectangle"}, {CARTA::ANNELLIPSE, "ann ellipse"}, {CARTA::ANNPOLYGON, "ann polygon"}, - {CARTA::ANNVECTOR, "vector"}, {CARTA::ANNRULER, "ruler"}, {CARTA::ANNTEXT, "text"}, {CARTA::ANNCOMPASS, "compass"}}; - return region_names[type]; -} - -struct RegionState { - // struct used for region parameters - int reference_file_id; - CARTA::RegionType type; - std::vector control_points; - float rotation; - - RegionState() : reference_file_id(-1), type(CARTA::POINT), rotation(0) {} - RegionState(int ref_file_id_, CARTA::RegionType type_, const std::vector& control_points_, float rotation_) - : reference_file_id(ref_file_id_), type(type_), control_points(control_points_), rotation(rotation_) {} - - void operator=(const RegionState& other) { - reference_file_id = other.reference_file_id; - type = other.type; - control_points = other.control_points; - rotation = other.rotation; - } - bool operator==(const RegionState& rhs) { - return (reference_file_id == rhs.reference_file_id) && (type == rhs.type) && !RegionChanged(rhs); - } - bool operator!=(const RegionState& rhs) { - return (reference_file_id != rhs.reference_file_id) || (type != rhs.type) || RegionChanged(rhs); - } - - bool RegionDefined() { - return !control_points.empty(); - } - - bool RegionChanged(const RegionState& rhs) { - // Ignores annotation params (for interrupting region calculations) - return (rotation != rhs.rotation) || PointsChanged(rhs); - } - bool PointsChanged(const RegionState& rhs) { - // Points must be same size, order, and value to be unchanged - if (control_points.size() != rhs.control_points.size()) { - return true; - } - for (int i = 0; i < control_points.size(); ++i) { - float x(control_points[i].x()), y(control_points[i].y()); - float rhs_x(rhs.control_points[i].x()), rhs_y(rhs.control_points[i].y()); - if ((x != rhs_x) || (y != rhs_y)) { - return true; - } - } - return false; - } -}; - class Region { public: Region(const RegionState& state, std::shared_ptr csys); - inline bool IsValid() { // control points validated + inline bool IsValid() { return _valid; }; - // set new region parameters + inline bool RegionChanged() { + return _region_changed; + } + + // Set new region parameters bool UpdateRegion(const RegionState& state); - // state accessors + // RegionState accessors, including region type information inline RegionState GetRegionState() { std::lock_guard guard(_region_state_mutex); RegionState region_state = _region_state; return region_state; } - inline int GetReferenceFileId() { - return GetRegionState().reference_file_id; - } - - inline bool RegionChanged() { // reference image, type, points, or rotation changed - return _region_changed; - } - inline bool IsPoint() { - return GetRegionState().type == CARTA::POINT; + return GetRegionState().IsPoint(); } inline bool IsLineType() { // Not enclosed region defined by 2 or more points - std::vector line_types{ - CARTA::LINE, CARTA::POLYLINE, CARTA::ANNLINE, CARTA::ANNPOLYLINE, CARTA::ANNVECTOR, CARTA::ANNRULER}; - auto type = GetRegionState().type; - return std::find(line_types.begin(), line_types.end(), type) != line_types.end(); + return GetRegionState().IsLineType(); } - inline bool IsRotbox() { - RegionState rs = GetRegionState(); - return ((rs.type == CARTA::RECTANGLE || rs.type == CARTA::ANNRECTANGLE) && (rs.rotation != 0.0)); + inline bool IsInReferenceImage(int file_id) { + return file_id == GetRegionState().reference_file_id; } inline bool IsAnnotation() { - return GetRegionState().type > CARTA::POLYGON; + return GetRegionState().IsAnnotation(); } inline std::shared_ptr CoordinateSystem() { @@ -139,108 +71,53 @@ class Region { // Communication bool IsConnected(); void WaitForTaskCancellation(); + std::shared_mutex& GetActiveTaskMutex(); - // Converted region as approximate LCPolygon and its mask - std::shared_ptr GetImageRegion(int file_id, std::shared_ptr image_csys, - const casacore::IPosition& image_shape, const StokesSource& stokes_source = StokesSource(), bool report_error = true); + // LCRegion and mask for region applied to image. Must be a closed region (not line) and not annotation. + std::shared_ptr GetImageRegion(int file_id, std::shared_ptr csys, + const casacore::IPosition& shape, const StokesSource& stokes_source = StokesSource(), bool report_error = true); casacore::ArrayLattice GetImageRegionMask(int file_id); - // Converted region in Record for export + // Record for region applied to image, for export. Not for converting to LCRegion for analytics. casacore::TableRecord GetImageRegionRecord( - int file_id, std::shared_ptr output_csys, const casacore::IPosition& output_shape); - - std::shared_mutex& GetActiveTaskMutex(); + int file_id, std::shared_ptr csys, const casacore::IPosition& shape); private: - bool SetPoints(const std::vector& points); - - // check points: number required for region type, and values are finite + // Check points: number required for region type, and values are finite bool CheckPoints(const std::vector& points, CARTA::RegionType type); bool PointsFinite(const std::vector& points); - // Reset cache when region changes + // Cached LCRegion void ResetRegionCache(); + std::shared_ptr GetCachedLCRegion(int file_id, const StokesSource& stokes_source); - // Check if reference region is set successfully - bool ReferenceRegionValid(); - - // Apply region to reference image, set WCRegion and wcs control points. - void SetReferenceRegion(); - bool RectanglePointsToWorld(std::vector& pixel_points, std::vector& wcs_points); - void RectanglePointsToCorners(std::vector& pixel_points, float rotation, casacore::Vector& x, - casacore::Vector& y); - bool EllipsePointsToWorld(std::vector& pixel_points, std::vector& wcs_points, float& rotation); - - // Reference region as approximate polygon converted to image coordinates; used for data streams - bool UseApproximatePolygon(std::shared_ptr output_csys); - std::vector GetRectangleMidpoints(); - std::shared_ptr GetCachedPolygonRegion(int file_id); - std::shared_ptr GetAppliedPolygonRegion( - int file_id, std::shared_ptr output_csys, const casacore::IPosition& output_shape); - std::vector> GetReferencePolygonPoints(int num_vertices); - std::vector> GetApproximatePolygonPoints(int num_vertices); - std::vector GetApproximateEllipsePoints(int num_vertices); - double GetTotalSegmentLength(std::vector& points); - bool ConvertPointsToImagePixels(const std::vector& points, std::shared_ptr output_csys, - casacore::Vector& x, casacore::Vector& y); - void RemoveHorizontalPolygonPoints(casacore::Vector& x, casacore::Vector& y); - bool ValuesNear(float val1, float val2); - - // Region applied to any image; used for export - std::shared_ptr GetCachedLCRegion(int file_id); - std::shared_ptr GetConvertedLCRegion(int file_id, std::shared_ptr output_csys, - const casacore::IPosition& output_shape, const StokesSource& stokes_source = StokesSource(), bool report_error = true); - - // Control points converted to pixel coords in output image, returned in LCRegion Record format for export - casacore::TableRecord GetRegionPointsRecord( - int file_id, std::shared_ptr output_csys, const casacore::IPosition& output_shape); + // Record in pixel coordinates from control points, for reference image casacore::TableRecord GetControlPointsRecord(const casacore::IPosition& shape); - void CompleteLCRegionRecord(casacore::TableRecord& record, const casacore::IPosition& shape); - casacore::TableRecord GetPointRecord(std::shared_ptr output_csys, const casacore::IPosition& output_shape); - casacore::TableRecord GetLineRecord( - int file_id, std::shared_ptr image_csys, const casacore::IPosition& image_shape); - casacore::TableRecord GetPolygonRecord(std::shared_ptr output_csys); - casacore::TableRecord GetRotboxRecord(std::shared_ptr output_csys); - casacore::TableRecord GetEllipseRecord(std::shared_ptr output_csys); + casacore::TableRecord GetRotboxRecordForLCRegion(const casacore::IPosition& shape); void CompleteRegionRecord(casacore::TableRecord& record, const casacore::IPosition& image_shape); - // Utilities to convert control points - // Input: CARTA::Point. Returns: point (x, y) in reference world coords - bool ConvertCartaPointToWorld(const CARTA::Point& point, std::vector& world_point); - // Input: point (x,y) in reference world coords. Returns: point (x,y) in output pixel coords - bool ConvertWorldToPixel(std::vector& world_point, std::shared_ptr output_csys, - casacore::Vector& pixel_point); - - // region parameters struct - RegionState _region_state; - - // coord sys and shape of reference image + // Coordinate system of reference image std::shared_ptr _coord_sys; - // Reference region cache - std::mutex _region_mutex; // creation of casacore regions is not threadsafe - std::mutex _region_approx_mutex; + // Region parameters + RegionState _region_state; std::mutex _region_state_mutex; - // Use a shared lock for long time calculations, use an exclusive lock for the object destruction - mutable std::shared_mutex _active_task_mutex; - - // Region cached as original type - std::shared_ptr _reference_region; // 2D region applied to reference image - std::vector _wcs_control_points; // for manual region conversion + // Region flags + bool _valid; // RegionState set properly + bool _region_changed; // control points or rotation changed - // Converted regions - // Reference region converted to image; key is file_id - std::unordered_map> _applied_regions; - // Polygon approximation region converted to image; key is file_id - std::unordered_map> _polygon_regions; + // Region applied to reference image + std::shared_ptr _lcregion; + std::mutex _lcregion_mutex; + bool _lcregion_set; // may be nullptr if outside image - // region flags - bool _valid; // RegionState set properly - bool _region_changed; // control points or rotation changed - bool _reference_region_set; // indicates attempt was made; may be null wcregion outside image + // Converter to handle region applied to matched image + std::unique_ptr _region_converter; - // Communication + // Communication: + // Use a shared lock for long time calculations, use an exclusive lock for the object destruction + mutable std::shared_mutex _active_task_mutex; volatile bool _connected = true; }; diff --git a/src/Region/RegionConverter.cc b/src/Region/RegionConverter.cc new file mode 100644 index 000000000..2e3771851 --- /dev/null +++ b/src/Region/RegionConverter.cc @@ -0,0 +1,1012 @@ +/* This file is part of the CARTA Image Viewer: https://github.com/CARTAvis/carta-backend + Copyright 2018- Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), + Associated Universities, Inc. (AUI) and the Inter-University Institute for Data Intensive Astronomy (IDIA) + SPDX-License-Identifier: GPL-3.0-or-later +*/ + +//# RegionConverter.cc: implementation of class for converting a region + +#include "RegionConverter.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "Logger/Logger.h" +#include "Util/Image.h" +#include "Util/Message.h" + +using namespace carta; + +RegionConverter::RegionConverter(const RegionState& state, std::shared_ptr csys) + : _region_state(state), _reference_coord_sys(csys), _reference_region_set(false) {} + +// **************************************************************************** +// Apply region to reference image in world coordinates (WCRegion) and save wcs control points + +bool RegionConverter::ReferenceRegionValid() { + return _reference_region_set && bool(_reference_region); +} + +void RegionConverter::SetReferenceWCRegion() { + // Create WCRegion (world coordinate region) in the reference image according to type using wcs control points. + // Sets _reference_region (including to nullptr if fails) and _reference_region_set (attempted). + // Supports closed (not line-type) annotation regions, for conversion to matched image LCRegion then Record for export. + // Rotated box is rotated, do not use for Record! + casacore::WCRegion* region(nullptr); + std::vector world_points; // for converting one point at a time + casacore::IPosition pixel_axes(2, 0, 1); // first and second axes only + casacore::Vector abs_rel; + auto type(_region_state.type); + + try { + switch (type) { + case CARTA::POINT: + case CARTA::ANNPOINT: { + // Convert one point for wcs control points + if (CartaPointToWorld(_region_state.control_points[0], _wcs_control_points)) { + // WCBox blc and trc are same point + casacore::Vector const wcs_points(_wcs_control_points); + region = new casacore::WCBox(wcs_points, wcs_points, pixel_axes, *_reference_coord_sys, abs_rel); + } + break; + } + case CARTA::RECTANGLE: // 4 corners + case CARTA::POLYGON: // vertices + case CARTA::ANNRECTANGLE: // 4 corners + case CARTA::ANNPOLYGON: // vertices + case CARTA::ANNTEXT: { // 4 corners of text box + // Use corners/vertices for wcs control points: x0, y0, x1, y1, etc. + if (type == CARTA::RECTANGLE || type == CARTA::ANNRECTANGLE || type == CARTA::ANNTEXT) { + if (!RectangleControlPointsToWorld(_wcs_control_points)) { + _wcs_control_points.clear(); + } + } else { + for (auto& point : _region_state.control_points) { + if (CartaPointToWorld(point, world_points)) { + _wcs_control_points.push_back(world_points[0]); + _wcs_control_points.push_back(world_points[1]); + } else { + _wcs_control_points.clear(); + break; + } + } + } + + if (!_wcs_control_points.empty()) { + // Convert from Vector (wcs control points) to Quantum for WCPolygon + casacore::Quantum> qx, qy; + // Separate x and y + std::vector x, y; + for (size_t i = 0; i < _wcs_control_points.size(); i += 2) { + x.push_back(_wcs_control_points[i].getValue()); + y.push_back(_wcs_control_points[i + 1].getValue()); + } + casacore::Vector x_v(x), y_v(y); + casacore::Vector world_units = _reference_coord_sys->worldAxisUnits(); + qx = x_v; // set values + qx.setUnit(world_units(0)); // set unit + qy = y_v; // set values + qy.setUnit(world_units(1)); // set unit + + region = new casacore::WCPolygon(qx, qy, pixel_axes, *_reference_coord_sys); + } + break; + } + case CARTA::ELLIPSE: // [(cx, cy), (bmaj, bmin)] + case CARTA::ANNELLIPSE: // [(cx, cy), (bmaj, bmin)] + case CARTA::ANNCOMPASS: { // [(cx, cy), (length, length)} + float ellipse_rotation; + if (EllipseControlPointsToWorld(_wcs_control_points, ellipse_rotation)) { + // wcs control points order: xcenter, ycenter, major axis, minor axis + casacore::Quantity theta(ellipse_rotation, "deg"); + theta.convert("rad"); + region = new casacore::WCEllipsoid(_wcs_control_points[0], _wcs_control_points[1], _wcs_control_points[2], + _wcs_control_points[3], theta, 0, 1, *_reference_coord_sys); + } + break; + } + default: // no WCRegion for line-type regions + break; + } + } catch (casacore::AipsError& err) { // region failed + spdlog::error("region type {} failed: {}", type, err.getMesg()); + } + + std::shared_ptr shared_region = std::shared_ptr(region); + std::atomic_store(&_reference_region, shared_region); + + // Flag indicates that attempt was made, to avoid repeated attempts + _reference_region_set = true; +} + +bool RegionConverter::RectangleControlPointsToWorld(std::vector& world_corners) { + // Convert CARTA rectangle points (cx, cy), (width, height) to corners in world coordinates (reference image) + // Get 4 corner points in pixel coordinates from control points, applying rotation + casacore::Vector x, y; + if (!_region_state.GetRectangleCorners(x, y)) { + return false; + } + + // Convert corners to wcs in one call for efficiency, rather than one point at a time + size_t num_points(x.size()), num_axes(_reference_coord_sys->nPixelAxes()); + casacore::Matrix pixel_coords(num_axes, num_points); + casacore::Matrix world_coords(num_axes, num_points); + pixel_coords = 0.0; + pixel_coords.row(0) = x; + pixel_coords.row(1) = y; + casacore::Vector failures; + if (!_reference_coord_sys->toWorldMany(world_coords, pixel_coords, failures)) { + return false; + } + + // Save x and y values as Quantities + casacore::Vector world_units = _reference_coord_sys->worldAxisUnits(); + casacore::Vector x_wcs = world_coords.row(0); + casacore::Vector y_wcs = world_coords.row(1); + // reference points: corners (x0, y0, x1, y1, x2, y2, x3, y3) in world coordinates + world_corners.resize(num_points * 2); + for (int i = 0; i < num_points; ++i) { + world_corners[i * 2] = casacore::Quantity(x_wcs(i), world_units(0)); + world_corners[(i * 2) + 1] = casacore::Quantity(y_wcs(i), world_units(1)); + } + return true; +} + +bool RegionConverter::EllipseControlPointsToWorld(std::vector& wcs_points, float& ellipse_rotation) { + // Convert CARTA ellipse points (cx, cy), (bmaj, bmin) to world coordinates, adjust rotation + auto pixel_points = _region_state.control_points; + ellipse_rotation = _region_state.rotation; + + // Convert center and store in wcs points + std::vector center_world; + if (!CartaPointToWorld(pixel_points[0], center_world)) { + return false; + } + wcs_points = center_world; + + // Convert bmaj, bmin from pixel length to world length + float bmaj(pixel_points[1].x()), bmin(pixel_points[1].y()); + casacore::Quantity bmaj_world = _reference_coord_sys->toWorldLength(bmaj, 0); + casacore::Quantity bmin_world = _reference_coord_sys->toWorldLength(bmin, 1); + + // Check if bmaj/bmin units conform (false for PV image, in arcsec and Hz) + if (!bmaj_world.isConform(bmin_world.getUnit())) { + return false; + } + + // bmaj > bmin (world coords) required for WCEllipsoid; adjust rotation angle + if (bmaj_world > bmin_world) { + // carta rotation is from y-axis, ellipse rotation is from x-axis + ellipse_rotation += 90.0; + } else { + // swapping takes care of 90 deg adjustment + std::swap(bmaj_world, bmin_world); + } + + wcs_points.push_back(bmaj_world); + wcs_points.push_back(bmin_world); + return true; +} + +bool RegionConverter::CartaPointToWorld(const CARTA::Point& point, std::vector& world_point) { + // Converts a CARTA point(x, y) in pixel coordinates to a Quantity vector [x, y] in world coordinates. + // Returns whether conversion was successful + + // Vectors must be same number of axes as in coord system for conversion: + int naxes(_reference_coord_sys->nPixelAxes()); + casacore::Vector pixel_values(naxes), world_values(naxes); + pixel_values = 0.0; // set "extra" axes to 0, not needed + pixel_values(0) = point.x(); + pixel_values(1) = point.y(); + + // convert pixel vector to world vector + if (!_reference_coord_sys->toWorld(world_values, pixel_values)) { + return false; + } + + // Set Quantities from world values and units + casacore::Vector world_units = _reference_coord_sys->worldAxisUnits(); + world_point.clear(); + world_point.push_back(casacore::Quantity(world_values(0), world_units(0))); + world_point.push_back(casacore::Quantity(world_values(1), world_units(1))); + return true; +} + +// ************************************************************************* +// Convert region to any image + +std::shared_ptr RegionConverter::GetCachedLCRegion(int file_id, bool use_approx_polygon) { + // Return cached region applied to image with file_id, if cached. + std::lock_guard guard(_region_mutex); + if (_converted_regions.find(file_id) != _converted_regions.end()) { + return _converted_regions.at(file_id); + } else if (use_approx_polygon && _polygon_regions.find(file_id) != _polygon_regions.end()) { + // Return cached polygon-approximated region applied to image with file_id + return _polygon_regions.at(file_id); + } + + return std::shared_ptr(); +} + +std::shared_ptr RegionConverter::GetImageRegion(int file_id, std::shared_ptr output_csys, + const casacore::IPosition& output_shape, const StokesSource& stokes_source, bool report_error) { + // Apply region to non-reference image, possibly as an approximate polygon to avoid distortion. + std::shared_ptr lc_region; + + // Analytic, closed regions only + if (_region_state.IsLineType() || _region_state.IsAnnotation()) { + return lc_region; + } + + if (stokes_source.IsOriginalImage()) { + // The cache of converted LCRegions is only for the original image (not computed stokes image). In order to avoid the ambiguity + lc_region = GetCachedLCRegion(file_id); + } + + if (!lc_region) { + // Create converted LCRegion and cache it + bool cache_polygon(false); + bool use_polygon = UseApproximatePolygon(output_csys); + if (_region_state.IsRotbox()) { + // Rotbox is always converted from a polygon, either box corners only (use_polygon==false) or many points (use_polygon==true) + lc_region = GetAppliedPolygonRegion(file_id, output_csys, output_shape, use_polygon); + cache_polygon = true; + } else if (use_polygon) { // distortion + // Approximate region as polygon points then convert the points. + spdlog::debug("Using polygon approximation to avoid distortion in matched region"); + lc_region = GetAppliedPolygonRegion(file_id, output_csys, output_shape, use_polygon); + cache_polygon = true; + } else { // no distortion + // Do direct region conversion from reference WCRegion (no distortion detected) + lc_region = GetConvertedLCRegion(file_id, output_csys, output_shape, stokes_source, report_error); + if (lc_region) { + // Region conversion succeeded + spdlog::debug("Using direct region conversion for matched image"); + } else { + // Region conversion failed (e.g. outside image) but can convert points to world + spdlog::debug("Using polygon approximation for failed conversion of matched region"); + lc_region = GetAppliedPolygonRegion(file_id, output_csys, output_shape, use_polygon); + cache_polygon = true; + } + } + + if (cache_polygon && lc_region && stokes_source.IsOriginalImage()) { + // Cache converted polygon, only for the original image (not computed stokes). + std::lock_guard guard(_region_mutex); + _polygon_regions[file_id] = lc_region; + } + } + + return lc_region; +} + +std::shared_ptr RegionConverter::GetConvertedLCRegion(int file_id, + std::shared_ptr output_csys, const casacore::IPosition& output_shape, const StokesSource& stokes_source, + bool report_error) { + // Convert reference WCRegion to LCRegion in output coord_sys and shape, and cache converted region. + // Check cache before calling this else will needlessly create a new LCRegion and cache it. + std::shared_ptr lc_region; + + try { + // Convert reference WCRegion to LCRegion using output csys and shape + if (!_reference_region_set) { + SetReferenceWCRegion(); + } + + if (ReferenceRegionValid()) { + std::shared_ptr reference_region = std::atomic_load(&_reference_region); + lc_region.reset(reference_region->toLCRegion(*output_csys.get(), output_shape)); + } + } catch (const casacore::AipsError& err) { + if (report_error) { + spdlog::error("Error converting region type {} to file {}: {}", _region_state.type, file_id, err.getMesg()); + } + } + + if (lc_region && stokes_source.IsOriginalImage()) { + // Cache the lattice coordinate region only for the original image (not computed stokes image). + std::lock_guard guard(_region_mutex); + _converted_regions[file_id] = lc_region; + } + + return lc_region; +} + +// ************************************************************************* +// Region as polygon to avoid distortion in matched image + +bool RegionConverter::UseApproximatePolygon(std::shared_ptr output_csys) { + // Determine whether to convert region directly, or approximate it as a polygon in the output image. + // Closed region types: rectangle, ellipse, polygon. + // Check ellipse and rectangle distortion; always use polygon for polygon regions. + CARTA::RegionType region_type = _region_state.type; + if ((region_type != CARTA::RegionType::ELLIPSE) && (region_type != CARTA::RegionType::RECTANGLE)) { + return true; + } + + // Ratio of vector lengths in reference image region + double ref_length_ratio; + double x_length(_region_state.control_points[1].x()), y_length(_region_state.control_points[1].y()); + if (region_type == CARTA::RegionType::ELLIPSE) { + ref_length_ratio = x_length / y_length; + } else { + ref_length_ratio = y_length / x_length; + } + + // Make vector of endpoints and center, to check lengths against reference image lengths + std::vector points; // [p0, p1, p2, p3, center] + if (region_type == CARTA::RegionType::ELLIPSE) { + // Make "polygon" with only 4 points + points = GetApproximateEllipsePoints(4); + } else { + // Get midpoints of 4 sides of rectangle + points = GetRectangleMidpoints(); + } + points.push_back(_region_state.control_points[0]); + + // Convert reference pixel points to output pixel points, then check vector length ratio and dot product + casacore::Vector x, y; + if (PointsToImagePixels(points, output_csys, x, y)) { + // vector0 is (center, p0), vector1 is (center, p1) + auto v0_delta_x = x[0] - x[4]; + auto v0_delta_y = y[0] - y[4]; + auto v1_delta_x = x[1] - x[4]; + auto v1_delta_y = y[1] - y[4]; + + // Compare reference length ratio to converted length ratio + // Ratio of vector lengths in converted region + auto v0_length = sqrt((v0_delta_x * v0_delta_x) + (v0_delta_y * v0_delta_y)); + auto v1_length = sqrt((v1_delta_x * v1_delta_x) + (v1_delta_y * v1_delta_y)); + double converted_length_ratio = v1_length / v0_length; + double length_ratio_difference = fabs(ref_length_ratio - converted_length_ratio); + // spdlog::debug("{} distortion check: length ratio difference={:.3e}", length_ratio_difference); + if (length_ratio_difference > 1e-4) { + // Failed ratio check, use polygon + return true; + } + + // Passed ratio check; check dot product of converted region + double converted_dot_product = (v0_delta_x * v1_delta_x) + (v0_delta_y * v1_delta_y); + // spdlog::debug("{} distortion check: dot product={:.3e}", converted_dot_product); + if (fabs(converted_dot_product) > 1e-2) { + // failed dot product test, use polygon + return true; + } + } else { + spdlog::error("Error converting region points to matched image."); + return true; + } + + return false; +} + +std::vector RegionConverter::GetRectangleMidpoints() { + // Return midpoints of 4 sides of rectangle + // Find corners with rotation: blc, brc, trc, tlc + std::vector midpoints; + casacore::Vector x, y; + if (_region_state.GetRectangleCorners(x, y)) { + // start with right side brc, trc + midpoints.push_back(Message::Point((x[1] + x[2]) / 2.0, (y[1] + y[2]) / 2.0)); + midpoints.push_back(Message::Point((x[2] + x[3]) / 2.0, (y[2] + y[3]) / 2.0)); + midpoints.push_back(Message::Point((x[3] + x[0]) / 2.0, (y[3] + y[0]) / 2.0)); + midpoints.push_back(Message::Point((x[0] + x[1]) / 2.0, (y[0] + y[1]) / 2.0)); + } + return midpoints; +} + +std::shared_ptr RegionConverter::GetAppliedPolygonRegion( + int file_id, std::shared_ptr output_csys, const casacore::IPosition& output_shape, bool has_distortion) { + // Approximate region as polygon pixel vertices, and convert to given csys. + // If has distortion, use DEFAULT_VERTEX_COUNT vertices (else is rotbox and just use corners) + std::shared_ptr lc_region; + + bool is_point(_region_state.IsPoint()); + size_t nvertices(is_point ? 1 : DEFAULT_VERTEX_COUNT); + + // Set reference region as points along polygon segments + auto polygon_points = GetReferencePolygonPoints(nvertices, has_distortion); + if (polygon_points.empty()) { + return lc_region; + } + + // Convert polygon points to x and y pixel coords in matched image + casacore::Vector x, y; + if (polygon_points.size() == 1) { + // Point, ellipse, and rotbox with no distortion have one vector for all points + if (!PointsToImagePixels(polygon_points[0], output_csys, x, y)) { + spdlog::error("Error approximating region as polygon in matched image."); + return lc_region; + } + if (!is_point && has_distortion) { + // if ~horizontal then remove intermediate points to fix kinks + RemoveHorizontalPolygonPoints(x, y); + } + } else { + // Rectangle and polygon have one vector for each segment of original rectangle/polygon + for (auto& segment : polygon_points) { + casacore::Vector segment_x, segment_y; + if (!PointsToImagePixels(segment, output_csys, segment_x, segment_y)) { + spdlog::error("Error approximating region as polygon in matched image."); + return lc_region; + } + + if (has_distortion) { + // if ~horizontal then remove intermediate points to fix "kinks" in mask + RemoveHorizontalPolygonPoints(segment_x, segment_y); + } + + // Resize x and y to append selected segment points + auto old_size = x.size(); + x.resize(old_size + segment_x.size(), true); + y.resize(old_size + segment_y.size(), true); + for (auto i = 0; i < segment_x.size(); ++i) { + x[old_size + i] = segment_x[i]; + y[old_size + i] = segment_y[i]; + } + } + } + + // Use converted pixel points to create LCRegion (LCBox for point, else LCPolygon) + try { + if (is_point) { + // Point is not a polygon (needs at least 3 points), use LCBox instead with blc, trc = point + size_t ndim(output_shape.size()); + casacore::Vector blc(ndim, 0.0), trc(ndim); + blc(0) = x(0); + blc(1) = y(0); + trc(0) = x(0); + trc(1) = y(0); + + for (size_t i = 2; i < ndim; ++i) { + trc(i) = output_shape(i) - 1; + } + + lc_region.reset(new casacore::LCBox(blc, trc, output_shape)); + } else { + // Need 2D shape + casacore::IPosition keep_axes(2, 0, 1); + casacore::IPosition region_shape(output_shape.keepAxes(keep_axes)); + lc_region.reset(new casacore::LCPolygon(x, y, region_shape)); + } + } catch (const casacore::AipsError& err) { + spdlog::error("Cannot apply region type {} to file {}: {}", _region_state.type, file_id, err.getMesg()); + } + + return lc_region; +} + +std::vector> RegionConverter::GetReferencePolygonPoints(int num_vertices, bool has_distortion) { + // Approximate reference region as polygon with input number of vertices. + // Returns points for supported region types. + std::vector> points; + + switch (_region_state.type) { + case CARTA::POINT: { + points.push_back(_region_state.control_points); + } + case CARTA::RECTANGLE: + case CARTA::POLYGON: { + return GetApproximatePolygonPoints(num_vertices, has_distortion); + } + case CARTA::ELLIPSE: { + points.push_back(GetApproximateEllipsePoints(num_vertices)); + } + default: + return points; + } + return points; +} + +std::vector> RegionConverter::GetApproximatePolygonPoints(int num_vertices, bool has_distortion) { + // Approximate RECTANGLE or POLYGON region as polygon with num_vertices if has distortion. Else convert rotbox to polygon. + // Returns vector of points for each segment of polygon/rectangle, or empty vector for other region types. + std::vector> polygon_points; + std::vector region_vertices; + CARTA::RegionType region_type(_region_state.type); + + // Rectangle corners or polygon points as polygon vertices for segments. + if (region_type == CARTA::RegionType::RECTANGLE) { + casacore::Vector x, y; + _region_state.GetRectangleCorners(x, y); + + for (size_t i = 0; i < x.size(); ++i) { + region_vertices.push_back(Message::Point(x(i), y(i))); + } + } else if (region_type == CARTA::RegionType::POLYGON) { + region_vertices = _region_state.control_points; + } else { + spdlog::error("Error approximating region as polygon: type {} not supported", _region_state.type); + return polygon_points; + } + + // Close polygon + CARTA::Point first_point(region_vertices[0]); + region_vertices.push_back(first_point); + + if (has_distortion) { + // Divide each polygon/rectangle segment into target number of segments with target length + double total_length = GetTotalSegmentLength(region_vertices); + double min_length = 0.1; + double target_segment_length = total_length / num_vertices; + target_segment_length = std::max(target_segment_length, min_length); + + for (size_t i = 1; i < region_vertices.size(); ++i) { + // Handle segment from point[i-1] to point[i] + std::vector segment_points; + + auto delta_x = region_vertices[i].x() - region_vertices[i - 1].x(); + auto delta_y = region_vertices[i].y() - region_vertices[i - 1].y(); + auto segment_length = sqrt((delta_x * delta_x) + (delta_y * delta_y)); + auto dir_x = delta_x / segment_length; + auto dir_y = delta_y / segment_length; + auto target_nsegment = round(segment_length / target_segment_length); + auto target_length = segment_length / target_nsegment; + + auto first_segment_point(region_vertices[i - 1]); + segment_points.push_back(first_segment_point); + + auto first_x(first_segment_point.x()); + auto first_y(first_segment_point.y()); + + for (size_t j = 1; j < target_nsegment; ++j) { + auto length_from_first = j * target_length; + auto x_offset = dir_x * length_from_first; + auto y_offset = dir_y * length_from_first; + segment_points.push_back(Message::Point(first_x + x_offset, first_y + y_offset)); + } + + polygon_points.push_back(segment_points); + } + } else { + // No need to add extra points, no distortion + polygon_points.push_back(region_vertices); + } + + return polygon_points; +} + +std::vector RegionConverter::GetApproximateEllipsePoints(int num_vertices) { + // Approximate ELLIPSE region as polygon with num_vertices, return points + std::vector polygon_points; + + auto cx = _region_state.control_points[0].x(); + auto cy = _region_state.control_points[0].y(); + auto bmaj = _region_state.control_points[1].x(); + auto bmin = _region_state.control_points[1].y(); + + auto delta_theta = 2.0 * M_PI / num_vertices; + auto rotation = _region_state.rotation * M_PI / 180.0; + auto cos_rotation = cos(rotation); + auto sin_rotation = sin(rotation); + + for (int i = 0; i < num_vertices; ++i) { + auto theta = i * delta_theta; + auto rot_bmin = bmin * cos(theta); + auto rot_bmaj = bmaj * sin(theta); + + auto x_offset = (cos_rotation * rot_bmin) - (sin_rotation * rot_bmaj); + auto y_offset = (sin_rotation * rot_bmin) + (cos_rotation * rot_bmaj); + + polygon_points.push_back(Message::Point(cx + x_offset, cy + y_offset)); + } + + return polygon_points; +} + +double RegionConverter::GetTotalSegmentLength(std::vector& points) { + // Accumulate length of each point-to-point segment; returns total length. + double total_length(0.0); + + for (size_t i = 1; i < points.size(); ++i) { + auto delta_x = points[i].x() - points[i - 1].x(); + auto delta_y = points[i].y() - points[i - 1].y(); + total_length += sqrt((delta_x * delta_x) + (delta_y * delta_y)); + } + + return total_length; +} + +void RegionConverter::RemoveHorizontalPolygonPoints(casacore::Vector& x, casacore::Vector& y) { + // When polygon points have close y-points (horizontal segment), the x-range is masked only to the next point. + // Remove points not near integral pixel. + std::vector keep_x, keep_y; + size_t npoints(x.size()); + + for (int i = 0; i < npoints - 2; ++i) { + if (i == 0) { + // always include first point of segment + keep_x.push_back(x[i]); + keep_y.push_back(y[i]); + continue; + } + + float this_y = y[i]; + float next_y = y[i + 1]; + if (!ValuesNear(this_y, next_y)) { + // Line connecting points not ~horizontal - keep point + keep_x.push_back(x[i]); + keep_y.push_back(y[i]); + } + } + + if (keep_x.size() < npoints) { + // Set to new vector with points removed + x = casacore::Vector(keep_x); + y = casacore::Vector(keep_y); + } +} + +bool RegionConverter::ValuesNear(float val1, float val2) { + // near and nearAbs in casacore Math + return val1 == 0 || val2 == 0 ? casacore::nearAbs(val1, val2) : casacore::near(val1, val2); +} + +// *************************************************************** +// Apply region to any image and return LCRegion Record for export + +casacore::TableRecord RegionConverter::GetImageRegionRecord( + int file_id, std::shared_ptr output_csys, const casacore::IPosition& output_shape) { + // Return Record describing Region applied to output image in pixel coordinates + casacore::TableRecord record; + + // No LCRegion for lines. LCRegion for rotbox is a polygon; Record should be for unrotated box with rotation field. + if (!_region_state.IsLineType() && !_region_state.IsRotbox()) { + // Get record from converted LCRegion, for enclosed regions only. + // Check converted regions cache (but not polygon regions cache) + std::shared_ptr lc_region = GetCachedLCRegion(file_id, false); + + // Convert reference region to output image + if (!lc_region) { + lc_region = GetConvertedLCRegion(file_id, output_csys, output_shape); + } + + // Get LCRegion definition as Record + if (lc_region) { + record = lc_region->toRecord("region"); + if (record.isDefined("region")) { + record = record.asRecord("region"); + } + } + } + + if (record.empty()) { + // LCRegion failed: is outside the image or is a line or rotated rectangle. + // Manually convert control points and put in Record. + record = GetRegionPointsRecord(output_csys, output_shape); + } + return record; +} + +casacore::TableRecord RegionConverter::GetRegionPointsRecord( + std::shared_ptr output_csys, const casacore::IPosition& output_shape) { + // Convert control points to output coord sys if needed, and return completed record. + // Used when LCRegion::toRecord() fails, usually outside image. + casacore::TableRecord record; + if (!_reference_region_set) { + SetReferenceWCRegion(); // set wcs control points + } + + switch (_region_state.type) { + case CARTA::RegionType::POINT: + case CARTA::RegionType::ANNPOINT: { + record = GetPointRecord(output_csys, output_shape); + break; + } + case CARTA::RegionType::LINE: + case CARTA::RegionType::POLYLINE: + case CARTA::RegionType::ANNLINE: + case CARTA::RegionType::ANNPOLYLINE: + case CARTA::RegionType::ANNVECTOR: + case CARTA::RegionType::ANNRULER: { + record = GetLineRecord(output_csys); + break; + } + case CARTA::RegionType::RECTANGLE: + case CARTA::RegionType::POLYGON: + case CARTA::RegionType::ANNRECTANGLE: + case CARTA::RegionType::ANNPOLYGON: + case CARTA::RegionType::ANNTEXT: { + // Rectangle types are LCPolygon with 4 (unrotated) corners. + record = (_region_state.IsRotbox() ? GetRotboxRecord(output_csys) : GetPolygonRecord(output_csys)); + break; + } + case CARTA::RegionType::ELLIPSE: + case CARTA::RegionType::ANNELLIPSE: + case CARTA::RegionType::ANNCOMPASS: { + record = GetEllipseRecord(output_csys); + break; + } + default: + break; + } + return record; +} + +casacore::TableRecord RegionConverter::GetPointRecord( + std::shared_ptr output_csys, const casacore::IPosition& output_shape) { + // Convert wcs points to output image in format of LCBox::toRecord() + casacore::TableRecord record; + try { + // wcs control points is single point (x, y) + casacore::Vector pixel_point; + if (WorldPointToImagePixels(_wcs_control_points, output_csys, pixel_point)) { + auto ndim = output_shape.size(); + casacore::Vector blc(ndim), trc(ndim); + blc = 0.0; + blc(0) = pixel_point(0); + blc(1) = pixel_point(1); + trc(0) = pixel_point(0); + trc(1) = pixel_point(1); + + for (size_t i = 2; i < ndim; ++i) { + trc(i) = output_shape(i) - 1; + } + + record.define("name", "LCBox"); + record.define("blc", blc); + record.define("trc", trc); + } else { + spdlog::error("Error converting point to image."); + } + } catch (const casacore::AipsError& err) { + spdlog::error("Error converting point to image: {}", err.getMesg()); + } + + return record; +} + +casacore::TableRecord RegionConverter::GetLineRecord(std::shared_ptr image_csys) { + // Convert control points for line-type region to output image pixels in format of LCPolygon::toRecord() + casacore::TableRecord record; + casacore::Vector x, y; + if (PointsToImagePixels(_region_state.control_points, image_csys, x, y)) { + record.define("name", _region_state.GetLineRegionName()); + record.define("x", x); + record.define("y", y); + } + return record; +} + +casacore::TableRecord RegionConverter::GetPolygonRecord(std::shared_ptr output_csys) { + // Convert wcs points to output image in format of LCPolygon::toRecord() + // This is for POLYGON or RECTANGLE (points are four corners of box) + casacore::TableRecord record; + auto type = _region_state.type; + + try { + size_t npoints(_wcs_control_points.size() / 2); + casacore::Vector x(npoints), y(npoints); // Record fields + + // Convert each wcs control point to pixel coords in output csys + for (size_t i = 0; i < _wcs_control_points.size(); i += 2) { + std::vector world_point(2); + world_point[0] = _wcs_control_points[i]; + world_point[1] = _wcs_control_points[i + 1]; + casacore::Vector pixel_point; + if (WorldPointToImagePixels(world_point, output_csys, pixel_point)) { + // Add to x and y Vectors + int index(i / 2); + x(index) = pixel_point(0); + y(index) = pixel_point(1); + } else { + spdlog::error("Error converting region type {} to image pixels.", type); + return record; + } + } + + if (type == CARTA::RegionType::POLYGON) { + // LCPolygon::toRecord adds first point as last point to close region + x.resize(npoints + 1, true); + x(npoints) = x(0); + y.resize(npoints + 1, true); + y(npoints) = y(0); + } + + // Add fields for this region type + record.define("name", "LCPolygon"); + record.define("x", x); + record.define("y", y); + } catch (const casacore::AipsError& err) { + spdlog::error("Error converting region type () to image: {}", type, err.getMesg()); + } + + return record; +} + +casacore::TableRecord RegionConverter::GetRotboxRecord(std::shared_ptr output_csys) { + // Convert control points for rotated box-type region (ignore rotation) to output image pixels + // in format of LCPolygon::toRecord(). + casacore::TableRecord record; + try { + // Get 4 corner points (unrotated) in pixel coordinates + casacore::Vector x, y; + bool apply_rotation(false); + if (_region_state.GetRectangleCorners(x, y, apply_rotation)) { + // Convert corners to reference world coords. + // Cannot use wcs control points because rotation was applied. + int num_axes(_reference_coord_sys->nPixelAxes()); + auto num_points(x.size()); + casacore::Matrix pixel_coords(num_axes, num_points); + casacore::Matrix world_coords(num_axes, num_points); + pixel_coords = 0.0; + pixel_coords.row(0) = x; + pixel_coords.row(1) = y; + casacore::Vector failures; + if (!_reference_coord_sys->toWorldMany(world_coords, pixel_coords, failures)) { + spdlog::error("Error converting rectangle pixel coordinates to world."); + return record; + } + + // Convert reference world coord points to output pixel points + casacore::Vector ref_x_world = world_coords.row(0); + casacore::Vector ref_y_world = world_coords.row(1); + casacore::Vector ref_world_units = _reference_coord_sys->worldAxisUnits(); + casacore::Vector out_x_pix(num_points), out_y_pix(num_points); + for (size_t i = 0; i < num_points; i++) { + // Reference world point as Quantity + std::vector ref_world_point; + ref_world_point.push_back(casacore::Quantity(ref_x_world(i), ref_world_units(0))); + ref_world_point.push_back(casacore::Quantity(ref_y_world(i), ref_world_units(1))); + // Convert to output pixel point + casacore::Vector out_pixel_point; + if (WorldPointToImagePixels(ref_world_point, output_csys, out_pixel_point)) { + out_x_pix(i) = out_pixel_point(0); + out_y_pix(i) = out_pixel_point(1); + } else { + spdlog::error("Error converting rectangle coordinates to image."); + return record; + } + } + + // Add fields for this region type + record.define("name", "LCPolygon"); + record.define("x", out_x_pix); + record.define("y", out_y_pix); + } + } catch (const casacore::AipsError& err) { + spdlog::error("Error converting rectangle to image: {}", err.getMesg()); + } + return record; +} + +casacore::TableRecord RegionConverter::GetEllipseRecord(std::shared_ptr output_csys) { + // Convert wcs points to output image in format of LCEllipsoid::toRecord() + casacore::TableRecord record; + casacore::Vector center(2), radii(2); // for record + + // Center point + std::vector ref_world_point(2); + ref_world_point[0] = _wcs_control_points[0]; + ref_world_point[1] = _wcs_control_points[1]; + casacore::Vector out_pixel_point; + + try { + if (WorldPointToImagePixels(ref_world_point, output_csys, out_pixel_point)) { + center(0) = out_pixel_point(0); + center(1) = out_pixel_point(1); + + // Convert radii to output world units, then to pixels using increment + casacore::Quantity bmaj = _wcs_control_points[2]; + casacore::Quantity bmin = _wcs_control_points[3]; + casacore::Vector out_increments(output_csys->increment()); + casacore::Vector out_units(output_csys->worldAxisUnits()); + bmaj.convert(out_units(0)); + bmin.convert(out_units(1)); + radii(0) = fabs(bmaj.getValue() / out_increments(0)); + radii(1) = fabs(bmin.getValue() / out_increments(1)); + + // Add fields for this region type + record.define("name", "LCEllipsoid"); + record.define("center", center); + record.define("radii", radii); + + // LCEllipsoid measured from major (x) axis + casacore::Quantity theta = casacore::Quantity(_region_state.rotation + 90.0, "deg"); + theta.convert("rad"); + record.define("theta", theta.getValue()); + } else { + spdlog::error("Incompatible coordinate systems for ellipse conversion."); + } + } catch (const casacore::AipsError& err) { + spdlog::error("Error converting ellipse to image: {}", err.getMesg()); + } + + return record; +} + +// ************************************************************************* +// Utilities for pixel/world conversion + +bool RegionConverter::PointsToImagePixels(const std::vector& points, std::shared_ptr output_csys, + casacore::Vector& x, casacore::Vector& y) { + // Convert pixel coords in reference image (points) to pixel coords in output image coordinate system (x and y). + // ref pixels -> ref world -> output world -> output pixels + bool converted(true); + try { + // Convert each pixel point to output csys pixel + size_t npoints(points.size()); + x.resize(npoints); + y.resize(npoints); + + for (auto i = 0; i < npoints; ++i) { + // Convert pixel to world (reference image) [x, y] + std::vector world_point; + if (CartaPointToWorld(points[i], world_point)) { + // Convert world to pixel (output image) [x, y] + casacore::Vector pixel_point; + if (WorldPointToImagePixels(world_point, output_csys, pixel_point)) { + x(i) = pixel_point(0); + y(i) = pixel_point(1); + } else { // world to pixel failed + spdlog::error("Error converting region to output image pixel coords."); + converted = false; + break; + } + } else { // pixel to world failed + spdlog::error("Error converting region to reference image world coords."); + converted = false; + break; + } + } + } catch (const casacore::AipsError& err) { + spdlog::error("Error converting region to output image: {}", err.getMesg()); + converted = false; + } + + return converted; +} + +bool RegionConverter::WorldPointToImagePixels(std::vector& world_point, + std::shared_ptr output_csys, casacore::Vector& pixel_point) { + // Convert reference world-coord point to output pixel-coord point: ref world -> output world -> output pixels. + // Both images must have direction coordinates or linear coordinates. + // Returns pixel points with success or throws exception (catch in calling function). + bool success(false); + if (_reference_coord_sys->hasDirectionCoordinate() && output_csys->hasDirectionCoordinate()) { + // Input and output direction reference frames + casacore::MDirection::Types reference_dir_type = _reference_coord_sys->directionCoordinate().directionType(); + casacore::MDirection::Types output_dir_type = output_csys->directionCoordinate().directionType(); + + // Convert world point from reference to output coord sys + casacore::MDirection world_direction(world_point[0], world_point[1], reference_dir_type); + if (reference_dir_type != output_dir_type) { + world_direction = casacore::MDirection::Convert(world_direction, output_dir_type)(); + } + + // Convert output world point to pixel point + output_csys->directionCoordinate().toPixel(pixel_point, world_direction); + success = true; + } else if (_reference_coord_sys->hasLinearCoordinate() && output_csys->hasLinearCoordinate()) { + // Get linear axes indices + auto indices = output_csys->linearAxesNumbers(); + if (indices.size() != 2) { + return false; + } + // Input and output linear frames + casacore::Vector output_units = output_csys->worldAxisUnits(); + casacore::Vector world_point_value(output_csys->nWorldAxes(), 0); + world_point_value(indices(0)) = world_point[0].get(output_units(indices(0))).getValue(); + world_point_value(indices(1)) = world_point[1].get(output_units(indices(1))).getValue(); + + // Convert world point to output pixel point + casacore::Vector tmp_pixel_point; + output_csys->toPixel(tmp_pixel_point, world_point_value); + + // Only fill the pixel coordinate results + pixel_point.resize(2); + pixel_point(0) = tmp_pixel_point(indices(0)); + pixel_point(1) = tmp_pixel_point(indices(1)); + success = true; + } + return success; +} diff --git a/src/Region/RegionConverter.h b/src/Region/RegionConverter.h new file mode 100644 index 000000000..1d13d89ce --- /dev/null +++ b/src/Region/RegionConverter.h @@ -0,0 +1,103 @@ +/* This file is part of the CARTA Image Viewer: https://github.com/CARTAvis/carta-backend + Copyright 2018- Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), + Associated Universities, Inc. (AUI) and the Inter-University Institute for Data Intensive Astronomy (IDIA) + SPDX-License-Identifier: GPL-3.0-or-later +*/ + +//# RegionConverter.h: class for managing region conversion to matched image + +#ifndef CARTA_SRC_REGION_REGIONCONVERTER_H_ +#define CARTA_SRC_REGION_REGIONCONVERTER_H_ + +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include "RegionState.h" +#include "Util/Stokes.h" + +#define DEFAULT_VERTEX_COUNT 1000 + +namespace carta { + +class RegionConverter { +public: + RegionConverter(const RegionState& state, std::shared_ptr csys); + + // LCRegion and mask for region applied to image. Must be a closed region (not line) and not annotation. + std::shared_ptr GetCachedLCRegion(int file_id, bool use_approx_polygon = true); + std::shared_ptr GetImageRegion(int file_id, std::shared_ptr csys, + const casacore::IPosition& shape, const StokesSource& stokes_source = StokesSource(), bool report_error = true); + casacore::TableRecord GetImageRegionRecord( + int file_id, std::shared_ptr csys, const casacore::IPosition& shape); + +private: + // Reference region in world coordinates (reference coordinate system) + void SetReferenceWCRegion(); + bool ReferenceRegionValid(); + bool RectangleControlPointsToWorld(std::vector& wcs_points); + bool EllipseControlPointsToWorld(std::vector& wcs_points, float& ellipse_rotation); + bool CartaPointToWorld(const CARTA::Point& point, std::vector& world_point); + + // Region converted directly to image (not polygon approximation) + std::shared_ptr GetConvertedLCRegion(int file_id, std::shared_ptr output_csys, + const casacore::IPosition& output_shape, const StokesSource& stokes_source = StokesSource(), bool report_error = true); + + // Reference region converted to image as approximate polygon, cached with file_id + bool UseApproximatePolygon(std::shared_ptr output_csys); + std::vector GetRectangleMidpoints(); + std::shared_ptr GetAppliedPolygonRegion( + int file_id, std::shared_ptr output_csys, const casacore::IPosition& output_shape, bool has_distortion); + std::vector> GetReferencePolygonPoints(int num_vertices, bool has_distortion); + std::vector> GetApproximatePolygonPoints(int num_vertices, bool has_distortion); + std::vector GetApproximateEllipsePoints(int num_vertices); + double GetTotalSegmentLength(std::vector& points); + void RemoveHorizontalPolygonPoints(casacore::Vector& x, casacore::Vector& y); + bool ValuesNear(float val1, float val2); + + // Record for region converted to pixel coords in output image, returned in LCRegion Record format + casacore::TableRecord GetRegionPointsRecord( + std::shared_ptr output_csys, const casacore::IPosition& output_shape); + void CompleteLCRegionRecord(casacore::TableRecord& record, const casacore::IPosition& shape); + casacore::TableRecord GetPointRecord(std::shared_ptr output_csys, const casacore::IPosition& output_shape); + casacore::TableRecord GetLineRecord(std::shared_ptr image_csys); + casacore::TableRecord GetPolygonRecord(std::shared_ptr output_csys); + casacore::TableRecord GetRotboxRecord(std::shared_ptr output_csys); + casacore::TableRecord GetEllipseRecord(std::shared_ptr output_csys); + + // Utilities for pixel/world conversion + bool PointsToImagePixels(const std::vector& points, std::shared_ptr output_csys, + casacore::Vector& x, casacore::Vector& y); + // World point as (x,y) quantities to output pixel point as (x,y) vector + bool WorldPointToImagePixels(std::vector& world_point, std::shared_ptr output_csys, + casacore::Vector& pixel_point); + + // Reference image region parameters + RegionState _region_state; + std::shared_ptr _reference_coord_sys; + + // Reference region: control points converted to world coords, + // cached as WCRegion, and flag that WCRegion was attempted + // (may be null if outside coordinate system) + std::vector _wcs_control_points; + std::shared_ptr _reference_region; + bool _reference_region_set; + + // Converted regions: reference region applied to image, + // or as polygon approximation applied to image. + // Map key is file_id + std::mutex _region_mutex; + std::unordered_map> _converted_regions; + std::unordered_map> _polygon_regions; +}; + +} // namespace carta + +#endif // CARTA_SRC_REGION_REGIONCONVERTER_H_ diff --git a/src/Region/RegionHandler.cc b/src/Region/RegionHandler.cc index bf50dc4a7..609daf957 100644 --- a/src/Region/RegionHandler.cc +++ b/src/Region/RegionHandler.cc @@ -96,14 +96,6 @@ bool RegionHandler::SetRegion(int& region_id, RegionState& region_state, std::sh return valid_region; } -bool RegionHandler::RegionChanged(int region_id) { - // Used to trigger sending profiles etc., so not for annotation regions - if (!RegionSet(region_id, true)) { - return false; - } - return GetRegion(region_id)->RegionChanged(); -} - void RegionHandler::RemoveRegion(int region_id) { // Call destructor and erase from map if (!RegionSet(region_id)) { @@ -258,10 +250,9 @@ void RegionHandler::ExportRegion(int file_id, std::shared_ptr frame, CART } } - bool pixel_coord(coord_type == CARTA::CoordinateType::PIXEL); + bool export_pixel_coord(coord_type == CARTA::CoordinateType::PIXEL); auto output_csys = frame->CoordinateSystem(); - - if (!pixel_coord && !output_csys->hasDirectionCoordinate()) { + if (!export_pixel_coord && !output_csys->hasDirectionCoordinate()) { // Export fails, cannot convert to world coordinates export_ack.set_success(false); export_ack.set_message("Cannot export regions in world coordinates for linear coordinate system."); @@ -275,7 +266,7 @@ void RegionHandler::ExportRegion(int file_id, std::shared_ptr frame, CART exporter = std::unique_ptr(new CrtfImportExport(output_csys, output_shape, frame->StokesAxis())); break; case CARTA::FileType::DS9_REG: - exporter = std::unique_ptr(new Ds9ImportExport(output_csys, output_shape, pixel_coord)); + exporter = std::unique_ptr(new Ds9ImportExport(output_csys, output_shape, export_pixel_coord)); break; default: break; @@ -291,7 +282,7 @@ void RegionHandler::ExportRegion(int file_id, std::shared_ptr frame, CART auto region = GetRegion(region_id); auto region_state = region->GetRegionState(); - if ((region_state.reference_file_id == file_id) && pixel_coord) { + if ((region_state.reference_file_id == file_id) && export_pixel_coord) { // Use RegionState control points with reference file id for pixel export region_added = exporter->AddExportRegion(region_state, region_style); } else { @@ -299,7 +290,7 @@ void RegionHandler::ExportRegion(int file_id, std::shared_ptr frame, CART // Use Record containing pixel coords of region converted to output image casacore::TableRecord region_record = region->GetImageRegionRecord(file_id, output_csys, output_shape); if (!region_record.empty()) { - region_added = exporter->AddExportRegion(region_state, region_style, region_record, pixel_coord); + region_added = exporter->AddExportRegion(region_state, region_style, region_record, export_pixel_coord); } } catch (const casacore::AipsError& err) { spdlog::error("Export region record failed: {}", err.getMesg()); @@ -2284,6 +2275,7 @@ bool RegionHandler::FillPointSpatialProfileData(int file_id, int region_id, std: bool RegionHandler::FillLineSpatialProfileData(int file_id, int region_id, std::function cb) { // Line spatial profiles. Use callback to return each profile individually. + Timer t; if (!RegionFileIdsValid(region_id, file_id, true)) { return false; } @@ -2342,6 +2334,7 @@ bool RegionHandler::FillLineSpatialProfileData(int file_id, int region_id, std:: }); } + spdlog::performance("Line spatial data in {:.3f} ms", t.Elapsed().ms()); return profile_ok; } @@ -2380,13 +2373,15 @@ bool RegionHandler::GetLineSpatialData(int file_id, int region_id, const std::st } bool RegionHandler::IsPointRegion(int region_id) { + // Analytic region, not annotation if (RegionSet(region_id, true)) { - return GetRegion(region_id)->IsPoint(); + return GetRegion(region_id)->IsPoint() && !GetRegion(region_id)->IsAnnotation(); } return false; } bool RegionHandler::IsLineRegion(int region_id) { + // Analytic region, not annotation if (RegionSet(region_id, true)) { return GetRegion(region_id)->IsLineType() && !GetRegion(region_id)->IsAnnotation(); } @@ -2394,6 +2389,7 @@ bool RegionHandler::IsLineRegion(int region_id) { } bool RegionHandler::IsClosedRegion(int region_id) { + // Analytic region, not annotation if (RegionSet(region_id, true)) { auto type = GetRegion(region_id)->GetRegionState().type; return (type == CARTA::RegionType::RECTANGLE) || (type == CARTA::RegionType::ELLIPSE) || (type == CARTA::RegionType::POLYGON); diff --git a/src/Region/RegionHandler.h b/src/Region/RegionHandler.h index 931919c2d..529a01d4f 100644 --- a/src/Region/RegionHandler.h +++ b/src/Region/RegionHandler.h @@ -35,7 +35,6 @@ class RegionHandler { // Regions bool SetRegion(int& region_id, RegionState& region_state, std::shared_ptr csys); - bool RegionChanged(int region_id); void RemoveRegion(int region_id); std::shared_ptr GetRegion(int region_id); diff --git a/src/Region/RegionImportExport.cc b/src/Region/RegionImportExport.cc index a9d21dc85..cc6244341 100644 --- a/src/Region/RegionImportExport.cc +++ b/src/Region/RegionImportExport.cc @@ -380,13 +380,30 @@ bool RegionImportExport::ConvertRecordToRectangle( // Convert casacore Record to box Quantity control points. // Rectangles are exported to Record as LCPolygon with 4 points: blc, brc, trc, tlc. // The input Record for a rotbox must be the corners of an unrotated box (rotation in the region state) - casacore::Vector x = region_record.asArrayFloat("x"); - casacore::Vector y = region_record.asArrayFloat("y"); + casacore::Vector x, y; + + if (region_record.dataType("x") == casacore::TpArrayFloat) { + casacore::Vector xf, yf; + xf = region_record.asArrayFloat("x"); + yf = region_record.asArrayFloat("y"); + + // Convert to Double + auto xf_size(xf.size()); + x.resize(xf_size); + y.resize(xf_size); + for (auto i = 0; i < xf_size; ++i) { + x(i) = xf(i); + y(i) = yf(i); + } + } else { + x = region_record.asArrayDouble("x"); + y = region_record.asArrayDouble("y"); + } // Make zero-based if (region_record.asBool("oneRel")) { - x -= (float)1.0; - y -= (float)1.0; + x -= 1.0; + y -= 1.0; } double cx, cy, width, height; @@ -445,7 +462,7 @@ bool RegionImportExport::ConvertRecordToEllipse(const RegionState& region_state, // RegionState needed to check if bmaj/bmin swapped for LCEllipsoid casacore::Vector center = region_record.asArrayFloat("center"); casacore::Vector radii = region_record.asArrayFloat("radii"); - casacore::Float theta = region_record.asFloat("theta"); // radians + casacore::Double theta = region_record.asDouble("theta"); // radians rotation = casacore::Quantity(theta, "rad"); rotation.convert("deg"); // CASA rotang, from x-axis diff --git a/src/Region/RegionState.h b/src/Region/RegionState.h new file mode 100644 index 000000000..30495e039 --- /dev/null +++ b/src/Region/RegionState.h @@ -0,0 +1,149 @@ +/* This file is part of the CARTA Image Viewer: https://github.com/CARTAvis/carta-backend + Copyright 2018- Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), + Associated Universities, Inc. (AUI) and the Inter-University Institute for Data Intensive Astronomy (IDIA) + SPDX-License-Identifier: GPL-3.0-or-later +*/ + +//# RegionState.h: struct for holding region parameters + +#ifndef CARTA_SRC_REGION_REGIONSTATE_H_ +#define CARTA_SRC_REGION_REGIONSTATE_H_ + +#include + +namespace carta { + +struct RegionState { + // struct used for region parameters + int reference_file_id; + CARTA::RegionType type; + std::vector control_points; + float rotation; + + RegionState() : reference_file_id(-1), type(CARTA::POINT), rotation(0) {} + RegionState(int ref_file_id_, CARTA::RegionType type_, const std::vector& control_points_, float rotation_) + : reference_file_id(ref_file_id_), type(type_), control_points(control_points_), rotation(rotation_) {} + + bool operator==(const RegionState& rhs) { + return (reference_file_id == rhs.reference_file_id) && (type == rhs.type) && !RegionChanged(rhs); + } + + bool operator!=(const RegionState& rhs) { + return (reference_file_id != rhs.reference_file_id) || (type != rhs.type) || RegionChanged(rhs); + } + + bool RegionDefined() { + return !control_points.empty(); + } + + bool RegionChanged(const RegionState& rhs) { + // Ignores annotation params (for interrupting region calculations) + return (rotation != rhs.rotation) || PointsChanged(rhs); + } + bool PointsChanged(const RegionState& rhs) { + // Points must be same size, order, and value to be unchanged + if (control_points.size() != rhs.control_points.size()) { + return true; + } + for (int i = 0; i < control_points.size(); ++i) { + float x(control_points[i].x()), y(control_points[i].y()); + float rhs_x(rhs.control_points[i].x()), rhs_y(rhs.control_points[i].y()); + if ((x != rhs_x) || (y != rhs_y)) { + return true; + } + } + return false; + } + + bool IsPoint() { + // Includes annotation types. + return type == CARTA::POINT || type == CARTA::ANNPOINT; + } + + bool IsLineType() { + // Not enclosed region, defined by points. Includes annotation types. + std::vector line_types{ + CARTA::LINE, CARTA::POLYLINE, CARTA::ANNLINE, CARTA::ANNPOLYLINE, CARTA::ANNVECTOR, CARTA::ANNRULER}; + return std::find(line_types.begin(), line_types.end(), type) != line_types.end(); + } + + bool IsBox() { + // Rectangle-type regions. Includes annotation types. + std::vector rectangle_types{CARTA::RECTANGLE, CARTA::ANNRECTANGLE, CARTA::ANNTEXT}; + return std::find(rectangle_types.begin(), rectangle_types.end(), type) != rectangle_types.end(); + } + + bool IsRotbox() { + // Rectangle-type regions with rotation. Includes annotation types. + return IsBox() && (rotation != 0.0); + } + + bool IsAnnotation() { + return type > CARTA::POLYGON; + } + + bool GetRectangleCorners(casacore::Vector& x, casacore::Vector& y, bool apply_rotation = true) { + // Convert rectangle points [[cx, cy], [width, height]] to corner points. Optionally apply rotation. + if (type != CARTA::RECTANGLE && type != CARTA::ANNRECTANGLE && type != CARTA::ANNTEXT) { + return false; + } + float center_x(control_points[0].x()), center_y(control_points[0].y()); + float width(control_points[1].x()), height(control_points[1].y()); + x.resize(4); + y.resize(4); + + if (rotation == 0.0 || !apply_rotation) { + float x_min(center_x - width / 2.0f), x_max(center_x + width / 2.0f); + float y_min(center_y - height / 2.0f), y_max(center_y + height / 2.0f); + // Bottom left + x(0) = x_min; + y(0) = y_min; + // Bottom right + x(1) = x_max; + y(1) = y_min; + // Top right + x(2) = x_max; + y(2) = y_max; + // Top left + x(3) = x_min; + y(3) = y_max; + } else { + // Apply rotation matrix to get width and height vectors in rotated basis + float cos_x = cos(rotation * M_PI / 180.0f); + float sin_x = sin(rotation * M_PI / 180.0f); + float width_vector_x = cos_x * width; + float width_vector_y = sin_x * width; + float height_vector_x = -sin_x * height; + float height_vector_y = cos_x * height; + + // Bottom left + x(0) = center_x + (-width_vector_x - height_vector_x) / 2.0f; + y(0) = center_y + (-width_vector_y - height_vector_y) / 2.0f; + // Bottom right + x(1) = center_x + (width_vector_x - height_vector_x) / 2.0f; + y(1) = center_y + (width_vector_y - height_vector_y) / 2.0f; + // Top right + x(2) = center_x + (width_vector_x + height_vector_x) / 2.0f; + y(2) = center_y + (width_vector_y + height_vector_y) / 2.0f; + // Top left + x(3) = center_x + (-width_vector_x + height_vector_x) / 2.0f; + y(3) = center_y + (-width_vector_y + height_vector_y) / 2.0f; + } + return true; + } + + std::string GetLineRegionName() { + // Names not defined in casacore, for region record + std::string name; + if (IsLineType()) { + std::unordered_map line_region_names = {{CARTA::LINE, "line"}, {CARTA::POLYLINE, "polyline"}, + {CARTA::ANNLINE, "line"}, {CARTA::ANNPOLYLINE, "polyline"}, {CARTA::ANNVECTOR, "vector"}, {CARTA::ANNRULER, "ruler"}}; + name = line_region_names.at(type); + } + return name; + } +}; + +} // namespace carta + +#endif // CARTA_SRC_REGION_REGIONSTATE_H_ diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index d9f149279..be1972fa7 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -39,6 +39,7 @@ set(TEST_SOURCES CommonTestUtilities.cc TestBlockSmooth.cc TestContour.cc + TestCursorSpatialProfiles.cc TestExprImage.cc TestFileInfo.cc TestFileList.cc @@ -48,14 +49,19 @@ set(TEST_SOURCES TestHdf5Image.cc TestHistogram.cc TestImageFitting.cc - TestLineSpatialProfiles.cc TestMain.cc TestMoment.cc TestNormalizedUnits.cc TestProgramSettings.cc TestPvGenerator.cc + TestRegion.cc + TestRegionImportExport.cc + TestRegionHistogram.cc + TestRegionMatched.cc + TestRegionSpatialProfiles.cc + TestRegionSpectralProfiles.cc + TestRegionStats.cc TestRestApi.cc - TestSpatialProfiles.cc TestTileEncoding.cc TestUtil.cc TestVoTable.cc) diff --git a/test/CommonTestUtilities.cc b/test/CommonTestUtilities.cc index 315689423..768049d51 100644 --- a/test/CommonTestUtilities.cc +++ b/test/CommonTestUtilities.cc @@ -386,27 +386,8 @@ void CmpSpatialProfiles( const std::vector& data_vec, const std::pair, std::vector>& data_profiles) { EXPECT_EQ(data_vec.size(), 1); for (const auto& data : data_vec) { - CmpVectors(GetSpatialProfileValues(data.profiles(0)), data_profiles.first); - CmpVectors(GetSpatialProfileValues(data.profiles(1)), data_profiles.second); - } -} - -void CmpVectors(const std::vector& data1, const std::vector& data2, float abs_err) { - EXPECT_EQ(data1.size(), data2.size()); - if (data1.size() == data2.size()) { - for (int i = 0; i < data1.size(); ++i) { - CmpValues(data1[i], data2[i], abs_err); - } - } -} - -void CmpValues(float data1, float data2, float abs_err) { - if (!std::isnan(data1) || !std::isnan(data2)) { - if (abs_err > 0) { - EXPECT_NEAR(data1, data2, abs_err); - } else { - EXPECT_FLOAT_EQ(data1, data2); - } + CmpVectors(GetSpatialProfileValues(data.profiles(0)), data_profiles.first); + CmpVectors(GetSpatialProfileValues(data.profiles(1)), data_profiles.second); } } diff --git a/test/CommonTestUtilities.h b/test/CommonTestUtilities.h index 1944aee5a..b3f5c8804 100644 --- a/test/CommonTestUtilities.h +++ b/test/CommonTestUtilities.h @@ -92,12 +92,15 @@ void GetImageData(std::vector& data, std::shared_ptr std::vector GetSpectralProfileValues(const CARTA::SpectralProfile& profile); std::vector GetSpatialProfileValues(const CARTA::SpatialProfile& profile); -void CmpValues(float data1, float data2, float abs_err = 0); -void CmpVectors(const std::vector& data1, const std::vector& data2, float abs_err = 0); void CmpSpatialProfiles( const std::vector& data_vec, const std::pair, std::vector>& data_profiles); bool CmpHistograms(const carta::Histogram& hist1, const carta::Histogram& hist2); +template +void CmpValues(T data1, T data2, T abs_err = 0); +template +void CmpVectors(const std::vector& data1, const std::vector& data2, T abs_err = 0); + #include "CommonTestUtilities.tcc" #endif // CARTA_TEST_COMMONTESTUTILITIES_H_ diff --git a/test/CommonTestUtilities.tcc b/test/CommonTestUtilities.tcc index 31778eb89..a3fb02242 100644 --- a/test/CommonTestUtilities.tcc +++ b/test/CommonTestUtilities.tcc @@ -20,4 +20,21 @@ std::vector GetSpectralProfileValues(const CARTA::SpectralProfile& profile) { return values; } +template +void CmpVectors(const std::vector& data1, const std::vector& data2, T abs_err) { + EXPECT_EQ(data1.size(), data2.size()); + if (data1.size() == data2.size()) { + for (int i = 0; i < data1.size(); ++i) { + CmpValues(data1[i], data2[i], abs_err); + } + } +} + +template +void CmpValues(T data1, T data2, T abs_err) { + if (!std::isnan(data1) || !std::isnan(data2)) { + EXPECT_NEAR(data1, data2, abs_err); + } +} + #endif // CARTA_TEST_COMMONTESTUTILITIES_TCC_ diff --git a/test/TestSpatialProfiles.cc b/test/TestCursorSpatialProfiles.cc similarity index 88% rename from test/TestSpatialProfiles.cc rename to test/TestCursorSpatialProfiles.cc index 8a532bd1f..8aa36f078 100644 --- a/test/TestSpatialProfiles.cc +++ b/test/TestCursorSpatialProfiles.cc @@ -19,7 +19,7 @@ using ::testing::Pointwise; static const std::string IMAGE_OPTS = "-s 0 -n row column -d 10"; -class SpatialProfileTest : public ::testing::Test, public ImageGenerator { +class CursorSpatialProfileTest : public ::testing::Test, public ImageGenerator { public: static std::tuple GetProfiles(CARTA::SpatialProfileData& data) { if (data.profiles(0).coordinate().back() == 'x') { @@ -99,7 +99,7 @@ class SpatialProfileTest : public ::testing::Test, public ImageGenerator { } }; -TEST_F(SpatialProfileTest, SmallFitsProfile) { +TEST_F(CursorSpatialProfileTest, SmallFitsProfile) { auto path_string = GeneratedFitsImagePath("10 10", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -129,18 +129,18 @@ TEST_F(SpatialProfileTest, SmallFitsProfile) { EXPECT_EQ(x_profile.mip(), 0); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 10); - CmpVectors(x_vals, reader.ReadProfileX(5)); + CmpVectors(x_vals, reader.ReadProfileX(5)); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 10); EXPECT_EQ(y_profile.mip(), 0); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 10); - CmpVectors(y_vals, reader.ReadProfileY(5)); + CmpVectors(y_vals, reader.ReadProfileY(5)); } } -TEST_F(SpatialProfileTest, SmallHdf5Profile) { +TEST_F(CursorSpatialProfileTest, SmallHdf5Profile) { auto path_string = GeneratedHdf5ImagePath("10 10", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -170,18 +170,18 @@ TEST_F(SpatialProfileTest, SmallHdf5Profile) { EXPECT_EQ(x_profile.mip(), 0); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 10); - CmpVectors(x_vals, reader.ReadProfileX(5)); + CmpVectors(x_vals, reader.ReadProfileX(5)); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 10); EXPECT_EQ(y_profile.mip(), 0); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 10); - CmpVectors(y_vals, reader.ReadProfileY(5)); + CmpVectors(y_vals, reader.ReadProfileY(5)); } } -TEST_F(SpatialProfileTest, LowResFitsProfile) { +TEST_F(CursorSpatialProfileTest, LowResFitsProfile) { auto path_string = GeneratedFitsImagePath("130 100", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -205,18 +205,18 @@ TEST_F(SpatialProfileTest, LowResFitsProfile) { EXPECT_EQ(x_profile.mip(), 2); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 66); - CmpVectors(x_vals, Decimated(reader.ReadProfileX(50), 2)); + CmpVectors(x_vals, Decimated(reader.ReadProfileX(50), 2)); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 100); EXPECT_EQ(y_profile.mip(), 2); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 50); - CmpVectors(y_vals, Decimated(reader.ReadProfileY(50), 2)); + CmpVectors(y_vals, Decimated(reader.ReadProfileY(50), 2)); } } -TEST_F(SpatialProfileTest, LowResHdf5ProfileExactMipAvailable) { +TEST_F(CursorSpatialProfileTest, LowResHdf5ProfileExactMipAvailable) { auto path_string = GeneratedHdf5ImagePath("130 100", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -240,18 +240,18 @@ TEST_F(SpatialProfileTest, LowResHdf5ProfileExactMipAvailable) { EXPECT_EQ(x_profile.mip(), 2); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 65); - CmpVectors(x_vals, Downsampled({reader.ReadProfileX(50), reader.ReadProfileX(51)}), 1e-5); + CmpVectors(x_vals, Downsampled({reader.ReadProfileX(50), reader.ReadProfileX(51)}), 1e-5); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 100); EXPECT_EQ(y_profile.mip(), 2); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 50); - CmpVectors(y_vals, Downsampled({reader.ReadProfileY(50), reader.ReadProfileY(51)}), 1e-5); + CmpVectors(y_vals, Downsampled({reader.ReadProfileY(50), reader.ReadProfileY(51)}), 1e-5); } } -TEST_F(SpatialProfileTest, LowResHdf5ProfileLowerMipAvailable) { +TEST_F(CursorSpatialProfileTest, LowResHdf5ProfileLowerMipAvailable) { auto path_string = GeneratedHdf5ImagePath("130 100", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -277,18 +277,18 @@ TEST_F(SpatialProfileTest, LowResHdf5ProfileLowerMipAvailable) { EXPECT_EQ(x_profile.mip(), 2); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 65); - CmpVectors(x_vals, Downsampled({reader.ReadProfileX(50), reader.ReadProfileX(51)}), 1e-5); + CmpVectors(x_vals, Downsampled({reader.ReadProfileX(50), reader.ReadProfileX(51)}), 1e-5); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 100); EXPECT_EQ(y_profile.mip(), 2); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 50); - CmpVectors(y_vals, Downsampled({reader.ReadProfileY(50), reader.ReadProfileY(51)}), 1e-5); + CmpVectors(y_vals, Downsampled({reader.ReadProfileY(50), reader.ReadProfileY(51)}), 1e-5); } } -TEST_F(SpatialProfileTest, LowResHdf5ProfileNoMipAvailable) { +TEST_F(CursorSpatialProfileTest, LowResHdf5ProfileNoMipAvailable) { auto path_string = GeneratedHdf5ImagePath("120 100", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -314,18 +314,18 @@ TEST_F(SpatialProfileTest, LowResHdf5ProfileNoMipAvailable) { EXPECT_EQ(x_profile.mip(), 2); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 60); - CmpVectors(x_vals, Decimated(reader.ReadProfileX(50), 2)); + CmpVectors(x_vals, Decimated(reader.ReadProfileX(50), 2)); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 100); EXPECT_EQ(y_profile.mip(), 2); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 50); - CmpVectors(y_vals, Decimated(reader.ReadProfileY(50), 2)); + CmpVectors(y_vals, Decimated(reader.ReadProfileY(50), 2)); } } -TEST_F(SpatialProfileTest, FullResFitsStartEnd) { +TEST_F(CursorSpatialProfileTest, FullResFitsStartEnd) { auto path_string = GeneratedFitsImagePath("400 300", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -349,18 +349,18 @@ TEST_F(SpatialProfileTest, FullResFitsStartEnd) { EXPECT_EQ(x_profile.mip(), 0); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 100); - CmpVectors(x_vals, Segment(reader.ReadProfileX(150), 100, 200)); + CmpVectors(x_vals, Segment(reader.ReadProfileX(150), 100, 200)); EXPECT_EQ(y_profile.start(), 100); EXPECT_EQ(y_profile.end(), 200); EXPECT_EQ(y_profile.mip(), 0); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 100); - CmpVectors(y_vals, Segment(reader.ReadProfileY(150), 100, 200)); + CmpVectors(y_vals, Segment(reader.ReadProfileY(150), 100, 200)); } } -TEST_F(SpatialProfileTest, FullResHdf5StartEnd) { +TEST_F(CursorSpatialProfileTest, FullResHdf5StartEnd) { auto path_string = GeneratedHdf5ImagePath("400 300", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -384,18 +384,18 @@ TEST_F(SpatialProfileTest, FullResHdf5StartEnd) { EXPECT_EQ(x_profile.mip(), 0); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 100); - CmpVectors(x_vals, Segment(reader.ReadProfileX(150), 100, 200)); + CmpVectors(x_vals, Segment(reader.ReadProfileX(150), 100, 200)); EXPECT_EQ(y_profile.start(), 100); EXPECT_EQ(y_profile.end(), 200); EXPECT_EQ(y_profile.mip(), 0); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 100); - CmpVectors(y_vals, Segment(reader.ReadProfileY(150), 100, 200)); + CmpVectors(y_vals, Segment(reader.ReadProfileY(150), 100, 200)); } } -TEST_F(SpatialProfileTest, LowResFitsStartEnd) { +TEST_F(CursorSpatialProfileTest, LowResFitsStartEnd) { auto path_string = GeneratedFitsImagePath("400 300", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -420,7 +420,7 @@ TEST_F(SpatialProfileTest, LowResFitsStartEnd) { auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 24); // Data to decimate has endpoints rounded up to mip*2 - CmpVectors(x_vals, Decimated(Segment(reader.ReadProfileX(150), 104, 200), 4)); + CmpVectors(x_vals, Decimated(Segment(reader.ReadProfileX(150), 104, 200), 4)); EXPECT_EQ(y_profile.start(), 100); EXPECT_EQ(y_profile.end(), 200); @@ -428,11 +428,11 @@ TEST_F(SpatialProfileTest, LowResFitsStartEnd) { auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 24); // Data to decimate has endpoints rounded up to mip*2 - CmpVectors(y_vals, Decimated(Segment(reader.ReadProfileY(150), 104, 200), 4)); + CmpVectors(y_vals, Decimated(Segment(reader.ReadProfileY(150), 104, 200), 4)); } } -TEST_F(SpatialProfileTest, LowResHdf5StartEnd) { +TEST_F(CursorSpatialProfileTest, LowResHdf5StartEnd) { auto path_string = GeneratedHdf5ImagePath("400 300", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -457,7 +457,7 @@ TEST_F(SpatialProfileTest, LowResHdf5StartEnd) { auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 25); // Downsampled region is selected so that it includes the requested row - CmpVectors(x_vals, + CmpVectors(x_vals, Segment(Downsampled({reader.ReadProfileX(148), reader.ReadProfileX(149), reader.ReadProfileX(150), reader.ReadProfileX(151)}), 25, 50), 1e-5); @@ -468,14 +468,14 @@ TEST_F(SpatialProfileTest, LowResHdf5StartEnd) { auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 25); // Downsampled region is selected so that it includes the requested column - CmpVectors(y_vals, + CmpVectors(y_vals, Segment(Downsampled({reader.ReadProfileY(148), reader.ReadProfileY(149), reader.ReadProfileY(150), reader.ReadProfileY(151)}), 25, 50), 1e-5); } } -TEST_F(SpatialProfileTest, Hdf5MultipleChunkFullRes) { +TEST_F(CursorSpatialProfileTest, Hdf5MultipleChunkFullRes) { auto path_string = GeneratedHdf5ImagePath("3000 2000", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -498,18 +498,18 @@ TEST_F(SpatialProfileTest, Hdf5MultipleChunkFullRes) { EXPECT_EQ(x_profile.mip(), 0); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 3000); - CmpVectors(x_vals, reader.ReadProfileX(150)); + CmpVectors(x_vals, reader.ReadProfileX(150)); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 2000); EXPECT_EQ(y_profile.mip(), 0); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 2000); - CmpVectors(y_vals, reader.ReadProfileY(150)); + CmpVectors(y_vals, reader.ReadProfileY(150)); } } -TEST_F(SpatialProfileTest, Hdf5MultipleChunkFullResStartEnd) { +TEST_F(CursorSpatialProfileTest, Hdf5MultipleChunkFullResStartEnd) { auto path_string = GeneratedHdf5ImagePath("3000 2000", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -533,18 +533,18 @@ TEST_F(SpatialProfileTest, Hdf5MultipleChunkFullResStartEnd) { EXPECT_EQ(x_profile.mip(), 0); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 500); - CmpVectors(x_vals, Segment(reader.ReadProfileX(1250), 1000, 1500)); + CmpVectors(x_vals, Segment(reader.ReadProfileX(1250), 1000, 1500)); EXPECT_EQ(y_profile.start(), 1000); EXPECT_EQ(y_profile.end(), 1500); EXPECT_EQ(y_profile.mip(), 0); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 500); - CmpVectors(y_vals, Segment(reader.ReadProfileY(1250), 1000, 1500)); + CmpVectors(y_vals, Segment(reader.ReadProfileY(1250), 1000, 1500)); } } -TEST_F(SpatialProfileTest, FitsChannelChange) { +TEST_F(CursorSpatialProfileTest, FitsChannelChange) { auto path_string = GeneratedFitsImagePath("10 10 2", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -576,18 +576,18 @@ TEST_F(SpatialProfileTest, FitsChannelChange) { EXPECT_EQ(x_profile.mip(), 0); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 10); - CmpVectors(x_vals, reader.ReadProfileX(5, 1)); + CmpVectors(x_vals, reader.ReadProfileX(5, 1)); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 10); EXPECT_EQ(y_profile.mip(), 0); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 10); - CmpVectors(y_vals, reader.ReadProfileY(5, 1)); + CmpVectors(y_vals, reader.ReadProfileY(5, 1)); } } -TEST_F(SpatialProfileTest, FitsChannelStokesChange) { +TEST_F(CursorSpatialProfileTest, FitsChannelStokesChange) { auto path_string = GeneratedFitsImagePath("10 10 2 2", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -625,18 +625,18 @@ TEST_F(SpatialProfileTest, FitsChannelStokesChange) { EXPECT_EQ(x_profile.mip(), 0); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 10); - CmpVectors(x_vals, reader.ReadProfileX(y, channel, spatial_config_stokes)); + CmpVectors(x_vals, reader.ReadProfileX(y, channel, spatial_config_stokes)); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 10); EXPECT_EQ(y_profile.mip(), 0); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 10); - CmpVectors(y_vals, reader.ReadProfileY(x, channel, spatial_config_stokes)); + CmpVectors(y_vals, reader.ReadProfileY(x, channel, spatial_config_stokes)); } } -TEST_F(SpatialProfileTest, ContiguousHDF5ChannelChange) { +TEST_F(CursorSpatialProfileTest, ContiguousHDF5ChannelChange) { auto path_string = GeneratedHdf5ImagePath("10 10 2", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -668,18 +668,18 @@ TEST_F(SpatialProfileTest, ContiguousHDF5ChannelChange) { EXPECT_EQ(x_profile.mip(), 0); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 10); - CmpVectors(x_vals, reader.ReadProfileX(5, 1)); + CmpVectors(x_vals, reader.ReadProfileX(5, 1)); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 10); EXPECT_EQ(y_profile.mip(), 0); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 10); - CmpVectors(y_vals, reader.ReadProfileY(5, 1)); + CmpVectors(y_vals, reader.ReadProfileY(5, 1)); } } -TEST_F(SpatialProfileTest, ChunkedHDF5ChannelChange) { +TEST_F(CursorSpatialProfileTest, ChunkedHDF5ChannelChange) { auto path_string = GeneratedHdf5ImagePath("1000 1000 2", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -711,18 +711,18 @@ TEST_F(SpatialProfileTest, ChunkedHDF5ChannelChange) { EXPECT_EQ(x_profile.mip(), 0); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 1000); - CmpVectors(x_vals, reader.ReadProfileX(5, 1)); + CmpVectors(x_vals, reader.ReadProfileX(5, 1)); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 1000); EXPECT_EQ(y_profile.mip(), 0); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 1000); - CmpVectors(y_vals, reader.ReadProfileY(5, 1)); + CmpVectors(y_vals, reader.ReadProfileY(5, 1)); } } -TEST_F(SpatialProfileTest, ChunkedHDF5ChannelStokesChange) { +TEST_F(CursorSpatialProfileTest, ChunkedHDF5ChannelStokesChange) { auto path_string = GeneratedHdf5ImagePath("1000 1000 2 2", IMAGE_OPTS); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::unique_ptr frame(new Frame(0, loader, "0")); @@ -760,13 +760,13 @@ TEST_F(SpatialProfileTest, ChunkedHDF5ChannelStokesChange) { EXPECT_EQ(x_profile.mip(), 0); auto x_vals = ProfileValues(x_profile); EXPECT_EQ(x_vals.size(), 1000); - CmpVectors(x_vals, reader.ReadProfileX(y, channel, spatial_config_stokes)); + CmpVectors(x_vals, reader.ReadProfileX(y, channel, spatial_config_stokes)); EXPECT_EQ(y_profile.start(), 0); EXPECT_EQ(y_profile.end(), 1000); EXPECT_EQ(y_profile.mip(), 0); auto y_vals = ProfileValues(y_profile); EXPECT_EQ(y_vals.size(), 1000); - CmpVectors(y_vals, reader.ReadProfileY(x, channel, spatial_config_stokes)); + CmpVectors(y_vals, reader.ReadProfileY(x, channel, spatial_config_stokes)); } } diff --git a/test/TestExprImage.cc b/test/TestExprImage.cc index 420bff2e2..dd85e8fbe 100644 --- a/test/TestExprImage.cc +++ b/test/TestExprImage.cc @@ -74,10 +74,10 @@ class ImageExprTest : public ::testing::Test { EXPECT_EQ(image_shape, expr_shape); // Compare image xprofile * 2 to expr xprofile for_each(image_xprofile.begin(), image_xprofile.end(), [](float& a) { a *= 2; }); - CmpVectors(image_xprofile, expr_xprofile.tovector()); + CmpVectors(image_xprofile, expr_xprofile.tovector()); // Compare image yprofile * 2 to expr yprofile for_each(image_yprofile.begin(), image_yprofile.end(), [](float& a) { a *= 2; }); - CmpVectors(image_yprofile, expr_yprofile.tovector()); + CmpVectors(image_yprofile, expr_yprofile.tovector()); } void SaveImageExpr(const std::string& file_name, const std::string& hdu, CARTA::FileType file_type) { diff --git a/test/TestLineSpatialProfiles.cc b/test/TestLineSpatialProfiles.cc deleted file mode 100644 index 9da319f0a..000000000 --- a/test/TestLineSpatialProfiles.cc +++ /dev/null @@ -1,222 +0,0 @@ -/* This file is part of the CARTA Image Viewer: https://github.com/CARTAvis/carta-backend - Copyright 2018- Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), - Associated Universities, Inc. (AUI) and the Inter-University Institute for Data Intensive Astronomy (IDIA) - SPDX-License-Identifier: GPL-3.0-or-later -*/ - -#include -#include - -#include "CommonTestUtilities.h" -#include "ImageData/FileLoader.h" -#include "Region/Region.h" -#include "Region/RegionHandler.h" -#include "src/Frame/Frame.h" - -using namespace carta; - -using ::testing::FloatNear; -using ::testing::Pointwise; - -class LineSpatialProfileTest : public ::testing::Test { -public: - static bool SetLineRegion(carta::RegionHandler& region_handler, int file_id, int& region_id, const std::vector& endpoints, - std::shared_ptr csys, bool is_annotation) { - std::vector control_points; - for (auto i = 0; i < endpoints.size(); i += 2) { - control_points.push_back(Message::Point(endpoints[i], endpoints[i + 1])); - } - - // Define RegionState for line region and set region (region_id updated) - auto npoints(control_points.size()); - CARTA::RegionType region_type(CARTA::RegionType::LINE); - if (is_annotation) { - if (npoints > 2) { - region_type = CARTA::RegionType::ANNPOLYLINE; - } else { - region_type = CARTA::RegionType::ANNLINE; - } - } else if (npoints > 2) { - region_type = CARTA::RegionType::POLYLINE; - } - RegionState region_state(file_id, region_type, control_points, 0.0); - return region_handler.SetRegion(region_id, region_state, csys); - } - - static bool GetLineProfiles(const std::string& image_path, const std::vector& endpoints, - const std::vector& spatial_reqs, CARTA::SpatialProfileData& spatial_profile, - bool is_annotation = false) { - std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); - std::shared_ptr frame(new Frame(0, loader, "0")); - carta::RegionHandler region_handler; - - // Set line region - int file_id(0), region_id(-1); - auto csys = frame->CoordinateSystem(); - if (!SetLineRegion(region_handler, file_id, region_id, endpoints, csys, is_annotation)) { - return false; - } - - // Get spatial profiles - if (!region_handler.SetSpatialRequirements(region_id, file_id, frame, spatial_reqs)) { - return false; - } - - return region_handler.FillLineSpatialProfileData( - file_id, region_id, [&](CARTA::SpatialProfileData profile_data) { spatial_profile = profile_data; }); - } - - static std::vector ProfileValues(CARTA::SpatialProfile& profile) { - std::string buffer = profile.raw_values_fp32(); - std::vector values(buffer.size() / sizeof(float)); - memcpy(values.data(), buffer.data(), buffer.size()); - return values; - } - - void SetUp() { - setenv("HDF5_USE_FILE_LOCKING", "FALSE", 0); - } - - static void TestAveragingWidthRange(int width, bool expected_width_range) { - std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); - std::vector endpoints = {0.0, 0.0, 9.0, 9.0}; - int start(0), end(0), mip(0); - std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; - CARTA::SpatialProfileData spatial_profile; - - if (expected_width_range) { - ASSERT_TRUE(GetLineProfiles(image_path, endpoints, spatial_reqs, spatial_profile)); - ASSERT_EQ(spatial_profile.profiles_size(), 1); - } else { - ASSERT_FALSE(GetLineProfiles(image_path, endpoints, spatial_reqs, spatial_profile)); - ASSERT_EQ(spatial_profile.profiles_size(), 0); - } - } -}; - -TEST_F(LineSpatialProfileTest, FitsLineProfile) { - std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); - std::vector endpoints = {0.0, 0.0, 9.0, 9.0}; - int start(0), end(0), mip(0), width(3); - std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; - - CARTA::SpatialProfileData spatial_profile; - bool ok = GetLineProfiles(image_path, endpoints, spatial_reqs, spatial_profile); - - ASSERT_TRUE(ok); - ASSERT_EQ(spatial_profile.profiles_size(), 1); -} - -TEST_F(LineSpatialProfileTest, Hdf5LineProfile) { - std::string image_path = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); - std::vector endpoints = {0.0, 0.0, 9.0, 9.0}; - int start(0), end(0), mip(0), width(3); - std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; - - CARTA::SpatialProfileData spatial_profile; - bool ok = GetLineProfiles(image_path, endpoints, spatial_reqs, spatial_profile); - - ASSERT_TRUE(ok); - ASSERT_EQ(spatial_profile.profiles_size(), 1); -} - -TEST_F(LineSpatialProfileTest, FitsHorizontalCutProfile) { - std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); - std::vector endpoints = {9.0, 5.0, 1.0, 5.0}; // Set line region at y=5 - int start(0), end(0), mip(0), width(1); - std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; - - CARTA::SpatialProfileData spatial_profile; - bool ok = GetLineProfiles(image_path, endpoints, spatial_reqs, spatial_profile); - - ASSERT_TRUE(ok); - ASSERT_EQ(spatial_profile.profiles_size(), 1); - - // Profile data - auto profile = spatial_profile.profiles(0); - std::vector profile_data = ProfileValues(profile); - EXPECT_EQ(profile_data.size(), 9); - - // Read image data slice for first channel - FitsDataReader reader(image_path); - auto image_data = reader.ReadRegion({1, 5, 0}, {10, 6, 1}); - - // Profile data width=1 of horizontal line is same as slice - CmpVectors(profile_data, image_data); -} - -TEST_F(LineSpatialProfileTest, FitsVerticalCutProfile) { - std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); - std::vector endpoints = {5.0, 9.0, 5.0, 1.0}; // Set line region at x=5 - int start(0), end(0), mip(0), width(1); - std::vector spatial_reqs = {Message::SpatialConfig("y", start, end, mip, width)}; - - CARTA::SpatialProfileData spatial_profile; - bool ok = GetLineProfiles(image_path, endpoints, spatial_reqs, spatial_profile); - - ASSERT_TRUE(ok); - ASSERT_EQ(spatial_profile.profiles_size(), 1); - - // Profile data - auto profile = spatial_profile.profiles(0); - std::vector profile_data = ProfileValues(profile); - EXPECT_EQ(profile_data.size(), 9); - - // Read image data slice for first channel - FitsDataReader reader(image_path); - auto image_data = reader.ReadRegion({5, 1, 0}, {6, 10, 1}); - - // Profile data width=1 of horizontal line is same as slice - CmpVectors(profile_data, image_data); -} - -TEST_F(LineSpatialProfileTest, FitsPolylineProfile) { - std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); - std::vector endpoints = {1.0, 1.0, 9.0, 1.0, 9.0, 5.0}; - int start(0), end(0), mip(0), width(1); - std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; - std::vector spatial_profiles; - - CARTA::SpatialProfileData spatial_profile; - bool ok = GetLineProfiles(image_path, endpoints, spatial_reqs, spatial_profile); - ASSERT_TRUE(ok); - ASSERT_EQ(spatial_profile.profiles_size(), 1); - - // Profile data - auto profile = spatial_profile.profiles(0); - std::vector profile_data = ProfileValues(profile); - EXPECT_EQ(profile_data.size(), 13); - - // Read image data slice for first channel - FitsDataReader reader(image_path); - auto line0_data = reader.ReadRegion({1, 1, 0}, {10, 2, 1}); - // Trim line 1: [9, 1] already covered by line 0 - auto line1_data = reader.ReadRegion({9, 2, 0}, {10, 6, 1}); - auto image_data = line0_data; - for (size_t i = 0; i < line1_data.size(); ++i) { - image_data.push_back(line1_data[i]); - } - - // Profile data width=1 of polyline is same as slices - CmpVectors(profile_data, image_data); -} - -TEST_F(LineSpatialProfileTest, AveragingWidthRange) { - TestAveragingWidthRange(0, false); - TestAveragingWidthRange(1, true); - TestAveragingWidthRange(20, true); - TestAveragingWidthRange(21, false); -} - -TEST_F(LineSpatialProfileTest, FitsAnnotationLineProfile) { - std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); - std::vector endpoints = {0.0, 0.0, 9.0, 9.0}; - int start(0), end(0), mip(0), width(3); - std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; - - CARTA::SpatialProfileData spatial_profile; - bool is_annotation(true); - bool ok = GetLineProfiles(image_path, endpoints, spatial_reqs, spatial_profile, is_annotation); - - ASSERT_FALSE(ok); -} diff --git a/test/TestPvGenerator.cc b/test/TestPvGenerator.cc index 2af07a321..0e7bcb9fb 100644 --- a/test/TestPvGenerator.cc +++ b/test/TestPvGenerator.cc @@ -39,7 +39,7 @@ class PvGeneratorTest : public ::testing::Test, public ImageGenerator { } static void TestAveragingWidthRange(int width, bool expected_width_range) { - auto image_path = TestRoot() / "data/images/fits/noise_3d.fits"; // 10x10x10 image + auto image_path = FileFinder::FitsImagePath("noise_3d.fits"); // 10x10x10 image std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); std::shared_ptr frame(new Frame(0, loader, "0")); carta::RegionHandler region_handler; @@ -69,7 +69,7 @@ class PvGeneratorTest : public ::testing::Test, public ImageGenerator { }; TEST_F(PvGeneratorTest, FitsPvImage) { - auto image_path = TestRoot() / "data/images/fits/noise_3d.fits"; // 10x10x10 image + auto image_path = FileFinder::FitsImagePath("noise_3d.fits"); // 10x10x10 image std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); std::shared_ptr frame(new Frame(0, loader, "0")); @@ -139,7 +139,7 @@ TEST_F(PvGeneratorTest, FitsPvImage) { } TEST_F(PvGeneratorTest, FitsPvImageHorizontalCut) { - auto image_path = TestRoot() / "data/images/fits/noise_3d.fits"; // 10x10x10 image + auto image_path = FileFinder::FitsImagePath("noise_3d.fits"); // 10x10x10 image std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); std::shared_ptr frame(new Frame(0, loader, "0")); @@ -208,7 +208,7 @@ TEST_F(PvGeneratorTest, FitsPvImageHorizontalCut) { EXPECT_EQ(pv_data.shape()(1), frame->Depth()); // Read image data slice - FitsDataReader reader(image_path.string()); + FitsDataReader reader(image_path); auto image_data = reader.ReadRegion({1, 5, 0}, {10, 6, 10}); EXPECT_EQ(pv_data.size(), image_data.size()); @@ -216,7 +216,7 @@ TEST_F(PvGeneratorTest, FitsPvImageHorizontalCut) { } TEST_F(PvGeneratorTest, FitsPvImageVerticalCut) { - auto image_path = TestRoot() / "data/images/fits/noise_3d.fits"; // 10x10x10 image + auto image_path = FileFinder::FitsImagePath("noise_3d.fits"); // 10x10x10 image std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); std::shared_ptr frame(new Frame(0, loader, "0")); @@ -285,7 +285,7 @@ TEST_F(PvGeneratorTest, FitsPvImageVerticalCut) { EXPECT_EQ(pv_data.shape()(1), frame->Depth()); // Read image data slice - FitsDataReader reader(image_path.string()); + FitsDataReader reader(image_path); auto image_data = reader.ReadRegion({5, 1, 0}, {6, 10, 10}); EXPECT_EQ(pv_data.size(), image_data.size()); @@ -293,7 +293,6 @@ TEST_F(PvGeneratorTest, FitsPvImageVerticalCut) { } TEST_F(PvGeneratorTest, TestNoSpectralAxis) { - auto image_path = TestRoot() / "data/images/hdf5/noise_10px_10px.hdf5"; auto path_string = GeneratedHdf5ImagePath("10 10 10"); std::shared_ptr loader(carta::FileLoader::GetLoader(path_string)); std::shared_ptr frame(new Frame(0, loader, "0")); @@ -327,7 +326,7 @@ TEST_F(PvGeneratorTest, AveragingWidthRange) { TEST_F(PvGeneratorTest, PvImageSpectralRange) { // FITS - auto image_path = TestRoot() / "data/images/fits/noise_3d.fits"; // 10x10x10 image + auto image_path = FileFinder::FitsImagePath("noise_3d.fits"); // 10x10x10 image std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); std::shared_ptr frame(new Frame(0, loader, "0")); auto csys = frame->CoordinateSystem(); @@ -357,7 +356,7 @@ TEST_F(PvGeneratorTest, PvImageSpectralRange) { } TEST_F(PvGeneratorTest, PvImageReversedAxes) { - auto image_path = TestRoot() / "data/images/fits/noise_3d.fits"; // 10x10x10 image + auto image_path = FileFinder::FitsImagePath("noise_3d.fits"); // 10x10x10 image std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); std::shared_ptr frame(new Frame(0, loader, "0")); auto csys = frame->CoordinateSystem(); @@ -397,7 +396,7 @@ TEST_F(PvGeneratorTest, PvImageReversedAxes) { } TEST_F(PvGeneratorTest, PvImageKeep) { - auto image_path = TestRoot() / "data/images/fits/noise_3d.fits"; // 10x10x10 image + auto image_path = FileFinder::FitsImagePath("noise_3d.fits"); // 10x10x10 image std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); std::shared_ptr frame(new Frame(0, loader, "0")); auto csys = frame->CoordinateSystem(); @@ -445,7 +444,7 @@ TEST_F(PvGeneratorTest, PvImageKeep) { } TEST_F(PvGeneratorTest, FitsPvAnnotationLine) { - auto image_path = TestRoot() / "data/images/fits/noise_3d.fits"; // 10x10x10 image + auto image_path = FileFinder::FitsImagePath("noise_3d.fits"); // 10x10x10 image std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); std::shared_ptr frame(new Frame(0, loader, "0")); carta::RegionHandler region_handler; diff --git a/test/TestRegion.cc b/test/TestRegion.cc new file mode 100644 index 000000000..c2a6d7d0f --- /dev/null +++ b/test/TestRegion.cc @@ -0,0 +1,418 @@ +/* This file is part of the CARTA Image Viewer: https://github.com/CARTAvis/carta-backend + Copyright 2018- Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), + Associated Universities, Inc. (AUI) and the Inter-University Institute for Data Intensive Astronomy (IDIA) + SPDX-License-Identifier: GPL-3.0-or-later +*/ + +#include +#include + +#include "CommonTestUtilities.h" +#include "ImageData/FileLoader.h" +#include "Region/Region.h" +#include "Region/RegionHandler.h" +#include "src/Frame/Frame.h" + +using namespace carta; + +class RegionTest : public ::testing::Test { +public: + static bool SetRegion(carta::RegionHandler& region_handler, int file_id, int& region_id, CARTA::RegionType type, + const std::vector& points, float rotation, std::shared_ptr csys) { + std::vector control_points; + for (auto i = 0; i < points.size(); i += 2) { + control_points.push_back(Message::Point(points[i], points[i + 1])); + } + + // Define RegionState and set region (region_id updated) + RegionState region_state(file_id, type, control_points, rotation); + return region_handler.SetRegion(region_id, region_state, csys); + } +}; + +TEST_F(RegionTest, TestSetUpdateRemoveRegion) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::RECTANGLE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(0.0); + auto csys = frame->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + // RegionHandler checks + ASSERT_TRUE(ok); + ASSERT_FALSE(region_handler.IsPointRegion(region_id)); + ASSERT_FALSE(region_handler.IsLineRegion(region_id)); + ASSERT_TRUE(region_handler.IsClosedRegion(region_id)); + + // Region checks + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + ASSERT_TRUE(region->IsValid()); + ASSERT_FALSE(region->IsPoint()); + ASSERT_FALSE(region->IsLineType()); + ASSERT_FALSE(region->IsAnnotation()); + ASSERT_FALSE(region->RegionChanged()); + ASSERT_TRUE(region->IsConnected()); + ASSERT_EQ(region->CoordinateSystem(), csys); + auto region_state = region->GetRegionState(); + ASSERT_FALSE(region_state.IsRotbox()); + + // Update region + rotation = 30.0; + ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + ASSERT_TRUE(region->IsValid()); + ASSERT_TRUE(region->RegionChanged()); + auto new_region_state = region->GetRegionState(); + ASSERT_FALSE(region_state == new_region_state); + ASSERT_TRUE(new_region_state.IsRotbox()); + + // Remove region and frame (not set, should not cause error) + region_handler.RemoveRegion(region_id); + auto no_region = region_handler.GetRegion(region_id); + ASSERT_FALSE(no_region); +} + +TEST_F(RegionTest, TestReferenceImageRectangleLCRegion) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::RECTANGLE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(0.0); + auto csys = frame->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + // Get Region as 3D LCRegion + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + auto image_shape = frame->ImageShape(); + auto lc_region = region->GetImageRegion(file_id, csys, image_shape); + ASSERT_TRUE(lc_region); // shared_ptr + ASSERT_EQ(lc_region->ndim(), 2); + ASSERT_EQ(lc_region->latticeShape()(0), image_shape(0)); + ASSERT_EQ(lc_region->latticeShape()(1), image_shape(1)); + ASSERT_EQ(lc_region->shape(), casacore::IPosition(2, 5, 3)); +} + +TEST_F(RegionTest, TestReferenceImageRotboxLCRegion) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); // 10x10x10 + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::RECTANGLE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(30.0); + auto csys = frame->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + + // Get Region as 3D LCRegion + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + auto image_shape = frame->ImageShape(); + auto lc_region = region->GetImageRegion(file_id, csys, image_shape); + ASSERT_TRUE(lc_region); // shared_ptr + ASSERT_EQ(lc_region->ndim(), 2); + ASSERT_EQ(lc_region->latticeShape()(0), image_shape(0)); + ASSERT_EQ(lc_region->latticeShape()(1), image_shape(1)); + ASSERT_EQ(lc_region->shape(), casacore::IPosition(2, 5, 5)); +} + +TEST_F(RegionTest, TestReferenceImageEllipseLCRegion) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::ELLIPSE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(0.0); + auto csys = frame->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + + // Get Region as 3D LCRegion + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + auto image_shape = frame->ImageShape(); + auto lc_region = region->GetImageRegion(file_id, csys, image_shape); + ASSERT_TRUE(lc_region); // shared_ptr + ASSERT_EQ(lc_region->ndim(), 2); + ASSERT_EQ(lc_region->latticeShape()(0), image_shape(0)); + ASSERT_EQ(lc_region->latticeShape()(1), image_shape(1)); + ASSERT_EQ(lc_region->shape(), casacore::IPosition(2, 7, 9)); +} + +TEST_F(RegionTest, TestReferenceImagePolygonLCRegion) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::POLYGON; + std::vector points = {5.0, 5.0, 4.0, 3.0, 1.0, 6.0, 3.0, 8.0}; + float rotation(0.0); + auto csys = frame->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + + // Get Region as 3D LCRegion + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + auto image_shape = frame->ImageShape(); + auto lc_region = region->GetImageRegion(file_id, csys, image_shape); + ASSERT_TRUE(lc_region); // shared_ptr + ASSERT_EQ(lc_region->ndim(), 2); + ASSERT_EQ(lc_region->latticeShape()(0), image_shape(0)); + ASSERT_EQ(lc_region->latticeShape()(1), image_shape(1)); + ASSERT_EQ(lc_region->shape(), casacore::IPosition(2, 5, 6)); +} + +TEST_F(RegionTest, TestReferenceImagePointRecord) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::POINT; + std::vector points = {4.0, 2.0}; + float rotation(0.0); + auto csys = frame->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + + // Get Region as casacore::Record + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + auto image_shape = frame->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "LCBox"); + ASSERT_FALSE(region_record.asBool("oneRel")); + auto blc = region_record.asArrayFloat("blc").tovector(); + auto trc = region_record.asArrayFloat("trc").tovector(); + ASSERT_EQ(blc.size(), 3); + ASSERT_EQ(trc.size(), 3); + ASSERT_FLOAT_EQ(blc[0], points[0]); + ASSERT_FLOAT_EQ(blc[1], points[1]); + ASSERT_FLOAT_EQ(blc[2], 0.0); // min channel + ASSERT_FLOAT_EQ(trc[0], points[0]); + ASSERT_FLOAT_EQ(trc[1], points[1]); + ASSERT_FLOAT_EQ(trc[2], 9.0); // max channel +} + +TEST_F(RegionTest, TestReferenceImageLineRecord) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::LINE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(0.0); + auto csys = frame->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + + // Get Region as casacore::Record (from control points, no casacore Line region) + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + auto image_shape = frame->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "line"); + ASSERT_FALSE(region_record.asBool("oneRel")); + auto x = region_record.asArrayFloat("x").tovector(); + auto y = region_record.asArrayFloat("y").tovector(); + ASSERT_EQ(x.size(), 2); + ASSERT_EQ(y.size(), 2); + ASSERT_FLOAT_EQ(x[0], points[0]); + ASSERT_FLOAT_EQ(x[1], points[2]); + ASSERT_FLOAT_EQ(y[0], points[1]); + ASSERT_FLOAT_EQ(y[1], points[3]); +} + +TEST_F(RegionTest, TestReferenceImageRectangleRecord) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::RECTANGLE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(0.0); + auto csys = frame->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + + // Get Region as casacore::Record + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + auto image_shape = frame->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "LCPolygon"); // box corners set as polygon + ASSERT_FALSE(region_record.asBool("oneRel")); + // x, y order is [blc, brc, trc, tlc, blc] + auto x = region_record.asArrayDouble("x").tovector(); + auto y = region_record.asArrayDouble("y").tovector(); + float left_x = points[0] - (points[2] / 2.0); + float right_x = points[0] + (points[2] / 2.0); + float bottom_y = points[1] - (points[3] / 2.0); + float top_y = points[1] + (points[3] / 2.0); + ASSERT_EQ(x.size(), 5); // casacore repeats first point to close polygon + ASSERT_EQ(y.size(), 5); // casacore repeats first point to close polygon + ASSERT_FLOAT_EQ(x[0], left_x); + ASSERT_FLOAT_EQ(x[1], right_x); + ASSERT_FLOAT_EQ(x[2], right_x); + ASSERT_FLOAT_EQ(x[3], left_x); + ASSERT_FLOAT_EQ(x[4], left_x); + ASSERT_FLOAT_EQ(y[0], bottom_y); + ASSERT_FLOAT_EQ(y[1], bottom_y); + ASSERT_FLOAT_EQ(y[2], top_y); + ASSERT_FLOAT_EQ(y[3], top_y); + ASSERT_FLOAT_EQ(y[4], bottom_y); +} + +TEST_F(RegionTest, TestReferenceImageRotboxRecord) { + // Record is for unrotated rectangle; RegionState used for angle in export + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); // 10x10x10 + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::RECTANGLE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(30.0); + auto csys = frame->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + + // Get Region as casacore::Record + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + auto image_shape = frame->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "LCPolygon"); // box corners set as polygon + ASSERT_FALSE(region_record.asBool("oneRel")); + // x, y order is [blc, brc, trc, tlc, blc] + auto x = region_record.asArrayDouble("x").tovector(); + auto y = region_record.asArrayDouble("y").tovector(); + float left_x = points[0] - (points[2] / 2.0); + float right_x = points[0] + (points[2] / 2.0); + float bottom_y = points[1] - (points[3] / 2.0); + float top_y = points[1] + (points[3] / 2.0); + ASSERT_EQ(x.size(), 5); // repeats first point to close polygon + ASSERT_EQ(y.size(), 5); // repeats first point to close polygon + ASSERT_FLOAT_EQ(x[0], left_x); + ASSERT_FLOAT_EQ(x[1], right_x); + ASSERT_FLOAT_EQ(x[2], right_x); + ASSERT_FLOAT_EQ(x[3], left_x); + ASSERT_FLOAT_EQ(x[4], left_x); + ASSERT_FLOAT_EQ(y[0], bottom_y); + ASSERT_FLOAT_EQ(y[1], bottom_y); + ASSERT_FLOAT_EQ(y[2], top_y); + ASSERT_FLOAT_EQ(y[3], top_y); + ASSERT_FLOAT_EQ(y[4], bottom_y); +} + +TEST_F(RegionTest, TestReferenceImageEllipseRecord) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::ELLIPSE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(0.0); + auto csys = frame->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + + // Get Region as casacore::Record + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + auto image_shape = frame->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "LCEllipsoid"); + ASSERT_FALSE(region_record.asBool("oneRel")); + auto center = region_record.asArrayFloat("center").tovector(); + ASSERT_FLOAT_EQ(center[0], points[0]); + ASSERT_FLOAT_EQ(center[1], points[1]); + auto radii = region_record.asArrayFloat("radii").tovector(); + ASSERT_FLOAT_EQ(radii[0], points[2]); + ASSERT_FLOAT_EQ(radii[1], points[3]); +} + +TEST_F(RegionTest, TestReferenceImagePolygonRecord) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::POLYGON; + std::vector points = {5.0, 5.0, 4.0, 3.0, 1.0, 6.0, 3.0, 8.0}; // 4 points + float rotation(0.0); + auto csys = frame->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + + // Get Region as casacore::Record + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + auto image_shape = frame->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "LCPolygon"); + ASSERT_FALSE(region_record.asBool("oneRel")); + // x, y points + auto x = region_record.asArrayFloat("x").tovector(); + auto y = region_record.asArrayFloat("y").tovector(); + ASSERT_EQ(x.size(), 5); // casacore repeats first point to close polygon + ASSERT_EQ(y.size(), 5); // casacore repeats first point to close polygon + ASSERT_FLOAT_EQ(x[0], points[0]); + ASSERT_FLOAT_EQ(x[1], points[2]); + ASSERT_FLOAT_EQ(x[2], points[4]); + ASSERT_FLOAT_EQ(x[3], points[6]); + ASSERT_FLOAT_EQ(y[0], points[1]); + ASSERT_FLOAT_EQ(y[1], points[3]); + ASSERT_FLOAT_EQ(y[2], points[5]); + ASSERT_FLOAT_EQ(y[3], points[7]); +} diff --git a/test/TestRegionHistogram.cc b/test/TestRegionHistogram.cc new file mode 100644 index 000000000..3ac57d5b1 --- /dev/null +++ b/test/TestRegionHistogram.cc @@ -0,0 +1,94 @@ +/* This file is part of the CARTA Image Viewer: https://github.com/CARTAvis/carta-backend + Copyright 2018- Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), + Associated Universities, Inc. (AUI) and the Inter-University Institute for Data Intensive Astronomy (IDIA) + SPDX-License-Identifier: GPL-3.0-or-later +*/ + +#include +#include + +#include "CommonTestUtilities.h" +#include "ImageData/FileLoader.h" +#include "Region/Region.h" +#include "Region/RegionHandler.h" +#include "src/Frame/Frame.h" + +using namespace carta; + +class RegionHistogramTest : public ::testing::Test { +public: + static bool SetRegion(carta::RegionHandler& region_handler, int file_id, int& region_id, const std::vector& points, + std::shared_ptr csys, bool is_annotation) { + std::vector control_points; + for (auto i = 0; i < points.size(); i += 2) { + control_points.push_back(Message::Point(points[i], points[i + 1])); + } + + // Define RegionState for line region and set region (region_id updated) + auto npoints(control_points.size()); + CARTA::RegionType region_type = CARTA::RegionType::POLYGON; + if (is_annotation) { + region_type = CARTA::RegionType::ANNPOLYGON; + } + RegionState region_state(file_id, region_type, control_points, 0.0); + return region_handler.SetRegion(region_id, region_state, csys); + } + + static bool RegionHistogram(const std::string& image_path, const std::vector& endpoints, + CARTA::RegionHistogramData& region_histogram, bool is_annotation = false) { + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + carta::RegionHandler region_handler; + + // Set polygon region + int file_id(0), region_id(-1); + auto csys = frame->CoordinateSystem(); + if (!SetRegion(region_handler, file_id, region_id, endpoints, csys, is_annotation)) { + return false; + } + + // Set histogram requirements + auto histogram_req_message = Message::SetHistogramRequirements(file_id, region_id); + std::vector histogram_configs = { + histogram_req_message.histograms().begin(), histogram_req_message.histograms().end()}; + if (!region_handler.SetHistogramRequirements(region_id, file_id, frame, histogram_configs)) { + return false; + } + + // Get histogram + return region_handler.FillRegionHistogramData( + [&](CARTA::RegionHistogramData histogram_data) { region_histogram = histogram_data; }, region_id, file_id); + } +}; + +TEST_F(RegionHistogramTest, TestFitsRegionHistogram) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {1.0, 1.0, 1.0, 4.0, 4.0, 4.0, 4.0, 1.0}; + CARTA::RegionHistogramData histogram_data; + bool ok = RegionHistogram(image_path, endpoints, histogram_data); + + // Check histogram data fields + ASSERT_TRUE(ok); + ASSERT_EQ(histogram_data.file_id(), 0); + ASSERT_EQ(histogram_data.region_id(), 1); + ASSERT_EQ(histogram_data.channel(), 0); + ASSERT_EQ(histogram_data.stokes(), 0); + ASSERT_TRUE(histogram_data.has_histograms()); + ASSERT_EQ(histogram_data.progress(), 1.0); + ASSERT_TRUE(histogram_data.has_config()); + int expected_num_bins = sqrt(4 * 4); // region bounding box is 4x4 + ASSERT_EQ(histogram_data.histograms().num_bins(), expected_num_bins); + + FitsDataReader reader(image_path); + auto image_data = reader.ReadRegion({1, 1, 0}, {5, 5, 1}); + double expected_mean = std::accumulate(image_data.begin(), image_data.end(), 0.0) / image_data.size(); + ASSERT_DOUBLE_EQ(histogram_data.histograms().mean(), expected_mean); +} + +TEST_F(RegionHistogramTest, TestFitsAnnotationRegionHistogram) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {0.0, 0.0, 0.0, 3.0, 3.0, 3.0, 3.0, 0.0}; + CARTA::RegionHistogramData histogram_data; + bool ok = RegionHistogram(image_path, endpoints, histogram_data, true); + ASSERT_FALSE(ok); +} diff --git a/test/TestRegionImportExport.cc b/test/TestRegionImportExport.cc new file mode 100644 index 000000000..31bae0f9b --- /dev/null +++ b/test/TestRegionImportExport.cc @@ -0,0 +1,362 @@ +/* This file is part of the CARTA Image Viewer: https://github.com/CARTAvis/carta-backend + Copyright 2018- Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), + Associated Universities, Inc. (AUI) and the Inter-University Institute for Data Intensive Astronomy (IDIA) + SPDX-License-Identifier: GPL-3.0-or-later +*/ + +#include +#include + +#include + +#include "CommonTestUtilities.h" +#include "ImageData/FileLoader.h" +#include "Region/Region.h" +#include "Region/RegionHandler.h" +#include "src/Frame/Frame.h" + +using namespace carta; + +class RegionImportExportTest : public ::testing::Test { +public: + static bool SetRegion(carta::RegionHandler& region_handler, int file_id, int& region_id, CARTA::RegionType type, + const std::vector& points, float rotation, std::shared_ptr csys) { + std::vector control_points; + for (auto i = 0; i < points.size(); i += 2) { + control_points.push_back(Message::Point(points[i], points[i + 1])); + } + + // Define RegionState and set region (region_id updated) + RegionState region_state(file_id, type, control_points, rotation); + return region_handler.SetRegion(region_id, region_state, csys); + } + + static int SetAllRegions(carta::RegionHandler& region_handler, int file_id, std::shared_ptr csys) { + std::vector point_points = {5.0, 5.0}; + std::vector line_rectangle_ellipse_points = {5.0, 5.0, 4.0, 3.0}; + std::vector circle_points = {5.0, 5.0, 3.0, 3.0}; + std::vector poly_points = {5.0, 5.0, 4.0, 3.0, 1.0, 6.0, 3.0, 8.0}; + std::unordered_map> region_points = {{CARTA::POINT, point_points}, + {CARTA::LINE, line_rectangle_ellipse_points}, {CARTA::POLYLINE, poly_points}, {CARTA::RECTANGLE, line_rectangle_ellipse_points}, + {CARTA::ELLIPSE, line_rectangle_ellipse_points}, {CARTA::POLYGON, poly_points}, {CARTA::ANNPOINT, point_points}, + {CARTA::ANNLINE, line_rectangle_ellipse_points}, {CARTA::ANNPOLYLINE, poly_points}, + {CARTA::ANNRECTANGLE, line_rectangle_ellipse_points}, {CARTA::ANNELLIPSE, line_rectangle_ellipse_points}, + {CARTA::ANNPOLYGON, poly_points}, {CARTA::ANNVECTOR, line_rectangle_ellipse_points}, + {CARTA::ANNRULER, line_rectangle_ellipse_points}, {CARTA::ANNTEXT, line_rectangle_ellipse_points}, + {CARTA::ANNCOMPASS, circle_points}}; + float rotation(0.0); + int num_regions(0); + + // Add all region types + int region_id(-1); + for (int i = 0; i < CARTA::RegionType_ARRAYSIZE; ++i) { + CARTA::RegionType type = static_cast(i); + if (type == CARTA::RegionType::ANNULUS) { + continue; + } + auto points = region_points[type]; + if (!SetRegion(region_handler, file_id, region_id, type, points, rotation, csys)) { + return 0; + } + num_regions++; + region_id = -1; + } + + // Add special region types: rotbox, circle (analytical and annotation) + rotation = 30.0; + if (!SetRegion(region_handler, file_id, region_id, CARTA::RECTANGLE, line_rectangle_ellipse_points, rotation, csys)) { + return 0; + } + num_regions++; + + region_id = -1; + if (!SetRegion(region_handler, file_id, region_id, CARTA::ANNRECTANGLE, line_rectangle_ellipse_points, rotation, csys)) { + return 0; + } + num_regions++; + + region_id = -1; + rotation = 0.0; + if (!SetRegion(region_handler, file_id, region_id, CARTA::ELLIPSE, circle_points, 0.0, csys)) { + return 0; + } + num_regions++; + + region_id = -1; + if (!SetRegion(region_handler, file_id, region_id, CARTA::ANNELLIPSE, circle_points, 0.0, csys)) { + return 0; + } + num_regions++; + + return num_regions; + } + + static CARTA::RegionStyle GetRegionStyle(CARTA::RegionType type) { + bool is_annotation = type > CARTA::POLYGON; + std::string color = is_annotation ? "#FFBA01" : "#2EE6D6"; + + CARTA::RegionStyle region_style; + region_style.set_color(color); + region_style.set_line_width(2); + + if (is_annotation) { + // Default fields set by frontend + auto annotation_style = region_style.mutable_annotation_style(); + if (type == CARTA::ANNPOINT) { + annotation_style->set_point_shape(CARTA::PointAnnotationShape::SQUARE); + annotation_style->set_point_width(6); + } else if (type == CARTA::ANNTEXT || type == CARTA::ANNCOMPASS || type == CARTA::ANNRULER) { + annotation_style->set_font("Helvetica"); + annotation_style->set_font_size(20); + annotation_style->set_font_style("Normal"); + if (type == CARTA::ANNTEXT) { + annotation_style->set_text_label0("Text"); + annotation_style->set_text_position(CARTA::TextAnnotationPosition::CENTER); + } else if (type == CARTA::ANNCOMPASS) { + annotation_style->set_coordinate_system("PIXEL"); + annotation_style->set_is_east_arrow(true); + annotation_style->set_is_north_arrow(true); + annotation_style->set_text_label0("N"); + annotation_style->set_text_label1("E"); + } else if (type == CARTA::ANNRULER) { + annotation_style->set_coordinate_system("PIXEL"); + } + } + } + return region_style; + } + + static std::string ConcatContents(const std::vector& string_vector) { + std::string one_string; + for (auto& item : string_vector) { + one_string.append(item + "\n"); + } + return one_string; + } +}; + +TEST_F(RegionImportExportTest, TestCrtfPixExportImport) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set all region types in frame0 + carta::RegionHandler region_handler; + int file_id(0); + int num_regions = SetAllRegions(region_handler, file_id, frame0->CoordinateSystem()); + // All CARTA regions except ANNULUS (- 1) plus 2 rotbox and 2 circle (+ 4) + ASSERT_EQ(num_regions, CARTA::RegionType_ARRAYSIZE - 1 + 4); + + // Set RegionStyle map for export + std::map region_style_map; + for (int i = 0; i < num_regions; ++i) { + int region_id = i + 1; // region 0 is cursor + CARTA::RegionType region_type = region_handler.GetRegion(region_id)->GetRegionState().type; + CARTA::RegionStyle style = GetRegionStyle(region_type); + region_style_map[region_id] = style; + } + std::string filename; // do not export to file + + // Export all regions in frame0 (reference image) + CARTA::ExportRegionAck export_ack0; + region_handler.ExportRegion(file_id, frame0, CARTA::CRTF, CARTA::PIXEL, region_style_map, filename, export_ack0); + // Check that all regions were exported + ASSERT_EQ(export_ack0.contents_size(), num_regions + 2); // header, textbox for text + + // Import all regions in frame0 (reference image) + std::vector export_contents = {export_ack0.contents().begin(), export_ack0.contents().end()}; + auto contents_string = ConcatContents(export_contents); + bool file_is_filename(false); + CARTA::ImportRegionAck import_ack0; + region_handler.ImportRegion(file_id, frame0, CARTA::CRTF, contents_string, file_is_filename, import_ack0); + // Check that all regions were imported + ASSERT_EQ(import_ack0.regions_size(), num_regions); + + // Export all regions in frame1 (matched image) + file_id = 1; + CARTA::ExportRegionAck export_ack1; + region_handler.ExportRegion(file_id, frame1, CARTA::CRTF, CARTA::PIXEL, region_style_map, filename, export_ack1); + // Check that all regions were exported + ASSERT_EQ(export_ack1.contents_size(), num_regions + 2); // header, textbox for text + + // Import all regions in frame1 (matched image) + export_contents = {export_ack1.contents().begin(), export_ack1.contents().end()}; + contents_string = ConcatContents(export_contents); + CARTA::ImportRegionAck import_ack1; + region_handler.ImportRegion(file_id, frame1, CARTA::CRTF, contents_string, file_is_filename, import_ack1); + // Check that all regions were imported + ASSERT_EQ(import_ack1.regions_size(), num_regions); +} + +TEST_F(RegionImportExportTest, TestCrtfWorldExportImport) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set all region types in frame0 + carta::RegionHandler region_handler; + int file_id(0); + int num_regions = SetAllRegions(region_handler, file_id, frame0->CoordinateSystem()); + // All CARTA regions except ANNULUS (- 1) plus 2 rotbox and 2 circle (+ 4) + ASSERT_EQ(num_regions, CARTA::RegionType_ARRAYSIZE - 1 + 4); + + // Export all regions in frame0 (reference image) + std::map region_style_map; + for (int i = 0; i < num_regions; ++i) { + int region_id = i + 1; // region 0 is cursor + CARTA::RegionType region_type = region_handler.GetRegion(region_id)->GetRegionState().type; + CARTA::RegionStyle style = GetRegionStyle(region_type); + region_style_map[region_id] = style; + } + std::string filename; // do not export to file + CARTA::ExportRegionAck export_ack0; + region_handler.ExportRegion(file_id, frame0, CARTA::CRTF, CARTA::WORLD, region_style_map, filename, export_ack0); + // Check that all regions were exported + ASSERT_EQ(export_ack0.contents_size(), num_regions + 2); // header, textbox for text + + // Import all regions in frame0 (reference image) + std::vector export_contents = {export_ack0.contents().begin(), export_ack0.contents().end()}; + auto contents_string = ConcatContents(export_contents); + bool file_is_filename(false); + CARTA::ImportRegionAck import_ack0; + region_handler.ImportRegion(file_id, frame0, CARTA::CRTF, contents_string, file_is_filename, import_ack0); + // Check that all regions were imported + ASSERT_EQ(import_ack0.regions_size(), num_regions); + + // Export all regions in frame1 (matched image) + file_id = 1; + CARTA::ExportRegionAck export_ack1; + region_handler.ExportRegion(file_id, frame1, CARTA::CRTF, CARTA::WORLD, region_style_map, filename, export_ack1); + // Check that all regions were exported + ASSERT_EQ(export_ack1.contents_size(), num_regions + 2); // header, textbox for text + + // Import all regions in frame1 (matched image) + export_contents = {export_ack1.contents().begin(), export_ack1.contents().end()}; + contents_string = ConcatContents(export_contents); + CARTA::ImportRegionAck import_ack1; + region_handler.ImportRegion(file_id, frame1, CARTA::CRTF, contents_string, file_is_filename, import_ack1); + // Check that all regions were imported + ASSERT_EQ(import_ack1.regions_size(), num_regions); +} + +TEST_F(RegionImportExportTest, TestDs9PixExportImport) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set all region types in frame0 + carta::RegionHandler region_handler; + int file_id(0); + int num_regions = SetAllRegions(region_handler, file_id, frame0->CoordinateSystem()); + // All CARTA regions except ANNULUS (- 1) plus 2 rotbox and 2 circle (+ 4) + ASSERT_EQ(num_regions, CARTA::RegionType_ARRAYSIZE - 1 + 4); + + // Export all regions in frame0 (reference image) + std::map region_style_map; + for (int i = 0; i < num_regions; ++i) { + int region_id = i + 1; // region 0 is cursor + CARTA::RegionType region_type = region_handler.GetRegion(region_id)->GetRegionState().type; + CARTA::RegionStyle style = GetRegionStyle(region_type); + region_style_map[region_id] = style; + } + std::string filename; // do not export to file + CARTA::ExportRegionAck export_ack0; + region_handler.ExportRegion(file_id, frame0, CARTA::DS9_REG, CARTA::PIXEL, region_style_map, filename, export_ack0); + // Check that all regions were exported + ASSERT_EQ(export_ack0.contents_size(), num_regions + 3); // header + globals, coord sys, textbox for text + + // Import all regions in frame0 (reference image) + std::vector export_contents = {export_ack0.contents().begin(), export_ack0.contents().end()}; + auto contents_string = ConcatContents(export_contents); + bool file_is_filename(false); + CARTA::ImportRegionAck import_ack0; + region_handler.ImportRegion(file_id, frame0, CARTA::DS9_REG, contents_string, file_is_filename, import_ack0); + // Check that all regions were imported + ASSERT_EQ(import_ack0.regions_size(), num_regions); + + // Export all regions in frame1 (matched image) + file_id = 1; + CARTA::ExportRegionAck export_ack1; + region_handler.ExportRegion(file_id, frame1, CARTA::DS9_REG, CARTA::PIXEL, region_style_map, filename, export_ack1); + // Check that all regions were exported + ASSERT_EQ(export_ack1.contents_size(), num_regions + 2); // header + globals, coord sys (textbox + text in same string) + + // Import all regions in frame1 (matched image) + export_contents = {export_ack1.contents().begin(), export_ack1.contents().end()}; + contents_string = ConcatContents(export_contents); + CARTA::ImportRegionAck import_ack1; + region_handler.ImportRegion(file_id, frame1, CARTA::DS9_REG, contents_string, file_is_filename, import_ack1); + // Check that all regions were imported + ASSERT_EQ(import_ack1.regions_size(), num_regions); +} + +TEST_F(RegionImportExportTest, TestDs9WorldExportImport) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set all region types in frame0 + carta::RegionHandler region_handler; + int file_id(0); + int num_regions = SetAllRegions(region_handler, file_id, frame0->CoordinateSystem()); + // All CARTA regions except ANNULUS (- 1) plus 2 rotbox and 2 circle (+ 4) + ASSERT_EQ(num_regions, CARTA::RegionType_ARRAYSIZE - 1 + 4); + + // Export all regions in frame0 (reference image) + std::map region_style_map; + for (int i = 0; i < num_regions; ++i) { + int region_id = i + 1; // region 0 is cursor + CARTA::RegionType region_type = region_handler.GetRegion(region_id)->GetRegionState().type; + CARTA::RegionStyle style = GetRegionStyle(region_type); + region_style_map[region_id] = style; + } + std::string filename; // do not export to file + CARTA::ExportRegionAck export_ack0; + region_handler.ExportRegion(file_id, frame0, CARTA::DS9_REG, CARTA::WORLD, region_style_map, filename, export_ack0); + // Check that all regions were exported + ASSERT_EQ(export_ack0.contents_size(), num_regions + 2); // header + globals, coord sys (textbox + text in same string) + + // Import all regions in frame0 (reference image) + std::vector export_contents = {export_ack0.contents().begin(), export_ack0.contents().end()}; + auto contents_string = ConcatContents(export_contents); + bool file_is_filename(false); + CARTA::ImportRegionAck import_ack0; + region_handler.ImportRegion(file_id, frame0, CARTA::DS9_REG, contents_string, file_is_filename, import_ack0); + // Check that all regions were imported + ASSERT_EQ(import_ack0.regions_size(), num_regions); + + // Export all regions in frame1 (matched image) + file_id = 1; + CARTA::ExportRegionAck export_ack1; + region_handler.ExportRegion(file_id, frame1, CARTA::DS9_REG, CARTA::WORLD, region_style_map, filename, export_ack1); + // Check that all regions were exported + ASSERT_EQ(export_ack1.contents_size(), num_regions + 2); // header + globals, coord sys (textbox + text in same string) + + // Import all regions in frame1 (matched image) + export_contents = {export_ack1.contents().begin(), export_ack1.contents().end()}; + contents_string = ConcatContents(export_contents); + CARTA::ImportRegionAck import_ack1; + region_handler.ImportRegion(file_id, frame1, CARTA::DS9_REG, contents_string, file_is_filename, import_ack1); + // Check that all regions were imported + ASSERT_EQ(import_ack1.regions_size(), num_regions); +} diff --git a/test/TestRegionMatched.cc b/test/TestRegionMatched.cc new file mode 100644 index 000000000..5da734e3e --- /dev/null +++ b/test/TestRegionMatched.cc @@ -0,0 +1,447 @@ +/* This file is part of the CARTA Image Viewer: https://github.com/CARTAvis/carta-backend + Copyright 2018- Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), + Associated Universities, Inc. (AUI) and the Inter-University Institute for Data Intensive Astronomy (IDIA) + SPDX-License-Identifier: GPL-3.0-or-later +*/ + +#include +#include + +#include "CommonTestUtilities.h" +#include "ImageData/FileLoader.h" +#include "Region/Region.h" +#include "Region/RegionHandler.h" +#include "src/Frame/Frame.h" + +using namespace carta; + +class RegionMatchedTest : public ::testing::Test { +public: + static bool SetRegion(carta::RegionHandler& region_handler, int file_id, int& region_id, CARTA::RegionType type, + const std::vector& points, float rotation, std::shared_ptr csys) { + std::vector control_points; + for (auto i = 0; i < points.size(); i += 2) { + control_points.push_back(Message::Point(points[i], points[i + 1])); + } + + // Define RegionState and set region (region_id updated) + RegionState region_state(file_id, type, control_points, rotation); + return region_handler.SetRegion(region_id, region_state, csys); + } +}; + +TEST_F(RegionMatchedTest, TestMatchedImageRectangleLCRegion) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set rectangle in frame 0 + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::RECTANGLE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(0.0); + auto csys = frame0->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + + // Get Region as 2D LCRegion in frame1 + file_id = 1; + csys = frame1->CoordinateSystem(); + auto image_shape = frame1->ImageShape(); + auto lc_region = region->GetImageRegion(file_id, csys, image_shape); + + // Check LCRegion + ASSERT_TRUE(lc_region); // shared_ptr + ASSERT_EQ(lc_region->ndim(), image_shape.size()); + ASSERT_EQ(lc_region->latticeShape(), image_shape); + ASSERT_EQ(lc_region->shape(), casacore::IPosition(2, 5, 3)); +} + +TEST_F(RegionMatchedTest, TestMatchedImageRotboxLCRegion) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set rectangle in frame 0 + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::RECTANGLE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(30.0); + auto csys = frame0->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + + // Get Region as 2D LCRegion in frame1 + file_id = 1; + csys = frame1->CoordinateSystem(); + auto image_shape = frame1->ImageShape(); + auto lc_region = region->GetImageRegion(file_id, csys, image_shape); + + // Check LCRegion + ASSERT_TRUE(lc_region); // shared_ptr + ASSERT_EQ(lc_region->ndim(), image_shape.size()); + ASSERT_EQ(lc_region->latticeShape(), image_shape); + ASSERT_EQ(lc_region->shape(), casacore::IPosition(2, 5, 5)); +} + +TEST_F(RegionMatchedTest, TestMatchedImageEllipseLCRegion) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set region in frame0 + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::ELLIPSE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(0.0); + auto csys = frame0->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + + // Get Region as 2D LCRegion in frame1 + file_id = 1; + csys = frame1->CoordinateSystem(); + auto image_shape = frame1->ImageShape(); + auto lc_region = region->GetImageRegion(file_id, csys, image_shape); + + // Check LCRegion + ASSERT_TRUE(lc_region); // shared_ptr + ASSERT_EQ(lc_region->ndim(), image_shape.size()); + ASSERT_EQ(lc_region->latticeShape(), image_shape); + ASSERT_EQ(lc_region->shape(), casacore::IPosition(2, 7, 9)); +} + +TEST_F(RegionMatchedTest, TestMatchedImagePolygonLCRegion) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set region in frame0 + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::POLYGON; + std::vector points = {5.0, 5.0, 4.0, 3.0, 1.0, 6.0, 3.0, 8.0}; + float rotation(0.0); + auto csys = frame0->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + + // Get Region as 2D LCRegion in frame1 + file_id = 1; + csys = frame1->CoordinateSystem(); + auto image_shape = frame1->ImageShape(); + auto lc_region = region->GetImageRegion(file_id, csys, image_shape); + + // Check LCRegion + ASSERT_TRUE(lc_region); // shared_ptr + ASSERT_EQ(lc_region->ndim(), image_shape.size()); + ASSERT_EQ(lc_region->latticeShape(), image_shape); + ASSERT_EQ(lc_region->shape(), casacore::IPosition(2, 5, 6)); +} + +TEST_F(RegionMatchedTest, TestMatchedImagePointRecord) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set region in frame0 + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::POINT; + std::vector points = {4.0, 2.0}; + float rotation(0.0); + auto csys = frame0->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + + // Get Region as casacore::Record in frame1 + file_id = 1; + csys = frame1->CoordinateSystem(); + auto image_shape = frame1->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + + // Check record + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "LCBox"); // box with blc = trc + ASSERT_TRUE(region_record.asBool("oneRel")); // 1-based pixels + auto blc = region_record.asArrayFloat("blc").tovector(); + auto trc = region_record.asArrayFloat("trc").tovector(); + ASSERT_EQ(blc.size(), 2); + ASSERT_EQ(trc.size(), 2); + ASSERT_FLOAT_EQ(blc[0], points[0] + 1.0); + ASSERT_FLOAT_EQ(blc[1], points[1] + 1.0); + ASSERT_FLOAT_EQ(trc[0], points[0] + 1.0); + ASSERT_FLOAT_EQ(trc[1], points[1] + 1.0); +} + +TEST_F(RegionMatchedTest, TestMatchedImageLineRecord) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set region in frame0 + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::LINE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(0.0); + auto csys = frame0->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + + // Get Region as casacore::Record in frame1 + file_id = 1; + csys = frame1->CoordinateSystem(); + auto image_shape = frame1->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + + // Check record + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "line"); + ASSERT_FALSE(region_record.asBool("oneRel")); + auto x = region_record.asArrayDouble("x").tovector(); + auto y = region_record.asArrayDouble("y").tovector(); + ASSERT_EQ(x.size(), 2); + ASSERT_EQ(y.size(), 2); + ASSERT_FLOAT_EQ(x[0], points[0]); + ASSERT_FLOAT_EQ(x[1], points[2]); + ASSERT_FLOAT_EQ(y[0], points[1]); + ASSERT_FLOAT_EQ(y[1], points[3]); +} + +TEST_F(RegionMatchedTest, TestMatchedImageRectangleRecord) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set region in frame0 + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::RECTANGLE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(0.0); + auto csys = frame0->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + + // Get Region as casacore::Record in frame1 + file_id = 1; + csys = frame1->CoordinateSystem(); + auto image_shape = frame1->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + + // Check record + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "LCPolygon"); // box corners set as polygon + ASSERT_TRUE(region_record.asBool("oneRel")); // 1-based pixels + // x, y order is [blc, brc, trc, tlc, blc] + auto x = region_record.asArrayFloat("x").tovector(); + auto y = region_record.asArrayFloat("y").tovector(); + float left_x = points[0] - (points[2] / 2.0) + 1.0; + float right_x = points[0] + (points[2] / 2.0) + 1.0; + float bottom_y = points[1] - (points[3] / 2.0) + 1.0; + float top_y = points[1] + (points[3] / 2.0) + 1.0; + ASSERT_EQ(x.size(), 5); // casacore repeats first point to close polygon + ASSERT_EQ(y.size(), 5); // casacore repeats first point to close polygon + ASSERT_FLOAT_EQ(x[0], left_x); + ASSERT_FLOAT_EQ(x[1], right_x); + ASSERT_FLOAT_EQ(y[0], bottom_y); + ASSERT_FLOAT_EQ(y[2], top_y); +} + +TEST_F(RegionMatchedTest, TestMatchedImageRotboxRecord) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set region in frame0 + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::RECTANGLE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(30.0); + auto csys = frame0->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + + // Get Region as casacore::Record in frame1 + file_id = 1; + csys = frame1->CoordinateSystem(); + auto image_shape = frame1->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + + // Check record + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "LCPolygon"); // box corners set as polygon + ASSERT_FALSE(region_record.asBool("oneRel")); + // x, y order is [blc, brc, trc, tlc, blc] + auto x = region_record.asArrayFloat("x").tovector(); + auto y = region_record.asArrayFloat("y").tovector(); + // keep original rectangle pixel points with rotation for export + float left_x = points[0] - (points[2] / 2.0); + float right_x = points[0] + (points[2] / 2.0); + float bottom_y = points[1] - (points[3] / 2.0); + float top_y = points[1] + (points[3] / 2.0); + ASSERT_EQ(x.size(), 4); + ASSERT_EQ(y.size(), 4); + ASSERT_FLOAT_EQ(x[0], left_x); + ASSERT_FLOAT_EQ(x[1], right_x); + ASSERT_FLOAT_EQ(x[2], right_x); + ASSERT_FLOAT_EQ(x[3], left_x); + ASSERT_FLOAT_EQ(y[0], bottom_y); + ASSERT_FLOAT_EQ(y[1], bottom_y); + ASSERT_FLOAT_EQ(y[2], top_y); + ASSERT_FLOAT_EQ(y[3], top_y); +} + +TEST_F(RegionMatchedTest, TestMatchedImageEllipseRecord) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set region in frame0 + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::ELLIPSE; + std::vector points = {5.0, 5.0, 4.0, 3.0}; + float rotation(0.0); + auto csys = frame0->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + + // Get Region as casacore::Record in frame1 + file_id = 1; + csys = frame1->CoordinateSystem(); + auto image_shape = frame1->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + + // Check record + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "LCEllipsoid"); + ASSERT_TRUE(region_record.asBool("oneRel")); // 1-based pixels + auto center = region_record.asArrayFloat("center").tovector(); + ASSERT_FLOAT_EQ(center[0], points[0] + 1.0); + ASSERT_FLOAT_EQ(center[1], points[1] + 1.0); + auto radii = region_record.asArrayFloat("radii").tovector(); + ASSERT_FLOAT_EQ(radii[0], points[3]); + ASSERT_FLOAT_EQ(radii[1], points[2]); +} + +TEST_F(RegionMatchedTest, TestMatchedImagePolygonRecord) { + // frame 0 + std::string image_path0 = FileFinder::FitsImagePath("noise_10px_10px.fits"); + std::shared_ptr loader0(carta::FileLoader::GetLoader(image_path0)); + std::shared_ptr frame0(new Frame(0, loader0, "0")); + // frame 1 + std::string image_path1 = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::shared_ptr loader1(carta::FileLoader::GetLoader(image_path1)); + std::shared_ptr frame1(new Frame(0, loader1, "0")); + + // Set region in frame0 + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + CARTA::RegionType region_type = CARTA::RegionType::POLYGON; + std::vector points = {5.0, 5.0, 4.0, 3.0, 1.0, 6.0, 3.0, 8.0}; // 4 points + float rotation(0.0); + auto csys = frame0->CoordinateSystem(); + bool ok = SetRegion(region_handler, file_id, region_id, region_type, points, rotation, csys); + ASSERT_TRUE(ok); + auto region = region_handler.GetRegion(region_id); + ASSERT_TRUE(region); // shared_ptr + + // Get Region as casacore::Record in frame1 + file_id = 1; + csys = frame1->CoordinateSystem(); + auto image_shape = frame1->ImageShape(); + auto region_record = region->GetImageRegionRecord(file_id, csys, image_shape); + + // Check record + ASSERT_GT(region_record.nfields(), 0); + ASSERT_EQ(region_record.asInt("isRegion"), 1); + ASSERT_EQ(region_record.asString("name"), "LCPolygon"); + ASSERT_TRUE(region_record.asBool("oneRel")); // 1-based pixels + // x, y points + auto x = region_record.asArrayFloat("x").tovector(); + auto y = region_record.asArrayFloat("y").tovector(); + ASSERT_EQ(x.size(), 5); // casacore repeats first point to close polygon + ASSERT_EQ(y.size(), 5); // casacore repeats first point to close polygon + ASSERT_FLOAT_EQ(x[0], points[0] + 1.0); + ASSERT_FLOAT_EQ(x[1], points[2] + 1.0); + ASSERT_FLOAT_EQ(x[2], points[4] + 1.0); + ASSERT_FLOAT_EQ(x[3], points[6] + 1.0); + ASSERT_FLOAT_EQ(y[0], points[1] + 1.0); + ASSERT_FLOAT_EQ(y[1], points[3] + 1.0); + ASSERT_FLOAT_EQ(y[2], points[5] + 1.0); + ASSERT_FLOAT_EQ(y[3], points[7] + 1.0); +} diff --git a/test/TestRegionSpatialProfiles.cc b/test/TestRegionSpatialProfiles.cc new file mode 100644 index 000000000..849fc0cab --- /dev/null +++ b/test/TestRegionSpatialProfiles.cc @@ -0,0 +1,346 @@ +/* This file is part of the CARTA Image Viewer: https://github.com/CARTAvis/carta-backend + Copyright 2018- Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), + Associated Universities, Inc. (AUI) and the Inter-University Institute for Data Intensive Astronomy (IDIA) + SPDX-License-Identifier: GPL-3.0-or-later +*/ + +#include +#include + +#include "CommonTestUtilities.h" +#include "ImageData/FileLoader.h" +#include "Region/Region.h" +#include "Region/RegionHandler.h" +#include "src/Frame/Frame.h" + +using namespace carta; + +class RegionSpatialProfileTest : public ::testing::Test { +public: + static bool SetRegion(carta::RegionHandler& region_handler, int file_id, int& region_id, const std::vector& points, + std::shared_ptr csys, bool is_annotation) { + std::vector control_points; + for (auto i = 0; i < points.size(); i += 2) { + control_points.push_back(Message::Point(points[i], points[i + 1])); + } + + // Define RegionState for line region and set region (region_id updated) + auto npoints(control_points.size()); + CARTA::RegionType region_type; + if (npoints == 1) { + if (is_annotation) { + region_type = CARTA::RegionType::ANNPOINT; + } else { + region_type = CARTA::RegionType::POINT; + } + } else { + region_type = CARTA::RegionType::LINE; + if (is_annotation) { + if (npoints > 2) { + region_type = CARTA::RegionType::ANNPOLYLINE; + } else { + region_type = CARTA::RegionType::ANNLINE; + } + } else { + if (npoints > 2) { + region_type = CARTA::RegionType::POLYLINE; + } else { + region_type = CARTA::RegionType::LINE; + } + } + } + RegionState region_state(file_id, region_type, control_points, 0.0); + return region_handler.SetRegion(region_id, region_state, csys); + } + + static bool RegionSpatialProfile(const std::string& image_path, const std::vector& endpoints, + const std::vector& spatial_reqs, CARTA::SpatialProfileData& spatial_profile, + bool is_annotation = false) { + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + carta::RegionHandler region_handler; + + // Set line region + int file_id(0), region_id(-1); + auto csys = frame->CoordinateSystem(); + if (!SetRegion(region_handler, file_id, region_id, endpoints, csys, is_annotation)) { + return false; + } + + // Set spatial requirements + if (!region_handler.SetSpatialRequirements(region_id, file_id, frame, spatial_reqs)) { + return false; + } + + // Get spatial profiles + if (endpoints.size() == 2) { + std::vector profiles; + bool ok = region_handler.FillPointSpatialProfileData(file_id, region_id, profiles); + if (ok) { + spatial_profile = profiles[0]; + } + return ok; + } else { + return region_handler.FillLineSpatialProfileData( + file_id, region_id, [&](CARTA::SpatialProfileData profile_data) { spatial_profile = profile_data; }); + } + } + + static std::vector ProfileValues(CARTA::SpatialProfile& profile) { + std::string buffer = profile.raw_values_fp32(); + std::vector values(buffer.size() / sizeof(float)); + memcpy(values.data(), buffer.data(), buffer.size()); + return values; + } + + void SetUp() { + setenv("HDF5_USE_FILE_LOCKING", "FALSE", 0); + } + + static void TestAveragingWidthRange(int width, bool expected_width_range) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {0.0, 0.0, 9.0, 9.0}; + int start(0), end(0), mip(0); + std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; + CARTA::SpatialProfileData spatial_profile; + + if (expected_width_range) { + ASSERT_TRUE(RegionSpatialProfile(image_path, endpoints, spatial_reqs, spatial_profile)); + ASSERT_EQ(spatial_profile.profiles_size(), 1); + } else { + ASSERT_FALSE(RegionSpatialProfile(image_path, endpoints, spatial_reqs, spatial_profile)); + ASSERT_EQ(spatial_profile.profiles_size(), 0); + } + } +}; + +TEST_F(RegionSpatialProfileTest, TestSpatialRequirements) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + + // Set line region + carta::RegionHandler region_handler; + int file_id(0), region_id(-1); + std::vector endpoints = {0.0, 0.0, 9.0, 9.0}; + auto csys = frame->CoordinateSystem(); + bool is_annotation(false); + bool ok = SetRegion(region_handler, file_id, region_id, endpoints, csys, is_annotation); + ASSERT_TRUE(ok); + + // Set spatial requirements + int start(0), end(0), mip(0), width(3); + std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; + ok = region_handler.SetSpatialRequirements(region_id, file_id, frame, spatial_reqs); + ASSERT_TRUE(ok); + + // Get regions with spatial requirements for file id + auto region_ids = region_handler.GetSpatialReqRegionsForFile(file_id); + ASSERT_EQ(region_ids.size(), 1); + ASSERT_EQ(region_ids[0], region_id); + + // Get files with spatial requirements for region id + auto file_ids = region_handler.GetSpatialReqFilesForRegion(region_id); + ASSERT_EQ(file_ids.size(), 1); + ASSERT_EQ(file_ids[0], file_id); +} + +TEST_F(RegionSpatialProfileTest, FitsLineProfile) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {0.0, 0.0, 9.0, 9.0}; + int start(0), end(0), mip(0), width(3); + std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; + + CARTA::SpatialProfileData spatial_profile; + bool ok = RegionSpatialProfile(image_path, endpoints, spatial_reqs, spatial_profile); + + ASSERT_TRUE(ok); + ASSERT_EQ(spatial_profile.profiles_size(), 1); + + // Check file/region ids for requirements +} + +TEST_F(RegionSpatialProfileTest, Hdf5LineProfile) { + std::string image_path = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::vector endpoints = {0.0, 0.0, 9.0, 9.0}; + int start(0), end(0), mip(0), width(3); + std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; + + CARTA::SpatialProfileData spatial_profile; + bool ok = RegionSpatialProfile(image_path, endpoints, spatial_reqs, spatial_profile); + + ASSERT_TRUE(ok); + ASSERT_EQ(spatial_profile.profiles_size(), 1); +} + +TEST_F(RegionSpatialProfileTest, FitsHorizontalCutProfile) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {9.0, 5.0, 1.0, 5.0}; // Set line region at y=5 + int start(0), end(0), mip(0), width(1); + std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; + + CARTA::SpatialProfileData spatial_profile; + bool ok = RegionSpatialProfile(image_path, endpoints, spatial_reqs, spatial_profile); + + ASSERT_TRUE(ok); + ASSERT_EQ(spatial_profile.profiles_size(), 1); + + // Profile data + auto profile = spatial_profile.profiles(0); + std::vector profile_data = ProfileValues(profile); + EXPECT_EQ(profile_data.size(), 9); + + // Read image data slice for first channel + // Profile data of horizontal line with width=1 is same as slice + FitsDataReader reader(image_path); + auto image_data = reader.ReadRegion({1, 5, 0}, {10, 6, 1}); + CmpVectors(profile_data, image_data); +} + +TEST_F(RegionSpatialProfileTest, FitsVerticalCutProfile) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {5.0, 9.0, 5.0, 1.0}; // Set line region at x=5 + int start(0), end(0), mip(0), width(1); + std::vector spatial_reqs = {Message::SpatialConfig("y", start, end, mip, width)}; + + CARTA::SpatialProfileData spatial_profile; + bool ok = RegionSpatialProfile(image_path, endpoints, spatial_reqs, spatial_profile); + + ASSERT_TRUE(ok); + ASSERT_EQ(spatial_profile.profiles_size(), 1); + + // Profile data + auto profile = spatial_profile.profiles(0); + std::vector profile_data = ProfileValues(profile); + EXPECT_EQ(profile_data.size(), 9); + + // Read image data slice for first channel + // Profile data of vertical line with width=1 is same as slice + FitsDataReader reader(image_path); + auto image_data = reader.ReadRegion({5, 1, 0}, {6, 10, 1}); + CmpVectors(profile_data, image_data); +} + +TEST_F(RegionSpatialProfileTest, FitsPolylineProfile) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {1.0, 1.0, 9.0, 1.0, 9.0, 5.0}; + int start(0), end(0), mip(0), width(1); + std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; + std::vector spatial_profiles; + + CARTA::SpatialProfileData spatial_profile; + bool ok = RegionSpatialProfile(image_path, endpoints, spatial_reqs, spatial_profile); + ASSERT_TRUE(ok); + ASSERT_EQ(spatial_profile.profiles_size(), 1); + + // Profile data + auto profile = spatial_profile.profiles(0); + std::vector profile_data = ProfileValues(profile); + EXPECT_EQ(profile_data.size(), 13); + + // Read image data slice for first channel + FitsDataReader reader(image_path); + auto line0_data = reader.ReadRegion({1, 1, 0}, {10, 2, 1}); + // Trim line 1: [9, 1] already covered by line 0 + auto line1_data = reader.ReadRegion({9, 2, 0}, {10, 6, 1}); + auto image_data = line0_data; + for (size_t i = 0; i < line1_data.size(); ++i) { + image_data.push_back(line1_data[i]); + } + + // Profile data of polyline with width=1 is same as image data slices + CmpVectors(profile_data, image_data); +} + +TEST_F(RegionSpatialProfileTest, AveragingWidthRange) { + TestAveragingWidthRange(0, false); + TestAveragingWidthRange(1, true); + TestAveragingWidthRange(20, true); + TestAveragingWidthRange(21, false); +} + +TEST_F(RegionSpatialProfileTest, FitsAnnotationLineProfile) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {0.0, 0.0, 9.0, 9.0}; + int start(0), end(0), mip(0), width(3); + std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; + + CARTA::SpatialProfileData spatial_profile; + bool is_annotation(true); + bool ok = RegionSpatialProfile(image_path, endpoints, spatial_reqs, spatial_profile, is_annotation); + + ASSERT_FALSE(ok); +} + +TEST_F(RegionSpatialProfileTest, FitsPointProfile) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {0.0, 0.0}; + int start(0), end(0), mip(0), width(1); + std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; + + CARTA::SpatialProfileData spatial_profile; + bool ok = RegionSpatialProfile(image_path, endpoints, spatial_reqs, spatial_profile); + + ASSERT_TRUE(ok); + ASSERT_EQ(spatial_profile.profiles_size(), 1); + + // Profile data for 10 x 10 image frame + auto profile = spatial_profile.profiles(0); + std::vector profile_data = ProfileValues(profile); + EXPECT_EQ(profile_data.size(), 10); + + // Read image data slice for first channel + // Profile data of point is same as slice + FitsDataReader reader(image_path); + auto image_data = reader.ReadRegion({0, 0, 0}, {10, 1, 1}); + CmpVectors(profile_data, image_data); +} + +TEST_F(RegionSpatialProfileTest, Hdf5PointProfile) { + std::string image_path = FileFinder::Hdf5ImagePath("noise_10px_10px.hdf5"); + std::vector points = {0.0, 0.0}; + int start(0), end(0), mip(0), width(1); + std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; + + CARTA::SpatialProfileData spatial_profile; + bool ok = RegionSpatialProfile(image_path, points, spatial_reqs, spatial_profile); + + ASSERT_TRUE(ok); + ASSERT_EQ(spatial_profile.profiles_size(), 1); + + // Profile data for 10 x 10 image frame + auto profile = spatial_profile.profiles(0); + std::vector profile_data = ProfileValues(profile); + EXPECT_EQ(profile_data.size(), 10); + + // Read image data slice for first channel + // Profile data of point is same as slice + Hdf5DataReader reader(image_path); + auto image_data = reader.ReadRegion({0, 0, 0}, {10, 1, 1}); + CmpVectors(profile_data, image_data); +} + +TEST_F(RegionSpatialProfileTest, FitsPointProfileOutsideImage) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {-2.0, -2.0}; + int start(0), end(0), mip(0), width(1); + std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; + + CARTA::SpatialProfileData spatial_profile; + bool ok = RegionSpatialProfile(image_path, endpoints, spatial_reqs, spatial_profile); + + ASSERT_FALSE(ok); +} + +TEST_F(RegionSpatialProfileTest, FitsAnnotationPointProfile) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector point = {0.0, 0.0}; + int start(0), end(0), mip(0), width(3); + std::vector spatial_reqs = {Message::SpatialConfig("x", start, end, mip, width)}; + + CARTA::SpatialProfileData spatial_profile; + bool is_annotation(true); + bool ok = RegionSpatialProfile(image_path, point, spatial_reqs, spatial_profile, is_annotation); + + ASSERT_FALSE(ok); +} diff --git a/test/TestRegionSpectralProfiles.cc b/test/TestRegionSpectralProfiles.cc new file mode 100644 index 000000000..441ae5f4f --- /dev/null +++ b/test/TestRegionSpectralProfiles.cc @@ -0,0 +1,173 @@ +/* This file is part of the CARTA Image Viewer: https://github.com/CARTAvis/carta-backend + Copyright 2018- Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), + Associated Universities, Inc. (AUI) and the Inter-University Institute for Data Intensive Astronomy (IDIA) + SPDX-License-Identifier: GPL-3.0-or-later +*/ + +#include +#include + +#include "CommonTestUtilities.h" +#include "ImageData/FileLoader.h" +#include "Region/Region.h" +#include "Region/RegionHandler.h" +#include "src/Frame/Frame.h" + +using namespace carta; + +class RegionSpectralProfileTest : public ::testing::Test { +public: + static bool SetRegion(carta::RegionHandler& region_handler, int file_id, int& region_id, const std::vector& points, + std::shared_ptr csys, bool is_annotation) { + std::vector control_points; + for (auto i = 0; i < points.size(); i += 2) { + control_points.push_back(Message::Point(points[i], points[i + 1])); + } + + // Define RegionState and set region + CARTA::RegionType region_type; + if (control_points.size() > 1) { + region_type = is_annotation ? CARTA::ANNPOLYGON : CARTA::POLYGON; + } else { + region_type = is_annotation ? CARTA::ANNPOINT : CARTA::POINT; + } + RegionState region_state(file_id, region_type, control_points, 0.0); + return region_handler.SetRegion(region_id, region_state, csys); + } + + static bool SpectralProfile(const std::string& image_path, const std::vector& points, CARTA::SpectralProfileData& spectral_data, + bool is_annotation = false) { + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + carta::RegionHandler region_handler; + + // Set polygon or point region + int file_id(0), region_id(-1); + auto csys = frame->CoordinateSystem(); + if (!SetRegion(region_handler, file_id, region_id, points, csys, is_annotation)) { + return false; + } + + // Set spectral requirements (Message requests 10 stats types) + auto spectral_req_message = Message::SetSpectralRequirements(file_id, region_id, "z"); + std::vector spectral_requirements = { + spectral_req_message.spectral_profiles().begin(), spectral_req_message.spectral_profiles().end()}; + if (!region_handler.SetSpectralRequirements(region_id, file_id, frame, spectral_requirements)) { + return false; + } + + // Get spectral profile + return region_handler.FillSpectralProfileData( + [&](CARTA::SpectralProfileData profile_data) { spectral_data = profile_data; }, region_id, file_id, false); + } + + static std::vector GetExpectedMeanProfile(std::string& image_path, int num_channels, CARTA::RegionType type) { + // Read image for profile for box region blc (0,0) trc (3,3) or point (3,3) + FitsDataReader reader(image_path); + std::vector profile; + if (type == CARTA::POLYGON) { + for (hsize_t i = 0; i < num_channels; ++i) { + auto channel_data = reader.ReadRegion({0, 0, i}, {4, 4, i + 1}); + double sum = std::accumulate(channel_data.begin(), channel_data.end(), 0.0); + double mean = sum / channel_data.size(); + profile.push_back(mean); + } + } else { + hsize_t chan = num_channels; + auto channel_data = reader.ReadRegion({3, 3, 0}, {4, 4, chan}); + for (float data : channel_data) { + profile.push_back(data); + } + } + return profile; + } +}; + +TEST_F(RegionSpectralProfileTest, TestPolygonSpectralProfile) { + // Box described as 4-corner polygon + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + int num_channels = 10; + std::vector points = {0.0, 0.0, 0.0, 3.0, 3.0, 3.0, 3.0, 0.0}; + CARTA::SpectralProfileData spectral_data; + bool ok = SpectralProfile(image_path, points, spectral_data); + + // Check spectral profiles + ASSERT_TRUE(ok); + ASSERT_EQ(spectral_data.file_id(), 0); + ASSERT_EQ(spectral_data.region_id(), 1); + ASSERT_EQ(spectral_data.stokes(), 0); + ASSERT_EQ(spectral_data.progress(), 1.0); + + auto num_profiles = spectral_data.profiles_size(); + // Expected types set in Message::SetSpectralRequirements + std::vector expected_types = {CARTA::StatsType::NumPixels, CARTA::StatsType::Sum, CARTA::StatsType::FluxDensity, + CARTA::StatsType::Mean, CARTA::StatsType::RMS, CARTA::StatsType::Sigma, CARTA::StatsType::SumSq, CARTA::StatsType::Min, + CARTA::StatsType::Max, CARTA::StatsType::Extrema}; + ASSERT_EQ(num_profiles, expected_types.size()); + + for (int i = 0; i < num_profiles; ++i) { + auto profile = spectral_data.profiles(i); + ASSERT_EQ(profile.coordinate(), "z"); + ASSERT_EQ(profile.stats_type(), expected_types[i]); + auto values = profile.raw_values_fp64(); // string (bytes) + ASSERT_EQ(values.size(), num_channels * 8); // double (8-byte) values + + if (profile.stats_type() == CARTA::StatsType::Mean) { + auto expected_profile = GetExpectedMeanProfile(image_path, num_channels, CARTA::POLYGON); + auto actual_profile = GetSpectralProfileValues(profile); + CmpVectors(actual_profile, expected_profile); + } + } +} + +TEST_F(RegionSpectralProfileTest, TestAnnPolygonSpectralProfile) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector points = {0.0, 0.0, 0.0, 3.0, 3.0, 3.0, 3.0, 0.0}; + CARTA::SpectralProfileData spectral_data; + bool ok = SpectralProfile(image_path, points, spectral_data, true); + ASSERT_FALSE(ok); +} + +TEST_F(RegionSpectralProfileTest, TestPointSpectralProfile) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + int num_channels = 10; + std::vector points = {3.0, 3.0}; + CARTA::SpectralProfileData spectral_data; + bool ok = SpectralProfile(image_path, points, spectral_data); + + // Check spectral profiles + ASSERT_TRUE(ok); + ASSERT_EQ(spectral_data.file_id(), 0); + ASSERT_EQ(spectral_data.region_id(), 1); + ASSERT_EQ(spectral_data.stokes(), 0); + ASSERT_EQ(spectral_data.progress(), 1.0); + + auto num_profiles = spectral_data.profiles_size(); + // Expected types set in Message::SetSpectralRequirements + std::vector expected_types = {CARTA::StatsType::NumPixels, CARTA::StatsType::Sum, CARTA::StatsType::FluxDensity, + CARTA::StatsType::Mean, CARTA::StatsType::RMS, CARTA::StatsType::Sigma, CARTA::StatsType::SumSq, CARTA::StatsType::Min, + CARTA::StatsType::Max, CARTA::StatsType::Extrema}; + ASSERT_EQ(num_profiles, expected_types.size()); + + for (int i = 0; i < num_profiles; ++i) { + auto profile = spectral_data.profiles(i); + ASSERT_EQ(profile.coordinate(), "z"); + ASSERT_EQ(profile.stats_type(), expected_types[i]); + auto values = profile.raw_values_fp64(); // string (bytes) + ASSERT_EQ(values.size(), num_channels * 8); // double (8-byte) values + + if (profile.stats_type() == CARTA::StatsType::Mean) { + auto expected_profile = GetExpectedMeanProfile(image_path, num_channels, CARTA::POINT); + auto actual_profile = GetSpectralProfileValues(profile); + CmpVectors(actual_profile, expected_profile); + } + } +} + +TEST_F(RegionSpectralProfileTest, TestAnnPointSpectralProfile) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector points = {3.0, 3.0}; + CARTA::SpectralProfileData spectral_data; + bool ok = SpectralProfile(image_path, points, spectral_data, true); + ASSERT_FALSE(ok); +} diff --git a/test/TestRegionStats.cc b/test/TestRegionStats.cc new file mode 100644 index 000000000..c3b6c7399 --- /dev/null +++ b/test/TestRegionStats.cc @@ -0,0 +1,108 @@ +/* This file is part of the CARTA Image Viewer: https://github.com/CARTAvis/carta-backend + Copyright 2018- Academia Sinica Institute of Astronomy and Astrophysics (ASIAA), + Associated Universities, Inc. (AUI) and the Inter-University Institute for Data Intensive Astronomy (IDIA) + SPDX-License-Identifier: GPL-3.0-or-later +*/ + +#include +#include + +#include "CommonTestUtilities.h" +#include "ImageData/FileLoader.h" +#include "Region/Region.h" +#include "Region/RegionHandler.h" +#include "src/Frame/Frame.h" + +using namespace carta; + +class RegionStatsTest : public ::testing::Test { +public: + static bool SetRegion(carta::RegionHandler& region_handler, int file_id, int& region_id, const std::vector& points, + std::shared_ptr csys, bool is_annotation) { + std::vector control_points; + for (auto i = 0; i < points.size(); i += 2) { + control_points.push_back(Message::Point(points[i], points[i + 1])); + } + + // Define RegionState for line region and set region (region_id updated) + auto npoints(control_points.size()); + CARTA::RegionType region_type = CARTA::RegionType::POLYGON; + if (is_annotation) { + region_type = CARTA::RegionType::ANNPOLYGON; + } + RegionState region_state(file_id, region_type, control_points, 0.0); + return region_handler.SetRegion(region_id, region_state, csys); + } + + static bool RegionStats(const std::string& image_path, const std::vector& endpoints, CARTA::RegionStatsData& region_stats, + bool is_annotation = false) { + std::shared_ptr loader(carta::FileLoader::GetLoader(image_path)); + std::shared_ptr frame(new Frame(0, loader, "0")); + carta::RegionHandler region_handler; + + // Set polygon region + int file_id(0), region_id(-1); + auto csys = frame->CoordinateSystem(); + if (!SetRegion(region_handler, file_id, region_id, endpoints, csys, is_annotation)) { + return false; + } + + // Set stats requirements + auto stats_req_message = Message::SetStatsRequirements(file_id, region_id); + std::vector stats_configs = { + stats_req_message.stats_configs().begin(), stats_req_message.stats_configs().end()}; + if (!region_handler.SetStatsRequirements(region_id, file_id, frame, stats_configs)) { + return false; + } + + // Get stats + return region_handler.FillRegionStatsData( + [&](CARTA::RegionStatsData stats_data) { region_stats = stats_data; }, region_id, file_id); + } +}; + +TEST_F(RegionStatsTest, TestFitsRegionStats) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {1.0, 1.0, 1.0, 4.0, 4.0, 4.0, 4.0, 1.0}; + CARTA::RegionStatsData stats_data; + bool ok = RegionStats(image_path, endpoints, stats_data); + + // Check stats fields + ASSERT_TRUE(ok); + ASSERT_EQ(stats_data.file_id(), 0); + ASSERT_EQ(stats_data.region_id(), 1); + ASSERT_EQ(stats_data.channel(), 0); + ASSERT_EQ(stats_data.stokes(), 0); + ASSERT_GT(stats_data.statistics_size(), 0); + + // Calc expected stats from image data + FitsDataReader reader(image_path); + auto image_data = reader.ReadRegion({1, 1, 0}, {5, 5, 1}); + double expected_sum = std::accumulate(image_data.begin(), image_data.end(), 0.0); + double expected_mean = expected_sum / image_data.size(); + double expected_min = *std::min_element(image_data.begin(), image_data.end()); + double expected_max = *std::max_element(image_data.begin(), image_data.end()); + + // Check some stats + for (size_t i = 0; i < stats_data.statistics_size(); ++i) { + if (stats_data.statistics(i).stats_type() == CARTA::StatsType::NumPixels) { + ASSERT_DOUBLE_EQ(stats_data.statistics(i).value(), (double)image_data.size()); + } else if (stats_data.statistics(i).stats_type() == CARTA::StatsType::Sum) { + ASSERT_DOUBLE_EQ(stats_data.statistics(i).value(), expected_sum); + } else if (stats_data.statistics(i).stats_type() == CARTA::StatsType::Mean) { + ASSERT_DOUBLE_EQ(stats_data.statistics(i).value(), expected_mean); + } else if (stats_data.statistics(i).stats_type() == CARTA::StatsType::Min) { + ASSERT_DOUBLE_EQ(stats_data.statistics(i).value(), expected_min); + } else if (stats_data.statistics(i).stats_type() == CARTA::StatsType::Max) { + ASSERT_DOUBLE_EQ(stats_data.statistics(i).value(), expected_max); + } + } +} + +TEST_F(RegionStatsTest, TestFitsAnnotationRegionStats) { + std::string image_path = FileFinder::FitsImagePath("noise_3d.fits"); + std::vector endpoints = {0.0, 0.0, 0.0, 3.0, 3.0, 3.0, 3.0, 0.0}; + CARTA::RegionStatsData stats_data; + bool ok = RegionStats(image_path, endpoints, stats_data, true); + ASSERT_FALSE(ok); +} From f17756b72ee89426f30ea619b671d77d6deaac31 Mon Sep 17 00:00:00 2001 From: pford Date: Tue, 9 Apr 2024 15:01:43 -0600 Subject: [PATCH 3/5] Improve region spatial profile and PV performance (#1353) * Use cached image data for region data when current channel+stokes * Use openmp for RegionHandler line profiles and calculate stats locally * Fix slicer data from image cache for >2D image * Fix image cache slicer and line profiles return status * Fix formatting * Update changelog * Use counter instead of for loop variable for pv image progress * Convert rotbox corners to LCPolygon if no distortion * Set min length for approximated polygon segments between points * Improve LCRegion for rotbox in reference image * Fix caching rotbox polygon region * Remove unused variable * Fix bugs from rotbox LCPolygon implementation * Fix polygon approximation to decrease error * Formatting pass --------- Co-authored-by: Pam Harris --- CHANGELOG.md | 1 + src/Frame/Frame.cc | 152 ++++++++++++++++++++++++------------ src/Frame/Frame.h | 2 +- src/Region/Region.cc | 1 - src/Region/RegionHandler.cc | 50 ++++++++---- 5 files changed, 141 insertions(+), 65 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 3f7162879..31cc83ffb 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -14,6 +14,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0 * Improved file IDs for generated images ([#1224](https://github.com/CARTAvis/carta-frontend/issues/1224)). * Store global settings in a singleton class ([#1302](https://github.com/CARTAvis/carta-backend/issues/1302)). * Move Region code to RegionConverter for maintenance and performance ([#1347](https://github.com/CARTAvis/carta-backend/issues/1347)). +* Improve performance of region spatial profiles and PV image generation ([#1339](https://github.com/CARTAvis/carta-backend/issues/1339)). ## [4.1.0] diff --git a/src/Frame/Frame.cc b/src/Frame/Frame.cc index 8bf80d34a..bfdfb1acc 100644 --- a/src/Frame/Frame.cc +++ b/src/Frame/Frame.cc @@ -182,6 +182,11 @@ int Frame::StokesAxis() { return _stokes_axis; } +bool Frame::IsCurrentZStokes(const StokesSource& stokes_source) { + return (stokes_source.z_range.from == stokes_source.z_range.to) && (stokes_source.z_range.from == CurrentZ()) && + (stokes_source.stokes == CurrentStokes()); +} + bool Frame::GetBeams(std::vector& beams) { std::string error; bool beams_ok = _loader->GetBeams(beams, error); @@ -331,12 +336,12 @@ bool Frame::SetImageChannels(int new_z, int new_stokes, std::string& message) { bool z_ok(CheckZ(new_z)); bool stokes_ok(CheckStokes(new_stokes)); if (z_ok && stokes_ok) { - _z_index = new_z; - _stokes_index = new_stokes; - // invalidate the image cache InvalidateImageCache(); + _z_index = new_z; + _stokes_index = new_stokes; + if (!(_loader->UseTileCache() && _loader->HasMip(2)) || IsComputedStokes(_stokes_index)) { // Reload the full channel cache for loaders which use it FillImageCache(); @@ -366,7 +371,6 @@ bool Frame::SetCursor(float x, float y) { bool Frame::FillImageCache() { // get image data for z, stokes - bool write_lock(true); queuing_rw_mutex_scoped cache_lock(&_cache_mutex, write_lock); @@ -1605,68 +1609,120 @@ bool Frame::GetSlicerSubImage(const StokesSlicer& stokes_slicer, casacore::SubIm bool Frame::GetRegionData(const StokesRegion& stokes_region, std::vector& data, bool report_performance) { // Get image data with a region applied Timer t; - casacore::SubImage sub_image; - bool subimage_ok = GetRegionSubImage(stokes_region, sub_image); + std::vector region_mask; - if (!subimage_ok) { - return false; - } - - casacore::IPosition subimage_shape = sub_image.shape(); - if (subimage_shape.empty()) { - return false; + if (IsCurrentZStokes(stokes_region.stokes_source)) { + try { + // Slice cached image data using LCRegion bounding box + casacore::Slicer bounding_box = stokes_region.image_region.asLCRegion().boundingBox(); + StokesSlicer stokes_slicer(stokes_region.stokes_source, bounding_box); + data.resize(bounding_box.length().product()); + + if (GetSlicerData(stokes_slicer, data.data())) { + // Next get the LCRegion as a mask (LCRegion is a Lattice) + casacore::Array tmpmask = stokes_region.image_region.asLCRegion().get(); + region_mask = tmpmask.tovector(); + } else { + data.clear(); + } + } catch (const casacore::AipsError& err) { + // ImageRegion underlying region was not LCRegion + data.clear(); + } } - try { - casacore::IPosition start(subimage_shape.size(), 0); - casacore::IPosition count(subimage_shape); - casacore::Slicer slicer(start, count); // entire subimage - bool is_computed_stokes(!stokes_region.stokes_source.IsOriginalImage()); + if (data.empty()) { + // Apply region to image to get SubImage data + casacore::SubImage sub_image; + bool subimage_ok = GetRegionSubImage(stokes_region, sub_image); - // Get image data - std::unique_lock ulock(_image_mutex); - if (_loader->IsGenerated() || is_computed_stokes) { // For the image in memory - casacore::Array tmp; - sub_image.doGetSlice(tmp, slicer); - data = tmp.tovector(); - } else { - data.resize(subimage_shape.product()); // must size correctly before sharing - casacore::Array tmp(subimage_shape, data.data(), casacore::StorageInitPolicy::SHARE); - sub_image.doGetSlice(tmp, slicer); + if (!subimage_ok) { + return false; } - // Get mask that defines region in subimage bounding box - casacore::Array tmpmask; - sub_image.doGetMaskSlice(tmpmask, slicer); - ulock.unlock(); + casacore::IPosition subimage_shape = sub_image.shape(); + if (subimage_shape.empty()) { + return false; + } - // Apply mask to data - std::vector datamask = tmpmask.tovector(); - for (size_t i = 0; i < data.size(); ++i) { - if (!datamask[i]) { - data[i] = NAN; + try { + casacore::IPosition start(subimage_shape.size(), 0); + casacore::IPosition count(subimage_shape); + casacore::Slicer slicer(start, count); // entire subimage + bool is_computed_stokes(!stokes_region.stokes_source.IsOriginalImage()); + + // Get image data and mask, with image mutex locked + std::unique_lock ulock(_image_mutex); + casacore::Array tmpdata; + if (_loader->IsGenerated() || is_computed_stokes) { // For the image in memory + sub_image.doGetSlice(tmpdata, slicer); + data = tmpdata.tovector(); + } else { + data.resize(subimage_shape.product()); // must size correctly before sharing + tmpdata = casacore::Array(subimage_shape, data.data(), casacore::StorageInitPolicy::SHARE); + sub_image.doGetSlice(tmpdata, slicer); } + + // Get mask that defines region in subimage bounding box + casacore::Array tmpmask; + sub_image.doGetMaskSlice(tmpmask, slicer); + ulock.unlock(); + region_mask = tmpmask.tovector(); + } catch (const casacore::AipsError& err) { + data.clear(); + return false; } + } - if (report_performance) { - spdlog::performance("Get region subimage data in {:.3f} ms", t.Elapsed().ms()); + // Apply mask to data + for (size_t i = 0; i < data.size(); ++i) { + if (!region_mask[i]) { + data[i] = NAN; } + } - return true; - } catch (casacore::AipsError& err) { - data.clear(); + if (report_performance) { + spdlog::performance("Get region subimage data in {:.3f} ms", t.Elapsed().ms()); } - return false; + return true; } bool Frame::GetSlicerData(const StokesSlicer& stokes_slicer, float* data) { - // Get image data with a slicer applied + // Get image data with a slicer applied; data must be correctly resized + bool data_ok(false); casacore::Array tmp(stokes_slicer.slicer.length(), data, casacore::StorageInitPolicy::SHARE); - std::unique_lock ulock(_image_mutex); - bool data_ok = _loader->GetSlice(tmp, stokes_slicer); - _loader->CloseImageIfUpdated(); - ulock.unlock(); + + if (_image_cache_valid && IsCurrentZStokes(stokes_slicer.stokes_source)) { + // Slice image cache + auto cache_shape = ImageShape(); + auto slicer_start = stokes_slicer.slicer.start(); + auto slicer_end = stokes_slicer.slicer.end(); + + // Adjust cache shape and slicer for single channel and stokes + if (_spectral_axis >= 0) { + cache_shape(_spectral_axis) = 1; + slicer_start(_spectral_axis) = 0; + slicer_end(_spectral_axis) = 0; + } + if (_stokes_axis >= 0) { + cache_shape(_stokes_axis) = 1; + slicer_start(_stokes_axis) = 0; + slicer_end(_stokes_axis) = 0; + } + casacore::Slicer cache_slicer(slicer_start, slicer_end, casacore::Slicer::endIsLast); + + queuing_rw_mutex_scoped cache_lock(&_cache_mutex, false); // read lock + casacore::Array image_cache_as_array(cache_shape, _image_cache.get(), casacore::StorageInitPolicy::SHARE); + tmp = image_cache_as_array(cache_slicer); + data_ok = true; + } else { + // Use loader to slice image + std::unique_lock ulock(_image_mutex); + data_ok = _loader->GetSlice(tmp, stokes_slicer); + _loader->CloseImageIfUpdated(); + ulock.unlock(); + } return data_ok; } diff --git a/src/Frame/Frame.h b/src/Frame/Frame.h index 0580f0496..7d56cce54 100644 --- a/src/Frame/Frame.h +++ b/src/Frame/Frame.h @@ -105,6 +105,7 @@ class Frame { size_t NumStokes(); // if no stokes axis, nstokes=1 int CurrentZ(); int CurrentStokes(); + bool IsCurrentZStokes(const StokesSource& stokes_source); int SpectralAxis(); int StokesAxis(); bool GetBeams(std::vector& beams); @@ -296,7 +297,6 @@ class Frame { ContourSettings _contour_settings; // Image data cache and mutex - // std::vector _image_cache; // image data for current z, stokes long long int _image_cache_size; std::unique_ptr _image_cache; bool _image_cache_valid; // cached image data is valid for current z and stokes diff --git a/src/Region/Region.cc b/src/Region/Region.cc index 157d7334f..daf47ffd3 100644 --- a/src/Region/Region.cc +++ b/src/Region/Region.cc @@ -12,7 +12,6 @@ #include #include "RegionConverter.h" -// #include "Util/Image.h" using namespace carta; diff --git a/src/Region/RegionHandler.cc b/src/Region/RegionHandler.cc index 609daf957..6087aa2f3 100644 --- a/src/Region/RegionHandler.cc +++ b/src/Region/RegionHandler.cc @@ -1248,7 +1248,6 @@ bool RegionHandler::CalculatePvPreviewImage(int frame_id, int preview_id, bool q // Get box region LCRegion and mask bool cancel(false); - casacore::Slicer box_bounding_box; auto box_lc_region = ApplyRegionToFile(box_region_id, frame_id); if (!box_lc_region) { @@ -2332,6 +2331,7 @@ bool RegionHandler::FillLineSpatialProfileData(int file_id, int region_id, std:: profile, coordinate, mip, axis_type, crpix, crval, cdelt, unit); cb(profile_message); }); + spdlog::performance("Fill line spatial profile in {:.3f} ms", t.Elapsed().ms()); } spdlog::performance("Line spatial data in {:.3f} ms", t.Elapsed().ms()); @@ -2467,25 +2467,37 @@ bool RegionHandler::GetLineProfiles(int file_id, int region_id, int width, const if (line_box_regions.GetLineBoxRegions(line_region_state, line_coord_sys, width, increment, box_regions, message)) { auto t_start = std::chrono::high_resolution_clock::now(); auto num_profiles = box_regions.size(); + size_t iprofile; + // Use completed profiles (not iprofile) for progress. + // iprofile not in order so progress is uneven. + size_t completed_profiles(0); + + // Return this to column 0 when uncomment (moved for format check): + // #pragma omp parallel for private(iprofile) shared(progress, t_start, completed_profiles) + for (iprofile = 0; iprofile < num_profiles; ++iprofile) { + if (cancelled) { + continue; + } - for (size_t iprofile = 0; iprofile < num_profiles; ++iprofile) { // Frame/region closing, or line changed if (CancelLineProfiles(region_id, file_id, line_region_state)) { cancelled = true; - return false; } // PV generator: check if user canceled if (per_z && _stop_pv[file_id]) { spdlog::debug("Stopping line profiles: PV generator cancelled"); cancelled = true; - return false; } // Line spatial profile: check if requirements removed if (!per_z && !HasSpatialRequirements(region_id, file_id, coordinate, width)) { cancelled = true; - return false; + } + + if (cancelled) { + profiles.resize(); + continue; } // Get mean profile for requested file_id and log number of pixels in region @@ -2509,7 +2521,7 @@ bool RegionHandler::GetLineProfiles(int file_id, int region_id, int width, const profiles.row(iprofile) = region_profile; } - progress = float(iprofile + 1) / float(num_profiles); + progress = float(++completed_profiles) / float(num_profiles); if (per_z) { // Update progress if time interval elapsed @@ -2524,7 +2536,7 @@ bool RegionHandler::GetLineProfiles(int file_id, int region_id, int width, const } } - return (progress >= 1.0) && !allEQ(profiles, NAN); + return (!cancelled) && (progress >= 1.0) && !allEQ(profiles, NAN); } bool RegionHandler::CancelLineProfiles(int region_id, int file_id, RegionState& region_state) { @@ -2559,9 +2571,9 @@ casacore::Vector RegionHandler::GetTemporaryRegionProfile(int region_idx, } std::lock_guard guard(_line_profile_mutex); + // Set temporary region int region_id(TEMP_REGION_ID); SetRegion(region_id, region_state, reference_csys); - if (!RegionSet(region_id, true)) { return profile; } @@ -2602,19 +2614,27 @@ casacore::Vector RegionHandler::GetTemporaryRegionProfile(int region_idx, }); } else { // Use BasicStats to get num_pixels and mean for current channel and stokes - // Get region + // Get region as LCRegion StokesRegion stokes_region; std::shared_ptr lc_region; std::shared_lock frame_lock(_frames.at(file_id)->GetActiveTaskMutex()); if (ApplyRegionToFile(region_id, file_id, z_range, stokes_index, lc_region, stokes_region)) { - // Get region data (report_performance = false, too much output) + // Get region data by applying LCRegion to image std::vector region_data; if (_frames.at(file_id)->GetRegionData(stokes_region, region_data, false)) { - // Get BasicStats - BasicStats basic_stats; - CalcBasicStats(basic_stats, region_data.data(), region_data.size()); - num_pixels = basic_stats.num_pixels; - profile[0] = basic_stats.mean; + // Very small region, just calc needed stats here + num_pixels = 0; + double sum(0.0); + for (size_t i = 0; i < region_data.size(); ++i) { + float val(region_data[i]); + if (std::isfinite(val)) { + num_pixels++; + sum += (double)val; + } + } + if (num_pixels > 0) { + profile[0] = sum / num_pixels; + } } } frame_lock.unlock(); From be2999f5921b39d2abba2ac4812b92da487b586b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Adrianna=20Pi=C5=84ska?= Date: Thu, 18 Apr 2024 11:29:58 +0200 Subject: [PATCH 4/5] Move to GitHub project (#1372) * Updated checkbox in PR template * Bumped protobuf --- .github/pull_request_template.md | 2 +- carta-protobuf | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/.github/pull_request_template.md b/.github/pull_request_template.md index b378bba4e..515217888 100644 --- a/.github/pull_request_template.md +++ b/.github/pull_request_template.md @@ -13,4 +13,4 @@ - [ ] protobuf updated to the latest dev commit / no protobuf update needed - [ ] protobuf version bumped / protobuf version not bumped - [ ] added reviewers and assignee -- [ ] added ZenHub estimate, milestone, and release +- [ ] GitHub Project estimate added diff --git a/carta-protobuf b/carta-protobuf index 701220443..3219acafa 160000 --- a/carta-protobuf +++ b/carta-protobuf @@ -1 +1 @@ -Subproject commit 701220443c977641b1884d2fed97a5b71ed46b9d +Subproject commit 3219acafa7bc454fe8d448f4dd99cef29f5ac5e8 From 848e98281d92783c1cac1ef2d26206656084fa93 Mon Sep 17 00:00:00 2001 From: Cheng-Chin Chiang Date: Tue, 23 Apr 2024 16:15:50 +0800 Subject: [PATCH 5/5] Bump the submodule carta-protobuf (#1367) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit * Include absl libraries for the new version of protobuf * bump the submodule carta-protobuf * Correct the identification of the new protobuf version * Bumped protobuf version --------- Co-authored-by: Adrianna Pińska --- CMakeLists.txt | 11 +++++++++++ carta-protobuf | 2 +- 2 files changed, 12 insertions(+), 1 deletion(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index 5b8c7d2de..ae8a92aa3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -109,6 +109,17 @@ endif () FIND_PACKAGE(HDF5 REQUIRED COMPONENTS CXX) FIND_PACKAGE(Protobuf REQUIRED) INCLUDE_DIRECTORIES(${PROTOBUF_INCLUDE_DIR}) + +if (Protobuf_VERSION VERSION_GREATER_EQUAL "4.25.3") + FIND_PACKAGE(absl REQUIRED) + set(PROTOBUF_LIBRARY + ${PROTOBUF_LIBRARY} + absl_hash + absl_log_internal_message + absl_log_internal_nullguard) + message(STATUS "Newer protobuf version (${Protobuf_VERSION}) includes abseil libraries") +endif () + FIND_PACKAGE(Threads) INCLUDE_DIRECTORIES(${HDF5_INCLUDE_DIR}) diff --git a/carta-protobuf b/carta-protobuf index 3219acafa..202c31d25 160000 --- a/carta-protobuf +++ b/carta-protobuf @@ -1 +1 @@ -Subproject commit 3219acafa7bc454fe8d448f4dd99cef29f5ac5e8 +Subproject commit 202c31d25ec4bda97cbb9bd71fee589541872847