Day3 homework

day3/slice_time_image.py

@@ -21,10 +21,12 @@
 """
 
 # Add any extra imports here
+import os
 import sys
-
+import numpy as np
+import matplotlib.pyplot as plt
 import nibabel as nib
-
+from scipy.interpolate import InterpolatedUnivariateSpline as IUS
 
 def slice_time_image(img, slice_times, TR):
     """Run slice-timing correction on nibabel image 'img'.
@@ -47,8 +49,32 @@ def slice_time_image(img, slice_times, TR):
         A new copy of the input image with slice-time interpolation applied
     """
     # Get the image data as an array;
+    data = img.get_data()
     # Make a new empty array "interp_data" to hold the interpolated data;
+    interp_data = np.zeros(data.shape)
+
     # Do the interpolation;
+    time_offset = slice_times/slice_times.shape * TR
+    print "Time Offsets: %s" % time_offset
+    slice_0_times = TR * np.arange(data.shape[-1])
+
+    # For each z coordinate (slice index):
+    for z in range(data.shape[2]):
+        print "Starting Slice %s of %s" % (z+1, data.shape[2])
+        slice_z_times = slice_0_times + time_offset[z]
+        # For each x coordinate:
+        for x in range(data.shape[0]):
+            # For each y coordinate:
+            for y in range(data.shape[1]):
+                # extract the time series at this x, y, z coordinate;
+                current_timecourse = data[x,y,z,:]
+                # make a linear interpolator object with the `slice_z_times` and the extracted time series;
+                current_interp = IUS(slice_z_times, current_timecourse, k=1)
+                # resample this interpolator at the slice 0 times;
+                current_timecourse_estimate = current_interp(slice_0_times)
+                # put this new resampled time series into the new data at the same position
+                interp_data[x,y,z,:] = current_timecourse_estimate
+
     # This is how to make a new image with the interpolated data
     new_img = nib.Nifti1Image(interp_data, img.affine, img.header)
     return new_img
@@ -71,7 +97,13 @@ def slice_time_file(fname, slice_times, TR):
     TR: double
         Repetition time in seconds.
     """
-    # Hint: use os.path.split and os.path.join to make the new filename
+    img = nib.load(fname)
+    interp_img = slice_time_image(img, slice_times, TR)
+
+    path, tail = os.path.split(fname)
+    new_tail = 'Interp_' + tail  # I decided to use a different new filename
+    # I thought it would provide more information than just adding 'a'
+    new_fname = os.path.join(path,new_tail)
     nib.save(interp_img, new_fname)
 
 
@@ -82,12 +114,32 @@ def main():
     fname = sys.argv[1]
     # Assume the TR
     TR = 2.0
+
     # Assume the slices were acquired even slices first, inferior to
     # superior, then odd slices, inferior to superior, where the most inferior
     # slice is index 0 and 0 is an even number.  What are the slice acquisition
     # times in seconds, where the first value is the acquisition time of slice
     # 0, the second is acquisition time of slice 1, etc?
-    slice_times = ?
+
+    # have to load in image to know the number of slices there are
+    img = nib.load(fname)
+    num_slices = img.shape[2]
+    print "Number of slices: %s" % num_slices
+
+    # use ceil in the next two lines in case of an odd number of slices
+    first_half = np.arange(0, np.ceil(num_slices/2.))
+    second_half = np.arange(np.ceil(num_slices/2.), num_slices)
+    slice_times = np.array([])
+
+    for i in first_half:
+        slice_times = np.append(slice_times,first_half[i])
+        if i < len(second_half):
+            # if there are an odd number of slices
+            # len(second_half)<len(first_half)
+            # so we have to make sure that there are more slices to append
+            slice_times = np.append(slice_times,second_half[i])
+
+    print "Slice Ordering: %s" % slice_times
     slice_time_file(fname, slice_times, TR)