diff --git a/PIL/ImageMorph.py b/PIL/ImageMorph.py
new file mode 100644
index 000000000..a78e341ed
--- /dev/null
+++ b/PIL/ImageMorph.py
@@ -0,0 +1,234 @@
+# A binary morphology add-on for the Python Imaging Library
+#
+# History:
+#   2014-06-04 Initial version.
+#
+# Copyright (c) 2014 Dov Grobgeld <dov.grobgeld@gmail.com>
+
+from PIL import Image
+from PIL import _imagingmorph
+import re
+
+LUT_SIZE = 1<<9
+class LutBuilder:
+    """A class for building MorphLut's from a descriptive language
+
+      The input patterns is a list of a strings sequences like these:
+
+      4:(...
+         .1.
+         111)->1
+
+      (whitespaces including linebreaks are ignored). The option 4
+      descibes a series of symmetry operations (in this case a
+      4-rotation), the pattern is decribed by:
+
+         . or X - Ignore
+         1 - Pixel is on
+         0 - Pixel is off
+
+      The result of the operation is described after "->" string.
+
+      The default is to return the current pixel value, which is
+      returned if no other match is found.
+
+      Operations:
+         4 - 4 way rotation
+         N - Negate
+         1 - Dummy op for no other operation (an op must always be given)
+         M - Mirroring
+
+      Example:
+
+      lb = LutBuilder(patterns = ["4:(... .1. 111)->1"])
+      lut = lb.build_lut()
+         
+    """
+    def __init__(self,patterns = None,op_name=None):
+        if patterns is not None:
+            self.patterns = patterns
+        else:
+            self.patterns = []
+        self.lut = None
+        if op_name is not None:
+            known_patterns = {
+                'corner' : ['1:(... ... ...)->0',
+                            '4:(00. 01. ...)->1'],
+                'dilation4' : ['4:(... .0. .1.)->1'],
+                'dilation8' : ['4:(... .0. .1.)->1',
+                               '4:(... .0. ..1)->1'],
+                'erosion4' : ['4:(... .1. .0.)->0'],
+                'erosion8' : ['4:(... .1. .0.)->0',
+                              '4:(... .1. ..0)->0'],
+                'edge' : ['1:(... ... ...)->0',
+                          '4:(.0. .1. ...)->1',
+                          '4:(01. .1. ...)->1']
+            }
+            if not op_name in known_patterns:
+                raise Exception('Unknown pattern '+op_name+'!')
+
+            self.patterns = known_patterns[op_name]
+
+    def add_patterns(self, patterns):
+        self.patterns += patterns
+
+    def build_default_lut(self):
+        symbols = ['\0','\1']
+        m = 1 << 4  # pos of current pixel
+        self.lut = bytearray(''.join([symbols[(i & m)>0] for i in range(LUT_SIZE)]))
+
+    def get_lut(self):
+        return self.lut
+
+    def _string_permute(self, pattern, permutation):
+        """string_permute takes a pattern and a permutation and returns the
+        string permuted accordinging to the permutation list.
+        """
+        assert(len(permutation)==9)
+        return ''.join([pattern[p] for p in permutation])
+
+    def _pattern_permute(self, basic_pattern, options, basic_result):
+        """pattern_permute takes a basic pattern and its result and clones
+        the mattern according to the modifications described in the $options
+        parameter. It returns a list of all cloned patterns."""
+        patterns = [(basic_pattern, basic_result)]
+
+        # rotations
+        if '4' in options:
+            res = patterns[-1][1]
+            for i in range(4):
+                patterns.append(
+                    (self._string_permute(patterns[-1][0],
+                                    [6,3,0,
+                                     7,4,1,
+                                     8,5,2]), res))
+        # mirror
+        if 'M' in options:
+            n = len(patterns)
+            for pattern,res in patterns[0:n]:
+                patterns.append(
+                    (self._string_permute(pattern, [2,1,0,
+					      5,4,3,
+					      8,7,6]), res))
+
+        # negate
+        if 'N' in options:
+            n = len(patterns)
+            for pattern,res in patterns[0:n]:
+                # Swap 0 and 1
+                pattern = (pattern
+                           .replace('0','Z')
+                           .replace('1','0')
+                           .replace('Z','1'))
+                res = '%d'%(1-int(res))
+                patterns.append((pattern, res))
+
+        return patterns
+
+    def build_lut(self):
+        """Compile all patterns into a morphology lut.
+
+        TBD :Build based on (file) morphlut:modify_lut
+        """
+        self.build_default_lut()
+        patterns = []
+
+        # Parse and create symmetries of the patterns strings
+        for p in self.patterns:
+            m = re.search(r'(\w*):?\s*\((.+?)\)\s*->\s*(\d)', p.replace('\n',''))
+            if not m:
+                raise Exception('Syntax error in pattern "'+p+'"')
+            options = m.group(1)
+            pattern = m.group(2)
+            result = int(m.group(3))
+
+            # Get rid of spaces
+            pattern= pattern.replace(' ','').replace('\n','')
+
+            patterns += self._pattern_permute(pattern, options, result)
+
+#        # Debugging
+#        for p,r in patterns:
+#            print p,r
+#        print '--'
+
+        # compile the patterns into regular expressions for speed
+        for i in range(len(patterns)):
+            p = patterns[i][0].replace('.','X').replace('X','[01]')
+            p = re.compile(p)
+            patterns[i] = (p, patterns[i][1])
+
+        # Step through table and find patterns that match.
+        # Note that all the patterns are searched. The last one
+        # caught overrides
+        for i in range(LUT_SIZE):
+            # Build the bit pattern
+            bitpattern = bin(i)[2:]
+            bitpattern = ('0'*(9-len(bitpattern)) + bitpattern)[::-1]
+
+            for p,r in patterns:
+                if p.match(bitpattern):
+                    self.lut[i] = ['\0','\1'][r]
+
+        return self.lut
+        
+class MorphOp:
+    """A class for binary morphological operators"""
+
+    def __init__(self,
+                 lut=None,
+                 op_name = None,
+                 patterns = None):
+        """Create a binary morphological operator"""
+        self.lut = lut
+        if op_name is not None:
+            self.lut = LutBuilder(op_name = op_name).build_lut()
+        elif patterns is not None:
+            self.lut = LutBuilder(patterns = patterns).build_lut()
+
+    def apply(self, image):
+        """Run a single morphological operation on an image
+
+        Returns a tuple of the number of changed pixels and the
+        morphed image"""
+        if self.lut is None:
+            raise Exception('No operator loaded')
+        
+        outimage = Image.new(image.mode, image.size, None)
+        count = _imagingmorph.apply(str(self.lut), image.im.id, outimage.im.id)
+        return count, outimage
+    
+    def match(self, image):
+        """Get a list of coordinates matching the morphological operation on an image
+
+        Returns a list of tuples of (x,y) coordinates of all matching pixels."""
+        if self.lut is None:
+            raise Exception('No operator loaded')
+        
+        return  _imagingmorph.match(str(self.lut), image.im.id)
+
+    def get_on_pixels(self, image):
+        """Get a list of all turned on pixels in a binary image
+
+        Returns a list of tuples of (x,y) coordinates of all matching pixels."""
+        
+        return  _imagingmorph.get_on_pixels(image.im.id)
+
+    def load_lut(self, filename):
+        """Load an operator from an mrl file"""
+        self.lut = bytearray(open(filename,'rb').read())
+        if len(self.lut)!= 8192:
+            self.lut = None
+            raise Exception('Wrong size operator file!')
+
+    def save_lut(self, filename):
+        """Load an operator save mrl file"""
+        if self.lut is None:
+            raise Exception('No operator loaded')
+        open(filename,'wb').write(self.lut)
+
+    def set_lut(self, lut):
+        """Set the lut from an external source"""
+        self.lut = lut
+
+    
diff --git a/Tests/test_imagemorph.py b/Tests/test_imagemorph.py
new file mode 100644
index 000000000..7d0ce8796
--- /dev/null
+++ b/Tests/test_imagemorph.py
@@ -0,0 +1,138 @@
+# Test the ImageMorphology functionality
+from tester import *
+
+from PIL import Image
+from PIL import ImageMorph
+
+def img_to_string(im):
+    """Turn a (small) binary image into a string representation"""
+    chars = '.1'
+    width, height = im.size
+    return '\n'.join(
+        [''.join([chars[im.getpixel((c,r))>0] for c in range(width)])
+         for r in range(height)])
+
+def string_to_img(image_string):
+    """Turn a string image representation into a binary image"""
+    rows = [s for s in image_string.replace(' ','').split('\n')
+            if len(s)]
+    height = len(rows)
+    width = len(rows[0])
+    im = Image.new('L',(width,height))
+    for i in range(width):
+        for j in range(height):
+            c = rows[j][i]
+            v = c in 'X1'
+            im.putpixel((i,j),v)
+
+    return im
+
+def img_string_normalize(im):
+  return img_to_string(string_to_img(im))
+
+def assert_img_equal(A,B):
+  assert_equal(img_to_string(A), img_to_string(B))
+
+def assert_img_equal_img_string(A,Bstring):
+  assert_equal(img_to_string(A), img_string_normalize(Bstring))
+
+A = string_to_img(
+"""    
+.......
+.......
+..111..
+..111..
+..111..
+.......
+.......
+"""
+)    
+
+# Test the named patterns
+
+# erosion8
+mop = ImageMorph.MorphOp(op_name='erosion8')
+count,Aout = mop.apply(A)
+assert_equal(count,8)
+assert_img_equal_img_string(Aout,
+"""
+.......
+.......
+.......
+...1...
+.......
+.......
+.......
+""")
+
+# erosion8
+mop = ImageMorph.MorphOp(op_name='dilation8')
+count,Aout = mop.apply(A)
+assert_equal(count,16)
+assert_img_equal_img_string(Aout,
+"""
+.......
+.11111.
+.11111.
+.11111.
+.11111.
+.11111.
+.......
+""")
+
+# erosion4
+mop = ImageMorph.MorphOp(op_name='dilation4')
+count,Aout = mop.apply(A)
+assert_equal(count,12)
+assert_img_equal_img_string(Aout,
+"""
+.......
+..111..
+.11111.
+.11111.
+.11111.
+..111..
+.......
+""")
+
+# edge 
+mop = ImageMorph.MorphOp(op_name='edge')
+count,Aout = mop.apply(A)
+assert_equal(count,1)
+assert_img_equal_img_string(Aout,
+"""
+.......
+.......
+..111..
+..1.1..
+..111..
+.......
+.......
+""")
+
+# Create a corner detector pattern
+mop = ImageMorph.MorphOp(patterns = ['1:(... ... ...)->0',
+                                     '4:(00. 01. ...)->1'])
+count,Aout = mop.apply(A)
+assert_equal(count,5)
+assert_img_equal_img_string(Aout,
+"""
+.......
+.......
+..1.1..
+.......
+..1.1..
+.......
+.......
+""")
+
+# Test the coordinate counting with the same operator
+coords = mop.match(A)
+assert_equal(len(coords), 4)
+assert_equal(tuple(coords),
+             ((2,2),(4,2),(2,4),(4,4)))
+
+coords = mop.get_on_pixels(Aout)
+assert_equal(len(coords), 4)
+assert_equal(tuple(coords),
+             ((2,2),(4,2),(2,4),(4,4)))
diff --git a/_imagingmorph.c b/_imagingmorph.c
new file mode 100644
index 000000000..f80124022
--- /dev/null
+++ b/_imagingmorph.c
@@ -0,0 +1,286 @@
+/*
+ * The Python Imaging Library
+ *
+ * A binary morphology add-on for the Python Imaging Library
+ *
+ * History:
+ *   2014-06-04 Initial version.
+ *
+ * Copyright (c) 2014 Dov Grobgeld <dov.grobgeld@gmail.com>
+ *
+ * See the README file for information on usage and redistribution.
+ */
+
+#include "Python.h"
+#include "Imaging.h"
+#include "py3.h"
+
+#define LUT_SIZE (1<<9)
+
+/* Apply a morphologic LUT to a binary image. Outputs a
+   a new binary image.
+
+   Expected parameters:
+
+      1. a LUT - a 512 byte size lookup table.
+      2. an input Imaging image id.
+      3. an output Imaging image id
+
+   Returns number of changed pixels.
+*/
+static PyObject*
+apply(PyObject *self, PyObject* args)
+{
+    const char *lut;
+    Py_ssize_t lut_len, i0, i1;
+    Imaging imgin, imgout;
+    int width, height;
+    int row_idx, col_idx;
+    UINT8 **inrows, **outrows;
+    int num_changed_pixels = 0;
+
+    if (!PyArg_ParseTuple(args, "s#nn", &lut, &lut_len, &i0, &i1)) {
+        PyErr_SetString(PyExc_RuntimeError, "Argument parsing problem");
+
+        return NULL;
+    }
+
+    if (lut_len < LUT_SIZE) {
+        PyErr_SetString(PyExc_RuntimeError, "The morphology LUT has the wrong size");
+        return NULL;
+    }
+
+    imgin = (Imaging) i0;
+    imgout = (Imaging) i1;
+    width = imgin->xsize;
+    height = imgin->ysize;
+
+    if (imgin->type != IMAGING_TYPE_UINT8 &&
+        imgin->bands != 1) {
+        PyErr_SetString(PyExc_RuntimeError, "Unsupported image type");
+        return NULL;
+    }
+    if (imgout->type != IMAGING_TYPE_UINT8 &&
+        imgout->bands != 1) {
+        PyErr_SetString(PyExc_RuntimeError, "Unsupported image type");
+        return NULL;
+    }
+
+    inrows = imgin->image8;
+    outrows = imgout->image8;
+
+    for (row_idx=0; row_idx < height; row_idx++) {
+        UINT8 *outrow = outrows[row_idx];
+        UINT8 *inrow = inrows[row_idx];
+        UINT8 *prow, *nrow; /* Previous and next row */
+
+        /* zero boundary conditions. TBD support other modes */
+        outrow[0] = outrow[width-1] = 0;
+        if (row_idx==0 || row_idx == height-1) {
+            for(col_idx=0; col_idx<width; col_idx++)
+                outrow[col_idx] = 0;
+            continue;
+        }
+
+        prow = inrows[row_idx-1];
+        nrow = inrows[row_idx+1];
+
+        for (col_idx=1; col_idx<width-1; col_idx++) {
+            int cim = col_idx-1;
+            int cip = col_idx+1;
+            unsigned char b0 = prow[cim] &1;
+            unsigned char b1 = prow[col_idx]&1;
+            unsigned char b2 = prow[cip]&1;
+
+            unsigned char b3 = inrow[cim]&1;
+            unsigned char b4 = inrow[col_idx]&1;
+            unsigned char b5 = inrow[cip]&1;
+
+            unsigned char b6 = nrow[cim]&1;
+            unsigned char b7 = nrow[col_idx]&1;
+            unsigned char b8 = nrow[cip]&1;
+
+            int lut_idx = (b0 
+                           |(b1 << 1)
+                           |(b2 << 2)
+                           |(b3 << 3)
+                           |(b4 << 4)
+                           |(b5 << 5)
+                           |(b6 << 6)
+                           |(b7 << 7)
+                           |(b8 << 8));
+            outrow[col_idx] = 255*(lut[lut_idx]&1);
+            num_changed_pixels += ((b4&1)!=(outrow[col_idx]&1));
+        }
+    }
+    return Py_BuildValue("i",num_changed_pixels);
+}
+
+/* Match a morphologic LUT to a binary image and return a list
+   of the coordinates of all all matching pixels.
+
+   Expected parameters:
+
+      1. a LUT - a 512 byte size lookup table.
+      2. an input Imaging image id.
+
+   Returns list of matching pixels.
+*/
+static PyObject*
+match(PyObject *self, PyObject* args)
+{
+    const char *lut;
+    Py_ssize_t lut_len, i0;
+    Imaging imgin;
+    int width, height;
+    int row_idx, col_idx;
+    UINT8 **inrows;
+    PyObject *ret = PyList_New(0);
+
+    if (!PyArg_ParseTuple(args, "s#n", &lut, &lut_len, &i0)) {
+        PyErr_SetString(PyExc_RuntimeError, "Argument parsing problem");
+
+        return NULL;
+    }
+
+    if (lut_len < LUT_SIZE) {
+        PyErr_SetString(PyExc_RuntimeError, "The morphology LUT has the wrong size");
+        return NULL;
+    }
+
+    imgin = (Imaging) i0;
+
+    if (imgin->type != IMAGING_TYPE_UINT8 &&
+        imgin->bands != 1) {
+        PyErr_SetString(PyExc_RuntimeError, "Unsupported image type");
+        return NULL;
+    }
+
+    inrows = imgin->image8;
+    width = imgin->xsize;
+    height = imgin->ysize;
+
+    for (row_idx=1; row_idx < height-1; row_idx++) {
+        UINT8 *inrow = inrows[row_idx];
+        UINT8 *prow, *nrow;
+
+        prow = inrows[row_idx-1];
+        nrow = inrows[row_idx+1];
+
+        for (col_idx=1; col_idx<width-1; col_idx++) {
+            int cim = col_idx-1;
+            int cip = col_idx+1;
+            unsigned char b0 = prow[cim] &1;
+            unsigned char b1 = prow[col_idx]&1;
+            unsigned char b2 = prow[cip]&1;
+
+            unsigned char b3 = inrow[cim]&1;
+            unsigned char b4 = inrow[col_idx]&1;
+            unsigned char b5 = inrow[cip]&1;
+
+            unsigned char b6 = nrow[cim]&1;
+            unsigned char b7 = nrow[col_idx]&1;
+            unsigned char b8 = nrow[cip]&1;
+
+            int lut_idx = (b0 
+                           |(b1 << 1)
+                           |(b2 << 2)
+                           |(b3 << 3)
+                           |(b4 << 4)
+                           |(b5 << 5)
+                           |(b6 << 6)
+                           |(b7 << 7)
+                           |(b8 << 8));
+            if (lut[lut_idx]) {
+                PyObject *coordObj = Py_BuildValue("(nn)",col_idx,row_idx);
+                PyList_Append(ret, coordObj);
+            }
+        }
+    }
+
+    return ret;
+}
+
+/* Return a list of the coordinates of all turned on pixels in an image.
+   May be used to extract features after a sequence of MorphOp's were applied.
+   This is faster than match as only 1x1 lookup is made.
+*/
+static PyObject*
+get_on_pixels(PyObject *self, PyObject* args)
+{
+    Py_ssize_t i0;
+    Imaging img;
+    UINT8 **rows;
+    int row_idx, col_idx;
+    int width, height;
+    PyObject *ret = PyList_New(0);
+  
+    if (!PyArg_ParseTuple(args, "n", &i0)) {
+        PyErr_SetString(PyExc_RuntimeError, "Argument parsing problem");
+
+        return NULL;
+    }
+    img = (Imaging) i0;
+    rows = img->image8;
+    width = img->xsize;
+    height = img->ysize;
+  
+    for (row_idx=0; row_idx < height; row_idx++) {
+        UINT8 *row = rows[row_idx];
+        for (col_idx=0; col_idx<width; col_idx++) {
+            if (row[col_idx]) {
+                PyObject *coordObj = Py_BuildValue("(nn)",col_idx,row_idx);
+                PyList_Append(ret, coordObj);
+            }
+        }
+    }
+    return ret;
+}
+
+
+static int
+setup_module(PyObject* m)
+{
+    PyObject* d = PyModule_GetDict(m);
+
+    PyDict_SetItemString(d, "__version", PyUnicode_FromString("0.1"));
+
+    return 0;
+}
+
+static PyMethodDef functions[] = {
+    /* Functions */
+    {"apply", (PyCFunction)apply, 1},
+    {"get_on_pixels", (PyCFunction)get_on_pixels, 1},
+    {"match", (PyCFunction)match, 1},
+};
+
+#if PY_VERSION_HEX >= 0x03000000
+PyMODINIT_FUNC
+PyInit__imagingmorph(void) {
+    PyObject* m;
+
+    static PyModuleDef module_def = {
+        PyModuleDef_HEAD_INIT,
+        "_imagingmorph",         /* m_name */
+        "A module for doing image morphology",               /* m_doc */
+        -1,                 /* m_size */
+        functions,          /* m_methods */
+    };
+
+    m = PyModule_Create(&module_def);
+
+    if (setup_module(m) < 0)
+        return NULL;
+
+    return m;
+}
+#else
+PyMODINIT_FUNC
+init_imagingmorph(void)
+{
+    PyObject* m = Py_InitModule("_imagingmorph", functions);
+    setup_module(m);
+}
+#endif
+
diff --git a/setup.py b/setup.py
index 50ca985e3..cee913637 100644
--- a/setup.py
+++ b/setup.py
@@ -543,6 +543,9 @@ class pil_build_ext(build_ext):
         if os.path.isfile("_imagingmath.c"):
             exts.append(Extension("PIL._imagingmath", ["_imagingmath.c"]))
 
+        if os.path.isfile("_imagingmorph.c"):
+            exts.append(Extension("PIL._imagingmorph", ["_imagingmorph.c"]))
+
         self.extensions[:] = exts
 
         build_ext.build_extensions(self)