Search code examples
pythonpython-imaging-library16-bitpgm

Python and 16-bit PGM

I have 16-bit PGM images that I am trying to read in Python. It seems (?) like PIL does not support this format?

import Image
im = Image.open('test.pgm')
im.show()

Shows roughly the image, but it isn't right. There are dark bands throughout and img is reported to have mode=L. I think this is related to an early question I had about 16-bit TIFF files. Is 16-bit that rare that PIL just does not support it? Any advice how I can read 16-bit PGM files in Python, using PIL or another standard library, or home-grown code?

Solution

  • The following depends only on numpy to load the image, which can be 8-bit or 16-bit raw PGM/PPM. I also show a couple different ways to view the image. The one that uses PIL (import Image) requires that the data first be converted to 8-bit.

    #!/usr/bin/python2 -u

from __future__ import print_function
import sys, numpy

def read_pnm_from_stream( fd ):
   pnm = type('pnm',(object,),{}) ## create an empty container
   pnm.header = fd.readline()
   pnm.magic = pnm.header.split()[0]
   pnm.maxsample = 1 if ( pnm.magic == 'P4' ) else 0
   while ( len(pnm.header.split()) < 3+(1,0)[pnm.maxsample] ): s = fd.readline() ; pnm.header += s if ( len(s) and s[0] != '#' ) else ''
   pnm.width, pnm.height = [int(item) for item in pnm.header.split()[1:3]]
   pnm.samples = 3 if ( pnm.magic == 'P6' ) else 1
   if ( pnm.maxsample == 0 ): pnm.maxsample = int(pnm.header.split()[3])
   pnm.pixels = numpy.fromfile( fd, count=pnm.width*pnm.height*pnm.samples, dtype='u1' if pnm.maxsample < 256 else '>u2' )
   pnm.pixels = pnm.pixels.reshape(pnm.height,pnm.width) if pnm.samples==1 else pnm.pixels.reshape(pnm.height,pnm.width,pnm.samples)
   return pnm

if __name__ == '__main__':

## read image
 # src = read_pnm_from_stream( open(filename) )
   src = read_pnm_from_stream( sys.stdin )
 # print("src.header="+src.header.strip(), file=sys.stderr )
 # print("src.pixels="+repr(src.pixels), file=sys.stderr )

## write image
   dst=src
   dst.pixels = numpy.array([ dst.maxsample-i for i in src.pixels ],dtype=dst.pixels.dtype) ## example image processing
 # print("dst shape: "+str(dst.pixels.shape), file=sys.stderr )
   sys.stdout.write(("P5" if dst.samples==1 else "P6")+"\n"+str(dst.width)+" "+str(dst.height)+"\n"+str(dst.maxsample)+"\n");
   dst.pixels.tofile( sys.stdout ) ## seems to work, I'm not sure how it decides about endianness

## view using Image
   import Image
   viewable = dst.pixels if dst.pixels.dtype == numpy.dtype('u1') else numpy.array([ x>>8 for x in dst.pixels],dtype='u1')
   Image.fromarray(viewable).show()

## view using scipy
   import scipy.misc
   scipy.misc.toimage(dst.pixels).show()

    Usage notes

    • I eventually figured out "how it decides about endianness" -- it is actually storing the image in memory as big-endian (rather than native). This scheme might slow down any non-trivial image processing -- although other performance issues with Python may relegate this concern to a triviality (see below).

    • I asked a question related to the endianness concern here. I also ran into some interesting confusion related to endianness with this because I was testing by preprocessing the image with pnmdepth 65535 which is not good (by itself) for testing endianness since the low and high bytes might end up being the same (I didn't notice right away because print(array) outputs decimal). I should have also applied pnmgamma to save myself some confusion.

    • Because Python is so slow, numpy tries to be sneakyclever about how it applies certain operations (see broadcasting). The first rule of thumb for efficiency with numpy is let numpy handle iteration for you (or put another way don't write your own for loops). The funny thing in the code above is that it only partially follows this rule when doing the "example image processing", and therefore the performance of that line has an extreme dependency on the parameters that were given to reshape.

    • The next big numpy endianness mystery: Why does newbyteorder() seem to return an array, when it's documented to return a dtype. This is relevant if you want to convert to native endian with dst.pixels=dst.pixels.byteswap(True).newbyteorder().

    • Hints on porting to Python 3: binary input with an ASCII text header, read from stdin