diff --git a/PIL/Image.py b/PIL/Image.py
index 03f3973ee..7e02c49bb 100644
--- a/PIL/Image.py
+++ b/PIL/Image.py
@@ -1555,7 +1555,7 @@ class Image(object):
 
         return self._new(self.im.resize(size, resample))
 
-    def rotate(self, angle, resample=NEAREST, expand=0):
+    def rotate(self, angle, resample=NEAREST, expand=0, center=None, translate=None):
         """
         Returns a rotated copy of this image.  This method returns a
         copy of this image, rotated the given number of degrees counter
@@ -1569,10 +1569,14 @@ class Image(object):
            (cubic spline interpolation in a 4x4 environment).
            If omitted, or if the image has mode "1" or "P", it is
            set :py:attr:`PIL.Image.NEAREST`. See :ref:`concept-filters`.
+        :param center: Optional center of rotation (a 2-tuple).  Origin is
+           the upper left corner.  Default is the center of the image.
+        :param translate: An optional final translation.
         :param expand: Optional expansion flag.  If true, expands the output
            image to make it large enough to hold the entire rotated image.
            If false or omitted, make the output image the same size as the
-           input image.
+           input image.  Note that the expand flag assumes rotation around
+           the center and no translation.
         :returns: An :py:class:`~PIL.Image.Image` object.
         """
 
@@ -1588,32 +1592,59 @@ class Image(object):
         if angle == 270 and expand:
             return self.transpose(ROTATE_270)
 
+        # Calculate the affine matrix.  Note that this is the reverse
+        # transformation (from destination image to source) because we
+        # want to interpolate the (discrete) destination pixel from
+        # the local area around the (floating) source pixel.
+
+        # The matrix we actually want (note that it operates from the right):
+        # (1, 0, tx)   (1, 0, cx)   ( cos a, sin a, 0)   (1, 0, -cx)
+        # (0, 1, ty) * (0, 1, cy) * (-sin a, cos a, 0) * (0, 1, -cy)
+        # (0, 0,  1)   (0, 0,  1)   (     0,     0, 1)   (0, 0,   1)
+
+        # The reverse matrix is thus:
+        # (1, 0, cx)   ( cos -a, sin -a, 0)   (1, 0, -cx)   (1, 0, -tx)
+        # (0, 1, cy) * (-sin -a, cos -a, 0) * (0, 1, -cy) * (0, 1, -ty)
+        # (0, 0,  1)   (      0,      0, 1)   (0, 0,   1)   (0, 0,   1)
+
+        # In any case, the final translation may be updated at the end to
+        # compensate for the expand flag.
+
+        w, h = self.size
+
+        if translate is None:
+            translate = [0, 0]
+        if center is None:
+            center = [w / 2.0, h / 2.0]
+
         angle = - math.radians(angle)
         matrix = [
             round(math.cos(angle), 15), round(math.sin(angle), 15), 0.0,
             round(-math.sin(angle), 15), round(math.cos(angle), 15), 0.0
-            ]
-
-        def transform(x, y, matrix=matrix):
+        ]
+        def transform(x, y, matrix):
             (a, b, c, d, e, f) = matrix
             return a*x + b*y + c, d*x + e*y + f
+        matrix[2], matrix[5] = transform(-center[0] - translate[0], -center[1] - translate[1], matrix)
+        matrix[2] += center[0]
+        matrix[5] += center[1]
 
-        w, h = self.size
         if expand:
             # calculate output size
             xx = []
             yy = []
             for x, y in ((0, 0), (w, 0), (w, h), (0, h)):
-                x, y = transform(x, y)
+                x, y = transform(x, y, matrix)
                 xx.append(x)
                 yy.append(y)
-            w = int(math.ceil(max(xx)) - math.floor(min(xx)))
-            h = int(math.ceil(max(yy)) - math.floor(min(yy)))
+            nw = int(math.ceil(max(xx)) - math.floor(min(xx)))
+            nh = int(math.ceil(max(yy)) - math.floor(min(yy)))
 
-        # adjust center
-        x, y = transform(w / 2.0, h / 2.0)
-        matrix[2] = self.size[0] / 2.0 - x
-        matrix[5] = self.size[1] / 2.0 - y
+            # We multiply a translation matrix from the right.  Because of its
+            # special form, this is the same as taking the image of the translation vector
+            # as new translation vector.
+            matrix[2], matrix[5] = transform(-(nw - w) / 2.0, -(nh - h) / 2.0, matrix)
+            w, h = nw, nh
 
         return self.transform((w, h), AFFINE, matrix, resample)
 
diff --git a/Tests/test_image_rotate.py b/Tests/test_image_rotate.py
index dfe5a0731..f2a3c4fa8 100644
--- a/Tests/test_image_rotate.py
+++ b/Tests/test_image_rotate.py
@@ -5,11 +5,11 @@ from PIL import Image
 class TestImageRotate(PillowTestCase):
 
     def test_rotate(self):
-        def rotate(im, mode, angle):
-            out = im.rotate(angle)
+        def rotate(im, mode, angle, center=None, translate=None):
+            out = im.rotate(angle, center=center, translate=translate)
             self.assertEqual(out.mode, mode)
             self.assertEqual(out.size, im.size)  # default rotate clips output
-            out = im.rotate(angle, expand=1)
+            out = im.rotate(angle, center=center, translate=translate, expand=1)
             self.assertEqual(out.mode, mode)
             if angle % 180 == 0:
                 self.assertEqual(out.size, im.size)
@@ -31,6 +31,10 @@ class TestImageRotate(PillowTestCase):
             im = Image.new('RGB',(0,0))
             rotate(im, im.mode, angle)
 
+        rotate(im, im.mode, 45, center=(0, 0))
+        rotate(im, im.mode, 45, translate=(im.size[0]/2, 0))
+        rotate(im, im.mode, 45, center=(0, 0), translate=(im.size[0]/2, 0))
+
 
 if __name__ == '__main__':
     unittest.main()