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:

Eyebright

can you attach an example img.pgm? Off topic: checked your site; you might like to see this: seems you are not the only ones looking for warmer water around the arctic... (supporting evidence for your (coleages) thesis perhaps?) – Instinctive

PGM here: db.tt/phaR587 P.S. One doesn't have to look very hard to find these things... :(. – Eyebright

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()

Georginageorgine

Very nice, but in the case of this test file, >u2 produces the wrong values (range 4098 to 65287) while u2 produces the correct values (528 to 2047). You mention big-endian in another comment. The data was produced on and I am reading on an Intel (little endian) chip. I presume it was written in native format. – Eyebright

The specification says "The most significant byte is first.", which is big endian. See also en.wikipedia.org/wiki/Netpbm_format#16-bit_extensions. – Georginageorgine

Matlab reads the data as big endian, so the image shown in your question would be wrong(?). You can always swap the bytes later if you are reading a non-standard file. – Georginageorgine

Added a byteorder parameter to the function. – Georginageorgine

Yup you are right. The data is written in little-endian, contrary to the spec. My image looks correct but with a colorbar it is clearly scaled wrong. – Eyebright

Thanks for this valuable script. I have a bug to report though: I am using the script to parse a binary 16 bit pgm whose data section starts with the value 2851 (alternatively an 8 bit pgm that starts with 11 35). This translates to the hex values 0B 23 which as characters are a vertical tab (interpreted as whitespace by the regex) and a # (interpreted as the beginning of a comment). This causes a crash because the first part of the data section (up until the next 0A or 0D) is interpreted as part of the header which in turn causes a buffer overflow due to the offset being too great. – Waterlogged

This string will break your regex: "P5\n120 120\n2866\n #\n " where the header ends after "2866\n". Concretly this occurs in a binary pgm if the data section starts with 0x20 0x23 0x0a 0x20 or in other words if the pgm image starts with the grey values 32 35 10 32. As I see this, errors of this type cannot be avoided since you can represent any comment by a carefully selected byte string. As I see it the only way to make this work would be to forbid multi line comments in the last header line and adapting the regex accordingly. – Waterlogged

I just played around some more with the problem but couldn't find a good solution. You could try and request that the last comment cannot have more than one line, but that didn't work for me either. I think the most robust would be to disallow comments after the last header value (maxval). Then you could use this regex (excuse the horrible formatting): 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)", buffer).groups() – Waterlogged

Unfortunately the standard over at netpbm.sourceforge.net is not very helpful either: "Note also that this means if you have a comment right before the raster, the newline at the end of the comment is not sufficient to delimit the raster." That does not offer a solution, does it? – Waterlogged

In contrast to the netpbm web pages, the netpbm man pages disallow comments after the last header value: Characters from a "#" to the next end-of-line, before the maxval line, are comments and are ignored. – Georginageorgine

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))

Merca

quite elegant, and works for mankoffs binary! I only got weird output when testing it on a standard string formatted pgm file... – Instinctive

@Instinctive - Yeah, I didn't intend it for an ascii pgm file. However, it's quite simple to use np.loadtxt or something similar in that case. – Merca

Close but still a bug. The file is 614417 bytes long which is equal to 640*480*2 + 17, which is a 17 byte header and a 640x480 two-byte (16 bit) data. The image displays properly decoded as such manually in other languages (IDL), and using built-in routines elsewhere (GIMP, MATLAB). I will post a version of the image in the question shortly. Sorry for not providing all this info initially, I'm figuring it out too as I go... – Eyebright

OK I got it. Change next(infile) to infile.read(17). But what if I don't want to hard-code this? It will work, all my PGMs are the same, but it would be nice to get it right. Anyway, thank you for the solution so far. – Eyebright

More details: f = open('foo.pgm'); h=next(f); print f.tell() prints 8192, while f = open('foo.pgm'); h=f.read(17); print f.tell() prints 17. – Eyebright

Just solved this part in an edit to my first answer: open as text just to know the width and height and calculate i_header_end: filesize-(width*height*2); close and open as binary and read from i_header_end onworth. (a bit odd but it seems to work) – Instinctive

@Eyebright - I wasn't thinking and used next(infile) instead of infile.readline(). There's no need to manually calculate the size of the header. – Merca

Answer 3

from here I understand that the header information can be separated by either spaces, carriage returns or others. If yours is separated by spaces (inform me if otherwise) you can do:

with open('img.pgm') as f:
    lines = f.readlines()
    data = np.array([line.split() for line in lines[1:]], dtype=np.int16).T

your data is now an array in int16 format!

Suppose you are still interested in the header information, you can do:

class Header(object):
    def __init__(self, type, width, height, maxval):
        self.type = type
        self.width = int(width)
        self.height = int(height)
        self.maxval = int(maxval)

h = Header(*lines[0].split()[:4])

so that you can check the image data against the read lines:

assert (h.width, h.height) == data.shape    
assert h.maxval >= data.max()

Edit: with the image data being binary, the file has to be opened as 'rb' and read from after the header information:

import numpy as np

def as_array(filepath):
    f = open(filepath, 'r')
    w, h = size = tuple(int(v) for v in next(f).split()[1:3])
    data_size = w * h * 2

    f.seek(0, 2)
    filesize = f.tell()
    f.close()
    i_header_end = filesize - (data_size)

    f = open(filepath, 'rb')
    f.seek(i_header_end)
    buffer = f.read()
    f.close()

    # convert binary data to an array of the right shape
    data = np.frombuffer(buffer, dtype=np.uint16).reshape((w, h))

    return data

a = as_array('foo.pgm')

Instinctive

I think the link you attached describes my format correctly. However, I have the P5 "raw" format (the more common one, described first). The header is ASCII, but the data below is binary, and it seems that readlines() is failing because of this. – Eyebright

Right. readlines() reads one line, but the interpretion of that line has to be via np.fromstring(), or, like you and Joe Kington propose, directly with np.fromfile() since you know it is binary anyway. There is another problem however: see my second answer – Instinctive

Answer 4

Indeed, the 'string' after the header is a binary in your file. I solved that below (found the following: ndarray: [2047 2047 2047 ..., 540 539 539]) but there is another problem: the file is not long enough; counts only 289872 numbers instead of 640*480...

I am terribly sorry for my exageration by making a class for it...

import numpy as np
import Image

class PGM(object):
    def __init__(self, filepath):

        with open(filepath) as f:

            # suppose all header info in first line:
            info = f.readline().split()
            self.type = info[0]
            self.width, self.height, self.maxval = [int(v) for v in info[1:]]
            size = self.width * self.height

            lines = f.readlines()
            dt = [np.int8, np.int16][self.maxval > 255]
            try:
                # this will work if lines are integers separated by e.g. spaces
                self.data = np.array([l.split() for l in lines], dtype=dt).T
            except ValueError:
                # data is binary
                data = np.fromstring(lines[0], dtype=dt)
                if data.size < size:
                    # this is the case for the 'db.tt/phaR587 (foo.pgm)'
                    #raise ValueError('data binary string probably uncomplete')
                    data = np.hstack((data, np.zeros(size-data.size)))
                self.data = data[:size].reshape((self.width, self.height))

            assert (self.width, self.height) == self.data.shape
            assert self.maxval >= self.data.max()

        self._img = None

    def get_img(self):
        if self._img is None:
            # only executed once
            size = (self.width, self.height)
            mode = 'L'
            data = self.data
            self.img = Image.frombuffer(mode, size, data)

        return self.img

    Image = property(get_img)

mypgm = PGM('foo.pgm')

mypgm.Image

edit: great Idea from Joe Kington to fill image with zeros!

Instinctive

File is long enough. I think readline() is reading too much. Perhaps some of the binary is on the first line too? – Eyebright

Answer 5

Thanks to the answer by @joe-kington for helping figure this out. The solution follows.

There is a little bit of extra work to not hard-code the known header length (17 bytes in this case), but to determine it from the header. The PGM standard says that the header usually ends with a newline but can end with any whitespace. I think this code will break on a PGM that uses non-newline whitespace for the end-of-header delimeter. Header size in this case would be determined by the size of variables holding width, height, and maxsize, plus two bytes for 'P5', plus 4 bytes of whitespace.

Other cases where this might break are if the width or height are larger than an int (very big image). Or if the PGM is 8-bit rather than 16-bit (which can be determined from maxval, and possible width, height, and the filesize).

#!/usr/bin/python
import numpy as np
import matplotlib.pyplot as plt

file='foo.pgm'
infile = open(file,'r')
header = next(infile)
width, height, maxval = [int(item) for item in header.split()[1:]]
infile.seek(len(header))
image = np.fromfile(infile, dtype=np.uint16).reshape((height, width))
print width, height, maxval
plt.figimage(image)

Eyebright

The dtype should be big endian. – Georginageorgine