Image reconstruction

Having described an image by a set of moments, it may prove useful to investigate which moments give rise to which characteristics of the image, or vice versa. This can be achieved by reconstructing the original image from the calculated moments. Moment reconstruction, for moments with orthogonal basis functions (such as Legendre and Zernike moments which will be discussed later) has been developed extensively, [11,12,18,19]. However, where the basis set is non-orthogonal (such as Cartesian and centralised moments), only one method has appeared (although, non-orthogonal transform methods exist). This is the method of moment matching for non-orthogonal moment reconstruction [18]. The method is based upon creating a continuous function that has identical moments to that of the original function. In this section it has been applied first to Cartesian moments and then to the centralised moments. It must be noted that in applying the theory to sampled images, the continuous conditions are replaced by discrete versions, reducing the accuracy of the final function.

Assuming that all moments $M_{pq}$ of a function $f(x,y)$ and of order $N=(p+q)$ are known from zero through to order $N_{max}$, it is then possible to obtain the continuous function $g(x,y)$ whose moments match those of the original function $f(x,y)$, up to order $N_{max}$. (With reference to the Taylor series expansion in Section 1.1), assuming that the given continuous function can be defined as:

$$g(x,y) = g_{00} + g_{10}x + g_{01}y + g_{20}x^2 + g_{11}xy +~... g_{pq}x^py^q$$ (22)

which reduces to:

$$g(x,y) = \sum_{p=0} \sum_{q=0} g_{pq}x^py^q ~~~;~~~N_{max}=p+q$$ (23)

then the constant coefficients $g_{pq}$, are calculated, so that the moments of $g(x,y)$ match those of $f(x,y)$, assuming that the image is a continuous function bounded by:

$$x~\epsilon~[-1,1]~~,~~y~\epsilon~[-1,1]$$ (24)

These limits can be achieved by normalising the pixel range over which the Cartesian moments are calculated, thus:

$$\int_{-1}^{1}\int_{-1}^{1}g(x,y)x^py^q~dx~dy \equiv M_{pq}$$ (25)

Substituting Equation 1.22 into Equation 1.25 and then solving the integration produces a set of Linear Equations (LE), the number of which is determined by the order $(p+q)$ of reconstruction. These can then be solved for the coefficients $g_{pq}$ (in terms of the moments $M_{pq}$) by using matrix inversion. For order three ($(p+q)\leq 3$), the LEs in matrix form are:

$$\left[ \begin{array}{*{3}{c}} \hspace*{0.5cm} & \hspace{0.5... ...{*{1}{c}} M_{00} \\ M_{20} \\ M_{02} \\ \end{array} \right]$$ (26)

$$\left[ \begin{array}{*{3}{c}} \hspace*{0.5cm} & \hspace{0.5... ...{*{1}{c}} M_{10} \\ M_{30} \\ M_{12} \\ \end{array} \right]$$ (27)

$$\left[ \begin{array}{*{3}{c}} \hspace*{0.5cm} & \hspace{0.5... ...{*{1}{c}} M_{01} \\ M_{03} \\ M_{21} \\ \end{array} \right]$$ (28)

and finally:

$$g_{11} = \frac{9}{4}M_{11}$$ (29)

Applying matrix inversion to the first matrix, Equation 1.26 produces:

$$\left[ \begin{array}{*{3}{c}} \hspace*{0.5cm} & \hspace{0.5... ...{*{1}{c}} g_{00} \\ g_{20} \\ g_{02} \\ \end{array} \right]$$ (30)

By repeating this for all the matrices, it is possible to calculate all the coefficients. If they are then substituted back into Equation 1.22 an expression for $g(x,y)$ is produced. This expression can then used to reconstruct an approximation of the original image. The reconstruction function $g(x,y)$ is now in terms of weighted sums of the moments $M_{pq}$, which have been previously calculated from the original function $f(x,y)$. The resultant function $g(x,y)$ for order three is:

$$16g(x,y) = (14M_{00} -15M_{20} - 15M_{02})$$

$$+ \mbox{} (90M_{10} - 105m_{30} - 45M_{12})x$$

$$+ \mbox{} (90M_{01} - 105M_{03} - 45M_{21})y$$ (31)

$$+ \mbox{} (-15M_{00} + 45M_{20})x^2 + 36M_{11}xy$$

$$+ \mbox{} (-15M_{00} + 45M_{02})y^2 + (-105M_{10} + 175M_{30})x^3$$

$$+ \mbox{} (-45M_{01} + 135M_{21})x^2y + (-45M_{10} + 135M_{12})xy^2$$

$$+ \mbox{} (-105M_{01} + 175M_{03})y^3$$

Implementing this method to order $(p+q)=8~$ for binary images of simple shapes produces recognisable results, as shown in Figures 1.3 and 1.4. Figure 1.3a is the original image from which the moments were calculated and Figure 1.3b is the image reconstructed from the moments. The borders of the shape appear unclear, but they appear when the reconstructed image is thresholded, Figure 1.3c. Here the level of the applied threshold was adjusted by visual comparison with the original image. Due to the nature of the continuous function, the final shape is dependent on the threshold level, as is apparent in Figure 1.3c, as compared with the original image, Figure 1.3a.

Figure 1.3: Order $8$ Cartesian reconstruction of an ellipse.

$\parbox{4.8cm}{ \center{ \fbox{\rotatebox{-0}{\scalebox{0.25}{\includegraphics{images/theory/original_ellipse.ps}}}} } }$$\parbox{4.8cm}{ \center{ \fbox{\rotatebox{-0}{\scalebox{0.25}{\includegraphics{images/theory/reconst_ellipse.ps}}}} } }$$\parbox{4.8cm}{ \center{ \fbox{\rotatebox{-0}{\scalebox{0.25}{\includegraphics{images/theory/threshold_ellipse_reconst.ps}}}} } }$
$\parbox{4.8cm}{\center (a) Original Image}$$\parbox{4.8cm}{\center (b) Reconstructed continuous function}$$\parbox{4.8cm}{\center (c) Thresholded}$

This analysis is then repeated for the rectangle in Figure 1.4. The corners of the rectangle in Figure 1.4c are missing. The corners represent the high frequency content in the image, thus will be described fully by higher order moments. So the thresholded shape will converge to the original shape as the number of moments (and thus the order) increases. However for more complex shapes, higher accuracy $(p+q)\gg8$ is needed. This is analogous to the high frequency information needed to reconstruct pulsed time domain waveforms, using methods like Fourier series. As the order (and accuracy) increases, so does the number of LEs that need to be solved (reconstruction for order eight resulted in forty five LE's). Further, if it was required to increase the order of reconstruction (using Equation 1.31), then all coefficients need to be re-calculated. This is due to the correlated nature of the Cartesian moments, each moment does not simply provide its own individual contribution, (unlike the orthogonal case which will be discussed later in this chapter). It is interesting to note the effects of the Gibbs phenomena [16] which are more evident in the reconstructed ellipse - Figure 1.3b. The Gibbs phenomena (explained in terms of Fourier series) is the inability for a continuous function to recreate a step function - no matter how many finite high order terms are used, an overshoot of the function will occur. Here the discontinuous edge of the original intensity function of the ellipse (between the ellipse and the background) appear unclear in the reconstruction. While outside of the original area of the ellipse, `ripples' of overshoot of the continuous function are visible.

Figure 1.4: Order $8$ Cartesian reconstruction of a rectangle.

$\parbox{4.8cm}{ \center{ \fbox{\rotatebox{-0}{\scalebox{0.25}{\includegraphics{images/theory/original_square.ps}}}} } }$$\parbox{4.8cm}{ \center{ \fbox{\rotatebox{-0}{\scalebox{0.25}{\includegraphics{images/theory/reconst_square2.ps}}}} } }$$\parbox{4.8cm}{ \center{ \fbox{\rotatebox{-0}{\scalebox{0.25}{\includegraphics{images/theory/threshold_square_reconst.ps}}}} } }$
$\parbox{4.8cm}{\center (a) Original Image}$$\parbox{4.8cm}{\center (b) Reconstructed continuous function}$$\parbox{4.8cm}{\center (c) Thresholded}$

By assuming the same constraints as for Cartesian moment matching, the theory can be extended to centralised moments. The continuous function $g(x,y)$ is now defined as:

$$g(x,y) = \sum_{p=0} \sum_{q=0} g_{pq}(x-\overline{x})^p(y-\overline{y})^q ~~~;~~~N_{max}=p+q$$ (32)

similarly, Equation 1.25 becomes:

$$\int_{-1}^{1}\int_{-1}^{1}g(x,y)(x-\overline{x})^p(y-\overline{y})^q~dx~dy \equiv M_{pq}$$ (33)

where $\overline{x}$ and $\overline{y}$ are the $x$ and $y$ COM's, respectively. Solving for $g(x,y)$ is then achieved in the same manner as already described for the Cartesian case.