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 <[email protected]>
CARTAvis · Apr 9, 2024 · f17756b · f17756b
1 parent de4da52
commit f17756b
Show file tree

Hide file tree

Showing 5 changed files with 141 additions and 65 deletions.
diff --git 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
@@ -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<CARTA::Beam>& 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<float>& data, bool report_performance) {
     // Get image data with a region applied
     Timer t;
-    casacore::SubImage<float> sub_image;
-    bool subimage_ok = GetRegionSubImage(stokes_region, sub_image);
+    std::vector<bool> 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<bool>)
+                casacore::Array<bool> 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<float> sub_image;
+        bool subimage_ok = GetRegionSubImage(stokes_region, sub_image);
 
-        // Get image data
-        std::unique_lock<std::mutex> ulock(_image_mutex);
-        if (_loader->IsGenerated() || is_computed_stokes) { // For the image in memory
-            casacore::Array<float> tmp;
-            sub_image.doGetSlice(tmp, slicer);
-            data = tmp.tovector();
-        } else {
-            data.resize(subimage_shape.product()); // must size correctly before sharing
-            casacore::Array<float> 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<bool> 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<bool> 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<std::mutex> ulock(_image_mutex);
+            casacore::Array<float> 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<float>(subimage_shape, data.data(), casacore::StorageInitPolicy::SHARE);
+                sub_image.doGetSlice(tmpdata, slicer);
             }
+
+            // Get mask that defines region in subimage bounding box
+            casacore::Array<bool> 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<float> tmp(stokes_slicer.slicer.length(), data, casacore::StorageInitPolicy::SHARE);
-    std::unique_lock<std::mutex> 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<float> 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<std::mutex> 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
@@ -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<CARTA::Beam>& beams);
@@ -296,7 +297,6 @@ class Frame {
     ContourSettings _contour_settings;
 
     // Image data cache and mutex
-    //    std::vector<float> _image_cache; // image data for current z, stokes
     long long int _image_cache_size;
     std::unique_ptr<float[]> _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
@@ -12,7 +12,6 @@
 #include <casacore/lattices/LRegions/RegionType.h>
 
 #include "RegionConverter.h"
-// #include "Util/Image.h"
 
 using namespace carta;
 

diff --git 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<float> RegionHandler::GetTemporaryRegionProfile(int region_idx,
     }
 
     std::lock_guard<std::mutex> 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<float> 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<casacore::LCRegion> 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<float> region_data;
             if (_frames.at(file_id)->GetRegionData(stokes_region, region_data, false)) {
-                // Get BasicStats
-                BasicStats<float> 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();