From 8a0b743686426bc0c5fac203c01df90e5bef2809 Mon Sep 17 00:00:00 2001 From: MDS Date: Tue, 20 Sep 2016 10:32:09 -0400 Subject: [PATCH 1/4] advanced_tools.py volumes.py Explicitly account for data type after loading nifti images. Due to different standards at different sites an overflow issue occured, messing up the resulting image. resampling.py Added bicubic upsampling. niitools.py Adjusted which upsampling/downsampling method to use --- niitools.py | 3 +- niitools/advanced_tools.py | 13 ++- niitools/resampling.py | 193 +++++++++++++++++++++++++++++++++++++ niitools/volumes.py | 25 +++-- 4 files changed, 223 insertions(+), 11 deletions(-) mode change 100644 => 100755 niitools/resampling.py diff --git a/niitools.py b/niitools.py index 8c195cc..ed7254a 100755 --- a/niitools.py +++ b/niitools.py @@ -20,7 +20,8 @@ '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, diff --git a/niitools/advanced_tools.py b/niitools/advanced_tools.py index 6e851bd..b72cbd6 100755 --- a/niitools/advanced_tools.py +++ b/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) diff --git a/niitools/resampling.py b/niitools/resampling.py old mode 100644 new mode 100755 index 7b00785..622f3bc --- a/niitools/resampling.py +++ b/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,15 @@ from numpy import newaxis as nax import nibabel as nib +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 +140,186 @@ 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) + print('Upsampling with scales: ' + str(all_upsampling)) + 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 + old_shape = old_volume.shape + new_shape = tuple([old_shape[ii]*usampling for ii, usampling in enumerate(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 + + 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 + with open(zoom_values_file, 'w') 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'): + """ + 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,'r') as zfile: + [a, b, c, all_upsampling, polynomial, old_shape] = pickle.load(zfile) + + print('Downsampling with scales: ' + str(1./np.asarray(all_upsampling))) + 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] + + # 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]))] + 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) diff --git a/niitools/volumes.py b/niitools/volumes.py index 2256f21..6c2d777 100644 --- a/niitools/volumes.py +++ b/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,11 +243,13 @@ 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 = header['dim'][1:4] assert np.all(dims.shape == padAmount.shape) newDims = dims + padAmount * 2 @@ -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) From 474a891422c578a50e785cdc8186f15f76e91c1b Mon Sep 17 00:00:00 2001 From: MDS Date: Tue, 12 Mar 2019 11:12:56 -0400 Subject: [PATCH 2/4] updated pyniftitools based on WMH pipeline --- .gitignore | 0 LICENSE | 0 README.md | 0 niitools/__init__.py | 0 niitools/headers.py | 0 niitools/labels.py | 0 niitools/resampling.py | 87 +++++++++++++++++++++++++----------------- niitools/util.py | 0 niitools/volumes.py | 4 +- 9 files changed, 55 insertions(+), 36 deletions(-) mode change 100644 => 100755 .gitignore mode change 100644 => 100755 LICENSE mode change 100644 => 100755 README.md mode change 100644 => 100755 niitools/__init__.py mode change 100644 => 100755 niitools/headers.py mode change 100644 => 100755 niitools/labels.py mode change 100644 => 100755 niitools/util.py mode change 100644 => 100755 niitools/volumes.py diff --git a/.gitignore b/.gitignore old mode 100644 new mode 100755 diff --git a/LICENSE b/LICENSE old mode 100644 new mode 100755 diff --git a/README.md b/README.md old mode 100644 new mode 100755 diff --git a/niitools/__init__.py b/niitools/__init__.py old mode 100644 new mode 100755 diff --git a/niitools/headers.py b/niitools/headers.py old mode 100644 new mode 100755 diff --git a/niitools/labels.py b/niitools/labels.py old mode 100644 new mode 100755 diff --git a/niitools/resampling.py b/niitools/resampling.py index 622f3bc..035fdf2 100755 --- a/niitools/resampling.py +++ b/niitools/resampling.py @@ -8,6 +8,7 @@ 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 @@ -180,7 +181,6 @@ def upsample(niiFileName, upsampledFile, zoom_values_file='upsampling_log.pickle all_upsampling[idx]= 1. else: all_upsampling[idx]= np.round(zoom) - print('Upsampling with scales: ' + str(all_upsampling)) else: upsample_factor = float(upsample_factor) all_upsampling = [upsample_factor for zoom in header.get_zooms()[0:3]] @@ -189,26 +189,33 @@ def upsample(niiFileName, upsampledFile, zoom_values_file='upsampling_log.pickle polynomial = int(polynomial) old_volume = np.squeeze(nii.get_data().astype(float)) - # get new volume shape + # get new volume shape and update upsampling based on rounded numbers old_shape = old_volume.shape - new_shape = tuple([old_shape[ii]*usampling for ii, usampling in enumerate(all_upsampling)]) - # get indices - a,b,c = np.indices(new_shape) - a = a.ravel() - b = b.ravel() - c = c.ravel() + 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) + # 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) + # 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[vol<=0.] = 0 # nonsensical values + + vol = sn.zoom(old_volume, all_upsampling) print('Done.') @@ -238,13 +245,16 @@ def upsample(niiFileName, upsampledFile, zoom_values_file='upsampling_log.pickle newNii.to_filename(upsampledFile) # save upsampling factors + a=0 + b=0 + c=0 with open(zoom_values_file, 'w') 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'): +def downsample(niiFileName, downsampled_file, zoom_values_file='upsampling_log.pickle', order=3): """ downsample a nifti which has been upsampled with the function above. @@ -268,30 +278,39 @@ def downsample(niiFileName, downsampled_file, zoom_values_file='upsampling_log.p [a, b, c, all_upsampling, polynomial, old_shape] = pickle.load(zfile) print('Downsampling with scales: ' + str(1./np.asarray(all_upsampling))) - if len(all_upsampling) == 1: - downsample_values = 1./np.asarray(3*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: - downsample_values = 1./np.asarray(all_upsampling) + 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] + # #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 = [] + # 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]))] + # 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 diff --git a/niitools/util.py b/niitools/util.py old mode 100644 new mode 100755 diff --git a/niitools/volumes.py b/niitools/volumes.py old mode 100644 new mode 100755 index 6c2d777..184afd8 --- a/niitools/volumes.py +++ b/niitools/volumes.py @@ -249,12 +249,12 @@ def pad(niiFileName, paddedNiiFileName, maskNiiFileName, padAmountMM='30'): # get the amount of padding in voxels pixdim = np.asarray(header.get_zooms()) padAmount = np.ceil(padAmountMM / pixdim) - dims = 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 From 6fc098b722cc044cd2765a4c3f0284633d9d1c39 Mon Sep 17 00:00:00 2001 From: MDS Date: Mon, 29 Apr 2019 12:35:13 -0400 Subject: [PATCH 3/4] update after ncerebro integration --- niitools.py | 1 + niitools/util.py | 16 ++++++++++++++++ 2 files changed, 17 insertions(+) diff --git a/niitools.py b/niitools.py index ed7254a..6201128 100755 --- a/niitools.py +++ b/niitools.py @@ -15,6 +15,7 @@ '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, diff --git a/niitools/util.py b/niitools/util.py index d3dc74e..b0e4e9a 100755 --- a/niitools/util.py +++ b/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) \ No newline at end of file From 8e98ea1e5733cd85b53369b0e00d51b5134c57dd Mon Sep 17 00:00:00 2001 From: MDS Date: Mon, 29 Apr 2019 14:38:11 -0400 Subject: [PATCH 4/4] conversion to python 3 --- niitools/resampling.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/niitools/resampling.py b/niitools/resampling.py index 035fdf2..189d7e2 100755 --- a/niitools/resampling.py +++ b/niitools/resampling.py @@ -248,7 +248,7 @@ def upsample(niiFileName, upsampledFile, zoom_values_file='upsampling_log.pickle a=0 b=0 c=0 - with open(zoom_values_file, 'w') as outfile: + 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) @@ -274,7 +274,7 @@ def downsample(niiFileName, downsampled_file, zoom_values_file='upsampling_log.p header['bitpix'] = 32 # corresponds to float32 downsample_factor=[] - with open(zoom_values_file,'r') as zfile: + 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)))