{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Cartesian Calibrationless Reconstruction using GroupLASSO Regularizer\n\nAuthor: Chaithya G R\n\nIn this tutorial we will reconstruct an MRI image from cartesian kspace\nmeasurements.\n\n## Import neuroimaging data\n\nWe use the toy datasets available in pysap, more specifically a 2D parallel MRI\nbrain slice on 32 channels and the acquisition cartesian 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 CalibrationlessReconstructor\nimport pysap\nfrom pysap.data import get_sample_data\n\n# Third party import\nfrom modopt.math.metrics import ssim\nfrom modopt.opt.proximity import GroupLASSO\nimport numpy as np\n\n\n# Loading input data\ncartesian_ref_image = get_sample_data('2d-pmri')\nimage = pysap.Image(data=np.sqrt(np.sum(cartesian_ref_image.data**2, axis=0)))\n# Obtain MRI cartesian mask\nmask = get_sample_data(\"cartesian-mri-mask\")\nkspace_loc = convert_mask_to_locations(mask.data)\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 and the associated observations\nfourier_op = FFT(samples=kspace_loc, shape=image.shape,\n                 n_coils=cartesian_ref_image.shape[0])\nkspace_obs = fourier_op.op(cartesian_ref_image)\n\n# Zero Filled reconstruction\nzero_filled = fourier_op.adj_op(kspace_obs)\nimage_rec0 = pysap.Image(data=np.sqrt(np.sum(np.abs(zero_filled)**2, axis=0)))\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 = WaveletN(\n    wavelet_name='sym8',\n    nb_scale=4,\n    n_coils=cartesian_ref_image.shape[0],\n)\nregularizer_op = GroupLASSO(weights=6e-8)\n# Setup Reconstructor\nreconstructor = CalibrationlessReconstructor(\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=300,\n)\nimage_rec = pysap.Image(data=np.sqrt(np.sum(np.abs(x_final)**2, axis=0)))\nrecon_ssim = ssim(image_rec, image)\nprint('The Reconstruction SSIM is : ' + str(recon_ssim))\n# image_rec.show()"
      ]
    }
  ],
  "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
}