partial_fourier.ipynb 7.54 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Imports"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import h5py\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
25
26
27
28
29
30
31
32
33
34
    "# Load data\n",
    "\n",
    "Load complex image data from MATLAB mat-file (v7.3 or later), which is actually an HDF5 format\n",
    "\n",
    "Complex data is loaded as a (real, imag) tuple, so it neds to be explicitly converted to complex double\n",
    "\n",
    "In this section:\n",
    "- using h5py module\n",
    "- np.transpose\n",
    "- 1j as imaginary constant\n"
35
36
37
38
39
40
41
42
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
43
    "# get hdf5 object for the mat-file\n",
44
    "h = h5py.File('data.mat','r')\n",
45
46
47
48
49
50
51
    "\n",
    "# get img variable from the mat-file\n",
    "dat = h.get('img')\n",
    "\n",
    "# turn array of (real, imag) tuples into an array of complex doubles\n",
    "# transpose to keep data in same orientation as MATLAB\n",
    "img = np.transpose(dat['real'] + 1j*dat['imag'])"
52
53
54
55
56
57
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
58
59
60
61
62
63
64
65
66
    "# 6/8 Partial Fourier sampling\n",
    "\n",
    "Fourier transform the image to get k-space data, and add complex Gaussian noise\n",
    "\n",
    "To simulate 6/8 Partial Fourier sampling, zero out the bottom 1/4 of k-space\n",
    "\n",
    "In this section:\n",
    "- np.random.randn\n",
    "- np.fft\n",
67
68
    "- 0-based indexing\n",
    "- image plotting"
69
70
71
72
73
74
75
76
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
77
    "# generate normally-distributed complex noise\n",
78
    "n = np.random.randn(96,96) + 1j*np.random.randn(96,96)\n",
79
80
81
82
83
    "\n",
    "# Fourier transform the image and add noise\n",
    "y = np.fft.fftshift(np.fft.fft2(img), axes=0) + n\n",
    "\n",
    "# set bottom 24/96 lines to 0\n",
84
85
86
87
88
    "y[72:,:] = 0\n",
    "\n",
    "# show sampling\n",
    "_, ax = plt.subplots()\n",
    "ax.imshow(np.log(np.abs(np.fft.fftshift(y, axes=1))))"
89
90
91
92
93
94
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
95
96
97
98
99
100
101
102
103
    "# Estimate phase\n",
    "\n",
    "Filter the k-space data and extract a low-resolution phase estimate\n",
    "\n",
    "Filtering can help reduce ringing in the phase image\n",
    "\n",
    "In this section:\n",
    "- np.pad\n",
    "- np.hanning\n",
104
105
    "- reshaping 1D array to 2D array using np.newaxis (or None)\n",
    "- subplots with titles"
106
107
108
109
110
111
112
113
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
114
115
116
117
118
119
120
121
122
123
124
125
    "# create zero-padded hanning filter for ky-filtering\n",
    "filt = np.pad(np.hanning(48),24,'constant')\n",
    "\n",
    "# reshape 1D array into 2D array\n",
    "filt = filt[:,np.newaxis]\n",
    "# or\n",
    "# filt = filt[:,None]\n",
    "\n",
    "# generate low-res image with inverse Fourier transform\n",
    "low = np.fft.ifft2(np.fft.ifftshift(y*filt, axes=0))\n",
    "\n",
    "# get phase image\n",
126
127
128
129
130
131
132
133
    "phs = np.exp(1j*np.angle(low))\n",
    "\n",
    "# show phase estimate alongside true phase\n",
    "_, ax = plt.subplots(1,2)\n",
    "ax[0].imshow(np.angle(img))\n",
    "ax[0].set_title('True image phase')\n",
    "ax[1].imshow(np.angle(phs))\n",
    "ax[1].set_title('Estimated phase')"
134
135
136
137
138
139
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
    "# POCS reconstruction\n",
    "\n",
    "Perform the projection-onto-convex-sets (POCS) partial Fourier reconstruction method.\n",
    "\n",
    "POCS is an iterative scheme estimates the reconstructed image as any element in the intersection of the following two (convex) sets:\n",
    "1. Set of images consistent with the measured data\n",
    "2. Set of images that are non-negative real\n",
    "\n",
    "This requires prior knowledge of the image phase (hence the estimate above), and it works because although we have less than a full k-space of measurements, we're now only estimating half the number of free parameters (real values only, instead of real + imag), and we're no longer under-determined. Equivalently, consider the fact that real-valued images have conjugate symmetric k-spaces, so we only require half of k-space to reconstruct our image.\n",
    "\n",
    "In this section:\n",
    "- np.zeros\n",
    "- range() builtin\n",
    "- point-wise multiplication (*)\n",
    "- np.fft operations default to last axis, not first\n",
    "- np.maximum vs np.max\n"
156
157
158
159
160
161
162
163
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
164
    "# initialise image estimate to be zeros\n",
165
    "est = np.zeros((96,96))\n",
166
167
    "\n",
    "# set the number of iterations \n",
168
    "iters = 10\n",
169
170
    "\n",
    "# each iteration cycles between projections\n",
171
    "for i in range(iters):\n",
172
173
174
175
176
177
178
179
    "# projection onto data-consistent set:\n",
    "    # use a-priori phase to get complex image\n",
    "    est = est*phs\n",
    "    \n",
    "    # Fourier transform to get k-space\n",
    "    est = np.fft.fftshift(np.fft.fft2(est), axes=0)\n",
    "    \n",
    "    # replace data with measured lines\n",
180
    "    est[:72,:] = y[:72,:]\n",
181
182
183
184
185
186
187
188
189
190
191
192
193
    "    \n",
    "    # inverse Fourier transform to get back to image space\n",
    "    est = np.fft.ifft2(np.fft.ifftshift(est, axes=0))\n",
    "\n",
    "# projection onto non-negative reals:\n",
    "    # remove a-priori phase\n",
    "    est = est*np.conj(phs)\n",
    "    \n",
    "    # get real part\n",
    "    est = np.real(est)\n",
    "\n",
    "    # ensure output is non-negative\n",
    "    est = np.maximum(est, 0)"
194
195
196
197
198
199
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
200
201
202
203
204
205
    "# Display error and plot reconstruction\n",
    "\n",
    "The POCS reconstruction is compared to a zero-filled reconstruction (i.e., where the missing data is zeroed prior to inverse Fourier Transform)\n",
    "\n",
    "In this section:\n",
    "- print formatted strings to standard output\n",
206
    "- 2D subplots with min/max scales, figure size\n",
207
    "- np.sum sums over all elements by default"
208
209
210
211
212
213
214
215
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
216
217
218
219
220
221
222
223
224
225
226
227
    "# compute zero-filled recon\n",
    "zf = np.fft.ifft2(np.fft.ifftshift(y, axes=0))\n",
    "\n",
    "# compute rmse for zero-filled and POCS recon\n",
    "err_zf = np.sqrt(np.sum(np.abs(zf - img)**2))\n",
    "err_pocs = np.sqrt(np.sum(np.abs(est*phs - img)**2))\n",
    "\n",
    "# print errors\n",
    "print(f'RMSE for zero-filled recon: {err_zf}')\n",
    "print(f'RMSE for POCS recon: {err_pocs}')\n",
    "\n",
    "# plot both recons side-by-side\n",
228
    "_, ax = plt.subplots(2,2,figsize=(16,16))\n",
229
230
    "\n",
    "# plot zero-filled\n",
231
232
233
    "ax[0,0].imshow(np.abs(zf), vmin=0, vmax=1)\n",
    "ax[0,0].set_title('Zero-filled')\n",
    "ax[1,0].plot(np.abs(zf[:,47]))\n",
234
235
    "\n",
    "# plot POCS\n",
236
237
238
    "ax[0,1].imshow(est, vmin=0, vmax=1)\n",
    "ax[0,1].set_title('POCS recon')\n",
    "ax[1,1].plot(np.abs(est[:,47]))"
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
   ]
  }
 ],
 "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",
258
   "version": "3.7.3"
259
260
261
262
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
263
}