Geometry based three-dimensional image processing method for electronic cluster eye

Abstract

In recent years, much attention has been paid to the electronic cluster eye (eCley), a new type of artificial compound eyes, because of its small size, wide field of view (FOV) and sensitivity to motion objects. An eCley is composed of a certain number of optical channels organized as an array. Each optical channel spans a small and fixed field of view (FOV). To obtain a complete image with a full FOV, the images from all the optical channels are required to be fused together. The parallax from unparallel neighboring optical channels in eCley may lead to reconstructed image blurring and incorrectly estimated depth. To solve this problem, this paper proposes a geometry based three-dimensional image processing method (G3D) for eCley to obtain a complete focused image and dense depth map. In G3D, we derive the geometry relationship of optical channels in eCley to obtain the mathematical relation between the parallax and depth among unparallel neighboring optical channels. Based on the geometry relationship, all of the optical channels are used to estimate the depth map and reconstruct a focused image. Subsequently, by using an edge-aware interpolation method, we can further gain a sharply focused image and a depth map. The effectiveness of the proposed method is verified by the experimental results.

Keywords

Image processing electronic cluster eye artificial eyes mobile robotics geometric algorithms

1. Introduction

Modern robots are progressively becoming more realistic and similar to humans [23] being asked to solve tasks such as path planning [58], path following [18], trajectory tracking [49], reproducing the behavior of human limbs [24].

A fundamental role in robotics is played by Artificial Intelligence [55], Computational Modelling [32, 57], and Metaheuristic Optimisation. Algorithms belonging to the last category are usually inspired by nature. Some recently proposed examples include inspiration by music [41, 40], physics [43, 42, 44], swarm behavior [34, 37, 36]. Artificial Intelligence is often integrated within robot controllers [45, 28] and combined to neural systems for Machine Intelligence [19, 2, 13].

The eye of the robot is very important and complex device since it is one of the sensors that allow a robot to make decisions and thus being autonomous. Several problems arise in relation to the eyes of the robots. Some are related to design of the eye itself, e.g. retina modelling [25] while others are related to the interpretation of what the cameras incorporated the eyes capture. Some examples include image segmentation [3], feature extraction [47, 35], image recognition [7, 26], and image processing in medical applications [16, 17].

The present paper focuses on robots’ eyes and, more specifically on artificial compound eyes. Inspired by natural compound eyes, an artificial compound eye is manufactured to be an alternative imaging system to miniaturize an optical system [48]. Due to many merits, such as small size, wide field of view (FOV) and sensitivity to motion objects, much attention has been paid to artificial compound eye systems in recent years [51]. Inspired by a wasp parasite called Xenos Peckii, an electronic cluster eye (eCley) [5] was introduced by Brückner et al. in 2010 as a new type of artificial compound eyes. The details about the optical module, the working principle and the manufacturing process can be referred to [5]. The camera module of eCley is shown in Fig. 1 and its basic system parameters are listed in Table 1.

Table 1
Parameters of the eCley system in [5]

Property	Value
Size of eCley ( $L\times W\times H$ )	6.8 mm $\times$ 5.2 mm $\times$ 1.4 mm
FOV ( $H\times V$ )	58 ${}^{\circ}\times$ 46 ${}^{\circ}$
F-number	3.7
Focal length	0.778 mm
No. of channels	13 $\times$ 17
Pixel pitch	3.2 $\mu$ m
Diameter of microlens	0.375–0.386 mm
Angular shift	4 ${}^{\circ}$

Figure 1.

The eCley module (the one cent coin is used for the size comparison [5]).

An eCley is composed of a certain number of optical channels organized as an array. The array of microlenses is used for imaging. As shown in Fig. 2, each imaging channel can capture only a small part of the full FOV, so the FOV of each imaging channel is small and the FOVs of neighboring channels are partial overlapped. By using this structure, the total track length is only about 1.4 mm, which is much shorter than single-aperture camera [1]. In addition, since the focal length of eCley is only 0.778 mm, it has a merit that the eCley focuses on very closed objects easily. Therefore, eCley can be used for the navigation of micro-robot and the internal examination based on micro capsule robot. While the captured image is composed of a series of separated, small FOV and low-resolution images. It is hard for robots to make decision based on the small FOV and low-resolution image of a single channel or even several channels, thus the post-process to obtain a complete high-resolution image is required.

In the eCley system, a complete image with a full FOV is obtained by fusing the images from all the optical channels. The principal optic axes of optical channels are considered to be unparallel so as to gain a larger FOV than the parallel case on condition that the overall size of eCley is not changed. Each channel has its inherent offset against its neighboring channel. Thus, the parallaxes of objects at the same distance are not identical in different channels. The parallax from unparallel neighboring optical channels in eCley may lead to reconstructed image blurring and incorrectly estimated depth.

Figure 2.

(a) The captured image. (b) The extracted images after flat field correction. (c) Reconstructed image focused on a single distance. (d) The reconstructed image of our method. (e) The depth map of our method.

To solve this problem, a geometry based three-dimensional image processing method (G3D) for eCley is proposed to obtain a complete focused image and dense depth map. The main idea of G3D is the introduction of the geometry relationship of optical channels for eCley. Based on the geometry relationship, the mathematical relation between the parallax and depth among unparallel neighboring optical channels can be derived. Further, a sharply focused image and a depth map can be gained. Compared with the methods in [30, 31, 59], G3D has three main advantages:

G3D can reconstruct a complete all-in-focus image.

Based on the geometry relationship derived, all channel images can be directly used to estimate the depth of an object.

G3D can reach sub-pixel precision in principle, because the matching is based on the space coordinates, instead of pixel shift.

The rest of this paper is organized as follows. Section 2 reviews the related work. Section 3 presents the geometry based three-dimensional image processing method. Experimental results are shown in Section 4. Finally, conclusions are drawn in Section 5.

2. Related work

This study discusses the three-dimensional image processing method with respect to eCley covering object depth estimation and image reconstruction. So this section starts from brief introductions of object depth estimation and image reconstruction in the area of computer vision, and then moves to the survey of the related work on artificial compound eyes.

Object depth estimation is one of the most active research topics in computer vision, see [56, 6]. According to the different setting of system parameters, two different methods, depth from defocus (DFD) [52] and stereo matching [38, 11], are widely studied. DFD attempts to extract depth information from multiple images (at least two) captured at different focus setting [4], while the stereo matching method tries to obtain depth from multiple images captured at different positions. The performance analysis and comparisons between DFD and stereo matching algorithms was presented in [39]. Moreover, many public test images can be found in the web site by the link: http://vision.middlebury.edu/, and plenty of stereo matching and multi-view stereo algorithms have been listed and evaluated.

In the meantime, super-resolution (SR) image reconstruction is also a very attractive research area over the past two decades. An SR reconstruction algorithm is to obtain one or multiple high-resolution images from one or multiple low-resolution images. In general, SR algorithms consist of four main categories: nonuniform interpolation approach, frequency domain approach, regularized SR reconstruction approach and projection onto convex sets approach [33, 27].

In the area of artificial compound eyes, except for the system design, image reconstruction and depth estimation are two active research areas [51]. In SR reconstruction, Kitamura et al. proposed a pixel rearrange method to reconstruct the high-resolution image based on Thin Observation Module by Bound Optics (TOMBO) [21]. In [29], an iterative back projection method (IBP) was presented to further improve the reconstructed results. El-Sallam et al. introduced a spectral-based blind method for TOMBO [9, 10], but the point spread function (PSF) of sub-images should be known before. As for depth estimation, most algorithms were designed for TOMBO. In [53], a multiple baselines method was discussed. In [14], a sum of squared difference (SSD) evaluation based back-project method was described based on the fact that a back-project method will result in poor reconstruction performance when the object is located at the wrong position. Especially for eCley, a braiding algorithm was introduced in [30, 31] to obtain the reconstructed image with the final resolution 640 $*$ 480. The final image can only focus on a specified distance at a time. If the object is not in the specified distance, the image blurring appears. In [59], Ziegler et al. presented an algorithm based on stereo matching and fusion method to obtain the final depth map and complete image. The depth of each channel image was obtained by fusing the results of stereo matching of four neighbouring images. In [20], an approach based on the intensity transitional area was discussed. The approach can achieve subpixel precision, but the depth information can be obtained only at the object edge.

Figure 3.

Schematic diagram of G3D.

In summary, the work on image processing for eCley stays at the initial stage and lots of problems are required to be solved before it goes to practical applications.

3. G3D

To eliminate the reconstructed image blurring and incorrectly estimated depth resulting from the parallax coming from unparallel neighboring optical channels in eCley, this section discusses the proposed geometry based three-dimensional image processing method (G3D) for eCley. The main idea of G3D is that the geometry relationship of optical channels in eCley is derived to obtain the mathematical relation between the parallax and depth among unparallel neighboring optical channels. The geometry relationship is very useful for gaining a sharply focused image and a depth map.

The main steps of G3D are shown in Fig. 3. In the eCley system, the captured image goes through several processes: image extraction, flat field correction, spatial coordinate calculation, matching cost computation, matching cost aggregation, initial depth estimation, initial image reconstruction, and depth & image refinement. Finally, the dense depth map and focused image can be obtained. Each step will be described in detail as follows.

Figure 4.

The flat field correction of a channel image. (a) Before the flat field correction. (b) After the flat field correction.

3.1 Image extraction and flat field correction

An eCley system consists of multiple optical channels, but the images of all the channels are captured by one image sensor. The image of each channel occupies a small circular region on the image sensor. Assume that the eCley system has $m\times n$ channels, where $m$ and $n$ are odd numbers, and the resolution of the image sensor is $\textit{IW}\times\textit{IH}$ pixels. The central pixel $(\frac{IW}{2},\frac{IH}{2})$ of the central channel $(\frac{m+1}{2},\frac{n+1}{2})$ is located at the center of the image sensor. The intervals between images of adjacent channels are equivalent. In our eCley system, there are 13 $\times$ 17 channels, and the image captured by each channel is a circular region with a diameter of $79$ pixels. The total resolution of a raw image captured by the image sensor is 1536 $\times$ 2048 pixels. The interval between images of adjacent channels with regard to the central pixels is 111 pixels. Therefore, the channel image is extracted as a squared image centered at the central pixel of channel with a resolution of 55 $\times$ 55 pixels. Because the large oblique angle of incidence leads to high distortion, in this paper, only 13 $\times$ 13 channels are considered.

Due to the border effect, the amount of incident light at the pixels near the border of image channel is less than the rest pixels’, so the border of each image is darker than the rest of the image, which is shown on the left-up corner of the left image (pixels in the red box) in Fig. 4. The pixels located in the left-up corner of this red box are darker than the rest pixels in the red box. With the effect of intensity attenuation, the depth of each pixel can not be accurate estimate because the brightness consistency assumption does not hold and the influence is more severe than the random sampling noise and the system noise of eCley. To reduce these effects, the flat field correction [21, 50] is performed. The steps are as follows:

1.
Reference image capture. Two types of reference images are captured by using the eCley. One image is the white reference image, the other is the black reference. The black reference image is the Dark Current (DC) offset signal. To reduce the sampling noise, multiple reference images are captured continuously, in our case, 5 images of each reference type are enough to obtain a smooth reference image, then the average intensity of each pixel is calculated as the pixel intensity. $W(x,y)$ and $B(x,y)$ are the averaged responsivity of reference white image and black image at pixel $(x,y)$ , respectively.
2.
Pixel response calculation. The pixel responsivity of camera sensor $\Pr(x,y)$ is calculated according to equation $\Pr(x,y)=W(x,y)-B(x,y)$ . Due to the white reference image is the image of a uniform white object. The $\Pr(x,y)$ can reflect the distortion ratio of camera sensor. To further reduce the noise effect, a Least Square (LS) polynomial fitting method is used to fit the curve of pixel responsivity in each channel. The pixel responsivity after LS fitting is $P(x,y)$ .
3.
Flat field correction. For a captured object image $T(x,y)$ , the corrected image $T^{\prime}(x,y)$ is computed as follows:

$T^{\prime}(x,y)=\frac{T(x,y)-B(x,y)}{P(x,y)}$ (1)

The result of the flat field correction is shown in Fig. 4b. Comparing with the pixels in the red box of both left and right images, the light intensity of the border area is reduced in the original image after flat field correction. So the effect of intensity attenuation in the red box is reduced significantly.

Figure 5.
The simplified graph of the principal optic axis in eCley (only one dimension). The principal optic axes are not parallel with each other, and the angular deviation between adjacent channels is about $4^{\circ}$ .

3.2 Spatial coordinate calculation

This step is the key point of G3D and is also the foundation of the subsequent steps of G3D. The aim of spatial coordinate calculation is to handle the problem of nonuniform disparities at the same distance due to the nonuniform spatial sampling interval in the eCley system. Concretely, on the basis of several depth plane candidates, the pixel is back-projected to each depth plane to obtain its corresponding spatial coordinates, thus, several processes in the local stereo matching method [38], such as matching cost computation, matching cost aggregation and depth selection, can be used to estimate the correct depth by using the spatial coordinates.

Figure 6.

Interleaving pixels from neighboring channels in the eCley. Left: two marginal channels. Right: two central channels. For simplicity, each channel has nine pixels in one dimension. The central channels nearly have the uniform spatial sampling interval, while the marginal channels do not have. At the same distance, the central channels have exactly the same parallax for all pixels, but the parallax of the pixels in the marginal channels is not uniform.

Figure 7.

The simplified graph showing the relationship of camera system parameters (only in one dimension).

Differing from a single traditional lens camera, an array of microlenses is placed on the eCley, where each microlens has a small diameter, up to 400 $\mu$ m, and the spacing between neighboring microlens is less than 10 $\mu$ m. With such a small diameter, the FOV of a microlens is quite small. Thus, all microlenses are considered to jointly gain a large FOV. However, the principal optic axes of adjacent channels in eCley are not parallel and there is a certain degree of angular shift, as shown in Fig. 5. The large oblique angle of incidence results in the image astigmatism. To alleviate the astigmatism, some measures, such as the introduction of a well designed chirped array in [8], were studied. But the oblique incidence will lead to the nonuniform spatial sampling of the object. As shown in Fig. 6, the sampling interval of the central channel is nearly uniform at the same distance, but the sampling interval of the marginal channel with an oblique incidence is nonuniform and larger than the central channel’s. Even if the distances to the object are the same, the disparities of center image pair and marginal image pair are not identical. Consequently, the stereo matching method based on the disparity of pixels [38] can not yield accurate results.

Differing from the spatial sampling interval, the angular sampling interval of channel images is uniform. So this study proposes to directly use the spatial coordinates of pixels to obtain the corresponding pixels among neighboring channels so as to estimate the correct depth of object. The pixel is back-projected into space based on the geometry relationship of optical channels and the incidence angle of the pixel, thus, the spatial coordinates of pixels captured the same object will be nearly the same at the correct depth plane. The calculation of spatial coordinates is described in detail as follows.

In a camera system, suppose that $W$ , $\Delta x$ and $f$ are the width of an image sensor, the pixel pitch and the focal length, respectively, which are shown in Fig. 7. The angular resolution $\Delta\theta$ of this camera system is

$\Delta\theta=\arctan\frac{\Delta x}{f}$ (2)

and the FOV $\theta$ is

$\theta=\arctan\frac{W}{f}$ (3)

In our eCley system, the pixel pitch is 3.2 $\mu$ m; the channel image has a size of 55 $\times$ 55 pixels; and the focal length $f$ is 778 $\mu$ m. According to Eqs (2) and (3), the angular resolution $\Delta\theta=0.235^{\circ}$ and the FOV $\theta=6.345^{\circ}$ are obtained.

Figure 8.

Schematic graph for imaging process. Brown point: optical center; red point: image center; green point: object position; blue point: the projection position of the object on the image plane. The principal optic axis of the image plane and the image plane are not orthogonal and it has a horizontal offset angle $\Delta\varphi_{x}$ and a vertical offset angle $\Delta\varphi_{y}$ .

For a micro-lens channel shown in Fig. 8, the principal optic axis has a horizontal offset angle $\Delta\varphi_{x}$ and a vertical offset angle $\Delta\varphi_{y}$ . Let $Ic(Ic_{x},Ic_{y})$ be the center pixel of an image. An object point $T(D_{x},D_{y},Z)$ is projected into the image of this channel at $I(I_{x},I_{y})$ . Thus, the horizontal and vertical incident angles, $\Delta\psi_{h}$ and $\Delta\psi_{v}$ of object $T$ , are

$\displaystyle\Delta\psi_{h}=\Delta\varphi_{x}+(I_{x}-Ic_{x})*\Delta\theta$ $\displaystyle\Delta\psi_{v}=\Delta\varphi_{y}+(I_{y}-Ic_{y})*\Delta\theta$

The spatial coordinates of the object have the following relationship

$\displaystyle D_{x}=\tan(\Delta\psi_{h})*Z$ $\displaystyle D_{y}=\tan(\Delta\psi_{v})*Z$

Once the depth $Z$ is known, the spatial coordinate of an object can be obtained. But the spatial coordinate of each channel image is not based on a unified coordinate system. The origins of the coordinate systems of different channel images are distinct, so a coordinate transformation is required.

In a unified coordinate system, let the optical center of the central channel $c(c_{x},c_{y})$ be the original point $O(O_{x},O_{y})=(0,0)$ , and $B$ be the interval between adjacent channels. The unified coordinate $O_{i}(O_{ix},O_{iy})$ of the optical center in the $i^{\rm th}$ channel can be computed as

$\displaystyle O_{ix}=O_{x}+(i_{x}-c_{x})*B$ $\displaystyle O_{iy}=O_{y}+(i_{y}-c_{y})*B$

where $i_{x}\in[-N,N]$ and $i_{y}\in[-M,M]$ indicate the horizontal and vertical positions of channels. The total number of channels is $(2*N+1)\times(2*M+1)$ , and the central channel $c(c_{x},c_{y})$ is located at $(i_{x},i_{y})=(0,0)$ .

Figure 9.

Schematic graph showing the object captured by 4 adjacent channels.

Consequently, for a specified depth $Z$ , according to Eqs (5) and (6), the unified spatial coordinate $T^{\prime}(D^{\prime}_{x},D^{\prime}_{y},Z)$ of a pixel $I(I_{x},I_{y})$ in the $i^{\rm th}$ channel is calculated as

$\displaystyle D^{\prime}_{x}=D_{x}^{i}+O_{ix}$ $\displaystyle D^{\prime}_{y}=D_{y}^{i}+O_{iy}$

where the superscript $i$ is added to Eq. (7) to indicate the $i^{\rm th}$ channel.

3.3 Matching cost computation

Matching cost computation is to calculate the intensity difference between a pixel and its corresponding pixels at adjacent channels. Like stereo matching, the depth is estimated based on several depth plane candidates. With the spatial coordinates of each pixel in each depth plane obtained in Eq. (7), the matching cost of these pixels with similar spatial coordinates can be obtained. If the pixel is placed at the proper distance, the matching cost would be small due to the similar intensities of its corresponding pixels.

According to the eCley parameters, there is a slight offset between the principal optic axes of adjacent channel images, though an object in the FOV of eCley can still be captured by at least four channels (except for the object located at the marginal of the FOV), as show in Fig. 9. Due to the fact that the direction of FOV between adjacent channels has a shift and the width of each image is only $55$ pixels, which lead to an inherent parallax between two adjacent channel images and little usefully contextual information for depth matching, respectively, the traditional two-image stereo correspondence algorithms easily cause mismatching. Although the method proposed by Ziegler et al. [59] estimates the disparity of one channel image with four neighboring images, essentially, this method is based on two image stereo corresponding algorithm.

Figure 10.

An example for showing the pixels projected to the spatial plane, where the black points denote the projected point of pixels; blue and green boxes indicate the sampling regions of one pixel. (a) The spatial plane with several pixels projected to this plane. (b) A magnified figure of the green box in (a).

In order to gather as much information as possible for depth estimation, this study uses the pixels from different images located at similar spatial coordinate to compute the matching cost simultaneously. For a specified depth plane, all pixels are projected to this plane based on Eq. (7) in Section 3.2. As shown in Fig. 10, the black dots are the spatial points at this depth plane. The spatial coordinates of several pixels from different images may not be consistent and have a small deviation. According to the expected resolution of the reconstructed image, the depth plane is divided into multiple small and non-overlapped windows, each of which indicates a sampling area of a re-sampling point. Thus, the pixel points within the small window are considered as the supporting points of the re-sampling point. As shown in Fig. 10b, the orange dot is a re-sampling point; the pixel points located in the blue box are the supporting points of the re-sampling point; the intensity of the orange dot is the average intensity of the supporting points. Moreover, the sampling interval is related to the depth. The farther the depth plane is, the larger the sampling interval will be. For simplicity, those points in all candidate depth planes are re-projected to a fixed depth plane $Z_{f}$ , thus the sampling intervals at different distances are consistent. The coordinates of each point are obtained by using Eq. (8)

$\displaystyle X_{zf}=D^{\prime}_{x}*\frac{Z_{f}}{Z}$ $\displaystyle Y_{zf}=D^{\prime}_{y}*\frac{Z_{f}}{Z}$

where $Z$ is a known depth and $S(X_{zf},Y_{zf},d)$ is the spatial coordinate of a pixel, where we suppose all pixels are projected to the plane at depth $d$ . The corresponding intensity is a $3$ -dimensional vector $I(X_{zf},Y_{zf},d)$ of the RGB space. The depth plane is re-sampled with a resolution of $W*H$ . The re-sampling interval $R w$ is

$\displaystyle Rw=$ $\hskip 10.0pt\min\left\{\frac{\max(X_{zf})-\min(X_{zf})}{W},\right.$ $\hskip 10.0pt\left.\frac{\max(Y_{zf})-\min(Y_{zf})}{H}\right\}$

The sampling area of a re-sampling point $S^{\prime}(rx,ry,$ $d)$ is a small window centered at $S^{\prime}(rx,ry,d)$ with a size of $Rw\times Rw$ . Thus the intensity $I^{\prime}(rx,ry,d)$ of the re-sampling point and the matching cost $C(rx,ry,d)$ of this point is

$\displaystyle I^{\prime}(rx,ry,d)=\frac{1}{Ni}\sum_{k}I(k)$ $\displaystyle C(rx,ry,d)=$ $\displaystyle\sqrt{\frac{1}{Ni}\sum_{k}\sum_{m\in{r,g,b}}(I_{m}(k)-I^{\prime}_% {m}(rx,ry,d))^{2}}$

where $k$ is a 3-dimensional vector that indicates the supporting pixel points of $S^{\prime}(rx,ry,d)$ ; $I(\cdot)$ is a intensity function; $N i$ is the number of supporting points; and $m$ indicates the color channels.

3.4 Matching cost aggregation

Matching cost aggregation is widely used in local stereo matching methods to aggregate the local matching information and improve the matching results. In this study, due to the matching ambiguity, the depth estimation results based on the point-wise matching cost may be noisy and less accurate. So a local window-based edge preserved method based on the bilateral filter is used for cost aggregation [54, 12]. Assume that a window $P$ with a size of $[\textit{WL},\textit{WL}]$ is centered at point $S^{\prime}(rx,ry,d)$ . Thus its aggregated cost $\textit{Cs}(rx,ry,d)$ is computed as

$\displaystyle\textit{Cs}(rx,ry,d)=$ $\displaystyle\sum_{(qx,qy)\in P}\frac{\omega(qx,qy,d)\times C(qx,qy,d)}{\omega% (qx,qy,d)}$

where $(qx,qy)$ indicates the coordinate of the pixel points located in the window $P$ ; $\omega(qx,qy,d)$ is the weight of point $S^{\prime}(qx,qy,d)$ . The weight is based on an edge-preserved filter and can be written as

$\displaystyle\omega(qx,qy,d)=\exp$ $\displaystyle∼{}∼{}\left(-\left(\frac{\triangle cl(qx,qy,d)}{\gamma_{a}}+\frac% {\triangle gl(qx,qy,d)}{\gamma_{b}}\right)\right)$ $\displaystyle\triangle cl(qx,qy,d)=$ $\displaystyle\sqrt{\sum_{m\in{r,g,b}}(I^{\prime}_{m}(qx,qy,d)-I^{\prime}_{m}(% rx,ry,d))^{2}}$ $\displaystyle\triangle gl(qx,qy,d)=$ $\displaystyle\sqrt{\sum_{cg\in{x,y}}(S^{\prime}_{cg}(qx,qy,d)-S^{\prime}_{cg}(% rx,ry,d))^{2}}$

where $\triangle cl$ measures the color difference between point $S^{\prime}(qx,qy,d)$ and point $S^{\prime}(rx,ry,d)$ ; $\triangle gl$ measures the Euclidean distance between two points; $\gamma_{a}$ and $\gamma_{b}$ are weighting parameters and determined empirically; $c g$ indicates the horizontal and vertical coordinates.

3.5 Initial depth estimation and image reconstruction

In this section, the depth with the minimal matching cost is selected as the initial estimated depth. For different values of distance $d$ , the depth information $d^{\prime}(rx,ry)$ of a point can be selected by using Winner-Takes-All (WTA) method [38].

$\displaystyle d^{\prime}(rx,ry)=\arg\min_{d}\textit{Cs}(rx,ry,d)$

In addition, it is well known that the parallax of an object captured by neighboring images is related to the depth of the object. If an object is placed in an improper depth plane, the obtained image of the object will be blurred. So the image reconstruction is based on depth information. Concretely, the intensity $I^{\prime}(rx,ry,d)$ of a point $S^{\prime}(rx,ry,d)$ at specified distance $d$ is the average intensity of pixels in the sampling area. Multiple intensities of $S^{\prime}(rx,ry,d)$ are obtained based on multiple values of $d$ . Among these intensities, the intensity corresponding to the minimal matching cost is considered as the intensity of this point. If the difference of these intensities is small, the matching cost will be low, which indicates that these points corresponding to these intensities are more likely to be the same spatial point or they are different points with similar intensities. While the different points will not lead to the reduction of the image quality. So this method is simple and effective. The intensity of the reconstructed image is

$\displaystyle\textit{In}(rx,ry)=I^{\prime}(rx,ry,d^{\prime}(rx,ry))$

3.6 Depth and image refinement

However, the depth obtained may still suffer mismatching. Therefore, a reliable depth map is constructed by comparing the minimal cost $C(rx,ry,d^{\prime})$ with the second minimal cost $C(rx,ry,d^{\prime}_{1})$ . If the remainder that the difference between the two costs divided by $C(rx,ry,d^{\prime}_{1})$ is larger than a threshold $\lambda$ , the depth is marked as a reliable depth, otherwise, it is not a reliable depth. The $\lambda$ is determined empirically, the larger value indicates a more reliable depth. Thus, the reliable depth map $D(rx,ry)$ is computed by Eq. (3.6).

$\displaystyle D(rx,ry)=d^{\prime}(rx,ry),\ re(rx,ry)>\lambda$ $\displaystyle D(rx,ry)=0,\ \ \ re(rx,ry)\leqslant\lambda$ (15) $\displaystyle re(rx,ry)=\frac{C(rx,ry,d^{\prime})-C(rx,ry,d^{\prime}_{1})}{C(% rx,ry,d^{\prime}_{1})}$

Now a sparse reliable depth map and an initial reconstructed image $\textit{In}(px,py)$ are obtained. For simplicity, two dimensional $p=(px,py)$ is used. To fill the depth in those invalid pixels, based on the initial reconstructed image, an edge-aware interpolation method called matting Laplacian [22] is performed. The depth refinement problem of $D$ can be formulated as minimizing the follow function

$\displaystyle E(d_{f})\!=\!d_{fv}^{\top}Md_{fv}{+}\lambda_{d}(d_{fv}{-}D_{v})^% {\top}O(d_{fv}{-}D_{v})$

where $d_{fv}$ and $D_{v}$ are the vector forms of the final depth map $d_{f}$ and $D$ , respectively. $O$ is a diagonal matrix whose $O_{pp}$ is $1$ if pixel $p$ is a valid pixel, and $0$ otherwise. $\lambda_{d}$ balance between fidelity to the consistent depth map and smoothness. $M$ is the matting Laplacian matrix and the ( $p, q$ ) element of $M$ is

$\displaystyle\sum\limits_{g|s\in\vartheta_{g}}\left(\delta_{pq}-\frac{1}{|% \vartheta_{g}|}(1+(\textit{In}(p)-\mu_{g})\right.$ $\displaystyle\quad∼{}\left.\left(\Sigma_{g}+\frac{\epsilon_{g}}{|\vartheta_{g}% |}U_{3}\right)^{-1}(\textit{In}(q)-\mu_{g}))\right)$

where $U_{3}$ is a $3\times 3$ identity matrix, $\mu_{g}$ and $\Sigma_{g}$ are $3\times 1$ mean vector and $3\times 3$ covariance matrix of colors in a window $\vartheta_{g}$ , respectively. $|\vartheta_{g}|$ is the size of $\vartheta_{g}$ and $\epsilon_{g}$ is a regulation parameter. $\textit{In}(p)$ and $\textit{In}(p)$ are the colors of input image In at pixel $p$ and $q$ respectively.

Therefore, the final depth map $d_{fv}$ of image In can be obtained by solving the follow linear system,

$(M+\lambda_{d}O)d_{fv}=\lambda_{d}OD_{v}$ (18)

Then based on the final depth map $d_{f}$ , a more smooth image $I_{s}(p)$ is obtained,

$I_{s}(p)=I^{\prime}(p,d_{f}(p))$ (19)

4. Experimental results

In this section, the feasibility of the proposed method is tested and the comparisons with the two methods in [59, 46] are made in terms of the quality of reconstructed images and the depth map. Actually, the depth estimation method in [59] is replaced by the cost filter technique [15] in this paper. Both the proposed method and the other algorithms for comparison were implemented in MATLAB 2016b.

Table 2
Parameters used in this study

Reconstructed image resolution	640 $\times$ 640
Angular resolution $\Delta\theta$	0.235 ${}^{\circ}$
Pixel pitch	3.2 $\mu$ m
$\gamma_{a}$	0.1
$\gamma_{b}$	3
Window size of $P$	11 $\times$ 11
$\lambda$	0.05
$\lambda_{d}$	0.005

Figure 11.

Input image captured by the eCley camera system.(a) Original image before flat field correction. (b) Image after flat field correction.

Since there is no ground truth for our image sets, the quantitative comparisons are performed by comparing the reconstruction errors between the final image and the original channel images. Suppose a reconstructed image $I_{s}(p)$ with a resolution [ $W_{s},H_{s}$ ], pixel $p$ is obtained from a supporting pixel sets $I^{s}(p)$ of original channels images. Then the reconstructed mean squared error (MSE) is,

$\displaystyle\textit{mse}=\frac{1}{W_{s}{\times}H_{s}}\sum\limits_{p}\!\left(% \frac{1}{N_{s}}\sum\limits_{k\in N_{s}}\!(I_{s}(p){-}I^{s}(k)^{2})\right)$

where $N_{s}$ is the number of supporting pixels.

The captured image is extracted from all channels and shown in Fig. 11, where the left is the original image without flat field correction and the right is the image after flat field correction. In this paper, due to the large oblique distortion, only $13\times 13$ channels are considered. The most time-consuming step is the cost aggregation. The larger the window size of P is, the more time is needed for cost aggregation. Figure 12 shows the executing time at the cost aggregation step with different window sizes and the MSE of the final reconstructed image. As shown in the figure, the computing time increases as the window size goes larger. In the meantime, the MSE decreases. Therefore, to balance the time consumption and the result accuracy, the window size with $11\times 11$ is used in the study. The parameters used in this paper are listed in Table 2.

Table 3

MSE of the reconstructed image

Image	G3D	[59]	[46]
Image 1	0.048	0.0823	0.0843
Image 2	0.0536	0.0593	0.0654
Image 3	0.0535	0.0766	0.0829
Image 4	0.045	0.0464	0.0523
Image 5	0.0374	0.0465	0.0439

Figure 12.

The executing time at the cost aggregation step with different window sizes $P$ . The X-axis is the width of the windows.

Figure 13 shows five test images and the corresponding reconstructed images. The first, second and fifth images both have three fronto-parallel objects placed in front of the eCley system, the third image has two objects, different from the above two images, they are placed on a textured slanted plane, and the fourth image is a pure slanted plane, all of the objects shown in image are placed up to 200 mm away from the eCley system. As shown in Fig. 13, the proposed method has a better performance than the one gained by method in [59] and by the method in [46], especially for the third and the fourth image, where the words on the slanted plane have less artifact. In addition, the MSE listed in Table 3 also indicates that the proposed method has the minimal reconstructed error among the three methods.

Figure 13.

Test images and the reconstructed image. (a) Four test images. (b) Reconstructed image based on the depth of proposed method. (c) Reconstructed image based on [46]. (d) Reconstructed image based on [59].

Figure 14.

The final depth map. (a) The proposed method. (b) Depth in [46]. (c) Depth in [59].

Figure 15.

The 3D point cloud of image 1. (a) The proposed method. (b) Result in [46]. (c) Result in [59].

Figure 16.

The 3D point cloud of image 2. (a) The proposed method. (b) Result in [46]. (c) Result in [59].

Figure 17.

The 3D point cloud of image 3. (a) The proposed method. (b) Result in [46]. (c) Result in [59].

Figure 18.

The 3D point cloud of image 4. (a) The proposed method. (b) Result in [46]. (c) Result in [59].

Figure 19.

The 3D point cloud of image 5. (a) The proposed method. (b) Result in [46]. (c) Result in [59].

Figure 14 shows the final depth map obtained by the proposed method and the methods in [46, 59]. Since the principal optical axis among eCley apertures is not parallel, therefore, the estimated depth obtained by directly using the stereo method on the eCley images will suffer from the problem that the depth of an object will decrease as the object approaches to the image border. As shown in the right two columns in Fig. 14, especially for the images of the first two rows, the objects are fronto-parallel to the eCley system. While in the proposed method, this problem is largely avoided, and the depth of an object is more accurate.

Figures 15–19 show the 3D point clouds of the reconstructed images based on the depth maps. As shown in Figs 15, 16 and 19, the objects obtained by the proposed method are located at the same plane, and the objects are accurately recovered, while based on the methods in [46, 59], the objects near the border of the image are not located at the same plane. In addition, for the slanted plane shown in Figs 17 and 18, the proposed method can also accurately recover the patterns on the slanted plane, while in the methods [46, 59], the reconstructed images are much worse and the patterns on the slanted plane are hardly recognised. Therefore, based on the depth maps and the 3D point cloud maps, the reconstructed images obtained by the proposed method and the methods in [46, 59] are not much different, while the stereo images have much difference. The distance of an object in the proposed method has not dramatic change, but for the methods in [46, 59], the distance of the object located in the corner of the reconstructed image always has a dramatic change.

It is worth pointing out that there is not much difference in the aspect of computing time, several minutes on Matlab 2016b, between the proposed method and the one in [59], the computing times are listed in Table 4. The two methods use the same cost aggregation approach. As for the method in [46], since this is a global method, and a lot of iterations is needed to obtain a good result, the computing time is about two time longer than the proposed method to obtain a sufficient result. In the proposed method, all neighboring images of the reference image are used for cost matching, while the methods in [59, 46] use a two-view stereo matching method and fusion. In this study, no fusion process is required.

Table 4

Average computing time

Method	G3D	[59]	[46]
Time (sec)	169.53	143.34	302.78

5. Conclusions

In this work, a three-dimensional image processing method is presented by using the multiple apertures images of eCley to reconstruct a focused image and to estimate depth information. Based on the projection geometry, the proposed method directly obtains an initial depth map and a reconstructed image according to the spatial coordinate of pixels. The reconstructed image and the initial depth map are further used to obtain a dense and consistent depth map by using an edge-aware interpolation method. Subsequently, the reconstructed image is further refined based on the final dense and consistent depth map. Therefore, the three-dimensional structure of the scene can be recovered. Differing from the depth map obtained by stereo matching and fusion method, the proposed method back-projects the pixel into space and directly uses the images from all channels to estimate the depth map. The information from multiple channels is merged into a single matching step. The reconstructed image and the depth map are compared with the methods in literature. Experimental results show that the proposed method can correctly estimate the reconstructed image and the depth map, and achieve a better reconstructed image and depth results than two state-of-the-art stereo methods. Future work will focus on the auto calibration method to rectify the channel image and then further improve the depth result and the use of super-resolution image reconstruction because the results still contain some errors.

Footnotes

Acknowledgments

This work is supported by the National Natural Science Foundation of China (61373047, 61672437, 61702428 and 51641506), Sichuan Science and Technology Program (18ZDYF2877, 18ZDYF1985), the State Key Laboratory of Robotics Program (2014-O09), Scientific Research Foundation of CUIT (J201508), Scientific Research Fund of Sichuan Provincial Science & Technology Department (2015GZ0304); Fund of Robot Technology Used for Special Environment Key Laboratory of Sichuan Province (14zxtk04) and China Scholarship Council.

References

Adams

. Hand held dyson relay lens for anti-counterfeiting. In Proc IEEE International Conference on Image Systems and Technologies 2010; 273-278.

Aydin

Demirtas

Ates

Tunga

. Emotion recognition with eigen features of frequency band activities embedded in induced brain oscillations mediated by affective pictures. Int J Neural Syst 2016; 26(3).

Baumgartner

Flesia

Gimenez

Pucheta

. A new image segmentation framework based on two-dimensional hidden markov models. Integrated Computer-Aided Engineering 2016; 23(1): 1-13.

Ben-Ari

. A unified approach for registration and depth in depth from defocus. IEEE Transactions on Pattern Analysis and Machine Intelligence 2014; 36(6): 1041-1055.

Brückner

Duparré

Leitel

Dannberg

Bräuer

Tünnermann

. Thin wafer-level camera lenses inspired by insect compound eyes. Optics Express 2010; 18(24): 24379-24394.

Curkovic

Vucina

Curkovic

. Enhanced 3d parameterization for integrated shape synthesis by fitting parameter values to point sets. Integrated Computer-Aided Engineering 2017; 24(3): 241-260.

Delibasis

Georgakopoulos

Kottari

Plagianakos

Maglogiannis

. Geodesically-corrected zernike descriptors for pose recognition in omni-directional images. Integrated Computer-Aided Engineering 2016; 23(2): 185-199.

Duparré

Wippermann

Dannberg

Reimann

. Chirped arrays of refractive ellipsoidal microlenses for aberration correction under oblique incidence. Optics Express 2005; 13(26): 10539-10551.

El-Sallam

Boussaid

. Spectral-based blind image restoration method for thin tombo imagers. Sensors 2008; 8(9): 6108-6124.

10.

EI-Sallam

Boussaid

. A high resolution color image restoration algorithm for thin tombo imaging systems. Sensors 2009; 9(6): 4649-4668.

11.

Furukawa

Hernández

. Multi-view stereo: A tutorial. Found Trends Comput Graph Vis Jun 2015; 9(1-2): 1-148.

12.

Sun

Tang

. Guided image filtering. IEEE Transactions on Pattern Analysis and Machine Intelligence 2013; 35(6): 1397-1409.

13.

Hirschauer

Adeli

Buford

. Computer-aided diagnosis of parkinson’s disease using enhanced probabilistic neural network. J Medical Systems 2015; 39(11): 179.

14.

Horisaki

Irie

Nakao

Ogura

Toyoda

Masaki

Tanida

. 3d information acquisition using a compound imaging system. In Proc SPIE 2007; 66950F.

15.

Hosni

Rhemann

Bleyer

Rother

Gelautz

. Fast cost-volume filtering for visual correspondence and beyond. IEEE Transactions on Pattern Analysis and Machine Intelligence 2013; 35(2): 504-511.

16.

Hsu

. Automatic atrium contour tracking in ultrasound imaging. Integrated Computer-Aided Engineering 2016; 23(4): 401-411. Online Available: https://doi.org/10.3233/ICA-160517

17.

Hsu

. A hybrid approach for brain image registration with local constraints. Integrated Computer-Aided Engineering 2017; 24(1): 73-85.

18.

Iacca

Caraffini

Neri

. Memory-saving memetic computing for path-following mobile robots. Appl Soft Comput 2013; 13(4): 2003-2016.

19.

Iliya

Neri

. Towards artificial speech therapy: A neural system for impaired speech segmentation. Int J Neural Syst 2016; 26(6): 1-16.

20.

Jiang

Zhu

Kuhnert

K-D

Kuhnert

. Distance measuring using calibrating subpixel distances of stereo pixel pairs in artificial compound eye. In Proc International Conference on Informative and Cybernetics for Computational Social Systems 2014; 118-122.

21.

Kitamura

Shogenji

Yamada

Miyatake

Miyamoto

Morimoto

Masaki

Kondou

Miyazaki

Tanida

, et al. Reconstruction of a high-resolution image on a compound-eye image-capturing system. Applied Optics 2004; 43(8): 1719-1727.

22.

Levin

Lischinski

Weiss

. A closed-form solution to natural image matting. IEEE Transactions on Pattern Analysis and Machine Intelligence 2008; 30(2): 228-242.

23.

Zhou

. Increasing N200 potentials via visual stimulus depicting humanoid robot behavior. Int J Neural Syst 2016; 26(1).

24.

Chen

Zhang

Tong

Zhou

. Real-time control of an exoskeleton hand robot with myoelectric pattern recognition. Int J Neural Syst 2017; 27(5): 1-11.

25.

Martínez-Cañada

Morillas

Pino

Ros

Pelayo

. A computational framework for realistic retina modeling. Int J Neural Syst 2016; 26(7): 1-18.

26.

Mesquita

Mello

CAB

. Object recognition using saliency guided searching. Integrated Computer-Aided Engineering 2016; 23(4): 385-400.

27.

Nasrollahi

Moeslund

. Super-resolution: A comprehensive survey. Machine Vision and Applications 2014; 25(6): 1423-1468.

28.

Neri

Mininno

. Memetic compact differential evolution for cartesian robot control. IEEE Com Int Mag 2010; 5(2): 54-65.

29.

Nitta

Shogenji

Miyatake

Tanida

. Image reconstruction for thin observation module by bound optics by using the iterative backprojection method. Applied Optics 2006; 45(13): 2893-2900.

30.

Oberdörster

Brückner

Wippermann

Bräuer

. Correcting distortion and braiding of micro-images from multi-aperture imaging systems. In Proc SPIE 2011; 78750B.

31.

Oberdörster

Brückner

Wippermann

Bräuer

Lensch

. Digital focusing and refocusing with thin multi-aperture cameras. In Proc SPIE 2012; 829907.

32.

Pan

Paun

Zhang

Neri

. Spiking neural p systems with communication on request. Int J Neural Syst, vol 27, 2017, to appear.

33.

Park

Kang

. Super-resolution image reconstruction: A technical overview. IEEE Signal Processing Magazine 2003; 20(3): 21-36.

34.

Poikolainen

Neri

Caraffini

. Cluster-based population initialization for differential evolution frameworks. Inf Sci 2015; 297: 216-235.

35.

Rigos

Tsekouras

Vousdoukas

Chatzipavlis

Velegrakis

. A chebyshev polynomial radial basis function neural network for automated shoreline extraction from coastal imagery. Integrated Computer-Aided Engineering 2016; 23(2): 141-160.

36.

Rostami

Neri

Epitropakis

. Progressive preference articulation for decision making in multi-objective optimisation problems. Integrated Computer-Aided Engineering, vol 24, 2017, to appear.

37.

Rostami

Neri

. Covariance matrix adaptation pareto archived evolution strategy with hypervolume-sorted adaptive grid algorithm. Integrated Computer-Aided Engineering 2016; 23(4): 313-329.

38.

Scharstein

Szeliski

. A taxonomy and evaluation of dense two-frame stereo correspondence algorithms. International Journal of Computer Vision 2002; 47(1-3): 7-42.

39.

Schechner

Kiryati

. Depth from defocus vs. stereo: How different really are they? International Journal of Computer Vision 2000; 39(2): 141-162.

40.

Siddique

Adeli

. Harmony search algorithm and its variants. IJPRAI 2015; 29(8).

41.

Siddique

Adeli

. Hybrid harmony search algorithms. International Journal on Artificial Intelligence Tools 2015; 24(6).

42.

Siddique

Adeli

. Gravitational search algorithm and itsvariants. IJPRAI 2016; 30(8): 1-22. Online Available: https://doi.org/10.1142/S0218001416390018

43.

Siddique

Adeli

. Physics-based search and optimization: Inspirations from nature. Expert Systems 2016; 33(6): 607-623.

44.

Siddique

Adeli

. Simulated annealing, its variants and engineering applications. International Journal on Artificial Intelligence Tools 2016; 25(6): 1-24.

45.

Soto

Adeli

. Recent advances in control algorithms for smart structures and machines. Expert Systems 2017; 34(2).

46.

Sun

Zheng

N-N

Shum

H-Y

. Stereo matching using belief propagation. IEEE Transactions on Pattern Analysis and Machine Intelligence 2003; 25(7): 787-800.

47.

Tirronen

Neri

Kärkkäinen

Majava

Rossi

. An enhanced memetic differential evolution in filter design for defect detection in paper production. Evolutionary Computation 2008; 16(4): 529-555.

48.

Volkel

Herzig

Nussbaum

Dandliker

Hugle

. Microlens array imaging system for photolithography. Optical Engineering 1996; 35(11): 3323-3330.

49.

Wang

Zhang

Neri

Jiang

Zhao

Gheorghe

Ipate

Lefticaru

. Design and implementation of membrane controllers for trajectory tracking of nonholonomic wheeled mobile robots. Integrated Computer-Aided Engineering 2016; 23(1): 15-30.

50.

Williams

Shaddix

. Simultaneous correction of flat field and nonlinearity response of intensified charge-coupled devices. Review of Scientific Instruments 2007; 78(12): 123702.

51.

Jiang

Zhang

Schoenemann

Neri

Zhu

Han

Kuhnert

K-D

. Artificial compound eye: A survey of the state of the art. Artificial Intelligence Review 2017; 48(4): 573-603.

52.

Xian

Subbarao

. Performance evaluation of different depth from defocus (dfd) techniques. In Proc SPIE 2005; 600009.

53.

Yamada

Mitsui

Asano

Tanida

Takahashi

. Development of ultrathin three-dimensional image capturing system. In Proc SPIE 2006; 60560V.

54.

Yoon

K-J

Kweon

I-S

. Locally adaptive support-weight approach for visual correspondence search. In Proc CVPR IEEE, June 2005; 924-931.

55.

Zeinali

Story

. Competitive probabilistic neural network. Integrated Computer-Aided Engineering 2017; 24(2): 105-118.

56.

Zhang

Han

Zou

Chen

. An efficient approach to directly compute the exact hausdorff distance for 3d point sets. Integrated Computer-Aided Engineering 2017; 24(3): 261-277.

57.

Zhang

Rong

Neri

Pérez-Jiménez

. An optimization spiking neural P system for approximately solving combinatorial optimization problems. Int J Neural Syst 2014; 24(5).

58.

Zhu

Xiao

Wang

Zhang

. Global path planning of wheeled robots using multi-objective memetic algorithms. Integrated Computer-Aided Engineering 2015; 22(4): 387-404.

59.

Ziegler

Zilly

Schaefer

Keinert

Schöberl

Foessel

. Dense lightfield reconstruction from multi aperture cameras. In Proc International Conference on Image Processing, Oct 2014; 1937-1941.

Geometry based three-dimensional image processing method for electronic cluster eye

Abstract

Keywords

1. Introduction

Table 1 Parameters of the eCley system in [5]

3.5 Initial depth estimation and image reconstruction

3.6 Depth and image refinement

Table 2 Parameters used in this study

Footnotes

Acknowledgments

References

Table 1
Parameters of the eCley system in [5]

Table 2
Parameters used in this study