scikit-image advanced panorama tutorial
Enhanced from the original demo as featured in the scikit-image paper.
Multiple overlapping images of the same scene, combined into a single image, can yield amazing results. This tutorial will illustrate how to accomplish panorama stitching using scikit-image, from loading the images to cleverly stitching them together.
First things first
Import NumPy and matplotlib, then define a utility function to compare multiple images
Load data
The ImageCollection
class provides an easy and efficient way to load and represent multiple images. Images in the ImageCollection
are not only read from disk when accessed.
Load a series of images into an ImageCollection
with a wildcard, as they share similar names.
Inspect these images using the convenience function compare()
defined earlier
Credit: Images of Private Arch and the trail to Delicate Arch in Arches National Park, USA, taken by Joshua D. Warner.
License: CC-BY 4.0
0. Pre-processing
This stage usually involves one or more of the following:
Resizing, often downscaling with fixed aspect ratio
Conversion to grayscale, as some feature descriptors are not defined for color images
Cropping to region(s) of interest
For convenience our example data is already resized smaller, and we won't bother cropping. However, they are presently in color so coversion to grayscale with skimage.color.rgb2gray
is appropriate.
1. Feature detection and matching
We need to estimate a projective transformation that relates these images together. The steps will be
Define one image as a target or destination image, which will remain anchored while the others are warped
Detect features in all three images
Match features from left and right images against the features in the center, anchored image.
In this three-shot series, the middle image pano1
is the logical anchor point.
We detect "Oriented FAST and rotated BRIEF" (ORB) features in both images.
Note: For efficiency, in this tutorial we're finding 800 keypoints. The results are good but small variations are expected. If you need a more robust estimate in practice, run multiple times and pick the best result or generate additional keypoints.
Match features from images 0 <-> 1 and 1 <-> 2.
Inspect these matched features side-by-side using the convenience function skimage.feature.plot_matches
.
Most of these line up similarly, but it isn't perfect. There are a number of obvious outliers or false matches.
Similar to above, decent signal but numerous false matches.
2. Transform estimation
To filter out the false matches, we apply RANdom SAmple Consensus (RANSAC), a powerful method of rejecting outliers available in skimage.transform.ransac
. The transformation is estimated using an iterative process based on randomly chosen subsets, finally selecting the model which corresponds best with the majority of matches.
We need to do this twice, once each for the transforms left -> center and right -> center.
The inliers
returned from RANSAC select the best subset of matches. How do they look?
Most of the false matches are rejected!
3. Warping
Next, we produce the panorama itself. We must warp, or transform, two of the three images so they will properly align with the stationary image.
Extent of output image
The first step is to find the shape of the output image to contain all three transformed images. To do this we consider the extents of all warped images.
Apply estimated transforms
Warp the images with skimage.transform.warp
according to the estimated models. A shift, or translation is needed to place as our middle image in the middle - it isn't truly stationary.
Values outside the input images are initially set to -1 to distinguish the "background", which is identified for later use.
Note: warp
takes the inverse mapping as an input.
Warp left panel into place
Warp right panel into place
Inspect the warped images:
4. Combining images the easy (and bad) way
This method simply
sums the warped images
tracks how many images overlapped to create each point
normalizes the result.
Finally, view the results!
What happened?! Why are there nasty dark lines at boundaries, and why does the middle look so blurry?
The lines are artifacts (boundary effect) from the warping method. When the image is warped with interpolation, edge pixels containing part image and part background combine these values. We would have bright lines if we'd chosen cval=2
in the warp
calls (try it!), but regardless of choice there will always be discontinuities.
...Unless you use order=0
in warp
, which is nearest neighbor. Then edges are perfect (try it!). But who wants to be limited to an inferior interpolation method?
Even then, it's blurry! Is there a better way?
5. Stitching images along a minimum-cost path
Let's step back a moment and consider: Is it even reasonable to blend pixels?
Take a look at a difference image, which is just one image subtracted from the other.
The surrounding flat gray is zero. A perfect overlap would show no structure!
Instead, the overlap region matches fairly well in the middle... but off to the sides where things start to look a little embossed, a simple average blurs the result. This caused the blurring in the previous, method (look again). Unfortunately, this is almost always the case for panoramas!
How can we fix this?
Let's attempt to find a vertical path through this difference image which stays as close to zero as possible. If we use that to build a mask, defining a transition between images, the result should appear seamless.
Seamless image stitching with Minimum-Cost Paths and skimage.graph
Among other things, skimage.graph
allows you to
start at any point on an array
find the path to any other point in the array
the path found minimizes the sum of values on the path.
The array is called a cost array, while the path found is a minimum-cost path or MCP.
To accomplish this we need
Starting and ending points for the path
A cost array (a modified difference image)
This method is so powerful that, with a carefully constructed cost array, the seed points are essentially irrelevant. It just works!
Define seed points
Construct cost array
This utility function exists to give a "cost break" for paths from the edge to the overlap region.
We will visually explore the results shortly. Examine the code later - for now, just use it.
Use this function to generate the cost array.
Allow the path to "slide" along top and bottom edges to the optimal horizontal position by setting top and bottom edges to zero cost.
Our cost array now looks like this
The tweak we made with generate_costs
is subtle but important. Can you see it?
Find the minimum-cost path (MCP)
Use skimage.graph.route_through_array
to find an optimal path through the cost array
Did it work?
That looks like a great seam to stitch these images together - the path looks very close to zero.
Irregularities
Due to the random element in the RANSAC transform estimation, everyone will have a slightly different blue path. Your path will look different from mine, and different from your neighbor's. That's expected! The awesome thing about MCP is that everyone just calculated the best possible path to stitch together their unique transforms!
Filling the mask
Turn that path into a mask, which will be 1 where we want the left image to show through and zero elsewhere. We need to fill the left side of the mask with ones over to our path.
Note: This is the inverse of NumPy masked array conventions (numpy.ma
), which specify a negative mask (mask == bad/missing) rather than a positive mask as used here (mask == good/selected).
Place the path into a new, empty array.
Ensure the path appears as expected
Label the various contiguous regions in the image using skimage.measure.label
Looks great!
Rinse and repeat
Apply the same principles to images 1 and 2: first, build the cost array
Add an additional constraint this time, to prevent this path crossing the prior one!
Check the result
Your results may look slightly different.
Compute the minimal cost path
Verify a reasonable result
Initialize the mask by placing the path in a new array
Fill the right side this time, again using skimage.measure.label
- the label of interest is 3
Final mask
The last mask for the middle image is one of exclusion - it will be displayed everywhere mask0
and mask2
are not.
Define a convenience function to place masks in alpha channels
Obtain final, alpha blended individual images and inspect them
What we have here is the world's most complicated and precisely-fitting jigsaw puzzle...
Plot all three together and view the results!
Fantastic! Without the black borders, you'd never know this was composed of separate images!
Bonus round: now, in color!
We converted to grayscale for ORB feature detection, back in the initial preprocessing steps. Since we stored our transforms and masks, adding color is straightforward!
Transform the colored images
Apply the custom alpha channel masks
View the result!
Save the combined, color panorama locally as './pano-advanced-output.png'
Once more, from the top
I hear what you're saying. "But Josh, those were too easy! The panoramas had too much overlap! Does this still work in the real world?"
Go back to the top. Under "Load Data" replace the string 'data/JDW_03*'
with 'data/JDW_9*'
, and re-run all of the cells in order.