{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# k-Space lab in Python\n",
    "\n",
    "> **_NOTE:_**  Using generative AI to solve these problems, or to write the lab report, is _not_ allowed.\n",
    "\n",
    "This is an interactive Jupyter Notebook. The notebook consists of markdown cells (like this) and code cells. To edit a code cell, you just click in it and start typing. To execute code you can either click the Run button, or use the keyboard shortcut Shift + Enter.\n",
    "\n",
    "This lab is making use of the NumPy (<a href=\"https://numpy.org/devdocs/user/quickstart.html\">tutorial</a>) and PyPlot (<a href=\"https://matplotlib.org/3.1.1/tutorials/introductory/pyplot.html\">tutorial</a>) packages.\n",
    "\n",
    "A few tips for numpy:\n",
    "* Numpy uses zero-based indexing, e.g. ``my_nparray[0]`` is the first element.\n",
    "* Arrays can be sliced using the ``[start:stop:step]`` notation, e.g. ``my_nparray[::2]`` for every second element or ``my_nparray[0:20]`` for the 20 first elements.\n",
    "* Before modifying a variable, it is a good idea to make a copy, e.g. ``my_kspace = kspace_res512_FOV24.copy()``, otherwise you might end up with unwanted results.\n",
    "* If you want to initialize a zero matrix with ``np.zeros``, you need to specify that it's a complex datatype using ``dtype=np.complex64``\n",
    "* The command ``np.log`` may give undesired results for arrays containing zeros. You can use ``np.log(x+1)`` for those instances.\n",
    "* ``ij`` or ``0+ij`` can be used to represent complex numbers. ``np.pi`` can be used for pi.\n",
    "* In numpy convention, the first dimension (axis=0) corresponds to rows and the second dimension (axis=1) corresponds to columns.\n",
    "\n",
    "A few tips for matplotlib:\n",
    "* Use ``plt.imshow`` to display images, and use the option ``cmap='gray'`` to display them in black and white.\n",
    "* When the question calls for multiple plots, you can either use ``plt.figure()`` to create multiple plots or use ``fig, axes = plt.subplots(2,2)`` to create a 2 by 2 grid.\n",
    "* If your figures appear too small, you can use ``figsize=(20,20)`` to make the plotting area larger, e.g. ``plt.figure(figsize=(20,20))`` or ``fig, axes = plt.subplots(2,2, figsize=(20,20))``\n",
    "\n",
    "Finally, if you need to reset your workspace, run the first cell again!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lab report\n",
    "Submit a written lab report no later than 3 weeks after the lab. The report should answer the questions posed for each task, and should include images and necessary code snippets. Each student should write and submit an individual lab report, even if you are working together during the session. As noted above, using generative AI (ChatGPT, Claude, Gemini, Co-Pilot, Llama etc.) is not allowed for writing the lab report or for solving the programming tasks. If you use a grammar and/or spell checker that has generative AI features (such as Grammarly or Writefull), please declare this in your lab report.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"RUN THIS CELL TO SETUP OR RESET THE ENVIRONMENT \n",
    "\"\"\"\n",
    "!wget -nc https://raw.githubusercontent.com/fyrdahl/kspace-lab/master/data.h5\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import h5py\n",
    "\n",
    "from scipy.signal import convolve2d\n",
    "\n",
    "with h5py.File('data.h5', 'r') as F:\n",
    "    kspace_res512_FOV24 = np.array(F['kspace_res512_FOV24'])\n",
    "    kspace_res512_FOV48 = np.array(F['kspace_res512_FOV48'])\n",
    "\n",
    "fftc = lambda f: np.fft.fftshift(np.fft.fft(np.fft.ifftshift(f, axes=0), axis=0), axes=0)\n",
    "ifftc = lambda F: np.fft.fftshift(np.fft.ifft(np.fft.ifftshift(F, axes=0), axis=0), axes=0)\n",
    "fft2c = lambda f: np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(f)))\n",
    "ifft2c = lambda F: np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(F)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1\n",
    "\n",
    "Use ``kspace_res512_FOV24``\n",
    "\n",
    "1. Show the magnitude, phase, real, and imaginary part of the k-space. To better visualize the k-space, you can use ``np.log``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2,2)\n",
    "# Plot magnitude here\n",
    "axes[0][0].imshow(np.abs(np.log(kspace_res512_FOV24)), cmap='gray')\n",
    "axes[0][0].set_title('Magnitude')\n",
    "axes[0][0].set_axis_off()\n",
    "\n",
    "# Plot phase here\n",
    "\n",
    "# Plot real part here\n",
    "\n",
    "# Plot imaginary part here\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2.  Take the inverse Fourier transform of the k-space data along one dimension (ifftc).\n",
    "    * Display the magnitude and phase of the result.\n",
    "    * Describe what happened, which domain is the data in – spatial or frequency?\n",
    "    * What are the unit on the x and y axis?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. Do an inverse FT of the k-space data in both the frequency and phase encoding direction.\n",
    "    * Display the magnitude and phase of the result.\n",
    "    * Which domain is the data in – spatial or frequency?\n",
    "    * Which dimension is most likely to be the phase encoding direction?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2\n",
    "\n",
    "Use ``kspace_res512_FOV48``\n",
    "\n",
    "1. Decrease the FOV in the spatial domain to 24x24 cm, and 16x16 cm _by removing parts of k-space_ \n",
    "    * Show the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3\n",
    "\n",
    "Use ``kspace_res512_FOV24``\n",
    "\n",
    "1. Decrease the spatial resolution by replacing k-space data with zeros before performing the ifft2c.\n",
    "    * The new images should have a 512x256, 256x256, 64x64, 512x64 and 64x512 resolution.\n",
    "    * Display the magnitude of k-space and the magnitude image in the spatial domain.\n",
    "    * Remember that you should use ``np.log(x+1)`` when taking the \"log\" a matrix containing zeros.\n",
    "2. What is the difference between the 512x64 and the 64x512 image?\n",
    "    * Describe the artifact that is appearing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4\n",
    "\n",
    "Use ``kspace_res512_FOV24``\n",
    "\n",
    "1. Set the k-space element ``[kx,ky]`` to 200000 for _each of_ these elements: ``[100,100]``, ``[200,200]``, ``[250,250]``, ``[255,255]``, ``[257,257]``.\n",
    "    * Display one image for each k-space adjustment.\n",
    "    * Describe the artifact.\n",
    "2. What is the correlation between location of the element you changed and the apperance of the artifact.\n",
    "3. What kind of artifact do you think this could simulate?\n",
    "\n",
    "> **_NOTE:_** Numpy arrays are mutable. To create an image for _each_ k-space location, you need to use ``.copy()`` to create copies of the original array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 5\n",
    "\n",
    "Use ``kspace_res512_FOV24``\n",
    "\n",
    "Change the k-space phase without altering the magnitude. Remember that complex\n",
    "numbers can be written as ``z = mag * np.exp(1j*phase)``\n",
    "\n",
    "1. _Increase_ the phase by 2 radians in _each_ of the following k-space rows:\n",
    "    * 1–20\n",
    "    * 201–220\n",
    "    * 237–256\n",
    "    * 247–266\n",
    "2. Display the magnitude image in the spatial domain for _each_ k-space adjustment.\n",
    "    * How does the appearance of the image correlate with the rows you adjusted?\n",
    "3. What kind of artifact do you think this could simulate?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 6\n",
    "\n",
    "Use ``kspace_res512_FOV24``\n",
    "\n",
    "1. Filter the spatial image with a boxcar filter using two kernel sizes: 7x7 and 4x4.\n",
    "3. Do an inverse Fourier transform and display the magnitude of k-space for each kernel size. Describe the results.\n",
    "\n",
    "> **_NOTE:_** To define a boxcar filter, you can use ``filter = np.ones((n,n))``, then use ``convolve2d(image,filter)`` to perform a 2D convolution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 7\n",
    "In your written report, reflect on the results of this exerice. Was there something in particular that you found particularly interesting about the connection between image space and k-space. Did any of the tasks help clarify a concept that was discussed in previous lectures or exercises?"
   ]
  }
 ],
 "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.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}