{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Neuroimaging non-cartesian reconstruction\n\nAuthor: Chaithya G R\n\nIn this tutorial we will reconstruct an MRI image from non-cartesian kspace\nmeasurements.\n\n## Import neuroimaging data\n\nWe use the toy datasets available in pysap, more specifically a 2D brain slice\nand the acquisition cartesian scheme.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Package import\nfrom mri.operators import NonCartesianFFT, WaveletUD2\nfrom mri.operators.utils import convert_locations_to_mask, \\\n    gridded_inverse_fourier_transform_nd\nfrom mri.reconstructors import SingleChannelReconstructor\nimport pysap\nfrom pysap.data import get_sample_data\n\n# Third party import\nfrom modopt.math.metrics import ssim\nfrom modopt.opt.linear import Identity\nfrom modopt.opt.proximity import SparseThreshold\nimport numpy as np\n\n# Loading input data\nimage = get_sample_data('2d-mri')\n\n# Obtain MRI non-cartesian mask\nradial_mask = get_sample_data(\"mri-radial-samples\")\nkspace_loc = radial_mask.data\nmask = pysap.Image(data=convert_locations_to_mask(kspace_loc, image.shape))\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 radial 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 and the associated observations\nfourier_op = NonCartesianFFT(samples=kspace_loc, shape=image.shape,\n                             implementation='cpu')\nkspace_obs = fourier_op.op(image.data)\n\n# Gridded solution\ngrid_space = np.linspace(-0.5, 0.5, num=image.shape[0])\ngrid2D = np.meshgrid(grid_space, grid_space)\ngrid_soln = gridded_inverse_fourier_transform_nd(kspace_loc, kspace_obs,\n                                                 tuple(grid2D), 'linear')\nimage_rec0 = pysap.Image(data=grid_soln)\n# image_rec0.show()\nbase_ssim = ssim(image_rec0, image)\nprint('The Base SSIM is : ' + str(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 = WaveletUD2(\n    wavelet_id=24,\n    nb_scale=4,\n)\nregularizer_op = SparseThreshold(Identity(), 6 * 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)\nx_final, costs, metrics = reconstructor.reconstruct(\n    kspace_data=kspace_obs,\n    optimization_alg='fista',\n    num_iterations=200,\n)\nimage_rec = pysap.Image(data=np.abs(x_final))\n# image_rec.show()\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
}