esmrmb-educational/ESMRMB_Educational.ipynb

{
 "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.. 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": "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/jupyter-jupfi/esmrmb-educational/bart//bart\n",
      "# BART version: \n",
      "v0.9.00-466-gb3be61a\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "echo \"# The BART used in this notebook:\"\n",
    "which bart\n",
    "echo \"# BART version: \"\n",
    "bart version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "4bfab7fe-902a-4aa4-8748-757fa25b269e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Virtual environment already exists.\n",
      "Activating virtual environment and installing packages...\n",
      "Requirement already satisfied: numpy in ./venv/lib/python3.10/site-packages (2.1.0)\n",
      "Requirement already satisfied: scipy in ./venv/lib/python3.10/site-packages (1.14.0)\n",
      "Requirement already satisfied: matplotlib in ./venv/lib/python3.10/site-packages (3.9.2)\n",
      "Requirement already satisfied: pyparsing>=2.3.1 in ./venv/lib/python3.10/site-packages (from matplotlib) (3.1.2)\n",
      "Requirement already satisfied: cycler>=0.10 in ./venv/lib/python3.10/site-packages (from matplotlib) (0.12.1)\n",
      "Requirement already satisfied: pillow>=8 in ./venv/lib/python3.10/site-packages (from matplotlib) (10.4.0)\n",
      "Requirement already satisfied: kiwisolver>=1.3.1 in ./venv/lib/python3.10/site-packages (from matplotlib) (1.4.5)\n",
      "Requirement already satisfied: packaging>=20.0 in ./venv/lib/python3.10/site-packages (from matplotlib) (24.1)\n",
      "Requirement already satisfied: python-dateutil>=2.7 in ./venv/lib/python3.10/site-packages (from matplotlib) (2.9.0.post0)\n",
      "Requirement already satisfied: contourpy>=1.0.1 in ./venv/lib/python3.10/site-packages (from matplotlib) (1.2.1)\n",
      "Requirement already satisfied: fonttools>=4.22.0 in ./venv/lib/python3.10/site-packages (from matplotlib) (4.53.1)\n",
      "Requirement already satisfied: six>=1.5 in ./venv/lib/python3.10/site-packages (from python-dateutil>=2.7->matplotlib) (1.16.0)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import subprocess\n",
    "\n",
    "# Function to execute shell command within Jupyter notebook\n",
    "def run_command(command):\n",
    "    process = subprocess.Popen(command, shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE)\n",
    "    stdout, stderr = process.communicate()\n",
    "    return process.returncode, stdout.decode(), stderr.decode()\n",
    "\n",
    "# Check if venv exists\n",
    "venv_path = './venv'\n",
    "if not os.path.exists(venv_path):\n",
    "    # Create venv if it doesn't exist\n",
    "    print(\"Creating virtual environment...\")\n",
    "    result, out, err = run_command('python3 -m venv venv')\n",
    "    if result != 0:\n",
    "        print(f\"Error creating virtual environment: {err}\")\n",
    "\n",
    "else:\n",
    "    print(\"Virtual environment already exists.\")\n",
    "\n",
    "# Activate the virtual environment and install packages\n",
    "activate_command = '. venv/bin/activate && pip install numpy scipy matplotlib'\n",
    "\n",
    "print(\"Activating virtual environment and installing packages...\")\n",
    "result, out, err = run_command(activate_command)\n",
    "if result != 0:\n",
    "    print(f\"Error activating virtual environment or installing packages: {err}\")\n",
    "else:\n",
    "    print(out)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "1575b526-7b4e-43f7-85d8-9c7d2b2a30e9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Python Stuff\n",
    "from bart import bart\n",
    "import cfl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "237a793d-f70e-430c-a471-38d3b427ee26",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Bart viewer\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib import cm\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import ipywidgets as widgets\n",
    "from IPython.display import display\n",
    "from traitlets import List, Bool\n",
    "\n",
    "class DimensionSelector(widgets.HBox):\n",
    "    value = List(Bool()).tag(sync=True)\n",
    "\n",
    "    def __init__(self, n_dims):\n",
    "        self.checkboxes = [widgets.Checkbox(value=i<2, description=f'Dim {i}') for i in range(n_dims)]\n",
    "        super().__init__(self.checkboxes)\n",
    "        for cb in self.checkboxes:\n",
    "            cb.observe(self.on_change, 'value')\n",
    "        self.value = [cb.value for cb in self.checkboxes]\n",
    "\n",
    "    def on_change(self, change):\n",
    "        selected = [cb for cb in self.checkboxes if cb.value]\n",
    "        if len(selected) > 2:\n",
    "            selected[0].value = False\n",
    "        self.value = [cb.value for cb in self.checkboxes]\n",
    "\n",
    "def interactive_bart_plot(data, cmap='gray', cbar_label='', mag=True, fsize=10):\n",
    "    data = np.array(data)\n",
    "    if mag:\n",
    "        data = np.abs(data)\n",
    "\n",
    "    def update_plot(plot_dims, interpolation, **slice_values):\n",
    "        plt.clf()\n",
    "        selected_dims = [i for i, sel in enumerate(plot_dims) if sel]\n",
    "        if len(selected_dims) != 2:\n",
    "            plt.text(0.5, 0.5, 'Please select exactly two dimensions for plotting', \n",
    "                     ha='center', va='center', fontsize=12)\n",
    "            plt.axis('off')\n",
    "        else:\n",
    "            # Create a slice object for indexing\n",
    "            index = [slice(None)] * data.ndim\n",
    "            for i in range(data.ndim):\n",
    "                if i not in selected_dims:\n",
    "                    index[i] = slice_values[f'slice_dim{i}']\n",
    "            \n",
    "            # Extract the 2D slice to plot\n",
    "            plot_data = data[tuple(index)]\n",
    "            \n",
    "            fig, ax1 = plt.subplots(figsize=(fsize, fsize))\n",
    "            im = ax1.imshow(plot_data, interpolation=interpolation, cmap=cmap)\n",
    "            \n",
    "            divider = make_axes_locatable(ax1)\n",
    "            cax = divider.append_axes(\"right\", size=\"5%\", pad=0.05)\n",
    "            cbar = plt.colorbar(im, cax=cax)\n",
    "            cbar.set_label(cbar_label)\n",
    "            \n",
    "            ax1.set_title(f'Dimensions: {selected_dims}, Interpolation: {interpolation}')\n",
    "            ax1.set_axis_off()\n",
    "        plt.show()\n",
    "\n",
    "    dim_selector = DimensionSelector(data.ndim)\n",
    "    \n",
    "    def create_slice_widget(i):\n",
    "        widget = widgets.IntText(value=0, description=f'Slice Dim {i}:', \n",
    "                                 layout=widgets.Layout(width='200px'))\n",
    "        \n",
    "        def on_value_change(change):\n",
    "            if change['new'] < 0:\n",
    "                widget.value = 0\n",
    "            elif change['new'] >= data.shape[i]:\n",
    "                widget.value = data.shape[i] - 1\n",
    "        \n",
    "        widget.observe(on_value_change, names='value')\n",
    "        return widget\n",
    "\n",
    "    slice_widgets = {f'slice_dim{i}': create_slice_widget(i) for i in range(data.ndim)}\n",
    "    \n",
    "    interpolation_options = ['nearest', 'bilinear', 'bicubic', 'spline16', 'spline36', 'hanning', 'hamming', 'hermite', 'kaiser', 'quadric', 'catrom', 'gaussian', 'bessel', 'mitchell', 'sinc', 'lanczos']\n",
    "    interpolation_widget = widgets.Dropdown(\n",
    "        options=interpolation_options,\n",
    "        value='nearest',\n",
    "        description='Interpolation:',\n",
    "        layout=widgets.Layout(width='200px')\n",
    "    )\n",
    "    \n",
    "    slice_ui = widgets.VBox(list(slice_widgets.values()))\n",
    "    ui = widgets.VBox([dim_selector, interpolation_widget, slice_ui])\n",
    "    \n",
    "    controls = {'plot_dims': dim_selector, 'interpolation': interpolation_widget, ** slice_widgets}\n",
    "    out = widgets.interactive_output(update_plot, controls)\n",
    "    \n",
    "    display(ui, out)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "cb006433-1224-4cf6-9cc3-40ae0fdd03fe",
   "metadata": {},
   "outputs": [],
   "source": [
    "# BartIO (used for bitmask)\n",
    "import io\n",
    "import contextlib\n",
    "\n",
    "def bartIO(*args, **kwargs):\n",
    "    buffer = io.StringIO()\n",
    "\n",
    "    with contextlib.redirect_stdout(buffer):\n",
    "        bart(*args, **kwargs)\n",
    "\n",
    "    return buffer.getvalue()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "9953cfd1-750d-4c44-aa50-bd3c59d4c276",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f5012800346c43b19e6fd17e6d9918e0",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "01519e2b671a47a780900969da42362e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Testing the viewer\n",
    "bart_phantom = bart(1, 'phantom -x 256 -B').T\n",
    "interactive_bart_plot(bart_phantom)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77d0c22c-bb16-4a5d-916b-0d4263fc8c50",
   "metadata": {},
   "source": [
    "## 2. Reading in the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "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": 11,
   "id": "845409eb-b38f-4b1f-9cbc-5fef3b16b778",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7886bfd3e2f243a2b51ef26822aa93a4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0e47ed3ea7244e17a7500836cacf3320",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interactive_bart_plot(kSpace)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "150faf49-acd6-4ee3-ad36-3a745e8e88d9",
   "metadata": {},
   "source": [
    "### 2.1 Data Exploration"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5f93fcca-48ab-4d42-9061-91361aadfb35",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "kSpace dimensions:\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": [
    "print(\"kSpace dimensions:\")\n",
    "bart(0, 'show -m', kSpace)\n",
    "\n",
    "DIM_X = 80\n",
    "DIM_Y = 80\n",
    "DIM_Z = 16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "77d36186-04bb-422d-9ec3-a2032fcbceda",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "New Shape:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "(80, 80, 16, 1, 1, 20)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Correct dimension\n",
    "kSpace_transposed = bart(1, 'transpose 0 2', kSpace)\n",
    "kSpace_transposed = bart(1, 'transpose 3 5', kSpace_transposed)\n",
    "print(\"New Shape:\")\n",
    "kSpace_transposed.shape"
   ]
  },
  {
   "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": 14,
   "id": "44062fa0-0d88-4957-a4a0-1e3c36a3048b",
   "metadata": {},
   "outputs": [],
   "source": [
    "bitmask= bartIO(0, \"bitmask 0 1 2\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "c074e817-5989-4e83-a108-79834c7b0bdf",
   "metadata": {},
   "outputs": [],
   "source": [
    "fft = bart(1, f'fft -i {int(bitmask)}', kSpace_transposed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "fcfc949d-8046-4868-bdfe-810b5f27454c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "644f112a16be41b9b53171303f0d47e8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cef42e6a597744c2a9d60d86f4778475",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interactive_bart_plot(fft)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "1e3cdf7c-8b6c-43f9-91d2-abe353b47ab7",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Echo times:\n",
      " [0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1  0.11 0.12 0.13 0.14\n",
      " 0.15 0.16 0.17 0.18 0.19 0.2 ]\n"
     ]
    }
   ],
   "source": [
    "echo_times = np.linspace(10, 200, num=echo_train_length) \n",
    "echo_times *= 0.001\n",
    "print(f'Echo times:\\n {echo_times}')\n",
    "echo_times_final = bart(1, 'transpose 0 5', echo_times)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "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 <file>] [--init [f:]*f] [--scale [f:]*f] [--min-flag d] [--max-flag d] [--max-mag-flag d] [--min [f:]*f] [--max [f:]*f] <enc> <echo/contrast images> [<coefficients>] \n",
      "\n",
      "Pixel-wise fitting of physical signal models.\n",
      "\n",
      "-T                      Multi-Echo Spin Echo: f(M0, R2) = M0 * exp(-t * R2)\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": [
    "! bart mobafit -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e0581e77-f92f-4b2f-8e06-d9eb8f1604a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "R2 = bart(1, 'mobafit -T', echo_times_final, fft)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "dce7fefb-88cd-443f-a85a-37ac7ec374d7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(80, 80, 16, 1, 1, 1, 2)"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "R2.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "0b59f1a1-60fe-4c97-9d7d-35e5180216fd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a2b8770756584b68afda965bde8cae0a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0101c630434d44a5b2047e1dc5631893",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interactive_bart_plot(R2)\n",
    "# Oh noes! The fit  doesn't work :("
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "8d35a117-a8f7-4a82-beb0-1dc29cb58e2a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Usage: threshold [-H] [-W] [-L] [-D] [-B] [-M] [-j d] [-b d] lambda <input> <output> \n",
      "\n",
      "Perform (soft) thresholding with parameter lambda.\n",
      "\n",
      "-H              hard thresholding\n",
      "-W              daubechies wavelet soft-thresholding\n",
      "-L              locally low rank soft-thresholding\n",
      "-D              divergence-free wavelet soft-thresholding\n",
      "-B              thresholding with binary output where (val>lambda)\n",
      "-M              thresholding with binary output where (val<lambda)\n",
      "-j bitmask      joint soft-thresholding\n",
      "-b blocksize    locally low rank block size\n",
      "-h              help\n"
     ]
    }
   ],
   "source": [
    "! bart threshold -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "1e1083db-f618-4aca-9f83-6f7b44b50f9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "MAX_T2 = 1\n",
    "\n",
    "T2 = bart(1, 'invert', R2)\n",
    "\n",
    "mask = bart(1, f'threshold -M {MAX_T2}', T2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "e7e30fdf-dd59-4bb6-8b86-4c0618dc0003",
   "metadata": {},
   "outputs": [],
   "source": [
    "T2_masked = T2 * mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "82477c08-8ad9-44eb-aa1f-5e95e9790142",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "933878e7ef0445978c79e0a2c133d6b6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7cfc58d2a91d4b80ad124ae90c4508e8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interactive_bart_plot(T2_masked, cmap='viridis') #(yay)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9b49f422-7373-4f2e-8994-872c4bcd1622",
   "metadata": {},
   "source": [
    "## 2.3 PICS Reco Fully Sampled Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "id": "9fd005c5-1ce2-4de6-a80d-1cbec18856b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "ecalib_full = bart(1, 'ecalib -m1', kSpace_transposed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 76,
   "id": "d51dc743-7697-4228-b593-11471e1375b1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "78e86e4d76a54be382e20d3fbafac351",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "565d2706a76645b3815576c68f77d98f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interactive_bart_plot(ecalib_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "id": "54508b05-8865-4713-990f-8565e624f559",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Usage: pics [-l ...] [-r f] [-R ...] [-c] [-s f] [-i d] [-t <file>] [-n] [-N] [-g] [--gpu-gridding] [-p <file>] [--precond] [-b d] [-e] [-W <file>] [-d d] [-u f] [-C d] [-f f] [-I,--ist] [--fista] [--eulermaruyama] [-m,--admm] [-a,--pridu] [-w f] [-S] [--shared-img-dims d] [-K] [-B <file>] [-P f] [-M] [-U,--lowmem] [--no-toeplitz] [--psf_export <file>] [--psf_import <file>] [--wavelet <string>] [--mpi d] [--fista_pqr f:f:f] [--fista_last] <kspace> <sensitivities> <output> \n",
      "\n",
      "Parallel-imaging compressed-sensing reconstruction.\n",
      "\n",
      "\n",
      "-l1/-l2                    toggle l1-wavelet or l2 regularization.\n",
      "-r lambda                  regularization parameter\n",
      "-R <T>:A:B:C               generalized regularization options (-Rh for help)\n",
      "-c                         real-value constraint\n",
      "-s step                    iteration stepsize\n",
      "-i iter                    max. number of iterations\n",
      "-t file                    k-space trajectory\n",
      "-n                         disable random wavelet cycle spinning\n",
      "-N                         do fully overlapping LLR blocks\n",
      "-g                         use GPU\n",
      "--gpu-gridding             use GPU for gridding\n",
      "-p file                    pattern or weights\n",
      "--precond                  interprete weights as preconditioner\n",
      "-b blk                     Lowrank block size\n",
      "-e                         Scale stepsize based on max. eigenvalue\n",
      "-W <img>                   Warm start with <img>\n",
      "-d level                   Debug level\n",
      "-u rho                     ADMM rho\n",
      "-C iter                    ADMM max. CG iterations\n",
      "-f rfov                    restrict FOV\n",
      "-I,--ist                   select IST\n",
      "--fista                    select FISTA\n",
      "--eulermaruyama            select Euler Maruyama\n",
      "-m,--admm                  select ADMM\n",
      "-a,--pridu                 select Primal Dual\n",
      "-w f                       inverse scaling of the data\n",
      "-S                         re-scale the image after reconstruction\n",
      "--shared-img-dims flags    deselect image dims with flags\n",
      "-K                         randshift for NUFFT\n",
      "-B file                    temporal (or other) basis\n",
      "-P eps                     Basis Pursuit formulation, || y- Ax ||_2 <= eps\n",
      "-M                         Simultaneous Multi-Slice reconstruction\n",
      "-U,--lowmem                Use low-mem mode of the nuFFT\n",
      "--no-toeplitz              Turn off Toeplitz mode of nuFFT\n",
      "--psf_export file          Export PSF to file\n",
      "--psf_import file          Import PSF from file\n",
      "--wavelet name             wavelet type (haar,dau2,cdf44)\n",
      "--mpi flags                distribute over this dimensions with use of MPI\n",
      "--fista_pqr p:q:r          parameters for FISTA acceleration\n",
      "--fista_last               end iteration with call to proximal op\n",
      "-h                         help\n"
     ]
    }
   ],
   "source": [
    "! bart pics -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "id": "5b1ce368-549c-4b78-bbae-d2392f9fdf57",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Size: 2048000 Samples: 2047999 Acc: 1.00\n",
      "l1 regularization: 0.030000\n",
      "Regularization terms: 1, Supporting variables: 0\n",
      "FISTA\n",
      "Total Time: 2.629528\n"
     ]
    }
   ],
   "source": [
    "bitmask = bartIO(0, 'bitmask 5')\n",
    "\n",
    "REGULARIZATION = 3e-2\n",
    "ITERATIONS=50\n",
    "\n",
    "pics_fully = bart(1, f'pics -RI:{int(bitmask)}:{REGULARIZATION} -i{ITERATIONS} -S', kSpace_transposed, ecalib_full)\n",
    "# pics_fully = bart(1, f'pics -l1 -r{REGULARIZATION} -i{ITERATIONS} -S', kSpace_transposed, ecalib_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "id": "eb3bec34-95a0-4419-8675-c04c2e16e8c1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1701098fdbf74bab9e8ff53378a6efb2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "8a74c508e5494351bdbfa116c327b69a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the pics reco\n",
    "interactive_bart_plot(pics_fully)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "id": "5c4b28cf-efaf-4535-ba12-a1301bca327d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "201b5688d03a428e94f657cae9b627ce",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "65f1ea0d04e54ab4ad568a3ea5e8dc3c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "R2_pics = bart(1, 'mobafit -T', echo_times_final, pics_fully)\n",
    "\n",
    "MAX_T2 = 1\n",
    "\n",
    "T2_pics = bart(1, 'invert', R2_pics)\n",
    "\n",
    "mask_pics = bart(1, f'threshold -M {MAX_T2}', T2_pics)\n",
    "\n",
    "T2_pics_masked = T2_pics * mask_pics\n",
    "\n",
    "interactive_bart_plot(T2_pics_masked, cmap='viridis')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f8d7d303-2747-4f41-bf5d-153bc83f2382",
   "metadata": {},
   "source": [
    "## 3.1 Undersampling Poisson ky - kz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "id": "46bb3f9a-aab7-4f9f-ab61-8c44856968b1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "points: 451, grid size: 80x16x(pi/4) = 1005 (R = 2.229068)\n"
     ]
    }
   ],
   "source": [
    "ACC_Y=2\n",
    "ACC_Z=2\n",
    "CENTER=16\n",
    "\n",
    "poisson_mask = bart(1, f'poisson -Y {kSpace_transposed.shape[1]} -Z {kSpace_transposed.shape[2]} -y {ACC_Y} -z {ACC_Z} -C {CENTER} -e')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "id": "1ca63cbe-772f-4799-8c19-d3e376bc1809",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "66fb3a2c1ffc4c7596fd98a1cf9e279e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "322c5e2785cd49d580c93160932b2026",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "interactive_bart_plot(poisson_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "id": "1eefa822-903f-47b5-8f75-ea857760a683",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dac913f003284bfbb617efdcba0fdde7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6850d70d2bec4b33836dee9c304646f4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Apply the mask\n",
    "kSpace_poisson = bart(1, 'fmac', poisson_mask, kSpace_transposed)\n",
    "interactive_bart_plot(kSpace_poisson)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "id": "e41fc27c-c713-4b26-ab0a-8841297cea04",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Usage: nlinv [-i d] [-d d] [-c] [-N] [-m d] [-U] [-f f] [-p <file>] [-t <file>] [-B <file>] [-I <file>] [-g] [-S] [--lowmem] [-x,--dims d:d:d] [--sens-dims d:d:d] [--real-time] [--fast] <kspace> <output> [<sensitivities>] \n",
      "\n",
      "Jointly estimate image and sensitivities with nonlinear\n",
      "inversion using {iter} iteration steps. Optionally outputs\n",
      "the sensitivities.\n",
      "\n",
      "-i iter              Number of Newton steps\n",
      "-d level             Debug level\n",
      "-c                   Real-value constraint\n",
      "-N                   Do not normalize image with coil sensitivities\n",
      "-m nmaps             Number of ENLIVE maps to use in reconstruction\n",
      "-U                   Do not combine ENLIVE maps in output\n",
      "-f FOV               restrict FOV\n",
      "-p file              pattern / transfer function\n",
      "-t file              kspace trajectory\n",
      "-B file              temporal (or other) basis\n",
      "-I file              File for initialization\n",
      "-g                   use gpu\n",
      "-S                   Re-scale image after reconstruction\n",
      "--lowmem             Use low-mem mode of the nuFFT\n",
      "-x,--dims x:y:z      Explicitly specify image dimensions\n",
      "--sens-dims x:y:z    Explicitly specify sens dimensions\n",
      "--real-time          Use real-time (temporal l2) regularization\n",
      "--fast               Use tuned but less generic model\n",
      "-h                   help\n"
     ]
    }
   ],
   "source": [
    "! bart nlinv -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "id": "a0727645-10b5-4c0b-b42b-560057531a83",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Model created (Cartesian):\n",
      "kspace:  [ 80  80  16   1   1  20   1   1   1   1   1   1   1   1   1   1 ]\n",
      "images:  [ 80  80  16   1   1  20   1   1   1   1   1   1   1   1   1   1 ]\n",
      "coils:   [ 80  80  16   1   1  20   1   1   1   1   1   1   1   1   1   1 ]\n",
      "pattern: [ 80  80  16   1   1  20   1   1   1   1   1   1   1   1   1   1 ]\n",
      "coilimg: [ 80  80  16   1   1  20   1   1   1   1   1   1   1   1   1   1 ]\n",
      "\n",
      "Scaling: 0.349091\n",
      "Step: 0, Res: 100.000001\n",
      "\t cg:   0\n",
      "Step: 1, Res: 99.959831\n",
      "\t cg:   0\n",
      "Step: 2, Res: 99.998280\n",
      "\t cg:   0\n",
      "Step: 3, Res: 99.999653\n",
      "\t cg:   0\n",
      "Step: 4, Res: 99.999687\n",
      "\t cg:   0\n",
      "Step: 5, Res: 99.998778\n",
      "\t cg:   0\n",
      "Step: 6, Res: 99.977362\n",
      "\t cg:   1\n",
      "Step: 7, Res: 99.069636\n",
      "\t cg:   2\n",
      "Step: 8, Res: 72.770714\n",
      "\t cg:   8\n",
      "Step: 9, Res: 60.890033\n",
      "\t cg:   1\n",
      "Step: 10, Res: 43.049279\n",
      "\t cg:  33\n",
      "Step: 11, Res: 136.236493\n",
      "\t cg:   0\n",
      "Step: 12, Res: 44.616604\n",
      "\t cg:   1\n",
      "Step: 13, Res: 14.514742\n",
      "\t cg:  16\n",
      "Step: 14, Res: 23.674170\n",
      "\t cg:   1\n",
      "Step: 15, Res: 4.198526\n",
      "\t cg:   6\n",
      "Step: 16, Res: 0.311999\n",
      "\t cg:  64\n",
      "Step: 17, Res: 16.272082\n",
      "\t cg:   4\n",
      "Step: 18, Res: 2.207014\n",
      "\t cg:   7\n",
      "Step: 19, Res: 0.049158\n",
      "\t cg:  44\n",
      "Step: 20, Res: 5.464199\n",
      "\t cg:   4\n",
      "Step: 21, Res: 0.183185\n",
      "\t cg:  51\n",
      "Step: 22, Res: 4.384588\n",
      "\t cg:   4\n",
      "Step: 23, Res: 0.102021\n",
      "\t cg:  48\n",
      "Step: 24, Res: 3.743119\n",
      "\t cg:   4\n",
      "Total time: 127.88 s\n"
     ]
    }
   ],
   "source": [
    "#Nlinv\n",
    "ITERATIONS=25\n",
    "nlinv_poisson = bart(1, f'nlinv -i{ITERATIONS} -x 80:80:16 -S -d4', kSpace_poisson) # This takes some time (5min)..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "id": "4307a819-22a8-42c0-a217-9b080ab00242",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a12fd5ce9a0b41a1bc613bc4be9006d9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9cfba7a6a35141f282f49a22860bc906",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "R2_nlinv_poisson = bart(1, 'mobafit -T', echo_times_final, nlinv_poisson)\n",
    "\n",
    "MAX_T2 = 1\n",
    "\n",
    "T2_nlinv_poisson = bart(1, 'invert', R2_nlinv_poisson)\n",
    "\n",
    "mask_nlinv_poisson = bart(1, f'threshold -M {MAX_T2}', T2_nlinv_poisson)\n",
    "\n",
    "T2_nlinv_poisson_masked = T2_nlinv_poisson * mask_nlinv_poisson\n",
    "\n",
    "interactive_bart_plot(T2_nlinv_poisson_masked, cmap='viridis')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8da50e61-6996-44c9-85d0-e97cef04f4c3",
   "metadata": {},
   "source": [
    "## 3.2 Uniform Undersampling\n",
    "Using bart upat - same pattern for every echo."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "id": "f7ed59e6-3998-4a15-9c43-d18dfe0f9211",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c5f68e117ec84de4bc4fb4aed78e4444",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dc1e29a9c7904e559d3a9e7f2c73a445",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "ACC_Y=4\n",
    "ACC_Z=4\n",
    "CENTER=9\n",
    "\n",
    "upat_mask = bart(1, f'upat -Y {kSpace_transposed.shape[1]} -Z {kSpace_transposed.shape[2]} -y {ACC_Y} -z {ACC_Z} -c {CENTER}')\n",
    "\n",
    "interactive_bart_plot(upat_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "id": "6d4c9a4d-b0d5-427a-b0f8-2d153d564f23",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "028d4e9c75d545668e872c0681144c73",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "47a02ed8ba394f5293bf0f5adcf82209",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Apply the mask\n",
    "kSpace_upat = bart(1, 'fmac', upat_mask, kSpace_transposed)\n",
    "interactive_bart_plot(kSpace_upat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "id": "1848b13f-50fa-4fcf-b480-ac3a441eb72f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Size: 2048000 Samples: 531200 Acc: 3.86\n",
      "l1 regularization: 0.100000\n",
      "Regularization terms: 1, Supporting variables: 0\n",
      "FISTA\n",
      "Total Time: 2.628610\n"
     ]
    }
   ],
   "source": [
    "# PICS reco because its faster - using ecalib from fully sampled data?\n",
    "\n",
    "bitmask = bartIO(0, 'bitmask 5')\n",
    "\n",
    "REGULARIZATION = 10e-2\n",
    "ITERATIONS=50\n",
    "\n",
    "pics_upat = bart(1, f'pics -RI:{int(bitmask)}:{REGULARIZATION} -i{ITERATIONS} -S', kSpace_upat, ecalib_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 90,
   "id": "8434c6e3-6a66-430f-95d6-6d86e57ba6ae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "d59d8a6a78ed4c7693c2e2e7360194fe",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "422d5da8e41e4b8f84981033361e4b81",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "R2_pics_upat = bart(1, 'mobafit -T', echo_times_final, pics_upat)\n",
    "\n",
    "MAX_T2 = 1\n",
    "\n",
    "T2_pics_upat = bart(1, 'invert', R2_pics_upat)\n",
    "\n",
    "mask_pics_upat = bart(1, f'threshold -M {MAX_T2}', T2_pics_upat)\n",
    "\n",
    "T2_pics_upat_masked = T2_pics_upat * mask_pics_upat\n",
    "\n",
    "interactive_bart_plot(T2_pics_upat_masked, cmap='viridis')"
   ]
  },
  {
   "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": 91,
   "id": "9c6fc6bd-d9b1-43f5-81e6-f8732d995668",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "cat << 'EOF' > create_pattern.sh\n",
    "#!/bin/bash\n",
    "\n",
    "create_pattern() {\n",
    "    # Default values for parameters\n",
    "    local R_Y=4\n",
    "    local R_Z=1\n",
    "    local Y=80\n",
    "    local Z=16\n",
    "    local N_ECHO=20\n",
    "    local CENTER_SIZE=16\n",
    "    local OUTPUT_NAME='upat_combined_new'\n",
    "\n",
    "    # Parse named arguments\n",
    "    while [[ \"$#\" -gt 0 ]]; do\n",
    "        case $1 in\n",
    "            -R_Y) R_Y=\"$2\"; shift 2 ;;\n",
    "            -R_Z) R_Z=\"$2\"; shift 2 ;;\n",
    "            -Y) Y=\"$2\"; shift 2 ;;\n",
    "            -Z) Z=\"$2\"; shift 2 ;;\n",
    "            -N_ECHO) N_ECHO=\"$2\"; shift 2 ;;\n",
    "            -CENTER_SIZE) CENTER_SIZE=\"$2\"; shift 2 ;;\n",
    "            -OUTPUT_NAME) OUTPUT_NAME=\"$2\"; shift 2 ;;\n",
    "            *) echo \"Unknown parameter passed: $1\"; return 1 ;;\n",
    "        esac\n",
    "    done\n",
    "\n",
    "    # Check if all required parameters are provided\n",
    "    if [ -z \"$R_Y\" ] || [ -z \"$R_Z\" ] || [ -z \"$Y\" ] || [ -z \"$Z\" ] || [ -z \"$N_ECHO\" ] || [ -z \"$CENTER_SIZE\" ] || [ -z \"$OUTPUT_NAME\" ]; then\n",
    "        echo \"Usage: create_pattern -R_Y <value> -R_Z <value> -Y <value> -Z <value> -N_ECHO <value> -CENTER_SIZE <value> -OUTPUT_NAME <value>\"\n",
    "        return 1\n",
    "    fi\n",
    "\n",
    "    # Create a temporary directory\n",
    "    local TMP_DIR=$(mktemp -d -t bart-XXXXXXX)\n",
    "\n",
    "    mkdir -p data\n",
    "\n",
    "    # Create the center\n",
    "    bart ones 2 $CENTER_SIZE $CENTER_SIZE $TMP_DIR/mask_center_new\n",
    "    bart transpose 0 2 $TMP_DIR/mask_center_new $TMP_DIR/mask_center_new_transposed\n",
    "\n",
    "    # Zero-pad the second (index 1) dimension to the full k-space dimension \n",
    "    bart resize -c 1 $Y 2 $Z $TMP_DIR/mask_center_new_transposed $TMP_DIR/mask_full_new\n",
    "\n",
    "    # Create pattern with undersampling in phase and slice directions by R_Y and R_Z\n",
    "    bart upat -Y $Y -Z $Z -y $R_Y -z $R_Z -c 1 $TMP_DIR/upat_new\n",
    "\n",
    "    # Initialize the final combined pattern file\n",
    "    bart saxpy 1 $TMP_DIR/mask_full_new $TMP_DIR/upat_new $TMP_DIR/upat_combined_new\n",
    "\n",
    "    # Loop over N_ECHO times to apply shifts and join them\n",
    "    for ((i=1; i<N_ECHO; i++)); do\n",
    "        # Apply a circular shift of one along the phase direction\n",
    "        bart circshift 1 $i $TMP_DIR/upat_new $TMP_DIR/upat_shifted_new\n",
    "\n",
    "        # Combine the center mask and the shifted pattern\n",
    "        bart saxpy 1 $TMP_DIR/mask_full_new $TMP_DIR/upat_shifted_new $TMP_DIR/tmp\n",
    "\n",
    "        # Join the shifted pattern along the 5th dimension \n",
    "        bart join 5 $TMP_DIR/upat_combined_new $TMP_DIR/tmp $TMP_DIR/upat_combined_new\n",
    "    done\n",
    "\n",
    "    bart threshold -B 0.1 $TMP_DIR/upat_combined_new data/\"$OUTPUT_NAME\"\n",
    "\n",
    "    # Clean up temporary directory\n",
    "    rm -rf $TMP_DIR\n",
    "\n",
    "    echo \"Pattern creation completed. Output file: $OUTPUT_NAME\"\n",
    "}\n",
    "EOF\n",
    "\n",
    "# Make the script executable\n",
    "chmod +x create_pattern.sh\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "id": "6d3961d2-9825-49a5-9d25-f44f0710d132",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pattern creation completed. Output file: upat_var\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "# Source the script to make the function available\n",
    "source ./create_pattern.sh\n",
    "\n",
    "# Parameters\n",
    "R_Y=4        # Undersampling factor 4 leads to bad results\n",
    "R_Z=1\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_pattern -R_Y $R_Y -R_Z $R_Z -Y $Y -Z $Z -N_ECHO $N_ECHO -CENTER_SIZE $CENTER_SIZE -OUTPUT_NAME \"upat_var\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "id": "bc553987-83fa-41c5-a4ca-54562477b616",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5d143663b26344eea1b952024af17332",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "46bd6b3df2064237a1c1f2f82f2b621d",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "upat_var = cfl.readcfl(\"data/upat_var\")\n",
    "interactive_bart_plot(upat_var) # it's a little hard to see ... echos are in dimension 5 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "id": "8985435a-d56b-48f8-98c7-b12f23ff8d5e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0731a18282cc4488bc12c98a63987f07",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "12970750b2cb410e8a72fab221f93837",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Apply the mask\n",
    "kSpace_upat_var = bart(1, 'fmac', upat_var, kSpace_transposed)\n",
    "interactive_bart_plot(kSpace_upat_var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "id": "c4cf0414-5f49-4ecd-9918-283e8d29152c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Size: 2048000 Samples: 819199 Acc: 2.50\n",
      "l1 regularization: 0.100000\n",
      "Regularization terms: 1, Supporting variables: 0\n",
      "FISTA\n",
      "Total Time: 2.655738\n"
     ]
    }
   ],
   "source": [
    "bitmask = bartIO(0, 'bitmask 5')\n",
    "\n",
    "REGULARIZATION = 10e-2\n",
    "ITERATIONS=50\n",
    "\n",
    "pics_upat_var = bart(1, f'pics -RI:{int(bitmask)}:{REGULARIZATION} -i{ITERATIONS} -S', kSpace_upat_var, ecalib_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "id": "ad8c6af0-3f84-4b96-a803-60ea3788c97c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3a864bdcbc494eda8db6b8312bf32f48",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "73929e445ee54a79aa2e6cae730e0435",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "R2_pics_upat_var = bart(1, 'mobafit -T', echo_times_final, pics_upat_var)\n",
    "\n",
    "MAX_T2 = 1\n",
    "\n",
    "T2_pics_upat_var = bart(1, 'invert', R2_pics_upat_var)\n",
    "\n",
    "mask_pics_upat_var = bart(1, f'threshold -M {MAX_T2}', T2_pics_upat_var)\n",
    "\n",
    "T2_pics_upat_var_masked = T2_pics_upat_var * mask_pics_upat_var\n",
    "\n",
    "interactive_bart_plot(T2_pics_upat_var_masked, cmap='viridis') # Undersampling artifacts :("
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16448609-399c-4967-896c-126819be8fec",
   "metadata": {},
   "source": [
    "## 4.0 Subspace Reco\n",
    "First fully for testing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "id": "bc87333c-2620-4f3e-b301-7ced06e19944",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Usage: signal [-F] [-B] [-T] [-S] [-M] [-G] [-C] [--fat] [-I] [-s] [--short-TR-LL-approx] [-0 f:f:f] [-1 f:f:f] [-2 f:f:f] [-3 f:f:f] [-4 f:f:f] [-5 f:f:f] [-r f] [-e f] [-i f] [-f f] [-d f] [-t f] [-n d] [-b d] [--av-spokes d] [-m d] <basis-functions> \n",
      "\n",
      "Analytical simulation tool.\n",
      "\n",
      "-F                      FLASH\n",
      "-B                      bSSFP\n",
      "-T                      TSE\n",
      "-S                      SE\n",
      "-M                      MOLLI\n",
      "-G                      MGRE\n",
      "-C                      IR MGRE\n",
      "--fat                   Simulate additional fat component.\n",
      "-I                      inversion recovery\n",
      "-s                      inversion recovery starting from steady state\n",
      "--short-TR-LL-approx    Short TR approximation for analytical LL model.\n",
      "-0 min:max:N            range of off-resonance frequency [Hz]\n",
      "-1 min:max:N            range of T1s [s]\n",
      "-2 min:max:N            range of T2s [s]\n",
      "-3 min:max:N            range of Mss\n",
      "-4 min:max:N            range of T1 values for fat [s]\n",
      "-5 min:max:N            range of FA values [°]\n",
      "-r TR                   repetition time\n",
      "-e TE                   echo time\n",
      "-i TI                   inversion time\n",
      "-f FA                   flip ange\n",
      "-d FF                   fat fraction\n",
      "-t T1 relax             T1 relax period (second) for MOLLI\n",
      "-n n                    number of measurements\n",
      "-b heart beats          number of heart beats for MOLLI\n",
      "--av-spokes d           Number of averaged consecutive spokes\n",
      "-m multi echos          number of multi gradient echos\n",
      "-h                      help\n"
     ]
    }
   ],
   "source": [
    "!bart signal -h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "id": "f985bfa4-6258-425d-9fbe-b0ef0cc059ad",
   "metadata": {},
   "outputs": [],
   "source": [
    "N_MEASUREMENTS = 20\n",
    "MIN_T2 = 10e-3\n",
    "MAX_T2 = 0.5\n",
    "N_T2_POINTS = 1000\n",
    "ECHO_TIME=10e-3\n",
    "\n",
    "signal = bart(1, f'signal -n{N_MEASUREMENTS} -2 {MIN_T2}:{MAX_T2}:{N_T2_POINTS} -e{ECHO_TIME}')\n",
    "#Should add a function for td plotting\n",
    "\n",
    "#Formatting\n",
    "signal = bart(1, f'transpose 2 5', signal)\n",
    "\n",
    "signal = bart(1, f'squeeze', signal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "id": "ad1ca17d-61cc-4cd8-acf0-682d954c6801",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(20, 1000)\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "c85e5f2495d0406aa65a1af52968edf2",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "59cfaf3dd53b42a195432929cc8a4702",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "print(signal.shape)\n",
    "interactive_bart_plot(signal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 103,
   "id": "36011621-5d6a-4877-b2bd-3f9f3e45b638",
   "metadata": {},
   "outputs": [],
   "source": [
    "U, S, V = bart(3, 'svd -e', signal)\n",
    "\n",
    "# Extract coefficients\n",
    "basis_tmp = bart(1, 'extract 1 0 3', U)\n",
    "basis_tmp = bart(1, 'transpose 1 6', basis_tmp)\n",
    "basis = bart(1, 'transpose 0 5', basis_tmp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 104,
   "id": "df8761cc-1ec3-487a-87e8-3f0142d3a3ed",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Basis: [ 80  80  16   1   1  20   3   1   1   1   1   1   1   1   1   1 ]\n",
      "Max:   [ 80  80  16   1   1   1   3   1   1   1   1   1   1   1   1   1 ]\n",
      "Size: 2048000 Samples: 2047999 Acc: 1.00\n",
      "l1 regularization: 0.030000\n",
      "Regularization terms: 1, Supporting variables: 0\n",
      "FISTA\n",
      "Total Time: 3.235034\n"
     ]
    }
   ],
   "source": [
    "REGULARIZATION=3e-2\n",
    "ITERATIONS=50\n",
    "\n",
    "bitmask = bartIO(0, 'bitmask 5')\n",
    "\n",
    "# Now we get a coefficient map\n",
    "cfl.writecfl(\"data/basis\", basis)\n",
    "coeff_full = bart(1, f'pics -B data/basis -RI:{int(bitmask)}:{REGULARIZATION} -i{ITERATIONS} -S', kSpace_transposed, ecalib_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "id": "5845b065-8226-450a-8fb5-849333bc2ee0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dd2191713095458f95be26f82fe7735e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7a2a987ac2214d7a80b5094fbb9c1dc8",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# From the coefficients and the basis we can get the image\n",
    "bitmask = int(bartIO(0, 'bitmask 6'))\n",
    "              \n",
    "img_subspace_full = bart(1, f'fmac -s {bitmask}', basis, coeff_full)\n",
    "\n",
    "interactive_bart_plot(img_subspace_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "id": "332d6b61-c35d-43b4-80cf-d52bfff3040a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b0201f1b764f4113a0b92c92ab08d010",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cde2edba4ea14a6cb745f3e812d44f8b",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Now the pixelwise fit\n",
    "R2_subspace_full = bart(1, 'mobafit -T', echo_times_final, img_subspace_full)\n",
    "\n",
    "MAX_T2 = 1\n",
    "\n",
    "T2_subspace_full = bart(1, 'invert', R2_subspace_full)\n",
    "\n",
    "mask_subspace_full = bart(1, f'threshold -M {MAX_T2}', T2_subspace_full)\n",
    "\n",
    "T2_subspace_full_masked = T2_subspace_full * mask_subspace_full\n",
    "\n",
    "interactive_bart_plot(T2_subspace_full_masked, cmap='viridis') # Something is off with the scaling"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb1744d3-9ee0-4af0-a3a3-09ddbea9db47",
   "metadata": {},
   "source": [
    "## Two dim USP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "id": "32dae0cd-423a-4a3d-b766-9cde90e961c1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Pattern creation completed. Output file: upat_2x2\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "# Source the script to make the function available\n",
    "source ./create_pattern.sh\n",
    "\n",
    "# Parameters\n",
    "R_Y=2      # Undersampling factor in the Y dimension\n",
    "R_Z=2      # Undersampling factor in the Z dimension\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_pattern -R_Y $R_Y -R_Z $R_Z -Y $Y -Z $Z -N_ECHO $N_ECHO -CENTER_SIZE $CENTER_SIZE -OUTPUT_NAME \"upat_2x2\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "id": "6a539cc0-d5aa-4c3c-be60-f415e1450a93",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "55739de9525c48e98a8e60e14baa66bb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1804c4a15d054991a4467159d24ed5c4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "upat_2x2 = cfl.readcfl(\"data/upat_2x2\")\n",
    "interactive_bart_plot(upat_2x2) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "id": "3912a3c6-eaa6-422f-8811-a6057aec2121",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "7599be5951fe4e6ab727aa0d9468dceb",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4a249ccdd6d040b1a343467ab93257cc",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Apply the pattern\n",
    "kSpace_2x2 = bart(1, 'fmac', upat_2x2, kSpace_transposed)\n",
    "interactive_bart_plot(kSpace_2x2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "id": "f8b24714-18b2-4bac-8578-16c133eeff97",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Basis: [ 80  80  16   1   1  20   3   1   1   1   1   1   1   1   1   1 ]\n",
      "Max:   [ 80  80  16   1   1   1   3   1   1   1   1   1   1   1   1   1 ]\n",
      "Size: 2048000 Samples: 819199 Acc: 2.50\n",
      "l1 regularization: 0.030000\n",
      "Regularization terms: 1, Supporting variables: 0\n",
      "FISTA\n",
      "Total Time: 3.139522\n"
     ]
    }
   ],
   "source": [
    "REGULARIZATION=3e-2\n",
    "ITERATIONS=50\n",
    "\n",
    "bitmask = bartIO(0, 'bitmask 5')\n",
    "\n",
    "coeff_2x2 = bart(1, f'pics -B data/basis -RI:{int(bitmask)}:{REGULARIZATION} -i{ITERATIONS} -S', kSpace_2x2, ecalib_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "id": "28d1efe2-1582-42f0-af6c-331b5f09d985",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fa610e35287443118f4135c3b306d65a",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "97d8c73e2dfe43ce8a9c57f287603acf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# From the coefficients and the basis we can get the image\n",
    "bitmask = int(bartIO(0, 'bitmask 6'))\n",
    "              \n",
    "img_2x2 = bart(1, f'fmac -s {bitmask}', basis, coeff_2x2)\n",
    "\n",
    "interactive_bart_plot(img_2x2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "5f141904-b972-4d26-9a0b-851d421c8a87",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "cabb35a96aa8458882241b511d71bf99",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "VBox(children=(DimensionSelector(children=(Checkbox(value=True, description='Dim 0'), Checkbox(value=True, des…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ca06a09f44e04bbb8b0930d619a14bc7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Output()"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Now the pixelwise fit\n",
    "R2_2x2 = bart(1, 'mobafit -T', echo_times_final, img_2x2)\n",
    "\n",
    "MAX_T2 = 1\n",
    "\n",
    "T2_2x2 = bart(1, 'invert', R2_2x2)\n",
    "\n",
    "mask_2x2 = bart(1, f'threshold -M {MAX_T2}', T2_2x2)\n",
    "\n",
    "T2_2x2_masked = T2_2x2 * mask_2x2\n",
    "\n",
    "interactive_bart_plot(T2_2x2_masked, cmap='viridis') # Something is off with the scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3754d681-ba36-4d7f-9c30-573af2cf624a",
   "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.10.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}