[ismrmrd] 226/281: adding crude Python example create/recon dataset scripts
Ghislain Vaillant
ghisvail-guest at moszumanska.debian.org
Wed Jan 14 20:01:18 UTC 2015
This is an automated email from the git hooks/post-receive script.
ghisvail-guest pushed a commit to annotated tag ismrmrd0.5
in repository ismrmrd.
commit db681da995282dc8764f9da88a386fa1d654420f
Author: Joseph Naegele <joseph.naegele at gmail.com>
Date: Thu Apr 24 10:35:38 2014 -0400
adding crude Python example create/recon dataset scripts
---
examples/python/ismrmrd_create_dataset.py | 165 ++++++++++++++++++++++++++++++
examples/python/ismrmrd_recon_dataset.py | 103 +++++++++++++++++++
examples/python/ismrmrd_test.py | 22 ----
3 files changed, 268 insertions(+), 22 deletions(-)
diff --git a/examples/python/ismrmrd_create_dataset.py b/examples/python/ismrmrd_create_dataset.py
new file mode 100644
index 0000000..1a14ce2
--- /dev/null
+++ b/examples/python/ismrmrd_create_dataset.py
@@ -0,0 +1,165 @@
+# coding: utf-8
+import os
+import ismrmrd
+import ismrmrd_xsd
+import numpy as np
+import matplotlib.pyplot as plt
+
+
+filename = 'testdata.h5'
+if os.path.isfile(filename):
+ os.remove(filename)
+# Create an empty ISMRMRD dataset
+dset = ismrmrd.IsmrmrdDataset(filename, 'dataset')
+
+# Synthesize the object
+nX, nY = 256, 256
+rho = np.zeros((nX, nY))
+x0, x1 = nX / 4, 3 * nX / 4
+y0, y1 = nY / 4, 3 * nY / 4
+rho[x0:x1, y0:y1] = 1
+
+plt.imshow(rho)
+
+# Synthesize some coil sensitivities
+X, Y = np.meshgrid(np.arange(nX) / nX / 2.0, np.arange(nY) / nY / 2.0)
+C = np.zeros((nX, nY, 4), dtype=np.complex64)
+C[:,:,0] = np.exp(-((X - 0.5) ** 2 + Y ** 2) + 1j * (X - 0.5))
+C[:,:,1] = np.exp(-((X + 0.5) ** 2 + Y ** 2) - 1j * (X + 0.5))
+C[:,:,2] = np.exp(-(X ** 2 + (Y - 0.5) ** 2) + 1j * (Y - 0.5))
+C[:,:,3] = np.exp(-(X ** 2 + (Y + 0.5) ** 2) - 1j * (Y + 0.5))
+ncoils = np.size(C, 2)
+
+# Synthesize the k-space data
+nreps = 5
+noiselevel = 0.05
+K = np.zeros((nX, nY, ncoils, nreps), dtype=np.complex64)
+for rep in range(nreps):
+ for coil in range(ncoils):
+ noise = noiselevel * (np.random.randn(nX, nY) + 1j * np.random.randn(nX, nY))
+ K[:,:,coil,rep] = np.fft.fftshift(np.fft.fft2(np.fft.fftshift(C[:,:,coil] * rho + noise)))
+
+rep0 = np.sqrt(np.sum(np.abs(K) ** 2, 2))
+plt.imshow(rep0[:,:,0])
+
+# Append a block of acquisitions, one repetition at a time (faster than one acquisition at a time)
+# Re-use the same Acquisition object because why not?
+acq = ismrmrd.Acquisition()
+for rep in range(nreps):
+ for line in range(nY):
+ head = ismrmrd.AcquisitionHeader()
+ head.version = 1
+ head.number_of_samples = nX
+ head.center_sample = nX / 2
+ head.active_channels = ncoils
+ rd, pd, sd = ismrmrd.floatArray(3), ismrmrd.floatArray(3), ismrmrd.floatArray(3)
+ rd[0] = 1; pd[1] = 1; sd[2] = 1
+ head.read_dir = rd
+ head.phase_dir = pd
+ head.slice_dir = sd
+ head.scan_counter = rep * nY + line
+ head.idx.kspace_encode_step_1 = line
+ head.idx.repetition = rep
+
+ # Note: the correct API for setting Acquisition flags looks like this:
+ # acq.setFlag(ismrmrd.FlagBit(ismrmrd.ACQ_FIRST_IN_ENCODE_STEP1))
+ # but using this API would require using only ugly "acq.setXXX" methods
+ # since the call to "acq.setHead()" below overwrites the Acquisition's header
+ head.flags = 0
+ if line == 0:
+ head.flags |= 1 << ismrmrd.ACQ_LAST_IN_ENCODE_STEP1
+ head.flags |= 1 << ismrmrd.ACQ_FIRST_IN_SLICE
+ head.flags |= 1 << ismrmrd.ACQ_FIRST_IN_REPETITION
+ elif line == nY - 1:
+ head.flags |= 1 << ismrmrd.ACQ_LAST_IN_ENCODE_STEP1
+ head.flags |= 1 << ismrmrd.ACQ_LAST_IN_SLICE
+ head.flags |= 1 << ismrmrd.ACQ_LAST_IN_REPETITION
+
+ acq.setHead(head)
+ # copy the array to ensure it's a continuous block of memory
+ data = np.zeros((2 * nX, ncoils), dtype=np.float32)
+ data[::2] = np.array([c.real for c in np.array(K[:,line,:,rep])])
+ data[1::2] = np.array([c.imag for c in np.array(K[:,line,:,rep])])
+ acq.setData(data.flatten())
+ dset.appendAcquisition(acq)
+
+
+# Fill the XML header
+header = ismrmrd_xsd.ismrmrdHeader()
+
+# Experimental Conditions
+exp = ismrmrd_xsd.experimentalConditionsType()
+exp.H1resonanceFrequency_Hz = 128000000
+header.experimentalConditions = exp
+
+# Acquisition System Information
+sys = ismrmrd_xsd.acquisitionSystemInformationType()
+sys.receiverChannels = ncoils
+header.acquisitionSystemInformation = sys
+
+# Encoding
+encoding = ismrmrd_xsd.encoding()
+encoding.trajectory = ismrmrd_xsd.trajectoryType.cartesian
+
+# Encoded Space
+fov = ismrmrd_xsd.fieldOfView_mm()
+fov.x = 256
+fov.y = 256
+fov.z = 5
+
+matrix = ismrmrd_xsd.matrixSize()
+matrix.x = np.size(K, 0)
+matrix.y = np.size(K, 1)
+matrix.z = 1
+
+space = ismrmrd_xsd.encodingSpaceType()
+space.matrixSize = matrix
+space.fieldOfView_mm = fov
+
+# Set encoded and recon space (same)
+encoding.encodedSpace = space
+encoding.reconSpace = space
+
+# Encoding limits
+limits = ismrmrd_xsd.encodingLimitsType()
+
+limits0 = ismrmrd_xsd.limitType()
+limits0.minimum = 0
+limits0.center = np.size(K, 0) / 2
+limits0.maximum = np.size(K, 0) - 1
+limits.kspaceEncodingStep0 = limits0
+
+limits1 = ismrmrd_xsd.limitType()
+limits1.minimum = 0
+limits1.center = np.size(K, 1) / 2
+limits1.maximum = np.size(K, 1) - 1
+limits.kspaceEncodingStep1 = limits1
+
+limits_rep = ismrmrd_xsd.limitType()
+limits_rep.minimum = 0
+limits_rep.center = nreps / 2
+limits_rep.maximum = nreps - 1
+limits.repetition = limits_rep
+
+limits_slice = ismrmrd_xsd.limitType()
+limits_slice.minimum = 0
+limits_slice.center = 0
+limits_slice.maximum = 0
+limits.slice = limits_slice
+
+limits_rest = ismrmrd_xsd.limitType()
+limits_rest.minimum = 0
+limits_rest.center = 0
+limits_rest.maximum = 0
+limits.average = limits_rest
+limits.contrast = limits_rest
+limits.kspaceEncodingStep2 = limits_rest
+limits.phase = limits_rest
+limits.segment = limits_rest
+limits.set = limits_rest
+
+encoding.encodingLimits = limits
+header.encoding.append(encoding)
+
+dset.writeHeader(header.toxml('utf-8'))
+dset.close()
diff --git a/examples/python/ismrmrd_recon_dataset.py b/examples/python/ismrmrd_recon_dataset.py
new file mode 100644
index 0000000..3d6c47d
--- /dev/null
+++ b/examples/python/ismrmrd_recon_dataset.py
@@ -0,0 +1,103 @@
+# coding: utf-8
+
+import os
+import ismrmrd
+import ismrmrd_xsd
+import numpy as np
+from numpy.fft import ifft, fftshift
+import matplotlib.pyplot as plt
+
+# Load file
+filename = 'testdata.h5'
+if not os.path.isfile(filename):
+ print("%s is not a valid file" % filename)
+ raise SystemExit
+dset = ismrmrd.IsmrmrdDataset(filename, 'dataset')
+
+header = ismrmrd_xsd.CreateFromDocument(dset.readHeader())
+enc = header.encoding[0]
+
+# Matrix size
+eNx = enc.encodedSpace.matrixSize.x
+eNy = enc.encodedSpace.matrixSize.y
+eNz = enc.encodedSpace.matrixSize.z
+rNx = enc.reconSpace.matrixSize.x
+rNy = enc.reconSpace.matrixSize.y
+rNz = enc.reconSpace.matrixSize.z
+
+# Field of View
+eFOVx = enc.encodedSpace.fieldOfView_mm.x
+eFOVy = enc.encodedSpace.fieldOfView_mm.y
+eFOVz = enc.encodedSpace.fieldOfView_mm.z
+rFOVx = enc.reconSpace.fieldOfView_mm.x
+rFOVy = enc.reconSpace.fieldOfView_mm.y
+rFOVz = enc.reconSpace.fieldOfView_mm.z
+
+# Number of Slices, Reps, Contrasts, etc.
+nslices = enc.encodingLimits.slice.maximum + 1
+ncoils = header.acquisitionSystemInformation.receiverChannels
+nreps = enc.encodingLimits.repetition.maximum + 1
+ncontrasts = enc.encodingLimits.contrast.maximum + 1
+
+# TODO: Ignore noise scans
+for acqnum in range(dset.getNumberOfAcquisitions()):
+ acq = dset.readAcquisition(acqnum)
+ if acq.getHead().flags & ismrmrd.ACQ_IS_NOISE_MEASUREMENT:
+ print("Found noise measurement @ %d" % acqnum)
+
+all_data = np.zeros((eNx, eNy, eNz, ncoils, nslices, ncontrasts, nreps), dtype=np.complex64)
+for acqnum in range(dset.getNumberOfAcquisitions()):
+ acq = dset.readAcquisition(acqnum)
+ head = acq.getHead()
+ rep = head.idx.repetition
+ contrast = head.idx.contrast
+ slice = head.idx.slice
+ y = head.idx.kspace_encode_step_1
+ z = head.idx.kspace_encode_step_2
+ floatdata = acq.getData()
+ data = np.array([complex(elem[0], elem[1]) for elem in zip(floatdata[::2], floatdata[1::2])])
+ data = data.reshape((eNx, ncoils))
+ all_data[:, y, z, :, slice, contrast, rep] = data
+
+fig = plt.figure()
+h, w = nreps * ncontrasts, eNz * nslices
+i = 0
+for rep in range(nreps):
+ for contrast in range(ncontrasts):
+ for slice in range(nslices):
+ for z in range(eNz):
+ K = all_data[:,:,z,:,slice, contrast, rep]
+ comb = np.sqrt(np.sum(np.abs(K) ** 2, 2))
+ a = fig.add_subplot(h, w, i)
+ plt.imshow(comb)
+ i += 1
+fig.set_size_inches(16, 16)
+
+images = []
+for rep in range(nreps):
+ for contrast in range(ncontrasts):
+ for slice in range(nslices):
+ K = all_data[:,:,:,:,slice, contrast, rep]
+ K = fftshift(ifft(fftshift(K, axes=0), axis=0), axes=0)
+
+ # chop if needed
+ if eNx != rNx:
+ i0 = (eNx - rNx) / 2
+ i1 = (eNx - rNx) / 2 + rNx
+ im = K[i0:i1,:,:,:]
+ else:
+ im = K
+
+ im = fftshift(ifft(fftshift(im, axes=1), axis=1), axes=1)
+ if np.size(im, 2) > 1:
+ im = fftshift(ifft(fftshift(im, axes=2), axis=2), axes=2)
+
+ im = np.squeeze(np.sqrt(np.sum(np.abs(im) ** 2, 3)))
+ images.append(im)
+
+l = len(images)
+fig = plt.figure()
+for n, im in enumerate(images):
+ a = fig.add_subplot(1, 5, n)
+ plt.imshow(im)
+fig.set_size_inches(16, 4)
diff --git a/examples/python/ismrmrd_test.py b/examples/python/ismrmrd_test.py
deleted file mode 100644
index ede39ab..0000000
--- a/examples/python/ismrmrd_test.py
+++ /dev/null
@@ -1,22 +0,0 @@
-import ismrmrd
-import ismrmrd_xsd
-import numpy as np
-
-f = ismrmrd.IsmrmrdDataset('./testdata.h5', '/dataset')
-
-a = f.readAcquisition(10)
-#print a.getData()
-#print(a.getPosition(0), a.getPosition(1), a.getPosition(2))
-#print(a.getReadDirection(0), a.getReadDirection(1), a.getReadDirection(2))
-
-d = a.getData()
-d2 = np.zeros(d.shape)
-a.setData(d2)
-
-i = f.readImage_ushort('the_square')
-print(i.getReadDirection(0), i.getReadDirection(1), i.getReadDirection(2))
-
-xml = f.readHeader()
-hdr = ismrmrd_xsd.CreateFromDocument(xml)
-
-print hdr.experimentalConditions.H1resonanceFrequency_Hz
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/ismrmrd.git
More information about the debian-science-commits
mailing list