diff --git a/_test/test_sym_op/test_cut_sqw_sym.m b/_test/test_sym_op/test_cut_sqw_sym.m
index ac842faf2c..e2c2ad030f 100644
--- a/_test/test_sym_op/test_cut_sqw_sym.m
+++ b/_test/test_sym_op/test_cut_sqw_sym.m
@@ -122,7 +122,6 @@ function test_cut_sym_reflect_half_to_whole_cut(obj)
 
         function test_cut_sym_with_pix(obj)
         % Test symmetrisation, keeping pixels
-            skipTest('Cut with pix disabled')
             clOb = set_temporary_config_options(hor_config, 'log_level', -1);
             w2sym = cut(obj.data, obj.proj, obj.bin,...
                 obj.width, obj.width, obj.ebins, obj.sym);
diff --git a/horace_core/sqw/@sqw/cut.m b/horace_core/sqw/@sqw/cut.m
index fc2d0d3134..70b160fee1 100644
--- a/horace_core/sqw/@sqw/cut.m
+++ b/horace_core/sqw/@sqw/cut.m
@@ -146,9 +146,9 @@
     [targ_ax_block, targ_proj] = obj.define_target_axes_block(targ_proj, pbin_tmp, sym);
 
     if return_cut
-        wout{cut_num} = cut_single_(obj, targ_proj, targ_ax_block, opt, log_level);
+        wout{cut_num} = cut_single_(obj, targ_proj, targ_ax_block, opt, log_level, sym);
     else
-        cut_single_(obj, targ_proj, targ_ax_block, opt, log_level);
+        cut_single_(obj, targ_proj, targ_ax_block, opt, log_level, sym);
     end
 end
 
diff --git a/horace_core/sqw/@sqw/private/cut_accumulate_data_.m b/horace_core/sqw/@sqw/private/cut_accumulate_data_.m
index 3402d27ac9..64e931b61f 100644
--- a/horace_core/sqw/@sqw/private/cut_accumulate_data_.m
+++ b/horace_core/sqw/@sqw/private/cut_accumulate_data_.m
@@ -1,5 +1,5 @@
 function [s, e, npix, pix_out, unique_runid] = ...
-    cut_accumulate_data_(obj, targ_proj, targ_axes, keep_pixels, log_level)
+    cut_accumulate_data_(obj, targ_proj, targ_axes, keep_pixels, log_level, sym)
 %%CUT_ACCUMULATE_DATA Accumulate image and pixel data for a cut
 %
 % Input:
@@ -74,7 +74,7 @@
 
 if keep_pixels
     [npix, s, e, pix_out, unique_runid] = cut_with_pixels(obj.pix, block_starts, block_sizes, targ_proj, ...
-        targ_axes, npix, s, e, log_level,keep_precision, pixel_contrib_name);
+        targ_axes, npix, s, e, log_level,keep_precision, pixel_contrib_name, sym);
 else
     [npix, s, e, pix_out, unique_runid] = cut_no_pixels(obj.pix, block_starts, block_sizes, targ_proj, ...
         targ_axes, npix, s, e, log_level, pixel_contrib_name);
@@ -166,11 +166,7 @@
 
 function [npix, s, e, pix_out, unique_runid] = cut_with_pixels(pix, block_starts, block_sizes, ...
     targ_proj, targ_axes, npix, s, e, ll, ...
-    keep_precision, pixel_contrib_name)
-
-if numel(targ_proj) > 1
-    error('HORACE:cut:not_implemented', 'Cannot cut sym and return pixels')
-end
+    keep_precision, pixel_contrib_name, sym)
 
 hc = hor_config;
 chunk_size = hc.mem_chunk_size;
@@ -225,6 +221,7 @@
         [npix, s, e, pix_ok, unique_runid_l, pix_indx, selected] = ...
             targ_proj(i).bin_pixels(targ_axes(i), candidate_pix, npix, s, e);
 
+        candidate_pix = sym{i}.transform_pix(candidate_pix, {}, selected);
         candidate_pix = candidate_pix.tag(selected);
 
         npix_step_retained = pix_ok.num_pixels; % just for logging the progress
diff --git a/horace_core/sqw/@sqw/private/cut_single_.m b/horace_core/sqw/@sqw/private/cut_single_.m
index 02c1df6d6d..c49590fe5a 100644
--- a/horace_core/sqw/@sqw/private/cut_single_.m
+++ b/horace_core/sqw/@sqw/private/cut_single_.m
@@ -1,4 +1,4 @@
-function wout = cut_single_(w, targ_proj, targ_axes, opt, log_level)
+function wout = cut_single_(w, targ_proj, targ_axes, opt, log_level, sym)
 %%CUT_SINGLE Perform a cut on a single sqw object
 %
 % Input:
@@ -27,7 +27,7 @@
 
 % Accumulate image and pixel data for cut
 [s, e, npix, pix_out,runid_contributed] = cut_accumulate_data_( ...
-    w, targ_proj, targ_axes, opt.keep_pix, log_level);
+    w, targ_proj, targ_axes, opt.keep_pix, log_level, sym);
 
 
 if isa(pix_out, 'MultipixBase')
@@ -51,7 +51,7 @@
 if opt.keep_pix
     wout = sqw();
     wout.main_header = w.main_header;
-    % NB detpar is no longer copied as detpar just exposes the detector_arrays 
+    % NB detpar is no longer copied as detpar just exposes the detector_arrays
     % already in experiment_info
     wout.data   = data_out.data;
     wout.pix  = data_out.pix;
diff --git a/horace_core/symop/@Symop/Symop.m b/horace_core/symop/@Symop/Symop.m
index 66f4c9d826..bd692a7d39 100644
--- a/horace_core/symop/@Symop/Symop.m
+++ b/horace_core/symop/@Symop/Symop.m
@@ -115,7 +115,7 @@ function disp(obj)
             end
         end
 
-        function pix = transform_pix(obj, pix, proj)
+        function pix = transform_pix(obj, pix, proj, selected)
             % Transform pixel coordinates into symmetry related coordinates
             %
             % The transformation converts the components of a vector which is
@@ -138,6 +138,7 @@ function disp(obj)
             %
             %   pix         PixelData object
             %
+            %   selected    Pixels to transform
             % Output:
             % -------
             %   pix         Transformed PixelData object
@@ -145,6 +146,9 @@ function disp(obj)
             if ~exist('proj', 'var')
                 proj = {};
             end
+            if ~exist('selected', 'var')
+                selected = 1:pix.num_pixels;
+            end
 
             % Check input
             if ~isa(pix, 'PixelDataBase')
@@ -156,6 +160,7 @@ function disp(obj)
             if isa(pix, 'PixelDataMemory')
                 for i = numel(obj):-1:1
                     sel = obj(i).in_irreducible(pix.q_coordinates, proj{:});
+                    sel(~selected) = false;
                     pix.q_coordinates(:, ~sel) = obj(i).transform_vec(pix.q_coordinates(:, ~sel));
                 end
             else