Seam line determination is important step of remote sensing image mosaicking. Automatic seam line selection problem using graph shortest path nding is discussed in this paper. Approaches to graph construction from image similarity cost function, graph memory layout and parallel seam line optimization using shortest path on graph are presented.
Acquisition of holistic image of a certain territory in the form of a mosaic from various (in the sense of geometry and brightness-contrast) separate images is an important problem when processing satellite and aerial photos. Wide range of remote sensing equipment with unique spectral sensitivity, variation of angle of view, position of the Sun and atmosphere condition, nonsimultaneous image acquisition cause unique deformation of the image of object of interest. Such a large number of negative factors means that it is practically impossible to compose the images, both from a geometric point of view, and from the point of view of radiometric conditions. One of the solutions of the described problem when seamless mosaic is created is optimal seam line detection within overlapping image regions.
To obtain a visually qualitative result that does
not introduce distortions into the information, strict
and, in some sense, contradictory restrictions on the
behavior of the seam line are superimposed in the
construction of photographic plans and mosaic of
images: a seam line must not cross high-altitude objects,
seam line must cross linear objects perpendicularly
and should be inside low-informativeness areas [
There are some di erent approaches for seam line creation:
Draw seam line along image borders (no optimization); Central image area maximization; Draw seam line inside areas with maximum similarity.
Draw seam line along image borders is used when there
are no quality demands to resulting image mosaic, or
there are strict time limitations. Since brightness and
geometrical distortions of remote sensing images are
minimal in their central areas we can draw seam line
as close to image center as possible using Voronoi
polygons [
Optimal seam line detection methods are based on
the search of the pixels with minimal di erences in
brightness, gradient and some other parameters which
are calculated from the image overlap regions
(intersections of image pairs in coordinate system of an
output image). Seam line being determined is a base
for united mosaic from several initial images. In this
case the visual divergence is signi cantly eliminated.
Many approaches bring optimal seam line
determination problem to energy function minimization. In
some articles normalized cross correlation (NCC) [
Even though in practice, graph cuts algorithms
have computational complexity close to linear [
Due to the fact that the solution of optimization
problems on graphs is associated with high
computational complexity, a number of authors try to simplify
the structure of the graph, leading to a reduction in
the time of solving the problem at the cost of a slight
degradation in quality [
The main contribution of this paper is a new approach to seam line optimization with low memory consumption and a high degree of parallelism. 3
To achieve a visual invisibility of the seam line the value of the pixels of overlapped images must coincide at points with the same coordinates in geographic coordinate system (or coordinate system of output image for non georeferenced mosaic). Seam line optimization is one of steps of remote sensing image mosaicking. Normally, it follows after geometry alignment and radiometric correction of input images. We expect that these two steps have already been correctly done, but not ideally due to some reasons and we need to eliminate the remaining visual inconsistencies. 3.1
In this article to create optimal seam line we propose to build energy cost function comprised of three functions: the image similarity matrix Ws, the image informativeness matrix Wi and the forbidden zone matrix (linear extensive, high-altitude and other visually obvious foreground objects). The similarity matrix Ws is calculated as follows:
Ws(x; y) = (I0(x; y)
I1(x; y))2 (1) where I0(x; y) and I1(x; y) denote the brightness on the rst and second overlapping images in the (x; y) position (here and later we assume that (x; y) is a pixel coordinates (row/column) of input images in coordinate system of output mosaic and all geometry transformations, including nonlinear ones, already have been done at geometry align step).
To de ne informativeness matrix Wi Moravec [
In order for the seam line to cross the line objects along a line close to the perpendicular and not intersect the high-altitude objects, we introduce the matrix of forbidden zones. To create forbidden zone matrix Wr(x; y) we rasterize vector lines and polygons, which represent linear extensive and high-altitude objects, using specially de ned penalty coe cients for crossing those objects. Due to the fact that the seam line receives an additional penalty when crossing objects marked in matrix Wr(x; y), the optimization algorithm will try to bypass such objects or cross them along the line of the minimum possible length, that for thin elongated objects (roads, rivers) will be perpendicular.
Resultant energy function de ned as a weighted sum of the particular matrices:
W (x; y) =
Ws(x; y) +
Wi(x; y) +
Wr(x; y); (2) where , , are the weighting coe cients, which dened each matrix contribution. In most cases, we can choose all weights to be the same, but if one wants to amplify one of the constraints, the corresponding coe cient can be increased. For example, if the visual similarity of images is more important to us, then we increase , and if it is more important for us to draw the seam line in places with low informativeness, then we increase the coe cient, if we want to completely prohibit the intersection of linear and high-altitude objects, we can set the coe cient close to in nity.
Thereby, seam line, which was delineated along the path with minimal energy cost sum, meet the seam line principal criteria: do not intersect high-altitude objects, linear objects should be intersected along the shortest distance (perpendicularly), and the other object should be intersected in the areas with maximum similarity and minimum informativeness. For optimal seam line determination we de ne energy function as weighted graph. Each node of the graph represents corresponding energy function matrix element. Let us connect all graph nodes to their neighbors, here we have two options: four- and eight- connected graph (Figure 1). We de ne the weight of the edge between energy matrix cells (x1; y1) and (x2; y2) as the sum of weights in the corresponded nodes multiplied by the distance between them: [W (x1; y1) + W (x2; y2)]p(x1 Thereby energy function matrix of M N size is represented by at graph consisted of M N nodes and PconnM N edges (where Pconn is a connection amount of energy matrix elements). Such a construction of the graph allows us not to store the graph explicitly, but to build it as necessary from the energy matrix on the y.
b) Eight-connectivity where Se is the size of one matrix element. Majority of the Earth remote sensing photos are 16 bit per channel images, therefore it is reasonable to use similar approach to energy matrix storage data type (Se = 2). The amount of memory required to store the graph explicitly in the form of compressed sparse row matrix (CSR) is: Scsr = SP trM N + SP trPconnM N + SePconnM N; (5) where rst term is a memory amount to store graph nodes, second term is memory amount to store graph edge pointers and third term is memory for edge weights, SP tr is node/edge pointer size (usually 4 or 8 byte), Se is size of edge weight.
When we store and use the graph in the described implicit way, it is possible to reduce memory usage by SP tr+SP trPScoenn+SePconn times. When 8 byte pointers and eight- connected graph are used, it is possible to reduce the RAM usage by a factor of 44 in comparison with the CSR graph storage type.
With this approach to storing the graph and when combined with the block method of storing the energy matrix W (x; y) in memory, the processor cache is more e ciently used, since when a graph is traversed, most nodes connected by edges are in the adjacent memory cells.
Moreover, the block storage of the energy matrix W (x; y) allows us to calculate and load matrix blocks as we use at the optimization algorithm, which further reduces the algorithm's requirement for the amount of memory and allows us to process images larger than the available RAM. 3.3
Typical satellite image size is about 50000 100000 pixels. When seam line between two images is created, we deal with the graph of 109 nodes and 1010 edges.
To avoid processing the entire array of data, we apply a hierarchical approach to image processing.
In the rst step, we construct optimal seam line at the overview level of the original images (Figure 2). To do this, we calculate the coordinates of the points of intersection of the image frames in the coordinate system of the output image (Vb and Ve), read input images from the overview levels of detail, build a coarse energy matrix and run the optimization process between the graph nodes corresponding to the points of intersection of the image frames. As a result, we obtain a coarse approximation of the seam line (blue polyline with black vertices). In the case where the image frames are not convex polygons and the intersection result is a multipolygon, each part of the multipolygon is processed independently.
In the second step, we will re ne the seam line between the vertices of the seam line obtained on an overview scale. To do this, we construct the ne energy matrix using fragments of input images from the detailed level and perform the optimization process again, taking as the initial and nal nodes of the graph, those that correspond to the vertices of the coarse seam line. Note that to carry out the optimization, it is not necessary to build the entire energy matrix W (x; y), but only the part that is near the coarse approximaImage blocks that is being processed
Vb Ve Refinement area around overview seam-line Seam-line obtained on an overview scale Seam-line obtained on a detailed scale Seam-line refined on a detailed scale tion of the seam line. Moreover, each image fragment can be processed in parallel.
But the new seam lime will have fractures near the vertices of the coarse seam line, because end point vertices of seam lime, are rigidly xed. To eliminate unwanted fractures, we will carry out the third stage of the seam line re nement on a detailed scale. To do this, we will take as a initial and nal nodes of the graph, those that correspond to the midpoints of the segments, obtained in the second step, and we will perform optimization again. Thus, we obtain a seam line, unnoticeable both at the detailed and at the overview scale (Figure 2). 3.4
For optimization we use the shortest path Dijkstra`s
algorithm [
Data: W (x; y), Vb, Ve
Result: S - optimal seam line vertices 1 D:f ill(inf); P:f ill(inf); 2 D[Vb] 0; 3 QMLB:insert(0; Vb) 4 while QMLB not empty do 5 c q:extractM in() 6 for p = 0; p < Pconn; ++p do 7 n c + Ap 8 // Compute edge weight according to (3) 9 ! (W (cx; cy) + W (nx; ny))length(Ap) 10 d D(cx; cy) + ! 11 if d < D(nx; ny) then 12 if P (nx; ny) is not set then 13 q:insert(d; n) 14 else 15 q:decreaseKey(d; n) 16 end 17 D(nx; ny) d 18 P (nx; ny) p 19 end 20 end 21 end 22 S turnBack(P; Ve)
Algorithm 1: Pseudo code of the optimization
algorithm taking into account structure of the graph
Implementation details taking into account the
proposed structure of the graph are given in the Algorithm
1. The input data for the algorithm are energy
function matrix W (x; y) and Vb, Ve { start and end points.
Note that the coordinates of the corresponding
mosaic pixels are used as pointers to the nodes of the
graph. In addition to the input data, the algorithm
uses several important internal structures. D(x; y) is
distance the matrix from the start node to each other.
The matrix D(x; y) must have a data type that allows
storing the sums of the elements of the matrix W (x; y)
without over owing. If we assume that the size of
the element of matrix W (x; y) is equal to two bytes,
then the size of the elements of matrix D(x; y) must
be equal to four or eight bytes. A is a possible arc list,
for four-connectivity case A4 = [ ; !; "; #], and for
eight-connectivity case A8 = [ ; !; "; #; %; &; -; .],
according to Figure 1. P (x; y) is node's parent
matrix. Since the graph has an ordered structure and
each node has four or eight neighbors, we can use the
arc number in array A as a pointer to the adjacent
node. Thus, no more than one byte is needed to store
the element of the matrix P . QMLB is a multi-level
bucket data structure. According to the calculations
in the paper [
We divide the problem into a number of blocks, while
re ning the seam line on detailed level, this allow us
to process unlimited amount of data and perform
parallel optimization in each individual blocks of
information. OpenMP [
Result: S - optimal seam line vertites 2 O1; O2 readImageOverview(I1; I2); 3 Vb; Ve computeIntersectionP oints(I1; I2); 4 Wcoarse computeEnergy(O1; O2); 5 Scoarse shortestP ath(Wcoarse; Vb; Ve); 7 B splitP athT oBlocks(Scoarse); 9 for step = 2; step <= 3; ++step do 11 #pragma omp parallel for 12 for i = 0; i < B:size(); ++i do 13 box computeBoundingBox(Bi); 14 Vb; Ve computeEndP oints(Bi); 16 begin #pragma omp critical 17 F1; F2 readF ineImage(I1; I2; box); 18 end 19 Wfine computeEnergy(F1; F2); 20 Sfine;i shortestP ath(Wfine; Vb; Ve); 21 end 23 S constructP athF romSegments(Sfine); 24 if step == 3 then 25 break; 26 end 27 B splitP athT oBlocks(S); 28 end 29 return S;
Algorithm 2: Simpli ed pseudo code of parallel hierarchical seam line optimization 4
Several datasets were used to test hierarchical seam line optimization algorithm:
Colored (RGB) 8-bit per channel middle resolution images acquired by Landsat-8 satellite system 16000 16000 pixels each (Figure 3); Colored (RGB) 8-bit per channel high resolution images acquired by aerial camera 4368 2912 pixels each with roads digitized (Figure 5); Panchromatic 16-bit high resolution images acquired by Resurs-P satellite system 36000 105000 pixels size.
The rst and second datasets were used to study the visual quality of the algorithm, and the latter was used to estimate the performance.
Algorithm was tested on a workstation equipped with a quad-core Intel Core i7 processor and 8 GB RAM.
Figure 4 shows the output mosaic obtained using the hierarchical seam line optimization algorithm, it can be seen that the seam line is practically indistinguishable both on a coarse scale and on a ne scale. The result of the algorithm on high resolution images using a vector layer of roads as forbidden zones is shown in Figure 6. It can be seen that the seam line crosses the roads as perpendicular as possible. But additional restrictions on the seam line in the form of forbidden zones lead to the fact that near the zone, the visual similarity of images on the di erent sides of the line is lost.
The dependence of the time for optimal seam line constructing on the number of vertices in the graph is shown in the Figure 7.
The construction of the seam line for maximum size graph, with one-processor version with full data loading into memory was not performed due to high consumption of memory. As studies have shown, when constructing an eight-connected graph, the visual quality of the seam line is somewhat higher, because it includes diagonally oriented edges, but the processing time increases due to the increase in the number of edges in the graph and the non-integer weight of the diagonal edges. Also on the Figure 7b it is seen that for the size of graph 109 nodes the processing time sharply increases, which is due to the predominance of the time of data loading from the hard drive over the optimization time. 5
In this paper the approach to construct an optimal seam line for the Earth remote sensing image mosaic, based on the graph shortest path nding is presented. The construction of an energy function matrix taking into account both informativeness and similarity of input images is discussed. Graph data structure design from an energy matrix is presented, which makes it possible to reduce the amount of RAM. The hierarchical parallel implementation of the optimization algorithm is presented with regard to the proposed graph data structure. The amount of RAM required for the optimization process is estimated.
The proposed hierarchical image processing approach not only allows us to reduce a seam line construction time, but also to make seam line unnoticeable both at the detailed and at the overview scale of the output mosaic. Using proposed approach, the time for optimal seam line creation was reduced by more than 10 times, which make it possible to produce seamless image mosaic with minimal CPU time usage.