# How can I sum pixel values over a rotated rectangle?

I have an optimization problem in which I need to sum pixel values in an image over a rectangular region. This is a core component of the optimization so it will be done often and the naive solution is not fast enough.

The region can be rotated arbitrarily so that I can't use an integral image; integral images only work for axis-aligned rectangles.

Does someone know a generalization of the integral images for that problem or has a different idea?

Update:

• The rectangles have different rotations, so that I can't rotate the rect.
• The rectangles have a size similar to the one of the full image (each ~20%)
• You're going to do this many times, for many different rotated rectangles. Will all of the rectangles be rotated the same amount? Or will they be different rotations each time? If the former, this is easy (rotate the entire image, then calculate an integral image for its rotation). The latter seems much harder.
– D.W.
Commented Jan 19, 2016 at 17:29
• @EvilJS, nice find! Write an answer? (It looks like their technique is similar to my answer, but they have developed it much further -- and it's great to see a fully-worked out paper on the subject.)
– D.W.
Commented Jan 19, 2016 at 19:50
• How large are your rectangles ?
– user16034
Commented Jan 19, 2016 at 20:49

One simple technique that's probably not optimal, but is better than naively enumerating all pixels in the rotated rectangle, is to use "integral columns" (thanks to Yves Daoust for the name):

Preprocessing: Let $I[\cdot,\cdot]$ be the image. For each $x$ value, build a 1D integral image for $I[x,\cdot]$. In other words, we build a table $S[\cdot,\cdot]$ where $S[x,y] = S[x,0] + S[x,1] + \dots + S[x,y]$. This can be done in linear time.

Answering a query: Given a rotated rectangle $R$, we find the range of $x$-values of pixels within $R$ (i.e., the leftmost part of $R$ and rightmost part of $R$). This can be done in $O(1)$ time. Next, for each $x$ in that range, we find the range $[y_\ell,y_u]$ of $y$-values such that $(x,y)$ is in $R$ for each $y_\ell \le y \le y_u$. We sum the pixel values of those by computing $S[x,y_u] - S[x,y_{\ell}-1]$. This takes $O(1)$ time per $x$-value, so the total running time to sum the pixels in a rectangle is proportional to the width of $R$ (in contrast, the naive method's running time is proportional to the area of $R$).

As a trivial optimization, you can precompute sum tables for both $x$ and $y$; for each $R$, you can check which is smaller, the width or height of $R$, and then decide whether to enumerate $x$-values or $y$-values.

• Integral columns ?
– user16034
Commented Jan 19, 2016 at 20:49
• I am not so sure that this isn't optimal. Integral images are efficient on axis-aligned rectangles, which are separable shapes (Cartesian product of two intervals). Rotated rectangles do not have AA-sides (just short runs of aligned pixels), and they have up to $2L$ "corners", where $L$ is the longest dimension in pixels, so that any decomposition in AA-rectangles requires $O(L)$ of them.
– user16034
Commented Jan 19, 2016 at 21:05
• @YvesDaoust, excellent name!
– D.W.
Commented Jan 19, 2016 at 23:56
• Thank you for the approach. It's easily understandable, looks implementable and can be parallelized perfectly. Commented Jan 20, 2016 at 17:03

Continuing my comment on the excellent answer of D.W.:

The integral images allow summation in constant time over an axis-aligned rectangle because they are the Cartesian product of two intervals.

You can speedup the summation over an arbitrary shape by partitioning it in a minimum number of axis-aligned rectangles. As shown on the figure, the possibilities are quite limited on a rotated rectangle.

In the worst case, rotation by $45°$, rectangles larger than one pixel are just impossible. There can be larger rectangles for small rotations, but in practice they are not really worth the effort. Note that it is not the area of the rectangles that matters, but their number.

Generally speaking, the total workload must be proportional to the length of the outline in the worst case. A comparable result (if not the same in disguise) is developed in "Yang, Luren, Albregtsen, Fritz, Fast and exact computation of moments using discrete Green's theorem. NOBIM-konferansen. 1994".

Integral images can be adapted to other angles, by summing in non axis-parallel directions, allowing to sum inside parallelograms.

Anyway, this is of little use as it doesn't lower the number of pieces required in the decomposition of arbitrarily rotated rectangles.

There is a quite fresh publication about Line Histogram which is beneficial to use, as this technique is less expensive than rotating image and uses simple line based histograms - which gives computational boost while still being flexible with given angles.

• Thanks. I just had a look an the paper and it also looks promising. Only strange thing is that they tested the algorithm on a Win XP(!) machine in 2012... Commented Jan 20, 2016 at 17:11

a quite old article

Exact Integral Images at Generic Angles for 2D Barcode Detection

abstract: We propose a novel solution to elegantly compute integral images at generic angles. Our method is exact in the sense that no approximations are used to derive it and it is vulnerable only to the unavoidable aliasing effects of discretization.

In their paper, the determination of a block sum for an arbitrary angle, is as same as the conventional one

block_sum = A - B - C + D

which is elegant, in my point of view.

• Could you write some abstract or summarize what this article describes? Could you properly attribute the article in case the link rots?
– Evil
Commented Nov 22, 2017 at 3:46

An approach to calculate approximate image integrals of arbitrary rotated rectangles up to desired accuracy is as follows (more accurate->slower):

An arbitrary rotated rectangle (defined by clockwise points p1 to p4) is subdivided into a maximum enclosed normal rectangle (R5) and four axis-aligned right triangles as shown below (R1 to R4).

One can use the ordinary integral image to compute (approximations of) all of these.

Having the normal integral image $$\text{sat}[x,y]$$ (eg as used by Viola,Jones) the algorithms to compute the rectangle, the triangles and rotated rectangle sums are as follows (pseudocode adapted from my Image Processing Library FILTER.js)

function satsum(sat, x0, y0, x1, y1)
{
// exact sat sum of axis-aligned rectangle defined by p0, p1 (top left, bottom right)
return sat[x1,y1] - sat[x0-1,y1] - sat[x1,y0-1] + sat[x0-1,y0-1];
}
function satsumt(sat, x0, y0, x1, y1, x2, y2, k)
{
// approximate sat sum of axis-aligned right triangle defined by p0, p1, p2
// k is number of subdivisions up to desired accuracy, k=Infinity gives maximum accuracy
// we approximate the triangle area by subdividing it into smaller rectangle sat sums
xm = min(x0, x1, x2);
ym = min(y0, y1, y2);
xM = max(x0, x1, x2);
yM = max(y0, y1, y2);
dx = xM - xm;
dy = yM - ym;
s = 0;
if (0 === dx || 0 === dy)
{
// zero area
s = 0;
}
else if (k <= 1)
{
//most crude approximation, half of enclosing rectangle sum
s = satsum(sat, xm, ym, xM, yM)/2;
}
else if (dx > dy)
{
//better approximation, subdivide along y
if (y0 === ym)
{
if (y1 === ym) xym = x2 === xm ? xM : xm;
else if (y2 === ym) xym = x1 === xm ? xM : xm;
else xym = x0;
}
else if (y1 === ym)
{
if (y2 === ym) xym = x0 === xm ? xM : xm;
else xym = x1;
}
else
{
xym = x2;
}
if (y0 === yM)
{
if (y1 === yM) xyM = x2 === xm ? xM : xm;
else if (y2 === yM) xyM = x1 === xm ? xM : xm;
else xyM = x0;
}
else if (y1 === yM)
{
if (y2 === yM) xyM = x0 === xm ? xM : xm;
else xyM = x1;
}
else
{
xyM = x2;
}
d = max(1, round(dy/k));
for (i=1,p=ym,y=ym+d; i<=k; ++i,y+=d)
{
if (y > yM) y = yM;
xi = round(xym + ((y-ym)/dy)*(xyM-xym));
s += satsum(sat, min(xym, xi), p, max(xym, xi), y);
if (y >= yM) break;
p = min(y+1, yM);
}
if (y < yM) s += satsum(sat, min(xym, xyM), y+1, max(xym, xyM), yM);
}
else
{
//better approximation, subdivide along x
if (x0 === xm)
{
if (x1 === xm) yxm = y2 === ym ? yM : ym;
else if (x2 === xm) yxm = y1 === ym ? yM : ym;
else yxm = y0;
}
else if (x1 === xm)
{
if (x2 === xm) yxm = y0 === ym ? yM : ym;
else yxm = y1;
}
else
{
yxm = y2;
}
if (x0 === xM)
{
if (x1 === xM) yxM = y2 === ym ? yM : ym;
else if (x2 === xM) yxM = y1 === ym ? yM : ym;
else yxM = y0;
}
else if (x1 === xM)
{
if (x2 === xM) yxM = y0 === ym ? yM : ym;
else yxM = y1;
}
else
{
yxM = y2;
}
d = max(1, round(dx/k));
for (i=1,p=xm,x=xm+d; i<=k; ++i,x+=d)
{
if (x > xM) x = xM;
yi = round(yxm + ((x-xm)/dx)*(yxM-yxm));
s += satsum(sat, p, min(yxm, yi), x, max(yxm, yi));
if (x >= xM) break;
p = min(x+1, xM);
}
if (x < xM) s += satsum(sat, x+1, min(yxm, yxM), xM, max(yxm, yxM));
}
return s;
}
function satsumr(sat, x1, y1, x2, y2, x3, y3, x4, y4, k)
{
// approximate sat sum for arbitrary rotated rectangle defined (clockwise) by p1 to p4
xm = min(x1, x2, x3, x4);
ym = min(y1, y2, y3, y4);
xM = max(x1, x2, x3, x4);
yM = max(y1, y2, y3, y4);
// (xm,ym), (xM,yM) is the normal rectangle enclosing the rotated rectangle
// (min(xi1, xi2),min(yi1, yi2)), (max(xi1, xi2),max(yi1, yi2)) is the maximum normal rectangle enclosed by the rotated rectangle computed by satsum
// the rest of the rotated rectangle are 4 axis-aligned right triangles computed approximately by satsumt

if (y1 === y2 || y2 === y3 || y3 === y4 || y4 === y1) return satsum(sat, xm, ym, xM, yM); // axis-aligned unrotated rectangle

if (y1 === ym) xi1 = x1;
else if (y2 === ym) xi1 = x2;
else if (y3 === ym) xi1 = x3;
else xi1 = x4;
if (y1 === yM) xi2 = x1;
else if (y2 === yM) xi2 = x2;
else if (y3 === yM) xi2 = x3;
else xi2 = x4;
if (x1 === xm) yi1 = y1;
else if (x2 === xm) yi1 = y2;
else if (x3 === xm) yi1 = y3;
else yi1 = y4;
if (x1 === xM) yi2 = y1;
else if (x2 === xM) yi2 = y2;
else if (x3 === xM) yi2 = y3;
else yi2 = y4;
return (
satsumt(sat, xm, yi1, xi1, ym, xi1, yi1, k) + // top left right triagle
satsumt(sat, xM, yi2, xi1, ym, xi1, yi2, k) + // top right right triagle
satsumt(sat, xm, yi1, xi2, yM, xi2, yi1, k) + // bottom left right triagle
satsumt(sat, xM, yi2, xi2, yM, xi2, yi2, k) + // bottom right right triagle
satsum (sat, min(xi1, xi2), min(yi1, yi2), max(xi1, xi2), max(yi1, yi2)) // center rectangle
);
}


References: