Gets difficult to optimize from here. It spends most of the time inside png_filters.c and the builtin zlib decode function.

png_decode.py

import struct
import zlib

png_header = "89504E470D0A1A0A".decode('hex')

def png_chunks(file):
    header = file.read(8)
    assert header == png_header
    length = 1
    while length > 0:
        length, chunk_type = struct.unpack('>I4s', file.read(8))
        data = file.read(length)
        crc  = struct.unpack('>I', file.read(4))[0]
        yield chunk_type, data, crc

def decode_IDAT(z, data, crc):
    ver = zlib.crc32('IDAT')
    ver = zlib.crc32(data, ver)
    ver &= 2**32 - 1
    assert crc == ver
    data = z.decompress(data)
    #hexdump(data[0:8])
    return data
    #data = zlib.decompress(data)

def hexdump(data, n=8):
    for i in range(0, len(data), n):
        print ' '.join(x.encode('hex') for x in data[i:i+n])

def read_unfiltered(chunks):
    unfiltered = ''
    z = zlib.decompressobj()
    for chunk_type, data, crc in chunks:
        if chunk_type == 'IDAT':
            unfiltered += decode_IDAT(z, data, crc)
        elif chunk_type == 'IEND':
            break
        elif chunk_type[0].isupper():
            print "{} {} crc={}".format(chunk_type, len(data), hex(crc))
            if len(data) < 100:
                hexdump(data)
            raise Exception("critical chunk")
    return unfiltered

from PIL import Image
from ctypes import CDLL, c_char, c_void_p, c_int

png_filters = CDLL('./png_filters.so')

png_unfilter = png_filters.png_unfilter
png_unfilter.restype = c_int
png_unfilter.argtypes = [c_void_p, c_void_p, c_int, c_int, c_int]

def unfilter_image(width, height, unfiltered):
    bpp = 3 # bytes per pixel
    pixels = (c_char*(width*height*bpp))()
    if png_unfilter(pixels, unfiltered, width, height, bpp) != 0:
        raise Exception("unknown filtering")
    return pixels

def read_file():
    with open('grumpycat.png', 'rb') as file:
        chunks = png_chunks(file)

        chunk_type, data, crc = chunks.next()
        assert chunk_type == 'IHDR'
        ihdr = struct.unpack('>II5b', data)
        width, height = ihdr[:2]
        bitdepth, colortype = ihdr[2:4]
        assert bitdepth  == 8
        assert colortype == 2
        compression_method, filter_method, interlace_method = ihdr[4:]
        assert compression_method == 0
        assert filter_method      == 0
        assert interlace_method   == 0

        unfiltered = read_unfiltered(chunks)
        assert len(unfiltered) == width*height*3 + height
        pixels = unfilter_image(width, height, unfiltered)
        assert len(pixels) == width*height*3

#    image = Image.fromstring('RGB', (width, height), pixels)
#    image.show()

read_file()
for i in range(100):
    read_file()

png_filters.c

// gcc --std=c99 png_filters.c -o png_filters.so -O2 -shared -fPIC
#include <stdlib.h>
#include <string.h>

typedef unsigned char ubyte;

static inline
int png_filter_pixel(int mode, int a, int b, int c)
{
    int p, pa, pb, pc;
    switch (mode) {
        case 0: return 0;
        case 1: return a;
        case 2: return b;
        case 3: return (a + b) / 2;
        case 4: // paeth
        p = a + b - c;
        pa = abs(p - a);
        pb = abs(p - b);
        pc = abs(p - c);
        if ((pa <= pb) && (pa <= pc))
            return a;
        else if (pb <= pc)
            return b;
        else
            return c;
        default: return 0;
    }
}

static inline
int sample(ubyte* pixels, int stride, int x, int y)
{
    if (x < 0) return 0;
    if (y < 0) return 0;
    return pixels[stride*y+x];
}

int png_unfilter(ubyte* pixels, ubyte* unfiltered, int width, int height, int bpp)
{
    const int stride = width*bpp;
    int scanline;
    int a, b, c;
    for (int y = 0; y < height; y++) {
        scanline = y*stride;

        int    mode  =  unfiltered[y*(stride+1)+0];
        ubyte* input = &unfiltered[y*(stride+1)+1];
        if (mode > 4) return 1;
        for (int x = 0; x < stride; x++) {
            a = sample(pixels, stride, x - bpp, y);
            b = sample(pixels, stride, x,       y - 1);
            c = sample(pixels, stride, x - bpp, y - 1);
            pixels[scanline+x] = input[x] + png_filter_pixel(mode, a, b, c);
        }
    }
    return 0;
}