Interactive Image Reconstruction from its Incomplete, Irregular and Imprecise Fragments

The aim of this paper is to present a novel method as an adjunct to a human/expert for assembling unknown images from their fragments. Reassembling an image from its partial and eroded or decayed constituent fragments is of great importance in variety of fields, such as forensics, failure analysis, anthropology, archaeology and art reconstruction. It is a common task, indeed, in archaeological research and artefact preservation. Reassembly from archaeological fragments is a much more involved problem in comparison to other fields due to the following reasons. The number of randomly mixed fragments is normally huge, the fragments are mostly corrupted, irregular, and frequently with uncertain and eroded boundaries. In archaeology the process of assembling may not be required to produce a complete solution. Indeed, it may be required to produce a series of possible options and even partial solutions which the archaeologist would appraise subject to a range of archaeological and often unquantifiable constraints. We propose such a method for archaeological images. This paper is organized as follows. Prior work related with our proposed method is shared in Section 2, while detailed explanation of the proposed method is presented in Section 3. Finally, this paper ends with the conclusion part in the Section 4.


Introduction
The aim of this paper is to present a novel method as an adjunct to a human/expert for assembling unknown images from their fragments. Reassembling an image from its partial and eroded or decayed constituent fragments is of great importance in variety of fields, such as forensics, failure analysis, anthropology, archaeology and art reconstruction. It is a common task, indeed, in archaeological research and artefact preservation. Reassembly from archaeological fragments is a much more involved problem in comparison to other fields due to the following reasons. The number of randomly mixed fragments is normally huge, the fragments are mostly corrupted, irregular, and frequently with uncertain and eroded boundaries. In archaeology the process of assembling may not be required to produce a complete solution. Indeed, it may be required to produce a series of possible options and even partial solutions which the archaeologist would appraise subject to a range of archaeological and often unquantifiable constraints. We propose such a method for archaeological images.
This paper is organized as follows. Prior work related with our proposed method is shared in Section 2, while detailed explanation of the proposed method is presented in Section 3. Finally, this paper ends with the conclusion part in the Section 4.

Related Work
The fact that reassembly of an unknown image manually from its irregular and very often eroded fragments is very time-consuming task, motivates researchers to design automated systems [1], [2], [3], [4]. The algorithms produced are usually built upon several subsystems which generate reliable or accurate scores between fragments. Although, these algorithms show great success in specific cases, they have, so far, limited success in reassembling incomplete images with vast numbers of eroded pieces automatically.
For 2D matching the algorithms proposed so far can be divided into two categories: a) colour based methods, and b) geometry based approaches. Colour based methods use the colour information present in the fragments and mainly at the boundaries [4], [5], [6], [7]. These algorithms are mostly efficient, but they fail when only the mean of the colour intensities belonging to a fragment is taken, or only when the boundary of the pair of fragments is considered. Geometry based approaches assemble the fragments by matching the boundary curves of the fragment pairs [8], [7], [9], [10]. The geometry based approaches can accurately match the adjacent pairs if there is no corrupted or eroded geometry at the broken parts of the fragments. Since the fragments of interest in our case have exposed to the natural elements, perhaps for thousands of years, it is inevitable that the edges of the pieces are eroded/corrupted, and their geometry is thereby destroyed. Moreover, the geometry based approaches are relatively slow as compared to other methods and they can fail in their search to lock at a local minimum of the objective function.
The latest research published focussed on initially finding matches between the fragments by checking the patterns at the boundaries of adjacent fragments. This pattern alignment is enhanced by several methods which can be listed briefly as follows, a) by using the colour information [11], b) by describing fractured edges using the new 'ribbon' idea [12], c) by normal maps [13]. These methods have great success when the features of the fragment pairs are significantly different. We construct an integrated approach with the above-mentioned methods, in which mostly the boundary correlations between the potential pairs dominate the matches. Matching scores are derived and the scores between the potential pairs are adjusted with each of the several algorithmic iterations in a graph based approach.

Proposed Method
Our proposed method is not meant to replace the experts' (in our case the archaeologist's) skills and knowledge but rather to be an aid and an additional tool. This aspect is incorporated in the structure of the algorithm. In our approach we support the most likely match with periodicity check, line continuity and inner layer correlation all of these terms are explained in this paper below. At the end, the most meaningful match is selected as decided by the expert/archaeologist. The method is explained in the following subsections which cover formulation of the problem, border and boundary extraction, score assignment, decision unit and score rearrangement.

Formulation of the Problem
Given the set of fragments G = G i of the image (I), we wish to find the most likely match between fragments. In the composition of the image I each fragment is identified by its pixel boundary vector B i , i = 1...N, where N is the total number of available fragments. The image (I), has more fragments than N, such that some of the fragments are lost, however, only N number of fragments are accessible. Because of the possibly of imprecise, eroded or corrupted boundaries we also check the inner layers. A layer is a sequence of pixel values once or several times removed off the boundary of a given fragment. Hence, each fragment is defined by its boundary matrix for each layer, B ij , where j is the layer number 1…ly, and ly is the maximum layer number. To locate the maximum correlation exactly in the whole boundary matrix of the fragment we also detect the corners (C) of a fragment and divide the boundary vector into borders (sub-boundaries) F k , where k = 1...C, F k �B ij � G i �I, as shown in the Figure 1.
Another reason for dividing the boundary matrix into borders is the fact that each boundary matrix is huge compared to the likely correct matched pairs. This result in vast and unrelated parts of the boundary matrix. These unrelated parts decrease the correlation between correctly matched pairs and give misleading results subsequently.
To avoid these problems, we prefer to check the relatively small but reasonable portions of the boundary vector/matrix borders as explained in the following subsection.

Border and Boundary Extraction
As already mentioned, each fragment is defined by its boundary and related borders for the layer of interest. To extract the boundary of the image we transform the image into binary (for instance the fragment is all white and the background is black), by knowing that the photographic image of each fragment is acquired at a constant colour background. The structuring element disk with the distance from its origin (r) proportional to the intended layer number is used for morphological erosion operation. We simply take the difference of the binary images after two consecutive erosion operations, with r set equal to the desired layer number, ly and ly -1. Then we obtain the raw boundary vector, without corner indicators. Commonly used and well-known functions for corner detection work well if the fragment is close to rectangular in shape with corners at approximately right angles. However, when we have a fragment with other form as in Figure 2a these methods fail and produce false corner alarms, such as in Figure 2b. This is essentially due the fact that we are working in the discrete domain where there are no continuous lines. As a result, when we zoom in the one portion of the fragment as in Figure 2b. Hence for the corner detection purposes we designed our own specific algorithm. We obtain a slope plot by using the boundary coordinates as in Figure 3. Corners are checked both from periodicity of the slope vector in Figure 3 within the boundary vector, and maximum-minimum slopes in a local area of corner.
In digital world, a line can be defined by its slope period shorter than its size, meaning that each line has a periodic behaviour that determines its starting point to end point slope. The line cannot follow same the same slope at each consecutive point as it can be seen at Figure 2(b), however it follows a period throughout the line. The proposed corner detection algorithm is designed based on this claim. If we assume that a boundary vector is taken, and its slope between two cascaded points are calculated, this brings a slope vector (s) of that boundary as it can be seen at Figure 3. Initially, enough data points (m1) depending on the least border size that is investigating, are taken from the boundary slope vector. After that, period of this data array (s(1)…s(m1)) is calculated, the calculated period of this line is less than its size (m1). Later on, the data is accumulated from "s" vector one at a time, and period of the accumulated data array (s(1)…s(m1), s(m1+1), s(m1+2)…s(m2)) is calculated. Moreover, if the calculated period is still less than the size of the array (m2), then this indicates that "s(m2)", belongs to the line. On the other hand, if the period is equal to the size of the array (m2) than this indicates that the line might be following another direction, starting from m2. This doubt about m2, that whether it is a corner point, or not, is overcome by accumulating the rest of the s vector size at a time, and its corresponding period, as it is done until m2 point. Furthermore, if they are equal in every time, then this indicates that slope vector is aperiodic after m2, and it is unique. After m2 is marked as a corner, the algorithm is iteratively computed for the s' vector, which corresponds to [s(m1+1), s(m1+2)….s(end)].
Within the local area next maximum-minimum and the previous maximumminimum slope values are checked and the detected corner is verified. This approach is also supported by the local standard deviation calculations of the differences at boundary's x values and y values. Peaks in the sliding standard deviation vector over the slope vector indicate a possible corner. These two checks, supports the proposed corner detection algorithm for the worst cases that it might fail.
Corner extraction results for an ideal fragment is shown in Figure 4(a), and corners for a fragment with curves can be found at Figure 4(b). We put an indicator (NaN in our case) for the detected corners in the boundary vector. Hence, we obtain a single vector to represent each G i .

Score Assignment
We used a dynamic graph based approach to track the relations between fragment pairs, and their corresponding borders. We basically use a cross correlation metric to determine the weights in the dynamic score graph.
Cross correlation calculated from the equation 1, where c y1y2 is the sample covariance function calculated from 2 and s y1 , s y2 are corresponding sample standard deviations of the time series y 1t and y 2t at lags k = 0, ±1, ±2, … calculated from equation 3. µ y1 and µ y2 are the sample means of the time series. Initially, the scores are calculated over the first layer. Elimination with a predetermined threshold is performed and most likely solutions are gathered. With this knowledge of the most possible matches, we supervise the inner layer correlations. However, the correlations at the inner layers are relatively small as it can be seen at the Figure 5. As a result, the threshold that eliminates the wrong matches and indicates frame borders of the image, (at the frame borders of the overall image I, we do not expect any match) is decreased. The whole idea is represented in the Figure 6.  At the inner layers the size of the corresponding border vector is smaller (fewer samples). This reduction in the length of the inner layers (pixel changes in the layer 1 to layer 3 cross correlation calculation) results in changes of the maximum correlation location within the two corresponding borders. To overcome this problem, the smaller border vector is linearly extrapolated |layer1 -layer2| times. Here, layer1 is the layer number of the first border and layer2 is the layer number of the second border.

Calculation of the Cross Correlation in 3D
Cross correlation is calculated between each fragment pairs at the desired layer (j). Especially, in each fragment pair (a1 and a2), all F k s of the related Ba 1 j and Ba 2 j. However, since RGB colour space intensities are used throughout the work presented, three dimensional vectors are obtained for the boundary vectors. The cross correlations between the corresponding channels of the border vectors are calculated. (red channel of a1 -red channel of a2) Then, each cross-correlation vector for red, green and blue for the related pair at the predetermined layer are multiplied and to yield a score vector. Since, the orientations of the fragments are unknown, the cross correlations should be controlled in two directions. First, the cross correlation is calculated in the same manner as the borders located in the boundary vector, and then one of them is reversed and then the cross-correlation calculation is performed again.

Storage of the Data
Five-dimensional arrays are used to store the scores, and for the graph theory applications, implication tables are constructed. The 3D illustration of the array is shared in Figure 7. In the score array, 2nd dimension is used to construct the implication table. In the 3rd dimension, the cross-correlation multiplication in R-G-B channels vector (score vectors) are stored. The 4th dimension is used to store the score vectors layer by layer. To be more specific, if ly is the last layer that is going to be used, then 4th dimension simply constructs the score vectors for the following combinations: 11, 12, 13… 1ly; 21, 22,...,2ly;ly1,ly2,...lyly. As mentioned earlier, the score vector is calculated first with the same locations of the borders vectors as they are in their boundary vectors and secondly, with one of the border vector is flipped. Hence, another 5th dimension is created for the same implication table, scores vector and layers in which one of the border vector is reversed. The periodicity of lines detection is also used as a metric which strengthens the score between two fragments. This signals that there is a tight correlation between the boundaries of two fragments. The detection is accomplished by checking the scores vector (product of cross correlations in R-G-B channels). When there are symmetrical peaks around the maximum point, as in Figure 8(a), it can be deduced that precise shifts (amount of the period) of the vectors results in again high correlation, which can be seen at the Figure 8(b). Then, after periodicity detection, the score between the related fragments is increased exponentially.

Decision Unit and Score Rearrangement
The score assignment initially starts automatically. First, the orientation of the vectors of interest are decided by checking the maximum values of the scores vector calculates between the 1st layers of the fragments. The correct orientation (whether the vector is reversed in the boundary vector, or not) although the fragments are not matched at all, yields higher acme among the score vector than the uncorrected orientation of the two border vectors. Subsequently, we also encode the correct orientation of the vectors to the scores array.
The designed algorithm iteratively goes over the layers up to the predetermined ly (maximum number of layers to be used) as it is illustrated in the Figure 8 for the most likely matches. The peak value of the score matrix in each layer is checked by the threshold and the shape of the score matrix is also considered as a decision metric. Mostly for the correct matches, the score vector shows an impulsive behaviour around the principally correlated locations of the border vectors. Moreover, at the inner layers of the score vector, it is expected that the impulsive behaviour will continue with a lower peak, but with the almost same principal location as it can be seen at Figure 5. However, for the wrong matches this impulsive behaviour is not observed.
The algorithm presents the most likely matches to the expert/archaeologist after the steps shown in Figure 6, starting from the furthest correlated fragments. Knowledge of the context, such as the archaeological excavation environment, historical cognisance and human intelligence here take place and decide whether the proposed match is correct or not. With this decision the scores are re-adapted and if the given match by the algorithm is rejected the corresponding score becomes NaN, and not used again.
(a) (b) Figure 8: Plot of the score vector (a) and the periodic lines that we detected by our algorithm(b)

Conclusion and Future Work
In this paper we presented a novel interactive image reconstruction algorithm. Our method consists mainly 3 parts. First, we extract the boundary vector and the corners of the fragments. Secondly, we construct a graph based on the correlation between the inner boundaries of the fragments. Afterwards, the expert/archaeologist is introduced to contribute to the process and to control the proposed algorithmic matches. Our results from initial trials show that the correct matches are determined.
Although our study is only a small step, we believe that this model will insert another insight, and will be useful for the construction of future image reassembly algorithms. In the future, the 3D space orientation of the fragments can be added as a metric for the matches of fractured pairs. Moreover, this method can be applied to the complex data-sets with imprecise boundaries and locally corrupted 3D fragments.