Skip to content

Commit

Permalink
Updated Runs
Browse files Browse the repository at this point in the history
  • Loading branch information
miramliu committed Sep 23, 2022
1 parent 4f815f5 commit b8fc832
Show file tree
Hide file tree
Showing 8 changed files with 708 additions and 283 deletions.
21 changes: 12 additions & 9 deletions MRI/.ipynb_checkpoints/PlanarWave_FFT_Scratch-checkpoint.ipynb
Original file line number Diff line number Diff line change
Expand Up @@ -82,7 +82,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -185,22 +185,25 @@
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"fig, (ax1, ax2) = pl.subplots(1,2)\n",
"angles = np.linspace(0,np.pi,10)\n",
"ax1.set_title('Planar Wave')\n",
"ax2.set_title('FT of Planar Wave')\n",
"#fig, (ax1, ax2) = pl.subplots(1,2)\n",
"angles = np.linspace(0,np.pi/2,10)\n",
"#ax1.set_title('Planar Wave')\n",
"#ax2.set_title('FT of Planar Wave')\n",
"for tilt in angles:\n",
" fig, (ax1, ax2) = pl.subplots(1,2)\n",
" #phi being the tilt angle\n",
" img = np.sin(2*np.pi*(X*np.cos(tilt) + Y*np.sin(tilt)) / 20)\n",
" ax1.imshow(img,cmap = 'gray')\n",
" ax1.imshow(img,cmap = 'gray',extent = [0,1,0,1])\n",
" ft = np.fft.fft2(img)\n",
" kspace = np.fft.fftshift(ft) #moves zero frequency component to center\n",
" ax2.imshow(np.abs(kspace),cmap = 'gray')\n",
"\n",
" #ax2.plot(img[:,51],x)\n",
" \n",
" pl.tight_layout()\n",
" pl.pause(.5)\n",
" \n",
" "
Expand Down Expand Up @@ -247,7 +250,7 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
Expand Down
315 changes: 315 additions & 0 deletions MRI/.ipynb_checkpoints/RunMe-Original-checkpoint.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,315 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"# Load packages \n",
"from PrepFunctions import * "
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Using matplotlib backend: MacOSX\n"
]
}
],
"source": [
"%matplotlib"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# How magnetic waves are used to construct non-invasive internal anatomic and physiologic images \n",
"\n",
"Mira Liu 9/09/2022\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# How the way we hear sounds is like MRI:\n",
"\n",
"### - When plucking a string, one can hear a note. \n",
"### - That note comes from oscillation at a particular frequency. \n",
"### - We then hear the sum of the frequencies, that's a chord.\n",
"\n",
"# BUT when we hear a chord, we don't just hear that one sum: \n",
"### - We hear the individual notes (waves) that are making up that chord! \n",
"### - Our brain can deconstruct a chord and hear the individual notes of the chord \n",
"### - Being able to \"hear\" individual frequencyes from a sum is how MRI works! \n",
"\n",
"## <span style=\"font-weight:bold;color:blue\">Run the following code (click the code below and press \"shift\" and \"return\") to see the C major chord as the individual notes, and their sum at the bottom is the sound that reaches our ears:</span>"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"# Three sound waves \n",
"\n",
"x = np.linspace(0, .1, 1000) #from 0 to .1 seconds, sampling rate of .1/1000\n",
"\n",
"C_wave = np.sin(2*np.pi*261*x) #frequency of middle C 261 Hz\n",
"E_wave = np.sin(2*np.pi*329*x) #frequency of E 329 Hz\n",
"G_wave = np.sin(2*np.pi*392*x) #frequency of G 392 Hz\n",
"\n",
"ShowThreeWaves(x,C_wave,E_wave,G_wave)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Fourier Transform: \n",
"### - Our brain is performing a \"fourier transform\", where signal is decomposed into the waves of varying frequencies. \n",
"### - This \"decomposition\" is called a fourier transform or \"FT\"\n",
"### - One wave with one frequency returns a line at the frequency it is oscillating at. \n",
"### - A chord, the sum of multiple waves, will show lines at the frequencies of the notes the chord is made of!\n",
"\n",
"## <span style=\"font-weight:bold;color:blue\">Run the following code to see how as we increase the frequency of the wave (higher pitch from C to E to G) it shifts frequency in fourier space, and how a chord is made of those three individual notes!</span>"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"# Load Data\n",
"with open('data/Fig1_1dChordWaves.pickle','rb') as handle: \n",
" wave1,ft_sin1,wave2,ft_sin2,wave3,ft_sin3,wavesum,ft_sinsum = pickle.load(handle)\n",
"\n",
"# Three sound waves and their fourier transform\n",
"x = np.linspace(0, .05, 10001) #from 0 to 1, sampling rate of 1/1000ft_freq = np.linspace(-500,500,1000) #kmax = 1/2Deltax, Delta kx = 1/FOVx\n",
"ft_freq = np.linspace(-100000,100000,10001) #kmax = 1/2Deltax, Delta kx = 1/FOVx\n",
"\n",
"ShowFTOfWaves(x,ft_freq,wave1,wave2,wave3,ft_sin1,ft_sin2,ft_sin3,wavesum,ft_sinsum)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Hearing sound waves make sense, but what about in images?\n",
"### - Make those 1D sound waves two dimensional (2D), they're planar waves now\n",
"### - this means they will be two dimensional frequency space\n",
"### - that doesn't make much sense intuitively, does it... so let's take a look:\n",
"\n",
"## <span style=\"font-weight:bold;color:blue\">Run the following code to see that the 2D planar wave of a single frequency, and its 2D fourier transform! Like before, it's just a point at the wave's frequency! </span>"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [],
"source": [
"with open('data/Fig2_2dWaves.pickle','rb') as handle: \n",
" img,kspace = pickle.load(handle)\n",
" \n",
"x = np.linspace(0, 1, 100)\n",
"ft_freq = np.linspace(-50,50,100)\n",
"\n",
"ShowFTof2Dwaves(x,ft_freq,img,kspace)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# What does 2D wave mean?\n",
"\n",
"### - A 1D wave doesn't have any direction, it only has frequency. (this can be hard to understand, but if you think about it, even for a 1D wave to travel in a direction, that direction has to be a direction in 2D space!)\n",
"\n",
"### - A 2D wave has both a frequency AND a direction in which it is going!\n",
"### - This means that the direction of the frequency must be included in the Fourier Transform.\n",
"\n",
"## <span style=\"font-weight:bold;color:blue\">Run the following code to see rotating planar waves means rotation of the fourier transform!</span>"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [],
"source": [
"with open('data/Fig3_Rotating.pickle','rb') as handle: \n",
" image_rotating,kspace_rotating = pickle.load(handle)\n",
"\n",
"fig, (ax1, ax2) = pl.subplots(1,2)\n",
"for j in range(np.shape(image_rotating)[2]):\n",
" ax1.imshow(image_rotating[:,:,j],cmap = 'gray',extent = [0,1,0,1])\n",
" ax1.set_title('Rotating\\nPlanar Wave')\n",
" ax2.imshow(np.abs(kspace_rotating[:,:,j]),cmap = 'gray',extent = [-50,50,-50,50])\n",
" ax2.set_title('Corresponding FT\\n\"k-space\"')\n",
"\n",
" pl.pause(.5)\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# What about the 2D frequency though? \n",
"\n",
"## Remember the 1D chord notes? Just like that...\n",
"### - frequency is from the distance of the point from the center\n",
"### - the direction of the wave is the angle of the delta function from the center! \n",
"\n",
"## <span style=\"font-weight:bold;color:blue\">Run the following code to see a wave that is increasing in frequency while rotating. We can see that the fourier transform is showing the points are moving away from the center and rotating as well!</span>"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"with open('data/Fig4_RotatingIncreasing.pickle','rb') as handle: \n",
" image_rotateincrease,kspace_rotateincrease = pickle.load(handle)\n",
"fig, (ax1, ax2) = pl.subplots(1,2)\n",
"for j in range(np.shape(image_rotating)[2]):\n",
" ax1.imshow(image_rotateincrease[:,:,j],cmap = 'gray',extent = [0,1,0,1])\n",
" ax1.set_title('Rotating & Increasing Hz\\nPlanar Wave')\n",
" ax2.imshow(np.abs(kspace_rotateincrease[:,:,j]),cmap = 'gray',extent = [-50,50,-50,50])\n",
" ax2.set_title('Corresponding FT\\n\"k-space\"')\n",
"\n",
" pl.pause(.5)\n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# From this MRI can make images! \n",
"\n",
"### - Just as a chord is formed by summing the 1D waves of individual notes, an image is formed by summing the 2D waves of individual frequencies! \n",
"\n",
"### - MRIs can use magnetic fields to measure the frequency, direction, and phase of the 2D wavepatterns of tissue. \n",
"### - By imaging this we can get the individual contribution of each type of wave, and from this we can get images of biologic anatomy and physiology. \n",
"\n",
"## <span style=\"font-weight:bold;color:blue\">Run the following code to see an example of an MRI brain scan and how it is formed from summing 2D waves</span>\n"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [],
"source": [
"name = 'Coronal' #Sagittal, Axial\n",
"\n",
"fov_img = 256\n",
"\n",
"with open('data/'+name + '_superposition.pickle','rb') as handle1:\n",
" stack_image = pickle.load(handle1)\n",
"with open('data/'+ name +'_superposition_masks.pickle','rb') as handle2:\n",
" stack_waves = pickle.load(handle2)\n",
"with open('data/'+ name + '_superposition_kspace.pickle','rb') as handle3:\n",
" stack_kspace = pickle.load(handle3)\n",
"with open('data/'+ name + '_FullImages.pickle','rb') as handle4:\n",
" img,kspace = pickle.load(handle4)\n",
" \n",
" \n",
"fig, (ax1, ax2, ax3,ax4) = pl.subplots(1,4)\n",
"ax4.imshow(img,cmap = 'gray')\n",
"ax4.set_title('Final Image')\n",
"for i in np.arange(0,np.shape(stack_image)[2]-65):\n",
" ax1.imshow(stack_kspace[:,:,i],vmin = 0,vmax = 8000,extent = [-.5, .5, -.5, .5],cmap = 'gray')\n",
" ax1.set_title('k-space')\n",
" ax2.imshow(stack_waves[:,:,i],cmap = 'gray')\n",
" ax2.set_title('planar\\nwaves')\n",
" ax3.imshow(stack_image[:,:,i],cmap = 'gray')\n",
" ax3.set_title('Image\\nForming')\n",
" pl.pause(.1)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# See what the MRI measures, and the image we form from it: \n",
"\n",
"<span style=\"font-weight:bold;color:blue\">Run the following code to see an example of an MRI brain scan and how it is formed from summing 2D waves</span>"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"fig, ((ax1,ax2), (ax3,ax4)) = pl.subplots(2,2)\n",
"\n",
"# MRI image\n",
"ax1.imshow(np.abs(kspace),vmin = 0,vmax = 8000,extent = [-.5, .5, -.5, .5],cmap = 'gray')\n",
"ax1.set_title('k-space collected \\nwith MRI')\n",
"ax2.imshow(img,cmap = 'gray')\n",
"ax2.set_title('Image from MRI')\n",
"\n",
"# Chord\n",
"x = np.linspace(0, .05, 10001) #from 0 to 1, sampling rate of 1/1000ft_freq = np.linspace(-500,500,1000) #kmax = 1/2Deltax, Delta kx = 1/FOVx\n",
"ft_freq = np.linspace(-100000,100000,10001) #kmax = 1/2Deltax, Delta kx = 1/FOVx\n",
"\n",
"ax4.plot(x,wavesum)\n",
"ax4.set_title('C Chord\\n(What our ears receive)')\n",
"ax4.set_xlabel('Time (s)')\n",
"ax3.plot(ft_freq,abs(ft_sinsum))\n",
"ax3.set_title('FT of C Chord\\n(What our brains process)')\n",
"ax3.set_xlim([0,500])\n",
"ax3.vlines(x=261,ymin =0,ymax=5000,color = 'red',label = ' C')\n",
"ax3.vlines(x=320,ymin =0,ymax=4000,color = 'purple', label = 'E')\n",
"ax3.vlines(x=400,ymin =0,ymax=4000,color = 'green',label = 'G')\n",
"ax3.set_xlabel('Hz (1/s)')\n",
"ax3.legend()\n",
"fig.tight_layout()\n",
"pl.show()\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"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.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Loading

0 comments on commit b8fc832

Please sign in to comment.