Merge pull request #2254 from uploadcare/resample-roi

Region of interest (box) for resampling

PIL/Image.py

@@ -1674,7 +1674,7 @@ class Image(object):
 
         return m_im
 
-    def resize(self, size, resample=NEAREST):
+    def resize(self, size, resample=NEAREST, box=None):
         """
         Returns a resized copy of this image.
 
@@ -1687,6 +1687,10 @@ class Image(object):
            If omitted, or if the image has mode "1" or "P", it is
            set :py:attr:`PIL.Image.NEAREST`.
            See: :ref:`concept-filters`.
+        :param box: An optional 4-tuple of floats giving the region
+           of the source image which should be scaled.
+           The values should be within (0, 0, width, height) rectangle.
+           If omitted or None, the entire source is used.
         :returns: An :py:class:`~PIL.Image.Image` object.
         """
 
@@ -1695,22 +1699,28 @@ class Image(object):
         ):
             raise ValueError("unknown resampling filter")
 
-        self.load()
-
         size = tuple(size)
-        if self.size == size:
+
+        if box is None:
+            box = (0, 0) + self.size
+        else:
+            box = tuple(box)
+
+        if self.size == size and box == (0, 0) + self.size:
             return self.copy()
 
         if self.mode in ("1", "P"):
             resample = NEAREST
 
         if self.mode == 'LA':
-            return self.convert('La').resize(size, resample).convert('LA')
+            return self.convert('La').resize(size, resample, box).convert('LA')
 
         if self.mode == 'RGBA':
-            return self.convert('RGBa').resize(size, resample).convert('RGBA')
+            return self.convert('RGBa').resize(size, resample, box).convert('RGBA')
 
-        return self._new(self.im.resize(size, resample))
+        self.load()
+
+        return self._new(self.im.resize(size, resample, box))
 
     def rotate(self, angle, resample=NEAREST, expand=0, center=None,
                translate=None):

Tests/test_image_resample.py

@@ -1,4 +1,7 @@
-from __future__ import print_function
+from __future__ import division, print_function
+
+from contextlib import contextmanager
+
 from helper import unittest, PillowTestCase, hopper
 from PIL import Image, ImageDraw
 
@@ -300,24 +303,46 @@ class CoreResampleAlphaCorrectTest(PillowTestCase):
 
 
 class CoreResamplePassesTest(PillowTestCase):
+    @contextmanager
+    def count(self, diff):
+        count = Image.core.getcount()
+        yield
+        self.assertEqual(Image.core.getcount() - count, diff)
+
     def test_horizontal(self):
         im = hopper('L')
-        count = Image.core.getcount()
+        with self.count(1):
-        im.resize((im.size[0] + 10, im.size[1]), Image.BILINEAR)
+            im.resize((im.size[0] - 10, im.size[1]), Image.BILINEAR)
-        self.assertEqual(Image.core.getcount(), count + 1)
 
     def test_vertical(self):
         im = hopper('L')
-        count = Image.core.getcount()
+        with self.count(1):
-        im.resize((im.size[0], im.size[1] + 10), Image.BILINEAR)
+            im.resize((im.size[0], im.size[1] - 10), Image.BILINEAR)
-        self.assertEqual(Image.core.getcount(), count + 1)
 
     def test_both(self):
         im = hopper('L')
-        count = Image.core.getcount()
+        with self.count(2):
-        im.resize((im.size[0] + 10, im.size[1] + 10), Image.BILINEAR)
+            im.resize((im.size[0] - 10, im.size[1] - 10), Image.BILINEAR)
-        self.assertEqual(Image.core.getcount(), count + 2)
 
+    def test_box_horizontal(self):
+        im = hopper('L')
+        box = (20, 0, im.size[0] - 20, im.size[1])
+        with self.count(1):
+            # the same size, but different box
+            with_box = im.resize(im.size, Image.BILINEAR, box)
+        with self.count(2):
+            cropped = im.crop(box).resize(im.size, Image.BILINEAR)
+        self.assert_image_similar(with_box, cropped, 0.1)
+
+    def test_box_vertical(self):
+        im = hopper('L')
+        box = (0, 20, im.size[0], im.size[1] - 20)
+        with self.count(1):
+            # the same size, but different box
+            with_box = im.resize(im.size, Image.BILINEAR, box)
+        with self.count(2):
+            cropped = im.crop(box).resize(im.size, Image.BILINEAR)
+        self.assert_image_similar(with_box, cropped, 0.1)
 
 class CoreResampleCoefficientsTest(PillowTestCase):
     def test_reduce(self):
@@ -347,5 +372,92 @@ class CoreResampleCoefficientsTest(PillowTestCase):
         self.assertEqual(histogram[0x100 * 3 + 0xff], 0x10000)  # fourth channel
 
 
+class CoreResampleBoxTest(PillowTestCase):
+    def test_wrong_arguments(self):
+        im = hopper()
+        for resample in (Image.NEAREST, Image.BOX, Image.BILINEAR, Image.HAMMING,
+                Image.BICUBIC, Image.LANCZOS):
+            im.resize((32, 32), resample, (0, 0, im.width, im.height))
+            im.resize((32, 32), resample, (20, 20, im.width, im.height))
+            im.resize((32, 32), resample, (20, 20, 20, 100))
+            im.resize((32, 32), resample, (20, 20, 100, 20))
+
+            with self.assertRaisesRegexp(TypeError, "must be sequence of length 4"):
+                im.resize((32, 32), resample, (im.width, im.height))
+
+            with self.assertRaisesRegexp(ValueError, "can't be negative"):
+                im.resize((32, 32), resample, (-20, 20, 100, 100))
+            with self.assertRaisesRegexp(ValueError, "can't be negative"):
+                im.resize((32, 32), resample, (20, -20, 100, 100))
+
+            with self.assertRaisesRegexp(ValueError, "can't be empty"):
+                im.resize((32, 32), resample, (20.1, 20, 20, 100))
+            with self.assertRaisesRegexp(ValueError, "can't be empty"):
+                im.resize((32, 32), resample, (20, 20.1, 100, 20))
+            with self.assertRaisesRegexp(ValueError, "can't be empty"):
+                im.resize((32, 32), resample, (20.1, 20.1, 20, 20))
+
+            with self.assertRaisesRegexp(ValueError, "can't exceed"):
+                im.resize((32, 32), resample, (0, 0, im.width + 1, im.height))
+            with self.assertRaisesRegexp(ValueError, "can't exceed"):
+                im.resize((32, 32), resample, (0, 0, im.width, im.height + 1))
+
+    def resize_tiled(self, im, dst_size, xtiles, ytiles):
+        def split_range(size, tiles):
+            scale = size / tiles
+            for i in range(tiles):
+                yield (int(round(scale * i)), int(round(scale * (i + 1))))
+
+        tiled = Image.new(im.mode, dst_size)
+        scale = (im.size[0] / tiled.size[0], im.size[1] / tiled.size[1])
+
+        for y0, y1 in split_range(dst_size[1], ytiles):
+            for x0, x1 in split_range(dst_size[0], xtiles):
+                box = (x0 * scale[0], y0 * scale[1],
+                       x1 * scale[0], y1 * scale[1])
+                tile = im.resize((x1 - x0, y1 - y0), Image.BICUBIC, box)
+                tiled.paste(tile, (x0, y0))
+        return tiled
+
+    def test_tiles(self):
+        im = Image.open("Tests/images/flower.jpg")
+        assert im.size == (480, 360)
+        dst_size = (251, 188)
+        reference = im.resize(dst_size, Image.BICUBIC)
+
+        for tiles in [(1, 1), (3, 3), (9, 7), (100, 100)]:
+            tiled = self.resize_tiled(im, dst_size, *tiles)
+            self.assert_image_similar(reference, tiled, 0.01)
+
+    def test_subsample(self):
+        # This test shows advantages of the subpixel resizing
+        # after supersampling (e.g. during JPEG decoding).
+        im = Image.open("Tests/images/flower.jpg")
+        assert im.size == (480, 360)
+        dst_size = (48, 36)
+        # Reference is cropped image resized to destination
+        reference = im.crop((0, 0, 473, 353)).resize(dst_size, Image.BICUBIC)
+        # Image.BOX emulates supersampling (480 / 8 = 60, 360 / 8 = 45)
+        supersampled = im.resize((60, 45), Image.BOX)
+
+        with_box = supersampled.resize(dst_size, Image.BICUBIC,
+                                       (0, 0, 59.125, 44.125))
+        without_box = supersampled.resize(dst_size, Image.BICUBIC)
+
+        # error with box should be much smaller than without
+        self.assert_image_similar(reference, with_box, 6)
+        with self.assertRaisesRegexp(AssertionError, "difference 29\."):
+            self.assert_image_similar(reference, without_box, 5)
+
+    def test_formats(self):
+        for resample in [Image.NEAREST, Image.BILINEAR]:
+            for mode in ['RGB', 'L', 'RGBA', 'LA', 'I', '']:
+                im = hopper(mode)
+                box = (20, 20, im.size[0] - 20, im.size[1] - 20)
+                with_box = im.resize((32, 32), resample, box)
+                cropped = im.crop(box).resize((32, 32), resample)
+                self.assert_image_similar(cropped, with_box, 0.4)
+
+
 if __name__ == '__main__':
     unittest.main()

_imaging.c

@@ -1503,24 +1503,43 @@ _resize(ImagingObject* self, PyObject* args)
 
     int xsize, ysize;
     int filter = IMAGING_TRANSFORM_NEAREST;
-    if (!PyArg_ParseTuple(args, "(ii)|i", &xsize, &ysize, &filter))
+    float box[4] = {0, 0, 0, 0};
-        return NULL;
     
     imIn = self->image;
+    box[2] = imIn->xsize;
+    box[3] = imIn->ysize;
+    
+    if (!PyArg_ParseTuple(args, "(ii)|i(ffff)", &xsize, &ysize, &filter,
+                          &box[0], &box[1], &box[2], &box[3]))
+        return NULL;
 
     if (xsize < 1 || ysize < 1) {
         return ImagingError_ValueError("height and width must be > 0");
     }
 
-    if (imIn->xsize == xsize && imIn->ysize == ysize) {
+    if (box[0] < 0 || box[1] < 0) {
+        return ImagingError_ValueError("box offset can't be negative");
+    }
+
+    if (box[2] > imIn->xsize || box[3] > imIn->ysize) {
+        return ImagingError_ValueError("box can't exceed original image size");
+    }
+
+    if (box[2] - box[0] < 0 || box[3] - box[1] < 0) {
+        return ImagingError_ValueError("box can't be empty");
+    }
+
+    if (box[0] == 0 && box[1] == 0 && box[2] == xsize && box[3] == ysize) {
         imOut = ImagingCopy(imIn);
     }
     else if (filter == IMAGING_TRANSFORM_NEAREST) {
         double a[6];
 
         memset(a, 0, sizeof a);
-        a[0] = (double) imIn->xsize / xsize;
+        a[0] = (double) (box[2] - box[0]) / xsize;
-        a[4] = (double) imIn->ysize / ysize;
+        a[4] = (double) (box[3] - box[1]) / ysize;
+        a[2] = box[0];
+        a[5] = box[1];
 
         imOut = ImagingNewDirty(imIn->mode, xsize, ysize);
 
@@ -1530,7 +1549,7 @@ _resize(ImagingObject* self, PyObject* args)
             a, filter, 1);
     }
     else {
-        imOut = ImagingResample(imIn, xsize, ysize, filter);
+        imOut = ImagingResample(imIn, xsize, ysize, filter, box);
     }
 
     return PyImagingNew(imOut);

libImaging/Imaging.h

@@ -291,7 +291,7 @@ extern Imaging ImagingRankFilter(Imaging im, int size, int rank);
 extern Imaging ImagingRotate90(Imaging imOut, Imaging imIn);
 extern Imaging ImagingRotate180(Imaging imOut, Imaging imIn);
 extern Imaging ImagingRotate270(Imaging imOut, Imaging imIn);
-extern Imaging ImagingResample(Imaging imIn, int xsize, int ysize, int filter);
+extern Imaging ImagingResample(Imaging imIn, int xsize, int ysize, int filter, float box[4]);
 extern Imaging ImagingTranspose(Imaging imOut, Imaging imIn);
 extern Imaging ImagingTransform(
     Imaging imOut, Imaging imIn, int method, int x0, int y0, int x1, int y1,

libImaging/Resample.c

@@ -134,16 +134,16 @@ static inline UINT8 clip8(int in)
 
 
 int
-precompute_coeffs(int inSize, int outSize, struct filter *filterp,
+precompute_coeffs(int inSize, float in0, float in1, int outSize,
-                  int **xboundsp, double **kkp) {
+                  struct filter *filterp, int **boundsp, double **kkp) {
     double support, scale, filterscale;
     double center, ww, ss;
-    int xx, x, kmax, xmin, xmax;
+    int xx, x, ksize, xmin, xmax;
-    int *xbounds;
+    int *bounds;
     double *kk, *k;
 
     /* prepare for horizontal stretch */
-    filterscale = scale = (double) inSize / outSize;
+    filterscale = scale = (double) (in1 - in0) / outSize;
     if (filterscale < 1.0) {
         filterscale = 1.0;
     }
@@ -152,27 +152,32 @@ precompute_coeffs(int inSize, int outSize, struct filter *filterp,
     support = filterp->support * filterscale;
 
     /* maximum number of coeffs */
-    kmax = (int) ceil(support) * 2 + 1;
+    ksize = (int) ceil(support) * 2 + 1;
 
     // check for overflow
-    if (outSize > INT_MAX / (kmax * sizeof(double)))
+    if (outSize > INT_MAX / (ksize * sizeof(double))) {
+        ImagingError_MemoryError();
         return 0;
+    }
 
     /* coefficient buffer */
     /* malloc check ok, overflow checked above */
-    kk = malloc(outSize * kmax * sizeof(double));
+    kk = malloc(outSize * ksize * sizeof(double));
-    if ( ! kk)
+    if ( ! kk) {
+        ImagingError_MemoryError();
         return 0;
+    }
 
-    /* malloc check ok, kmax*sizeof(double) > 2*sizeof(int) */
+    /* malloc check ok, ksize*sizeof(double) > 2*sizeof(int) */
-    xbounds = malloc(outSize * 2 * sizeof(int));
+    bounds = malloc(outSize * 2 * sizeof(int));
-    if ( ! xbounds) {
+    if ( ! bounds) {
         free(kk);
+        ImagingError_MemoryError();
         return 0;
     }
 
     for (xx = 0; xx < outSize; xx++) {
-        center = (xx + 0.5) * scale;
+        center = in0 + (xx + 0.5) * scale;
         ww = 0.0;
         ss = 1.0 / filterscale;
         // Round the value
@@ -184,7 +189,7 @@ precompute_coeffs(int inSize, int outSize, struct filter *filterp,
         if (xmax > inSize)
             xmax = inSize;
         xmax -= xmin;
-        k = &kk[xx * kmax];
+        k = &kk[xx * ksize];
         for (x = 0; x < xmax; x++) {
             double w = filterp->filter((x + xmin - center + 0.5) * ss);
             k[x] = w;
@@ -195,98 +200,75 @@ precompute_coeffs(int inSize, int outSize, struct filter *filterp,
                 k[x] /= ww;
         }
         // Remaining values should stay empty if they are used despite of xmax.
-        for (; x < kmax; x++) {
+        for (; x < ksize; x++) {
             k[x] = 0;
         }
-        xbounds[xx * 2 + 0] = xmin;
+        bounds[xx * 2 + 0] = xmin;
-        xbounds[xx * 2 + 1] = xmax;
+        bounds[xx * 2 + 1] = xmax;
     }
-    *xboundsp = xbounds;
+    *boundsp = bounds;
     *kkp = kk;
-    return kmax;
+    return ksize;
 }
 
 
-int
+void
-normalize_coeffs_8bpc(int outSize, int kmax, double *prekk, INT32 **kkp)
+normalize_coeffs_8bpc(int outSize, int ksize, double *prekk)
 {
     int x;
     INT32 *kk;
 
-    /* malloc check ok, overflow checked in precompute_coeffs */
+    // use the same buffer for normalized coefficients
-    kk = malloc(outSize * kmax * sizeof(INT32));
+    kk = (INT32 *) prekk;
-    if ( ! kk) {
-        return 0;
-    }
 
-    for (x = 0; x < outSize * kmax; x++) {
+    for (x = 0; x < outSize * ksize; x++) {
         if (prekk[x] < 0) {
             kk[x] = (int) (-0.5 + prekk[x] * (1 << PRECISION_BITS));
         } else {
             kk[x] = (int) (0.5 + prekk[x] * (1 << PRECISION_BITS));
         }
     }
-
-    *kkp = kk;
-    return kmax;
 }
 
 
 
-Imaging
+void    
-ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
+ImagingResampleHorizontal_8bpc(Imaging imOut, Imaging imIn, int offset,
+                               int ksize, int *bounds, double *prekk)
 {
     ImagingSectionCookie cookie;
-    Imaging imOut;
     int ss0, ss1, ss2, ss3;
-    int xx, yy, x, kmax, xmin, xmax;
+    int xx, yy, x, xmin, xmax;
-    int *xbounds;
     INT32 *k, *kk;
-    double *prekk;
 
-    kmax = precompute_coeffs(imIn->xsize, xsize, filterp, &xbounds, &prekk);
+    // use the same buffer for normalized coefficients
-    if ( ! kmax) {
+    kk = (INT32 *) prekk;
-        return (Imaging) ImagingError_MemoryError();
+    normalize_coeffs_8bpc(imOut->xsize, ksize, prekk);
-    }
-
-    kmax = normalize_coeffs_8bpc(xsize, kmax, prekk, &kk);
-    free(prekk);
-    if ( ! kmax) {
-        free(xbounds);
-        return (Imaging) ImagingError_MemoryError();
-    }
-
-    imOut = ImagingNewDirty(imIn->mode, xsize, imIn->ysize);
-    if ( ! imOut) {
-        free(kk);
-        free(xbounds);
-        return NULL;
-    }
 
     ImagingSectionEnter(&cookie);
     if (imIn->image8) {
         for (yy = 0; yy < imOut->ysize; yy++) {
-            for (xx = 0; xx < xsize; xx++) {
+            for (xx = 0; xx < imOut->xsize; xx++) {
-                xmin = xbounds[xx * 2 + 0];
+                xmin = bounds[xx * 2 + 0];
-                xmax = xbounds[xx * 2 + 1];
+                xmax = bounds[xx * 2 + 1];
-                k = &kk[xx * kmax];
+                k = &kk[xx * ksize];
                 ss0 = 1 << (PRECISION_BITS -1);
                 for (x = 0; x < xmax; x++)
-                    ss0 += ((UINT8) imIn->image8[yy][x + xmin]) * k[x];
+                    ss0 += ((UINT8) imIn->image8[yy + offset][x + xmin]) * k[x];
                 imOut->image8[yy][xx] = clip8(ss0);
             }
         }
     } else if (imIn->type == IMAGING_TYPE_UINT8) {
         if (imIn->bands == 2) {
             for (yy = 0; yy < imOut->ysize; yy++) {
-                for (xx = 0; xx < xsize; xx++) {
+                for (xx = 0; xx < imOut->xsize; xx++) {
-                    xmin = xbounds[xx * 2 + 0];
+                    xmin = bounds[xx * 2 + 0];
-                    xmax = xbounds[xx * 2 + 1];
+                    xmax = bounds[xx * 2 + 1];
-                    k = &kk[xx * kmax];
+                    k = &kk[xx * ksize];
                     ss0 = ss3 = 1 << (PRECISION_BITS -1);
                     for (x = 0; x < xmax; x++) {
-                        ss0 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 0]) * k[x];
+                        ss0 += ((UINT8) imIn->image[yy + offset][(x + xmin)*4 + 0]) * k[x];
-                        ss3 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 3]) * k[x];
+                        ss3 += ((UINT8) imIn->image[yy + offset][(x + xmin)*4 + 3]) * k[x];
                     }
                     ((UINT32 *) imOut->image[yy])[xx] = MAKE_UINT32(
                         clip8(ss0), 0, 0, clip8(ss3));
@@ -294,15 +276,15 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
             }
         } else if (imIn->bands == 3) {
             for (yy = 0; yy < imOut->ysize; yy++) {
-                for (xx = 0; xx < xsize; xx++) {
+                for (xx = 0; xx < imOut->xsize; xx++) {
-                    xmin = xbounds[xx * 2 + 0];
+                    xmin = bounds[xx * 2 + 0];
-                    xmax = xbounds[xx * 2 + 1];
+                    xmax = bounds[xx * 2 + 1];
-                    k = &kk[xx * kmax];
+                    k = &kk[xx * ksize];
                     ss0 = ss1 = ss2 = 1 << (PRECISION_BITS -1);
                     for (x = 0; x < xmax; x++) {
-                        ss0 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 0]) * k[x];
+                        ss0 += ((UINT8) imIn->image[yy + offset][(x + xmin)*4 + 0]) * k[x];
-                        ss1 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 1]) * k[x];
+                        ss1 += ((UINT8) imIn->image[yy + offset][(x + xmin)*4 + 1]) * k[x];
-                        ss2 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 2]) * k[x];
+                        ss2 += ((UINT8) imIn->image[yy + offset][(x + xmin)*4 + 2]) * k[x];
                     }
                     ((UINT32 *) imOut->image[yy])[xx] = MAKE_UINT32(
                         clip8(ss0), clip8(ss1), clip8(ss2), 0);
@@ -310,16 +292,16 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
             }
         } else {
             for (yy = 0; yy < imOut->ysize; yy++) {
-                for (xx = 0; xx < xsize; xx++) {
+                for (xx = 0; xx < imOut->xsize; xx++) {
-                    xmin = xbounds[xx * 2 + 0];
+                    xmin = bounds[xx * 2 + 0];
-                    xmax = xbounds[xx * 2 + 1];
+                    xmax = bounds[xx * 2 + 1];
-                    k = &kk[xx * kmax];
+                    k = &kk[xx * ksize];
                     ss0 = ss1 = ss2 = ss3 = 1 << (PRECISION_BITS -1);
                     for (x = 0; x < xmax; x++) {
-                        ss0 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 0]) * k[x];
+                        ss0 += ((UINT8) imIn->image[yy + offset][(x + xmin)*4 + 0]) * k[x];
-                        ss1 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 1]) * k[x];
+                        ss1 += ((UINT8) imIn->image[yy + offset][(x + xmin)*4 + 1]) * k[x];
-                        ss2 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 2]) * k[x];
+                        ss2 += ((UINT8) imIn->image[yy + offset][(x + xmin)*4 + 2]) * k[x];
-                        ss3 += ((UINT8) imIn->image[yy][(x + xmin)*4 + 3]) * k[x];
+                        ss3 += ((UINT8) imIn->image[yy + offset][(x + xmin)*4 + 3]) * k[x];
                     }
                     ((UINT32 *) imOut->image[yy])[xx] = MAKE_UINT32(
                         clip8(ss0), clip8(ss1), clip8(ss2), clip8(ss3));
@@ -327,50 +309,29 @@ ImagingResampleHorizontal_8bpc(Imaging imIn, int xsize, struct filter *filterp)
             }
         }
     }
-
     ImagingSectionLeave(&cookie);
-    free(kk);
-    free(xbounds);
-    return imOut;
 }
 
 
-Imaging
+void
-ImagingResampleVertical_8bpc(Imaging imIn, int ysize, struct filter *filterp)
+ImagingResampleVertical_8bpc(Imaging imOut, Imaging imIn, int offset,
+                             int ksize, int *bounds, double *prekk)
 {
     ImagingSectionCookie cookie;
-    Imaging imOut;
     int ss0, ss1, ss2, ss3;
-    int xx, yy, y, kmax, ymin, ymax;
+    int xx, yy, y, ymin, ymax;
-    int *xbounds;
     INT32 *k, *kk;
-    double *prekk;
 
-    kmax = precompute_coeffs(imIn->ysize, ysize, filterp, &xbounds, &prekk);
+    // use the same buffer for normalized coefficients
-    if ( ! kmax) {
+    kk = (INT32 *) prekk;
-        return (Imaging) ImagingError_MemoryError();
+    normalize_coeffs_8bpc(imOut->ysize, ksize, prekk);
-    }
-
-    kmax = normalize_coeffs_8bpc(ysize, kmax, prekk, &kk);
-    free(prekk);
-    if ( ! kmax) {
-        free(xbounds);
-        return (Imaging) ImagingError_MemoryError();
-    }
-
-    imOut = ImagingNewDirty(imIn->mode, imIn->xsize, ysize);
-    if ( ! imOut) {
-        free(kk);
-        free(xbounds);
-        return NULL;
-    }
 
     ImagingSectionEnter(&cookie);
     if (imIn->image8) {
-        for (yy = 0; yy < ysize; yy++) {
+        for (yy = 0; yy < imOut->ysize; yy++) {
-            k = &kk[yy * kmax];
+            k = &kk[yy * ksize];
-            ymin = xbounds[yy * 2 + 0];
+            ymin = bounds[yy * 2 + 0];
-            ymax = xbounds[yy * 2 + 1];
+            ymax = bounds[yy * 2 + 1];
             for (xx = 0; xx < imOut->xsize; xx++) {
                 ss0 = 1 << (PRECISION_BITS -1);
                 for (y = 0; y < ymax; y++)
@@ -380,10 +341,10 @@ ImagingResampleVertical_8bpc(Imaging imIn, int ysize, struct filter *filterp)
         }
     } else if (imIn->type == IMAGING_TYPE_UINT8) {
         if (imIn->bands == 2) {
-            for (yy = 0; yy < ysize; yy++) {
+            for (yy = 0; yy < imOut->ysize; yy++) {
-                k = &kk[yy * kmax];
+                k = &kk[yy * ksize];
-                ymin = xbounds[yy * 2 + 0];
+                ymin = bounds[yy * 2 + 0];
-                ymax = xbounds[yy * 2 + 1];
+                ymax = bounds[yy * 2 + 1];
                 for (xx = 0; xx < imOut->xsize; xx++) {
                     ss0 = ss3 = 1 << (PRECISION_BITS -1);
                     for (y = 0; y < ymax; y++) {
@@ -395,10 +356,10 @@ ImagingResampleVertical_8bpc(Imaging imIn, int ysize, struct filter *filterp)
                 }
             }
         } else if (imIn->bands == 3) {
-            for (yy = 0; yy < ysize; yy++) {
+            for (yy = 0; yy < imOut->ysize; yy++) {
-                k = &kk[yy * kmax];
+                k = &kk[yy * ksize];
-                ymin = xbounds[yy * 2 + 0];
+                ymin = bounds[yy * 2 + 0];
-                ymax = xbounds[yy * 2 + 1];
+                ymax = bounds[yy * 2 + 1];
                 for (xx = 0; xx < imOut->xsize; xx++) {
                     ss0 = ss1 = ss2 = 1 << (PRECISION_BITS -1);
                     for (y = 0; y < ymax; y++) {
@@ -411,10 +372,10 @@ ImagingResampleVertical_8bpc(Imaging imIn, int ysize, struct filter *filterp)
                 }
             }
         } else {
-            for (yy = 0; yy < ysize; yy++) {
+            for (yy = 0; yy < imOut->ysize; yy++) {
-                k = &kk[yy * kmax];
+                k = &kk[yy * ksize];
-                ymin = xbounds[yy * 2 + 0];
+                ymin = bounds[yy * 2 + 0];
-                ymax = xbounds[yy * 2 + 1];
+                ymax = bounds[yy * 2 + 1];
                 for (xx = 0; xx < imOut->xsize; xx++) {
                     ss0 = ss1 = ss2 = ss3 = 1 << (PRECISION_BITS -1);
                     for (y = 0; y < ymax; y++) {
@@ -429,47 +390,30 @@ ImagingResampleVertical_8bpc(Imaging imIn, int ysize, struct filter *filterp)
             }
         }
     }
-
     ImagingSectionLeave(&cookie);
-    free(kk);
-    free(xbounds);
-    return imOut;
 }
 
 
-Imaging
+void
-ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
+ImagingResampleHorizontal_32bpc(Imaging imOut, Imaging imIn, int offset,
+                                int ksize, int *bounds, double *kk)
 {
     ImagingSectionCookie cookie;
-    Imaging imOut;
     double ss;
-    int xx, yy, x, kmax, xmin, xmax;
+    int xx, yy, x, xmin, xmax;
-    int *xbounds;
+    double *k;
-    double *k, *kk;
-
-    kmax = precompute_coeffs(imIn->xsize, xsize, filterp, &xbounds, &kk);
-    if ( ! kmax) {
-        return (Imaging) ImagingError_MemoryError();
-    }
-
-    imOut = ImagingNewDirty(imIn->mode, xsize, imIn->ysize);
-    if ( ! imOut) {
-        free(kk);
-        free(xbounds);
-        return NULL;
-    }
 
     ImagingSectionEnter(&cookie);
     switch(imIn->type) {
         case IMAGING_TYPE_INT32:
             for (yy = 0; yy < imOut->ysize; yy++) {
-                for (xx = 0; xx < xsize; xx++) {
+                for (xx = 0; xx < imOut->xsize; xx++) {
-                    xmin = xbounds[xx * 2 + 0];
+                    xmin = bounds[xx * 2 + 0];
-                    xmax = xbounds[xx * 2 + 1];
+                    xmax = bounds[xx * 2 + 1];
-                    k = &kk[xx * kmax];
+                    k = &kk[xx * ksize];
                     ss = 0.0;
                     for (x = 0; x < xmax; x++)
-                        ss += IMAGING_PIXEL_I(imIn, x + xmin, yy) * k[x];
+                        ss += IMAGING_PIXEL_I(imIn, x + xmin, yy + offset) * k[x];
                     IMAGING_PIXEL_I(imOut, xx, yy) = ROUND_UP(ss);
                 }
             }
@@ -477,55 +421,38 @@ ImagingResampleHorizontal_32bpc(Imaging imIn, int xsize, struct filter *filterp)
 
         case IMAGING_TYPE_FLOAT32:
             for (yy = 0; yy < imOut->ysize; yy++) {
-                for (xx = 0; xx < xsize; xx++) {
+                for (xx = 0; xx < imOut->xsize; xx++) {
-                    xmin = xbounds[xx * 2 + 0];
+                    xmin = bounds[xx * 2 + 0];
-                    xmax = xbounds[xx * 2 + 1];
+                    xmax = bounds[xx * 2 + 1];
-                    k = &kk[xx * kmax];
+                    k = &kk[xx * ksize];
                     ss = 0.0;
                     for (x = 0; x < xmax; x++)
-                        ss += IMAGING_PIXEL_F(imIn, x + xmin, yy) * k[x];
+                        ss += IMAGING_PIXEL_F(imIn, x + xmin, yy + offset) * k[x];
                     IMAGING_PIXEL_F(imOut, xx, yy) = ss;
                 }
             }
             break;
     }
-
     ImagingSectionLeave(&cookie);
-    free(kk);
-    free(xbounds);
-    return imOut;
 }
 
 
-Imaging
+void
-ImagingResampleVertical_32bpc(Imaging imIn, int ysize, struct filter *filterp)
+ImagingResampleVertical_32bpc(Imaging imOut, Imaging imIn, int offset,
+                              int ksize, int *bounds, double *kk)
 {
     ImagingSectionCookie cookie;
-    Imaging imOut;
     double ss;
-    int xx, yy, y, kmax, ymin, ymax;
+    int xx, yy, y, ymin, ymax;
-    int *xbounds;
+    double *k;
-    double *k, *kk;
-
-    kmax = precompute_coeffs(imIn->ysize, ysize, filterp, &xbounds, &kk);
-    if ( ! kmax) {
-        return (Imaging) ImagingError_MemoryError();
-    }
-
-    imOut = ImagingNewDirty(imIn->mode, imIn->xsize, ysize);
-    if ( ! imOut) {
-        free(kk);
-        free(xbounds);
-        return NULL;
-    }
 
     ImagingSectionEnter(&cookie);
     switch(imIn->type) {
         case IMAGING_TYPE_INT32:
-            for (yy = 0; yy < ysize; yy++) {
+            for (yy = 0; yy < imOut->ysize; yy++) {
-                ymin = xbounds[yy * 2 + 0];
+                ymin = bounds[yy * 2 + 0];
-                ymax = xbounds[yy * 2 + 1];
+                ymax = bounds[yy * 2 + 1];
-                k = &kk[yy * kmax];
+                k = &kk[yy * ksize];
                 for (xx = 0; xx < imOut->xsize; xx++) {
                     ss = 0.0;
                     for (y = 0; y < ymax; y++)
@@ -536,10 +463,10 @@ ImagingResampleVertical_32bpc(Imaging imIn, int ysize, struct filter *filterp)
             break;
 
         case IMAGING_TYPE_FLOAT32:
-            for (yy = 0; yy < ysize; yy++) {
+            for (yy = 0; yy < imOut->ysize; yy++) {
-                ymin = xbounds[yy * 2 + 0];
+                ymin = bounds[yy * 2 + 0];
-                ymax = xbounds[yy * 2 + 1];
+                ymax = bounds[yy * 2 + 1];
-                k = &kk[yy * kmax];
+                k = &kk[yy * ksize];
                 for (xx = 0; xx < imOut->xsize; xx++) {
                     ss = 0.0;
                     for (y = 0; y < ymax; y++)
@@ -549,22 +476,27 @@ ImagingResampleVertical_32bpc(Imaging imIn, int ysize, struct filter *filterp)
             }
             break;
     }
-
     ImagingSectionLeave(&cookie);
-    free(kk);
-    free(xbounds);
-    return imOut;
 }
 
 
+typedef void (*ResampleFunction)(Imaging imOut, Imaging imIn, int offset,
+                                 int ksize, int *bounds, double *kk);
+
+
 Imaging
-ImagingResample(Imaging imIn, int xsize, int ysize, int filter)
+ImagingResampleInner(Imaging imIn, int xsize, int ysize,
+                     struct filter *filterp, float box[4],
+                     ResampleFunction ResampleHorizontal,
+                     ResampleFunction ResampleVertical);
+
+
+Imaging
+ImagingResample(Imaging imIn, int xsize, int ysize, int filter, float box[4])
 {
-    Imaging imTemp = NULL;
-    Imaging imOut = NULL;
     struct filter *filterp;
-    Imaging (*ResampleHorizontal)(Imaging imIn, int xsize, struct filter *filterp);
+    ResampleFunction ResampleHorizontal;
-    Imaging (*ResampleVertical)(Imaging imIn, int xsize, struct filter *filterp);
+    ResampleFunction ResampleVertical;
 
     if (strcmp(imIn->mode, "P") == 0 || strcmp(imIn->mode, "1") == 0)
         return (Imaging) ImagingError_ModeError();
@@ -613,24 +545,93 @@ ImagingResample(Imaging imIn, int xsize, int ysize, int filter)
             );
     }
 
-    /* two-pass resize, first pass */
+    return ImagingResampleInner(imIn, xsize, ysize, filterp, box,
-    if (imIn->xsize != xsize) {
+                                ResampleHorizontal, ResampleVertical);
-        imTemp = ResampleHorizontal(imIn, xsize, filterp);
+}
-        if ( ! imTemp)
+
+
+Imaging
+ImagingResampleInner(Imaging imIn, int xsize, int ysize,
+                     struct filter *filterp, float box[4],
+                     ResampleFunction ResampleHorizontal,
+                     ResampleFunction ResampleVertical)
+{
+    Imaging imTemp = NULL;
+    Imaging imOut = NULL;
+
+    int i;
+    int yroi_min, yroi_max;
+    int ksize_horiz, ksize_vert;
+    int *bounds_horiz, *bounds_vert;
+    double *kk_horiz, *kk_vert;
+
+    ksize_horiz = precompute_coeffs(imIn->xsize, box[0], box[2], xsize,
+                                    filterp, &bounds_horiz, &kk_horiz);
+    if ( ! ksize_horiz) {
         return NULL;
+    }
+
+    ksize_vert = precompute_coeffs(imIn->ysize, box[1], box[3], ysize,
+                                   filterp, &bounds_vert, &kk_vert);
+    if ( ! ksize_vert) {
+        free(bounds_horiz);
+        free(kk_horiz);
+        return NULL;
+    }
+
+    // First used row in the source image
+    yroi_min = bounds_vert[0];
+    // Last used row in the source image
+    yroi_max = bounds_vert[ysize*2 - 2] + bounds_vert[ysize*2 - 1];
+
+
+    /* two-pass resize, first pass */
+    if (box[0] || box[2] != xsize) {
+        // Shift bounds for vertical pass
+        for (i = 0; i < ysize; i++) {
+            bounds_vert[i * 2] -= yroi_min;
+        }
+
+        imTemp = ImagingNewDirty(imIn->mode, xsize, yroi_max - yroi_min);
+        if (imTemp) {
+            ResampleHorizontal(imTemp, imIn, yroi_min,
+                               ksize_horiz, bounds_horiz, kk_horiz);
+        }
+        free(bounds_horiz);
+        free(kk_horiz);
+        if ( ! imTemp) {
+            free(bounds_vert);
+            free(kk_vert);
+            return NULL;
+        }
         imOut = imIn = imTemp;
+    } else {
+        // Free in any case
+        free(bounds_horiz);
+        free(kk_horiz);
     }
 
     /* second pass */
-    if (imIn->ysize != ysize) {
+    if (box[1] || box[3] != ysize) {
+        imOut = ImagingNewDirty(imIn->mode, imIn->xsize, ysize);
+        if (imOut) {
             /* imIn can be the original image or horizontally resampled one */
-        imOut = ResampleVertical(imIn, ysize, filterp);
+            ResampleVertical(imOut, imIn, 0,
+                             ksize_vert, bounds_vert, kk_vert);
+        }
         /* it's safe to call ImagingDelete with empty value
-           if there was no previous step. */
+           if previous step was not performed. */
         ImagingDelete(imTemp);
-        if ( ! imOut)
+        free(bounds_vert);
+        free(kk_vert);
+        if ( ! imOut) {
             return NULL;
         }
+    } else {
+        // Free in any case
+        free(bounds_vert);
+        free(kk_vert);
+    }
 
     /* none of the previous steps are performed, copying */
     if ( ! imOut) {