-
Notifications
You must be signed in to change notification settings - Fork 2
/
Copy pathBeattyPracticalSession2013.py
246 lines (189 loc) · 9.59 KB
/
BeattyPracticalSession2013.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
# -*- coding: utf-8 -*-
# ISMRM Sunrise Practical Session, Parallel Imaging Part I
#
# This document contains the first set of practical exercises for the
# ISMRM course on parallel imaging.
## Exercise Data
#
# Exercise ground truth information is found in im1.mat (brain image), smaps_phantom.mat
# (sensitivity maps) and noise_covariances.mat (noise covariance matrices).
#
# We start by clearing the workspace, loading this ground truth information
# and setting a few parameters from this.
import numpy as np
import scipy.ndimage as spImage
import matplotlib as mpl
import Display
import IsmrmSunrise
# <codecell>
print("Loading Ground Truth Information")
im1 = np.load('im1.npy')
smaps = np.load('smaps.npy')
Rn_normal_8 = np.load('Rn_normal_8.npy')
Rn_broken_8 = np.load('Rn_broken_8.npy')
Display.ShowImage2D([im1], windowTitle='Ground Truth Image')
Display.ShowImage2D(smaps, windowTitle='Coil Sensitivities', titles='Sensitivity', maxNumInRow=4)
Display.BlockOnOpenWindow()
numChannels = smaps.shape[2]
Rn = Rn_normal_8
nx = im1.shape[0]
ny = im1.shape[1]
# <codecell>
# Data Simulation
#
# Let's use our ground truth data to simulate some data
#
print("Simulating unaccelerated data")
channelIm = smaps * im1[:,:,np.newaxis]
Display.ShowImage2D(channelIm, windowTitle= 'Image multiplied by coil sensitivities', titles='Channel', maxNumInRow=4)
Display.BlockOnOpenWindow()
noiseScale = 0.1
noise = noiseScale * np.max(im1) * IsmrmSunrise.GenerateCorrelatedNoise(im1.shape, Rn)
data = IsmrmSunrise.TransformImageToKspace(channelIm, [0, 1]) + noise
# <codecell>
# Prewhiten
#
# Let's start by prewhitening so that we don't have to contend with the
# noise covariance matrix.
#
print("Prewhiten data and coil sensitivities")
dmtx = IsmrmSunrise.ComputeNoiseDecorrelationMatrixFromCovarianceMatrix(Rn)
csmTrue = IsmrmSunrise.ApplyNoiseDecorrelationMatrix(smaps, dmtx)
data = IsmrmSunrise.ApplyNoiseDecorrelationMatrix(data, dmtx)
Rn = np.eye(numChannels);
# <codecell>
## Channel Combination
#
# Let's use our ground truth data to perform an SNR-optimal channel
# combination
print("Perform SNR-optimal channel combination")
imFull = IsmrmSunrise.TransformKspaceToImage(data, [0, 1])
ccmRoemerOptimal = IsmrmSunrise.ComputeChannelCombinationMaps(csmTrue, np.eye(numChannels))
imRoemerOptimal = np.abs(np.sum(imFull * ccmRoemerOptimal, 2))
# How does this compare to our ground-truth image?
Display.ShowImage2D([im1, imRoemerOptimal, imRoemerOptimal - im1], windowTitle='Comparison of Ground Truth image to SNR-optimal channel combination', titles=["Ground Truth", "Roemer Optimal", "Difference"])
##
# The noise level varies slightly from pixel-to-pixel. Let's make a map of
# this
noiseMap = IsmrmSunrise.ComputeNoiseAmplification(ccmRoemerOptimal, np.eye(numChannels))
Display.ShowImage2D(noiseMap, windowTitle="Noise Map", colormap=mpl.cm.jet)
##
# How does this compare to a root-sum-of-squares (SOS) reconstruction?
imSoS = IsmrmSunrise.ComputeRootSumOfSquaresChannelCombination(imFull)
Display.ShowImage2D([imRoemerOptimal, imSoS, imRoemerOptimal-imSoS], windowTitle='SoS compared to SNR-optimal channel combination', titles=["Roemer Optimal", "SoS", "Difference"])
Display.BlockOnOpenWindow()
# <codecell>
##
# Looks like we have a scaling and/or shading issue that makes it hard to
# do this comparison. Let's normalize to SOS shading
print("Compare after normalizing to SoS shading")
ccmRoemerOptimal, shadingCorrectionIm = IsmrmSunrise.NormalizeShadingToSoS(ccmRoemerOptimal)
imRoemerOptimal = np.abs(np.sum(imFull * ccmRoemerOptimal, 2));
imTrue = shadingCorrectionIm * im1;
Display.ShowImage2D([imTrue, imRoemerOptimal, imSoS], windowTitle="After correcting for shading differences", titles=["Ground Truth", "Roemer Optimal", "SoS"])
Display.ShowImage2D([imRoemerOptimal - imTrue, imSoS-imTrue], windowTitle="Difference Images after correcting for shading differences", titles=["Roemer Optimal", "SoS"])
Display.BlockOnOpenWindow()
# <codecell>
##
# SOS channel combination aligns the phase of the background noise before
# summing, leading to a DC ofset in these pixels. Let's try channel
# combination based on coil sensitivity estimates from some low resolution
# calibration data.
#
print("Try low res channel combination")
calx = 32
caly = 32
calShape = np.array([calx, caly])
calNoiseScale = 0.1;
calData = IsmrmSunrise.TransformImageToKspace(channelIm, dim=[0,1], kShape = calShape)
noise = calNoiseScale * np.max(im1) * IsmrmSunrise.GenerateCorrelatedNoise(calShape, Rn)
calData = calData + noise
calData = IsmrmSunrise.ApplyNoiseDecorrelationMatrix(calData, dmtx);
f = (np.mat(np.hamming(calShape[0])).T * np.mat(np.hamming(calShape[1]))).A
filteredCalData = calData * f[:,:,np.newaxis]
calIm = IsmrmSunrise.TransformKspaceToImage(filteredCalData, [0, 1], im1.shape);
# Use a circular region of support for pixels
pixelMask = ((np.sum(np.abs(smaps),2))>0).astype(float)
pixelMask = pixelMask[:,:, np.newaxis]
calIm = calIm * pixelMask
csmWalsh = IsmrmSunrise.EstimateCsmWalsh(calIm)
csmMckenzie = IsmrmSunrise.EstimateCsmMckenzie(calIm)
csmTrue = IsmrmSunrise.NormalizeShadingToSoS(csmTrue)[0]
csmWalsh = IsmrmSunrise.NormalizeShadingToSoS(csmWalsh)[0]
csmMckenzie = IsmrmSunrise.NormalizeShadingToSoS(csmMckenzie)[0]
ccmMckenzie = IsmrmSunrise.ComputeChannelCombinationMaps(csmMckenzie)
ccmWalsh = IsmrmSunrise.ComputeChannelCombinationMaps(csmWalsh)
imMckenzie = np.sum(imFull * ccmMckenzie, 2)
imWalsh = np.sum(imFull * ccmWalsh, 2)
Display.ShowImage2D([imTrue, imRoemerOptimal, imSoS, imWalsh, imMckenzie],
windowTitle='Comparison of channel combination methods', titles=["source image", "true csm", "SoS", "Walsh csm", "McKenzie csm"])
Display.ShowImage2D([imRoemerOptimal-imTrue, imSoS-imTrue, np.abs(imWalsh)-np.abs(imTrue), np.abs(imMckenzie)-np.abs(imTrue)],
windowTitle='Difference images for comparison of channel combination methods', titles=["true csm", "SoS", "Walsh csm", "McKenzie csm"])
#Display.ShowImage2D(csmWalsh, windowTitle="csmWalsh")
#Display.ShowImage2D(csmTrue, windowTitle="csmTrue")
Display.BlockOnOpenWindow()
# <codecell>
## Create Accelerated Data
#
print("Create Accelerated Data")
accFactor = 4;
samplingPattern = IsmrmSunrise.GenerateAcceleratedSamplingPattern(im1.shape, accFactor, 0)
samplingMask = np.logical_or(samplingPattern == 1,samplingPattern == 3).astype(float)
dataAccel = data * samplingMask[:,:,np.newaxis]
imAlias = IsmrmSunrise.TransformKspaceToImage(dataAccel, [0,1])
Display.ShowImage2D(imAlias, maxNumInRow=4, windowTitle='Aliased Images', titles='channel')
## Create Data Driven & Model Driven Joint Encoding Relations
#
print("Create Joint Encoding Relations")
kernelShape = [5, 7]
jerLookupDd = IsmrmSunrise.ComputeJerDataDriven(calData, kernelShape)
unfilteredCalIm = IsmrmSunrise.TransformKspaceToImage(calData, [0, 1], imShape = 2 * calData.shape)
jerLookupMd = IsmrmSunrise.ComputeJerModelDriven(unfilteredCalIm, kernelShape)
Display.BlockOnOpenWindow()
# <codecell>
# Use Various Calibration Approaches to Create Unmixing Images
#
print("Create Unmixing Images")
numRecons = 4;
titles = ['SENSE true csm', 'SENSE estimated csm', 'PARS', 'GRAPPA']
unmix = np.zeros([nx, ny, numChannels, numRecons], dtype=complex)
unmix[:,:,:,0] = IsmrmSunrise.ComputeSenseUnmixing(accFactor, csmTrue, np.eye(numChannels), 0) * accFactor
unmix[:,:,:,1] = IsmrmSunrise.ComputeSenseUnmixing(accFactor, csmMckenzie, np.eye(numChannels), 0) * accFactor
print("Done ComputeSenseUnmixing")
unmix[:,:,:,2] = IsmrmSunrise.ComputeJerUnmixing(jerLookupMd, accFactor, ccmMckenzie, 0.0, False)
unmix[:,:,:,3] = IsmrmSunrise.ComputeJerUnmixing(jerLookupDd, accFactor, ccmMckenzie, 0.0, False)
print("Done ComputeJerUnmixing")
#
# SENSE unmixing
Display.ShowImage2D(unmix[:,:,:,0], maxNumInRow=4, titles='Channel', windowTitle="Unmixing coefficients from true sensitivities")
# SENSE unmixing with estimated coil sensitivity maps
Display.ShowImage2D(unmix[:,:,:,1], maxNumInRow=4, titles='Channel', windowTitle="Unmixing coefficients from estimated sensitivities")
# PARS unmixing
Display.ShowImage2D(unmix[:,:,:,2], maxNumInRow=4, titles='Channel', windowTitle="Unmixing coefficients from PARS")
# GRAPPA unmixing
Display.ShowImage2D(unmix[:,:,:,3], maxNumInRow=4, titles='Channel', windowTitle="Unmixing coefficients from GRAPPA")
Display.BlockOnOpenWindow()
# Apply and Analyze Reconstructions
print("Apply and Analyze Reconstructions")
matrixShape = [nx, ny, numRecons]
aem = np.zeros(matrixShape, dtype=complex)
gmap = np.zeros(matrixShape, dtype=complex)
im_hat = np.zeros(matrixShape, dtype=complex)
im_diff = np.zeros(matrixShape, dtype=complex)
signalMask = spImage.binary_closing(im1>100.0, structure = np.ones((5,5))).astype(int)
for reconIndex in range(numRecons):
aem[:,:,reconIndex] = IsmrmSunrise.ComputeAliasingEnergyMap(signalMask, csmTrue, unmix[:,:,:,reconIndex], accFactor)
gmap[:,:,reconIndex] = IsmrmSunrise.ComputeGmap(unmix[:,:,:,reconIndex], ccmRoemerOptimal, accFactor, Rn)
im_hat[:,:,reconIndex] = np.sum(imAlias * unmix[:,:,:,reconIndex], 2)
im_diff[:,:,reconIndex] = np.abs(im_hat[:,:,reconIndex]) - np.abs(imTrue)
# aliasing energy maps
Display.ShowImage2D(aem, titles=titles, windowTitle="Aliasing Energy Maps, R=4", colormap=mpl.cm.jet)
#ismrm_imshow(aem, [0 0.1], [], titles); colormap(jet);
# g-factor (relative noise amplification) maps
Display.ShowImage2D(gmap, titles=titles, windowTitle="G-factor maps, R=4", colormap=mpl.cm.jet)
#ismrm_imshow(gmap, [0 6], [], titles); colormap(jet);
# reconstructed images
Display.ShowImage2D(im_hat, windowTitle='Reconstructed Images, R=4', titles=titles)
# difference images
Display.ShowImage2D(im_diff, windowTitle='Difference Images, R=4', titles=titles)
Display.BlockOnOpenWindow()