.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "auto_examples/anatomical/example_anat_motion_EPI.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_auto_examples_anatomical_example_anat_motion_EPI.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_auto_examples_anatomical_example_anat_motion_EPI.py:


Creating motion artifacts on an anatomical EPI with SNAKE-fMRI
==============================================================

This examples walks through the elementary components of SNAKE and demonstrates how
to add motion artifacts to the simulation.

.. GENERATED FROM PYTHON SOURCE LINES 12-21

.. code-block:: Python


    # Imports
    import numpy as np
    import matplotlib.pyplot as plt

    from snake.core.phantom import Phantom
    from snake.core.sampling import EPI3dAcquisitionSampler
    from snake.core.simulation import GreConfig, SimConfig, default_hardware








.. GENERATED FROM PYTHON SOURCE LINES 22-24

Setting up the base simulation Config.  This configuration holds all key
parameters for the simulation, describing the scanner parameters.

.. GENERATED FROM PYTHON SOURCE LINES 24-34

.. code-block:: Python


    sim_conf = SimConfig(
        max_sim_time=6,
        seq=GreConfig(TR=50, TE=30, FA=3),
        hardware=default_hardware,
    )
    sim_conf.hardware.n_coils = 8
    sim_conf.fov.res_mm = (3, 3, 3)
    sim_conf






.. raw:: html

    <div class="output_subarea output_html rendered_html output_result">
    <table style="border:1px solid lightgray;"><caption style="border:1px solid lightgray;"><strong>SimConfig</strong></caption>
    <tr><td>max_sim_time(<i>float</i>)</td><td>6</td></tr>
    <tr><td>seq(<i>GreConfig</i>)</td><td><table style="border:1px solid lightgray;"><caption style="border:1px solid lightgray;"><strong>GreConfig</strong></caption>
    <tr><td>TR (<i>float</i>)</td><td>TE (<i>float</i>)</td><td>FA (<i>float</i>)</td></tr>
    <tr><td>50</td><td>30</td><td>3</td></tr>
    </table></td></tr>
    <tr><td>hardware(<i>HardwareConfig</i>)</td><td><table style="border:1px solid lightgray;"><caption style="border:1px solid lightgray;"><strong>HardwareConfig</strong></caption>
    <tr><td>gmax (<i>float</i>)</td><td>smax (<i>float</i>)</td><td>n_coils (<i>int</i>)</td><td>dwell_time_ms (<i>float</i>)</td><td>raster_time_ms (<i>float</i>)</td><td>field (<i>float</i>)</td></tr>
    <tr><td>40</td><td>200</td><td>8</td><td>0.001</td><td>0.005</td><td>3.0</td></tr>
    </table></td></tr>
    <tr><td>fov(<i>FOVConfig</i>)</td><td><table style="border:1px solid lightgray;"><caption style="border:1px solid lightgray;"><strong>FOVConfig</strong></caption>
    <tr><td>size (<i>ThreeFloats</i>)</td><td>offset (<i>ThreeFloats</i>)</td><td>angles (<i>ThreeFloats</i>)</td><td>res_mm (<i>ThreeFloats</i>)</td></tr>
    <tr><td>(181, 217, 181)</td><td>(-90.25, -126.25, -72.25)</td><td>(0, 0, 0)</td><td>(3, 3, 3)</td></tr>
    </table></td></tr>
    <tr><td>rng_seed(<i>int</i>)</td><td>19290506</td></tr>
    </table>
    </div>
    <br />
    <br />

.. GENERATED FROM PYTHON SOURCE LINES 35-44

Creating the base Phantom
-------------------------

The simulation acquires
the data describe in a phantom. A phantom consists of fuzzy segmentation of
head tissue, and their MR intrinsic parameters
(density, T1, T2, T2*,  magnetic susceptibilities)

Here we use Brainweb reference mask and values for convenience.

.. GENERATED FROM PYTHON SOURCE LINES 44-53

.. code-block:: Python


    phantom = Phantom.from_brainweb(
        sub_id=4, sim_conf=sim_conf, output_res=1, tissue_file="tissue_7T"
    )

    # Here are the tissue availables and their parameters
    phantom.affine






.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    No matrix size found in the header. The header is probably missing.

    array([[1., 0., 0., 0.],
           [0., 1., 0., 0.],
           [0., 0., 1., 0.],
           [0., 0., 0., 1.]], dtype=float32)



.. GENERATED FROM PYTHON SOURCE LINES 54-58

Adding motion to the Phantom
----------------------------
The motion is added to the phantom by applying a transformation to the
simulation's FOV configuration (as if the phantom was moving in the scanner).

.. GENERATED FROM PYTHON SOURCE LINES 58-71

.. code-block:: Python


    from snake.core.handlers import RandomMotionImageHandler

    motion = RandomMotionImageHandler(
        ts_std_mms=[1, 1, 1],
        rs_std_degs=[0.1, 0.1, 0.1],
    )

    motion2 = RandomMotionImageHandler(
        ts_std_mms=[3, 3, 3],
        rs_std_degs=[0.5, 0.5, 0.5],
    )








.. GENERATED FROM PYTHON SOURCE LINES 72-74

Setting up Acquisition Pattern and Initializing Result file.
------------------------------------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 74-82

.. code-block:: Python


    # The next piece of simulation is the acquisition trajectory.
    # Here nothing fancy, we are using a EPI (fully sampled), that samples a 3D
    # k-space (this akin to the 3D EPI sequence of XXXX)

    sampler = EPI3dAcquisitionSampler(accelz=1, acsz=0.1, orderz="top-down")









.. GENERATED FROM PYTHON SOURCE LINES 83-97

Acquisition with Cartesian Engine
---------------------------------

The generated file ``example_EPI.mrd`` does not contains any k-space data for
now, only the sampling trajectory. let's put some in. In order to do so, we
need to setup the **acquisition engine** that models the MR physics, and get
sampled at the specified k-space trajectory.

SNAKE comes with two models for the MR Physics:

- model="simple" :: Each k-space shot acquires a constant signal, which is the
  image contrast at TE.
- model="T2s" :: Each k-space shot is degraded by the T2* decay induced by
  each tissue.

.. GENERATED FROM PYTHON SOURCE LINES 97-143

.. code-block:: Python


    # Here we will use the "simple" model, which is faster.
    #
    # SNAKE's Engine are capable of simulating the data in parallel, by distributing
    # the shots to be acquired to a set of processes. To do so , we need to specify
    # the number of jobs that will run in parallel, as well as the size of a job.
    # Setting the job size and the number of jobs can have a great impact on total
    # runtime and memory consumption.
    #
    # Here, we have a single frame to acquire with 60 frames (one EPI per slice), so
    # a single worker will do.

    from snake.core.engine import EPIAcquisitionEngine

    engine = EPIAcquisitionEngine(model="simple")

    engine(
        "example_nomotion.mrd",
        sampler=sampler,
        phantom=phantom,
        handlers=[],
        sim_conf=sim_conf,
        worker_chunk_size=16,
        n_workers=4,
    )

    engine(
        "example_motion.mrd",
        sampler=sampler,
        phantom=phantom,
        handlers=[motion],
        sim_conf=sim_conf,
        worker_chunk_size=16,
        n_workers=4,
    )

    engine(
        "example_motion2.mrd",
        sampler=sampler,
        phantom=phantom,
        handlers=[motion2],
        sim_conf=sim_conf,
        worker_chunk_size=16,
        n_workers=4,
    )



.. rst-class:: sphx-glr-script-out

.. code-block:: pytb

    Traceback (most recent call last):
      File "/volatile/github-ci-mind-inria/gpu_mind_runner/_work/snake-fmri/snake-fmri/examples/anatomical/example_anat_motion_EPI.py", line 113, in <module>
        engine(
      File "/volatile/github-ci-mind-inria/gpu_mind_runner/_work/snake-fmri/snake-fmri/src/snake/core/engine/base.py", line 299, in __call__
        raise exc
      File "/volatile/github-ci-mind-inria/gpu_mind_runner/_work/snake-fmri/snake-fmri/src/snake/core/engine/base.py", line 296, in __call__
        f_chunk = str(future.result())
      File "/volatile/github-ci-mind-inria/gpu_mind_runner/_work/_tool/Python/3.10.16/x64/lib/python3.10/concurrent/futures/_base.py", line 451, in result
        return self.__get_result()
      File "/volatile/github-ci-mind-inria/gpu_mind_runner/_work/_tool/Python/3.10.16/x64/lib/python3.10/concurrent/futures/_base.py", line 403, in __get_result
        raise self._exception
    IndexError: index 1 is out of bounds for axis 0 with size 1




.. GENERATED FROM PYTHON SOURCE LINES 144-146

Simple reconstruction
---------------------

.. GENERATED FROM PYTHON SOURCE LINES 146-166

.. code-block:: Python


    from snake.mrd_utils import CartesianFrameDataLoader
    from snake.toolkit.reconstructors import ZeroFilledReconstructor


    with CartesianFrameDataLoader("example_nomotion.mrd") as data_loader:
        rec = ZeroFilledReconstructor(n_jobs=1)
        rec_nomotion = rec.reconstruct(data_loader).squeeze()

    with CartesianFrameDataLoader("example_motion.mrd") as data_loader:
        rec = ZeroFilledReconstructor(n_jobs=1)
        rec_motion = rec.reconstruct(data_loader).squeeze()
        motion = data_loader.get_dynamic(0)

    with CartesianFrameDataLoader("example_motion2.mrd") as data_loader:
        rec = ZeroFilledReconstructor(n_jobs=1)
        rec_motion2 = rec.reconstruct(data_loader).squeeze()
        motion2 = data_loader.get_dynamic(0)



.. GENERATED FROM PYTHON SOURCE LINES 167-170

.. code-block:: Python

    with CartesianFrameDataLoader("example_motion2.mrd") as data_loader:
        sim_conf.max_sim_time


.. GENERATED FROM PYTHON SOURCE LINES 171-173

.. code-block:: Python

    rec_motion.shape


.. GENERATED FROM PYTHON SOURCE LINES 174-176

Show the motion
---------------

.. GENERATED FROM PYTHON SOURCE LINES 176-194

.. code-block:: Python

    fig, axs = plt.subplots(2, 1, figsize=(10, 5), sharex=True)
    time = np.arange(len(motion.data[0])) * sim_conf.seq.TR / 1000
    TR_vol = 3  # We generated 2 frames in 6 seconds
    axs[0].axvline(TR_vol, c="gray")
    axs[1].axvline(TR_vol, c="gray")


    axs[0].plot(time, motion.data[:3, :].T)
    axs[1].plot(time, motion.data[3:, :].T)
    axs[0].set_prop_cycle(None)  # reset color cycling
    axs[1].set_prop_cycle(None)
    axs[0].plot(time, motion2.data[:3, :].T, linestyle="dashed")
    axs[1].plot(time, motion2.data[3:, :].T, linestyle="dashed")

    axs[1].set_xlabel("time (s)")
    axs[0].set_ylabel("Translation (mm)")
    axs[1].set_ylabel("Rotation (deg)")


.. GENERATED FROM PYTHON SOURCE LINES 195-197

Visualizing the reconstructed data
----------------------------------

.. GENERATED FROM PYTHON SOURCE LINES 197-227

.. code-block:: Python


    import matplotlib.pyplot as plt

    from snake.toolkit.plotting import axis3dcut

    fig, axs = plt.subplots(1, 3, figsize=(26, 5))

    axis3dcut(
        rec_motion[0].T,
        None,
        None,
        cbar=False,
        cuts=(0.5, 0.5, 0.5),
        ax=axs[0],
    )
    axis3dcut(
        rec_motion2[0].T,
        None,
        None,
        cbar=False,
        cuts=(0.5, 0.5, 0.5),
        ax=axs[1],
    )
    axis3dcut(rec_nomotion[0].T, None, None, cbar=False, cuts=(0.5, 0.5, 0.5), ax=axs[2])
    axs[0].set_title("Small motion")
    axs[1].set_title("Big motion")
    axs[2].set_title("No motion")
    plt.show()



.. GENERATED FROM PYTHON SOURCE LINES 228-231

.. code-block:: Python

    sim_conf.seq.TR_eff




.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 19.578 seconds)


.. _sphx_glr_download_auto_examples_anatomical_example_anat_motion_EPI.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: example_anat_motion_EPI.ipynb <example_anat_motion_EPI.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: example_anat_motion_EPI.py <example_anat_motion_EPI.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: example_anat_motion_EPI.zip <example_anat_motion_EPI.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_