commit 6c719eaad277d3412445c824919fde49e6afa5e3 Author: jupfi Date: Tue Jun 25 15:59:36 2024 +0200 Different basic undersampling methods. diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..2462722 --- /dev/null +++ b/.gitignore @@ -0,0 +1,8 @@ +venv/ +*.cfl +*.hdr +fig/ +.ipynb_checkpoints/ +data/ +__pycache__/ +*.mat \ No newline at end of file diff --git a/ESMRMB_Educational.ipynb b/ESMRMB_Educational.ipynb new file mode 100644 index 0000000..34ee800 --- /dev/null +++ b/ESMRMB_Educational.ipynb @@ -0,0 +1,1355 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "ca29b946-be1d-4862-9e8d-7e89a0650348", + "metadata": {}, + "source": [ + "## 1 Setup BART\n", + "This notebook is designed to run on a local system. It uses the python kernel, however, almost all commands use the `%%bash` cell magic to be executed in a `bash` subshell.\n", + "\n", + "We will use BART version 0.9.00. In particular, we will take advantage of the newly added looping feature. For more information check the announcement on our [mailing list](https://lists.eecs.berkeley.edu/sympa/arc/mrirecon/2023-12/msg00000.html).\n", + "\n", + "\n", + "This version has been archived at CERN as: \n", + "\n", + "BART: version 0.9.00 (2023) DOI:10.5281/zenodo.10277939" + ] + }, + { + "cell_type": "markdown", + "id": "7602584e-b4dc-44cc-b5fc-67d92a9eac2b", + "metadata": {}, + "source": [ + "### 1.1 Local Usage\n", + "- Install bart from its [github repository](https://github.com/mrirecon/bart)\n", + "- Set the `BART_TOOLBOX_PATH` to the BART directory and add it to the `PATH`\n", + "\n", + "```bash\n", + "export BART_TOOLBOX_PATH=/path/to/bart \n", + "export PATH=$BART_TOOLBOX_PATH:$PATH\n", + "```\n", + "\n", + "Although the simplest way to call the BART CLI tools is from a terminal, there are also wrapper functions that allow the tools to be used from Matlab and Python. These are located under the `$BART_TOOLBOX_PATH/matlab` and `$BART_TOOLBOX_PATH/python` directories." + ] + }, + { + "cell_type": "markdown", + "id": "84076a2e-2983-44cb-81cc-0d50505ba1b2", + "metadata": {}, + "source": [ + "We will also use the [CFL viewer](https://github.com/mrirecon/view). Install it locally after configuring BART" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "d48523c5-e101-4589-94b6-4b7793572c97", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "os.environ['DEMO'] = '1'\n", + "os.environ['DEMO_INSTALL'] = '0'" + ] + }, + { + "cell_type": "markdown", + "id": "f62d94b7-07e9-4ae0-a3d7-3fa2b565895c", + "metadata": {}, + "source": [ + "### 1.2 Clone and compile BART v0.9.00\n", + "We clone BART into the current working directory of this notebook and delete any previous installation in this directory." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "72bc4d27-375b-405f-9089-c08e1c531096", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "skipped .. DEMO_INSTALL=0\n" + ] + } + ], + "source": [ + "%%bash\n", + "\n", + "# Clone Bart\n", + "if [ 1 -eq $DEMO_INSTALL ]; then\n", + "\n", + " [ -d bart ] && rm -r bart\n", + " git clone https://github.com/mrirecon/bart/ bart &> /dev/null\n", + "\n", + " cd bart\n", + "\n", + " make -j32 &> /dev/null\n", + "else\n", + " echo skipped .. DEMO_INSTALL=$DEMO_INSTALL\n", + "fi" + ] + }, + { + "cell_type": "markdown", + "id": "76468f36-d2d6-4fdb-ac53-b51ffce3ab48", + "metadata": {}, + "source": [ + "### 1.3 Add BART to PATH variable\n", + "\n", + "We add the BART directory to the PATH variable and include the python wrapper for reading *.cfl files:" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "a8854d78-bd89-4ea7-bc6f-eafa4a6f139e", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import sys\n", + "\n", + "os.environ['BART_TOOLBOX_PATH']=os.getcwd()+\"/bart/\"\n", + "os.environ['PATH'] = os.environ['BART_TOOLBOX_PATH'] + \":\" + os.environ['PATH']\n", + "sys.path.append(os.environ['BART_TOOLBOX_PATH'] + \"/python/\")\n", + "os.environ['DEBUG_LEVEL'] = '2'\n", + "os.environ['BART_DEBUG_LEVEL'] = '2'" + ] + }, + { + "cell_type": "markdown", + "id": "3d494c56-88ab-4440-8484-1968acbcd1cb", + "metadata": {}, + "source": [ + "Check BART setup:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "5abe217f-02ba-4ca4-9de2-7862af2d0ce6", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "# The BART used in this notebook:\n", + "/home/jpfitzer/bart-0.9.00/bart\n", + "# BART version: \n", + "v0.9.00\n" + ] + } + ], + "source": [ + "%%bash\n", + "\n", + "echo \"# The BART used in this notebook:\"\n", + "which bart\n", + "echo \"# BART version: \"\n", + "bart version" + ] + }, + { + "cell_type": "markdown", + "id": "e6803f18-fc73-4abe-941c-9898bb640ce7", + "metadata": {}, + "source": [ + "### 1.4 Install Interactive CFL Viewer\n" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "4476e853-4eb7-4937-84fc-82ed6e80d0a6", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "skipped .. DEMO_INSTALL=0\n" + ] + } + ], + "source": [ + "%%bash\n", + "\n", + "# Clone View\n", + "if [ 1 -eq $DEMO_INSTALL ]; then\n", + "\n", + " [ -d view ] && rm -r view\n", + " git clone https://github.com/mrirecon/view/ view &> /dev/null\n", + "\n", + " cd view\n", + "\n", + " make &> /dev/null\n", + "\n", + "else\n", + " echo skipped .. DEMO_INSTALL=$DEMO_INSTALL\n", + "fi" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "c1721b60-7cb3-4cc1-bee0-065352f8d55a", + "metadata": {}, + "outputs": [], + "source": [ + "os.environ['VIEW_TOOLBOX_PATH']=os.getcwd()+\"/view/\"\n", + "os.environ['PATH'] = os.environ['VIEW_TOOLBOX_PATH'] + \":\" + os.environ['PATH']" + ] + }, + { + "cell_type": "markdown", + "id": "79a7f963-a4fc-41d3-9998-5ba0a196a543", + "metadata": {}, + "source": [ + "Check view setup:" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "c60a8e73-8a0b-4cb8-9724-891d860307e4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "# The BART viewer used in this notebook:\n", + "/usr/bin/view\n" + ] + } + ], + "source": [ + "%%bash\n", + "\n", + "echo \"# The BART viewer used in this notebook:\"\n", + "which view" + ] + }, + { + "cell_type": "markdown", + "id": "e691812f-83fa-48eb-96b4-65ca85d21dbf", + "metadata": {}, + "source": [ + "### 1.5 Python Data Writer" + ] + }, + { + "cell_type": "code", + "execution_count": 151, + "id": "f7b7800e-8a4c-4852-b94c-2b47cd4a52de", + "metadata": {}, + "outputs": [], + "source": [ + "# Copyright 2013-2015. The Regents of the University of California.\n", + "# Copyright 2021. Uecker Lab. University Center Göttingen.\n", + "# All rights reserved. Use of this source code is governed by\n", + "# a BSD-style license which can be found in the LICENSE file.\n", + "#\n", + "# Authors:\n", + "# 2013 Martin Uecker \n", + "# 2015 Jonathan Tamir \n", + "\n", + "from __future__ import print_function\n", + "from __future__ import with_statement\n", + "\n", + "import numpy as np\n", + "import mmap\n", + "import os\n", + "\n", + "from IPython.display import Image\n", + "\n", + "\n", + "def readcfl(name):\n", + " # get dims from .hdr\n", + " with open(name + \".hdr\", \"rt\") as h:\n", + " h.readline() # skip\n", + " l = h.readline()\n", + " dims = [int(i) for i in l.split()]\n", + "\n", + " # remove singleton dimensions from the end\n", + " n = np.prod(dims)\n", + " dims_prod = np.cumprod(dims)\n", + " dims = dims[:np.searchsorted(dims_prod, n)+1]\n", + "\n", + " # load data and reshape into dims\n", + " with open(name + \".cfl\", \"rb\") as d:\n", + " a = np.fromfile(d, dtype=np.complex64, count=n);\n", + " return a.reshape(dims, order='F') # column-major\n", + "\n", + "def readmulticfl(name):\n", + " # get dims from .hdr\n", + " with open(name + \".hdr\", \"rt\") as h:\n", + " lines = h.read().splitlines()\n", + "\n", + " index_dim = 1 + lines.index('# Dimensions')\n", + " total_size = int(lines[index_dim])\n", + " index_sizes = 1 + lines.index('# SizesDimensions')\n", + " sizes = [int(i) for i in lines[index_sizes].split()]\n", + " index_dims = 1 + lines.index('# MultiDimensions')\n", + "\n", + " with open(name + \".cfl\", \"rb\") as d:\n", + " a = np.fromfile(d, dtype=np.complex64, count=total_size)\n", + "\n", + " offset = 0\n", + " result = []\n", + " for i in range(len(sizes)):\n", + " dims = ([int(i) for i in lines[index_dims + i].split()])\n", + " n = np.prod(dims)\n", + " result.append(a[offset:offset+n].reshape(dims, order='F'))\n", + " offset += n\n", + "\n", + " if total_size != offset:\n", + " print(\"Error\")\n", + "\n", + " return result\n", + "\n", + "\n", + "def writecfl(name, array):\n", + " with open(name + \".hdr\", \"wt\") as h:\n", + " h.write('# Dimensions\\n')\n", + " for i in (array.shape):\n", + " h.write(\"%d \" % i)\n", + " h.write('\\n')\n", + "\n", + " size = np.prod(array.shape) * np.dtype(np.complex64).itemsize\n", + "\n", + " with open(name + \".cfl\", \"a+b\") as d:\n", + " os.ftruncate(d.fileno(), size)\n", + " mm = mmap.mmap(d.fileno(), size, flags=mmap.MAP_SHARED, prot=mmap.PROT_WRITE)\n", + " if array.dtype != np.complex64:\n", + " array = array.astype(np.complex64)\n", + " mm.write(np.ascontiguousarray(array.T))\n", + " mm.close()\n", + " #with mmap.mmap(d.fileno(), size, flags=mmap.MAP_SHARED, prot=mmap.PROT_WRITE) as mm:\n", + " # mm.write(array.astype(np.complex64).tobytes(order='F'))\n", + "\n", + "def writemulticfl(name, arrays):\n", + " size = 0\n", + " dims = []\n", + "\n", + " for array in arrays:\n", + " size += array.size\n", + " dims.append(array.shape)\n", + "\n", + " with open(name + \".hdr\", \"wt\") as h:\n", + " h.write('# Dimensions\\n')\n", + " h.write(\"%d\\n\" % size)\n", + "\n", + " h.write('# SizesDimensions\\n')\n", + " for dim in dims:\n", + " h.write(\"%d \" % len(dim))\n", + " h.write('\\n')\n", + "\n", + " h.write('# MultiDimensions\\n')\n", + " for dim in dims:\n", + " for i in dim:\n", + " h.write(\"%d \" % i)\n", + " h.write('\\n')\n", + " \n", + " size = size * np.dtype(np.complex64).itemsize\n", + "\n", + " with open(name + \".cfl\", \"a+b\") as d:\n", + " os.ftruncate(d.fileno(), size)\n", + " mm = mmap.mmap(d.fileno(), size, flags=mmap.MAP_SHARED, prot=mmap.PROT_WRITE)\n", + " for array in arrays:\n", + " if array.dtype != np.complex64:\n", + " array = array.astype(np.complex64)\n", + " mm.write(np.ascontiguousarray(array.T))\n", + " mm.close()\n" + ] + }, + { + "cell_type": "markdown", + "id": "77d0c22c-bb16-4a5d-916b-0d4263fc8c50", + "metadata": {}, + "source": [ + "## 2. Reading in the data" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "f1883fa9-f260-47a3-9598-474905d10f89", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Echo train length: 20\n", + "Echo spacing: 10 ms\n", + "(16, 80, 80, 20)\n" + ] + } + ], + "source": [ + "import scipy.io as sio\n", + "\n", + "rawName = 'T2_CPMG.mat'\n", + "mat_data_0=sio.loadmat(rawName)\n", + "\n", + "# Format: sampled, # kSpace [kRd, kPh, kSl, kSpace_echo_1, kSpace_echo_2, ..., kSpace_echo_nETL] (102400, 23)\n", + "# So the first three are the coordinates of the kspace, and the rest are the echoes\n", + "\n", + "kSpaces3D = mat_data_0['kSpaces3D']\n", + "# self.mapVals['sampled'] = np.concatenate((kRD, kPH, kSL, dataAll_sampled), axis=1)\n", + "\n", + "# nReadout, nPhase, nSlice\n", + "nPoints = (80, 80, 16)\n", + "\n", + "echo_train_length = mat_data_0['kSpaces3D'].shape[1] - 3 # Because the first 3 are kRD, kPH, kSL -> should give 20\n", + "print(f\"Echo train length: {echo_train_length}\")\n", + "\n", + "echo_spacing = mat_data_0['echoSpacing'][0][0]\n", + "print(f\"Echo spacing: {echo_spacing} ms\")\n", + "\n", + "k_readout = kSpaces3D[:, 0]\n", + "k_phase = kSpaces3D[:, 1]\n", + "k_slice = kSpaces3D[:, 2]\n", + "\n", + "# The rest of the data is the echoes\n", + "echos = kSpaces3D[:, 3:]\n", + "\n", + "# Reshape the kspace data for bart\n", + "kSpace = echos.reshape(nPoints[2], nPoints[1], nPoints[0], echo_train_length)\n", + "\n", + "print(kSpace.shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 208, + "id": "96b28fec-ff5a-44a3-8a0c-9b01b60cff0c", + "metadata": {}, + "outputs": [], + "source": [ + "# Write the cfl file\n", + "cfl = writecfl('data/kSpace', kSpace)" + ] + }, + { + "cell_type": "markdown", + "id": "150faf49-acd6-4ee3-ad36-3a745e8e88d9", + "metadata": {}, + "source": [ + "### 2.1 Data Exploration" + ] + }, + { + "cell_type": "code", + "execution_count": 209, + "id": "5f93fcca-48ab-4d42-9061-91361aadfb35", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "ksp:\n", + "Type: complex float\n", + "Dimensions: 16\n", + "AoD:\t16\t80\t80\t20\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\n" + ] + } + ], + "source": [ + "%%bash\n", + "\n", + "echo \"ksp:\"\n", + "bart show -m data/kSpace\n", + "\n", + "DIM_X=80\n", + "DIM_Y=80\n", + "DIM_Z=16" + ] + }, + { + "cell_type": "markdown", + "id": "188249e4-dbbb-4f7b-9d28-e28338de14c3", + "metadata": {}, + "source": [ + "### 2.2 FFT\n", + "\n", + "We have three dimensional data so we perform the fft along the first three dimension. " + ] + }, + { + "cell_type": "code", + "execution_count": 210, + "id": "c074e817-5989-4e83-a108-79834c7b0bdf", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "\n", + "bart fft -i $(bart bitmask 0 1 2) data/kSpace data/fft" + ] + }, + { + "cell_type": "code", + "execution_count": 211, + "id": "fcfc949d-8046-4868-bdfe-810b5f27454c", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "\n", + "# view data/fft" + ] + }, + { + "cell_type": "code", + "execution_count": 212, + "id": "4c72186c-61e3-4e1d-9fe7-e95088ab8c38", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "\n", + "# Put the echos on the fifth dimension:\n", + "bart transpose 3 5 data/fft data/fft_transposed\n", + "# Put the slices on the correct dimension\n", + "bart transpose 0 2 data/fft_transposed data/fft_transposed" + ] + }, + { + "cell_type": "code", + "execution_count": 213, + "id": "0a9648a7-6c94-489e-b637-c4338b8cf113", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "fft_transposed:\n", + "Type: complex float\n", + "Dimensions: 16\n", + "AoD:\t80\t80\t16\t1\t1\t20\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\n" + ] + } + ], + "source": [ + "%%bash\n", + "\n", + "echo \"fft_transposed:\"\n", + "bart show -m data/fft_transposed" + ] + }, + { + "cell_type": "code", + "execution_count": 214, + "id": "b9dade46-455d-4816-9b35-1bf6abf53dcc", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "# Create the echotimes file:\n", + "bart vec 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 data/echo_times\n", + "bart scale 0.001 data/echo_times data/echo_times_scaled\n", + "# Move the echo_times to the correct dimension:\n", + "bart transpose 0 5 data/echo_times_scaled data/echo_times_final" + ] + }, + { + "cell_type": "code", + "execution_count": 215, + "id": "a119e34e-8f2f-42e6-81f9-4023f753fe81", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Usage: mobafit [-T] [-I] [-L] [-G] [-D] [-m d] [-a] [-i d] [-g] [-B ] [--init [f:]*f] [--scale [f:]*f] [--min-flag d] [--max-flag d] [--max-mag-flag d] [--min [f:]*f] [--max [f:]*f] [] \n", + "\n", + "Pixel-wise fitting of physical signal models.\n", + "\n", + "-T TSE\n", + "-I Inversion Recovery: f(M0, R1, c) = M0 * (1 - exp(-t * R1 + c))\n", + "-L Inversion Recovery Look-Locker\n", + "-G MGRE\n", + "-D diffusion\n", + "-m model Select the MGRE model from enum { WF = 0, WFR2S, WF2R2S, R2S, PHASEDIFF } [default: WFR2S]\n", + "-a fit magnitude of signal model to data\n", + "-i iter Number of IRGNM steps\n", + "-g use gpu\n", + "-B file temporal (or other) basis\n", + "--init [f:]*f Initial values of parameters in model-based reconstruction\n", + "--scale [f:]*f Scaling\n", + "--min-flag flags Apply minimum constraint on selected maps\n", + "--max-flag flags Apply maximum constraint on selected maps\n", + "--max-mag-flag flags Apply maximum magnitude constraint on selected maps\n", + "--min [f:]*f Min bound (map must be selected with \"min-flag\")\n", + "--max [f:]*f Max bound (map must be selected with \"max-flag\" or \"max-mag-flag\")\n", + "-h help\n" + ] + } + ], + "source": [ + "%%bash\n", + "bart mobafit -h" + ] + }, + { + "cell_type": "code", + "execution_count": 216, + "id": "66f1afe4-7598-47c0-ae99-581beca68a11", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "# Fit the model:\n", + "bart mobafit -T data/echo_times_final data/fft_transposed data/fit" + ] + }, + { + "cell_type": "code", + "execution_count": 217, + "id": "c1f0ece5-2118-4bbc-80c2-e692a01b5f6b", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "# Now some values will be very large so we can apply a threshold to obtain a mask\n", + "MAX_T2=1000\n", + "\n", + "bart threshold -M $MAX_T2 data/fit data/mask" + ] + }, + { + "cell_type": "code", + "execution_count": 218, + "id": "68fe0ec1-b8a8-453f-8eef-9c9645b7b064", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "# Multiply the fit with the mask\n", + "bart fmac data/fit data/mask data/fit_mask" + ] + }, + { + "cell_type": "code", + "execution_count": 219, + "id": "54dcbf68-b920-48d6-a072-fe25f52e712a", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "# Select slice\n", + "bart slice 6 1 data/fit_mask data/R2_map" + ] + }, + { + "cell_type": "code", + "execution_count": 220, + "id": "5e159890-56b0-4629-89f8-abc33eba42c5", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "# Invert the data to get T2\n", + "bart invert data/R2_map data/T2_map " + ] + }, + { + "cell_type": "markdown", + "id": "f8d7d303-2747-4f41-bf5d-153bc83f2382", + "metadata": {}, + "source": [ + "## 3.1 Undersampling Poisson ky - kz" + ] + }, + { + "cell_type": "code", + "execution_count": 221, + "id": "65569e35-233b-4ac9-baf4-6668fcc581a4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "points: 448, grid size: 80x16x(pi/4) = 1005 (R = 2.243995)\n", + "Type: complex float\n", + "Dimensions: 16\n", + "AoD:\t1\t80\t16\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\n", + "Type: complex float\n", + "Dimensions: 16\n", + "AoD:\t80\t80\t16\t20\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\n" + ] + } + ], + "source": [ + "%%bash\n", + "ACC_Y=2\n", + "ACC_Z=2\n", + "\n", + "bart poisson -Y 80 -Z 16 -y $ACC_Y -z $ACC_Z -C 16 -e data/vd_mask\n", + "\n", + "bart show -m data/vd_mask\n", + "\n", + "bart transpose 0 2 data/kSpace data/kSpace_transposed\n", + "\n", + "bart show -m data/kSpace_transposed" + ] + }, + { + "cell_type": "code", + "execution_count": 222, + "id": "822e7962-12d6-47a2-a33e-f3d31e34100b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Writing 1 image(s)...done.\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFAAAAAQCAIAAAAOK2+WAAABY0lEQVRIia1WS5YEIQgT73/nmoXTvHQIqFSzUrQIn4Bl46fyPI+Z0eLWyPqwNtuW+ebj8YnHo3Jv2m6ZGeZord8n0eUiYIJBV1Dj13puUVRoX4YtvSokrcOyKyl0pazhhUMBtPCEoPFOdj+tMGJI3hZKhzwOk6HRwtqSEnmEfVT486+PeDJPRfK228t4uche4RE6KNsSKfBoYsLcqG9pCNEIISfkaSNaCWdmVDSJi24Tpde1GadC3EaeLKNZqOhfTyj1A0oSlYUbsSnmSZUWBoadJXt8E6QdLXGYipaNa3rShirynvR4Yevi+0Edgyn6ljTFqadvIoz8bajfW6TArwZ1NimHoo8cLlHjrNTPaaxYXcOiGr0KZxDboy0jNnOlYKnMiISsIUjqnqI6y9ji0Zf9EyfqBpY345g9lHNmNdwbh//S1EjURbjAKdB+luQjPHKyXAH1n8pC5AA7FPlYXNWwlj8kWyUQianQ+gAAAABJRU5ErkJggg==", + "text/plain": [ + "" + ] + }, + "execution_count": 222, + "metadata": { + "image/png": { + "width": 1000 + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "!bart transpose 1 2 data/vd_mask data/img\n", + "!bart toimg -w 1.0 data/img fig/vd_mask.png\n", + "Image(\"fig/vd_mask.png\", width=1000)" + ] + }, + { + "cell_type": "code", + "execution_count": 223, + "id": "4693848f-fe88-42e5-9200-17d38dfe2765", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "bart fmac data/vd_mask data/kSpace_transposed data/kSpace_vd" + ] + }, + { + "cell_type": "code", + "execution_count": 224, + "id": "10817602-2bd3-40a1-950f-2e9adc7ce97f", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Scaling: 0.349166\n" + ] + } + ], + "source": [ + "%%bash\n", + "# For nlinv the echos can't be in the coil dimension (3)\n", + "bart transpose 3 5 data/kSpace_vd data/kSpace_vd_transposed\n", + "# The scaling seems to help\n", + "bart nlinv -i 30 -x 80:80:16 -S data/kSpace_vd_transposed data/nlinv_vd" + ] + }, + { + "cell_type": "code", + "execution_count": 225, + "id": "2448bf68-a7c0-43ee-a6b1-1c74e058308e", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash \n", + "bart pattern data/kSpace_vd data/pattern_vd" + ] + }, + { + "cell_type": "code", + "execution_count": 226, + "id": "41288837-b086-4ebb-9a71-0dbb7f722ef1", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "# Fit the model:\n", + "bart mobafit -T data/echo_times_final data/nlinv_vd data/fit_vd\n", + "\n", + "# Now some values will be very large so we can apply a threshold to obtain a mask\n", + "MAX_T2=1000\n", + "\n", + "bart threshold -M $MAX_T2 data/fit_vd data/mask_vd\n", + "\n", + "# Multiply the fit with the mask\n", + "bart fmac data/fit_vd data/mask_vd data/fit_mask_vd\n", + "\n", + "# Select slice\n", + "bart slice 6 1 data/fit_mask_vd data/R2_map_vd\n", + "\n", + "# Invert the data to get T2\n", + "bart invert data/R2_map_vd data/T2_map_vd " + ] + }, + { + "cell_type": "markdown", + "id": "6a1e8727-0677-4c9c-8e0e-e3622cef5c17", + "metadata": {}, + "source": [ + "### 3.2 Undersampling Poisson-Disc in-plane\n", + "This isn't useful for a real sequence because then you undersample in readout direction." + ] + }, + { + "cell_type": "code", + "execution_count": 227, + "id": "67334262-2ce1-4c36-8198-0235c7423a84", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "points: 1450, grid size: 80x80x(pi/4) = 5026 (R = 3.466585)\n", + "Type: complex float\n", + "Dimensions: 16\n", + "AoD:\t80\t80\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\n", + "Type: complex float\n", + "Dimensions: 16\n", + "AoD:\t80\t80\t16\t20\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\t1\n", + "Scaling: 0.351537\n" + ] + } + ], + "source": [ + "%%bash\n", + "\n", + "ACC_Y=2\n", + "ACC_Z=2\n", + "\n", + "bart poisson -Y 80 -Z 80 -y $ACC_Y -z $ACC_Z -C 16 -e data/ip_mask\n", + "\n", + "bart squeeze data/ip_mask data/ip_mask\n", + "\n", + "bart show -m data/ip_mask\n", + "\n", + "bart show -m data/kSpace_transposed\n", + "\n", + "bart fmac data/ip_mask data/kSpace_transposed data/kSpace_ip\n", + "\n", + "# For nlinv the echos can't be in the coil dimension (3)\n", + "bart transpose 3 5 data/kSpace_ip data/kSpace_ip_transposed\n", + "\n", + "bart nlinv -i 30 -x 80:80:16 -S data/kSpace_ip_transposed data/nlinv_ip\n", + "# bart fft -i $(bart bitmask 0 1 2) data/kSpace_ip data/fft_ip\n", + "\n", + "bart pattern data/kSpace_ip data/pattern_ip\n", + "\n", + "# bart transpose 3 5 data/fft_ip data/fft_ip_transposed\n", + "\n", + "# Fit the model:\n", + "bart mobafit -T data/echo_times_final data/nlinv_ip data/fit_ip\n", + "\n", + "# Now some values will be very large so we can apply a threshold to obtain a mask\n", + "MAX_T2=1000\n", + "\n", + "bart threshold -M $MAX_T2 data/fit_ip data/mask_ip\n", + "\n", + "# Multiply the fit with the mask\n", + "bart fmac data/fit_ip data/mask_ip data/fit_mask_ip\n", + "\n", + "# Select slice\n", + "bart slice 6 1 data/fit_mask_ip data/R2_map_ip\n", + "\n", + "# Invert the data to get T2\n", + "bart invert data/R2_map_ip data/T2_map_ip " + ] + }, + { + "cell_type": "code", + "execution_count": 228, + "id": "d53121cd-a255-492b-a037-e6b9c6fbf6f2", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Writing 1 image(s)...done.\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFAAAABQCAIAAAABc2X6AAAEzklEQVR4nN1XSZLjMAxT/P8/ew5dpVFIAAQV25MeHrpsmQu4QekxHpfzPMPf9eF3S0jjP8kKyppbTnu2l336dmGTKdLOatnJtw88TA+2MetDJ47+0wKb4yDTbewW5SFxJnbdWF8tR/FLeYvkVNlfaGuqwZNHMy/DsEKUtjrtf3y95fEzk3F+cuQMhf6NmZfkyeaN8VmZUlkI6P9Kabk2N21jzkVNrxRR2vFeCwel9hb8+MN/gZTl7E5ULk15h5stvSBzwY1+G8vdhn50V8PzjdsbTvwRFZ+66wDjtm4NhclcJBZjb82E1XM3sA/CWUVn7HXTyuutUZqSReCnTEIQqEZfapoObxQTq0NdjKscAiud12kEiOYoOuO0d8202N5NdXMHPIfQp4brsK7AbKXA5tBEU0IUOZdw2QXWHTEsJUmK0uRXk8k1+m4Dd4ja+WryhMn2pS3jtnBSrwPjA5NgtTcx9kzH3B2ITQDDUsaAmoLDy/B7EfN5j71YyWEDHVpirtgrk+58tZxj6Oxwe5zM7vmDrfHEk3Kv2IJpuFbsykSPTDeQG1VDKYdN8J/Zk3L/hdsYomQ8p+Qa7iBz4SyIT04lqjHGK+i9Xi/9Og/Xr+unbMViY0DvPhmYAHsaZlTh8Fgxzc8iwE9KQXnVD05a2U5XATcMB2GvmOfzqnCsOQS9NR6znyfBak1+Q3IOIZ8ZcZ6HAo1U9B+FIzgN7Z1Tmn3NlEq22Mt2ykQIeSgPVGhMMHzLSjckqzkmnycMYWzDO0YaP0boeUXDLjCTjWzhTK0K+ZMJr3eVQ4U8veG8KwJAGQvavtFN9phXKDOzuCEY55kCLzZ4HUIYTHl+PQIsyMA/ZpClQ2nxFDUl346i6CFEWM/s6ghlgHBL12ue7FZsSWgU3M8ZDm5l6Nl8PqCjnAasKOQnSBitVEdqFEsm9zxMHECYJzMwx0i9ZR4Zpj3GYrGcQ+bkPM/6Fg3107dcVtjo84v8LDfjBuRBDj0D+VfLkJX+fIHXbYLMAinzlX4dh9+Of1GJ2Oc7Q7ZKnoH6ku8kFgiiElDH+t9SXqHMkGvsqcxuKVEInS0jkRwoc9tIjAu3PboWTKBfmUIR1RPhOY+uClRqsJLnV41JBGV1KZkZQgrPNLWySGX/tfkgJWjVQvhpiCheWU5HzdF3cnCGUfxVCPwZYycbAyJedQnaoncG6pTsNc+3V0NDgng2R1230QTU9TPVGPFqVyzVIqJmQvHJ5G2mUPLcfNaEF6zcbAUm7bFkqZKBmZqmPc049WznAKKugti2d1Vn1ZjYUkwqcjyUUyrK1yIO5tAFr8dpA5zQ1HuxYf55w4BfZ291LT5nMm2io1Ava28dZi7D6NZp8xaD9kZaiEnCI7VIEBvE57faH59CnN6WdS3b1fokgF3f0tKps7R6DgdKxlnma4jKd90iFWbVJbxwd5hgGqKpP0M0aVbPv0NCV/ZTe9+LJDZwPjsFur29TiT/DiunWjjXJldKq70l6+jqmDA2zNui2YWhEQMM1cx63Z6tFuc2cibzQ6q7S8rFZsrra+4kO4dOrszHF72Z2yunb6l/k63eTKHPuF3vsPD8tJjkVJqb598iedIgaTkr8O2pZgnzuf2z4fdlvorJzOvX3yqCch+WP60uH9jRdpH8AAAAAElFTkSuQmCC", + "text/plain": [ + "" + ] + }, + "execution_count": 228, + "metadata": { + "image/png": { + "width": 500 + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "!bart transpose 1 2 data/ip_mask data/img\n", + "!bart toimg -w 1.0 data/img fig/ip_mask.png\n", + "Image(\"fig/ip_mask.png\", width=500)" + ] + }, + { + "cell_type": "markdown", + "id": "8da50e61-6996-44c9-85d0-e97cef04f4c3", + "metadata": {}, + "source": [ + "## 3.3 Uniform Undersampling\n", + "Using bart upat - same pattern for every echo." + ] + }, + { + "cell_type": "code", + "execution_count": 229, + "id": "a9aa3845-2dc3-4e83-b68f-49b965af11ab", + "metadata": {}, + "outputs": [], + "source": [ + "%%bash\n", + "ACC_Y=2\n", + "ACC_Z=2\n", + "CENTER_SIZE=8\n", + "\n", + "bart upat -Y 80 -Z 16 -y $ACC_Y -z $ACC_Z -c $CENTER_SIZE data/upat" + ] + }, + { + "cell_type": "code", + "execution_count": 230, + "id": "ad44e09a-a74e-493e-bb3f-3cd3f7607eb6", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Scaling: 0.349421\n" + ] + } + ], + "source": [ + "%%bash\n", + "bart fmac data/upat data/kSpace_transposed data/kSpace_upat\n", + "\n", + "# For nlinv the echos can't be in the coil dimension (3)\n", + "bart transpose 3 5 data/kSpace_upat data/kSpace_upat_transposed\n", + "\n", + "bart nlinv -i 30 -x 80:80:16 -S data/kSpace_upat_transposed data/nlinv_upat\n", + "# bart fft -i $(bart bitmask 0 1 2) data/kSpace_upat data/fft_upat\n", + "\n", + "# bart transpose 3 5 data/fft_upat data/fft_upat_transposed\n", + "\n", + "# Fit the model:\n", + "bart mobafit -T data/echo_times_final data/nlinv_upat data/fit_upat\n", + "\n", + "# Now some values will be very large so we can apply a threshold to obtain a mask\n", + "MAX_T2=1000\n", + "\n", + "bart threshold -M $MAX_T2 data/fit_upat data/mask_upat\n", + "\n", + "# Multiply the fit with the mask\n", + "bart fmac data/fit_upat data/mask_upat data/fit_mask_upat\n", + "\n", + "# Select slice\n", + "bart slice 6 1 data/fit_mask_upat data/R2_map_upat\n", + "\n", + "# Invert the data to get T2\n", + "bart invert data/R2_map_upat data/T2_map_upat " + ] + }, + { + "cell_type": "code", + "execution_count": 231, + "id": "5132b2af-0e1d-44b9-90ba-d6031ca4d10a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Writing 1 image(s)...done.\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAFAAAAAQCAIAAAAOK2+WAAAARklEQVRIie3RIRIAMAgDwdD//5mKeqYiqemtwMJw6m5J/8xS3tl0qSp70nrz13sUNqNwAIUnFDajcACFJxQ2o3AAhScUNtsCZEIXIDnDpQAAAABJRU5ErkJggg==", + "text/plain": [ + "" + ] + }, + "execution_count": 231, + "metadata": { + "image/png": { + "width": 1000 + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "!bart transpose 1 2 data/upat data/img\n", + "!bart toimg -w 1.0 data/img fig/upat_mask.png\n", + "Image(\"fig/upat_mask.png\", width=1000)" + ] + }, + { + "cell_type": "markdown", + "id": "4c359bc3-c6e2-425b-acb7-52b5894ad654", + "metadata": {}, + "source": [ + "### 3.3.1 Uniform Undersampling\n", + "With varying patterns along the echos" + ] + }, + { + "cell_type": "code", + "execution_count": 246, + "id": "0882a788-ecd2-4968-9715-a57816d66e9c", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Pattern creation completed. Output file: upat_combined\n" + ] + } + ], + "source": [ + "%%bash\n", + "# Parameters\n", + "R=2 # Undersampling factor 4 leads to bad results\n", + "Y=80 # Dimension Y\n", + "Z=16 # Dimension Z\n", + "N_ECHO=20 # Number of ECHO shifts\n", + "CENTER_SIZE=16\n", + "\n", + "# Create the center \n", + "bart ones 2 $CENTER_SIZE $CENTER_SIZE data/mask_center\n", + "bart transpose 0 2 data/mask_center data/mask_center\n", + "\n", + "# Zero-pad the second (index 1) dimension to the full k-space dimension \n", + "bart resize -c 1 $Y 2 $Z data/mask_center data/mask_full\n", + "\n", + "# Create pattern with undersampling in phase direction by R=4\n", + "bart upat -Y $Y -Z $Z -y $R -z 1 -c 1 data/upat\n", + "\n", + "# Initialize the final combined pattern file\n", + "bart saxpy 1 data/mask_full data/upat data/upat_combined\n", + "\n", + "# Loop over N_ECHO times to apply shifts and join them\n", + "for ((i=1; i" + ] + }, + "execution_count": 248, + "metadata": { + "image/png": { + "width": 1000 + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "!bart transpose 1 2 data/upat_combined data/img\n", + "!bart slice 5 1 data/img data/img\n", + "!bart toimg -w 1.0 data/img fig/upat_combined_mask.png\n", + "Image(\"fig/upat_combined_mask.png\", width=1000)" + ] + }, + { + "cell_type": "markdown", + "id": "d6fd6361-4ff8-4fc8-9b59-c2ccc75aa997", + "metadata": {}, + "source": [ + "### 3.3.1 Uniform Undersampling (in-plane)\n", + "With varying patterns along the echos in-plane (kx ky)" + ] + }, + { + "cell_type": "code", + "execution_count": 238, + "id": "b69cf32c-7712-4018-baa0-7f9646cb2502", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Pattern creation completed. Output file: upat_combined\n" + ] + } + ], + "source": [ + "%%bash\n", + "# Parameters\n", + "R=4 # Undersampling factor\n", + "Y=80 # Dimension Y\n", + "X=80 # Dimension Z\n", + "N_ECHO=20 # Number of ECHO shifts\n", + "CENTER_SIZE=16\n", + "\n", + "# Create the center \n", + "bart ones 2 $Y $CENTER_SIZE data/mask_center\n", + "\n", + "# Zero-pad the second (index 1) dimension to the full k-space dimension \n", + "bart resize -c 1 $Y data/mask_center data/mask_full\n", + "\n", + "# Now the mask needs to be in x-y plane\n", + "#bart transpose 1 2 data/mask_full data/mask_full\n", + "#bart transpose 0 1 data/mask_full data/mask_full\n", + "\n", + "# Create pattern with undersampling in phase direction by R=4\n", + "bart upat -Y $Y -Z $X -y $R -z 1 -c 1 data/upat\n", + "\n", + "# Transpose this so it fits the bart dataformat\n", + "bart transpose 0 2 data/upat data/upat\n", + "\n", + "# Initialize the final combined pattern file\n", + "bart saxpy 1 data/mask_full data/upat data/upat_combined_ip\n", + "\n", + "# Loop over N_ECHO times to apply shifts and join them\n", + "for ((i=1; i" + ] + }, + "execution_count": 240, + "metadata": { + "image/png": { + "width": 500 + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "!bart transpose 0 1 data/upat_combined_ip data/img\n", + "!bart slice 5 1 data/img data/img\n", + "!bart toimg -w 1.0 data/img fig/upat_combined_mask_ip.png\n", + "Image(\"fig/upat_combined_mask_ip.png\", width=500)" + ] + }, + { + "cell_type": "code", + "execution_count": 245, + "id": "fa611a8c-a923-4015-825a-71752aab022a", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Scaling: 0.348982\n" + ] + } + ], + "source": [ + "%%bash\n", + "# For nlinv the echos can't be in the coil dimension (3) also to make it work with the upat_combined\n", + "bart transpose 3 5 data/kSpace_transposed data/kSpace_upat_combined_transposed\n", + "\n", + "bart fmac data/upat_combined data/kSpace_upat_combined_transposed data/kSpace_upat_combined_ip\n", + "\n", + "# Compute pattern (just to check)\n", + "bart pattern data/kSpace_upat_combined_ip data/calc_pattern_ip\n", + "\n", + "bart nlinv -i 30 -x 80:80:16 -S data/kSpace_upat_combined_ip data/nlinv_upat_combined_ip\n", + "# bart fft -i $(bart bitmask 0 1 2) data/kSpace_upat data/fft_upat\n", + "\n", + "# bart transpose 3 5 data/fft_upat data/fft_upat_transposed\n", + "\n", + "# Fit the model:\n", + "bart mobafit -T data/echo_times_final data/nlinv_upat_combined_ip data/fit_upat_combined_ip\n", + "\n", + "# Now some values will be very large so we can apply a threshold to obtain a mask\n", + "MAX_T2=1000\n", + "\n", + "bart threshold -M $MAX_T2 data/fit_upat_combined_ip data/mask_upat_combined_ip\n", + "\n", + "# Multiply the fit with the mask\n", + "bart fmac data/fit_upat_combined_ip data/mask_upat_combined_ip data/fit_mask_upat_combined_ip\n", + "\n", + "# Select slice\n", + "bart slice 6 1 data/fit_mask_upat_combined_ip data/R2_map_upat_combined_ip\n", + "\n", + "# Invert the data to get T2\n", + "bart invert data/R2_map_upat_combined_ip data/T2_map_upat_combined_ip" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c3668de-d997-4a76-b711-c413c7378bca", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.2" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/README.md b/README.md new file mode 100644 index 0000000..481bee6 --- /dev/null +++ b/README.md @@ -0,0 +1,3 @@ +# ESMRMB Educational + +Testing different reconstruction methods. diff --git a/__init__.py b/__init__.py new file mode 100644 index 0000000..e69de29 diff --git a/cfl.py b/cfl.py new file mode 100644 index 0000000..265baa8 --- /dev/null +++ b/cfl.py @@ -0,0 +1,114 @@ +# Copyright 2013-2015. The Regents of the University of California. +# Copyright 2021. Uecker Lab. University Center Göttingen. +# All rights reserved. Use of this source code is governed by +# a BSD-style license which can be found in the LICENSE file. +# +# Authors: +# 2013 Martin Uecker +# 2015 Jonathan Tamir + +from __future__ import print_function +from __future__ import with_statement + +import numpy as np +import mmap +import os + + +def readcfl(name): + # get dims from .hdr + with open(name + ".hdr", "rt") as h: + h.readline() # skip + l = h.readline() + dims = [int(i) for i in l.split()] + + # remove singleton dimensions from the end + n = np.prod(dims) + dims_prod = np.cumprod(dims) + dims = dims[:np.searchsorted(dims_prod, n)+1] + + # load data and reshape into dims + with open(name + ".cfl", "rb") as d: + a = np.fromfile(d, dtype=np.complex64, count=n); + return a.reshape(dims, order='F') # column-major + +def readmulticfl(name): + # get dims from .hdr + with open(name + ".hdr", "rt") as h: + lines = h.read().splitlines() + + index_dim = 1 + lines.index('# Dimensions') + total_size = int(lines[index_dim]) + index_sizes = 1 + lines.index('# SizesDimensions') + sizes = [int(i) for i in lines[index_sizes].split()] + index_dims = 1 + lines.index('# MultiDimensions') + + with open(name + ".cfl", "rb") as d: + a = np.fromfile(d, dtype=np.complex64, count=total_size) + + offset = 0 + result = [] + for i in range(len(sizes)): + dims = ([int(i) for i in lines[index_dims + i].split()]) + n = np.prod(dims) + result.append(a[offset:offset+n].reshape(dims, order='F')) + offset += n + + if total_size != offset: + print("Error") + + return result + + +def writecfl(name, array): + with open(name + ".hdr", "wt") as h: + h.write('# Dimensions\n') + for i in (array.shape): + h.write("%d " % i) + h.write('\n') + + size = np.prod(array.shape) * np.dtype(np.complex64).itemsize + + with open(name + ".cfl", "a+b") as d: + os.ftruncate(d.fileno(), size) + mm = mmap.mmap(d.fileno(), size, flags=mmap.MAP_SHARED, prot=mmap.PROT_WRITE) + if array.dtype != np.complex64: + array = array.astype(np.complex64) + mm.write(np.ascontiguousarray(array.T)) + mm.close() + #with mmap.mmap(d.fileno(), size, flags=mmap.MAP_SHARED, prot=mmap.PROT_WRITE) as mm: + # mm.write(array.astype(np.complex64).tobytes(order='F')) + +def writemulticfl(name, arrays): + size = 0 + dims = [] + + for array in arrays: + size += array.size + dims.append(array.shape) + + with open(name + ".hdr", "wt") as h: + h.write('# Dimensions\n') + h.write("%d\n" % size) + + h.write('# SizesDimensions\n') + for dim in dims: + h.write("%d " % len(dim)) + h.write('\n') + + h.write('# MultiDimensions\n') + for dim in dims: + for i in dim: + h.write("%d " % i) + h.write('\n') + + size = size * np.dtype(np.complex64).itemsize + + with open(name + ".cfl", "a+b") as d: + os.ftruncate(d.fileno(), size) + mm = mmap.mmap(d.fileno(), size, flags=mmap.MAP_SHARED, prot=mmap.PROT_WRITE) + for array in arrays: + if array.dtype != np.complex64: + array = array.astype(np.complex64) + mm.write(np.ascontiguousarray(array.T)) + mm.close() diff --git a/plotBoMap.py b/plotBoMap.py new file mode 100644 index 0000000..72a63cd --- /dev/null +++ b/plotBoMap.py @@ -0,0 +1,48 @@ +import scipy.io as sio +import matplotlib.pyplot as plt +import numpy as np +import matplotlib.gridspec as gridspec +import matplotlib.cm as cm +from mpl_toolkits.axes_grid1 import make_axes_locatable + +# Bo Map +mat_data_0=sio.loadmat('bo.mat') +tmap = mat_data_0['boMap'] +print(np.min(tmap)) +print(np.max(tmap)) +vmin = -0.3 +vmax = 0.3 + +def multiSlicePlot(data,nRow, nCol, nSlrep, vmin, vmax, sliceInit = 0): + """"" + data = imageMatrix3D [sl,ph,rd] + nRow, nCol dimensions of the multislicePlot + nSlRep = number of slices to represent + sliceInit = slice in wihich we start representing + """"" + images = [] + fig = plt.figure(figsize=(nRow, nCol), dpi=500) + gs1 = gridspec.GridSpec(nRow, nCol) + gs1.update(wspace=0.020, hspace=0.020) # set the spacing between axes. + for i in range(nSlrep): + ii = i + sliceInit + # print(ii) + ax1 = plt.subplot(gs1[i]) + # plt.axis('off') + ax1.set_xticklabels([]) + ax1.set_yticklabels([]) + ax1.set_aspect('equal') + dataAux = data[int(ii), :, :] + imgPlot = ax1.imshow(dataAux,vmin=vmin, vmax=vmax) + images.append(imgPlot) + ax1.axis('off') + return fig + +fig = multiSlicePlot(tmap*1e3,3,6,16,vmin,vmax,0) + +cbar_ax = fig.add_axes([0.92, 0.15, 0.01, 0.7]) +norm = plt.Normalize(vmin=vmin, vmax=vmax) +cbar = plt.colorbar(cm.ScalarMappable(norm=norm, cmap=cm.viridis), cax=cbar_ax) +cbar.set_label('Bo [mT]', fontsize=2) +cbar.ax.tick_params(labelsize=2) +plt.show() \ No newline at end of file diff --git a/plotComparativeTimes.py b/plotComparativeTimes.py new file mode 100644 index 0000000..192a90d --- /dev/null +++ b/plotComparativeTimes.py @@ -0,0 +1,39 @@ +import scipy.io as sio +import matplotlib.pyplot as plt +import numpy as np +import matplotlib.gridspec as gridspec +import matplotlib.cm as cm +from mpl_toolkits.axes_grid1 import make_axes_locatable + +imagesPath=['T2_CPMG','T2_APCPMG', 'T2_APCP', 'T2_CP', 'T1'] +sliceRep = 12 +vmin = 0 +vmax = 200 + +nMaps = len(imagesPath) +fig = plt.figure(figsize=(1, 5), dpi=500) +gs1 = gridspec.GridSpec(1, nMaps) +gs1.update(wspace=0.020, hspace=0.020) +for i in range(nMaps): + if 'T1' in imagesPath[i]: + mat_data = sio.loadmat(imagesPath[i]+'.mat') + data = mat_data['t1map'] + print(data.shape) + else: + mat_data = sio.loadmat(imagesPath[i]+'.mat') + data = mat_data['t2map'] + print(data.shape) + ax1 = plt.subplot(gs1[i]) + ax1.set_xticklabels([]) + ax1.set_yticklabels([]) + ax1.set_aspect('equal') + dataAux = data[sliceRep, :, :] + imgPlot = ax1.imshow(abs(dataAux),vmin=vmin, vmax=vmax) + ax1.axis('off') + ax1.set_title(imagesPath[i], fontsize=4) +cbar_ax = fig.add_axes([0.92, 0.35, 0.01, 0.30]) +norm = plt.Normalize(vmin=vmin, vmax=vmax) +cbar = plt.colorbar(cm.ScalarMappable(norm=norm, cmap=cm.viridis), cax=cbar_ax) +cbar.set_label('Time [ms]', fontsize=2) +cbar.ax.tick_params(labelsize=2) +plt.show() \ No newline at end of file diff --git a/plotPDMap.py b/plotPDMap.py new file mode 100644 index 0000000..3a0e716 --- /dev/null +++ b/plotPDMap.py @@ -0,0 +1,46 @@ +import scipy.io as sio +import matplotlib.pyplot as plt +import numpy as np +import matplotlib.gridspec as gridspec +import matplotlib.cm as cm +from mpl_toolkits.axes_grid1 import make_axes_locatable + +# PD +mat_data_0=sio.loadmat('PD.mat') +img = mat_data_0['image3D'] +vmin = np.min(np.abs(img)) +vmax = np.max(np.abs(img)) + +def multiSlicePlot(data,nRow, nCol, nSlrep, vmin, vmax, sliceInit = 0): + """"" + data = imageMatrix3D [sl,ph,rd] + nRow, nCol dimensions of the multislicePlot + nSlRep = number of slices to represent + sliceInit = slice in wihich we start representing + """"" + images = [] + fig = plt.figure(figsize=(nRow, nCol), dpi=500) + gs1 = gridspec.GridSpec(nRow, nCol) + gs1.update(wspace=0.020, hspace=0.020) # set the spacing between axes. + for i in range(nSlrep): + ii = i + sliceInit + # print(ii) + ax1 = plt.subplot(gs1[i]) + # plt.axis('off') + ax1.set_xticklabels([]) + ax1.set_yticklabels([]) + ax1.set_aspect('equal') + dataAux = data[int(ii), :, :] + imgPlot = ax1.imshow(dataAux,vmin=vmin, vmax=vmax) + images.append(imgPlot) + ax1.axis('off') + return fig + +fig = multiSlicePlot(np.abs(img),3,6,16,vmin,vmax,0) + +cbar_ax = fig.add_axes([0.92, 0.15, 0.01, 0.7]) +norm = plt.Normalize(vmin=vmin, vmax=vmax) +cbar = plt.colorbar(cm.ScalarMappable(norm=norm, cmap=cm.viridis), cax=cbar_ax) +cbar.set_label('PD [a.u.]', fontsize=2) +cbar.ax.tick_params(labelsize=2) +plt.show() \ No newline at end of file diff --git a/plotT1Map.py b/plotT1Map.py new file mode 100644 index 0000000..42c7996 --- /dev/null +++ b/plotT1Map.py @@ -0,0 +1,12 @@ +import scipy.io as sio +import matplotlib.pyplot as plt +import numpy as np +from bart import bart + +# T1 dict +T1_dict = sio.loadmat('T1.mat') + +rawName = 'T2_CPMG.mat' +mat_data_0=sio.loadmat(rawName) + +print(mat_data_0.keys()) \ No newline at end of file diff --git a/plotTimesMap.py b/plotTimesMap.py new file mode 100644 index 0000000..db13876 --- /dev/null +++ b/plotTimesMap.py @@ -0,0 +1,52 @@ +import scipy.io as sio +import matplotlib.pyplot as plt +import numpy as np +import matplotlib.gridspec as gridspec +import matplotlib.cm as cm +from mpl_toolkits.axes_grid1 import make_axes_locatable + +# Map +rawName = 'T2_CPMG.mat' +mat_data_0=sio.loadmat(rawName) +if 'T1' in rawName: + tmap = mat_data_0['t1map'] +else: + tmap = mat_data_0['t2map'] +print(np.min(tmap)) +print(np.max(tmap)) +vmin = 0 +vmax = 200 + +def multiSlicePlot(data,nRow, nCol, nSlrep, vmin, vmax, sliceInit = 0): + """"" + data = imageMatrix3D [sl,ph,rd] + nRow, nCol dimensions of the multislicePlot + nSlRep = number of slices to represent + sliceInit = slice in wihich we start representing + """"" + images = [] + fig = plt.figure(figsize=(nRow, nCol), dpi=500) + gs1 = gridspec.GridSpec(nRow, nCol) + gs1.update(wspace=0.020, hspace=0.020) # set the spacing between axes. + for i in range(nSlrep): + ii = i + sliceInit + # print(ii) + ax1 = plt.subplot(gs1[i]) + # plt.axis('off') + ax1.set_xticklabels([]) + ax1.set_yticklabels([]) + ax1.set_aspect('equal') + dataAux = data[int(ii), :, :] + imgPlot = ax1.imshow(dataAux,vmin=vmin, vmax=vmax) + images.append(imgPlot) + ax1.axis('off') + return fig + +fig = multiSlicePlot(tmap,3,6,16,vmin,vmax,0) + +cbar_ax = fig.add_axes([0.92, 0.15, 0.01, 0.7]) +norm = plt.Normalize(vmin=vmin, vmax=vmax) +cbar = plt.colorbar(cm.ScalarMappable(norm=norm, cmap=cm.viridis), cax=cbar_ax) +cbar.set_label(rawName+ '[ms]', fontsize=2) +cbar.ax.tick_params(labelsize=2) +plt.show() \ No newline at end of file diff --git a/rawDatasExplanation.txt b/rawDatasExplanation.txt new file mode 100644 index 0000000..7b3cd1a --- /dev/null +++ b/rawDatasExplanation.txt @@ -0,0 +1,36 @@ + +# T2 Maps +t2maps_dict = { + 'kSpaces3D': sampled, # kSpace [kRd, kPh, kSl, kSpace_echo_1, kSpace_echo_2, ..., kSpace_echo_nETL] + 'images3D': imgMSE, # images [nSl, nPh, nETL, nRD] + 'echoSpacing': echoSpacing, # Echo spacing + 'FOV': fov, # FOV [sl,ph,rd] + 't2map': t2Map # t2Map (ms) [nSl, nPh, nRD] +} + +# T1 Map +t1map_dict = { + 'kSpaces3D': sampled, # kSpace [kRd, kPh, kSl, kSpace_echo_1, kSpace_echo_2, ..., kSpace_echo_nETL] + 'images3D': imgAll, # images [nSl, nPh, nImg, nRD] + 'inversionTimes': tI_vector, # Inversion times corresponding each image + 'FOV': fov, # FOV [sl,ph,rd] + 't1map': t1Map # t2Map (ms) [nSl, nPh, nRD] +} + +# Bo Map +bomap_dict = { + 'kSpace3D_0': sampled_0, # kSpace delay 0 [kRd, kPh, kSl, kSpace] + 'image3D_0': img0, # images delay 0 [nSl, nPh, nRD] + 'kSpace3D_delay': sampled_delay, # kSpace delay [kRd, kPh, kSl, kSpace] + 'image3D_delay': img80, # images delay [nSl, nPh, nRD] + 'delay': delay, # delay (s) + 'FOV': fov, # FOV [sl,ph,rd] + 'boMap': boCeros # boMap (T) [nSl, nPh, nRD] +} + +# PD +pd_dict = { + 'kSpace3D': sampled_0, # kSpace [kRd, kPh, kSl, kSpace] + 'image3D': img0, # images [nSl, nPh, nRD] + 'FOV': fov, # FOV [sl,ph,rd] +} \ No newline at end of file diff --git a/reco.py b/reco.py new file mode 100644 index 0000000..413f5bc --- /dev/null +++ b/reco.py @@ -0,0 +1,78 @@ +import scipy.io as sio +import matplotlib.pyplot as plt +import numpy as np +from scipy.optimize import curve_fit +import sys,os +from cfl import writecfl + +os.environ['TOOLBOX_PATH'] = '/home/jpfitzer/bart-0.9.00' +os.environ['BART_TOOLBOX_PATH'] = '/home/jpfitzer/bart-0.9.00' +sys.path.append('/home/jpfitzer/bart-0.9.00/bart/python') + +rawName = 'T2_CPMG.mat' +mat_data_0=sio.loadmat(rawName) + +# Format: sampled, # kSpace [kRd, kPh, kSl, kSpace_echo_1, kSpace_echo_2, ..., kSpace_echo_nETL] (102400, 23) +# So the first three are the coordinates of the kspace, and the rest are the echoes +print(mat_data_0['kSpaces3D'].shape) + +kSpaces3D = mat_data_0['kSpaces3D'] +# self.mapVals['sampled'] = np.concatenate((kRD, kPH, kSL, dataAll_sampled), axis=1) + +# nReadout, nPhase, nSlice +nPoints = (80, 80, 16) + +echo_train_length = mat_data_0['kSpaces3D'].shape[1] - 3 # Because the first 3 are kRD, kPH, kSL -> should give 20 +print(f"Echo train length: {echo_train_length}") + +echo_spacing = mat_data_0['echoSpacing'][0][0] +print(f"Echo spacing: {echo_spacing}") + +k_readout = kSpaces3D[:, 0] +k_phase = kSpaces3D[:, 1] +k_slice = kSpaces3D[:, 2] + +# The rest of the data is the echoes +echos = kSpaces3D[:, 3:] + +# Reshape the kspace data for bart +kSpace = echos.reshape(nPoints[2], nPoints[1], nPoints[0], echo_train_length) + + + +print(kSpace.shape) + +# cfl = writecfl('kSpace', kSpace) + +# Create the image with bart fft -i 7 kSpace fft -> three dimensional +# Put the echos on the fifth dimension: +# bart transpose 3 5 +# Put the slices on the correct dimension +# bart transpose 0 2 + +# traj = writecfl('traj', np.stack((k_readout, k_phase, k_slice), axis=1)) + +# Echo times with echo spacing +TE = np.linspace(echo_spacing, echo_spacing * echo_train_length, echo_train_length, endpoint=True) +print("TE: ", TE) + +# Create the echotimes file: +# bart vec ... echo_times +# bart scale 0.001 echo_times echo_times_scaled +# Move the echo_times to the correct dimension: +# bart transpose 0 5 echo_times_scaled echo_times_final + +# Fit the model: +# bart mobafit -T echo_times_final fft_transposed fit + +# Now some values will be very large so we can apply a threshold to obtain a mask +# bart threshold -M 1000 reco/fit reco/mask + +# Multiply the fit with the mask +# bart fmac reco/fit reco/fit reco/fit_mask + +# Select slice +# bart slice 6 1 reco/fit_mask reco/R2_map + +# Invert the data to get T2 +# bart invert reco/R2_map reco/T2_map \ No newline at end of file