Save numpy array as image with high precision (16 bits) with scikit-image

Question

I am working with 2D floating-point numpy arrays that I would like to save to greyscale .png files with high precision (e.g. 16 bits). I would like to do this using the scikit-image skimage.io package if possible.

Here's the main thing I've tried:

import numpy as np
from skimage import io, exposure, img_as_uint, img_as_float

im = np.array([[1., 2.], [3., 4.]], dtype='float64')
im = exposure.rescale_intensity(im, out_range='float')
im = img_as_uint(im)
im

produces:

array([[    0, 21845],
       [43690, 65535]], dtype=uint16)

First I tried saving this as an image then reloading using the Python Imaging Library:

# try with pil:
io.use_plugin('pil')
io.imsave('test_16bit.png', im)
im2 = io.imread('test_16bit.png')
im2

produces:

array([[  0,  85],
       [170, 255]], dtype=uint8)

So somewhere (in either the write or read) I have lost precision. I then tried with the matplotlib plugin:

# try with matplotlib:
io.use_plugin('matplotlib')
io.imsave('test_16bit.png', im)
im3 = io.imread('test_16bit.png')
im3

gives me a 32-bit float:

array([[ 0.        ,  0.33333334],
       [ 0.66666669,  1.        ]], dtype=float32)

but I doubt this is really 32-bits given that I saved a 16-bit uint to the file. It would be great if someone could point me to where I'm going wrong. I would like this to extend to 3D arrays too (i.e. saving 16 bits per colour channel, for 48 bits per image).

UPDATE:

The problem is with imsave. The images are 8 bits per channel. How can one use io.imsave to output a high bit-depth image?

Answer 1

You wanna use the freeimage library to do so:

import numpy as np
from skimage import io, exposure, img_as_uint, img_as_float

io.use_plugin('freeimage')

im = np.array([[1., 2.], [3., 4.]], dtype='float64')
im = exposure.rescale_intensity(im, out_range='float')
im = img_as_uint(im)

io.imsave('test_16bit.png', im)
im2 = io.imread('test_16bit.png')

Result:

[[    0 21845]
 [43690 65535]]

As for 3D arrays, you need to construct the array properly and then it'll work:

# im = np.array([[1, 2.], [3., 4.]], dtype='float64')
im = np.linspace(0, 1., 300).reshape(10, 10, 3)
im = exposure.rescale_intensity(im, out_range='float')
im = img_as_uint(im)

io.imsave('test_16bit.png', im)
im2 = io.imread('test_16bit.png')

Note that the read image is flipped, so something like np.fliplr(np.flipud(im2)) will bring it to original shape.