{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Neuroimaging cartesian reconstruction\n\nAuthor: Chaithya G R\n\nIn this tutorial we will reconstruct an MRI image from the sparse kspace\nmeasurements.\n\n## Import neuroimaging data\n\nWe use the toy datasets available in pysap, more specifically a 2D brain slice\nand the cartesian acquisition scheme.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Package import\nfrom mri.operators import FFT, WaveletN\nfrom mri.operators.utils import convert_mask_to_locations\nfrom mri.reconstructors import SingleChannelReconstructor\nimport pysap\nfrom pysap.data import get_sample_data\n\n# Third party import\nfrom modopt.opt.proximity import SparseThreshold\nfrom modopt.opt.linear import Identity\nfrom modopt.math.metrics import ssim\nimport numpy as np\n\n# Loading input data\nimage = get_sample_data('2d-mri')\n\n# Obtain K-Space Cartesian Mask\nmask = get_sample_data(\"cartesian-mri-mask\")\n\n# View Input\n# image.show()\n# mask.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## Generate the kspace\n\nFrom the 2D brain slice and the acquisition mask, we retrospectively\nundersample the k-space using a cartesian acquisition mask\nWe then reconstruct the zero order solution as a baseline\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Get the locations of the kspace samples\nkspace_loc = convert_mask_to_locations(mask.data)\n# Generate the subsampled kspace\nfourier_op = FFT(samples=kspace_loc, shape=image.shape)\nkspace_data = fourier_op.op(image)\n\n# Zero order solution\nimage_rec0 = pysap.Image(data=fourier_op.adj_op(kspace_data),\n                         metadata=image.metadata)\n# image_rec0.show()\n\n# Calculate SSIM\nbase_ssim = ssim(image_rec0, image)\nprint(base_ssim)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "## FISTA optimization\n\nWe now want to refine the zero order solution using a FISTA optimization.\nThe cost function is set to Proximity Cost + Gradient Cost\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Setup the operators\nlinear_op = WaveletN(wavelet_name=\"sym8\", nb_scales=4)\nregularizer_op = SparseThreshold(Identity(), 2 * 1e-7, thresh_type=\"soft\")\n# Setup Reconstructor\nreconstructor = SingleChannelReconstructor(\n    fourier_op=fourier_op,\n    linear_op=linear_op,\n    regularizer_op=regularizer_op,\n    gradient_formulation='synthesis',\n    verbose=1,\n)\n# Start Reconstruction\nx_final, costs, metrics = reconstructor.reconstruct(\n    kspace_data=kspace_data,\n    optimization_alg='fista',\n    num_iterations=200,\n)\nimage_rec = pysap.Image(data=np.abs(x_final))\n# image_rec.show()\n# Calculate SSIM\nrecon_ssim = ssim(image_rec, image)\nprint('The Reconstruction SSIM is : ' + str(recon_ssim))"
      ]
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}