diff --git a/.github/dependabot.yml b/.github/dependabot.yml index 11d0cb6..142397f 100644 --- a/.github/dependabot.yml +++ b/.github/dependabot.yml @@ -3,7 +3,7 @@ updates: - package-ecosystem: "github-actions" directory: "/" schedule: - interval: "daily" + interval: "monthly" groups: all-actions: patterns: ["*"] diff --git a/.github/workflows/examples.yml b/.github/workflows/examples.yml index 0832306..5b066a9 100644 --- a/.github/workflows/examples.yml +++ b/.github/workflows/examples.yml @@ -25,7 +25,7 @@ jobs: run: | pip install --index-url=https://download.pytorch.org/whl/cpu \ --extra-index-url https://pypi.org/simple/ \ - --upgrade --upgrade-strategy eager -e .[examples] + --upgrade --upgrade-strategy eager -e .[dev] # Build the book - name: Build the book diff --git a/examples/grpe_flash_dixon.ipynb b/examples/grpe_flash_dixon.ipynb index 6a49333..41abbb7 100644 --- a/examples/grpe_flash_dixon.ipynb +++ b/examples/grpe_flash_dixon.ipynb @@ -45,7 +45,7 @@ "from mrpro.operators.models import SpoiledGRE\n", "from scipy.interpolate import interp1d\n", "\n", - "from mrseq.scripts.grpe_flash_dixon import main as create_seq\n", + "from mrseq.sequences.grpe_flash_dixon import main as create_seq\n", "from mrseq.utils import combine_ismrmrd_files\n", "from mrseq.utils import sys_defaults" ] diff --git a/examples/radial_flash.ipynb b/examples/radial_flash.ipynb index e94c724..2912767 100644 --- a/examples/radial_flash.ipynb +++ b/examples/radial_flash.ipynb @@ -4,6 +4,14 @@ "cell_type": "markdown", "id": "0", "metadata": {}, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PTB-MR/mrseq/blob/main/examples/radial_flash.ipynb)" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, "source": [ "# Radial FLASH\n", "\n", @@ -12,7 +20,7 @@ }, { "cell_type": "markdown", - "id": "1", + "id": "2", "metadata": {}, "source": [ "### Imports" @@ -21,7 +29,20 @@ { "cell_type": "code", "execution_count": null, - "id": "2", + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('mrseq'):\n", + " %pip install mrseq[examples]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -37,14 +58,14 @@ "from mrpro.data.traj_calculators import KTrajectoryIsmrmrd\n", "from mrpro.operators.models import SpoiledGRE\n", "\n", - "from mrseq.scripts.radial_flash import main as create_seq\n", + "from mrseq.sequences.radial_flash import main as create_seq\n", "from mrseq.utils import combine_ismrmrd_files\n", "from mrseq.utils import sys_defaults" ] }, { "cell_type": "markdown", - "id": "3", + "id": "5", "metadata": {}, "source": [ "### Settings\n", @@ -55,7 +76,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -69,7 +90,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "### Create the digital phantom\n", @@ -80,7 +101,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -90,7 +111,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ "### Create the radial FLASH sequence\n", @@ -101,7 +122,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -119,7 +140,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "### Simulate the sequence\n", @@ -130,7 +151,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -142,7 +163,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "### Reconstruct the image\n", @@ -152,7 +173,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -163,7 +184,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "### Compare to theoretical model\n", @@ -174,7 +195,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -184,7 +205,7 @@ { "cell_type": "code", "execution_count": null, - "id": "15", + "id": "17", "metadata": {}, "outputs": [], "source": [ diff --git a/examples/spiral_flash.ipynb b/examples/spiral_flash.ipynb index 77dbad9..b7edade 100644 --- a/examples/spiral_flash.ipynb +++ b/examples/spiral_flash.ipynb @@ -4,6 +4,14 @@ "cell_type": "markdown", "id": "0", "metadata": {}, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PTB-MR/mrseq/blob/main/examples/spiral_flash.ipynb)" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, "source": [ "# Spiral FLASH\n", "\n", @@ -12,7 +20,7 @@ }, { "cell_type": "markdown", - "id": "1", + "id": "2", "metadata": {}, "source": [ "### Imports" @@ -21,7 +29,20 @@ { "cell_type": "code", "execution_count": null, - "id": "2", + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('mrseq'):\n", + " %pip install mrseq[examples]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -37,14 +58,14 @@ "from mrpro.data.traj_calculators import KTrajectoryIsmrmrd\n", "from mrpro.operators.models import SpoiledGRE\n", "\n", - "from mrseq.scripts.spiral_flash import main as create_seq\n", + "from mrseq.sequences.spiral_flash import main as create_seq\n", "from mrseq.utils import combine_ismrmrd_files\n", "from mrseq.utils import sys_defaults" ] }, { "cell_type": "markdown", - "id": "3", + "id": "5", "metadata": {}, "source": [ "### Settings\n", @@ -55,7 +76,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -69,7 +90,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "### Create the digital phantom\n", @@ -80,7 +101,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -90,7 +111,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ "### Create the spiral FLASH sequence\n", @@ -101,7 +122,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -119,7 +140,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "### Simulate the sequence\n", @@ -130,7 +151,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -142,7 +163,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "### Reconstruct the image\n", @@ -152,7 +173,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -163,7 +184,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "### Compare to theoretical model\n", @@ -174,7 +195,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -184,7 +205,7 @@ { "cell_type": "code", "execution_count": null, - "id": "15", + "id": "17", "metadata": {}, "outputs": [], "source": [ diff --git a/examples/t1_inv_rec_gre_single_line.ipynb b/examples/t1_inv_rec_gre_single_line.ipynb index 4ce6734..f1a4d67 100644 --- a/examples/t1_inv_rec_gre_single_line.ipynb +++ b/examples/t1_inv_rec_gre_single_line.ipynb @@ -4,6 +4,14 @@ "cell_type": "markdown", "id": "0", "metadata": {}, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PTB-MR/mrseq/blob/main/examples/t1_inv_rec_gre_single_line.ipynb)" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, "source": [ "# T1 Mapping - IR GRE\n", "\n", @@ -13,7 +21,7 @@ }, { "cell_type": "markdown", - "id": "1", + "id": "2", "metadata": {}, "source": [ "### Imports" @@ -22,7 +30,20 @@ { "cell_type": "code", "execution_count": null, - "id": "2", + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('mrseq'):\n", + " %pip install mrseq[examples]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -42,13 +63,13 @@ "from mrpro.operators import DictionaryMatchOp\n", "from mrpro.operators.models import InversionRecovery\n", "\n", - "from mrseq.scripts.t1_inv_rec_gre_single_line import main as create_seq\n", + "from mrseq.sequences.t1_inv_rec_gre_single_line import main as create_seq\n", "from mrseq.utils import sys_defaults" ] }, { "cell_type": "markdown", - "id": "3", + "id": "5", "metadata": {}, "source": [ "### Settings\n", @@ -58,7 +79,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -71,7 +92,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "### Create the digital phantom\n", @@ -82,7 +103,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -92,7 +113,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ "### Create the IR-GRE sequence\n", @@ -103,7 +124,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -120,7 +141,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "### Simulate the sequence\n", @@ -131,7 +152,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -142,7 +163,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "### Reconstruct the images at different inversion times\n", @@ -153,7 +174,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -166,7 +187,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "We can now plot the images at different inversion times." @@ -175,7 +196,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -190,7 +211,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "17", "metadata": {}, "source": [ "### Estimate the T1 maps\n", @@ -200,7 +221,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "18", "metadata": {}, "outputs": [], "source": [ diff --git a/examples/t1_inv_rec_se_single_line.ipynb b/examples/t1_inv_rec_se_single_line.ipynb index 3775785..fffac79 100644 --- a/examples/t1_inv_rec_se_single_line.ipynb +++ b/examples/t1_inv_rec_se_single_line.ipynb @@ -4,6 +4,14 @@ "cell_type": "markdown", "id": "0", "metadata": {}, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PTB-MR/mrseq/blob/main/examples/t1_inv_rec_se_single_line.ipynb)" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, "source": [ "# T1 Mapping - IR SE\n", "\n", @@ -13,7 +21,7 @@ }, { "cell_type": "markdown", - "id": "1", + "id": "2", "metadata": {}, "source": [ "### Imports" @@ -22,7 +30,20 @@ { "cell_type": "code", "execution_count": null, - "id": "2", + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('mrseq'):\n", + " %pip install mrseq[examples]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -42,13 +63,13 @@ "from mrpro.operators import DictionaryMatchOp\n", "from mrpro.operators.models import InversionRecovery\n", "\n", - "from mrseq.scripts.t1_inv_rec_se_single_line import main as create_seq\n", + "from mrseq.sequences.t1_inv_rec_se_single_line import main as create_seq\n", "from mrseq.utils import sys_defaults" ] }, { "cell_type": "markdown", - "id": "3", + "id": "5", "metadata": {}, "source": [ "### Settings\n", @@ -58,7 +79,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -71,7 +92,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "### Create the digital phantom\n", @@ -82,7 +103,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -92,7 +113,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ "### Create the IR-SE sequence\n", @@ -103,7 +124,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -120,7 +141,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "### Simulate the sequence\n", @@ -130,7 +151,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -141,7 +162,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "### Reconstruct the images at different inversion times\n", @@ -152,7 +173,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -165,7 +186,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "We can now plot the images at different inversion times." @@ -174,7 +195,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -189,7 +210,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "17", "metadata": {}, "source": [ "### Estimate the T1 maps\n", @@ -199,7 +220,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "18", "metadata": {}, "outputs": [], "source": [ diff --git a/examples/t1_molli_bssfp.ipynb b/examples/t1_molli_bssfp.ipynb index 199d4c3..219f064 100644 --- a/examples/t1_molli_bssfp.ipynb +++ b/examples/t1_molli_bssfp.ipynb @@ -4,6 +4,14 @@ "cell_type": "markdown", "id": "0", "metadata": {}, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PTB-MR/mrseq/blob/main/examples/t1_molli_bssfp.ipynb)" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, "source": [ "# T1 Mapping - MOLLI\n", "\n", @@ -13,7 +21,7 @@ }, { "cell_type": "markdown", - "id": "1", + "id": "2", "metadata": {}, "source": [ "### Imports" @@ -22,7 +30,20 @@ { "cell_type": "code", "execution_count": null, - "id": "2", + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('mrseq'):\n", + " %pip install mrseq[examples]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -45,13 +66,13 @@ "from mrpro.operators.models import MOLLI\n", "from mrpro.phantoms.coils import birdcage_2d\n", "\n", - "from mrseq.scripts.t1_molli_bssfp import main as create_seq\n", + "from mrseq.sequences.t1_molli_bssfp import main as create_seq\n", "from mrseq.utils import sys_defaults" ] }, { "cell_type": "markdown", - "id": "3", + "id": "5", "metadata": {}, "source": [ "### Settings\n", @@ -61,7 +82,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -73,7 +94,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "### Create the digital phantom\n", @@ -84,7 +105,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -98,7 +119,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ "### Create the MOLLI bSSFP sequence\n", @@ -109,7 +130,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -125,7 +146,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "### Simulate the sequence\n", @@ -135,7 +156,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -146,7 +167,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "### Reconstruct the images at different inversion times\n", @@ -157,7 +178,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -178,7 +199,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "We visualize the coil sensitivity maps which we used for the iterative SENSE reconstruction." @@ -187,7 +208,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -198,7 +219,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "17", "metadata": {}, "source": [ "We can now plot the images at different inversion times." @@ -207,7 +228,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -229,7 +250,7 @@ }, { "cell_type": "markdown", - "id": "17", + "id": "19", "metadata": {}, "source": [ "### Estimate the T1 maps\n", @@ -239,7 +260,7 @@ { "cell_type": "code", "execution_count": null, - "id": "18", + "id": "20", "metadata": {}, "outputs": [], "source": [ diff --git a/examples/t1_t2_spiral_cmrf.ipynb b/examples/t1_t2_spiral_cmrf.ipynb index b3dcaf6..4be7412 100644 --- a/examples/t1_t2_spiral_cmrf.ipynb +++ b/examples/t1_t2_spiral_cmrf.ipynb @@ -4,6 +4,14 @@ "cell_type": "markdown", "id": "0", "metadata": {}, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PTB-MR/mrseq/blob/main/examples/t1_t2_spiral_cmrf.ipynb)" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, "source": [ "# T1 and T2 Mapping - Spiral cMRF\n", "\n", @@ -12,7 +20,7 @@ }, { "cell_type": "markdown", - "id": "1", + "id": "2", "metadata": {}, "source": [ "### Imports" @@ -21,7 +29,20 @@ { "cell_type": "code", "execution_count": null, - "id": "2", + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('mrseq'):\n", + " %pip install mrseq[examples]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -42,14 +63,14 @@ "from mrpro.operators import DictionaryMatchOp\n", "from mrpro.operators.models import CardiacFingerprinting\n", "\n", - "from mrseq.scripts.t1_t2_spiral_cmrf import main as create_seq\n", + "from mrseq.sequences.t1_t2_spiral_cmrf import main as create_seq\n", "from mrseq.utils import combine_ismrmrd_files\n", "from mrseq.utils import sys_defaults" ] }, { "cell_type": "markdown", - "id": "3", + "id": "5", "metadata": {}, "source": [ "### Settings\n", @@ -61,7 +82,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -74,7 +95,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "### Create the digital phantom\n", @@ -87,7 +108,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -99,7 +120,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ "### Create the spiral cardiac MRF sequence\n", @@ -110,7 +131,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -127,7 +148,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "### Simulate the sequence\n", @@ -137,7 +158,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -149,7 +170,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "### Reconstruct series of images showing the dynamic contrast change\n", @@ -160,7 +181,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -182,7 +203,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "We can now plot the images at different time points along the acquisition." @@ -191,7 +212,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -206,7 +227,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "17", "metadata": {}, "source": [ "### Estimate the T1 and T2 maps\n", @@ -216,7 +237,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "18", "metadata": {}, "outputs": [], "source": [ diff --git a/examples/t1rho_se_single_line.ipynb b/examples/t1rho_se_single_line.ipynb index 72e4644..77faaa0 100644 --- a/examples/t1rho_se_single_line.ipynb +++ b/examples/t1rho_se_single_line.ipynb @@ -4,6 +4,14 @@ "cell_type": "markdown", "id": "0", "metadata": {}, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PTB-MR/mrseq/blob/main/examples/t1rho_se_single_line.ipynb)" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, "source": [ "# T1rho Mapping\n", "\n", @@ -13,7 +21,7 @@ }, { "cell_type": "markdown", - "id": "1", + "id": "2", "metadata": {}, "source": [ "### Imports" @@ -22,7 +30,20 @@ { "cell_type": "code", "execution_count": null, - "id": "2", + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('mrseq'):\n", + " %pip install mrseq[examples]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -42,13 +63,13 @@ "from mrpro.operators import DictionaryMatchOp\n", "from mrpro.operators.models import MonoExponentialDecay\n", "\n", - "from mrseq.scripts.t1rho_se_single_line import main as create_seq\n", + "from mrseq.sequences.t1rho_se_single_line import main as create_seq\n", "from mrseq.utils import sys_defaults" ] }, { "cell_type": "markdown", - "id": "3", + "id": "5", "metadata": {}, "source": [ "### Settings\n", @@ -58,7 +79,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -70,7 +91,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "### Create the digital phantom\n", @@ -81,7 +102,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -92,7 +113,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ "### Create the Spin-lock SE sequence\n", @@ -103,7 +124,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -119,7 +140,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "### Simulate the sequence\n", @@ -129,7 +150,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -140,7 +161,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "### Reconstruct the images at different spin-lock times\n", @@ -150,7 +171,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -163,7 +184,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "We can now plot the images at different inversion times." @@ -172,7 +193,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -188,7 +209,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "17", "metadata": {}, "source": [ "### Estimate the T1rho maps\n", @@ -199,7 +220,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "18", "metadata": {}, "outputs": [], "source": [ diff --git a/examples/t2_multi_echo_se_single_line.ipynb b/examples/t2_multi_echo_se_single_line.ipynb index 3d5b44a..98a1e37 100644 --- a/examples/t2_multi_echo_se_single_line.ipynb +++ b/examples/t2_multi_echo_se_single_line.ipynb @@ -4,6 +4,14 @@ "cell_type": "markdown", "id": "0", "metadata": {}, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PTB-MR/mrseq/blob/main/examples/t2_multi_echo_se_single_line.ipynb)" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, "source": [ "# T2 Mapping - Multi-echo SE\n", "\n", @@ -13,7 +21,7 @@ }, { "cell_type": "markdown", - "id": "1", + "id": "2", "metadata": {}, "source": [ "### Imports" @@ -22,7 +30,20 @@ { "cell_type": "code", "execution_count": null, - "id": "2", + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('mrseq'):\n", + " %pip install mrseq[examples]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -42,13 +63,13 @@ "from mrpro.operators import DictionaryMatchOp\n", "from mrpro.operators.models import MonoExponentialDecay\n", "\n", - "from mrseq.scripts.t2_multi_echo_se_single_line import main as create_seq\n", + "from mrseq.sequences.t2_multi_echo_se_single_line import main as create_seq\n", "from mrseq.utils import sys_defaults" ] }, { "cell_type": "markdown", - "id": "3", + "id": "5", "metadata": {}, "source": [ "### Settings\n", @@ -58,7 +79,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -71,7 +92,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "### Create the digital phantom\n", @@ -82,7 +103,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -92,7 +113,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ "### Create the multi-echo SE sequence\n", @@ -103,7 +124,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -120,7 +141,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "### Simulate the sequence\n", @@ -130,7 +151,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -141,7 +162,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "### Reconstruct the images at different echo times\n", @@ -152,7 +173,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -165,7 +186,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "We can now plot the images at different inversion times." @@ -174,7 +195,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -189,7 +210,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "17", "metadata": {}, "source": [ "### Estimate the T2 maps\n", @@ -199,7 +220,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "18", "metadata": {}, "outputs": [], "source": [ diff --git a/examples/t2_t2prep_flash.ipynb b/examples/t2_t2prep_flash.ipynb index e27cc64..9d0dcc0 100644 --- a/examples/t2_t2prep_flash.ipynb +++ b/examples/t2_t2prep_flash.ipynb @@ -4,6 +4,14 @@ "cell_type": "markdown", "id": "0", "metadata": {}, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PTB-MR/mrseq/blob/main/examples/t2_t2prep_flash.ipynb)" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, "source": [ "# T2 Mapping - T2prep FLASH\n", "\n", @@ -12,7 +20,7 @@ }, { "cell_type": "markdown", - "id": "1", + "id": "2", "metadata": {}, "source": [ "### Imports" @@ -21,7 +29,20 @@ { "cell_type": "code", "execution_count": null, - "id": "2", + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('mrseq'):\n", + " %pip install mrseq[examples]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -41,13 +62,13 @@ "from mrpro.operators import DictionaryMatchOp\n", "from mrpro.operators.models import MonoExponentialDecay\n", "\n", - "from mrseq.scripts.t2_t2prep_flash import main as create_seq\n", + "from mrseq.sequences.t2_t2prep_flash import main as create_seq\n", "from mrseq.utils import sys_defaults" ] }, { "cell_type": "markdown", - "id": "3", + "id": "5", "metadata": {}, "source": [ "### Settings\n", @@ -57,7 +78,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -70,7 +91,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "### Create the digital phantom\n", @@ -82,7 +103,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -95,7 +116,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ "### Create the T2prep FLASH sequence\n", @@ -105,7 +126,7 @@ }, { "cell_type": "markdown", - "id": "8", + "id": "10", "metadata": {}, "source": [ "For in-vivo applications we have to make sure the sequence can be run within a breathhold. This would require undersampling (acceleration > 1) and obtaining a high number of phase encoding points in each cardiac cycle. This will impair the accuracy of the obtained T2 maps. For the evaluation here we can make the sequence longer to increase accuracy. " @@ -114,7 +135,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9", + "id": "11", "metadata": {}, "outputs": [], "source": [ @@ -133,7 +154,7 @@ }, { "cell_type": "markdown", - "id": "10", + "id": "12", "metadata": {}, "source": [ "### Simulate the sequence\n", @@ -143,7 +164,7 @@ { "cell_type": "code", "execution_count": null, - "id": "11", + "id": "13", "metadata": {}, "outputs": [], "source": [ @@ -154,7 +175,7 @@ }, { "cell_type": "markdown", - "id": "12", + "id": "14", "metadata": {}, "source": [ "### Reconstruct the images with different T2-preparation pulses.\n", @@ -165,7 +186,7 @@ { "cell_type": "code", "execution_count": null, - "id": "13", + "id": "15", "metadata": {}, "outputs": [], "source": [ @@ -178,7 +199,7 @@ }, { "cell_type": "markdown", - "id": "14", + "id": "16", "metadata": {}, "source": [ "We can now plot the images with different T2-preparation times." @@ -187,7 +208,7 @@ { "cell_type": "code", "execution_count": null, - "id": "15", + "id": "17", "metadata": {}, "outputs": [], "source": [ @@ -202,7 +223,7 @@ }, { "cell_type": "markdown", - "id": "16", + "id": "18", "metadata": {}, "source": [ "### Estimate the T2 maps\n", @@ -212,7 +233,7 @@ { "cell_type": "code", "execution_count": null, - "id": "17", + "id": "19", "metadata": {}, "outputs": [], "source": [ diff --git a/examples/t2star_multi_echo_flash.ipynb b/examples/t2star_multi_echo_flash.ipynb index f41e875..09d33b7 100644 --- a/examples/t2star_multi_echo_flash.ipynb +++ b/examples/t2star_multi_echo_flash.ipynb @@ -4,6 +4,14 @@ "cell_type": "markdown", "id": "0", "metadata": {}, + "source": [ + "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/PTB-MR/mrseq/blob/main/examples/t2star_multi_echo_flash.ipynb)" + ] + }, + { + "cell_type": "markdown", + "id": "1", + "metadata": {}, "source": [ "# T2* Mapping - Multi-echo FLASH\n", "\n", @@ -12,7 +20,7 @@ }, { "cell_type": "markdown", - "id": "1", + "id": "2", "metadata": {}, "source": [ "### Imports" @@ -21,7 +29,20 @@ { "cell_type": "code", "execution_count": null, - "id": "2", + "id": "3", + "metadata": {}, + "outputs": [], + "source": [ + "import importlib\n", + "\n", + "if not importlib.util.find_spec('mrseq'):\n", + " %pip install mrseq[examples]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4", "metadata": {}, "outputs": [], "source": [ @@ -44,14 +65,14 @@ "from mrpro.operators.models import MonoExponentialDecay\n", "from mrpro.phantoms.coils import birdcage_2d\n", "\n", - "from mrseq.scripts.t2star_multi_echo_flash import main as create_seq\n", + "from mrseq.sequences.t2star_multi_echo_flash import main as create_seq\n", "from mrseq.utils import combine_ismrmrd_files\n", "from mrseq.utils import sys_defaults" ] }, { "cell_type": "markdown", - "id": "3", + "id": "5", "metadata": {}, "source": [ "### Settings\n", @@ -62,7 +83,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4", + "id": "6", "metadata": {}, "outputs": [], "source": [ @@ -74,7 +95,7 @@ }, { "cell_type": "markdown", - "id": "5", + "id": "7", "metadata": {}, "source": [ "### Create the digital phantom\n", @@ -85,7 +106,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6", + "id": "8", "metadata": {}, "outputs": [], "source": [ @@ -100,7 +121,7 @@ }, { "cell_type": "markdown", - "id": "7", + "id": "9", "metadata": {}, "source": [ "### Create the $T2^*$ mapping sequence\n", @@ -111,7 +132,7 @@ { "cell_type": "code", "execution_count": null, - "id": "8", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -129,7 +150,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "### Simulate the sequence\n", @@ -140,7 +161,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -152,7 +173,7 @@ }, { "cell_type": "markdown", - "id": "11", + "id": "13", "metadata": {}, "source": [ "### Reconstruct the image\n", @@ -162,7 +183,7 @@ { "cell_type": "code", "execution_count": null, - "id": "12", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -180,7 +201,7 @@ }, { "cell_type": "markdown", - "id": "13", + "id": "15", "metadata": {}, "source": [ "We can now visualize the images obtained at the different echo times." @@ -189,7 +210,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "16", "metadata": {}, "outputs": [], "source": [ @@ -204,7 +225,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "17", "metadata": {}, "source": [ "### Estimate $T2^*$\n", @@ -215,7 +236,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -225,7 +246,7 @@ { "cell_type": "code", "execution_count": null, - "id": "17", + "id": "19", "metadata": {}, "outputs": [], "source": [ diff --git a/pyproject.toml b/pyproject.toml index dc236f9..20c6f7a 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -54,8 +54,8 @@ tests = [ "pytest-cov", "pytest-xdist", ] -examples = ["mrzerocore>=0.4.3", "mrpro", "cmap", "jupyter-book<2.0"] -dev = ["mrseq[tests, examples]"] +examples = ["mrzerocore>=0.4.3", "mrpro", "cmap"] +dev = ["mrseq[tests, examples]", "jupyter-book<2.0"] [project.urls] "repository" = "https://github.com/PTB-MR/mrseq" diff --git a/src/mrseq/scripts/__init__.py b/src/mrseq/sequences/__init__.py similarity index 100% rename from src/mrseq/scripts/__init__.py rename to src/mrseq/sequences/__init__.py diff --git a/src/mrseq/scripts/grpe_flash_dixon.py b/src/mrseq/sequences/grpe_flash_dixon.py similarity index 99% rename from src/mrseq/scripts/grpe_flash_dixon.py rename to src/mrseq/sequences/grpe_flash_dixon.py index 2265a91..7b4d694 100644 --- a/src/mrseq/scripts/grpe_flash_dixon.py +++ b/src/mrseq/sequences/grpe_flash_dixon.py @@ -434,6 +434,15 @@ def grpe_flash_dixon_kernel( prot.close() delta_te_array = np.diff(time_to_echoes) + + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_x, fov_y, fov_z]) + seq.set_definition('ReconMatrix', (n_readout, n_rpe_points, n_rpe_points)) + seq.set_definition('SliceThickness', fov_z) + seq.set_definition('TE', [(te or float(min_te)) + idx * float(delta_te_array[0]) for idx in range(n_echoes)]) + seq.set_definition('TR', tr or float(min_tr)) + seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) + return seq, float(min_te), float(min_tr), float(delta_te_array[0]) @@ -537,7 +546,7 @@ def main( if (output_path / Path(filename + '_header.h5')).exists(): (output_path / Path(filename + '_header.h5')).unlink() - seq, min_te, min_tr, delta_te = grpe_flash_dixon_kernel( + seq, _min_te, min_tr, delta_te = grpe_flash_dixon_kernel( system=system, te=te, delta_te=delta_te, @@ -580,14 +589,6 @@ def main( print('\nCreating advanced test report...') print(seq.test_report()) - # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov_x, fov_y, fov_z]) - seq.set_definition('ReconMatrix', (n_readout, n_rpe_points, n_rpe_points)) - seq.set_definition('SliceThickness', fov_z) - seq.set_definition('TE', [(te or min_te) + idx * delta_te for idx in range(n_echoes)]) - seq.set_definition('TR', tr or min_tr) - seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) - # save seq-file to disk print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") write_sequence(seq, str(output_path / filename), create_signature=True, v141_compatibility=v141_compatibility) diff --git a/src/mrseq/scripts/radial_flash.py b/src/mrseq/sequences/radial_flash.py similarity index 99% rename from src/mrseq/scripts/radial_flash.py rename to src/mrseq/sequences/radial_flash.py index 316100c..e590eb2 100644 --- a/src/mrseq/scripts/radial_flash.py +++ b/src/mrseq/sequences/radial_flash.py @@ -279,6 +279,14 @@ def radial_flash_kernel( if mrd_header_file: prot.close() + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness * n_slices]) + seq.set_definition('ReconMatrix', (n_readout, n_readout, 1)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', te or min_te) + seq.set_definition('TR', tr or min_tr) + seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) + return seq, min_te, min_tr @@ -374,7 +382,7 @@ def main( if (output_path / Path(filename + '_header.h5')).exists(): (output_path / Path(filename + '_header.h5')).unlink() - seq, min_te, min_tr = radial_flash_kernel( + seq, _min_te, min_tr = radial_flash_kernel( system=system, te=te, tr=tr, @@ -412,14 +420,6 @@ def main( print('\nCreating advanced test report...') print(seq.test_report()) - # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness * n_slices]) - seq.set_definition('ReconMatrix', (n_readout, n_readout, 1)) - seq.set_definition('SliceThickness', slice_thickness) - seq.set_definition('TE', te or min_te) - seq.set_definition('TR', tr or min_tr) - seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) - # save seq-file to disk print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") write_sequence(seq, str(output_path / filename), create_signature=True, v141_compatibility=v141_compatibility) diff --git a/src/mrseq/scripts/spiral_flash.py b/src/mrseq/sequences/spiral_flash.py similarity index 99% rename from src/mrseq/scripts/spiral_flash.py rename to src/mrseq/sequences/spiral_flash.py index 8a03c9e..219adbd 100644 --- a/src/mrseq/scripts/spiral_flash.py +++ b/src/mrseq/sequences/spiral_flash.py @@ -247,6 +247,13 @@ def spiral_flash_kernel( if mrd_header_file: prot.close() + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness * n_slices]) + seq.set_definition('ReconMatrix', (n_readout, n_readout, 1)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', te or min_te) + seq.set_definition('TR', tr or min_tr) + return seq, min_te, min_tr @@ -336,7 +343,7 @@ def main( if (output_path / Path(filename + '_header.h5')).exists(): (output_path / Path(filename + '_header.h5')).unlink() - seq, min_te, min_tr = spiral_flash_kernel( + seq, _min_te, min_tr = spiral_flash_kernel( system=system, te=te, tr=tr, @@ -372,13 +379,6 @@ def main( print('\nCreating advanced test report...') print(seq.test_report()) - # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness * n_slices]) - seq.set_definition('ReconMatrix', (n_readout, n_readout, 1)) - seq.set_definition('SliceThickness', slice_thickness) - seq.set_definition('TE', te or min_te) - seq.set_definition('TR', tr or min_tr) - # save seq-file to disk print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") write_sequence(seq, str(output_path / filename), create_signature=True, v141_compatibility=v141_compatibility) diff --git a/src/mrseq/scripts/t1_inv_rec_gre_single_line.py b/src/mrseq/sequences/t1_inv_rec_gre_single_line.py similarity index 99% rename from src/mrseq/scripts/t1_inv_rec_gre_single_line.py rename to src/mrseq/sequences/t1_inv_rec_gre_single_line.py index 19402b6..2c5187d 100644 --- a/src/mrseq/scripts/t1_inv_rec_gre_single_line.py +++ b/src/mrseq/sequences/t1_inv_rec_gre_single_line.py @@ -197,6 +197,14 @@ def t1_inv_rec_gre_single_line_kernel( seq.add_block(pp.make_delay(tr_delay)) + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) + seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, 1)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', te or min_te) + seq.set_definition('TR', tr) + seq.set_definition('TI', inversion_times.tolist()) + return seq, time_to_first_tr_block, min_te @@ -273,7 +281,7 @@ def main( rf_bwt = 4 # bandwidth-time product of rf excitation pulse [Hz*s] rf_apodization = 0.5 # apodization factor of rf excitation pulse - seq, time_to_first_tr_block, min_te = t1_inv_rec_gre_single_line_kernel( + seq, time_to_first_tr_block, _min_te = t1_inv_rec_gre_single_line_kernel( system=system, inversion_times=inversion_times, te=te, @@ -312,14 +320,6 @@ def main( f'{Path(__file__).stem}_{int(fov_xy * 1000)}fov_{n_readout}nx_{n_phase_encoding}ny_{len(inversion_times)}TIs' ) - # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) - seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, 1)) - seq.set_definition('SliceThickness', slice_thickness) - seq.set_definition('TE', te or min_te) - seq.set_definition('TR', tr) - seq.set_definition('TI', inversion_times.tolist()) - # save seq-file to disk output_path = Path.cwd() / 'output' output_path.mkdir(parents=True, exist_ok=True) diff --git a/src/mrseq/scripts/t1_inv_rec_se_single_line.py b/src/mrseq/sequences/t1_inv_rec_se_single_line.py similarity index 99% rename from src/mrseq/scripts/t1_inv_rec_se_single_line.py rename to src/mrseq/sequences/t1_inv_rec_se_single_line.py index 290aaab..eabef5c 100644 --- a/src/mrseq/scripts/t1_inv_rec_se_single_line.py +++ b/src/mrseq/sequences/t1_inv_rec_se_single_line.py @@ -258,6 +258,14 @@ def t1_inv_rec_se_single_line_kernel( seq.add_block(pp.make_delay(tr_delay)) + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) + seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, 1)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', te or min_te) + seq.set_definition('TR', tr) + seq.set_definition('TI', inversion_times.tolist()) + return seq, time_to_first_tr_block, min_te @@ -344,7 +352,7 @@ def main( rf180_bwt = 4 # bandwidth-time product of rf refocusing pulse [Hz*s] rf180_apodization = 0.5 # apodization factor of rf refocusing pulse - seq, time_to_first_tr_block, min_te = t1_inv_rec_se_single_line_kernel( + seq, time_to_first_tr_block, _min_te = t1_inv_rec_se_single_line_kernel( system=system, inversion_times=inversion_times, te=te, @@ -389,14 +397,6 @@ def main( f'{Path(__file__).stem}_{int(fov_xy * 1000)}fov_{n_readout}nx_{n_phase_encoding}ny_{len(inversion_times)}TIs' ) - # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) - seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, 1)) - seq.set_definition('SliceThickness', slice_thickness) - seq.set_definition('TE', te or min_te) - seq.set_definition('TR', tr) - seq.set_definition('TI', inversion_times.tolist()) - # save seq-file to disk output_path = Path.cwd() / 'output' output_path.mkdir(parents=True, exist_ok=True) diff --git a/src/mrseq/scripts/t1_molli_bssfp.py b/src/mrseq/sequences/t1_molli_bssfp.py similarity index 99% rename from src/mrseq/scripts/t1_molli_bssfp.py rename to src/mrseq/sequences/t1_molli_bssfp.py index 5554778..45979a1 100644 --- a/src/mrseq/scripts/t1_molli_bssfp.py +++ b/src/mrseq/sequences/t1_molli_bssfp.py @@ -6,11 +6,11 @@ import pypulseq as pp from mrseq.preparations import add_t1_inv_prep +from mrseq.utils import cartesian_phase_encoding from mrseq.utils import find_gx_flat_time_on_adc_raster from mrseq.utils import round_to_raster from mrseq.utils import sys_defaults from mrseq.utils import write_sequence -from mrseq.utils.trajectory import cartesian_phase_encoding def t1_molli_bssfp_kernel( @@ -309,6 +309,15 @@ def t1_molli_bssfp_kernel( # add variable part of trigger delay (soft delay) seq.add_block(trig_soft_delay) + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) + seq.set_definition('ReconMatrix', (n_readout, n_readout, 1)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', te or min_te) + seq.set_definition('TR', tr or min_tr) + seq.set_definition('TI', inversion_times.tolist()) + seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) + return seq, min_te, min_tr @@ -401,7 +410,7 @@ def main( output_path = Path.cwd() / 'output' output_path.mkdir(parents=True, exist_ok=True) - seq, min_te, min_tr = t1_molli_bssfp_kernel( + seq, _min_te, _min_tr = t1_molli_bssfp_kernel( system=system, te=te, tr=tr, @@ -437,15 +446,6 @@ def main( print('\nCreating advanced test report...') print(seq.test_report()) - # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) - seq.set_definition('ReconMatrix', (n_readout, n_readout, 1)) - seq.set_definition('SliceThickness', slice_thickness) - seq.set_definition('TE', te or min_te) - seq.set_definition('TR', tr or min_tr) - seq.set_definition('TI', inversion_times.tolist()) - seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) - # save seq-file to disk print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") write_sequence(seq, str(output_path / filename), create_signature=True, v141_compatibility=v141_compatibility) diff --git a/src/mrseq/scripts/t1_t2_spiral_cmrf.py b/src/mrseq/sequences/t1_t2_spiral_cmrf.py similarity index 98% rename from src/mrseq/scripts/t1_t2_spiral_cmrf.py rename to src/mrseq/sequences/t1_t2_spiral_cmrf.py index 91cf143..8cc8a21 100644 --- a/src/mrseq/scripts/t1_t2_spiral_cmrf.py +++ b/src/mrseq/sequences/t1_t2_spiral_cmrf.py @@ -329,6 +329,18 @@ def t1_t2_spiral_cmrf_kernel( if mrd_header_file: prot.close() + # write all important parameters into the seq-file definitions + seq.set_definition('Name', 'cMRF_spiral') + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) + seq.set_definition('TE', min_te) + seq.set_definition('TI', time_since_inversion) + seq.set_definition('TR', tr) + seq.set_definition('t2prep_te', [0, 0, t2_prep_echo_times[0], t2_prep_echo_times[1], t2_prep_echo_times[2]]) + seq.set_definition('t1prep_ti', [time_since_inversion, 0, 0, 0, 0]) + seq.set_definition('slice_thickness', slice_thickness) + seq.set_definition('sampling_scheme', 'spiral') + seq.set_definition('number_of_readouts', int(n_readout)) + return seq, time_since_inversion, min_te @@ -408,7 +420,7 @@ def main( if (output_path / Path(filename + '_header.h5')).exists(): (output_path / Path(filename + '_header.h5')).unlink() - seq, inversion_time, te = t1_t2_spiral_cmrf_kernel( + seq, _inversion_time, _te = t1_t2_spiral_cmrf_kernel( system=system, t2_prep_echo_times=t2_prep_echo_times, tr=tr, @@ -441,18 +453,6 @@ def main( print('\nCreating advanced test report...') print(seq.test_report()) - # write all important parameters into the seq-file definitions - seq.set_definition('Name', 'cMRF_spiral') - seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) - seq.set_definition('TE', te) - seq.set_definition('TI', inversion_time) - seq.set_definition('TR', tr) - seq.set_definition('t2prep_te', [0, 0, t2_prep_echo_times[0], t2_prep_echo_times[1], t2_prep_echo_times[2]]) - seq.set_definition('t1prep_ti', [inversion_time, 0, 0, 0, 0]) - seq.set_definition('slice_thickness', slice_thickness) - seq.set_definition('sampling_scheme', 'spiral') - seq.set_definition('number_of_readouts', int(n_readout)) - # save seq-file to disk output_path = Path.cwd() / 'output' output_path.mkdir(parents=True, exist_ok=True) diff --git a/src/mrseq/scripts/t1rho_se_single_line.py b/src/mrseq/sequences/t1rho_se_single_line.py similarity index 99% rename from src/mrseq/scripts/t1rho_se_single_line.py rename to src/mrseq/sequences/t1rho_se_single_line.py index d396c84..6bf69d4 100644 --- a/src/mrseq/scripts/t1rho_se_single_line.py +++ b/src/mrseq/sequences/t1rho_se_single_line.py @@ -252,6 +252,14 @@ def t1rho_se_single_line_kernel( seq.add_block(pp.make_delay(tr_delay)) + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) + seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, 1)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', te or min_te) + seq.set_definition('TR', tr) + seq.set_definition('TSL', spin_lock_times.tolist()) + return seq, time_to_first_tr_block, min_te @@ -336,7 +344,7 @@ def main( rf180_bwt = 4 # bandwidth-time product of rf refocusing pulse [Hz*s] rf180_apodization = 0.5 # apodization factor of rf refocusing pulse - seq, time_to_first_tr_block, min_te = t1rho_se_single_line_kernel( + seq, time_to_first_tr_block, _min_te = t1rho_se_single_line_kernel( system=system, spin_lock_times=spin_lock_times, te=te, @@ -381,14 +389,6 @@ def main( f'{Path(__file__).stem}_{int(fov_xy * 1000)}fov_{n_readout}nx_{n_phase_encoding}ny_{len(spin_lock_times)}TIs' ) - # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) - seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, 1)) - seq.set_definition('SliceThickness', slice_thickness) - seq.set_definition('TE', te or min_te) - seq.set_definition('TR', tr) - seq.set_definition('TSL', spin_lock_times.tolist()) - # save seq-file to disk output_path = Path.cwd() / 'output' output_path.mkdir(parents=True, exist_ok=True) diff --git a/src/mrseq/scripts/t2_multi_echo_se_single_line.py b/src/mrseq/sequences/t2_multi_echo_se_single_line.py similarity index 100% rename from src/mrseq/scripts/t2_multi_echo_se_single_line.py rename to src/mrseq/sequences/t2_multi_echo_se_single_line.py index 835ed35..b3321fa 100644 --- a/src/mrseq/scripts/t2_multi_echo_se_single_line.py +++ b/src/mrseq/sequences/t2_multi_echo_se_single_line.py @@ -206,6 +206,13 @@ def t2_multi_echo_se_single_line_kernel( seq.add_block(pp.make_delay(tr_delay)) + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) + seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, 1)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', echo_times.tolist()) + seq.set_definition('TR', tr) + return seq, min_tr_first_echo_block @@ -326,13 +333,6 @@ def main( # define sequence filename filename = f'{Path(__file__).stem}_{int(fov_xy * 1000)}fov_{n_readout}nx_{n_phase_encoding}ny_{len(echo_times)}TEs' - # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) - seq.set_definition('ReconMatrix', (n_readout, n_phase_encoding, 1)) - seq.set_definition('SliceThickness', slice_thickness) - seq.set_definition('TE', echo_times.tolist()) - seq.set_definition('TR', tr) - # save seq-file to disk output_path = Path.cwd() / 'output' output_path.mkdir(parents=True, exist_ok=True) diff --git a/src/mrseq/scripts/t2_t2prep_flash.py b/src/mrseq/sequences/t2_t2prep_flash.py similarity index 99% rename from src/mrseq/scripts/t2_t2prep_flash.py rename to src/mrseq/sequences/t2_t2prep_flash.py index 8f69521..a22789c 100644 --- a/src/mrseq/scripts/t2_t2prep_flash.py +++ b/src/mrseq/sequences/t2_t2prep_flash.py @@ -6,11 +6,11 @@ import pypulseq as pp from mrseq.preparations.t2_prep import add_t2_prep +from mrseq.utils import cartesian_phase_encoding from mrseq.utils import find_gx_flat_time_on_adc_raster from mrseq.utils import round_to_raster from mrseq.utils import sys_defaults from mrseq.utils import write_sequence -from mrseq.utils.trajectory import cartesian_phase_encoding def t2_t2prep_flash_kernel( @@ -299,6 +299,14 @@ def t2_t2prep_flash_kernel( # add variable part of trigger delay (soft delay) seq.add_block(trig_soft_delay) + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) + seq.set_definition('ReconMatrix', (n_readout, n_readout, 1)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', te or min_te) + seq.set_definition('TR', tr or min_tr) + seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) + return seq, min_te, min_tr @@ -399,7 +407,7 @@ def main( output_path = Path.cwd() / 'output' output_path.mkdir(parents=True, exist_ok=True) - seq, min_te, min_tr = t2_t2prep_flash_kernel( + seq, _min_te, _min_tr = t2_t2prep_flash_kernel( system=system, te=te, tr=tr, @@ -438,14 +446,6 @@ def main( print('\nCreating advanced test report...') print(seq.test_report()) - # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) - seq.set_definition('ReconMatrix', (n_readout, n_readout, 1)) - seq.set_definition('SliceThickness', slice_thickness) - seq.set_definition('TE', te or min_te) - seq.set_definition('TR', tr or min_tr) - seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) - # save seq-file to disk print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") write_sequence(seq, str(output_path / filename), create_signature=True, v141_compatibility=v141_compatibility) diff --git a/src/mrseq/scripts/t2star_multi_echo_flash.py b/src/mrseq/sequences/t2star_multi_echo_flash.py similarity index 98% rename from src/mrseq/scripts/t2star_multi_echo_flash.py rename to src/mrseq/sequences/t2star_multi_echo_flash.py index e73638e..55a0719 100644 --- a/src/mrseq/scripts/t2star_multi_echo_flash.py +++ b/src/mrseq/sequences/t2star_multi_echo_flash.py @@ -6,6 +6,8 @@ import numpy as np import pypulseq as pp +from mrseq.utils import MultiEchoAcquisition +from mrseq.utils import cartesian_phase_encoding from mrseq.utils import find_gx_flat_time_on_adc_raster from mrseq.utils import round_to_raster from mrseq.utils import sys_defaults @@ -14,8 +16,6 @@ from mrseq.utils.ismrmrd import Limits from mrseq.utils.ismrmrd import MatrixSize from mrseq.utils.ismrmrd import create_header -from mrseq.utils.trajectory import MultiEchoAcquisition -from mrseq.utils.trajectory import cartesian_phase_encoding def t2star_multi_echo_flash_kernel( @@ -317,6 +317,15 @@ def t2star_multi_echo_flash_kernel( prot.close() delta_te_array = np.diff(time_to_echoes) + + # write all required parameters in the seq-file header/definitions + seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) + seq.set_definition('ReconMatrix', (n_readout, n_readout, 1)) + seq.set_definition('SliceThickness', slice_thickness) + seq.set_definition('TE', [(te or float(min_te)) + idx * float(delta_te_array[0]) for idx in range(n_echoes)]) + seq.set_definition('TR', tr or float(min_tr)) + seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) + return seq, float(min_te), float(min_tr), float(delta_te_array[0]) @@ -420,7 +429,7 @@ def main( output_path = Path.cwd() / 'output' output_path.mkdir(parents=True, exist_ok=True) - seq, min_te, min_tr, delta_te = t2star_multi_echo_flash_kernel( + seq, _min_te, _min_tr, delta_te = t2star_multi_echo_flash_kernel( system=system, te=te, delta_te=delta_te, @@ -461,14 +470,6 @@ def main( print('\nCreating advanced test report...') print(seq.test_report()) - # write all required parameters in the seq-file header/definitions - seq.set_definition('FOV', [fov_xy, fov_xy, slice_thickness]) - seq.set_definition('ReconMatrix', (n_readout, n_readout, 1)) - seq.set_definition('SliceThickness', slice_thickness) - seq.set_definition('TE', [(te or min_te) + idx * delta_te for idx in range(n_echoes)]) - seq.set_definition('TR', tr or min_tr) - seq.set_definition('ReadoutOversamplingFactor', readout_oversampling) - # save seq-file to disk print(f"\nSaving sequence file '{filename}.seq' into folder '{output_path}'.") write_sequence(seq, str(output_path / filename), create_signature=True, v141_compatibility=v141_compatibility) diff --git a/src/mrseq/utils/MultiEchoAcquisition.py b/src/mrseq/utils/MultiEchoAcquisition.py new file mode 100644 index 0000000..53d0bab --- /dev/null +++ b/src/mrseq/utils/MultiEchoAcquisition.py @@ -0,0 +1,209 @@ +"""Basic functionality for trajectory calculation.""" + +from typing import Literal + +import numpy as np +import pypulseq as pp + +from mrseq.utils import find_gx_flat_time_on_adc_raster +from mrseq.utils import round_to_raster +from mrseq.utils import sys_defaults + + +class MultiEchoAcquisition: + """ + Multi-echo gradient echo acquisition. + + Attributes + ---------- + system + PyPulseq system limits object. + n_readout_post_echo + Number of readout points after echo. + n_readout_pre_echo + Number of readout points before echo. + n_readout_with_partial_echo + Total number of readout points with partial echo. + te_delay + Additional delay after readout gradient gx to achieve desired delta echo time. + adc + ADC event object. + gx + Readout gradient object. + gx_pre + Pre-winder gradient object. + gx_post + Re-winder gradient object. + gx_between + Gradient between echoes. + """ + + def __init__( + self, + system: pp.Opts | None = None, + delta_te: float | None = None, + fov: float = 0.256, + n_readout: int = 128, + readout_oversampling: float = 2.0, + partial_echo_factor: float = 0.7, + gx_flat_time: float = 2.0e-3, + gx_pre_duration: float = 0.8e-3, + ): + """ + Initialize the MultiEchoAcquisition class and compute all required attributes. + + Parameters + ---------- + system + PyPulseq system limits object. + delta_te + Desired echo spacing (in seconds). Minimum echo spacing is used if set to None. + fov + Field of view in x direction (in meters). + n_readout + Number of frequency encoding steps. + readout_oversampling + Readout oversampling factor. + partial_echo_factor + Partial echo factor. + gx_flat_time + Flat time of the readout gradient. + gx_pre_duration + Duration of readout pre-winder gradient. + """ + # set system to default if not provided + self._system = sys_defaults if system is None else system + + delta_k = 1 / (fov * readout_oversampling) + self._n_readout_post_echo = int(n_readout * readout_oversampling / 2 - 1) + self._n_readout_post_echo += np.mod(self._n_readout_post_echo + 1, 2) # make odd + self._n_readout_pre_echo = int( + (n_readout * partial_echo_factor * readout_oversampling) - self._n_readout_post_echo - 1 + ) + self._n_readout_pre_echo += np.mod(self._n_readout_pre_echo, 2) # make even + + self._n_readout_with_partial_echo = self._n_readout_pre_echo + 1 + self._n_readout_post_echo + gx_flat_area = self._n_readout_with_partial_echo * delta_k + + # adc dwell time has to be on adc raster and gx flat time on gradient raster + self._gx_flat_time, _ = find_gx_flat_time_on_adc_raster( + self._n_readout_with_partial_echo, + gx_flat_time / self._n_readout_with_partial_echo, + self._system.grad_raster_time, + self._system.adc_raster_time, + ) + + self._gx = pp.make_trapezoid( + channel='x', flat_area=gx_flat_area, flat_time=self._gx_flat_time, system=self._system + ) + + self._adc = pp.make_adc( + num_samples=self._n_readout_with_partial_echo, + duration=self._gx.flat_time, + delay=self._gx.rise_time, + system=self._system, + ) + + self._gx_pre = pp.make_trapezoid( + channel='x', + area=-(self._gx.amplitude * self._gx.rise_time / 2 + delta_k * (self._n_readout_pre_echo + 0.5)), + duration=gx_pre_duration * partial_echo_factor, + system=self._system, + ) + self._gx_post = pp.make_trapezoid( + channel='x', + area=-(self._gx.amplitude * self._gx.fall_time / 2 + delta_k * (self._n_readout_post_echo + 0.5)), + duration=gx_pre_duration, + system=self._system, + ) + + self._gx_between = pp.make_trapezoid( + channel='x', + area=self._gx_pre.area - self._gx_post.area, + duration=gx_pre_duration, + system=self._system, + ) + + min_delta_te = pp.calc_duration(self._gx) + pp.calc_duration(self._gx_between) + if delta_te is None: + self._delta_te_delay = 0.0 + else: + self._delta_te_delay = round_to_raster(delta_te - min_delta_te, self._system.block_duration_raster) + if self._delta_te_delay < 0: + raise ValueError( + f'TE must be larger than {min_delta_te * 1000:.3f} ms. Current value is {delta_te * 1000:.3f} ms.' + ) + + def add_to_seq(self, seq: pp.Sequence, n_echoes: int, polarity: Literal['positive', 'negative'] = 'positive'): + """Add all gradients and adc to sequence. + + Parameters + ---------- + seq + PyPulseq Sequence object. + n_echoes + Number of echoes + polarity + Polarity of first readout gradient + + Returns + ------- + seq + PyPulseq Sequence object. + time_to_echoes + Time from beginning of sequence to echoes. + """ + # readout pre-winder + seq.add_block(self._gx_pre if polarity == 'positive' else pp.scale_grad(self._gx_pre, -1)) + + # add readout gradients and ADCs + seq, time_to_echoes = self.add_to_seq_without_pre_post_gradient(seq, n_echoes, polarity) + + # readout re-winder + seq.add_block(self._gx_post if polarity == 'positive' else pp.scale_grad(self._gx_post, -1)) + + return seq, time_to_echoes + + def add_to_seq_without_pre_post_gradient( + self, seq: pp.Sequence, n_echoes: int, polarity: Literal['positive', 'negative'] = 'positive' + ): + """Add readout gradients without pre- and re-winder gradients. + + Often the pre- and re-winder gradients are played out at the same time as phase encoding gradients or spoiler + gradients. + + Parameters + ---------- + seq + PyPulseq Sequence object. + n_echoes + Number of echoes + polarity + Polarity of first readout gradient + + Returns + ------- + seq + PyPulseq Sequence object. + time_to_echoes + Time from beginning of sequence to echoes. + """ + # add readout gradient and ADC + time_to_echoes = [] + for echo_ in range(n_echoes): + start_of_current_gx = sum(seq.block_durations.values()) + gx_sign = (-1) ** echo_ if polarity == 'positive' else (-1) ** (echo_ + 1) + labels = [] + labels.append(pp.make_label(label='REV', type='SET', value=gx_sign == -1)) + labels.append(pp.make_label(label='ECO', type='SET', value=echo_)) + seq.add_block(pp.scale_grad(self._gx, gx_sign), self._adc, *labels) + time_to_echoes.append( + start_of_current_gx + self._adc.delay + self._n_readout_pre_echo * self._adc.dwell + self._adc.dwell / 2 + ) + start_of_current_gx = sum(seq.block_durations.values()) + if echo_ < n_echoes - 1: + if self._delta_te_delay > 0: + seq.add_block(pp.make_delay(self._delta_te_delay)) + seq.add_block(pp.scale_grad(self._gx_between, -gx_sign)) + + return seq, time_to_echoes diff --git a/src/mrseq/utils/__init__.py b/src/mrseq/utils/__init__.py index 31bd0a0..25e483b 100644 --- a/src/mrseq/utils/__init__.py +++ b/src/mrseq/utils/__init__.py @@ -1,5 +1,7 @@ from mrseq.utils.sequence_helper import find_gx_flat_time_on_adc_raster, round_to_raster, write_sequence from mrseq.utils.system_defaults import sys_defaults -from mrseq.utils.trajectory import cartesian_phase_encoding, spiral_acquisition, MultiEchoAcquisition from mrseq.utils.vds import variable_density_spiral_trajectory +from mrseq.utils.spiral_sampling import spiral_acquisition +from mrseq.utils.cartesian_sampling import cartesian_phase_encoding +from mrseq.utils.MultiEchoAcquisition import MultiEchoAcquisition from mrseq.utils.ismrmrd import create_header, combine_ismrmrd_files diff --git a/src/mrseq/utils/cartesian_sampling.py b/src/mrseq/utils/cartesian_sampling.py new file mode 100644 index 0000000..a26c62d --- /dev/null +++ b/src/mrseq/utils/cartesian_sampling.py @@ -0,0 +1,82 @@ +"""Cartesian trajectory calculation.""" + +import warnings +from typing import Literal + +import numpy as np + + +def cartesian_phase_encoding( + n_phase_encoding: int, + acceleration: int = 1, + n_fully_sampled_center: int = 0, + sampling_order: Literal['linear', 'low_high', 'high_low', 'random'] = 'linear', + n_phase_encoding_per_shot: int | None = None, +) -> tuple[np.ndarray, np.ndarray]: + """Calculate Cartesian sampling trajectory. + + Parameters + ---------- + n_phase_encoding + number of phase encoding points before undersampling + acceleration + undersampling factor + n_fully_sampled_center + number of phsae encoding points in the fully sampled center. This will reduce the overall undersampling factor. + sampling_order + order how phase encoding points are sampled + n_phase_encoding_per_shot + used to ensure that all phase encoding points can be acquired in an integer number of shots. If None, this + parameter is ignored, i.e. equal to n_phase_encoding_per_shot = 1 + """ + if n_fully_sampled_center > n_phase_encoding: + warnings.warn( + 'Number of phase encoding steps in the fully sampled center will be reduced to the total number of phase ' + + 'encoding steps.', + stacklevel=2, + ) + n_fully_sampled_center = n_phase_encoding + + if sampling_order == 'random': + # Linear order of a fully sampled kpe dimension. Undersampling is done later. + kpe = np.arange(0, n_phase_encoding) + else: + # Always include k-space center and more points on the negative side of k-space + kpe_pos = np.arange(0, n_phase_encoding // 2, acceleration) + kpe_neg = -np.arange(acceleration, n_phase_encoding // 2 + 1, acceleration) + kpe = np.concatenate((kpe_neg, kpe_pos), axis=0) + + # Ensure fully sampled center + kpe_fully_sampled_center = np.arange( + -n_fully_sampled_center // 2, -n_fully_sampled_center // 2 + n_fully_sampled_center + ) + kpe = np.unique(np.concatenate((kpe, kpe_fully_sampled_center))) + + # Always acquire more to ensure desired resolution + if n_phase_encoding_per_shot and sampling_order != 'random': + kpe_extended = np.arange(-n_phase_encoding, n_phase_encoding) + kpe_extended = kpe_extended[np.argsort(np.abs(kpe_extended), kind='stable')] + kidx = 0 + while np.mod(len(kpe), n_phase_encoding_per_shot) > 0: + kpe = np.unique(np.concatenate((kpe, (kpe_extended[kidx],)))) + kidx += 1 + + # Different temporal orders of phase encoding points + if sampling_order == 'random': + perm = np.random.permutation(kpe) + npe = len(perm) // acceleration + if n_phase_encoding_per_shot: + npe += n_phase_encoding_per_shot - np.mod(npe, n_phase_encoding_per_shot) + kpe = kpe[perm[:npe]] + elif sampling_order == 'linear': + kpe = np.sort(kpe) + elif sampling_order == 'low_high': + idx = np.argsort(np.abs(kpe), kind='stable') + kpe = kpe[idx] + elif sampling_order == 'high_low': + idx = np.argsort(-np.abs(kpe), kind='stable') + kpe = kpe[idx] + else: + raise ValueError(f'sampling order {sampling_order} not supported.') + + return kpe, kpe_fully_sampled_center diff --git a/src/mrseq/utils/spiral_sampling.py b/src/mrseq/utils/spiral_sampling.py new file mode 100644 index 0000000..584858e --- /dev/null +++ b/src/mrseq/utils/spiral_sampling.py @@ -0,0 +1,255 @@ +"""Basic functionality for trajectory calculation.""" + +from typing import Any +from typing import Literal + +import numpy as np +import pypulseq as pp + +from mrseq.utils import variable_density_spiral_trajectory + + +def undersampled_variable_density_spiral( + system: pp.Opts, n_readout: int, fov: float, undersampling_factor: float +) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, int, float, float]: + """Create undersampled variable density spiral. + + The distribution of the k-space points of a spiral trajectory are restricted by the maximum gradient amplitude and + slew rate. This makes an analytic solution for a given undersampling factor challenging. Here we use an iterative + approach in order to achieve a variable density spiral with a certain number of readout samplings and undersampling + factor. + + During the iterative search, the undersampling for the edge of k-space is increased. If this is not enough, then we + also start to increase the undersampling in the k-space center. The field-of-view varies linearly bewtween the + k-space center and k-space edge. + + If the undersampling factor is to high, it might not be possible to find a suitable solution. + + Parameters + ---------- + system + PyPulseq system limits object. + n_readout + Number of readout points per spiral. + fov + Field of view (in meters). + undersampling_factor + Undersampling factor of spiral trajectory + + Returns + ------- + tuple(np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, int, float, float) + - k-space trajectory (traj) + - Gradient waveform (grad) + - Slew rate (slew) + - Time points for the trajectory (timing) + - Radius values (radius) + - Angular positions (theta) + - Number of spiral arms (n_spirals) + - Scaling of the field-of-view in the k-space center + - Scaling of the field-of-view in the k-space edge + + """ + # calculate single spiral trajectory + n_k0 = np.inf + fov_scaling_center = 1.0 + fov_scaling_edge = 1.0 + n_spirals = int(np.round(n_readout / undersampling_factor)) + while n_k0 > n_readout: + fov_coefficients = [fov * fov_scaling_center, -fov * (1 - fov_scaling_edge)] + + try: + traj, grad, slew, timing, radius, theta = variable_density_spiral_trajectory( + system=system, + sampling_period=system.grad_raster_time, + n_interleaves=n_spirals, + fov_coefficients=fov_coefficients, + max_kspace_radius=0.5 / fov * n_readout, + ) + n_k0 = len(grad) + fov_scaling_edge *= 0.95 + except ValueError: + # It is not possible to achieve the desired undersampling factor with the given system limits while keeping + # the full field-of-view in the k-space center. Reduce the field-of-view and try again. + n_k0 = np.inf + fov_scaling_center *= 0.95 + fov_scaling_edge = fov_scaling_center + + if fov_scaling_center < 0.1: + raise ValueError('Cannot find a suitable trajectory.') + + return traj, grad, slew, timing, radius, theta, n_spirals, fov_coefficients[0] / fov, fov_coefficients[1] / fov + 1 + + +def spiral_acquisition( + system: pp.Opts, + n_readout: int, + fov: float, + undersampling_factor: float, + readout_oversampling: Literal[1, 2, 4], + n_spirals: int | None, + max_pre_duration: float, + spiral_type: Literal['out', 'in-out'], +) -> tuple[list[Any], list[Any], Any, np.ndarray, float]: + """Generate a spiral acquisition sequence. + + Parameters + ---------- + system + PyPulseq system limits object. + n_readout + Number of readout points per spiral. + fov + Field of view (in meters). + undersampling_factor + Undersampling factor. + readout_oversampling + Oversampling factor for the readout trajectory. + n_spirals + Number of spirals to generate. If set to None, this value will be set based on the undersampling factor. + max_pre_duration : float + Maximum duration for pre-winder gradients (in seconds). + spiral_type + Type of spiral acquisition. 'out' for outward spirals, 'in-out' for spirals turning in and then out. + + Returns + ------- + gx_combined + List of combined gradient objects for the x-channel. + gy_combined + List of combined gradient objects for the y-channel. + adc + PyPulseq ADC object for the acquisition. + trajectory + K-space trajectory. + time_to_echo + Time to echo from beginning of gradients (in seconds). + """ + # calculate single spiral trajectory + traj, grad, _s, _timing, _r, _theta, n_spirals_undersampling, fov_scaling_center, fov_scaling_edge = ( + undersampled_variable_density_spiral(system, n_readout, fov, undersampling_factor) + ) + n_spirals = n_spirals_undersampling if n_spirals is None else n_spirals + print( + f'Target undersampling: {undersampling_factor} - ', + f'achieved undersampling: {n_readout**2 / (len(traj) * n_spirals_undersampling):.2f}', + f'FOV: {fov * fov_scaling_center:.3f} (k-sapce center) - {fov * fov_scaling_edge:.3f} (k-space edge)', + ) + + delta_angle = 2 * np.pi / n_spirals + n_samples_to_echo = 0.5 + if spiral_type == 'in-out': + n_samples_to_echo = len(grad) + grad = np.concatenate((-np.asarray(grad * np.exp(1j * np.pi))[::-1], grad)) + traj = np.concatenate((np.asarray(traj * np.exp(1j * np.pi))[::-1], traj)) + delta_angle = delta_angle / 2 + + # calculate ADC + n_readout_with_oversampling = len(grad) * readout_oversampling + adc_dwell_time = system.grad_raster_time / readout_oversampling + adc = pp.make_adc( + num_samples=n_readout_with_oversampling, dwell=adc_dwell_time, system=system, delay=system.adc_dead_time + ) + traj = np.interp( + np.linspace(0.5 / readout_oversampling, len(grad) - 0.5 / readout_oversampling, n_readout_with_oversampling), + np.linspace(0.5, len(grad) - 0.5, len(grad)), + traj, + ) + + print(f'Receiver bandwidth: {int(1.0 / (adc_dwell_time * n_readout_with_oversampling))} Hz/pixel') + + # Create gradient values and trajectory for different spirals + grad_list = [grad * np.exp(1j * delta_angle * idx) for idx in np.arange(n_spirals)] + traj_list = [traj * np.exp(1j * delta_angle * idx) for idx in np.arange(n_spirals)] + + # Create gradient objects + gx = [pp.make_arbitrary_grad(channel='x', waveform=g.real, delay=adc.delay, system=system) for g in grad_list] + gy = [pp.make_arbitrary_grad(channel='y', waveform=g.imag, delay=adc.delay, system=system) for g in grad_list] + + # Calculate pre- and re-winder gradients + gx_rew, gx_pre, gy_rew, gy_pre = [], [], [], [] + for gx_, gy_ in zip(gx, gy, strict=True): + gx_rew.append( + pp.make_extended_trapezoid_area( + area=-gx_.area if spiral_type == 'out' else -gx_.area / 2, + channel='x', + grad_start=gx_.last, + grad_end=0, + system=system, + convert_to_arbitrary=True, + )[0] + ) + gy_rew.append( + pp.make_extended_trapezoid_area( + area=-gy_.area if spiral_type == 'out' else -gy_.area / 2, + channel='y', + grad_start=gy_.last, + grad_end=0, + system=system, + convert_to_arbitrary=True, + )[0] + ) + + if spiral_type == 'in-out': + gx_pre.append( + pp.make_extended_trapezoid_area( + area=-gx_.area / 2, + channel='x', + grad_start=0, + grad_end=gx_.first, + system=system, + convert_to_arbitrary=True, + )[0] + ) + + gy_pre.append( + pp.make_extended_trapezoid_area( + area=-gy_.area / 2, + channel='y', + grad_start=0, + grad_end=gy_.first, + system=system, + convert_to_arbitrary=True, + )[0] + ) + else: + gx_pre.append(None) + gy_pre.append(None) + + if spiral_type == 'in-out': + adc.delay = max_pre_duration + + for i in range(len(gx_pre)): + gy_pre[i].delay = max_pre_duration - gy_pre[i].shape_dur + gx_pre[i].delay = max_pre_duration - gx_pre[i].shape_dur + else: + max_pre_duration = adc.delay + + def combine_gradients(*grad_objects, channel): + grad_list = [grad for grad in grad_objects if grad is not None] # Remove None + waveform_combined = np.concatenate([grad.waveform for grad in grad_list]) + + return pp.make_arbitrary_grad( + channel=channel, + waveform=waveform_combined, + first=0, + delay=grad_list[0].delay, + last=0, + system=system, + ) + + gx_combined = [ + combine_gradients(gx_pre, gx_in_out, gx_rew, channel='x') + for gx_pre, gx_in_out, gx_rew in zip(gx_pre, gx, gx_rew, strict=True) + ] + gy_combined = [ + combine_gradients(gy_pre, gy_in_out, gy_rew, channel='y') + for gy_pre, gy_in_out, gy_rew in zip(gy_pre, gy, gy_rew, strict=True) + ] + + # times -1 to match pulseq trajectory calculation + trajectory = -np.stack((np.asarray(traj_list).real, np.asarray(traj_list).imag), axis=-1) + + time_to_echo = max_pre_duration + n_samples_to_echo * readout_oversampling * adc.dwell + + return gx_combined, gy_combined, adc, trajectory, time_to_echo diff --git a/src/mrseq/utils/trajectory.py b/src/mrseq/utils/trajectory.py deleted file mode 100644 index c40b5ca..0000000 --- a/src/mrseq/utils/trajectory.py +++ /dev/null @@ -1,534 +0,0 @@ -"""Basic functionality for trajectory calculation.""" - -import warnings -from typing import Any -from typing import Literal - -import numpy as np -import pypulseq as pp - -from mrseq.utils.sequence_helper import find_gx_flat_time_on_adc_raster -from mrseq.utils.sequence_helper import round_to_raster -from mrseq.utils.system_defaults import sys_defaults -from mrseq.utils.vds import variable_density_spiral_trajectory - - -def cartesian_phase_encoding( - n_phase_encoding: int, - acceleration: int = 1, - n_fully_sampled_center: int = 0, - sampling_order: Literal['linear', 'low_high', 'high_low', 'random'] = 'linear', - n_phase_encoding_per_shot: int | None = None, -) -> tuple[np.ndarray, np.ndarray]: - """Calculate Cartesian sampling trajectory. - - Parameters - ---------- - n_phase_encoding - number of phase encoding points before undersampling - acceleration - undersampling factor - n_fully_sampled_center - number of phsae encoding points in the fully sampled center. This will reduce the overall undersampling factor. - sampling_order - order how phase encoding points are sampled - n_phase_encoding_per_shot - used to ensure that all phase encoding points can be acquired in an integer number of shots. If None, this - parameter is ignored, i.e. equal to n_phase_encoding_per_shot = 1 - """ - if n_fully_sampled_center > n_phase_encoding: - warnings.warn( - 'Number of phase encoding steps in the fully sampled center will be reduced to the total number of phase ' - + 'encoding steps.', - stacklevel=2, - ) - n_fully_sampled_center = n_phase_encoding - - if sampling_order == 'random': - # Linear order of a fully sampled kpe dimension. Undersampling is done later. - kpe = np.arange(0, n_phase_encoding) - else: - # Always include k-space center and more points on the negative side of k-space - kpe_pos = np.arange(0, n_phase_encoding // 2, acceleration) - kpe_neg = -np.arange(acceleration, n_phase_encoding // 2 + 1, acceleration) - kpe = np.concatenate((kpe_neg, kpe_pos), axis=0) - - # Ensure fully sampled center - kpe_fully_sampled_center = np.arange( - -n_fully_sampled_center // 2, -n_fully_sampled_center // 2 + n_fully_sampled_center - ) - kpe = np.unique(np.concatenate((kpe, kpe_fully_sampled_center))) - - # Always acquire more to ensure desired resolution - if n_phase_encoding_per_shot and sampling_order != 'random': - kpe_extended = np.arange(-n_phase_encoding, n_phase_encoding) - kpe_extended = kpe_extended[np.argsort(np.abs(kpe_extended), kind='stable')] - kidx = 0 - while np.mod(len(kpe), n_phase_encoding_per_shot) > 0: - kpe = np.unique(np.concatenate((kpe, (kpe_extended[kidx],)))) - kidx += 1 - - # Different temporal orders of phase encoding points - if sampling_order == 'random': - perm = np.random.permutation(kpe) - npe = len(perm) // acceleration - if n_phase_encoding_per_shot: - npe += n_phase_encoding_per_shot - np.mod(npe, n_phase_encoding_per_shot) - kpe = kpe[perm[:npe]] - elif sampling_order == 'linear': - kpe = np.sort(kpe) - elif sampling_order == 'low_high': - idx = np.argsort(np.abs(kpe), kind='stable') - kpe = kpe[idx] - elif sampling_order == 'high_low': - idx = np.argsort(-np.abs(kpe), kind='stable') - kpe = kpe[idx] - else: - raise ValueError(f'sampling order {sampling_order} not supported.') - - return kpe, kpe_fully_sampled_center - - -class MultiEchoAcquisition: - """ - Multi-echo gradient echo acquisition. - - Attributes - ---------- - system - PyPulseq system limits object. - n_readout_post_echo - Number of readout points after echo. - n_readout_pre_echo - Number of readout points before echo. - n_readout_with_partial_echo - Total number of readout points with partial echo. - te_delay - Additional delay after readout gradient gx to achieve desired delta echo time. - adc - ADC event object. - gx - Readout gradient object. - gx_pre - Pre-winder gradient object. - gx_post - Re-winder gradient object. - gx_between - Gradient between echoes. - """ - - def __init__( - self, - system: pp.Opts | None = None, - delta_te: float | None = None, - fov: float = 0.256, - n_readout: int = 128, - readout_oversampling: float = 2.0, - partial_echo_factor: float = 0.7, - gx_flat_time: float = 2.0e-3, - gx_pre_duration: float = 0.8e-3, - ): - """ - Initialize the MultiEchoAcquisition class and compute all required attributes. - - Parameters - ---------- - system - PyPulseq system limits object. - delta_te - Desired echo spacing (in seconds). Minimum echo spacing is used if set to None. - fov - Field of view in x direction (in meters). - n_readout - Number of frequency encoding steps. - readout_oversampling - Readout oversampling factor. - partial_echo_factor - Partial echo factor. - gx_flat_time - Flat time of the readout gradient. - gx_pre_duration - Duration of readout pre-winder gradient. - """ - # set system to default if not provided - self._system = sys_defaults if system is None else system - - delta_k = 1 / (fov * readout_oversampling) - self._n_readout_post_echo = int(n_readout * readout_oversampling / 2 - 1) - self._n_readout_post_echo += np.mod(self._n_readout_post_echo + 1, 2) # make odd - self._n_readout_pre_echo = int( - (n_readout * partial_echo_factor * readout_oversampling) - self._n_readout_post_echo - 1 - ) - self._n_readout_pre_echo += np.mod(self._n_readout_pre_echo, 2) # make even - - self._n_readout_with_partial_echo = self._n_readout_pre_echo + 1 + self._n_readout_post_echo - gx_flat_area = self._n_readout_with_partial_echo * delta_k - - # adc dwell time has to be on adc raster and gx flat time on gradient raster - self._gx_flat_time, _ = find_gx_flat_time_on_adc_raster( - self._n_readout_with_partial_echo, - gx_flat_time / self._n_readout_with_partial_echo, - self._system.grad_raster_time, - self._system.adc_raster_time, - ) - - self._gx = pp.make_trapezoid( - channel='x', flat_area=gx_flat_area, flat_time=self._gx_flat_time, system=self._system - ) - - self._adc = pp.make_adc( - num_samples=self._n_readout_with_partial_echo, - duration=self._gx.flat_time, - delay=self._gx.rise_time, - system=self._system, - ) - - self._gx_pre = pp.make_trapezoid( - channel='x', - area=-(self._gx.amplitude * self._gx.rise_time / 2 + delta_k * (self._n_readout_pre_echo + 0.5)), - duration=gx_pre_duration * partial_echo_factor, - system=self._system, - ) - self._gx_post = pp.make_trapezoid( - channel='x', - area=-(self._gx.amplitude * self._gx.fall_time / 2 + delta_k * (self._n_readout_post_echo + 0.5)), - duration=gx_pre_duration, - system=self._system, - ) - - self._gx_between = pp.make_trapezoid( - channel='x', - area=self._gx_pre.area - self._gx_post.area, - duration=gx_pre_duration, - system=self._system, - ) - - min_delta_te = pp.calc_duration(self._gx) + pp.calc_duration(self._gx_between) - if delta_te is None: - self._delta_te_delay = 0.0 - else: - self._delta_te_delay = round_to_raster(delta_te - min_delta_te, self._system.block_duration_raster) - if self._delta_te_delay < 0: - raise ValueError( - f'TE must be larger than {min_delta_te * 1000:.3f} ms. Current value is {delta_te * 1000:.3f} ms.' - ) - - def add_to_seq(self, seq: pp.Sequence, n_echoes: int, polarity: Literal['positive', 'negative'] = 'positive'): - """Add all gradients and adc to sequence. - - Parameters - ---------- - seq - PyPulseq Sequence object. - n_echoes - Number of echoes - polarity - Polarity of first readout gradient - - Returns - ------- - seq - PyPulseq Sequence object. - time_to_echoes - Time from beginning of sequence to echoes. - """ - # readout pre-winder - seq.add_block(self._gx_pre if polarity == 'positive' else pp.scale_grad(self._gx_pre, -1)) - - # add readout gradients and ADCs - seq, time_to_echoes = self.add_to_seq_without_pre_post_gradient(seq, n_echoes, polarity) - - # readout re-winder - seq.add_block(self._gx_post if polarity == 'positive' else pp.scale_grad(self._gx_post, -1)) - - return seq, time_to_echoes - - def add_to_seq_without_pre_post_gradient( - self, seq: pp.Sequence, n_echoes: int, polarity: Literal['positive', 'negative'] = 'positive' - ): - """Add readout gradients without pre- and re-winder gradients. - - Often the pre- and re-winder gradients are played out at the same time as phase encoding gradients or spoiler - gradients. - - Parameters - ---------- - seq - PyPulseq Sequence object. - n_echoes - Number of echoes - polarity - Polarity of first readout gradient - - Returns - ------- - seq - PyPulseq Sequence object. - time_to_echoes - Time from beginning of sequence to echoes. - """ - # add readout gradient and ADC - time_to_echoes = [] - for echo_ in range(n_echoes): - start_of_current_gx = sum(seq.block_durations.values()) - gx_sign = (-1) ** echo_ if polarity == 'positive' else (-1) ** (echo_ + 1) - labels = [] - labels.append(pp.make_label(label='REV', type='SET', value=gx_sign == -1)) - labels.append(pp.make_label(label='ECO', type='SET', value=echo_)) - seq.add_block(pp.scale_grad(self._gx, gx_sign), self._adc, *labels) - time_to_echoes.append( - start_of_current_gx + self._adc.delay + self._n_readout_pre_echo * self._adc.dwell + self._adc.dwell / 2 - ) - start_of_current_gx = sum(seq.block_durations.values()) - if echo_ < n_echoes - 1: - if self._delta_te_delay > 0: - seq.add_block(pp.make_delay(self._delta_te_delay)) - seq.add_block(pp.scale_grad(self._gx_between, -gx_sign)) - - return seq, time_to_echoes - - -def undersampled_variable_density_spiral( - system: pp.Opts, n_readout: int, fov: float, undersampling_factor: float -) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, int, float, float]: - """Create undersampled variable density spiral. - - The distribution of the k-space points of a spiral trajectory are restricted by the maximum gradient amplitude and - slew rate. This makes an analytic solution for a given undersampling factor challenging. Here we use an iterative - approach in order to achieve a variable density spiral with a certain number of readout samplings and undersampling - factor. - - During the iterative search, the undersampling for the edge of k-space is increased. If this is not enough, then we - also start to increase the undersampling in the k-space center. The field-of-view varies linearly bewtween the - k-space center and k-space edge. - - If the undersampling factor is to high, it might not be possible to find a suitable solution. - - Parameters - ---------- - system - PyPulseq system limits object. - n_readout - Number of readout points per spiral. - fov - Field of view (in meters). - undersampling_factor - Undersampling factor of spiral trajectory - - Returns - ------- - tuple(np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, np.ndarray, int, float, float) - - k-space trajectory (traj) - - Gradient waveform (grad) - - Slew rate (slew) - - Time points for the trajectory (timing) - - Radius values (radius) - - Angular positions (theta) - - Number of spiral arms (n_spirals) - - Scaling of the field-of-view in the k-space center - - Scaling of the field-of-view in the k-space edge - - """ - # calculate single spiral trajectory - n_k0 = np.inf - fov_scaling_center = 1.0 - fov_scaling_edge = 1.0 - n_spirals = int(np.round(n_readout / undersampling_factor)) - while n_k0 > n_readout: - fov_coefficients = [fov * fov_scaling_center, -fov * (1 - fov_scaling_edge)] - - try: - traj, grad, slew, timing, radius, theta = variable_density_spiral_trajectory( - system=system, - sampling_period=system.grad_raster_time, - n_interleaves=n_spirals, - fov_coefficients=fov_coefficients, - max_kspace_radius=0.5 / fov * n_readout, - ) - n_k0 = len(grad) - fov_scaling_edge *= 0.95 - except ValueError: - # It is not possible to achieve the desired undersampling factor with the given system limits while keeping - # the full field-of-view in the k-space center. Reduce the field-of-view and try again. - n_k0 = np.inf - fov_scaling_center *= 0.95 - fov_scaling_edge = fov_scaling_center - - if fov_scaling_center < 0.1: - raise ValueError('Cannot find a suitable trajectory.') - - return traj, grad, slew, timing, radius, theta, n_spirals, fov_coefficients[0] / fov, fov_coefficients[1] / fov + 1 - - -def spiral_acquisition( - system: pp.Opts, - n_readout: int, - fov: float, - undersampling_factor: float, - readout_oversampling: Literal[1, 2, 4], - n_spirals: int | None, - max_pre_duration: float, - spiral_type: Literal['out', 'in-out'], -) -> tuple[list[Any], list[Any], Any, np.ndarray, float]: - """Generate a spiral acquisition sequence. - - Parameters - ---------- - system - PyPulseq system limits object. - n_readout - Number of readout points per spiral. - fov - Field of view (in meters). - undersampling_factor - Undersampling factor. - readout_oversampling - Oversampling factor for the readout trajectory. - n_spirals - Number of spirals to generate. If set to None, this value will be set based on the undersampling factor. - max_pre_duration : float - Maximum duration for pre-winder gradients (in seconds). - spiral_type - Type of spiral acquisition. 'out' for outward spirals, 'in-out' for spirals turning in and then out. - - Returns - ------- - gx_combined - List of combined gradient objects for the x-channel. - gy_combined - List of combined gradient objects for the y-channel. - adc - PyPulseq ADC object for the acquisition. - trajectory - K-space trajectory. - time_to_echo - Time to echo from beginning of gradients (in seconds). - """ - # calculate single spiral trajectory - traj, grad, _s, _timing, _r, _theta, n_spirals_undersampling, fov_scaling_center, fov_scaling_edge = ( - undersampled_variable_density_spiral(system, n_readout, fov, undersampling_factor) - ) - n_spirals = n_spirals_undersampling if n_spirals is None else n_spirals - print( - f'Target undersampling: {undersampling_factor} - ', - f'achieved undersampling: {n_readout**2 / (len(traj) * n_spirals_undersampling):.2f}', - f'FOV: {fov * fov_scaling_center:.3f} (k-sapce center) - {fov * fov_scaling_edge:.3f} (k-space edge)', - ) - - delta_angle = 2 * np.pi / n_spirals - n_samples_to_echo = 0.5 - if spiral_type == 'in-out': - n_samples_to_echo = len(grad) - grad = np.concatenate((-np.asarray(grad * np.exp(1j * np.pi))[::-1], grad)) - traj = np.concatenate((np.asarray(traj * np.exp(1j * np.pi))[::-1], traj)) - delta_angle = delta_angle / 2 - - # calculate ADC - n_readout_with_oversampling = len(grad) * readout_oversampling - adc_dwell_time = system.grad_raster_time / readout_oversampling - adc = pp.make_adc( - num_samples=n_readout_with_oversampling, dwell=adc_dwell_time, system=system, delay=system.adc_dead_time - ) - traj = np.interp( - np.linspace(0.5 / readout_oversampling, len(grad) - 0.5 / readout_oversampling, n_readout_with_oversampling), - np.linspace(0.5, len(grad) - 0.5, len(grad)), - traj, - ) - - print(f'Receiver bandwidth: {int(1.0 / (adc_dwell_time * n_readout_with_oversampling))} Hz/pixel') - - # Create gradient values and trajectory for different spirals - grad_list = [grad * np.exp(1j * delta_angle * idx) for idx in np.arange(n_spirals)] - traj_list = [traj * np.exp(1j * delta_angle * idx) for idx in np.arange(n_spirals)] - - # Create gradient objects - gx = [pp.make_arbitrary_grad(channel='x', waveform=g.real, delay=adc.delay, system=system) for g in grad_list] - gy = [pp.make_arbitrary_grad(channel='y', waveform=g.imag, delay=adc.delay, system=system) for g in grad_list] - - # Calculate pre- and re-winder gradients - gx_rew, gx_pre, gy_rew, gy_pre = [], [], [], [] - for gx_, gy_ in zip(gx, gy, strict=True): - gx_rew.append( - pp.make_extended_trapezoid_area( - area=-gx_.area if spiral_type == 'out' else -gx_.area / 2, - channel='x', - grad_start=gx_.last, - grad_end=0, - system=system, - convert_to_arbitrary=True, - )[0] - ) - gy_rew.append( - pp.make_extended_trapezoid_area( - area=-gy_.area if spiral_type == 'out' else -gy_.area / 2, - channel='y', - grad_start=gy_.last, - grad_end=0, - system=system, - convert_to_arbitrary=True, - )[0] - ) - - if spiral_type == 'in-out': - gx_pre.append( - pp.make_extended_trapezoid_area( - area=-gx_.area / 2, - channel='x', - grad_start=0, - grad_end=gx_.first, - system=system, - convert_to_arbitrary=True, - )[0] - ) - - gy_pre.append( - pp.make_extended_trapezoid_area( - area=-gy_.area / 2, - channel='y', - grad_start=0, - grad_end=gy_.first, - system=system, - convert_to_arbitrary=True, - )[0] - ) - else: - gx_pre.append(None) - gy_pre.append(None) - - if spiral_type == 'in-out': - adc.delay = max_pre_duration - - for i in range(len(gx_pre)): - gy_pre[i].delay = max_pre_duration - gy_pre[i].shape_dur - gx_pre[i].delay = max_pre_duration - gx_pre[i].shape_dur - else: - max_pre_duration = adc.delay - - def combine_gradients(*grad_objects, channel): - grad_list = [grad for grad in grad_objects if grad is not None] # Remove None - waveform_combined = np.concatenate([grad.waveform for grad in grad_list]) - - return pp.make_arbitrary_grad( - channel=channel, - waveform=waveform_combined, - first=0, - delay=grad_list[0].delay, - last=0, - system=system, - ) - - gx_combined = [ - combine_gradients(gx_pre, gx_in_out, gx_rew, channel='x') - for gx_pre, gx_in_out, gx_rew in zip(gx_pre, gx, gx_rew, strict=True) - ] - gy_combined = [ - combine_gradients(gy_pre, gy_in_out, gy_rew, channel='y') - for gy_pre, gy_in_out, gy_rew in zip(gy_pre, gy, gy_rew, strict=True) - ] - - # times -1 to match pulseq trajectory calculation - trajectory = -np.stack((np.asarray(traj_list).real, np.asarray(traj_list).imag), axis=-1) - - time_to_echo = max_pre_duration + n_samples_to_echo * readout_oversampling * adc.dwell - - return gx_combined, gy_combined, adc, trajectory, time_to_echo diff --git a/tests/scripts/__init__.py b/tests/sequences/__init__.py similarity index 100% rename from tests/scripts/__init__.py rename to tests/sequences/__init__.py diff --git a/tests/scripts/test_grpe_flash_dixon.py b/tests/sequences/test_grpe_flash_dixon.py similarity index 96% rename from tests/scripts/test_grpe_flash_dixon.py rename to tests/sequences/test_grpe_flash_dixon.py index 34759eb..0dce989 100644 --- a/tests/scripts/test_grpe_flash_dixon.py +++ b/tests/sequences/test_grpe_flash_dixon.py @@ -1,7 +1,7 @@ """Tests for 3D FLASH sequence with golden radial phase encoding.""" import pytest -from mrseq.scripts.grpe_flash_dixon import main as create_seq +from mrseq.sequences.grpe_flash_dixon import main as create_seq EXPECTED_DUR = 17.84831 # defined 2025-11-09 diff --git a/tests/scripts/test_radial_flash.py b/tests/sequences/test_radial_flash.py similarity index 96% rename from tests/scripts/test_radial_flash.py rename to tests/sequences/test_radial_flash.py index 8d031c3..2616d3c 100644 --- a/tests/scripts/test_radial_flash.py +++ b/tests/sequences/test_radial_flash.py @@ -2,7 +2,7 @@ import numpy as np import pytest -from mrseq.scripts.radial_flash import main as create_seq +from mrseq.sequences.radial_flash import main as create_seq EXPECTED_DUR = 0.75043 # defined 2025-10-17 diff --git a/tests/scripts/test_spiral_flash.py b/tests/sequences/test_spiral_flash.py similarity index 96% rename from tests/scripts/test_spiral_flash.py rename to tests/sequences/test_spiral_flash.py index 42983b2..5a55e69 100644 --- a/tests/scripts/test_spiral_flash.py +++ b/tests/sequences/test_spiral_flash.py @@ -2,7 +2,7 @@ import numpy as np import pytest -from mrseq.scripts.spiral_flash import main as create_seq +from mrseq.sequences.spiral_flash import main as create_seq EXPECTED_DUR = 0.70588 # defined 2026-02-10 diff --git a/tests/scripts/test_t1_inv_rec_gre_single_line.py b/tests/sequences/test_t1_inv_rec_gre_single_line.py similarity index 94% rename from tests/scripts/test_t1_inv_rec_gre_single_line.py rename to tests/sequences/test_t1_inv_rec_gre_single_line.py index 26685c0..a4e2d5c 100644 --- a/tests/scripts/test_t1_inv_rec_gre_single_line.py +++ b/tests/sequences/test_t1_inv_rec_gre_single_line.py @@ -1,7 +1,7 @@ """Tests for Gold standard GRE-based inversion recovery sequence.""" import pytest -from mrseq.scripts.t1_inv_rec_gre_single_line import main as create_seq +from mrseq.sequences.t1_inv_rec_gre_single_line import main as create_seq EXPECTED_DUR = 7168.000320 # defined 2025-02-03 diff --git a/tests/scripts/test_t1_inv_rec_se_single_line.py b/tests/sequences/test_t1_inv_rec_se_single_line.py similarity index 94% rename from tests/scripts/test_t1_inv_rec_se_single_line.py rename to tests/sequences/test_t1_inv_rec_se_single_line.py index 5b628c4..4f7a5cc 100644 --- a/tests/scripts/test_t1_inv_rec_se_single_line.py +++ b/tests/sequences/test_t1_inv_rec_se_single_line.py @@ -1,7 +1,7 @@ """Tests for Gold standard SE-based inversion recovery sequence.""" import pytest -from mrseq.scripts.t1_inv_rec_se_single_line import main as create_seq +from mrseq.sequences.t1_inv_rec_se_single_line import main as create_seq EXPECTED_DUR = 7168.000320 # defined 2025-02-03 diff --git a/tests/scripts/test_t1_molli_bssfp.py b/tests/sequences/test_t1_molli_bssfp.py similarity index 93% rename from tests/scripts/test_t1_molli_bssfp.py rename to tests/sequences/test_t1_molli_bssfp.py index 0c084de..9de7b25 100644 --- a/tests/scripts/test_t1_molli_bssfp.py +++ b/tests/sequences/test_t1_molli_bssfp.py @@ -1,7 +1,7 @@ """Tests for 2D Cartesian FLASH with T2-preparation pulses for T2 mapping.""" import pytest -from mrseq.scripts.t1_molli_bssfp import main as create_seq +from mrseq.sequences.t1_molli_bssfp import main as create_seq EXPECTED_DUR = 10.95339 # defined 2025-11-24 diff --git a/tests/scripts/test_t1_t2_spiral_cmrf.py b/tests/sequences/test_t1_t2_spiral_cmrf.py similarity index 90% rename from tests/scripts/test_t1_t2_spiral_cmrf.py rename to tests/sequences/test_t1_t2_spiral_cmrf.py index b210ed3..be7f6ea 100644 --- a/tests/scripts/test_t1_t2_spiral_cmrf.py +++ b/tests/sequences/test_t1_t2_spiral_cmrf.py @@ -1,7 +1,7 @@ """Tests for cardiac MR Fingerprinting sequence with spiral readout.""" import pytest -from mrseq.scripts.t1_t2_spiral_cmrf import main as create_seq +from mrseq.sequences.t1_t2_spiral_cmrf import main as create_seq EXPECTED_DUR = 14.55141 # defined 2026-02-10 diff --git a/tests/scripts/test_t1rho_se_single_line.py b/tests/sequences/test_t1rho_se_single_line.py similarity index 95% rename from tests/scripts/test_t1rho_se_single_line.py rename to tests/sequences/test_t1rho_se_single_line.py index 30b9662..161c37b 100644 --- a/tests/scripts/test_t1rho_se_single_line.py +++ b/tests/sequences/test_t1rho_se_single_line.py @@ -1,7 +1,7 @@ """Tests for Gold standard SE-based T1rho mapping sequence.""" import pytest -from mrseq.scripts.t1rho_se_single_line import main as create_seq +from mrseq.sequences.t1rho_se_single_line import main as create_seq EXPECTED_DUR = 3072.000000 # defined 2025-02-20 diff --git a/tests/scripts/test_t2_multi_echo_se_single_line.py b/tests/sequences/test_t2_multi_echo_se_single_line.py similarity index 94% rename from tests/scripts/test_t2_multi_echo_se_single_line.py rename to tests/sequences/test_t2_multi_echo_se_single_line.py index 2967f33..ed4312f 100644 --- a/tests/scripts/test_t2_multi_echo_se_single_line.py +++ b/tests/sequences/test_t2_multi_echo_se_single_line.py @@ -2,7 +2,7 @@ import numpy as np import pytest -from mrseq.scripts.t2_multi_echo_se_single_line import main as create_seq +from mrseq.sequences.t2_multi_echo_se_single_line import main as create_seq EXPECTED_DUR = 5120.000970 # defined 2025-02-06 diff --git a/tests/scripts/test_t2_t2prep_flash.py b/tests/sequences/test_t2_t2prep_flash.py similarity index 95% rename from tests/scripts/test_t2_t2prep_flash.py rename to tests/sequences/test_t2_t2prep_flash.py index d0f356c..dca3aa3 100644 --- a/tests/scripts/test_t2_t2prep_flash.py +++ b/tests/sequences/test_t2_t2prep_flash.py @@ -2,7 +2,7 @@ import numpy as np import pytest -from mrseq.scripts.t2_t2prep_flash import main as create_seq +from mrseq.sequences.t2_t2prep_flash import main as create_seq EXPECTED_DUR = 8.17667 # defined 2026-01-26 diff --git a/tests/scripts/test_t2star_multi_echo_flash.py b/tests/sequences/test_t2star_multi_echo_flash.py similarity index 94% rename from tests/scripts/test_t2star_multi_echo_flash.py rename to tests/sequences/test_t2star_multi_echo_flash.py index c428cd7..436a0d1 100644 --- a/tests/scripts/test_t2star_multi_echo_flash.py +++ b/tests/sequences/test_t2star_multi_echo_flash.py @@ -1,7 +1,7 @@ """Tests for 2D Cartesian FLASH with T2-preparation pulses for T2 mapping.""" import pytest -from mrseq.scripts.t2star_multi_echo_flash import main as create_seq +from mrseq.sequences.t2star_multi_echo_flash import main as create_seq EXPECTED_DUR = 2.83989 # defined 2025-11-22 diff --git a/tests/utils/test_cartesian_sampling.py b/tests/utils/test_cartesian_sampling.py new file mode 100644 index 0000000..64c3d3b --- /dev/null +++ b/tests/utils/test_cartesian_sampling.py @@ -0,0 +1,67 @@ +"""Tests for Cartesian sampling.""" + +from typing import Literal + +import numpy as np +import pytest +from mrseq.utils.cartesian_sampling import cartesian_phase_encoding + + +@pytest.mark.parametrize('n_phase_encoding', [50, 51, 100]) +@pytest.mark.parametrize('acceleration', [1, 2, 3, 4, 6]) +@pytest.mark.parametrize('n_fully_sampled_center', [0, 8, 9]) +def test_cartesian_phase_encoding_identical_points( + n_phase_encoding: int, acceleration: int, n_fully_sampled_center: int +): + """Test that linear, low-high and high-low cover same phase encoding points.""" + pe_linear, pe_center_linear = cartesian_phase_encoding( + n_phase_encoding, acceleration, n_fully_sampled_center, sampling_order='linear' + ) + pe_low_high, pe_center_low_high = cartesian_phase_encoding( + n_phase_encoding, acceleration, n_fully_sampled_center, sampling_order='low_high' + ) + pe_high_low, pe_center_high_low = cartesian_phase_encoding( + n_phase_encoding, acceleration, n_fully_sampled_center, sampling_order='high_low' + ) + + np.testing.assert_allclose(pe_linear, np.sort(pe_low_high)) + np.testing.assert_allclose(pe_linear, np.sort(pe_high_low)) + np.testing.assert_allclose(pe_center_linear, np.sort(pe_center_low_high)) + np.testing.assert_allclose(pe_center_linear, np.sort(pe_center_high_low)) + + +@pytest.mark.parametrize('pattern', ['linear', 'low_high', 'high_low', 'random']) +def test_cartesian_phase_encoding_acceleration(pattern: Literal['linear', 'low_high', 'high_low', 'random']): + """Test correct undersampling factor.""" + n_pe_full = 100 + acceleration = 4 + + pe, _ = cartesian_phase_encoding(n_phase_encoding=n_pe_full, acceleration=acceleration, sampling_order=pattern) + assert len(pe) == n_pe_full // acceleration + + +@pytest.mark.parametrize('pattern', ['linear', 'low_high', 'high_low', 'random']) +@pytest.mark.parametrize('n_phase_encoding_per_shot', [3, 8, 11, 13]) +def test_cartesian_phase_encoding_integer_shots( + pattern: Literal['linear', 'low_high', 'high_low', 'random'], + n_phase_encoding_per_shot: int, +): + """Test that the total number of phase encoding points lead to an integer number.""" + n_pe_full = 100 + acceleration = 4 + + pe, _ = cartesian_phase_encoding( + n_phase_encoding=n_pe_full, + acceleration=acceleration, + sampling_order=pattern, + n_phase_encoding_per_shot=n_phase_encoding_per_shot, + ) + assert np.mod(len(pe), n_phase_encoding_per_shot) == 0 + + +def test_cartesian_phase_encoding_warning_fully_sampled_center(): + """Test if warning is raised for a fully sampled center which is too large.""" + with pytest.raises(Warning, match='Number of phase encoding steps in the fully sampled center will be reduced'): + cartesian_phase_encoding( + n_phase_encoding=10, acceleration=1, sampling_order='linear', n_fully_sampled_center=12 + ) diff --git a/tests/utils/test_multi_echo_acquisition.py b/tests/utils/test_multi_echo_acquisition.py new file mode 100644 index 0000000..f15a0a2 --- /dev/null +++ b/tests/utils/test_multi_echo_acquisition.py @@ -0,0 +1,87 @@ +"""Tests for MultiEchoAcquisition.""" + +import numpy as np +import pypulseq as pp +import pytest +from mrseq.utils import MultiEchoAcquisition +from scipy.signal import argrelextrema + + +def test_multi_gradient_echo(system_defaults): + """Test multi-echo gradient echo readout as part of a simple sequence.""" + seq = pp.Sequence(system=system_defaults) + rf = pp.make_block_pulse( + flip_angle=np.pi, + delay=system_defaults.rf_dead_time, + duration=2e-3, + phase_offset=0.0, + system=system_defaults, + ) + seq.add_block(rf) + mecho = MultiEchoAcquisition(system=seq.system) + seq, _ = mecho.add_to_seq(seq, n_echoes=3) + + +@pytest.mark.parametrize('delta_te', [3e-3, 5.34e-3]) +def test_multi_gradient_echo_set_delta_te(delta_te, system_defaults): + """Test pre-defined delta te.""" + seq = pp.Sequence(system=system_defaults) + mecho = MultiEchoAcquisition(system=seq.system, delta_te=delta_te) + seq, time_to_echoes = mecho.add_to_seq(seq, 6) + np.testing.assert_allclose(np.diff(time_to_echoes), delta_te) + + +def test_multi_gradient_echo_error_on_short_delta_te(system_defaults): + """Test if error is raised on too short delta echo time.""" + with pytest.raises(ValueError): + MultiEchoAcquisition(system=system_defaults, delta_te=1e-6) + + +@pytest.mark.parametrize('n_echoes', [1, 3, 8]) +@pytest.mark.parametrize('readout_oversampling', [1, 1.5, 2]) +@pytest.mark.parametrize('n_readout', [64, 128, 200]) +@pytest.mark.parametrize('partial_echo_factor', [1.0, 0.8, 0.7]) +@pytest.mark.parametrize('polarity', ['positive', 'negative']) +def test_multi_gradient_echo_timing( + n_echoes, readout_oversampling, n_readout, partial_echo_factor, polarity, system_defaults +): + """Test that zero crossing of gradient moment coincides with echo time and correct adc sample.""" + seq = pp.Sequence(system=system_defaults) + mecho = MultiEchoAcquisition( + system=seq.system, + n_readout=n_readout, + readout_oversampling=readout_oversampling, + partial_echo_factor=partial_echo_factor, + ) + seq, time_to_echoes = mecho.add_to_seq(seq, n_echoes, polarity) + + # Get full waveform for readout gradient + w = seq.waveforms_and_times() + gx_waveform = w[0][0][1] + gx_waveform_time = w[0][0][0] + + # Find k0-crossings + max_grad = np.max(np.abs(gx_waveform)) + dt = np.arange(gx_waveform_time[0], gx_waveform_time[-1], step=mecho._adc.dwell / 10) + gx_waveform_intp = np.interp(dt, gx_waveform_time, gx_waveform / max_grad) + m0_intp = np.cumsum(gx_waveform_intp) / len(gx_waveform_intp) + k0_idx = argrelextrema(np.abs(m0_intp), np.less, order=100)[0] + + # Remove k0-crossings at the beginning and end of the block + k0_idx = np.asarray([ki for ki in k0_idx if (ki > 100 and ki < len(dt) - 100)]) + + # Zero-crossing of the readout gradient should be within +/- .5 adc dwell time of the k-space center sample which is + # the (_n_readout_pre_echo + 1)th sample. This should also be the same as the time_to_echo - time + assert n_echoes == len(k0_idx) + current_time = pp.calc_duration(mecho._gx_pre) + for echo in range(n_echoes): + time_of_k0_adc_sample = ( + current_time + mecho._adc.delay + mecho._n_readout_pre_echo * mecho._adc.dwell + mecho._adc.dwell / 2 + ) + print(dt[k0_idx[echo]], time_of_k0_adc_sample, mecho._adc.dwell) + assert np.isclose(dt[k0_idx[echo]], time_of_k0_adc_sample, atol=mecho._adc.dwell / 2) + assert np.isclose(dt[k0_idx[echo]], time_to_echoes[echo], atol=mecho._adc.dwell / 2) + + current_time += pp.calc_duration(mecho._gx) + pp.calc_duration(mecho._gx_between) + + assert seq is not None diff --git a/tests/utils/test_spiral_sampling.py b/tests/utils/test_spiral_sampling.py new file mode 100644 index 0000000..26b3463 --- /dev/null +++ b/tests/utils/test_spiral_sampling.py @@ -0,0 +1,138 @@ +"""Tests for spiral sampling.""" + +from typing import Literal + +import numpy as np +import pypulseq as pp +import pytest +from mrseq.utils import spiral_acquisition +from mrseq.utils.spiral_sampling import undersampled_variable_density_spiral + + +def get_interp_waveform_for_gx_gy(seq: pp.Sequence, dt: np.ndarray | None = None, scale: float = 1.0): + """Interpolate gradient waveforms for the x and y axes. + + Parameters + ---------- + seq + The PyPulseq sequence object containing gradient waveforms. + dt + Desired time points for interpolation. If None, a default time array is generated. + scale + Scaling factor for the gradient waveforms. Default is 1. + + Returns + ------- + gx_waveform_intp + Interpolated gradient waveform for the x-axis. + gy_waveform_intp + Interpolated gradient waveform for the y-axis. + dt + Time points corresponding to the interpolated waveforms. + """ + w = seq.waveforms_and_times() + gx_waveform = w[0][0][1] * scale + gx_waveform_time = w[0][0][0] + + gy_waveform = w[0][1][1] * scale + gy_waveform_time = w[0][1][0] + + if dt is None: + dt = np.arange( + min(gx_waveform_time[0], gy_waveform_time[0]), max(gx_waveform_time[-1], gy_waveform_time[-1]), step=1e-7 + ) + gx_waveform_intp = np.interp(dt, gx_waveform_time, gx_waveform) + gy_waveform_intp = np.interp(dt, gy_waveform_time, gy_waveform) + + return gx_waveform_intp, gy_waveform_intp, dt + + +@pytest.mark.parametrize('n_readout', (128, 256)) +@pytest.mark.parametrize('fov', (128e-3, 320e-3)) +@pytest.mark.parametrize('undersampling_factor', (1, 2, 3, 4)) +def test_undersampled_variable_density_spiral( + system_defaults: pp.Opts, n_readout: int, fov: float, undersampling_factor: float +): + """Test spiral for different undersampling factors.""" + traj, _grad, _s, _timing, _r, _theta, n_spirals, _fov_scaling_center, _fov_scaling_edge = ( + undersampled_variable_density_spiral(system_defaults, n_readout, fov, undersampling_factor) + ) + total_number_of_points = len(traj) * n_spirals + assert np.round(n_readout**2 / total_number_of_points) == undersampling_factor + + +@pytest.mark.parametrize('n_readout', (64, 256)) +@pytest.mark.parametrize('fov', (128e-3, 320e-3)) +@pytest.mark.parametrize('n_spirals', (14, None)) +@pytest.mark.parametrize('undersampling_factor', (1, 3)) +@pytest.mark.parametrize('readout_oversampling', (1, 2, 4)) +@pytest.mark.parametrize('spiral_type', ('out', 'in-out')) +def test_spiral_acquisition( + system_defaults: pp.Opts, + n_readout: int, + fov: float, + undersampling_factor: float, + n_spirals: int, + readout_oversampling: Literal[1, 2, 4], + spiral_type: Literal['out', 'in-out'], +): + """Test spiral trajectories for different parameter combinations.""" + g_pre_duration = 2e-3 # make this duration long to work for all combinations + system_defaults.adc_dead_time = 100e-6 # make this longer to ensure it impacts the timing + + gx, gy, adc, trajectory, time_to_echo = spiral_acquisition( + system_defaults, + n_readout, + fov, + undersampling_factor, + readout_oversampling, + n_spirals, + g_pre_duration, + spiral_type=spiral_type, + ) + + # Verify timing for each spiral arm + for spiral_idx in range(len(gx)): + seq = pp.Sequence(system=system_defaults) + seq.add_block(gx[spiral_idx], gy[spiral_idx], adc) + + # Get full waveform for readout gradient + gx_waveform_intp, gy_waveform_intp, dt = get_interp_waveform_for_gx_gy(seq) + max_grad = np.max(np.abs(gx_waveform_intp)) + gx_waveform_intp /= max_grad + gy_waveform_intp /= max_grad + + m0_intp = (np.abs(np.cumsum(gx_waveform_intp)) + np.abs(np.cumsum(gy_waveform_intp))) / ( + 2 * len(gx_waveform_intp) + ) + + if spiral_type == 'out': + k0_idx = np.argmin(m0_intp[:-100]) + else: + k0_idx = np.argmin(m0_intp[100:-100]) + 100 + assert m0_intp[0] < 1e-3 + assert m0_intp[-1] < 1e-3 + assert np.isclose(dt[k0_idx], time_to_echo, atol=system_defaults.grad_raster_time) + + k_traj_adc = seq.calculate_kspace()[0] + # Ignore first and last elements because they are extrapolated for readout_oversampling > 1 + k_traj_spiral = trajectory[spiral_idx, :, :] + if readout_oversampling > 1: + k_traj_adc = k_traj_adc[:, readout_oversampling // 2 : -readout_oversampling // 2] + k_traj_spiral = k_traj_spiral[readout_oversampling // 2 : -readout_oversampling // 2, :] + k_traj_adc /= np.max(np.abs(k_traj_adc)) + k_traj_spiral /= np.max(np.abs(k_traj_spiral)) + np.testing.assert_allclose(k_traj_adc[0, :], k_traj_spiral[:, 0], atol=2.5e-3) + np.testing.assert_allclose(k_traj_adc[1, :], k_traj_spiral[:, 1], atol=2.5e-3) + + # Verify entire trajectory + seq = pp.Sequence(system=system_defaults) + for idx in range(len(gx)): + seq.add_block(gx[idx], gy[idx], adc) + seq.add_block(pp.make_delay(1e-3)) + + ok, error_report = seq.check_timing() + if not ok: + print('\nTiming check failed! Error listing follows\n') + print(error_report) + assert ok diff --git a/tests/utils/test_trajectory.py b/tests/utils/test_trajectory.py deleted file mode 100644 index 352c852..0000000 --- a/tests/utils/test_trajectory.py +++ /dev/null @@ -1,281 +0,0 @@ -"""Tests for sequence helper functions.""" - -from typing import Literal - -import numpy as np -import pypulseq as pp -import pytest -from mrseq.utils import spiral_acquisition -from mrseq.utils.trajectory import MultiEchoAcquisition -from mrseq.utils.trajectory import cartesian_phase_encoding -from mrseq.utils.trajectory import undersampled_variable_density_spiral -from scipy.signal import argrelextrema - - -@pytest.mark.parametrize('n_phase_encoding', [50, 51, 100]) -@pytest.mark.parametrize('acceleration', [1, 2, 3, 4, 6]) -@pytest.mark.parametrize('n_fully_sampled_center', [0, 8, 9]) -def test_cartesian_phase_encoding_identical_points( - n_phase_encoding: int, acceleration: int, n_fully_sampled_center: int -): - """Test that linear, low-high and high-low cover same phase encoding points.""" - pe_linear, pe_center_linear = cartesian_phase_encoding( - n_phase_encoding, acceleration, n_fully_sampled_center, sampling_order='linear' - ) - pe_low_high, pe_center_low_high = cartesian_phase_encoding( - n_phase_encoding, acceleration, n_fully_sampled_center, sampling_order='low_high' - ) - pe_high_low, pe_center_high_low = cartesian_phase_encoding( - n_phase_encoding, acceleration, n_fully_sampled_center, sampling_order='high_low' - ) - - np.testing.assert_allclose(pe_linear, np.sort(pe_low_high)) - np.testing.assert_allclose(pe_linear, np.sort(pe_high_low)) - np.testing.assert_allclose(pe_center_linear, np.sort(pe_center_low_high)) - np.testing.assert_allclose(pe_center_linear, np.sort(pe_center_high_low)) - - -@pytest.mark.parametrize('pattern', ['linear', 'low_high', 'high_low', 'random']) -def test_cartesian_phase_encoding_acceleration(pattern: Literal['linear', 'low_high', 'high_low', 'random']): - """Test correct undersampling factor.""" - n_pe_full = 100 - acceleration = 4 - - pe, _ = cartesian_phase_encoding(n_phase_encoding=n_pe_full, acceleration=acceleration, sampling_order=pattern) - assert len(pe) == n_pe_full // acceleration - - -@pytest.mark.parametrize('pattern', ['linear', 'low_high', 'high_low', 'random']) -@pytest.mark.parametrize('n_phase_encoding_per_shot', [3, 8, 11, 13]) -def test_cartesian_phase_encoding_integer_shots( - pattern: Literal['linear', 'low_high', 'high_low', 'random'], - n_phase_encoding_per_shot: int, -): - """Test that the total number of phase encoding points lead to an integer number.""" - n_pe_full = 100 - acceleration = 4 - - pe, _ = cartesian_phase_encoding( - n_phase_encoding=n_pe_full, - acceleration=acceleration, - sampling_order=pattern, - n_phase_encoding_per_shot=n_phase_encoding_per_shot, - ) - assert np.mod(len(pe), n_phase_encoding_per_shot) == 0 - - -def test_cartesian_phase_encoding_warning_fully_sampled_center(): - """Test if warning is raised for a fully sampled center which is too large.""" - with pytest.raises(Warning, match='Number of phase encoding steps in the fully sampled center will be reduced'): - cartesian_phase_encoding( - n_phase_encoding=10, acceleration=1, sampling_order='linear', n_fully_sampled_center=12 - ) - - -def test_multi_gradient_echo(system_defaults): - """Test multi-echo gradient echo readout as part of a simple sequence.""" - seq = pp.Sequence(system=system_defaults) - rf = pp.make_block_pulse( - flip_angle=np.pi, - delay=system_defaults.rf_dead_time, - duration=2e-3, - phase_offset=0.0, - system=system_defaults, - ) - seq.add_block(rf) - mecho = MultiEchoAcquisition(system=seq.system) - seq, _ = mecho.add_to_seq(seq, n_echoes=3) - - -def get_interp_waveform_for_gx_gy(seq: pp.Sequence, dt: np.ndarray | None = None, scale: float = 1.0): - """Interpolate gradient waveforms for the x and y axes. - - Parameters - ---------- - seq - The PyPulseq sequence object containing gradient waveforms. - dt - Desired time points for interpolation. If None, a default time array is generated. - scale - Scaling factor for the gradient waveforms. Default is 1. - - Returns - ------- - gx_waveform_intp - Interpolated gradient waveform for the x-axis. - gy_waveform_intp - Interpolated gradient waveform for the y-axis. - dt - Time points corresponding to the interpolated waveforms. - """ - w = seq.waveforms_and_times() - gx_waveform = w[0][0][1] * scale - gx_waveform_time = w[0][0][0] - - gy_waveform = w[0][1][1] * scale - gy_waveform_time = w[0][1][0] - - if dt is None: - dt = np.arange( - min(gx_waveform_time[0], gy_waveform_time[0]), max(gx_waveform_time[-1], gy_waveform_time[-1]), step=1e-7 - ) - gx_waveform_intp = np.interp(dt, gx_waveform_time, gx_waveform) - gy_waveform_intp = np.interp(dt, gy_waveform_time, gy_waveform) - - return gx_waveform_intp, gy_waveform_intp, dt - - -@pytest.mark.parametrize('n_readout', (128, 256)) -@pytest.mark.parametrize('fov', (128e-3, 320e-3)) -@pytest.mark.parametrize('undersampling_factor', (1, 2, 3, 4)) -def test_undersampled_variable_density_spiral( - system_defaults: pp.Opts, n_readout: int, fov: float, undersampling_factor: float -): - """Test spiral for different undersampling factors.""" - traj, _grad, _s, _timing, _r, _theta, n_spirals, _fov_scaling_center, _fov_scaling_edge = ( - undersampled_variable_density_spiral(system_defaults, n_readout, fov, undersampling_factor) - ) - total_number_of_points = len(traj) * n_spirals - assert np.round(n_readout**2 / total_number_of_points) == undersampling_factor - - -@pytest.mark.parametrize('n_readout', (64, 256)) -@pytest.mark.parametrize('fov', (128e-3, 320e-3)) -@pytest.mark.parametrize('n_spirals', (14, None)) -@pytest.mark.parametrize('undersampling_factor', (1, 3)) -@pytest.mark.parametrize('readout_oversampling', (1, 2, 4)) -@pytest.mark.parametrize('spiral_type', ('out', 'in-out')) -def test_spiral_acquisition( - system_defaults: pp.Opts, - n_readout: int, - fov: float, - undersampling_factor: float, - n_spirals: int, - readout_oversampling: Literal[1, 2, 4], - spiral_type: Literal['out', 'in-out'], -): - """Test spiral trajectories for different parameter combinations.""" - g_pre_duration = 2e-3 # make this duration long to work for all combinations - system_defaults.adc_dead_time = 100e-6 # make this longer to ensure it impacts the timing - - gx, gy, adc, trajectory, time_to_echo = spiral_acquisition( - system_defaults, - n_readout, - fov, - undersampling_factor, - readout_oversampling, - n_spirals, - g_pre_duration, - spiral_type=spiral_type, - ) - - # Verify timing for each spiral arm - for spiral_idx in range(len(gx)): - seq = pp.Sequence(system=system_defaults) - seq.add_block(gx[spiral_idx], gy[spiral_idx], adc) - - # Get full waveform for readout gradient - gx_waveform_intp, gy_waveform_intp, dt = get_interp_waveform_for_gx_gy(seq) - max_grad = np.max(np.abs(gx_waveform_intp)) - gx_waveform_intp /= max_grad - gy_waveform_intp /= max_grad - - m0_intp = (np.abs(np.cumsum(gx_waveform_intp)) + np.abs(np.cumsum(gy_waveform_intp))) / ( - 2 * len(gx_waveform_intp) - ) - - if spiral_type == 'out': - k0_idx = np.argmin(m0_intp[:-100]) - else: - k0_idx = np.argmin(m0_intp[100:-100]) + 100 - assert m0_intp[0] < 1e-3 - assert m0_intp[-1] < 1e-3 - assert np.isclose(dt[k0_idx], time_to_echo, atol=system_defaults.grad_raster_time) - - k_traj_adc = seq.calculate_kspace()[0] - # Ignore first and last elements because they are extrapolated for readout_oversampling > 1 - k_traj_spiral = trajectory[spiral_idx, :, :] - if readout_oversampling > 1: - k_traj_adc = k_traj_adc[:, readout_oversampling // 2 : -readout_oversampling // 2] - k_traj_spiral = k_traj_spiral[readout_oversampling // 2 : -readout_oversampling // 2, :] - k_traj_adc /= np.max(np.abs(k_traj_adc)) - k_traj_spiral /= np.max(np.abs(k_traj_spiral)) - np.testing.assert_allclose(k_traj_adc[0, :], k_traj_spiral[:, 0], atol=2.5e-3) - np.testing.assert_allclose(k_traj_adc[1, :], k_traj_spiral[:, 1], atol=2.5e-3) - - # Verify entire trajectory - seq = pp.Sequence(system=system_defaults) - for idx in range(len(gx)): - seq.add_block(gx[idx], gy[idx], adc) - seq.add_block(pp.make_delay(1e-3)) - - ok, error_report = seq.check_timing() - if not ok: - print('\nTiming check failed! Error listing follows\n') - print(error_report) - assert ok - - -@pytest.mark.parametrize('delta_te', [3e-3, 5.34e-3]) -def test_multi_gradient_echo_set_delta_te(delta_te, system_defaults): - """Test pre-defined delta te.""" - seq = pp.Sequence(system=system_defaults) - mecho = MultiEchoAcquisition(system=seq.system, delta_te=delta_te) - seq, time_to_echoes = mecho.add_to_seq(seq, 6) - np.testing.assert_allclose(np.diff(time_to_echoes), delta_te) - - -def test_multi_gradient_echo_error_on_short_delta_te(system_defaults): - """Test if error is raised on too short delta echo time.""" - with pytest.raises(ValueError): - MultiEchoAcquisition(system=system_defaults, delta_te=1e-6) - - -@pytest.mark.parametrize('n_echoes', [1, 3, 8]) -@pytest.mark.parametrize('readout_oversampling', [1, 1.5, 2]) -@pytest.mark.parametrize('n_readout', [64, 128, 200]) -@pytest.mark.parametrize('partial_echo_factor', [1.0, 0.8, 0.7]) -@pytest.mark.parametrize('polarity', ['positive', 'negative']) -def test_multi_gradient_echo_timing( - n_echoes, readout_oversampling, n_readout, partial_echo_factor, polarity, system_defaults -): - """Test that zero crossing of gradient moment coincides with echo time and correct adc sample.""" - seq = pp.Sequence(system=system_defaults) - mecho = MultiEchoAcquisition( - system=seq.system, - n_readout=n_readout, - readout_oversampling=readout_oversampling, - partial_echo_factor=partial_echo_factor, - ) - seq, time_to_echoes = mecho.add_to_seq(seq, n_echoes, polarity) - - # Get full waveform for readout gradient - w = seq.waveforms_and_times() - gx_waveform = w[0][0][1] - gx_waveform_time = w[0][0][0] - - # Find k0-crossings - max_grad = np.max(np.abs(gx_waveform)) - dt = np.arange(gx_waveform_time[0], gx_waveform_time[-1], step=mecho._adc.dwell / 10) - gx_waveform_intp = np.interp(dt, gx_waveform_time, gx_waveform / max_grad) - m0_intp = np.cumsum(gx_waveform_intp) / len(gx_waveform_intp) - k0_idx = argrelextrema(np.abs(m0_intp), np.less, order=100)[0] - - # Remove k0-crossings at the beginning and end of the block - k0_idx = np.asarray([ki for ki in k0_idx if (ki > 100 and ki < len(dt) - 100)]) - - # Zero-crossing of the readout gradient should be within +/- .5 adc dwell time of the k-space center sample which is - # the (_n_readout_pre_echo + 1)th sample. This should also be the same as the time_to_echo - time - assert n_echoes == len(k0_idx) - current_time = pp.calc_duration(mecho._gx_pre) - for echo in range(n_echoes): - time_of_k0_adc_sample = ( - current_time + mecho._adc.delay + mecho._n_readout_pre_echo * mecho._adc.dwell + mecho._adc.dwell / 2 - ) - print(dt[k0_idx[echo]], time_of_k0_adc_sample, mecho._adc.dwell) - assert np.isclose(dt[k0_idx[echo]], time_of_k0_adc_sample, atol=mecho._adc.dwell / 2) - assert np.isclose(dt[k0_idx[echo]], time_to_echoes[echo], atol=mecho._adc.dwell / 2) - - current_time += pp.calc_duration(mecho._gx) + pp.calc_duration(mecho._gx_between) - - assert seq is not None