Hi! When I worked at Synopsys, I optimized algorithms for a specialized embedded vision processor — the kind used in autonomous driving and ADAS. In this post I'll cover how to optimize an algorithm that estimates the distance to objects from a pair of images. It's called Stereo Block Matching, and the prefix sums method is what makes it fast. Done naively, the same algorithm can run more than 100x slower.
The stereo pair
Let's look at an algorithm that estimates the distance from a camera to objects from a stereo pair. A stereo pair can be captured with two cameras pointed at the same spot, like a pair of eyes. It's two images — one from the left camera and one from the right. You can overlay them on top of each other and get a 3D movie. Here's what a stereo pair looks like:
An example of two images. One was taken by the left camera, the other by the right.
We'll represent the distances as a depth map. A depth map is an image whose pixel values encode the distance from the camera to the object visible at that pixel:
A depth map built from a stereo pair
Let's look at the principle behind recovering a depth map from a stereo pair. You can think of the stereo pair as two eyes. Human eyes are also able to judge distances to objects. Let's look at the diagram and try to picture how a person sees.
Schematic view of a stereo pair. A1 — left camera. A2 — right camera. D — object. T — distance between the left and right cameras. Z — distance from the object to the camera plane. X0 — the intersection of the image plane with the line from the left camera to the object. (X0 − d) — the same for the right camera.
Here we have point D — the point whose distance we want to find. Point A1 is the left eye, A2 is the right one. For simplicity, imagine there's a pane of glass in front of our eyes with the image projected onto it. Point X0 is where the glass intersects the line from the left eye to the point, and point X0 − d is where the glass intersects the line from the right eye to the point. Here d is the shift of the object on the right image relative to the left one.
We need to find the distance Z. Let's apply the proportionality relation for the triangles:
The sides and heights of the red and green triangles are proportional.
We get the following equation:
Solving for Z:
The problem is that we somehow need to match what the left eye sees against what the right eye sees, so we can tell they're looking at the same point. In other words, we need to find the points X0 and X0 − d. We'll assume these points lie at the same height. We know the focal length and the distance between the cameras. So the task boils down to this: for every pixel in the left image, find d — how far that pixel is shifted in the right image.
You can read more about stereo vision in this article
Stereo Block Matching
There are several algorithms for solving this problem, and one of them is Stereo Block Matching.
It's based on the idea of finding similar neighborhoods around specific pixels. We'll consider the neighborhoods around two pixels as squares.
For each pixel, we'll go over all possible shifts, compute the sum of squared differences between the left and right images, and pick the square with the smallest sum. The shift of that most-similar square (d) reflects the relative distance to the object.
Pixel neighborhoods in the left- and right-camera images. The number is the light intensity at that pixel. As you can see, the intensities in this example are almost identical.
In this example we can see the two squares don't differ much. To measure how similar they are, we build a new square holding the element-wise difference. The smaller that difference, the better.
Element-wise difference of the pixel values
Once we have the square representing the difference between the two, we can sum up its elements:
0 + 1 + 0 + 1 + 1 + 0 + 2 + 0 + 0 = 5
For a given pixel X0 in the left image, we iterate over different values of d, which correspond to different pixels in the right image. We pick the d that yields the minimum sum, and that will be the distance to the object.
It's worth noting that besides the element-wise difference, there are many other cost functions. But we'll stick with this one.
The naive solution
In the naive solution we just do exactly what the algorithm describes, head-on. We go over every pixel in the left image, trying shifts from 0 up to some number, say 1000. For each pixel in the left image we walk over all the pixels in its neighborhood and, at the same time, walk over the pixels in the right image accounting for the shift. Then we sum up their differences. This is called the SAD (Sum of Absolute Differences) metric:
Remember that we can compute d once we can compute sad — for that we take the d that minimizes sad. And that d is the distance to the object.
So the formula for the distance to the object looks like this:
The disparity for x, y is the d that minimizes sad(x, y, d)
Let's write out the algorithm. For simplicity our squares will be of size 3.
Pseudocode
The code looks like this:
def stereo_block_matching(left_image, right_image):
height = len(left_image)
width = len(left_image[0])
# In the resulting image, pixel values encode the distance to the object
# (the brighter, the closer)
result = [[0 for x in range(width)] for y in range(height)]
for y in range(height):
for x in range(width):
best_disparity = None
best_difference = None
for disparity in range(1000):
difference = 0
for shift_y in range(-1, 2):
for shift_x in range(-1, 2):
difference += abs(left_image[y][x] - right_image[y + shift_y][x + shift_x])
if best_difference is None or difference < best_difference:
best_difference = difference
best_disparity = disparity
result[y][x] = best_disparity
return result
Let's look at the complexity of this algorithm. Assuming the image is square with side N, we search shifts from 0 to d, and the compared block is a square of size K, the total complexity is:
The prefix sums solution
If K equals 3, the complexity doesn't look too critical overall. But K is usually larger, since a bigger neighborhood accounts for more pixels and gives higher accuracy. Of course, if you make K too large, the algorithm gets less accurate again. A size of 21, for example, is commonly used. In that case K × K is more than four hundred. If we could get rid of that K × K factor in the asymptotics, we'd save a lot.
The prefix sums method can help us here. The idea is to reuse a sum we've already computed. Right now we compute the sum of differences over a pixel neighborhood, then move on to compute it for the next shift and recompute it from scratch, even though we're only adding a new column and dropping an old one. Instead of recomputing the whole sum, we can recompute only the new and old columns. The idea is the same one used for a box filter. The same thing happens when we move to the next row. We lose part of the already-computed sums; all we need is to add a new row and drop the old one, yet we recompute everything. But we can recompute just the new and old rows.
Column SAD
Now for a more detailed description. Let's introduce a variable column_sad — the summed pixel difference over a portion of a column.
Let's look at an example value of column_sad[x, y]:
Computing column_sad: take the columns in the left and right images
Computing column_sad: the column of pixel differences between the left and right images
In this example column_sad[x, y] equals 1 + 1 + 0 = 2.
The formula for column_sad looks like this:
The column_sad value
Computing column_sad quickly
This variable can be recomputed in constant time using its value from the previous row.
Here's an example of the fast update:
Add the element in the last row, subtract the one from the previous row
Add the difference in the last row, subtract it from the previous row
In this example column_sad[x, y] equals (1 + 1 + 0) + 8 − 1 = 9.
The same column_sad value, only computed faster
Computing SAD from column_sad
Using column_sad for a column, we can easily compute the SAD for the whole square:
Moreover, this sum can also be derived from the previous value:
Here's an example of the fast SAD update:
Add the last column, subtract the previous one
Add the difference in the last column, subtract the difference in the previous column
So, sad[x, y] = ((0 + 1 + 2) + (1 + 1 + 0) + (0 + 0 + 0)) + (3 + 2 + 4) − (0 + 1 + 2) = (1 + 1 + 0) + (0 + 0 + 0) + (3 + 2 + 4) = 2 + 0 + 9 = 11
By first computing the column_sad values for each y, and then, given those, computing sad for the current pixel and shift, we get the following complexity:
You can also notice that the sad and column_sad values for different shifts are independent of each other, so they can be computed separately and, accordingly, we only need to store w memory for column_sad.
The algorithm's code
from PIL import Image
def stereo_block_matching(left_image, right_image):
height = len(left_image)
width = len(left_image[0])
# In the resulting image, pixel values encode the distance to the object
# (the brighter, the closer)
result = [[0 for x in range(width)] for y in range(height)]
column_sad = [[[]]]
sad = [[[]]]
block_size = 21
s = block_size // 2
for y in range(height):
for x in range(width):
best_disparity = None
best_difference = None
for d in range(100):
sad[y][x][d] = column_sad[y][x + s][d] \
+ sad[y][x - 1][d] \
+ column_sad[y][x - s - 1][d]
column_sad[y][x][d] = abs(left_image[y + s][x] - right_image[y + s][d]) \
+ column_sad[y - 1][x][d] \
- abs(left_image[y + s][x] - right_image[y + s][d])
if best_difference is None or sad[y][x][d] < best_difference:
best_difference = sad[y][x][d]
best_disparity = d
result[y][x] = best_disparity
return result
Conclusion
That's it. We've looked at applying the prefix sums method in a rather unusual role — estimating the distance from a camera to objects.
In one of the next posts I'll show how to implement this algorithm even more efficiently on a processor that supports SIMD and VLIW, with a neat trick.
Top comments (0)