Skip to content

Commit

Permalink
Merge pull request #113 from StellarStorm/numpy_dtypes
Browse files Browse the repository at this point in the history
Update datatypes for Numpy 1.24

Numpy 1.24 deprecated

    np.bool
    np.int
    np.float

Attempting to run MedPy raised an AttributeError that these weren't installed (I found this for metrics.hd95) This replaces them with np.bool_, int, and float, respectively, so medpy can continue to work with the latest Numpy.

Casts to specific precision (e.g. np.float32) are unchanged
  • Loading branch information
StellarStorm authored Dec 1, 2023
2 parents 66265de + 00ebacb commit 3955fe2
Show file tree
Hide file tree
Showing 21 changed files with 1,122 additions and 1,123 deletions.
36 changes: 18 additions & 18 deletions bin/medpy_apparent_diffusion_coefficient.py
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@
Copyright (C) 2013 Oskar Maier
This program comes with ABSOLUTELY NO WARRANTY; This is free software,
and you are welcome to redistribute it under certain conditions; see
the LICENSE file or <http://www.gnu.org/licenses/> for details.
the LICENSE file or <http://www.gnu.org/licenses/> for details.
"""

# code
Expand All @@ -92,27 +92,27 @@ def main():
logger = Logger.getInstance()
if args.debug: logger.setLevel(logging.DEBUG)
elif args.verbose: logger.setLevel(logging.INFO)

# loading input images
b0img, b0hdr = load(args.b0image)
bximg, bxhdr = load(args.bximage)

# convert to float
b0img = b0img.astype(numpy.float)
bximg = bximg.astype(numpy.float)
b0img = b0img.astype(float)
bximg = bximg.astype(float)

# check if image are compatible
if not b0img.shape == bximg.shape:
raise ArgumentError('The input images shapes differ i.e. {} != {}.'.format(b0img.shape, bximg.shape))
if not header.get_pixel_spacing(b0hdr) == header.get_pixel_spacing(bxhdr):
raise ArgumentError('The input images voxel spacing differs i.e. {} != {}.'.format(header.get_pixel_spacing(b0hdr), header.get_pixel_spacing(bxhdr)))

# check if supplied threshold value as well as the b value is above 0
if args.threshold is not None and not args.threshold >= 0:
raise ArgumentError('The supplied threshold value must be greater than 0, otherwise a division through 0 might occur.')
if not args.b > 0:
raise ArgumentError('The supplied b-value must be greater than 0.')

# compute threshold value if not supplied
if args.threshold is None:
b0thr = otsu(b0img, 32) / 4. # divide by 4 to decrease impact
Expand All @@ -123,29 +123,29 @@ def main():
raise ArgumentError('The supplied bximage seems to contain negative values.')
else:
b0thr = bxthr = args.threshold

logger.debug('thresholds={}/{}, b-value={}'.format(b0thr, bxthr, args.b))

# threshold b0 + bx DW image to obtain a mask
# b0 mask avoid division through 0, bx mask avoids a zero in the ln(x) computation
mask = binary_fill_holes(b0img > b0thr) & binary_fill_holes(bximg > bxthr)

# perform a number of binary morphology steps to select the brain only
mask = binary_erosion(mask, iterations=1)
mask = largest_connected_component(mask)
mask = binary_dilation(mask, iterations=1)

logger.debug('excluding {} of {} voxels from the computation and setting them to zero'.format(numpy.count_nonzero(mask), numpy.prod(mask.shape)))

# compute the ADC
adc = numpy.zeros(b0img.shape, b0img.dtype)
adc[mask] = -1. * args.b * numpy.log(bximg[mask] / b0img[mask])
adc[adc < 0] = 0

# saving the resulting image
save(adc, args.output, b0hdr, args.force)


def getArguments(parser):
"Provides additional validation of the arguments collected by argparse."
return parser.parse_args()
Expand All @@ -157,13 +157,13 @@ def getParser():
parser.add_argument('bximage', help='the diffusion weighted image required with b=x')
parser.add_argument('b', type=int, help='the b-value used to acquire the bx-image (i.e. x)')
parser.add_argument('output', help='the computed apparent diffusion coefficient image')

parser.add_argument('-t', '--threshold', type=int, dest='threshold', help='set a fixed threshold for the input images to mask the computation')

parser.add_argument('-v', '--verbose', dest='verbose', action='store_true', help='verbose output')
parser.add_argument('-d', dest='debug', action='store_true', help='Display debug information.')
parser.add_argument('-f', '--force', dest='force', action='store_true', help='overwrite existing files')
return parser

if __name__ == "__main__":
main()
main()
78 changes: 39 additions & 39 deletions bin/medpy_binary_resampling.py
Original file line number Diff line number Diff line change
Expand Up @@ -58,7 +58,7 @@
Copyright (C) 2013 Oskar Maier
This program comes with ABSOLUTELY NO WARRANTY; This is free software,
and you are welcome to redistribute it under certain conditions; see
the LICENSE file or <http://www.gnu.org/licenses/> for details.
the LICENSE file or <http://www.gnu.org/licenses/> for details.
"""

# code
Expand All @@ -70,59 +70,59 @@ def main():
logger = Logger.getInstance()
if args.debug: logger.setLevel(logging.DEBUG)
elif args.verbose: logger.setLevel(logging.INFO)

# loading input images
img, hdr = load(args.input)
img = img.astype(numpy.bool)
img = img.astype(numpy.bool_)

# check spacing values
if not len(args.spacing) == img.ndim:
parser.error('The image has {} dimensions, but {} spacing parameters have been supplied.'.format(img.ndim, len(args.spacing)))

# check if output image exists
if not args.force:
if os.path.exists(args.output):
parser.error('The output image {} already exists.'.format(args.output))
parser.error('The output image {} already exists.'.format(args.output))

logger.debug('target voxel spacing: {}'.format(args.spacing))

# determine number of required complete slices for up-sampling
vs = header.get_pixel_spacing(hdr)
rcss = [int(y // x - 1) for x, y in zip(args.spacing, vs)] # TODO: For option b, remove the - 1; better: no option b, since I am rounding later anyway

# remove negatives and round up to next even number
rcss = [x if x > 0 else 0 for x in rcss]
rcss = [x if 0 == x % 2 else x + 1 for x in rcss]
logger.debug('intermediate slices to add per dimension: {}'.format(rcss))

# for each dimension requiring up-sampling, from the highest down, perform shape based slice interpolation
logger.info('Adding required slices using shape based interpolation.')
for dim, rcs in enumerate(rcss):
if rcs > 0:
logger.debug('adding {} intermediate slices to dimension {}'.format(rcs, dim))
img = shape_based_slice_interpolation(img, dim, rcs)
logger.debug('resulting new image shape: {}'.format(img.shape))

# compute and set new voxel spacing
nvs = [x / (y + 1.) for x, y in zip(vs, rcss)]
header.set_pixel_spacing(hdr, nvs)
logger.debug('intermediate voxel spacing: {}'.format(nvs))

# interpolate with nearest neighbour
logger.info('Re-sampling the image with a b-spline order of {}.'.format(args.order))
img, hdr = resample(img, hdr, args.spacing, args.order, mode='nearest')

# saving the resulting image
save(img, args.output, hdr, args.force)

def shape_based_slice_interpolation(img, dim, nslices):
"""
Adds `nslices` slices between all slices of the binary image `img` along dimension
`dim` respecting the original slice values to be situated in the middle of each
slice. Extrapolation situations are handled by simple repeating.
Interpolation of new slices is performed using shape based interpolation.
Parameters
----------
img : array_like
Expand All @@ -131,7 +131,7 @@ def shape_based_slice_interpolation(img, dim, nslices):
The dimension along which to add slices.
nslices : int
The number of slices to add. Must be an even number.
Returns
-------
out : ndarray
Expand All @@ -140,32 +140,32 @@ def shape_based_slice_interpolation(img, dim, nslices):
# check arguments
if not 0 == nslices % 2:
raise ValueError('nslices must be an even number')

out = None
slicer = [slice(None)] * img.ndim
chunk_full_shape = list(img.shape)
chunk_full_shape[dim] = nslices + 2

for sl1, sl2 in zip(numpy.rollaxis(img, dim)[:-1], numpy.rollaxis(img, dim)[1:]):
if 0 == numpy.count_nonzero(sl1) and 0 == numpy.count_nonzero(sl2):
chunk = numpy.zeros(chunk_full_shape, dtype=numpy.bool)
chunk = numpy.zeros(chunk_full_shape, dtype=numpy.bool_)
else:
chunk = shape_based_slice_insertation_object_wise(sl1, sl2, dim, nslices)
if out is None:
out = numpy.delete(chunk, -1, dim)
else:
out = numpy.concatenate((out, numpy.delete(chunk, -1, dim)), dim)
slicer[dim] = numpy.newaxis

slicer[dim] = numpy.newaxis
out = numpy.concatenate((out, sl2[slicer]), dim)

slicer[dim] = slice(0, 1)
for _ in range(nslices // 2):
out = numpy.concatenate((img[slicer], out), dim)
slicer[dim] = slice(-1, None)
for _ in range(nslices // 2):
out = numpy.concatenate((out, img[slicer]), dim)

return out

def shape_based_slice_insertation_object_wise(sl1, sl2, dim, nslices, order=3):
Expand All @@ -184,15 +184,15 @@ def shape_based_slice_insertation_object_wise(sl1, sl2, dim, nslices, order=3):
else:
out |= _out
return out

def shape_based_slice_insertation(sl1, sl2, dim, nslices, order=3):
"""
Insert `nslices` new slices between `sl1` and `sl2` along dimension `dim` using shape
based binary interpolation.
Extrapolation is handled adding `nslices`/2 step-wise eroded copies of the last slice
in each direction.
Parameters
----------
sl1 : array_like
Expand All @@ -205,16 +205,16 @@ def shape_based_slice_insertation(sl1, sl2, dim, nslices, order=3):
The number of slices to add.
order : int
The b-spline interpolation order for re-sampling the distance maps.
Returns
-------
out : ndarray
A binary image of size `sl1`.shape() extend by `nslices`+2 along the new
dimension `dim`. The border slices are the original slices `sl1` and `sl2`.
"""
sl1 = sl1.astype(numpy.bool)
sl2 = sl2.astype(numpy.bool)
sl1 = sl1.astype(numpy.bool_)
sl2 = sl2.astype(numpy.bool_)

# extrapolation through erosion
if 0 == numpy.count_nonzero(sl1):
slices = [sl1]
Expand All @@ -234,26 +234,26 @@ def shape_based_slice_insertation(sl1, sl2, dim, nslices, order=3):
slices.append(sl2)
return numpy.rollaxis(numpy.asarray(slices), 0, dim + 1)
#return numpy.asarray([sl.T for sl in slices]).T
# interpolation shape based

# interpolation shape based
# note: distance_transform_edt shows strange behaviour for ones-arrays
dt1 = distance_transform_edt(~sl1) - distance_transform_edt(sl1)
dt2 = distance_transform_edt(~sl2) - distance_transform_edt(sl2)

slicer = [slice(None)] * dt1.ndim
slicer = slicer[:dim] + [numpy.newaxis] + slicer[dim:]
out = numpy.concatenate((dt1[slicer], dt2[slicer]), axis=dim)
zoom_factors = [1] * dt1.ndim
zoom_factors = zoom_factors[:dim] + [(nslices + 2)/2.] + zoom_factors[dim:]
out = zoom(out, zoom_factors, order=order)

return out <= 0

def getArguments(parser):
"Provides additional validation of the arguments collected by argparse."
args = parser.parse_args()
if args.order < 0 or args.order > 5:
parser.error('The order has to be a number between 0 and 5.')
parser.error('The order has to be a number between 0 and 5.')
return args

def getParser():
Expand All @@ -268,6 +268,6 @@ def getParser():
parser.add_argument('-d', dest='debug', action='store_true', help='Display debug information.')
parser.add_argument('-f', '--force', dest='force', action='store_true', help='overwrite existing files')
return parser

if __name__ == "__main__":
main()
main()
30 changes: 15 additions & 15 deletions bin/medpy_extract_contour.py
Original file line number Diff line number Diff line change
Expand Up @@ -46,16 +46,16 @@
contour width, the surface of the volume will correspond with the
middle of the contour line. In the case of an odd contour width, the
contour will be shifted by one voxel towards the inside of the volume.
In the case of 3D volumes, the contours result in shells, which might
not be desired, as they do not visualize well in 2D views. With the
'--dimension' argument, a dimension along which to extract the contours
can be supplied.
Copyright (C) 2013 Oskar Maier
This program comes with ABSOLUTELY NO WARRANTY; This is free software,
and you are welcome to redistribute it under certain conditions; see
the LICENSE file or <http://www.gnu.org/licenses/> for details.
the LICENSE file or <http://www.gnu.org/licenses/> for details.
"""

# code
Expand All @@ -66,22 +66,22 @@ def main():
logger = Logger.getInstance()
if args.debug: logger.setLevel(logging.DEBUG)
elif args.verbose: logger.setLevel(logging.INFO)

# load input image
data_input, header_input = load(args.input)

# treat as binary
data_input = data_input.astype(numpy.bool)
data_input = data_input.astype(numpy.bool_)

# check dimension argument
if args.dimension and (not args.dimension >= 0 or not args.dimension < data_input.ndim):
argparse.ArgumentError(args.dimension, 'Invalid dimension of {} supplied. Image has only {} dimensions.'.format(args.dimension, data_input.ndim))

# compute erosion and dilation steps
erosions = int(math.ceil(args.width / 2.))
dilations = int(math.floor(args.width / 2.))
logger.debug("Performing {} erosions and {} dilations to achieve a contour of width {}.".format(erosions, dilations, args.width))

# erode, dilate and compute contour
if not args.dimension:
eroded = binary_erosion(data_input, iterations=erosions) if not 0 == erosions else data_input
Expand All @@ -95,17 +95,17 @@ def main():
slicer[args.dimension] = slice(sl, sl+1)
bs_slicer[args.dimension] = slice(1, 2)
bs = generate_binary_structure(data_input.ndim, 1)

eroded = binary_erosion(data_input[slicer], structure=bs[bs_slicer], iterations=erosions) if not 0 == erosions else data_input[slicer]
dilated = binary_dilation(data_input[slicer], structure=bs[bs_slicer], iterations=dilations) if not 0 == dilations else data_input[slicer]
data_output[slicer] = dilated - eroded
logger.debug("Contour image contains {} contour voxels.".format(numpy.count_nonzero(data_output)))

# save resulting volume
save(data_output, args.output, header_input, args.force)
logger.info("Successfully terminated.")

logger.info("Successfully terminated.")

def getArguments(parser):
"Provides additional validation of the arguments collected by argparse."
args = parser.parse_args()
Expand All @@ -123,7 +123,7 @@ def getParser():
parser.add_argument('-v', dest='verbose', action='store_true', help='Display more information.')
parser.add_argument('-d', dest='debug', action='store_true', help='Display debug information.')
parser.add_argument('-f', dest='force', action='store_true', help='Silently override existing output images.')
return parser
return parser

if __name__ == "__main__":
main()
main()
Loading

0 comments on commit 3955fe2

Please sign in to comment.