diff --git a/src/Frame/Frame.cc b/src/Frame/Frame.cc index 429e31d23..abab881b9 100644 --- a/src/Frame/Frame.cc +++ b/src/Frame/Frame.cc @@ -344,7 +344,9 @@ bool Frame::FillImageCache() { auto t_start_set_image_cache = std::chrono::high_resolution_clock::now(); casacore::Slicer section = GetImageSlicer(AxisRange(_z_index), _stokes_index); - if (!GetSlicerData(section, _image_cache)) { + _image_cache_size = section.length().product(); + _image_cache = std::shared_ptr(new float[_image_cache_size]); + if (!GetSlicerData(section, _image_cache.get())) { spdlog::error("Session {}: {}", _session_id, "Loading image cache failed."); return false; } @@ -368,7 +370,8 @@ void Frame::InvalidateImageCache() { void Frame::GetZMatrix(std::vector& z_matrix, size_t z, size_t stokes) { // fill matrix for given z and stokes casacore::Slicer section = GetImageSlicer(AxisRange(z), stokes); - GetSlicerData(section, z_matrix); + z_matrix.resize(section.length().product()); + GetSlicerData(section, z_matrix.data()); } // **************************************************** @@ -412,10 +415,10 @@ bool Frame::GetRasterData(std::vector& image_data, CARTA::ImageBounds& bo if (mean_filter && mip > 1) { // Perform down-sampling by calculating the mean for each MIPxMIP block BlockSmooth( - _image_cache.data(), image_data.data(), num_image_columns, num_image_rows, row_length_region, num_rows_region, x, y, mip); + _image_cache.get(), image_data.data(), num_image_columns, num_image_rows, row_length_region, num_rows_region, x, y, mip); } else { // Nearest neighbour filtering - NearestNeighbor(_image_cache.data(), image_data.data(), num_image_columns, row_length_region, num_rows_region, x, y, mip); + NearestNeighbor(_image_cache.get(), image_data.data(), num_image_columns, row_length_region, num_rows_region, x, y, mip); } auto t_end_raster_data_filter = std::chrono::high_resolution_clock::now(); @@ -591,7 +594,7 @@ bool Frame::ContourImage(ContourCallback& partial_contour_callback) { queuing_rw_mutex_scoped cache_lock(&_cache_mutex, false); if (_contour_settings.smoothing_mode == CARTA::SmoothingMode::NoSmoothing || _contour_settings.smoothing_factor <= 1) { - TraceContours(_image_cache.data(), _width, _height, scale, offset, _contour_settings.levels, vertex_data, index_data, + TraceContours(_image_cache.get(), _width, _height, scale, offset, _contour_settings.levels, vertex_data, index_data, _contour_settings.chunk_size, partial_contour_callback); return true; } else if (_contour_settings.smoothing_mode == CARTA::SmoothingMode::GaussianBlur) { @@ -604,8 +607,8 @@ bool Frame::ContourImage(ContourCallback& partial_contour_callback) { int64_t dest_width = _width - (2 * kernel_width); int64_t dest_height = _height - (2 * kernel_width); std::unique_ptr dest_array(new float[dest_width * dest_height]); - smooth_successful = GaussianSmooth(_image_cache.data(), dest_array.get(), source_width, source_height, dest_width, dest_height, - _contour_settings.smoothing_factor); + smooth_successful = GaussianSmooth( + _image_cache.get(), dest_array.get(), source_width, source_height, dest_width, dest_height, _contour_settings.smoothing_factor); // Can release lock early, as we're no longer using the image cache cache_lock.release(); if (smooth_successful) { @@ -829,11 +832,11 @@ bool Frame::GetBasicStats(int z, int stokes, BasicStats& stats) { if ((z == CurrentZ()) && (stokes == CurrentStokes())) { // calculate histogram from image cache - if (_image_cache.empty() && !FillImageCache()) { + if ((_image_cache_size == 0) && !FillImageCache()) { // cannot calculate return false; } - CalcBasicStats(_image_cache, stats); + CalcBasicStats(stats, _image_cache.get(), _image_cache_size); _image_basic_stats[cache_key] = stats; return true; } @@ -841,7 +844,7 @@ bool Frame::GetBasicStats(int z, int stokes, BasicStats& stats) { // calculate histogram from given z/stokes data std::vector data; GetZMatrix(data, z, stokes); - CalcBasicStats(data, stats); + CalcBasicStats(stats, data.data(), data.size()); // cache results _image_basic_stats[cache_key] = stats; @@ -897,17 +900,17 @@ bool Frame::CalculateHistogram(int region_id, int z, int stokes, int num_bins, B if ((z == CurrentZ()) && (stokes == CurrentStokes())) { // calculate histogram from current image cache - if (_image_cache.empty() && !FillImageCache()) { + if ((_image_cache_size == 0) && !FillImageCache()) { return false; } bool write_lock(false); queuing_rw_mutex_scoped cache_lock(&_cache_mutex, write_lock); - hist = CalcHistogram(num_bins, stats, _image_cache); + hist = CalcHistogram(num_bins, stats, _image_cache.get(), _image_cache_size); } else { // calculate histogram for z/stokes data std::vector data; GetZMatrix(data, z, stokes); - hist = CalcHistogram(num_bins, stats, data); + hist = CalcHistogram(num_bins, stats, data.data(), data.size()); } // cache image histogram @@ -1115,8 +1118,9 @@ bool Frame::FillSpatialProfileData(PointXy point, std::vector data; - if (GetSlicerData(section, data)) { + const auto N = section.length().product(); + std::shared_ptr data(new float[N]); // zero initialization + if (GetSlicerData(section, data.get())) { cursor_value = data[0]; } } @@ -1274,7 +1278,8 @@ bool Frame::FillSpatialProfileData(PointXy point, std::vector buffer; - if (!GetSlicerData(slicer, buffer)) { + const auto N = slicer.length().product(); + std::shared_ptr buffer(new float[N]); + if (!GetSlicerData(slicer, buffer.get())) { return false; } // copy buffer to spectral_data - memcpy(&spectral_data[start(_z_axis)], buffer.data(), nz * sizeof(float)); + memcpy(&spectral_data[start(_z_axis)], buffer.get(), nz * sizeof(float)); // update start z and determine progress start(_z_axis) += nz; @@ -1662,10 +1668,8 @@ bool Frame::GetRegionData(const casacore::LattRegionHolder& region, std::vector< return false; } -bool Frame::GetSlicerData(const casacore::Slicer& slicer, std::vector& data) { - // Get image data with a slicer applied - data.resize(slicer.length().product()); // must have vector the right size before share it with Array - casacore::Array tmp(slicer.length(), data.data(), casacore::StorageInitPolicy::SHARE); +bool Frame::GetSlicerData(const casacore::Slicer& slicer, float* data) { + casacore::Array tmp(slicer.length(), data, casacore::StorageInitPolicy::SHARE); std::unique_lock ulock(_image_mutex); bool data_ok = _loader->GetSlice(tmp, slicer); _loader->CloseImageIfUpdated(); diff --git a/src/Frame/Frame.h b/src/Frame/Frame.h index 1f9fd77ef..ae073970f 100644 --- a/src/Frame/Frame.h +++ b/src/Frame/Frame.h @@ -170,7 +170,7 @@ class Frame { casacore::IPosition GetRegionShape(const casacore::LattRegionHolder& region); // Returns data vector bool GetRegionData(const casacore::LattRegionHolder& region, std::vector& data); - bool GetSlicerData(const casacore::Slicer& slicer, std::vector& data); + bool GetSlicerData(const casacore::Slicer& slicer, float* data); // Returns stats_values map for spectral profiles and stats data bool GetRegionStats(const casacore::LattRegionHolder& region, const std::vector& required_stats, bool per_z, std::map>& stats_values); @@ -277,12 +277,14 @@ class Frame { ContourSettings _contour_settings; // Image data cache and mutex - std::vector _image_cache; // image data for current z, stokes - bool _image_cache_valid; // cached image data is valid for current z and stokes - queuing_rw_mutex _cache_mutex; // allow concurrent reads but lock for write - std::mutex _image_mutex; // only one disk access at a time - bool _cache_loaded; // channel cache is set - TileCache _tile_cache; // cache for full-resolution image tiles + // std::vector _image_cache; // image data for current z, stokes + long long int _image_cache_size; + std::shared_ptr _image_cache; + bool _image_cache_valid; // cached image data is valid for current z and stokes + queuing_rw_mutex _cache_mutex; // allow concurrent reads but lock for write + std::mutex _image_mutex; // only one disk access at a time + bool _cache_loaded; // channel cache is set + TileCache _tile_cache; // cache for full-resolution image tiles std::mutex _ignore_interrupt_X_mutex; std::mutex _ignore_interrupt_Y_mutex; diff --git a/src/ImageStats/BasicStatsCalculator.h b/src/ImageStats/BasicStatsCalculator.h index d49c9ec61..41a585903 100644 --- a/src/ImageStats/BasicStatsCalculator.h +++ b/src/ImageStats/BasicStatsCalculator.h @@ -32,13 +32,14 @@ class BasicStatsCalculator { T _min_val, _max_val; double _sum, _sum_squares; size_t _num_pixels; - const std::vector& _data; + const T* _data; + size_t _data_size; public: - BasicStatsCalculator(const std::vector& data); + BasicStatsCalculator(const T* data, size_t data_size); void join(BasicStatsCalculator& other); // NOLINT - void reduce(const size_t start, const size_t end); + void reduce(); BasicStats GetStats() const; }; diff --git a/src/ImageStats/BasicStatsCalculator.tcc b/src/ImageStats/BasicStatsCalculator.tcc index c7de755f1..3df265f8a 100644 --- a/src/ImageStats/BasicStatsCalculator.tcc +++ b/src/ImageStats/BasicStatsCalculator.tcc @@ -7,6 +7,8 @@ #ifndef CARTA_BACKEND_IMAGESTATS_BASICSTATSCALCULATOR_TCC_ #define CARTA_BACKEND_IMAGESTATS_BASICSTATSCALCULATOR_TCC_ +#include "Logger/Logger.h" + #include namespace carta { @@ -41,19 +43,20 @@ BasicStats::BasicStats() sumSq(0) {} template -BasicStatsCalculator::BasicStatsCalculator(const std::vector& data) +BasicStatsCalculator::BasicStatsCalculator(const T* data, size_t data_size) : _min_val(std::numeric_limits::max()), _max_val(std::numeric_limits::lowest()), _sum(0), _sum_squares(0), _num_pixels(0), - _data(data) {} + _data(data), + _data_size(data_size) {} template -void BasicStatsCalculator::reduce(const size_t start, const size_t end) { +void BasicStatsCalculator::reduce() { size_t i; #pragma omp parallel for private(i) shared(_data) reduction(min: _min_val) reduction(max:_max_val) reduction(+:_num_pixels) reduction(+:_sum) reduction(+:_sum_squares) - for (i = start; i < end; i++) { + for (i = 0; i < _data_size; i++) { T val = _data[i]; if (std::isfinite(val)) { if (val < _min_val) { diff --git a/src/ImageStats/Histogram.cc b/src/ImageStats/Histogram.cc index e621767bd..6d302e832 100644 --- a/src/ImageStats/Histogram.cc +++ b/src/ImageStats/Histogram.cc @@ -15,13 +15,13 @@ using namespace carta; -Histogram::Histogram(int num_bins, float min_value, float max_value, const std::vector& data) +Histogram::Histogram(int num_bins, float min_value, float max_value, const float* data, const size_t data_size) : _bin_width((max_value - min_value) / num_bins), _min_val(min_value), _max_val(max_value), _bin_center(min_value + (_bin_width * 0.5)), _histogram_bins(num_bins, 0) { - Fill(data); + Fill(data, data_size); } Histogram::Histogram(const Histogram& h) @@ -45,9 +45,9 @@ bool Histogram::Add(const Histogram& h) { return true; } -void Histogram::Fill(const std::vector& data) { +void Histogram::Fill(const float* data, const size_t data_size) { std::vector temp_bins; - const auto num_elements = data.size(); + const auto num_elements = data_size; const size_t num_bins = GetNbins(); ThreadManager::ApplyThreadLimit(); #pragma omp parallel diff --git a/src/ImageStats/Histogram.h b/src/ImageStats/Histogram.h index 8d74bd688..935208266 100644 --- a/src/ImageStats/Histogram.h +++ b/src/ImageStats/Histogram.h @@ -19,12 +19,12 @@ class Histogram { float _bin_center; // bin center std::vector _histogram_bins; // histogram bin counts - void Fill(const std::vector&); + void Fill(const float* data, const size_t data_size); static bool ConsistencyCheck(const Histogram&, const Histogram&); public: Histogram() = default; // required to create empty histograms used in references - Histogram(int num_bins, float min_value, float max_value, const std::vector& data); + Histogram(int num_bins, float min_value, float max_value, const float* data, const size_t data_size); Histogram(const Histogram& h); diff --git a/src/ImageStats/StatsCalculator.cc b/src/ImageStats/StatsCalculator.cc index 7d513a944..36c4df4ea 100644 --- a/src/ImageStats/StatsCalculator.cc +++ b/src/ImageStats/StatsCalculator.cc @@ -16,19 +16,19 @@ namespace carta { -void CalcBasicStats(const std::vector& data, BasicStats& stats) { +void CalcBasicStats(BasicStats& stats, const float* data, const size_t data_size) { // Calculate stats in BasicStats struct - BasicStatsCalculator mm(data); - mm.reduce(0, data.size()); + BasicStatsCalculator mm(data, data_size); + mm.reduce(); stats = mm.GetStats(); } -Histogram CalcHistogram(int num_bins, const BasicStats& stats, const std::vector& data) { - if ((stats.min_val == std::numeric_limits::max()) || (stats.max_val == std::numeric_limits::min()) || data.empty()) { +Histogram CalcHistogram(int num_bins, const BasicStats& stats, const float* data, const size_t data_size) { + if ((stats.min_val == std::numeric_limits::max()) || (stats.max_val == std::numeric_limits::min()) || data_size == 0) { // empty / NaN region - return Histogram(1, 0, 0, {}); + return Histogram(1, 0, 0, data, data_size); } else { - return Histogram(num_bins, stats.min_val, stats.max_val, data); + return Histogram(num_bins, stats.min_val, stats.max_val, data, data_size); } } diff --git a/src/ImageStats/StatsCalculator.h b/src/ImageStats/StatsCalculator.h index 09ab28945..d9412ca64 100644 --- a/src/ImageStats/StatsCalculator.h +++ b/src/ImageStats/StatsCalculator.h @@ -19,9 +19,9 @@ namespace carta { -void CalcBasicStats(const std::vector& data, BasicStats& stats); +void CalcBasicStats(BasicStats& stats, const float* data, const size_t data_size); -Histogram CalcHistogram(int num_bins, const BasicStats& stats, const std::vector& data); +Histogram CalcHistogram(int num_bins, const BasicStats& stats, const float* data, const size_t data_size); bool CalcStatsValues(std::map>& stats_values, const std::vector& requested_stats, const casacore::ImageInterface& image, bool per_channel = true); diff --git a/src/Region/RegionHandler.cc b/src/Region/RegionHandler.cc index ef8b444b9..149e18549 100644 --- a/src/Region/RegionHandler.cc +++ b/src/Region/RegionHandler.cc @@ -994,13 +994,13 @@ bool RegionHandler::GetRegionHistogramData( // Calculate and cache stats if (!have_basic_stats) { - CalcBasicStats(data[stokes], stats); + CalcBasicStats(stats, data[stokes].data(), data[stokes].size()); _histogram_cache[cache_id].SetBasicStats(stats); have_basic_stats = true; } // Calculate and cache histogram for number of bins - Histogram histo = CalcHistogram(num_bins, stats, data[stokes]); + Histogram histo = CalcHistogram(num_bins, stats, data[stokes].data(), data[stokes].size()); _histogram_cache[cache_id].SetHistogram(num_bins, histo); // Complete Histogram submessage diff --git a/src/Util/File.cc b/src/Util/File.cc index 817a5d0c8..345c01b2b 100644 --- a/src/Util/File.cc +++ b/src/Util/File.cc @@ -38,7 +38,7 @@ int GetNumItems(const std::string& path) { try { int counter = 0; auto it = fs::directory_iterator(path); - for (const auto f : it) { + for (const auto& f : it) { counter++; } return counter; diff --git a/test/TestHistogram.cc b/test/TestHistogram.cc index 7a16f6aa4..de69df9ea 100644 --- a/test/TestHistogram.cc +++ b/test/TestHistogram.cc @@ -71,7 +71,7 @@ TEST_F(HistogramTest, TestHistogramBehaviour) { data.push_back(00.0 - 1e-9); // should not appear data.push_back(10.0 + 1e+9); // should not appear data.push_back(11.0); // should not appear - carta::Histogram hist(10, 0.0f, 10.0f, data); + carta::Histogram hist(10, 0.0f, 10.0f, data.data(), data.size()); auto counts = accumulate(hist.GetHistogramBins().begin(), hist.GetHistogramBins().end(), 0); EXPECT_EQ(counts, 12); EXPECT_EQ(hist.GetHistogramBins()[0], 2); @@ -84,7 +84,7 @@ TEST_F(HistogramTest, TestHistogramBehaviour) { data.push_back(NAN); data.push_back(NAN); data.push_back(NAN); - carta::Histogram hist2(10, 0.0f, 10.0f, data); + carta::Histogram hist2(10, 0.0f, 10.0f, data.data(), data.size()); // expect 0 counts auto bins = hist2.GetHistogramBins(); EXPECT_EQ(accumulate(bins.begin(), bins.end(), 0), 0); @@ -93,7 +93,7 @@ TEST_F(HistogramTest, TestHistogramBehaviour) { TEST_F(HistogramTest, TestHistogramConstructor) { std::vector data(1024 * 1024); std::for_each(data.begin(), data.end(), [&](float& v) { v = float_random(mt); }); - carta::Histogram hist(1024, 0.0f, 1.0f, data); + carta::Histogram hist(1024, 0.0f, 1.0f, data.data(), data.size()); carta::Histogram hist2(hist); EXPECT_TRUE(CompareResults(hist, hist2)); } @@ -101,14 +101,14 @@ TEST_F(HistogramTest, TestHistogramConstructor) { TEST_F(HistogramTest, TestHistogramAdd) { std::vector data(1024 * 1024); std::for_each(data.begin(), data.end(), [&](float& v) { v = float_random(mt); }); - carta::Histogram hist(1024, 0.0f, 1.0f, data); + carta::Histogram hist(1024, 0.0f, 1.0f, data.data(), data.size()); const auto total_counts = accumulate(hist.GetHistogramBins().begin(), hist.GetHistogramBins().end(), 0); - carta::Histogram hist2(1024, 0.0f, 1.0f, data); + carta::Histogram hist2(1024, 0.0f, 1.0f, data.data(), data.size()); EXPECT_TRUE(CompareResults(hist, hist2)); // naive? hist.Add(hist2); const auto total_counts2 = accumulate(hist.GetHistogramBins().begin(), hist.GetHistogramBins().end(), 0); EXPECT_EQ(2 * total_counts, total_counts2); - carta::Histogram hist3(512, 0.0f, 1.0f, data); + carta::Histogram hist3(512, 0.0f, 1.0f, data.data(), data.size()); EXPECT_FALSE(hist.Add(hist3)); } @@ -118,9 +118,9 @@ TEST_F(HistogramTest, TestSingleThreading) { v = float_random(mt); } carta::ThreadManager::SetThreadLimit(1); - carta::Histogram hist_st(1024, 0.0f, 1.0f, data); + carta::Histogram hist_st(1024, 0.0f, 1.0f, data.data(), data.size()); for (auto i = 2; i < 24; i++) { - carta::Histogram hist_mt(1024, 0.0f, 1.0f, data); + carta::Histogram hist_mt(1024, 0.0f, 1.0f, data.data(), data.size()); EXPECT_TRUE(CompareResults(hist_st, hist_mt)); } } @@ -132,10 +132,10 @@ TEST_F(HistogramTest, TestMultithreading) { } carta::ThreadManager::SetThreadLimit(1); - carta::Histogram hist_st(1024, 0.0f, 1.0f, data); + carta::Histogram hist_st(1024, 0.0f, 1.0f, data.data(), data.size()); for (auto i = 2; i < 24; i++) { carta::ThreadManager::SetThreadLimit(i); - carta::Histogram hist_mt(1024, 0.0f, 1.0f, data); + carta::Histogram hist_mt(1024, 0.0f, 1.0f, data.data(), data.size()); EXPECT_TRUE(CompareResults(hist_st, hist_mt)); } } @@ -151,12 +151,12 @@ TEST_F(HistogramTest, TestMultithreadingPerformance) { carta::ThreadManager::SetThreadLimit(1); t.Start("single_threaded"); - carta::Histogram hist_st(1024, 0.0f, 1.0f, data); + carta::Histogram hist_st(1024, 0.0f, 1.0f, data.data(), data.size()); t.End("single_threaded"); carta::ThreadManager::SetThreadLimit(4); t.Start("multi_threaded"); - carta::Histogram hist_mt(1024, 0.0f, 1.0f, data); + carta::Histogram hist_mt(1024, 0.0f, 1.0f, data.data(), data.size()); t.End("multi_threaded"); auto st_time = t.GetMeasurement("single_threaded");