diff --git a/skimage/exposure/exposure.py b/skimage/exposure/exposure.py index 2d3b82f5..a0a576d6 100644 --- a/skimage/exposure/exposure.py +++ b/skimage/exposure/exposure.py @@ -174,6 +174,8 @@ def rescale_intensity(image, in_range=None, out_range=None): if out_range is None: omin, omax = dtype_range[dtype] + if imin >= 0: + omin = 0 else: omin, omax = out_range diff --git a/skimage/util/dtype.py b/skimage/util/dtype.py index e9bc3036..44476e44 100644 --- a/skimage/util/dtype.py +++ b/skimage/util/dtype.py @@ -10,8 +10,8 @@ dtype_range = {np.uint8: (0, 255), np.uint16: (0, 65535), np.int8: (-128, 127), np.int16: (-32768, 32767), - np.float32: (0, 1), - np.float64: (0, 1)} + np.float32: (-1, 1), + np.float64: (-1, 1)} integer_types = (np.uint8, np.uint16, np.int8, np.int16) @@ -20,19 +20,24 @@ _supported_types = (np.uint8, np.uint16, np.uint32, np.float32, np.float64) if np.__version__ >= "1.6.0": - dtype_range[np.float16] = (0, 1) + dtype_range[np.float16] = (-1, 1) _supported_types += (np.float16, ) -def convert(image, dtype, force_copy=False): +def convert(image, dtype, force_copy=False, uniform=False): """ Convert an image to the requested data-type. - Warnings are issued in case of precision loss, or when - negative values have to be scaled into the positive domain. - Floating point values must be in the range [0.0, 1.0]. + Warnings are issued in case of precision loss, or when negative values + are clipped during conversion to unsigned integer types (sign loss). + + Floating point values are expected to be normalized and will be clipped + to the range [0.0, 1.0] or [-1.0, 1.0] when converting to unsigned or + signed integers respectively. + Numbers are not shifted to the negative side when converting from - floating point or unsigned integer types to signed integer types. + unsigned to signed integer types. Negative values will be clipped when + converting to unsigned integers. Parameters ---------- @@ -42,17 +47,29 @@ def convert(image, dtype, force_copy=False): Target data-type. force_copy : bool Force a copy of the data, irrespective of its current dtype. + uniform : bool + Uniformly quantize the floating point range to the integer range. + By default (uniform=False) floating point values are scaled and + rounded to the nearest integers, which minimizes back and forth + conversion errors. + + References + ---------- + (1) DirectX data conversion rules. + http://msdn.microsoft.com/en-us/library/windows/desktop/dd607323%28v=vs.85%29.aspx + (2) Data Conversions. + In "OpenGL ES 2.0 Specification v2.0.25", pp 7-8. Khronos Group, 2010. + (3) Proper treatment of pixels as integers. A.W. Paeth. + In "Graphics Gems I", pp 249-256. Morgan Kaufmann, 1990. + (4) Dirty Pixels. J. Blinn. + In "Jim Blinn's corner: Dirty Pixels", pp 47-57. Morgan Kaufmann, 1998. """ image = np.asarray(image) - dtype = np.dtype(dtype).type - dtype_in = image.dtype.type dtypeobj = np.dtype(dtype) - dtypeobj_in = np.dtype(dtype_in) - kind = dtypeobj.kind - kind_in = dtypeobj_in.kind - itemsize = dtypeobj.itemsize - itemsize_in = dtypeobj_in.itemsize + dtypeobj_in = image.dtype + dtype = dtypeobj.type + dtype_in = dtypeobj_in.type if dtype_in == dtype: if force_copy: @@ -74,93 +91,142 @@ def convert(image, dtype, force_copy=False): # Return first of `dtypes` with itemsize greater than `itemsize` return next(dt for dt in dtypes if itemsize < np.dtype(dt).itemsize) + def _dtype2(kind, bits, itemsize=1): + # Return dtype of `kind` that can store a `bits` wide unsigned int + c = lambda x, y: x <= y if kind == 'u' else x < y + s = next(i for i in (itemsize, ) + (2, 4, 8) if c(bits, i*8)) + return np.dtype(kind + str(s)) + + def _scale(a, n, m, copy=True): + # Scale unsigned/positive integers from n to m bits + # Numbers can be represented exactly only if m is a multiple of n + # Output array is of same kind as input. + kind = a.dtype.kind + if n == m: + return a.copy() if copy else a + elif n > m: + # downscale with precision loss + prec_loss() + if copy: + b = np.empty(a.shape, _dtype2(kind, m)) + np.divide(a, 2**(n - m), out=b, dtype=a.dtype, + casting='unsafe') + return b + else: + a //= 2**(n - m) + return a + elif m % n == 0: + # exact upscale to a multiple of n bits + if copy: + b = np.empty(a.shape, _dtype2(kind, m)) + np.multiply(a, (2**m - 1) // (2**n - 1), out=b, dtype=b.dtype) + return b + else: + a = np.array(a, _dtype2(kind, m, a.dtype.itemsize), copy=False) + a *= (2**m - 1) // (2**n - 1) + return a + else: + # upscale to a multiple of n bits, + # then downscale with precision loss + prec_loss() + o = (m // n + 1) * n + if copy: + b = np.empty(a.shape, _dtype2(kind, o)) + np.multiply(a, (2**o - 1) // (2**n - 1), out=b, dtype=b.dtype) + b //= 2**(o - m) + return b + else: + a = np.array(a, _dtype2(kind, o, a.dtype.itemsize), copy=False) + a *= (2**o - 1) // (2**n - 1) + a //= 2**(o - m) + return a + + kind = dtypeobj.kind + kind_in = dtypeobj_in.kind + itemsize = dtypeobj.itemsize + itemsize_in = dtypeobj_in.itemsize + if kind in 'ui': + imin = np.iinfo(dtype).min + imax = np.iinfo(dtype).max + if kind_in in 'ui': + imin_in = np.iinfo(dtype_in).min + imax_in = np.iinfo(dtype_in).max + if kind_in == 'f': - if np.min(image) < 0 or np.max(image) > 1: - raise ValueError("Images of type float must be between 0 and 1") + if np.min(image) < -1.0 or np.max(image) > 1.0: + raise ValueError("Images of type float must be between -1 and 1.") if kind == 'f': # floating point -> floating point if itemsize_in > itemsize: prec_loss() return dtype(image) + # floating point -> integer prec_loss() # use float type that can represent output integer type image = np.array(image, _dtype(itemsize, dtype_in, np.float32, np.float64)) - image *= np.iinfo(dtype).max + 1 - np.clip(image, 0, np.iinfo(dtype).max, out=image) + if not uniform: + if kind == 'u': + image *= imax + else: + image *= imax - imin + image -= 1.0 + image /= 2.0 + np.rint(image, out=image) + np.clip(image, imin, imax, out=image) + elif kind == 'u': + image *= imax + 1 + np.clip(image, 0, imax, out=image) + else: + image *= (imax - imin + 1.0) / 2.0 + np.floor(image, out=image) + np.clip(image, imin, imax, out=image) return dtype(image) + if kind == 'f': # integer -> floating point if itemsize_in >= itemsize: prec_loss() - # use float type that can represent input integers + # use float type that can exactly represent input integers image = np.array(image, _dtype(itemsize_in, dtype, np.float32, np.float64)) - if np.iinfo(dtype_in).min: - sign_loss() - image -= np.iinfo(dtype_in).min - image /= np.iinfo(dtype_in).max - np.iinfo(dtype_in).min - return dtype(image) - if kind_in == 'u': - # unsigned integer -> integer - shift = 1 if kind == 'i' else 0 - if itemsize_in > itemsize: - prec_loss() - image = image >> 8 * (itemsize_in - itemsize) + shift - return dtype(image) - result = dtype(image) - result <<= 8 * (itemsize - itemsize_in) - shift - if itemsize - itemsize_in == 3: - # uint8 -> (u)int32 - # hint: 4294967295 == (255 << 24) + (255 << 16) + (255 << 8) + 255 - image = dtype(image) - image *= 2**16 + 2**8 + 1 - if shift: - result += image >> shift + if kind_in == 'u': + image /= imax_in + # DirectX uses this conversion also for signed ints + #if imin_in: + # np.maximum(image, -1.0, out=image) else: - result += image - return dtype(result) + image *= 2.0 + image += 1.0 + image /= imax_in - imin_in + return dtype(image) + + if kind_in == 'u': + if kind == 'i': + # unsigned integer -> signed integer + image = _scale(image, 8*itemsize_in, 8*itemsize-1) + return image.view(dtype) + else: + # unsigned integer -> unsigned integer + return _scale(image, 8*itemsize_in, 8*itemsize) + if kind == 'u': # signed integer -> unsigned integer sign_loss() - # use next higher precision signed integer type - image = np.array(image, _dtype(itemsize_in, - np.int16, np.int32, np.int64)) - image -= np.iinfo(dtype_in).min - if itemsize_in == itemsize: - return dtype(image) - if itemsize_in > itemsize: - prec_loss() - image >>= 8 * (itemsize_in - itemsize) - return dtype(image) - result = dtype(image) - result <<= 8 * (itemsize - itemsize_in) - if itemsize - itemsize_in == 3: - # int8 -> uint32 - image = dtype(image) - image *= 2**16 + 2**8 + 1 - result += dtype(image) + image = _scale(image, 8*itemsize_in-1, 8*itemsize) + result = np.empty(image.shape, dtype) + np.maximum(image, 0, out=result, dtype=image.dtype, casting='unsafe') return result - if kind == 'i': - # signed integer -> signed integer - if itemsize_in > itemsize: - prec_loss() - return dtype(image // 2**(8 * (itemsize_in - itemsize))) - # use next higher precision signed integer type - image = np.array(image, _dtype(itemsize_in, - np.int16, np.int32, np.int64)) - image -= np.iinfo(dtype_in).min - # use next higher precision signed integer type - result = np.array(image, _dtype(image.itemsize, np.int32, np.int64)) - result *= 2**(8 * (itemsize - itemsize_in)) - if itemsize - itemsize_in == 3: - # int8 -> int32 - image = dtype(image) - image *= 2**16 + 2**8 + 1 - result += image - result += np.iinfo(dtype).min - return dtype(result) + + # signed integer -> signed integer + if itemsize_in > itemsize: + return _scale(image, 8*itemsize_in-1, 8*itemsize-1) + image = image.astype(_dtype2('i', itemsize*8)) + image -= imin_in + image = _scale(image, 8*itemsize_in, 8*itemsize, copy=False) + image += imin + return dtype(image) def img_as_float(image, force_copy=False): diff --git a/skimage/util/tests/test_dtype.py b/skimage/util/tests/test_dtype.py index d5067981..2ef0044f 100644 --- a/skimage/util/tests/test_dtype.py +++ b/skimage/util/tests/test_dtype.py @@ -9,12 +9,13 @@ dtype_range = {np.uint8: (0, 255), np.uint16: (0, 65535), np.int8: (-128, 127), np.int16: (-32768, 32767), - np.float32: (0, 1), - np.float64: (0, 1)} + np.float32: (-1.0, 1.0), + np.float64: (-1.0, 1.0)} -def _verify_range(msg, x, vmin, vmax): +def _verify_range(msg, x, vmin, vmax, dtype): assert_equal(x[0], vmin) assert_equal(x[-1], vmax) + assert x.dtype == dtype def test_range(): for dtype in dtype_range: @@ -29,12 +30,13 @@ def test_range(): omin, omax = dtype_range[dt] - if imin == 0: + if imin == 0 or omin == 0: omin = 0 + imin = 0 - yield _verify_range, \ - "From %s to %s" % (np.dtype(dtype), np.dtype(dt)), \ - y, omin, omax + yield (_verify_range, + "From %s to %s" % (np.dtype(dtype), np.dtype(dt)), + y, omin, omax, np.dtype(dt)) def test_range_extra_dtypes(): @@ -57,9 +59,9 @@ def test_range_extra_dtypes(): x = np.linspace(imin, imax, 10).astype(dtype_in) y = convert(x, dt) omin, omax = dtype_range_extra[dt] - yield _verify_range, \ - "From %s to %s" % (np.dtype(dtype_in), np.dtype(dt)), \ - y, omin, omax + yield (_verify_range, + "From %s to %s" % (np.dtype(dtype_in), np.dtype(dt)), + y, omin, omax, np.dtype(dt)) def test_unsupported_dtype(): @@ -70,7 +72,7 @@ def test_unsupported_dtype(): def test_float_out_of_range(): too_high = np.array([2], dtype=np.float32) assert_raises(ValueError, img_as_int, too_high) - too_low = np.array([-1], dtype=np.float32) + too_low = np.array([-2], dtype=np.float32) assert_raises(ValueError, img_as_int, too_low) @@ -84,4 +86,3 @@ def test_copy(): if __name__ == '__main__': np.testing.run_module_suite() -