webinar3/subspace_T1_exercise/exercise_subspace_T1_bart.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: \n",
    "## Subspace-constrained Reconstruction for Inversion-prepared T1 Mapping in BART\n",
    "\n",
    "This exercise uses the BART command-line inteface (CLI) (http://mrirecon.github.io/bart/) and challenges you to apply the previously presented subspace-constrained reconstruction for inversion-prepared T1 Mapping using BART.\n",
    "\n",
    "**Main Reference**\n",
    "\n",
    "Wang X, Tan Z, Scholand N, Roeloffs V, Uecker M. Physics-based reconstruction methods for magnetic resonance imaging. accepted to Philos. Trans. R. Soc. A. arXiv:2010.01403.\n",
    "\n",
    "**Authors**: [Xiaoqing Wang](mailto:xiaoqing.wang@med.uni-goettingen.de), [Nick Scholand](mailto:nick.scholand@med.uni-goettingen.de)\n",
    "\n",
    "**Institution**: University Medical Center Göttingen"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Prerequisites\n",
    "\n",
    "Make sure the BART configuration is completed and these paths are set correctly.\n",
    "\n",
    "`TOOLBOX_PATH=/Your/path/to/bart`\n",
    "\n",
    "`PATH=$TOOLBOX_PATH:$PATH`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we need the visualization functions already defined in the main tutorials."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "echo \"BART version\"\n",
    "\n",
    "bart version\n",
    "\n",
    "# BASH function to write BART file to png,\n",
    "# display the png, and then remove it\n",
    "function imshow () {\n",
    "    bart toimg -W $1 $1.png > /dev/null\n",
    "    cat $1.png | display\n",
    "}\n",
    "export -f imshow\n",
    "\n",
    "function imshow_colormap () {\n",
    "    python3 python/save_maps.py $1 viridis 0 2.0 $1.png > /dev/null\n",
    "    cat $1.png | display\n",
    "}\n",
    "export -f imshow_colormap\n",
    "\n",
    "# BASH function to visualize timesteps of BART cfl file\n",
    "function show_timesteps () {\n",
    "\n",
    "    data=$1\n",
    "    shift\n",
    "    pos=(\"$@\")\n",
    "\n",
    "    ind=0\n",
    "\n",
    "    for t in \"${pos[@]}\"\n",
    "    do\n",
    "       bart slice 5 $t $data _slice$ind\n",
    "       ind=$((ind+1))\n",
    "    done\n",
    "\n",
    "    bart join 6 `seq -f \"_slice%g\" 0 $((ind-1))` _slices\n",
    "    \n",
    "    DIM=`bart show -d 0 _slices`\n",
    "\n",
    "    bart reshape $(bart bitmask 1 6) $((DIM*ind)) 1 {_,}slices\n",
    "\n",
    "    imshow slices\n",
    "}\n",
    "export -f show_timesteps"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Download Dataset\n",
    "\n",
    "Many publications using the BART toolbox provide tutorials, code and data. Check out the [BART homepage](https://mrirecon.github.io/bart/) for further details.\n",
    "Here, we make use one of the published datasets and download an inversion-recovery FLASH one already converted to a *.cfl* and *.hdr* file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wget -q https://zenodo.org/record/4060287/files/IR-FLASH.cfl\n",
    "wget -q https://zenodo.org/record/4060287/files/IR-FLASH.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us check the dimensions of the dataset by printing out the first two lines of the provided *.hdr* file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head -n2 IR-FLASH.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the dimensionality definitions in BART, we can read out and store the information about the acquisition for later use in the processing and reconstruction."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "READ=`bart show -d 0 IR-FLASH`\n",
    "SPOKES=`bart show -d 2 IR-FLASH`\n",
    "COILS=`bart show -d 3 IR-FLASH`\n",
    "TIME=`bart show -d 10 IR-FLASH`\n",
    "BR=$((READ/2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Dictionary Generation, SVD and Temporal Basis\n",
    "\n",
    "The exercise starts by creating the dictionary. \n",
    "\n",
    "First, we specify the sequence characteristics necessary for the simulation of the signal curves."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define Sequence Characteristics\n",
    "TR=0.0041 # [s]\n",
    "REP=$TIME\n",
    "\n",
    "# Dictionary characteristics\n",
    "## R1s\n",
    "NUM_R1S=1000\n",
    "MIN_R1S=5e-3\n",
    "MAX_R1S=5\n",
    "\n",
    "## Mss\n",
    "NUM_MSS=100\n",
    "MIN_MSS=1e-2\n",
    "MAX_MSS=1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use the `signal` tool in BART to simulate an analytical **inversion-prepared FLASH** signal by passing the `-F` and `-I` flag. Pass the number of **repetitions** with `-n` and the **repetition time** with `-r`. \n",
    "\n",
    "Pass the desired number of $R_1^*$ values with their minimum and maximum values exploiting the `-1 min:max:N` flag and repeat the same for $M_{ss}$ and the `-3 min:max:N` flag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Simulate dictionary based on the `signal` tool\n",
    "bart signal ...your-code-here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Have a look at the dimensions of the created dictionary by printing the corresponding *.hdr* file or using BARTs `show` tool."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print dimensions of created dictionary\n",
    "head -n2 dicc.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The range of simulated $R_1^*$ values is stored in the 6th (coeff1) dimension while the $M_{ss}$ values are in the 7th (coeff2) dimension.\n",
    "\n",
    "To process the dictionary all elements need to be stored in the 6th dimension. Therefore, use the `reshape` command to concentrate all simulated signals in the coeff1 dimension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# reshape the dicc 6th and 7th dimension to have all the elements \n",
    "# concentrated in the 6th dimension\n",
    "bart reshape ...your-code-here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Print the dimensions of the reshaped dictionary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the dimensions of the reshaped dictionary\n",
    "head -n2 dicc_reshape.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now squeeze all non-zero dimensions with the `squeeze` command and check you progress by printing the dimensions again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# squeeze the reshaped dictionary to remove non-zero dimensions\n",
    "bart squeeze dicc_reshape dicc_squeeze"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform an SVD on the squeezed dictionary with the `svd` command. Specify to perform an economy-size decomposition by adding the `-e` flag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform an SVD of the squeezed dictionary and output and\n",
    "# economy-size decomposition of the resulting matrix\n",
    "bart svd ...your-code-here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can create the temporal basis by extracting the desired number of coefficients from the columns in the 1st dimension of $U$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define number of desired temporal basis coefficients\n",
    "NUM_COE=4\n",
    "\n",
    "# Extract desired number of orthonormal columns from U\n",
    "bart extract ...your-code-here...\n",
    "\n",
    "# Print the basis dimensions\n",
    "head -n2 basis.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "lines_to_next_cell": 0
   },
   "source": [
    "As you might have seen in the dimensions of our createed basis, the coefficients are currently in the 1st and the time is in the 0th dimension.\n",
    "\n",
    "For compatibility with BART we need to move the dimensions to their natural position within the 16 dimensions of a *.cfl* file:\n",
    "\n",
    "    time          -> 5th (echo time) dimension\n",
    "    coeffiecients -> 6th (coeff1) dimension.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# transpose the basis to have time in the 5th dimension \n",
    "# and coefficients in the 6th dimension\n",
    "bart transpose 1 6 basis basis1\n",
    "bart transpose 0 5 basis1 basis_$NUM_COE\n",
    "\n",
    "# Print the transposed basis dimensions\n",
    "head -n2 basis_$NUM_COE.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coil Compression, Generate Golden-Angle-based Radial Trajectory and Gradient Delay Correction \n",
    "\n",
    "Print the dimensions of the previously downloaded dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "head -n2 IR-FLASH.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the non-Cartesian data handling expects cooridnate positions in the 0th dimension and samples in the first dimension, we need to transpose the 0th and 1st dimension with the `transpose` command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transpose the 0th and 1st dimension of the downloaded data\n",
    "# to ensure compatibility with BARTs non-Cartesian tools\n",
    "bart transpose 0 1 IR-FLASH ksp\n",
    "\n",
    "# Print dimensions\n",
    "head -n2 ksp.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Coil Compression\n",
    "To reduce the size of our dataset and therefore also decrease the computational complexity, we perform a coil compression with the `cc` command. By passing `-A` we choose to use all possible data and want to reduce the dataset to *NUM_VCOILS* virtual coils with `-p`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NUM_VCOILS=8\n",
    "\n",
    "# coil compression\n",
    "bart cc ...your-code-here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Generate the Trajectory\n",
    "In the next step we generate a trajectory with the `traj` tool. To match the acquisition of the downloaded data, we need to specify a radial `-r`, centered `-c`, double-angle `-D`, 7th tiny golden-angle `-G -s7` sampling. The timesteps are passed using `-t`, the spokes by `-y` and the samples are specified with `-x`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "NUM_TGA=7\n",
    "\n",
    "# Create the trajectory using the `traj` tool\n",
    "bart traj ...your-code-here...\n",
    "\n",
    "# Print out its dimensions\n",
    "head -n2 traj.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here the 3 in the zeroth dimensions includes the coordinates in (x,y,z) and is the reasion for the previous transpose of the downloaded dataset!\n",
    "\n",
    "\n",
    "#### Gradient Delay Correction\n",
    "Because the dataset is following an IR FLASH signal relaxation, the gradient delay correction should be applied to data in the last repetitions which are in a steady-state. Therefore, we extract some final repetitions from the trajectory and the dataset using the `extract` command. Have in mind that the time dimension is the 10th here!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the number of timesteps for gradient delay correction\n",
    "FRAMES_FOR_GDC=100\n",
    "\n",
    "# Extract the DATA_GDC last time frames from the trajectory and the dataset\n",
    "bart extract 10 $((REP-FRAMES_FOR_GDC)) $REP traj traj_extract\n",
    "bart extract 10 $((REP-FRAMES_FOR_GDC)) $REP ksp_cc ksp_extract"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The gradient delay tool `estdelay` only estimates delays between spokes and does not support a time dimension. Assuming that all extracted spokes were acquired during the steady-state, all extracted timesteps can be transposed into the 2nd dimension for spokes (phase2) using the `transpose` command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Transpose the 2nd and 10th dimension for later use with the `estdelay` tool\n",
    "bart transpose 10 2 traj_extract traj_extract1\n",
    "bart transpose 10 2 ksp_extract ksp_extract1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Afterwards the gradient delays can be estimated using the `estdelay` tool. It is recommended to use the RING technique by adding the `-R` flag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate and store the gradient delays usign RING\n",
    "GDELAY=$(bart estdelay -R traj_extract1 ksp_extract1)\n",
    "\n",
    "echo \"Gradient Delays: \"$GDELAY"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To correct the trajectory by this delay, we use the same `traj` command call as before, but pass the known delays with the `-q`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the trajectory with known gradient delays\n",
    "bart traj ...your-code-here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you do not want to run the trajectory creation, the following cell will load the precomputed ouput."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load precomputed results if necessary\n",
    "[ ! -f \"trajn.cfl\" ] && bart copy ref/trajn_precomp trajn"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coil Sensitivity Estimation\n",
    "\n",
    "The coil profile estimation is similar to the gradient delay estimation performed on some of the last timesteps of the IR FLASH dataset. Therefore, *FRAMES_FOR_CSE* spokes from the last timesteps are extracted using the `extract` command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "FRAMES_FOR_CSE=300\n",
    "\n",
    "# Extract last time frames (10th dim) from trajectory and k-space data\n",
    "bart extract ...your-code-here...\n",
    "bart extract ...your-code-here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By assuming all extracted frames to be acquired during the steady-state, we again can transpose the 10th dimension (time) with the spoke dimension (2nd/phase2)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bart transpose 2 10 traj_ss traj_ss2\n",
    "bart transpose 2 10 ksp_cc_ss ksp_cc_ss2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the coil profile estimation we want to use **ESPIRiT**. Because its BART implementation `ecalib` requires Cartesian k-space data, we need to grid our radial data first.\n",
    "\n",
    "A handy solution in BART exploits the gridding step during an nuFFT call. So we start with performing an inverse `nufft` on our extracted steady-state data. Constrain the output dimension using the `-d` flag to the readout samples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# apply an inverse nufft of the extracted steady-state dataset\n",
    "bart nufft ...your-code-here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting reconstruction can now be transformed back into the frequency domain with an unitary (`-u`) Fourier tranformation with the `fft` command."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# transform reconstruction in image space back to k-space\n",
    "# to create gridded k-space for the ESPIRiT implementation\n",
    "bart fft ...your-code-here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, use BARTs ESPIRiT implementation `ecalib` to estimate the coil sensitivities from the gridded k-space. We recommend using a threshold `-t`, just one ESPIRit map `-m` and smooth map transitions `-S`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "THRESHOLD=0.01\n",
    "NUM_ESPIRIT_MAP=1\n",
    "\n",
    "# Estimate coil sensitivities from gridded, steady-state k-space using `ecalib`\n",
    "bart ecalib ...your-code-here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you do not want to run the coil sensitivity estimation, the following line will load the precomputed sensitivities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load precomputed results if necessary\n",
    "[ ! -f \"sens_invivo.cfl\" ] && bart copy ref/sens_invivo_precomp sens_invivo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Subspace-Constrained Reconstruction\n",
    "To start the subspace-constrained reconstruction we need to verify the data and trajectory dimension again."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "echo \"Trajectory:\"\n",
    "head -n2 trajn.hdr\n",
    "\n",
    "echo \"Data:\"\n",
    "head -n2 ksp_cc.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `pics` tool in BART requires all timesteps to be in the 5th (echo) dimension. Use the `reshape` or `transpose` command to do this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bart transpose 5 10 trajn traj_final\n",
    "bart transpose 5 10 ksp_cc ksp_final"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we can perform the subspace-constrained reconstruction with BARTs **P**arallel **I**maging **C**ompressed **S**ensing tool: `pics`.\n",
    "\n",
    "We recommend scaling on maximum eigenvalues `-e` and $l_1$-wavelet regularization `-R W:A:B:C `in the image dimensions (A) and jointly (B) on the basis coefficients. You also might use a debugg value `-d` and pass a previously defined number of iterations `-i`.\n",
    "\n",
    "Most important for a subspace-constrained reconstruction is to pass the temporal basis with `-B`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DEBUG=4\n",
    "ITER=100\n",
    "REG=0.0015\n",
    "\n",
    "bart pics ...your-code-here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you do not want to run the time consuming reconstruction, the following cell will load the precomputed results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load precomputed results if necessary\n",
    "[ ! -f \"subspace_reco_invivo.cfl\" ] && \\\n",
    "bart copy ref/subspace_reco_invivo_precomp subspace_reco_invivo"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you can see from the dimensions of the reconstruction"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print dimensions of reconstruction\n",
    "head -n2 subspace_reco_invivo.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "we got a two times larger FOV through the 2-fold oversampling. The images can be resized to the desired base resolution *BR* using the `resize` command. Ensure to use a centered crop `-c` of the first two dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Resize reconstructions to remove 2-fold oversampling effects\n",
    "bart resize ...your-code-here..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can visualize the reconstructed maps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reshape and flip coefficient maps for improved visualization\n",
    "\n",
    "## Concentrate all coefficients in the column dimension (1st/phase1)\n",
    "bart reshape $(bart bitmask 1 6) $((BR*NUM_COE)) 1 coeff_maps subspace_maps\n",
    "\n",
    "## Flip the map in row dimension to have the forhead pointing to the top of the page\n",
    "bart flip $(bart bitmask 0) subspace_maps subspace_maps1\n",
    "\n",
    "imshow subspace_maps1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Projection to Time Series and T1 Map Estimation\n",
    "\n",
    "The reconstructed maps are in a subspace. To obtain the full time series from these coefficients we multiply the coefficient maps with the previously estimated basis using `fmac`. After summing the result using the `-s` flag in the coefficient dimension, we have recovered a full time series from the subspace coefficients."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiply the reconstructed coeff maps with the optained time basis\n",
    "# Squash the coeffiencient dimension (6) to sum all coeff timeseries up\n",
    "bart fmac ...your-code-here...\n",
    "\n",
    "# Print dimensions of time series\n",
    "head -n2 time_series.hdr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we visualize some exemplaric steps of the time series."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Flip the datas row dimension again to have the forehead pointing up\n",
    "bart flip $(bart bitmask 0) time_series time_series1\n",
    "\n",
    "# Visualize some exemplaric timesteps\n",
    "t=( 0, 100, 200, 300, 900 )\n",
    "show_timesteps time_series1 \"${t[@]}\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To estimate T1 maps from these timesteps we will perform a pixelwise fitting using the pyhton script provided in *mapping_piecewise.py*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform a pixelwise fitting of the time series to obtain the T1 maps\n",
    "MTR=4100 #TR in [us]\n",
    "M0_THRESHOLD=0.1\n",
    "DELAY=15.3e-3\n",
    "\n",
    "## Prepare time vector\n",
    "### Create vector from 0 to REP\n",
    "bart index 5 $REP tmp1\n",
    "### Multiply vector with TR in [us]\n",
    "bart scale $MTR tmp1 tmp2\n",
    "### Scale time vector from [us] to [s]\n",
    "bart scale 0.000001 tmp2 TI\n",
    "    \n",
    "# Perform pixelwise fitting\n",
    "python3 python/mapping_piecewise.py time_series1 T1 TI maps\n",
    "\n",
    "# Extract first three maps (Mss, M0, R1s)...\n",
    "bart extract 2 0 3 maps tmp1\n",
    "# ... and move them to the coeff dimension\n",
    "bart transpose 2 6 tmp1 tmp2\n",
    "\n",
    "# Estimate T1 from (Mss, M0, R1s)\n",
    "bart looklocker -t $M0_THRESHOLD -D $DELAY tmp2 tmp3\n",
    "\n",
    "# Correct T1 for reconstruction scaling factor\n",
    "bart scale 0.5 tmp3 subspace_invivo_T1map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you do not want to run the pixel-wise fitting the following cell will just load the precomputed result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load precomputed results if necessary\n",
    "[ ! -f \"subspace_invivo_T1map.cfl\" ] && \\\n",
    "bart copy ref/subspace_invivo_T1map_precomp subspace_invivo_T1map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fianlly, we can visualize the estimated T1 map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imshow_colormap subspace_invivo_T1map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thank you for your participation! Feel free to ask any questions:)"
   ]
  }
 ],
 "metadata": {
  "jupytext": {
   "formats": "ipynb,md"
  },
  "kernelspec": {
   "display_name": "Bash",
   "language": "bash",
   "name": "bash"
  },
  "language_info": {
   "codemirror_mode": "shell",
   "file_extension": ".sh",
   "mimetype": "text/x-sh",
   "name": "bash"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}