couple bugfixes + bicubic up/downsampling

niitools.py

@@ -15,12 +15,14 @@
         'ssd': niitools.volumes.ssd,
         'trim': niitools.volumes.trim_bounding_box,
         'pad': niitools.volumes.pad,
+        'binarize': niitools.util.binarize,
         'warp_ssd': niitools.volumes.warp_ssd,
         'gaussian_blur': niitools.volumes.gaussian_blur,
         'neutralize_affine': niitools.headers.neutralize_affine,
         'copy_header': niitools.headers.copy_header,
         'convert_type': niitools.headers.convert_type,
-        'upsample': niitools.resampling.upsample_axis,
+        'upsample': niitools.resampling.upsample,
+        'downsample': niitools.resampling.downsample,
         'masked_threshold': niitools.volumes.masked_threshold,
         'masked_threshold_count': niitools.volumes.masked_threshold_count,
         'scale_intensity': niitools.advanced_tools.scale_intensity,

niitools/advanced_tools.py

@@ -28,7 +28,9 @@ def scale_intensity(input, mask, new_intensity, output):
     mask_img = mask_nii.get_data()
     assert np.sum(np.abs(in_nii.get_affine() - mask_nii.get_affine())) < 1e-2
 
-    in_img = in_nii.get_data()
+    # reading in image and converting to float; this ensures that no overflow happens during processing
+    in_img = in_nii.get_data().astype(np.float32)
+
     out_img = in_img.copy()
     if (in_img.mean() < 0):
         print("%0.1f%% of voxels were below 0" % ())
@@ -43,8 +45,13 @@ def scale_intensity(input, mask, new_intensity, output):
 
     out_img[mask_img > 0] = voxels * (new_intensity / mode)
 
-    out_nii = nib.Nifti1Image(out_img, header=in_nii.get_header(),
-            affine=in_nii.get_affine())
+    # make data type of image to float 
+    h = in_nii.get_header()
+    out_dtype = np.float32
+    h['datatype'] = 16 # corresponds to float32
+    h['bitpix'] = 32 # corresponds to float32
+
+    out_nii = nib.Nifti1Image(out_img.astype(out_dtype), header=h, affine=in_nii.get_affine())
 
     out_nii.to_filename(output)
 

niitools/resampling.py

@@ -1,3 +1,4 @@
+#!/data/porpoise/Imaging/Packages/virtualenvs/WMH/bin/python
 from __future__ import print_function
 from __future__ import division
 import sys
@@ -7,6 +8,16 @@
 from numpy import newaxis as nax
 import nibabel as nib
 
+import scipy.ndimage as sn
+import scipy.ndimage.interpolation as sni
+import pickle
+import gc
+
+###
+# ToDo: 
+#  - Fix affine transformation of individual axis resampling (see e.g. upsample())
+###
+
 def downsample_axis(infile, outfile, axis, new_pixdim, method='linear'):
     """
     Downsamples a volume along a specified axis.
@@ -130,3 +141,204 @@ def upsample_axis(infile, outfile, outmask, axis, pixdim_ratio, method='linear')
     out_mask = nib.Nifti1Image(mask, header=hdr.copy(), affine=aff)
     out_mask.update_header()
     out_mask.to_filename(outmask)
+
+
+## for all axes at the same time
+
+def upsample(niiFileName, upsampledFile, zoom_values_file='upsampling_log.pickle', isotrop_res=True, upsample_factor=None, polynomial='3'):
+    """
+    Upsample a nifti and save the upsampled image.
+    The upsampling procedure has been implemented in a way that it is easily revertable.
+
+    Example arguments:
+    niiFileName = '10529_t1.nii.gz'
+    upsampled_file = '10529_t1_upsampled.nii.gz'
+    """
+
+    # load the nifti
+    nii = nib.load(niiFileName)
+    header = nii.get_header()
+    affine = nii.get_affine()
+
+    # make data type of image to float 
+    out_dtype = np.float32
+    header['datatype'] = 16 # corresponds to float32
+    header['bitpix'] = 32 # corresponds to float32
+
+    # in case nothing should be done
+    isotrop_res = bool(int(isotrop_res))
+    if ((not isotrop_res) and (upsample_factor is None)):
+        print('Uspampling not requested. Skipping...')
+        nii.to_filename(upsampledFile)
+        return nii
+
+    # convert input to number
+    if isotrop_res:
+        isotrop_res = float(np.min(header.get_zooms()[0:3]))
+        all_upsampling = [float(zoom)/isotrop_res for zoom in header.get_zooms()[0:3]]
+        for idx, zoom in enumerate(all_upsampling):
+            if zoom<1:
+                all_upsampling[idx]= 1.
+            else:
+                all_upsampling[idx]= np.round(zoom)
+    else:
+        upsample_factor = float(upsample_factor)
+        all_upsampling = [upsample_factor for zoom in header.get_zooms()[0:3]]
+
+
+    polynomial = int(polynomial)
+
+    old_volume = np.squeeze(nii.get_data().astype(float))
+    # get new volume shape and update upsampling based on rounded numbers
+    old_shape = old_volume.shape
+    print(old_shape)
+    new_shape = tuple([np.round(old_shape[ii]*usampling).astype(int) for ii, usampling in enumerate(all_upsampling)])
+    print(new_shape)
+    # all_upsampling = [float(new_shape[ii])/float(old_shape[ii]) for ii in np.arange(len(old_shape))]
+    # print('Upsampling with scales: ' + str(all_upsampling))
+
+    # # get indices
+    # a,b,c = np.indices(new_shape)
+    # a = a.ravel()
+    # b = b.ravel()
+    # c = c.ravel()
+
+    # upsample image
+    print('Upsampling volume...')
+    # u_values = sni.map_coordinates(old_volume, [a/all_upsampling[0], b/all_upsampling[1], c/all_upsampling[2]])
+    # del nii
+    # gc.collect()
+    # vol = np.zeros(new_shape, dtype = out_dtype)
+
+    # for jj in np.arange(len(a)):
+    #     vol[a[jj], b[jj], c[jj]] = u_values[jj].astype(out_dtype)
+
+    # # vol[vol<=0.] = 0 # nonsensical values
+
+    vol = sn.zoom(old_volume, all_upsampling)
+
+    print('Done.')
+
+    # update voxel sizes in header
+    if len(header.get_zooms())==3:
+        new_zooms = tuple( [header.get_zooms()[ii]/float(all_upsampling[ii]) for ii in np.arange(3)] ) # 3 spatial dimensions
+    elif len(header.get_zooms())>3:
+        tmp = [header.get_zooms()[ii]/float(all_upsampling[ii]) for ii in np.arange(3)]
+        tmp.extend(list(header.get_zooms()[3:]))
+        new_zooms = tuple(tmp) # 3 spatial dimensions + 1 time
+    else:
+        print('Cannot handle this stuff... ')
+        print(header.get_zooms())
+        raise Exception('Header has less than 2 entries. 2D?')
+
+    header.set_zooms(new_zooms)
+
+    # adapt affine according to scaling
+    all_upsampling.extend([1.]) # time
+    scaling = np.diag(1./np.asarray(all_upsampling))
+    affine = np.dot(affine, scaling)
+
+    # create new NII
+    newNii = nib.Nifti1Image(vol.astype(out_dtype), header=header, affine=affine)
+
+    # save niftis
+    newNii.to_filename(upsampledFile)
+
+    # save upsampling factors
+    a=0
+    b=0
+    c=0
+    with open(zoom_values_file, 'wb') as outfile:
+        pickle.dump([np.unique(a),np.unique(b),np.unique(c),all_upsampling[:-1],polynomial, old_shape], outfile)
+
+
+    return (newNii)
+
+def downsample(niiFileName, downsampled_file, zoom_values_file='upsampling_log.pickle', order=3):
+    """
+    downsample a nifti which has been upsampled with the function above.
+
+    Example arguments:
+    niiFileName = '10529_t1_upsampled.nii.gz'
+    downsample_file = '10529_t1_downsample.nii.gz'
+    zoom_values_file = 'upsampling_log.pickle'
+    """
+
+    # load the nifti
+    nii = nib.load(niiFileName)
+    header = nii.get_header()
+
+    # make data type of image to float 
+    out_dtype = np.float32
+    header['datatype'] = 16 # corresponds to float32
+    header['bitpix'] = 32 # corresponds to float32
+
+    downsample_factor=[]
+    with open(zoom_values_file,'rb') as zfile:
+        [a, b, c, all_upsampling, polynomial, old_shape] = pickle.load(zfile)
+
+    print('Downsampling with scales: ' + str(1./np.asarray(all_upsampling)))
+    if old_shape:
+        current_shape = nii.get_data().shape
+        print(old_shape)
+        downsample_values = np.asarray([float(old_shape[ii])/float(current_shape[ii]) for ii in np.arange(len(old_shape))])
+    else:
+        if len(all_upsampling) == 1:
+            downsample_values = 1./np.asarray(3*all_upsampling)
+        else:
+            downsample_values = 1./np.asarray(all_upsampling)
+
+    # #prepping for loop
+    # all_coords = [a,b,c]
+    # # print(np.arange(np.round(np.max(coords))))
+
+    # downsampling image
+    print('Downsampling volume...')
+    # downsample_indices = []
+    # for idx, factor in enumerate(downsample_values):
+    #     print('%f, %f'%(idx, factor))
+    #     coords = (all_coords[idx]*factor)
+    #     downsample_idx = []
+    #     for jj in np.arange(np.round(np.max(coords))):
+    #         downsample_idx.append(np.where(coords==jj)[0])
+
+    #     downsample_indices.append(downsample_idx)
+
+    # # TODO: change this to use slices notation
+    # vol = nii.get_data().astype(out_dtype)[np.squeeze(np.asarray(downsample_indices[0])),:,:]
+    # vol = vol[:,np.squeeze(np.asarray(downsample_indices[1])),:]
+    # vol = vol[:,:,np.squeeze(np.asarray(downsample_indices[2]))]
+    # if np.sum(np.abs([vol.shape[ii] - old_shape[ii] for ii in np.arange(len(old_shape))])) != 0:
+    #     print('Downsampled output shape not the same as target.')
+    vol = sn.zoom(nii.get_data(), downsample_values, order=int(order))
+    print('Done.')
+
+    # update voxel sizes in header
+    if len(header.get_zooms())==3:
+        new_zooms = tuple( [header.get_zooms()[ii]/float(downsample_values[ii]) for ii in np.arange(3)] ) # 3 spatial dimensions
+    elif len(header.get_zooms())>3:
+        tmp = [header.get_zooms()[ii]/float(downsample_values[ii]) for ii in np.arange(3)]
+        tmp.extend(list(header.get_zooms()[3:]))
+        new_zooms = tuple(tmp) # 3 spatial dimensions + 1 time
+    else:
+        print('Cannot handle this stuff... ')
+        print(header.get_zooms())
+        raise Exception('Header has less than 2 entries. 2D?')
+
+    # new_zooms = tuple( [header.get_zooms()[ii]/float(downsample_values[ii]) for ii in np.arange(len(header.get_zooms()))] )
+    header.set_zooms(new_zooms)
+
+    # adapt affine according to scaling
+    affine = nii.get_affine()
+    downsample_values = downsample_values.tolist()
+    downsample_values.extend([1.]) # time
+    scaling = np.diag(1./np.asarray(downsample_values))
+    affine = np.dot(affine, scaling)
+
+    # create new NII
+    newNii = nib.Nifti1Image(vol.astype(out_dtype), header=header, affine=affine)
+
+    # save niftis
+    newNii.to_filename(downsampled_file)
+
+    return (newNii)

niitools/util.py

@@ -1,5 +1,21 @@
+import nibabel as nib
+
 class SliceMaker(object):
     def __getitem__(self, slic):
         return slic
 
 slicer = SliceMaker()
+
+def binarize(infile, outfile, threshold):
+
+    # load file
+    nii = nib.load(infile)
+    vol = nii.get_data()
+
+    # create new NII
+    newNii = nib.Nifti1Image((vol>=float(threshold)).astype(int), header=nii.get_header(), affine=nii.get_affine())
+
+    # save niftis
+    newNii.to_filename(outfile)
+
+    return(newNii)

niitools/volumes.py

@@ -95,7 +95,10 @@ def _masked_threshold(infile, threshold, outfile, mode='scalar', excludefile='-'
     if type(threshold) is str:
         threshold = ast.literal_eval(threshold)
     nii = nib.load(infile)
-    data = nii.get_data()
+    data = nii.get_data().astype(np.float32)
+    header = nii.get_header()
+    header['datatype'] = 16 # corresponds to float32
+    header['bitpix'] = 32 # corresponds to float32
 
     if labelfile == '-':
         label_arr = np.ones(data.shape)
@@ -129,7 +132,7 @@ def _masked_threshold(infile, threshold, outfile, mode='scalar', excludefile='-'
         with open(outfile, 'w') as f:
             f.write(str(count))
     else:
-        out = nib.Nifti1Image(vol, header=nii.get_header(),
+        out = nib.Nifti1Image(vol, header=header,
                 affine=nii.get_affine())
         out.to_filename(outfile)
 
@@ -240,16 +243,18 @@ def pad(niiFileName, paddedNiiFileName, maskNiiFileName, padAmountMM='30'):
 
     # load the nifti
     nii = nib.load(niiFileName)
+    header = nii.get_header()
+    out_dtype = np.float32
 
     # get the amount of padding in voxels
-    pixdim = nii.get_header()['pixdim'][1:4]
+    pixdim = np.asarray(header.get_zooms())
     padAmount = np.ceil(padAmountMM / pixdim)
-    dims = nii.get_header()['dim'][1:4]
+    dims = np.asarray(header['dim'][1:4]).astype(int)
     assert np.all(dims.shape == padAmount.shape)
     newDims = dims + padAmount * 2
 
     # compute where the center is for padding
-    center = newDims/2
+    center = (newDims/2.).astype(np.int)
     starts = np.round(center - dims/2)
     ends = starts + dims
 
@@ -258,16 +263,22 @@ def pad(niiFileName, paddedNiiFileName, maskNiiFileName, padAmountMM='30'):
 
     # set the subvolume in the center of the image w/the padding around it
     vol = np.zeros(newDims)
-    vol[slicer] = nii.get_data()
+    vol[slicer] = nii.get_data().astype(out_dtype)
     volMask = np.zeros(newDims)
     volMask[slicer] = np.ones(dims)
 
     # update affine
     affine = nii.get_affine()
     affine[:3, 3] -= padAmountMM
+
+    # update header info
+    header['dim'][1:4] = newDims
+    header['datatype'] = 16 # corresponds to float32
+    header['bitpix'] = 32 # corresponds to float32
+
     # create niftis
-    newNii = nib.Nifti1Image(vol, header=nii.get_header(), affine=affine)
-    newNiiMask = nib.Nifti1Image(volMask, header=nii.get_header(), affine=affine)
+    newNii = nib.Nifti1Image(vol.astype(out_dtype), header=header, affine=affine)
+    newNiiMask = nib.Nifti1Image(volMask.astype(out_dtype), header=header, affine=affine)
 
     # save niftis
     newNii.to_filename(paddedNiiFileName)