Merge pull request #99 from cgohlke/patch-2

ENH: More efficient dtype conversion.

skimage/util/dtype.py

@@ -15,13 +15,20 @@ dtype_range = {np.uint8: (0, 255),
 
 integer_types = (np.uint8, np.uint16, np.int8, np.int16)
 
+_supported_types = (np.uint8, np.uint16, np.uint32,
+                    np.int8, np.int16, np.int32,
+                    np.float16, np.float32, np.float64)
 
-def _convert(image, dtype, prec_loss):
+
+def _convert(image, dtype):
     """
     Convert an image to the requested data-type.
 
-    Warnings are issues in case of precision loss, or when
+    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].
+    Numbers are not shifted to the negative side when converting from
+    floating point or unsigned integer types to signed integer types.
 
     Parameters
     ----------
@@ -29,59 +36,114 @@ def _convert(image, dtype, prec_loss):
         Input image.
     dtype : dtype
         Target data-type.
-    prec_loss : tuple
-        List of input data-types that, when converted to `dtype`,
-        would lose precision.
 
     """
     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
 
     if dtype_in == dtype:
         return image
 
-    if dtype_in in prec_loss:
-        log.warn('Possible precision loss, converting from '
-                 '%s to %s' % (np.dtype(dtype_in), np.dtype(dtype)))
+    if not (dtype_in in _supported_types and dtype in _supported_types):
+        raise ValueError("can not convert %s to %s." % (dtypeobj_in, dtypeobj))
 
-    try:
-        imin, imax = dtype_range[dtype_in]
-        omin, omax = dtype_range[dtype]
-    except KeyError:
-        raise ValueError("Unsure how to convert %s to %s." % \
-                         (np.dtype(dtype_in), np.dtype(dtype)))
+    def sign_loss():
+        log.warn("Possible sign loss when converting negative image of type "
+                 "%s to positive image of type %s." % (dtypeobj_in, dtypeobj))
 
-    sign_loss = (np.sign(imin) == -1) and (np.sign(omin) != -1)
+    def prec_loss():
+        log.warn("Possible precision loss when converting from "
+                 "%s to %s" % (dtypeobj_in, dtypeobj))
 
-    if sign_loss:
-        log.warn('Possible sign loss when converting '
-                 'negative image of type %s to positive '
-                 'image of type %s.' % (np.dtype(dtype_in), np.dtype(dtype)))
-
-    # If input type is non-negative, or if
-    # converting to a positive-only type, then we
-    # there's no need to shift numbers to the negative side
-    if sign_loss or np.sign(imin) != -1:
-        shift = 0
-        omin = 0
-    else:
-        shift = omin
-
-    scale = (omax - omin) / (imax - imin)
-
-    if dtype in integer_types:
-        round_fn = np.round
-    else:
-        round_fn = lambda x: x
-
-    # Do scaling/shifting calculations in floating point
-    image = image.astype(np.float64)
-    out = image - imin
-    out *= scale
-    out += shift
-    out = round_fn(out).astype(dtype)
-
-    return out
+    if kind_in == 'f':
+        if kind == 'f':
+            # floating point -> floating point
+            if itemsize_in > itemsize:
+                prec_loss()
+            return dtype(image)
+        # floating point -> integer
+        prec_loss()
+        image = np.array(image, dtype=np.promote_types(dtype_in, dtype))
+        image *= np.iinfo(dtype).max + 1
+        np.clip(image, 0, np.iinfo(dtype).max, out=image)
+        return dtype(image)
+    if kind == 'f':
+        # integer -> floating point
+        if itemsize_in >= itemsize:
+            prec_loss()
+        image = np.array(image, dtype=np.promote_types(dtype_in, dtype))
+        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
+        else:
+            result += image
+        return dtype(result)
+    if kind == 'u':
+        # signed integer -> unsigned integer
+        sign_loss()
+        image = np.array(image, dtype=np.promote_types(dtype_in, dtype))
+        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 += image
+        return result
+    if kind == 'i':
+        # signed integer -> signed integer
+        if itemsize_in > itemsize:
+            prec_loss()
+            return dtype(image // 2**(8 * (itemsize_in - itemsize)))
+        # upcast to next higher precision signed integer type
+        dt = next(dt for dt in (np.int16, np.int32, np.int64)
+                  if image.itemsize < np.dtype(dt).itemsize)
+        image = np.array(image, dtype=dt)
+        image -= np.iinfo(dtype_in).min
+        # upcast to next higher precision signed integer type
+        dt = next(dt for dt in (np.int32, np.int64)
+                  if image.itemsize < np.dtype(dt).itemsize)
+        result = np.array(image, dtype=dt)
+        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)
 
 
 def img_as_float(image):
@@ -103,8 +165,7 @@ def img_as_float(image):
     Negative input values will be shifted to the positive domain.
 
     """
-    prec_loss = ()
-    return _convert(image, np.float64, prec_loss)
+    return _convert(image, np.float64)
 
 
 def img_as_uint(image):
@@ -125,9 +186,7 @@ def img_as_uint(image):
     Negative input values will be shifted to the positive domain.
 
     """
-
-    prec_loss = (np.float32, np.float64)
-    return _convert(image, np.uint16, prec_loss)
+    return _convert(image, np.uint16)
 
 
 def img_as_int(image):
@@ -149,8 +208,7 @@ def img_as_int(image):
     the output image will still only have positive values.
 
     """
-    prec_loss = (np.float32, np.float64, np.uint16)
-    return _convert(image, np.int16, prec_loss)
+    return _convert(image, np.int16)
 
 
 def img_as_ubyte(image):
@@ -172,5 +230,4 @@ def img_as_ubyte(image):
     the output image will still only have positive values.
 
     """
-    prec_loss = (np.float32, np.float64, np.uint16, np.int16, np.int8)
-    return _convert(image, np.ubyte, prec_loss)
+    return _convert(image, np.uint8)