diff --git a/examples/Python/PET/osem_gated_w_motion.py b/examples/Python/PET/osem_gated_w_motion.py
new file mode 100644
index 000000000..6205c1afc
--- /dev/null
+++ b/examples/Python/PET/osem_gated_w_motion.py
@@ -0,0 +1,384 @@
+"""
+OSEM reconstruction demo for gated data with motion.
+
+We actually use the OSMAPOSL reconstructor in this demo. This reconstructor
+implements an Ordered Subsets (OS) version of the One Step Late algorithm (OSL)
+from Green et al for Maximum a Posteriori (MAP) maximisation. Here we use it
+for Maximum Likelihood (ML) in which case it is equivalent to OSEM.
+
+Usage:
+  osem_reconstruction_gated_w_motion [--help | options]
+
+Options:
+  -S <file>, --sino=<file>          sinogram pattern, * or % wildcard
+                                    (e.g., sino_ms*.hs). Enclose in quotations.
+  -a <pattern>, --attn=<pattern>    attenuation pattern, * or % wildcard
+                                    (e.g., attn_ms*.hv). Enclose in quotations.
+  -R <pattern>, --rand=<pattern>    randoms pattern, * or % wildcard
+                                    (e.g., rand_ms*.hs). Enclose in quotations.
+  -n <norm>, --norm=<norm>          ECAT8 bin normalization file
+  -T <pattern>, --trans=<pattern>   transformation pattern, * or % wildcard
+                                    (e.g., tm_ms*.txt). Enclose in quotations.
+  -t <str>, --trans_type=<str>      transformation type (tm, disp, def)
+                                    [default: tm]
+  -s <subs>, --subs=<subs>          number of subsets [default: 21]
+  -i <siter>, --subiter=<siter>     number of sub-iterations [default: 2]
+  -e <engn>, --engine=<engn>        reconstruction engine [default: STIR]
+  -o <outp>, --outp=<outp>          output file prefix [default: recon]
+  --initial=<str>                   initial estimate
+  --nxny=<nxny>                     image x and y dimension [default: 127]
+  --dxdy=<dxdy>                     image x and y spacing
+                                    (default: determined by scanner)
+  --b_spline_order=<val>            b-spline order for resampler [default: 1]
+  --verbosity=<int>                 Verbosity [default: 0]
+  --nifti                           save output as nifti
+  --gpu                             use gpu
+"""
+
+# SyneRBI Synergistic Image Reconstruction Framework (SIRF)
+# Copyright 2015 - 2018 Rutherford Appleton Laboratory STFC
+# Copyright 2015 - 2020 University College London.
+#
+# This is software developed for the Collaborative Computational
+# Project in Synergistic Reconstruction for Biomedical Imaging
+# (formerly CCP PETMR)
+# (http://www.ccpsynerbi.ac.uk/).
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#       http://www.apache.org/licenses/LICENSE-2.0
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+import os
+from glob import glob
+from docopt import docopt
+import sirf.Reg as reg
+
+
+__version__ = '0.1.0'
+args = docopt(__doc__, version=__version__)
+
+# import engine module
+if args['--engine'] == 'STIR':
+    import sirf.STIR as pet
+
+pet.AcquisitionData.set_storage_scheme('memory')
+
+# Multiple files
+sino_pattern = str(args['--sino']).replace('%', '*')
+attn_pattern = str(args['--attn']).replace('%', '*')
+rand_pattern = str(args['--rand']).replace('%', '*')
+trans_pattern = str(args['--trans']).replace('%', '*')
+trans_type = str(args['--trans_type'])
+
+if attn_pattern is None:
+    attn_pattern = ""
+if rand_pattern is None:
+    rand_pattern = ""
+
+# Norm
+norm_file = None
+if args['--norm']:
+    norm_file = str(args['--norm'])
+    if not os.path.isfile(norm_file):
+        raise pet.error("Norm file not found: " + norm_file)
+
+# Initial estimate
+initial_estimate = None
+if args['--initial']:
+    initial_estimate = str(args['--initial'])
+
+# Verbosity
+verbosity = int(args['--verbosity'])
+pet.set_verbosity(verbosity)
+if verbosity == 0:
+    msg_red = pet.MessageRedirector(None, None, None)
+
+# Rest of the input options
+num_subsets = int(args['--subs'])
+num_subiterations = int(args['--subiter'])
+nxny = int(args['--nxny'])
+outp_file = args['--outp']
+b_spline_order = args['--b_spline_order']
+nifti = True if args['--nifti'] else False
+use_gpu = True if args['--gpu'] else False
+
+
+def get_resampler(image, ref=None, trans=None):
+    """Returns a NiftyResample object for the specified transform and image."""
+    if ref is None:
+        ref = image
+    resampler = reg.NiftyResample()
+    resampler.set_reference_image(ref)
+    resampler.set_floating_image(image)
+    resampler.set_padding_value(0)
+    resampler.set_interpolation_type_to_linear()
+    if trans is not None:
+        resampler.add_transformation(trans)
+    return resampler
+
+
+def get_asm_attn(sino, attn, acq_model):
+    """Get attn ASM from sino, attn image and acq model."""
+    asm_attn = pet.AcquisitionSensitivityModel(attn, acq_model)
+    # temporary fix pending attenuation offset fix in STIR:
+    # converting attenuation into 'bin efficiency'
+    asm_attn.set_up(sino)
+    bin_eff = pet.AcquisitionData(sino)
+    bin_eff.fill(1.0)
+    asm_attn.unnormalise(bin_eff)
+    asm_attn = pet.AcquisitionSensitivityModel(bin_eff)
+    return asm_attn
+
+
+def parse_input_files():
+    """Parse input files."""
+    if trans_pattern is None:
+        raise AssertionError("--trans missing")
+    if sino_pattern is None:
+        raise AssertionError("--sino missing")
+    trans_files = sorted(glob(trans_pattern))
+    sino_files = sorted(glob(sino_pattern))
+    attn_files = sorted(glob(attn_pattern))
+    rand_files = sorted(glob(rand_pattern))
+
+    num_ms = len(sino_files)
+    # Check some sinograms found
+    if num_ms == 0:
+        raise AssertionError("No sinograms found!")
+    # Should have as many trans as sinos
+    if num_ms != len(trans_files):
+        raise AssertionError("#trans should match #sinos. "
+                             "#sinos = " + str(num_ms) +
+                             ", #trans = " + str(len(trans_files)))
+    # If any rand, check num == num_ms
+    if len(rand_files) > 0 and len(rand_files) != num_ms:
+        raise AssertionError("#rand should match #sinos. "
+                             "#sinos = " + str(num_ms) +
+                             ", #rand = " + str(len(rand_files)))
+
+    # For attn, there should be 0, 1 or num_ms images
+    if len(attn_files) > 1 and len(attn_files) != num_ms:
+        raise AssertionError("#attn should be 0, 1 or #sinos")
+
+    return num_ms, trans_files, sino_files, attn_files, rand_files
+
+
+def read_trans(num_ms, trans_files):
+    """Read transformations."""
+    trans = num_ms*[None]
+    for i in range(num_ms):
+        if trans_type == "disp":
+            trans[i] = reg.NiftiImageData3DDisplacement(trans_files[i])
+        elif trans_type == "tm":
+            trans[i] = reg.NiftiImageData3DDisplacement(
+                reg.AffineTransformation(trans_files[i]))
+        elif trans_type == "def":
+            trans[i] = reg.NiftiImageData3DDisplacement(
+                reg.NiftiImageData3DDeformation(trans_files[i]))
+        else:
+            raise pet.error("Unknown transformation type")
+    return trans
+
+
+def get_initial_estimate(sinos):
+    """Get initial estimate."""
+    if initial_estimate:
+        image = pet.ImageData(initial_estimate)
+    else:
+        # Create image based on ProjData
+        image = sinos[0].create_uniform_image(0.0, (nxny, nxny))
+        # If using GPU, need to make sure that image is right size.
+        if use_gpu:
+            dim = (127, 320, 320)
+            spacing = (2.03125, 2.08626, 2.08626)
+        # elif non-default spacing desired
+        elif args['--dxdy']:
+            dim = image.dimensions()
+            dxdy = float(args['--dxdy'])
+            spacing = (image.voxel_sizes()[0], dxdy, dxdy)
+        if use_gpu or args['--dxdy']:
+            image.initialise(dim=dim,
+                             vsize=spacing)
+    return image
+
+
+def get_attn_images(num_ms, attns, trans, image):
+    """
+    Get attn images.
+
+    If none supplied, return None.
+    If 1 supplied, resample to each motion state.
+    If num_ms supplied return them.
+    If GPU projector, resample attn images for correct dimensions.
+    """
+    resampled_attns = None
+    if len(attns) > 0:
+        resampled_attns = [0]*num_ms
+        # if using GPU, dimensions of attn and recon images have to match
+        ref = image if use_gpu else None
+        for i in range(num_ms):
+            # if we only have 1 attn image, then we need to resample into
+            # space of each gate. However, if we have num_ms attn images, then
+            # assume they are already in the correct position, so use None as
+            # transformation.
+            tran = trans[i] if len(attns) == 1 else None
+            # If only 1 attn image, then resample that. If we have num_ms attn
+            # images, then use each attn image of each frame.
+            attn = attns[0] if len(attns) == 1 else attns[i]
+            resam = get_resampler(attn, ref=ref, trans=tran)
+            resampled_attns[i] = resam.forward(attn)
+    return resampled_attns
+
+
+def get_acq_models(num_ms, resampled_attns, sinos, rands, image):
+    """Get acquisition models."""
+    print("Setting up acquisition models...")
+    if not use_gpu:
+        acq_models = num_ms * [pet.AcquisitionModelUsingRayTracingMatrix()]
+    else:
+        acq_models = num_ms * [pet.AcquisitionModelUsingNiftyPET()]
+        for acq_model in acq_models:
+            acq_model.set_use_truncation(True)
+            acq_model.set_cuda_verbosity(verbosity)
+
+    # If present, create ASM from ECAT8 normalisation data
+    asm_norm = None
+    if norm_file:
+        asm_norm = pet.AcquisitionSensitivityModel(norm_file)
+
+    # Loop over each motion state
+    for ind in range(num_ms):
+        # Create attn ASM if necessary
+        asm_attn = None
+        if resampled_attns:
+            asm_attn = get_asm_attn(sinos[ind], resampled_attns[ind],
+                                    acq_models[ind])
+
+        # Get ASM dependent on attn and/or norm
+        asm = None
+        if asm_norm and asm_attn:
+            if ind == 0:
+                print("ASM contains norm and attenuation...")
+            asm = pet.AcquisitionSensitivityModel(asm_norm, asm_attn)
+        elif asm_norm:
+            if ind == 0:
+                print("ASM contains norm...")
+            asm = asm_norm
+        elif asm_attn:
+            if ind == 0:
+                print("ASM contains attenuation...")
+            asm = asm_attn
+        if asm:
+            acq_models[ind].set_acquisition_sensitivity(asm)
+
+        if len(rands) > 0:
+            acq_models[ind].set_background_term(rands[ind])
+
+        # Set up
+        acq_models[ind].set_up(sinos[ind], image)
+        return acq_models
+
+
+def set_up_reconstructor(num_ms, sinos, acq_models, trans, image):
+    """Set up reconstructor."""
+    print("Setting up reconstructor...")
+
+    # define objective function to be maximized as
+    # Poisson logarithmic likelihood (with linear model for mean)
+    obj_fun = pet.PoissonLogLhLinModMeanGatedProjDataWMotion()
+
+    # Loop over each motion state and set
+    for i in range(num_ms):
+        obj_fun.add_gate(sinos[i], acq_models[i], trans[i])
+
+    # select Ordered Subsets Maximum A-Posteriori One Step Late as the
+    # reconstruction algorithm (since we are not using a penalty, or prior, in
+    # this example, we actually run OSEM);
+    # this algorithm does not converge to the maximum of the objective function
+    # but is used in practice to speed-up calculations
+    recon = pet.OSMAPOSLReconstructor()
+    recon.set_objective_function(obj_fun)
+    recon.set_num_subsets(num_subsets)
+    recon.set_num_subiterations(num_subiterations)
+
+    # set up the reconstructor based on a sample image
+    # (checks the validity of parameters, sets up objective function
+    # and other objects involved in the reconstruction, which involves
+    # computing/reading sensitivity image etc etc.)
+    print('setting up, please wait...')
+    recon.set_up(image)
+
+    # set the initial image estimate
+    recon.set_current_estimate(image)
+    return recon
+
+
+def main():
+
+    ###########################################################################
+    # Parse input files
+    ###########################################################################
+
+    num_ms, trans_files, sino_files, attn_files, rand_files = \
+        parse_input_files()
+
+    ###########################################################################
+    # Read input
+    ###########################################################################
+
+    trans = read_trans(num_ms, trans_files)
+    sinos = [pet.AcquisitionData(file) for file in sino_files]
+    attns = [pet.ImageData(file) for file in attn_files]
+    rands = [pet.AcquisitionData(file) for file in rand_files]
+
+    ###########################################################################
+    # Initialise recon image
+    ###########################################################################
+
+    image = get_initial_estimate(sinos)
+
+    ###########################################################################
+    # Resample attenuation images (if necessary)
+    ###########################################################################
+
+    resampled_attns = get_attn_images(num_ms, attns, trans, image)
+
+    ###########################################################################
+    # Set up acquisition models
+    ###########################################################################
+
+    acq_models = get_acq_models(num_ms, resampled_attns, sinos, rands, image)
+
+    ###########################################################################
+    # Set up reconstructor
+    ###########################################################################
+
+    recon = set_up_reconstructor(num_ms, sinos, acq_models, trans, image)
+
+    ###########################################################################
+    # Reconstruct
+    ###########################################################################
+
+    print('reconstructing, please wait...')
+    recon.process()
+    out = recon.get_output()
+    if not args['--nifti']:
+        out.write(outp_file)
+    else:
+        reg.NiftiImageData(out).write(outp_file)
+
+
+# if anything goes wrong, an exception will be thrown
+# (cf. Error Handling section in the spec)
+try:
+    main()
+    print('done')
+except pet.error as err:
+    # display error information
+    print('%s' % err.value)
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index bb162df38..ca5a91097 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -79,6 +79,7 @@ if (DISABLE_STIR)
   message(STATUS "STIR support disabled.")
 else()
   find_package(STIR 4.0.0 REQUIRED)
+  ## NiftyPET
   if (STIR_WITH_NiftyPET_PROJECTOR)
     set(NiftyPET_BOOL_STR "1")
   	message(STATUS "STIR was built with NiftyPET, GPU projectors will be enabled.")
@@ -91,6 +92,22 @@ else()
     MESSAGE(STATUS "STIR not built with NiftyPET.")
     set(NiftyPET_BOOL_STR "0")
   endif()
+
+  ## Gated recon w/ motion
+  if (STIR_WITH_EXPERIMENTAL AND STIR_BUILT_WITH_ITK)
+    set(STIR_GATED_MOTION "1")
+    message(STATUS "STIR built with EXPERIMENTAL and ITK, gated motion recon enabled.")
+    if(${CMAKE_VERSION} VERSION_LESS "3.12.0")
+      add_definitions(-DSTIR_GATED_MOTION)
+    else()
+      add_compile_definitions(STIR_GATED_MOTION)
+    endif()
+  else()
+    MESSAGE(STATUS "STIR not built with EXPERIMENTAL or ITK, gated motion recon disabled.")
+    set(STIR_GATED_MOTION "0")
+    remove_definitions(-DSTIR_GATED_MOTION)
+  endif()
+
   ADD_SUBDIRECTORY(xSTIR)
 endif()
 
diff --git a/src/Synergistic/tests/CMakeLists.txt b/src/Synergistic/tests/CMakeLists.txt
index a9fb52aa6..749e59e5d 100644
--- a/src/Synergistic/tests/CMakeLists.txt
+++ b/src/Synergistic/tests/CMakeLists.txt
@@ -56,6 +56,13 @@ ADD_TEST(NAME SYN_TEST_CPLUSPLUS
 	COMMAND SYN_TEST_CPLUSPLUS ${test_data}
 	WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
 
+# Poisson likelihood with gated proj data and motion
+if (STIR_WITH_EXPERIMENTAL AND STIR_BUILT_WITH_ITK)
+	ADD_EXECUTABLE (test_PoissonGatedMotion test_PoissonGatedMotion.cpp ${STIR_REGISTRIES})
+	TARGET_LINK_LIBRARIES(test_PoissonGatedMotion PUBLIC cstir Reg)
+	ADD_TEST(NAME SYN_TEST_CPLUSPLUS_POISSONWMOTION COMMAND test_PoissonGatedMotion WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR})
+endif()
+
 # If python
 if (BUILD_PYTHON)
 	# Make into test
diff --git a/src/Synergistic/tests/test_PoissonGatedMotion.cpp b/src/Synergistic/tests/test_PoissonGatedMotion.cpp
new file mode 100644
index 000000000..c9f5ac172
--- /dev/null
+++ b/src/Synergistic/tests/test_PoissonGatedMotion.cpp
@@ -0,0 +1,113 @@
+/*
+SyneRBI Synergistic Image Reconstruction Framework (SIRF)
+Copyright 2020 University College London
+
+This is software developed for the Collaborative Computational
+Project in Synergistic Reconstruction for Biomedical Imaging (formerly CCP PETMR)
+(http://www.ccpsynerbi.ac.uk/).
+
+Licensed under the Apache License, Version 2.0 (the "License");
+you may not use this file except in compliance with the License.
+You may obtain a copy of the License at
+http://www.apache.org/licenses/LICENSE-2.0
+Unless required by applicable law or agreed to in writing, software
+distributed under the License is distributed on an "AS IS" BASIS,
+WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+See the License for the specific language governing permissions and
+limitations under the License.
+
+*/
+
+/*!
+\file
+\ingroup STIR Tests
+
+\author Richard Brown
+\author SyneRBI
+*/
+
+// STIR stuff
+#include "sirf/STIR/stir_x.h"
+#include "sirf/common/getenv.h"
+
+// Reg stuff
+#include "sirf/Reg/NiftiImageData3DDisplacement.h"
+
+using namespace sirf;
+using namespace stir;
+
+int test4()
+{
+	try {
+
+		std::string SIRF_path = sirf::getenv("SIRF_PATH");
+		if (SIRF_path.length() < 1) {
+			std::cout << "SIRF_PATH not defined, cannot find data" << std::endl;
+			return 1;
+		}
+
+        // Open acquisition data
+		std::string path = SIRF_path + "/data/examples/PET/";
+		std::string f_acq_data = path + "Utahscat600k_ca_seg4.hs";
+        auto acq_data_sptr = std::make_shared<PETAcquisitionDataInMemory>(f_acq_data.c_str());
+
+        // Recon image
+        auto im_sptr = std::make_shared<STIRImageData>(*acq_data_sptr);
+        im_sptr->fill(0.f);
+
+        // Create blank displacement field and convert to STIR format
+        NiftiImageData3DDisplacement<float> disp;
+        disp.create_from_3D_image(*im_sptr);
+        const std::string temp_fname = "temp_disp.nii";
+        disp.write(temp_fname);
+        auto disp_sptr = MAKE_SHARED<STIRDisp3DF>(temp_fname, 3 /*bspline_order*/);
+
+        // create additive term
+		auto a_sptr = acq_data_sptr->new_acquisition_data();
+		a_sptr->fill(0.05f);
+		// create background term
+		auto b_sptr = acq_data_sptr->new_acquisition_data();
+		b_sptr->fill(0.1f);
+
+        // Create acquisition model
+        auto matrix_sptr = MAKE_SHARED<RayTracingMatrix>();
+        matrix_sptr->set_num_tangential_LORs(2);
+        auto am_sptr = MAKE_SHARED<AcqModUsingMatrix3DF>();
+        am_sptr->set_matrix(matrix_sptr);
+        am_sptr->set_additive_term(a_sptr);
+        am_sptr->set_background_term(b_sptr);
+        am_sptr->set_up(acq_data_sptr, im_sptr);
+
+        // Create poisson loglikelihood
+        PoissonLogLhLinModMeanGatedProjDataWMotion3DF obj_fn;
+
+        // Add motion gates
+        const unsigned num_motion_gates = 4;
+        for (unsigned i=0; i<num_motion_gates; ++i)
+            obj_fn.add_gate(acq_data_sptr, am_sptr, disp_sptr);
+
+        obj_fn.set_up(im_sptr->data_sptr());
+
+        // Recon
+        OSMAPOSLReconstruction<Image3DF> recon;
+        recon.set_up(im_sptr->data_sptr());
+        recon.set_num_subiterations(2);
+        recon.reconstruct(im_sptr->data_sptr());
+
+        std::cout << "\nSuccess!\n";
+
+        return 0;
+	}
+	catch (...)
+	{
+		return 1;
+	}
+}
+
+//int test5();
+
+int main()
+{
+	return test4();
+	//return test5();
+}
diff --git a/src/xSTIR/cSTIR/cstir.cpp b/src/xSTIR/cSTIR/cstir.cpp
index 172035553..2c29f1ce2 100644
--- a/src/xSTIR/cSTIR/cstir.cpp
+++ b/src/xSTIR/cSTIR/cstir.cpp
@@ -89,6 +89,12 @@ void* cSTIR_newObject(const char* name)
 			"PoissonLogLikelihoodWithLinearModelForMeanAndProjData"))
 			return NEW_OBJECT_HANDLE
 			(xSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndProjData3DF);
+#ifdef STIR_GATED_MOTION
+        if (boost::iequals(name,
+			"PoissonLogLhLinModMeanGatedProjDataWMotion3DF"))
+			return NEW_OBJECT_HANDLE
+			(PoissonLogLhLinModMeanGatedProjDataWMotion3DF);
+#endif
 		if (boost::iequals(name, "AcqModUsingMatrix"))
 			return NEW_OBJECT_HANDLE(AcqModUsingMatrix3DF);
 #ifdef STIR_WITH_NiftyPET_PROJECTOR
@@ -884,6 +890,45 @@ cSTIR_objectiveFunctionGradientNotDivided(void* ptr_f, void* ptr_i, int subset)
 	CATCH;
 }
 
+extern "C"
+void*
+cSTIR_PoissonGatedWMotion_add_gate(
+        void* poisson_ptr, void* ad_ptr, void* am_ptr,
+        const char* disp_fname, const int b_spline_order)
+{
+	try {
+#ifdef STIR_GATED_MOTION
+        SPTR_FROM_HANDLE(PoissonLogLhLinModMeanGatedProjDataWMotion3DF,
+                         poisson_sptr, poisson_ptr);
+        SPTR_FROM_HANDLE(PETAcquisitionData, ad_sptr, ad_ptr);
+        SPTR_FROM_HANDLE(PETAcquisitionModel, am_sptr, am_ptr);
+        // set it
+        poisson_sptr->add_gate(ad_sptr, am_sptr, disp_fname, b_spline_order);
+#else
+		throw std::runtime_error("cSTIR_PoissonGatedWMotion_add_gate: shouldn't be here");
+#endif
+		return new DataHandle;
+	}
+	CATCH;
+}
+
+extern "C"
+void*
+cSTIR_PoissonGatedWMotion_clear_gates(void* poisson_ptr)
+{
+	try {
+#ifdef STIR_GATED_MOTION
+        SPTR_FROM_HANDLE(PoissonLogLhLinModMeanGatedProjDataWMotion3DF,
+                         poisson_sptr, poisson_ptr);
+        poisson_sptr->clear_gates();
+#else
+		throw std::runtime_error("cSTIR_PoissonGatedWMotion_clear_gates: shouldn't be here");
+#endif
+		return new DataHandle;
+	}
+	CATCH;
+}
+
 extern "C"
 void*
 cSTIR_setupPrior(void* ptr_p, void* ptr_i)
diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h
index 5969ca83d..069903d95 100644
--- a/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h
+++ b/src/xSTIR/cSTIR/include/sirf/STIR/cstir.h
@@ -119,6 +119,12 @@ extern "C" {
 	void* cSTIR_objectiveFunctionGradientNotDivided
 		(void* ptr_f, void* ptr_i, int subset);
 
+    // Poisson gated proj data w motion
+	void* cSTIR_PoissonGatedWMotion_add_gate(
+            void* poisson_ptr, void* ad_ptr, void* am_ptr,
+            const char* disp_fname, const int b_spline_order);
+    void* cSTIR_PoissonGatedWMotion_clear_gates(void* poisson_ptr);
+
 	// Prior methods
 	void* cSTIR_setupPrior(void* ptr_p, void* ptr_i);
 	void* cSTIR_priorGradient(void* ptr_p, void* ptr_i);
diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h
index 885606221..5036a4180 100644
--- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h
+++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_types.h
@@ -67,6 +67,11 @@ limitations under the License.
 #ifdef STIR_WITH_NiftyPET_PROJECTOR
 #include "stir/recon_buildblock/NiftyPET_projector/ProjectorByBinPairUsingNiftyPET.h"
 #endif
+#ifdef STIR_GATED_MOTION
+#include "stir_experimental/motion/PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotionUsingAdjoint.h"
+#include "stir_experimental/motion/NonRigidObjectTransformationUsingBSplines.h"
+#include "stir_experimental/motion/Transform3DObjectImageProcessor.h"
+#endif
 
 #include "stir/StirException.h"
 #include "stir/TextWriter.h"
@@ -98,6 +103,12 @@ namespace sirf {
 	typedef stir::QuadraticPrior<float> QuadPrior3DF;
 	typedef stir::DataProcessor<Image3DF> DataProcessor3DF;
 	typedef stir::TruncateToCylindricalFOVImageProcessor<float> CylindricFilter3DF;
+#ifdef STIR_GATED_MOTION
+    typedef stir::ObjectTransformation<3,float> STIRTrans3DF;
+    typedef stir::NonRigidObjectTransformationUsingBSplines<3,float> STIRDisp3DF;
+    typedef stir::PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotionUsingAdjoint < Image3DF >
+        stirPoissonLogLhLinModMeanGatedWMotion3DF;
+#endif
 
 }
 
diff --git a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h
index 51d4225c8..324bda7ab 100644
--- a/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h
+++ b/src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h
@@ -400,8 +400,8 @@ The actual algorithm is described in
 		}
 
 		virtual stir::Succeeded set_up(
-			stir::shared_ptr<PETAcquisitionData> sptr_acq,
-			stir::shared_ptr<STIRImageData> sptr_image);
+			const stir::shared_ptr<PETAcquisitionData> sptr_acq,
+			const stir::shared_ptr<STIRImageData> sptr_image);
 
 		// computes and returns a subset of forward-projected data 
 		stir::shared_ptr<PETAcquisitionData>
@@ -456,8 +456,8 @@ The actual algorithm is described in
 				get_proj_matrix_sptr();
 		}
 		virtual stir::Succeeded set_up(
-			stir::shared_ptr<PETAcquisitionData> sptr_acq,
-			stir::shared_ptr<STIRImageData> sptr_image)
+			const stir::shared_ptr<PETAcquisitionData> sptr_acq,
+			const stir::shared_ptr<STIRImageData> sptr_image)
 		{
 			if (!sptr_matrix_.get())
 				return stir::Succeeded::no;
@@ -591,6 +591,43 @@ The actual algorithm is described in
 	typedef xSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndProjData3DF
 		PoissonLogLhLinModMeanProjData3DF;
 
+#ifdef STIR_GATED_MOTION
+    class xSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion3DF :
+		public stirPoissonLogLhLinModMeanGatedWMotion3DF {
+    private:
+        typedef stirPoissonLogLhLinModMeanGatedWMotion3DF base_type;
+	public:
+
+        /// Add a motion gate. Required: acq data, acq model and disp field
+        void add_gate(
+                stir::shared_ptr<PETAcquisitionData> ad_sptr,
+                stir::shared_ptr<AcqMod3DF> am_sptr,
+                stir::shared_ptr<STIRTrans3DF> trans_sptr);
+
+        /// Add a motion gate. Required: acq data, acq model,
+        /// disp field (as string) and the b-spline order
+        void add_gate(
+                stir::shared_ptr<PETAcquisitionData> ad_sptr,
+                stir::shared_ptr<AcqMod3DF> am_sptr,
+                const std::string &disp_fname,
+                const int b_spline_order);
+
+        /// Set up
+        virtual stir::Succeeded set_up(stir::shared_ptr<Image3DF> const& target_sptr);
+
+        /// Clear all gates
+        void clear_gates();
+
+	private:
+		std::vector<stir::shared_ptr<PETAcquisitionData> > _ad_vec;
+		std::vector<stir::shared_ptr<AcqMod3DF> > _am_vec;
+        std::vector<stir::shared_ptr<STIRTrans3DF> > _trans_vec;
+	};
+
+	typedef xSTIR_PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion3DF
+		PoissonLogLhLinModMeanGatedProjDataWMotion3DF;
+#endif
+
 	class xSTIR_IterativeReconstruction3DF :
 		public stir::IterativeReconstruction < Image3DF > {
 	public:
@@ -661,7 +698,7 @@ The actual algorithm is described in
 				("wrong frequency cut-off", __FILE__, __LINE__);
 			fc_ramp = fc;
 		}
-		stir::Succeeded set_up(stir::shared_ptr<STIRImageData> sptr_id)
+		stir::Succeeded set_up(const stir::shared_ptr<STIRImageData> sptr_id)
 		{
 			_sptr_image_data.reset(new STIRImageData(*sptr_id));
 			_is_set_up = true;
diff --git a/src/xSTIR/cSTIR/stir_x.cpp b/src/xSTIR/cSTIR/stir_x.cpp
index 0d3dda717..f512c1194 100644
--- a/src/xSTIR/cSTIR/stir_x.cpp
+++ b/src/xSTIR/cSTIR/stir_x.cpp
@@ -523,8 +523,8 @@ PETAttenuationModel::normalise(PETAcquisitionData& ad) const
 
 Succeeded 
 PETAcquisitionModel::set_up(
-	shared_ptr<PETAcquisitionData> sptr_acq,
-	shared_ptr<STIRImageData> sptr_image)
+	const shared_ptr<PETAcquisitionData> sptr_acq,
+	const shared_ptr<STIRImageData> sptr_image)
 {
 	Succeeded s = Succeeded::no;
 	if (sptr_projectors_.get()) {
@@ -668,3 +668,113 @@ PETAcquisitionModel::backward(PETAcquisitionData& ad,
 
 	return sptr_id;
 }
+
+#ifdef STIR_GATED_MOTION
+void
+PoissonLogLhLinModMeanGatedProjDataWMotion3DF::
+add_gate(
+        stir::shared_ptr<PETAcquisitionData> ad_sptr,
+        stir::shared_ptr<AcqMod3DF> am_sptr,
+        stir::shared_ptr<STIRTrans3DF> trans_sptr)
+{
+    _ad_vec.push_back(ad_sptr);
+    _am_vec.push_back(am_sptr);
+    _trans_vec.push_back(trans_sptr);
+}
+
+void
+PoissonLogLhLinModMeanGatedProjDataWMotion3DF::
+add_gate(
+        stir::shared_ptr<PETAcquisitionData> ad_sptr,
+        stir::shared_ptr<AcqMod3DF> am_sptr,
+        const std::string &disp_fname,
+        const int b_spline_order)
+{
+    auto disp_sptr = MAKE_SHARED<STIRDisp3DF>(disp_fname, b_spline_order);
+    add_gate(ad_sptr, am_sptr, disp_sptr);
+}
+
+static stir::shared_ptr<GatedProjData>
+vec_AcqData_to_GatedProjData(const std::vector<stir::shared_ptr<PETAcquisitionData> > &vec_ad)
+{
+    if (vec_ad.empty() || vec_ad.at(0) == nullptr)
+        return nullptr;
+
+    auto sino_sptr = MAKE_SHARED<stir::GatedProjData>();
+    sino_sptr->set_exam_info(*vec_ad.at(0)->get_exam_info_sptr());
+    sino_sptr->resize(vec_ad.size());
+    for (unsigned i=0; i<vec_ad.size(); ++i)
+        sino_sptr->set_proj_data_sptr(vec_ad.at(i)->data(), i+1);
+    return sino_sptr;
+}
+
+template<class Type>
+void check_all_or_none(const std::vector<stir::shared_ptr<Type> > &vec, const std::string &explanation)
+{
+    const bool first_is_null_ptr = vec.at(0) == nullptr;
+    for (unsigned i=1; i<vec.size(); ++i)
+        if ((vec.at(i) == nullptr) != first_is_null_ptr)
+            throw std::runtime_error("PoissonLogLhLinModMeanGatedProjDataWMotion3DF::set_up "
+                                     + explanation + ": expected all or none");
+}
+
+Succeeded
+PoissonLogLhLinModMeanGatedProjDataWMotion3DF::
+set_up(const stir::shared_ptr<Image3DF> &target_sptr)
+{
+    const unsigned num_gates = _ad_vec.size();
+    if (num_gates==0)
+        throw std::runtime_error("PoissonLogLhLinModMeanGatedProjDataWMotion3DF::set_up: no gates added.");
+
+    // Projector pair
+    set_projector_pair_sptr(_am_vec.at(0)->projectors_sptr());
+
+    // Raw data: vector<sinograms> -> GatedProjData
+    auto sino_sptr = vec_AcqData_to_GatedProjData(_ad_vec);
+    set_proj_data_sptr(sino_sptr);
+
+    // Additives
+    std::vector<stir::shared_ptr<PETAcquisitionData> > add_vec;
+    for (unsigned i=0; i<num_gates; ++i)
+        add_vec.push_back(_am_vec.at(i)->additive_term_sptr());
+    check_all_or_none(add_vec, "additive");
+    auto add_sptr = vec_AcqData_to_GatedProjData(add_vec);
+    set_additive_proj_data_sptr(add_sptr);
+
+    // Norms
+    std::vector<stir::shared_ptr<stir::BinNormalisation> > norm_vec;
+    for (unsigned i=0; i<num_gates; ++i)
+        norm_vec.push_back(_am_vec.at(i)->normalisation_sptr());
+    check_all_or_none(norm_vec, "norms");
+    if (norm_vec.at(0)) {
+        stir::VectorWithOffset<stir::shared_ptr<stir::BinNormalisation> > norm_stir_vec(1, num_gates);
+        for (unsigned i=1; i<=num_gates; ++i)
+            norm_stir_vec.at(i) = norm_vec.at(i-1);
+        set_normalisations(norm_stir_vec);
+    }
+
+    // Transformations
+    stir::VectorWithOffset<stir::shared_ptr<DataProcessor3DF> > trans_stir_vec(1,num_gates);
+    for (unsigned i=0; i<num_gates; ++i)
+        trans_stir_vec.at(i+1) = MAKE_SHARED<stir::Transform3DObjectImageProcessor<float> >(_trans_vec.at(i));
+    set_forward_transformations(trans_stir_vec);
+
+    // Blank TimeFrameDefinitions
+    std::vector<std::pair<double, double> > frame_times(num_gates);
+    for (unsigned i=0; i<num_gates; ++i)
+        frame_times[i] = std::pair<double,double>(i,i+1);
+    set_frame_definitions(frame_times);
+
+    // Call the base type set up
+    return base_type::set_up(target_sptr);
+}
+
+void
+PoissonLogLhLinModMeanGatedProjDataWMotion3DF::
+clear_gates()
+{
+    _ad_vec.clear();
+    _am_vec.clear();
+    _trans_vec.clear();
+}
+#endif
diff --git a/src/xSTIR/pSTIR/STIR.py.in b/src/xSTIR/pSTIR/STIR.py.in
index 5e9d2412f..66887e03a 100644
--- a/src/xSTIR/pSTIR/STIR.py.in
+++ b/src/xSTIR/pSTIR/STIR.py.in
@@ -42,6 +42,10 @@ import sirf.pystir as pystir
 
 import sirf.STIR_params as parms
 
+if @STIR_GATED_MOTION@:
+    import sirf.Reg as reg
+    import uuid
+
 try:
     input = raw_input
 except NameError:
@@ -1789,6 +1793,40 @@ class PoissonLogLikelihoodWithLinearModelForMeanAndProjData\
         parms.set_parameter\
             (self.handle, self.name, 'acquisition_data', ad.handle)
 
+
+if @STIR_GATED_MOTION@:
+    class PoissonLogLhLinModMeanGatedProjDataWMotion(ObjectiveFunction):
+        """
+        Class for STIR PoissonLogLikelihoodWithLinearModelForMeanAndGatedProjDataWithMotion object.
+        There are two versions of this class in the STIR code base. The one used here lives in the
+        experimental part, and was chosen becase it uses the adjoint for the motion, as opposed to the
+        version in the main code base, which uses the inverse.
+        """
+        def __init__(self):
+            """Init"""
+            self.handle = None
+            self.name = 'PoissonLogLhLinModMeanGatedProjDataWMotion3DF'
+            self.handle = pystir.cSTIR_newObject(self.name)
+            check_status(self.handle)
+        def __del__(self):
+            if self.handle is not None:
+                pyiutil.deleteDataHandle(self.handle)
+        def add_gate(self, acq_data, acq_model, disp, b_spline_order=1):
+            """Add a gate by giving an acquisition data, acquisition model, 
+            displacement filename, and optionally b-spline order."""
+            assert_validity(acq_data, AcquisitionData)
+            assert_validity(acq_model, AcquisitionModel)
+            assert_validity(disp, reg.NiftiImageData3DDisplacement)
+            disp_fname = "temp_disp_" + uuid.uuid4().hex + ".nii"
+            disp.write(disp_fname)
+            try_calling(pystir.cSTIR_PoissonGatedWMotion_add_gate(
+                self.handle, acq_data.handle, acq_model.handle, disp_fname, b_spline_order))
+            os.remove(disp_fname)
+        def clear_gates(self):
+            """Clear all gates."""
+            try_calling(pystir.cSTIR_PoissonGatedWMotion_clear_gates(self.handle))
+    
+
 class Reconstructor(object):
     '''
     Class for a generic PET reconstructor.