The imos method

Published on December 23, 2012

The imos method

The imos method extends the prefix-sum (cumulative-sum) algorithm to multiple dimensions and multiple degrees. In competitive programming, only the two-dimensional, first-degree variant tends to appear, but as a 2012 research result it turns out that applying the imos method in higher-dimensional spaces and with higher degrees is useful in the fields of signal processing and image processing.

The basics: the 1-dimensional, 0-degree imos method

The simplest form of the imos method adds a 0-degree function (a function whose top is flat, such as a rectangular function or a step function) on a one-dimensional line.
Figure 1: An image of adding a 0-degree function on a one-dimensional line

Example problem

You run a coffee shop. For each customer i\ (0 \leq i \lt C) who visits your shop, you are given an arrival time S_i and a departure time E_i (0 \leq S_i \lt E_i \leq T). What is the maximum number of customers M who were in the shop at the same time? When a departure and an arrival happen at the same instant, assume the departure happens first. (As for the solution, assume that even when C is very large you are expected to solve it in O(C+T) time.)

Naive solution

A naive solution is to add 1 to the count of every time step during which each customer stays. However, this takes O(CT) time. 
for (int i = 0; i < T; i++) table[i] = 0;
for (int i = 0; i < C; i++) {
  // For each time from S[i] to E[i] - 1, increment the count by 1
  for (int j = S[i]; j < E[i]; j++) {
    table[j]++;
  }
}
// Find the maximum
M = 0;
for (int i = 0; i < T; i++) {
  if (M < table[i]) M = table[i];
}

Solution using the imos method

The idea of the imos method is to record only what happens at the entrance and exit: record +1 on arrival and -1 on departure, and simulate the whole thing just before computing the answer. Recording takes O(C) and the simulation takes O(T), so the overall time complexity is O(C+T). 
for (int i = 0; i < T; i++) table[i] = 0;
for (int i = 0; i < C; i++) {
  table[S[i]]++;  // Arrival: increment the count by 1
  table[E[i]]--;  // Departure: decrement the count by 1
}
// Simulate
for (int i = 0; i < T; i++) {
  if (0 < i) table[i] += table[i - 1];
}
// Find the maximum
M = 0;
for (int i = 0; i < T; i++) {
  if (M < table[i]) M = table[i];
}

Extension to more dimensions

The greatest feature of the imos method is that, compared with the naive approach, it scales well as the number of dimensions grows. You cannot escape the curse of dimensionality, so applying it in very high dimensions is difficult; but two- and three-dimensional problem settings are extremely common, and the imos method is useful in those situations.

Example problem

You are playing a game in which you catch various kinds of monsters. Right now you are in a field of grass made up of W\times H tiles. In this grass, N kinds of monsters appear. Monster i appears only in the range A_i ≦ x < B_i,\ C_i ≦ y < D_i. To play the game efficiently, you want to know the maximum number of monster kinds that appear on any single tile. (Assume W\times H is small enough to iterate over, while N can be large.)
Figure 2: The rectangles (in red) where monsters appear, and
the number of monster kinds appearing on each tile

Naive solution

A naive solution is to add 1 to every tile contained in the rectangle where each monster appears. However, this takes O(NWH) time. 
for (int y = 0; y < H; y++) {
  for (int x = 0; x < W; x++) {
    tiles[y][x] = 0;
  }
}
for (int i = 0; i < N; i++) {
  // Add 1 to the range [(A[i],C[i]), (B[i],D[i])) where monster i appears
  for (int y = C[i]; y < D[i]; y++) {
    for (int x = A[i]; x < B[i]; x++) {
      tiles[y][x]++;
    }
  }
}
// Find the maximum
int tile_max = 0;
for (int y = 0; y < H; y++) {
  for (int x = 0; x < W; x++) {
    if (tile_max < tiles[y][x]) {
      tile_max = tiles[y][x];
    }
  }
}

Solution using the imos method

In the imos method, add +1 at the top-left corner (A[i],C[i]) of the rectangle, -1 at the top-right (A[i],D[i]), -1 at the bottom-left (B[i],C[i]), and +1 at the bottom-right (B[i],D[i]), then take prefix sums just before computing the answer. Recording takes O(N) and computing the prefix sums takes O(WH), so the overall time complexity is O(N+WH). 
for (int y = 0; y < H; y++) {
  for (int x = 0; x < W; x++) {
    tiles[y][x] = 0;
  }
}
// Adding the weights (Figure 3)
for (int i = 0; i < N; i++) {
  tiles[C[i]][A[i]]++;
  tiles[C[i]][B[i]]--;
  tiles[D[i]][A[i]]--;
  tiles[D[i]][B[i]]++;
}
// Prefix sum in the horizontal direction (Figures 4, 5)
for (int y = 0; y < H; y++) {
  for (int x = 1; x < W; x++) {
    tiles[y][x] += tiles[y][x - 1];
  }
}
// Prefix sum in the vertical direction (Figures 6, 7)
for (int y = 1; y < H; y++) {
  for (int x = 0; x < W; x++) {
    tiles[y][x] += tiles[y - 1][x];
  }
}
// Find the maximum
int tile_max = 0;
for (int y = 0; y < H; y++) {
  for (int x = 0; x < W; x++) {
    if (tile_max < tiles[y][x]) {
      tile_max = tiles[y][x];
    }
  }
}
Figure 3: The result of adding the weights
Figure 4: Accumulating in the horizontal direction
Figure 5: The result of accumulating horizontally
Figure 6: Accumulating in the vertical direction
Figure 7: The result of accumulating vertically

Extension to special coordinate systems

The imos method can also be applied to triangular and hexagonal coordinate systems. To decide where to add the weights and in which directions to accumulate, take differences in various directions with respect to the target shape so that the number of weighted cells becomes as small as possible.
Figure 8: A triangular coordinate system

The differencing technique

Figure 9 shows the weights to be filled in on a triangular coordinate system. The green lines indicate that a portion in between is omitted (for a triangle whose side length is 4, ignoring the lines gives the same result in these figures). To perform such weighting with the imos method, you need to take differences until the number of weighted cells becomes small enough, so as to learn the placement of the cells to weight and the directions in which to accumulate.
Figures 9 through 12 show the process of taking differences. Between them we took the "difference from top-left to bottom-right," the "difference from top-right to bottom-left," and the "difference from left to right." So, to compute a triangular prefix sum with the imos method, we place the values as in Figure 12 and perform the reverse operations: "prefix sum from left to right," "prefix sum from top-right to bottom-left," and "prefix sum from top-left to bottom-right."
Figure 9: The weights to be filled in
Figure 10: The top-left-to-bottom-right difference of Figure 9
Figure 11: The top-right-to-bottom-left difference of Figure 10
Figure 12: The left-to-right difference of Figure 11
The same method applies not only to triangles but also to hexagons, and to cases where the weight added is not a constant but has a slope (as long as it is a polynomial).

Extension to higher degrees

Problems posed in competitive programming are probably at most first degree, but the imos method can also handle high-degree piecewise-polynomial functions. Functions that are usually thought not to admit a good polynomial approximation can sometimes be expressed as piecewise polynomials, and exploiting this makes it possible to speed up the continuous wavelet transform.

Differences of a quadratic function

As a trial, let us perform the imos method for a quadratic function on a one-dimensional line.
01234 56789
00014 916253649
↓difference　　　↑prefix sum
01234 56789
00013 5791113
↓difference　　　↑prefix sum
01234 56789
00012 22222
↓difference　　　↑prefix sum
01234 56789
00011 00000
Figure 13: Differences of the quadratic function (x-2)_+^2
From Figure 13, we see that to add the quadratic function (x-a)_+^2, we add 1 to positions a+1 and a+2 of the table and finally take the prefix sum three times, which lets us apply the imos method using a quadratic function (where (･)_+ = {\rm max}(0,･)). By the same method, the imos method can be applied even to higher-degree functions.

Gaussian functions and the imos method

As an example of the higher-degree extension, let us look at the Gaussian function. To apply the imos method, we approximate the Gaussian function on the left-hand side of Equation 1 by a piecewise polynomial like the one on the right-hand side. As shown in Figure 14, this approximating function fits extremely well.
Equation 1: The Gaussian function (left) and its approximation (right)
Figure 14: The Gaussian function (blue) and its polynomial approximation (red)
-12-11-10 -9-8-7-6-5 -4-3-2-10 12
003 30000 00-11-110 00
↓prefix sum　　　↑difference
-12-11-10 -9-8-7-6-5 -4-3-2-10 12
003 66666 66-5-16-16 -16-16
↓prefix sum　　　↑difference
-12-11-10 -9-8-7-6-5 -4-3-2-10 12
003 915212733 394540248 -8-24
↓prefix sum　　　↑difference
-12-11-10 -9-8-7-6-5 -4-3-2-10 12
003 12274875108 147192232256264 256232
Figure 15: The Gaussian function approximated by the imos method
-12-11-10 -9-8-7-6-5 -4-3-2-10 12
1
.49		3
.42		7
.28		 14
.41		26
.57		45
.58		72
.76		108
.08		 149
.41		192
.20		230
.09		256
.31		265
.70		 256
.31		230
.09
Figure 16: The target Gaussian function to be approximated
Furthermore, since a d-dimensional Gaussian function can be expressed as a product of one-dimensional Gaussian functions, the d-dimensional imos method can be applied. And because the weight of a Gaussian function is determined by the Euclidean distance from its center (a circle in two dimensions, a sphere in three dimensions), it is extremely useful when you want a rotation-invariant value.
If you would like to know more about this, please refer to the presentation material from ICPR 2012.

Appendix

Problems where the imos method can be used

  • [Easy] Osaki … 1-dimensional 0-degree imos method, ACM-ICPC Japan Alumni Group Practice Contest for Japan Domestic 2007 - [AOJ]
  • [Easy] Nails … 0-degree imos method on triangular coordinates, Japanese Olympiad in Informatics 2012 Final Round Problem 4 - [AOJ]
  • [Medium] Paint Colors … coordinate compression + 2-dimensional 0-degree imos method + breadth-first search (depth-first search is also possible depending on stack capacity), Japanese Olympiad in Informatics 2008 Final Round Problem 5 - [AOJ]
  • [Medium] Plugs … 2-dimensional 0-degree imos method + ad hoc, Japanese Olympiad in Informatics 2010 Training Camp Day 4 Problem 4
  • [Hard] Pyramid … first-degree imos method on a special coordinate system, Autumn Fest 2012: I