Use Z axis rather than spectral axis to offset slicer before accessing channel cache

CHANGELOG.md

@@ -18,6 +18,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 * Fix hdf5 image distortion after animation stops ([#1368](https://github.com/CARTAvis/carta-backend/issues/1368)).
 * Fix matched polygon region approximation crash ([#1383](https://github.com/CARTAvis/carta-backend/issues/1383)).
 * Fix save image/export regions bug which could cause directory overwrite or deletion ([#1377](https://github.com/CARTAvis/carta-backend/issues/1377)).
+* Fix bug in cache slicer transformation which affects some images with rotated axes ([#1389](https://github.com/CARTAvis/carta-backend/pull/1389)).
 
 ### Changed
 * Move the loader cache to separate files ([#1021](https://github.com/CARTAvis/carta-backend/issues/1021)).

src/Frame/Frame.cc

@@ -193,6 +193,18 @@ int Frame::StokesAxis() {
     return _stokes_axis;
 }
 
+int Frame::XAxis() {
+    return _x_axis;
+}
+
+int Frame::YAxis() {
+    return _y_axis;
+}
+
+int Frame::ZAxis() {
+    return _z_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());
@@ -1694,10 +1706,10 @@ bool Frame::GetSlicerData(const StokesSlicer& stokes_slicer, float* data) {
         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 (_z_axis >= 0) {
+            cache_shape(_z_axis) = 1;
+            slicer_start(_z_axis) = 0;
+            slicer_end(_z_axis) = 0;
         }
         if (_stokes_axis >= 0) {
             cache_shape(_stokes_axis) = 1;

src/Frame/Frame.h

@@ -102,6 +102,9 @@ class Frame {
     bool IsCurrentZStokes(const StokesSource& stokes_source);
     int SpectralAxis();
     int StokesAxis();
+    int XAxis();
+    int YAxis();
+    int ZAxis();
     bool GetBeams(std::vector<CARTA::Beam>& beams);
 
     // Slicer to set z and stokes ranges with full xy plane

src/ImageData/FileLoader.cc

@@ -232,6 +232,8 @@ bool FileLoader::FindCoordinateAxes(casacore::IPosition& shape, std::vector<int>
         }
     }
 
+    _x_axis = render_axes[0];
+    _y_axis = render_axes[1];
     _width = shape(render_axes[0]);
     _height = shape(render_axes[1]);
     _image_plane_size = _width * _height;
@@ -938,8 +940,8 @@ typename FileLoader::ImageRef FileLoader::GetStokesImage(const StokesSource& sto
 
     if (_stokes_source != stokes_source) {
         // compute new stokes image with respect to the channel range
-        carta::PolarizationCalculator polarization_calculator(
-            GetImage(), AxisRange(stokes_source.z_range), AxisRange(stokes_source.x_range), AxisRange(stokes_source.y_range));
+        carta::PolarizationCalculator polarization_calculator(GetImage(), {_x_axis, _y_axis, _z_axis, _stokes_axis},
+            AxisRange(stokes_source.z_range), AxisRange(stokes_source.x_range), AxisRange(stokes_source.y_range));
 
         if (stokes_source.stokes == CARTA::PolarizationType::Ptotal) {
             _computed_stokes_image = polarization_calculator.ComputeTotalPolarizedIntensity();

src/ImageData/FileLoader.h

@@ -169,7 +169,7 @@ class FileLoader {
     casacore::IPosition _image_shape;
     size_t _num_dims, _image_plane_size;
     size_t _width, _height, _depth, _num_stokes;
-    int _z_axis, _stokes_axis;
+    int _x_axis, _y_axis, _z_axis, _stokes_axis;
     std::shared_ptr<casacore::CoordinateSystem> _coord_sys;
     bool _has_pixel_mask;
     casacore::DataType _data_type;

src/ImageData/PolarizationCalculator.cc

@@ -10,7 +10,7 @@
 using namespace carta;
 
 PolarizationCalculator::PolarizationCalculator(
-    std::shared_ptr<casacore::ImageInterface<float>> image, AxisRange z_range, AxisRange x_range, AxisRange y_range)
+    std::shared_ptr<casacore::ImageInterface<float>> image, std::vector<int> axes, AxisRange z_range, AxisRange x_range, AxisRange y_range)
     : _image(image), _image_valid(true) {
     const auto ndim = _image->ndim();
     if (ndim < 4) {
@@ -20,56 +20,45 @@ PolarizationCalculator::PolarizationCalculator(
     }
 
     const auto& coord_sys = _image->coordinates();
-    std::vector<int> dir_axes = {0, 1}; // By default, the spatial axes numbers are 0 and 1
-    if (coord_sys.hasDirectionCoordinate()) {
-        casacore::Vector<casacore::Int> tmp_axes = coord_sys.directionAxesNumbers();
-        dir_axes[0] = tmp_axes[0];
-        dir_axes[1] = tmp_axes[1];
-    }
-
-    auto spectral_axis = coord_sys.spectralAxisNumber();
-    if (spectral_axis < 0) {
-        spectral_axis = 2; // Assume spectral axis number = 2
-    }
 
-    auto stokes_axis = coord_sys.polarizationAxisNumber();
-    if (stokes_axis < 0) {
-        stokes_axis = 3; // Assume stokes axis number = 3
-    }
+    int x_axis = axes[0];
+    int y_axis = axes[1];
+    int z_axis = axes[2];
+    int stokes_axis = axes[3];
 
     const auto shape = _image->shape();
     casacore::IPosition blc(ndim, 0);
     casacore::IPosition trc = shape - 1;
 
     if (x_range.to == ALL_X) {
         x_range.from = 0;
-        x_range.to = shape(dir_axes[0]) - 1;
+        x_range.to = shape(x_axis) - 1;
     }
 
     if (y_range.to == ALL_Y) {
         y_range.from = 0;
-        y_range.to = shape(dir_axes[1]) - 1;
+        y_range.to = shape(y_axis) - 1;
     }
 
     if (z_range.to == ALL_Z) {
         z_range.from = 0;
-        z_range.to = shape(spectral_axis) - 1;
+        z_range.to = shape(z_axis) - 1;
     }
 
-    if (x_range.from < 0 || x_range.to >= shape(dir_axes[0]) || y_range.from < 0 || y_range.to >= shape(dir_axes[1]) || z_range.from < 0 ||
-        z_range.to >= shape(spectral_axis)) {
+    if (x_range.from < 0 || x_range.to >= shape(x_axis) || y_range.from < 0 || y_range.to >= shape(y_axis) || z_range.from < 0 ||
+        z_range.to >= shape(z_axis)) {
         spdlog::error("Invalid selection region.");
         _image_valid = false;
         return;
     }
 
     // Make a region
-    blc(dir_axes[0]) = x_range.from;
-    trc(dir_axes[0]) = x_range.to;
-    blc(dir_axes[1]) = y_range.from;
-    trc(dir_axes[1]) = y_range.to;
-    blc(spectral_axis) = z_range.from;
-    trc(spectral_axis) = z_range.to;
+    blc(x_axis) = x_range.from;
+    trc(x_axis) = x_range.to;
+    blc(y_axis) = y_range.from;
+    trc(y_axis) = y_range.to;
+    blc(z_axis) = z_range.from;
+    trc(z_axis) = z_range.to;
 
     // Get stokes indices and make stokes regions
     if (coord_sys.hasPolarizationCoordinate()) {

src/ImageData/PolarizationCalculator.h

@@ -23,8 +23,8 @@ class PolarizationCalculator {
     enum StokesTypes { I, Q, U, V };
 
 public:
-    PolarizationCalculator(std::shared_ptr<casacore::ImageInterface<float>> image, AxisRange z_range = AxisRange(ALL_Z),
-        AxisRange x_range = AxisRange(ALL_X), AxisRange y_range = AxisRange(ALL_Y));
+    PolarizationCalculator(std::shared_ptr<casacore::ImageInterface<float>> image, std::vector<int> axes, AxisRange z_range,
+        AxisRange x_range, AxisRange y_range);
     ~PolarizationCalculator() = default;
 
     std::shared_ptr<casacore::ImageInterface<float>> ComputeTotalPolarizedIntensity();

src/Region/RegionHandler.cc

@@ -913,11 +913,10 @@ bool RegionHandler::ApplyRegionToFile(int region_id, int file_id, const AxisRang
         } else {
             // Extension extends applied_region in xy axes by z/stokes axes only
             // Remove xy axes from z/stokes box
-            casacore::IPosition remove_xy(2, 0, 1);
+            casacore::IPosition remove_xy(2, _frames.at(file_id)->XAxis(), _frames.at(file_id)->YAxis());
             z_stokes_slicer =
                 casacore::Slicer(z_stokes_slicer.start().removeAxes(remove_xy), z_stokes_slicer.length().removeAxes(remove_xy));
             casacore::LCBox z_stokes_box(z_stokes_slicer, image_shape.removeAxes(remove_xy));
-
             casacore::IPosition extend_axes = casacore::IPosition::makeAxisPath(image_shape.size()).removeAxes(remove_xy);
             casacore::LCExtension final_region(*applied_region, extend_axes, z_stokes_box);
             stokes_region.image_region = casacore::ImageRegion(final_region);