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.py b/niitools.py index 8c195cc..6201128 100755 --- a/niitools.py +++ b/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, diff --git a/niitools/__init__.py b/niitools/__init__.py old mode 100644 new mode 100755 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/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 old mode 100644 new mode 100755 index 7b00785..189d7e2 --- 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,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) diff --git a/niitools/util.py b/niitools/util.py old mode 100644 new mode 100755 index d3dc74e..b0e4e9a --- 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 diff --git a/niitools/volumes.py b/niitools/volumes.py old mode 100644 new mode 100755 index 2256f21..184afd8 --- 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,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)