Migrate 2D image cache away from std::vector<float> (#967) (#1034)

* replace vector with pointers

* fix style

* remove additional debug printing

* fix compilation

* add range check to BasicStatsCalculator

* simplify BasicStatsCalculator

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<float[]>(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<float>& 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<float>& 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<float[]> 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,19 +832,19 @@ bool Frame::GetBasicStats(int z, int stokes, BasicStats<float>& 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;
         }
 
         // calculate histogram from given z/stokes data
         std::vector<float> 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<float> 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<CARTA::SetSpatialR
             cursor_value = cursor_value_with_current_stokes;
         } else {
             casacore::Slicer section = GetImageSlicer(AxisRange(x), AxisRange(y), AxisRange(CurrentZ()), stokes);
-            std::vector<float> data;
-            if (GetSlicerData(section, data)) {
+            const auto N = section.length().product();
+            std::shared_ptr<float[]> 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<CARTA::SetSpatialR
                         section = GetImageSlicer(AxisRange(x), AxisRange(start, end - 1), AxisRange(CurrentZ()), stokes);
                     }
 
-                    have_profile = GetSlicerData(section, profile);
+                    profile.resize(section.length().product());
+                    have_profile = GetSlicerData(section, profile.data());
                 }
             }
 
@@ -1483,13 +1488,14 @@ bool Frame::FillSpectralProfileData(std::function<void(CARTA::SpectralProfileDat
                     size_t nz = (start(_z_axis) + delta_z < profile_size ? delta_z : profile_size - start(_z_axis));
                     count(_z_axis) = nz;
                     casacore::Slicer slicer(start, count);
-                    std::vector<float> buffer;
-                    if (!GetSlicerData(slicer, buffer)) {
+                    const auto N = slicer.length().product();
+                    std::shared_ptr<float[]> 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<float>& 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<float> tmp(slicer.length(), data.data(), casacore::StorageInitPolicy::SHARE);
+bool Frame::GetSlicerData(const casacore::Slicer& slicer, float* data) {
+    casacore::Array<float> tmp(slicer.length(), data, casacore::StorageInitPolicy::SHARE);
     std::unique_lock<std::mutex> ulock(_image_mutex);
     bool data_ok = _loader->GetSlice(tmp, slicer);
     _loader->CloseImageIfUpdated();

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<float>& data);
-    bool GetSlicerData(const casacore::Slicer& slicer, std::vector<float>& 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<CARTA::StatsType>& required_stats, bool per_z,
         std::map<CARTA::StatsType, std::vector<double>>& stats_values);
@@ -277,12 +277,14 @@ class Frame {
     ContourSettings _contour_settings;
 
     // Image data cache and mutex
-    std::vector<float> _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<float> _image_cache; // image data for current z, stokes
+    long long int _image_cache_size;
+    std::shared_ptr<float[]> _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;
 

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<T>& _data;
+    const T* _data;
+    size_t _data_size;
 
 public:
-    BasicStatsCalculator(const std::vector<T>& 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<T> GetStats() const;
 };

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 <cmath>
 
 namespace carta {
@@ -41,19 +43,20 @@ BasicStats<T>::BasicStats()
       sumSq(0) {}
 
 template <typename T>
-BasicStatsCalculator<T>::BasicStatsCalculator(const std::vector<T>& data)
+BasicStatsCalculator<T>::BasicStatsCalculator(const T* data, size_t data_size)
     : _min_val(std::numeric_limits<T>::max()),
       _max_val(std::numeric_limits<T>::lowest()),
       _sum(0),
       _sum_squares(0),
       _num_pixels(0),
-      _data(data) {}
+      _data(data),
+      _data_size(data_size) {}
 
 template <typename T>
-void BasicStatsCalculator<T>::reduce(const size_t start, const size_t end) {
+void BasicStatsCalculator<T>::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) {

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<float>& 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<float>& data) {
+void Histogram::Fill(const float* data, const size_t data_size) {
     std::vector<int64_t> 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

src/ImageStats/Histogram.h

@@ -19,12 +19,12 @@ class Histogram {
     float _bin_center;                // bin center
     std::vector<int> _histogram_bins; // histogram bin counts
 
-    void Fill(const std::vector<float>&);
+    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<float>& data);
+    Histogram(int num_bins, float min_value, float max_value, const float* data, const size_t data_size);
 
     Histogram(const Histogram& h);
 

src/ImageStats/StatsCalculator.cc

@@ -16,19 +16,19 @@
 
 namespace carta {
 
-void CalcBasicStats(const std::vector<float>& data, BasicStats<float>& stats) {
+void CalcBasicStats(BasicStats<float>& stats, const float* data, const size_t data_size) {
     // Calculate stats in BasicStats struct
-    BasicStatsCalculator<float> mm(data);
-    mm.reduce(0, data.size());
+    BasicStatsCalculator<float> mm(data, data_size);
+    mm.reduce();
     stats = mm.GetStats();
 }
 
-Histogram CalcHistogram(int num_bins, const BasicStats<float>& stats, const std::vector<float>& data) {
-    if ((stats.min_val == std::numeric_limits<float>::max()) || (stats.max_val == std::numeric_limits<float>::min()) || data.empty()) {
+Histogram CalcHistogram(int num_bins, const BasicStats<float>& stats, const float* data, const size_t data_size) {
+    if ((stats.min_val == std::numeric_limits<float>::max()) || (stats.max_val == std::numeric_limits<float>::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);
     }
 }
 

src/ImageStats/StatsCalculator.h

@@ -19,9 +19,9 @@
 
 namespace carta {
 
-void CalcBasicStats(const std::vector<float>& data, BasicStats<float>& stats);
+void CalcBasicStats(BasicStats<float>& stats, const float* data, const size_t data_size);
 
-Histogram CalcHistogram(int num_bins, const BasicStats<float>& stats, const std::vector<float>& data);
+Histogram CalcHistogram(int num_bins, const BasicStats<float>& stats, const float* data, const size_t data_size);
 
 bool CalcStatsValues(std::map<CARTA::StatsType, std::vector<double>>& stats_values, const std::vector<CARTA::StatsType>& requested_stats,
     const casacore::ImageInterface<float>& image, bool per_channel = true);

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

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;