Numpy and 16-bit PGM

Question

What is an efficient and clear way to read 16-bit PGM images in Python with numpy?

I cannot use PIL to load 16-bit PGM images due to a PIL bug. I can read in the header with the following code:

dt = np.dtype([('type', 'a2'),
               ('space_0', 'a1', ),
               ('x', 'a3', ),
               ('space_1', 'a1', ),
               ('y', 'a3', ),
               ('space_2', 'a1', ),
               ('maxval', 'a5')])
header = np.fromfile( 'img.pgm', dtype=dt )
print header

This prints the correct data: ('P5', ' ', '640', ' ', '480', ' ', '65535') But I have a feeling that is not quite the best way. And beyond that, I'm having trouble how to figure out how to read in the following data of x by y (in this case 640x480) by 16-bit with the offset of size(header).

EDIT: IMAGE ADDED

MATLAB code to read and display the image is:

I = imread('foo.pgm'); 
imagesc(I);

And looks like this:

Answer 1

import re
import numpy

def read_pgm(filename, byteorder='>'):
    """Return image data from a raw PGM file as numpy array.

    Format specification: http://netpbm.sourceforge.net/doc/pgm.html

    """
    with open(filename, 'rb') as f:
        buffer = f.read()
    try:
        header, width, height, maxval = re.search(
            b"(^P5\s(?:\s*#.*[\r\n])*"
            b"(\d+)\s(?:\s*#.*[\r\n])*"
            b"(\d+)\s(?:\s*#.*[\r\n])*"
            b"(\d+)\s(?:\s*#.*[\r\n]\s)*)", buffer).groups()
    except AttributeError:
        raise ValueError("Not a raw PGM file: '%s'" % filename)
    return numpy.frombuffer(buffer,
                            dtype='u1' if int(maxval) < 256 else byteorder+'u2',
                            count=int(width)*int(height),
                            offset=len(header)
                            ).reshape((int(height), int(width)))


if __name__ == "__main__":
    from matplotlib import pyplot
    image = read_pgm("foo.pgm", byteorder='<')
    pyplot.imshow(image, pyplot.cm.gray)
    pyplot.show()

Answer 2

I'm not terribly familar with the PGM format, but generally speaking you'd just use numpy.fromfile. fromfile will start at whatever position the file pointer you pass to it is at, so you can simply seek (or read) to the end of the header, and then use fromfile to read the rest in.

You'll need to use infile.readline() instead of next(infile).

import numpy as np

with open('foo.pgm', 'r') as infile:
    header = infile.readline()
    width, height, maxval = [int(item) for item in header.split()[1:]]
    image = np.fromfile(infile, dtype=np.uint16).reshape((height, width))

On a side note, the "foo.pgm" file you pointed to in your comment appears to specify the wrong number of rows in the header.

If you're going to be reading in a lot of files that potentially have that problem, you can just pad the array with zeros or truncate it, like this.

import numpy as np

with open('foo.pgm', 'r') as infile:
    header = next(infile)
    width, height, maxval = [int(item) for item in header.split()[1:]]
    image = np.fromfile(infile, dtype=np.uint16)
    if image.size < width * height:
        pad = np.zeros(width * height - image.size, dtype=np.uint16)
        image = np.hstack([image, pad])
    if image.size > width * height:
        image = image[:width * height]
    image = image.reshape((height, width))