3D Scanning System for Automatic High-Resolution Plant Phenotyping

TL;DR

This paper introduces a fast, low-cost 3D scanning system for plants that achieves high accuracy in reconstructing complex plant structures, facilitating automated phenotyping.

Contribution

The authors developed a novel 3D scanning platform with improved calibration and background removal, enabling high-resolution plant imaging with minimal cost and time.

Findings

Reconstructed plant models show 3-4% measurement accuracy compared to manual methods.

System effectively captures thin leaves and fine stems in 3D models.

Scan times range from 3 to 30 minutes depending on image set size.

Abstract

Thin leaves, fine stems, self-occlusion, non-rigid and slowly changing structures make plants difficult for three-dimensional (3D) scanning and reconstruction -- two critical steps in automated visual phenotyping. Many current solutions such as laser scanning, structured light, and multiview stereo can struggle to acquire usable 3D models because of limitations in scanning resolution and calibration accuracy. In response, we have developed a fast, low-cost, 3D scanning platform to image plants on a rotating stage with two tilting DSLR cameras centred on the plant. This uses new methods of camera calibration and background removal to achieve high-accuracy 3D reconstruction. We assessed the system's accuracy using a 3D visual hull reconstruction algorithm applied on 2 plastic models of dicotyledonous plants, 2 sorghum plants and 2 wheat plants across different sets of tilt angles. Scan times ranged from 3 minutes (to capture 72 images using 2 tilt angles), to 30 minutes (to capture 360 images using 10 tilt angles). The leaf lengths, widths, areas and perimeters of the plastic models were measured manually and compared to measurements from the scanning system: results were within 3-4% of each other. The 3D reconstructions obtained with the scanning system show excellent geometric agreement with all six plant specimens, even plants with thin leaves and fine stems.

Tables

Table 1. TABLE I: Manual leaf measurements of large plastic plant.

Leaf	Length (mm)	Width (mm)	Perim. (mm)	Area (mm2)
B	152.6	74.5	350.89	8535.6
C	186.9	96.6	436.16	13303.4
D	221.0	115.6	519.1	18868.0

Equations

$${}^{K}p={}^K{\xi }_T\cdot {}^{T}p=[R|t]\cdot [X,Y,Z,1{]}^T$$

$$r=\sqrt{{\left(\frac{x}{z}\right)}^2+{\left(\frac{y}{z}\right)}^2}$$

$$u=\frac{fx}{z}(1+d_1{r}^2+d_2{r}^4)+c_u$$

$$v=\frac{fy}{z}(1+d_1{r}^2+d_2{r}^4)+c_v$$

$${}^{K^{\prime }}p={}^{K^{\prime }}{\xi }_K\cdot {}^{K}p={}^{K^{\prime }}{\xi }_K\otimes {}^K{\xi }_T\cdot {}^{T}p$$

$$a(t_x-x_m)+b(t_y-y_m)+c(t_z-z_m)=0$$

$${}^L{\xi }_W={}^L{\xi }_K\otimes {}^K{\xi }_W$$

$$\omega =\arctan \left(\frac{|\mathbf{n}\times \mathbf{y}|}{\mathbf{n}\cdot \mathbf{y}}\right)$$

$$\mathbf{w}=\frac{\mathbf{n}\times \mathbf{y}}{|\mathbf{n}\times \mathbf{y}|}$$

$$R_{\omega }=\cos \omega I+(1-\cos \omega )\mathbf{w}{\mathbf{w}}^T+\sin \omega [\mathbf{w}{]}_{\times }$$

$${}^L{\xi }_W=[I|t_0]$$

$${}^K{\xi }_W$$

$$=$$

$${}^{K^{\prime }}{\xi }_T={}^{K^{\prime }}{\xi }_K\otimes {}^K{\xi }_W$$

$${}^T{\xi }_W={}^T{\xi }_K\otimes {}^K{\xi }_W={\left({}^K{\xi }_T\right)}^{-1}\otimes {}^K{\xi }_W$$

$${}^W{\xi }_T$$

$$=$$

$${}^W{\xi }_{T0}$$

$$=$$

$${}^{K{\theta }_i}{\xi }_W$$

$${}^{K^{\prime }{\theta }_i}{\xi }_W$$

$$=$$

$${}^K{p}_{estimated}={}^{K{\theta }_i}{\xi }_W\otimes {}^W{\xi }_{T0}\cdot {}^{T}p$$

$$\Delta (u,v)=\left|L(u,v)-L^{\prime }(u,v)\right|$$

$$\begin{gathered} \Theta (u,v)=\left|A(u,v)-A^{\prime }(u,v)\right|+ \\ \left|B(u,v)-B^{\prime }(u,v)\right| \end{gathered}$$

$$\begin{gathered} \Psi (u,v)=\left|\frac{L(u,v)}{L(u+1,v)}-\frac{L^{\prime }(u,v)}{L^{\prime }(u+1,v)}\right|+ \\ \left|\frac{L(u,v)}{L(u,v+1)}-\frac{L^{\prime }(u,v)}{L^{\prime }(u,v+1)}\right| \end{gathered}$$

$$\Omega (u,v)=\frac{\alpha \Delta (u,v)+\beta \Theta (u,v)+\gamma \Psi (u,v)}{\alpha +\beta +\gamma }$$

Full text

1 Introduction

The structures of different plant species pose a range of challenges for 3D scanning and reconstruction. Various solutions to the issue of digitally imaging plants have been reported. For plants with larger leaves and simple structures (such as maize, sorghum, cereal seedlings) it is possible to capture a small number of digital images from various viewing angles (typically 3), analyse these in 2D, then develop a relationship between these 2D poses and the (destructively measured) leaf area and biomass of the species [1]. However, commercial systems using this approach are relatively expensive, with closed and proprietary analysis software; these have generally been deployed in large phenomics centres with conveyor based systems [2]. The 2D approach has difficulty in resolving concavities, leaf overlap and other occlusions; many of the powerful image analysis tools which can be applied to 3D meshes (e.g., volumetric based shape recognition and organ tracking) are more difficult in 2D [3].

Laser scanning, e.g., light detection and ranging (LIDAR), has been applied to plant digitisation but reconstructing a mesh from a pointcloud sufficiently dense to capture thin narrow leaves is computationally intensive. While this approach has been applied successfully to forestry [4, 5] and to statistical analysis of canopies, it is not well suited to extracting single plant attributes [4, 6]. Full waveform LIDAR is extremely expensive; simpler machine vision LIDAR systems of sufficient resolution can cost tens of thousands of dollars. Structured light approaches using affordable sensors such as the Kinect gaming sensor or camera-video projector setups do not offer the resolution or spatial repeatability to cope with complex plant structures [7, 8].

Recently, approaches using multiple images from a larger number of viewing angles have yielded promising results [3, 9, 10, 11, 12, 13], either by using a silhouette-based 3D mesh reconstruction method or patch based stereo reconstruction transforming the images to point clouds. Silhouette-based reconstruction is prone to errors in camera calibration and silhouette extraction. Image acquisition for 3D reconstruction remains largely manual [14], limiting the speed and accuracy of 3D reconstruction.

Existing background removal techniques for silhouette extraction [15, 16] are not reliable for complex plant structures. For current patch based methods, reconstruction quality is usually poor due to both weak-textures (patterns on leaves) and thin structures. In both cases, high accuracy 3D reconstruction requires a very rigid imaging station and the engineering required for sensor integration is costly. Recent work has included semi-automated approaches to modelling the plants [17, 18] increasing completeness and visual quality at the expense of throughput.

To address the shortcomings of existing scanning systems, we describe “PlantScan Lite”, an affordable automated 3D imaging platform system to accurately digitise complex plant structures. This can be readily built for less than AU$8000 of components (including DSLR cameras). The system has a number of novel elements including the use of high-resolution digital single-lens reflex (DSLR) cameras synchronised with a turntable; a high-accuracy, easy-to-setup camera calibration; and an accurate background removal method for 3D reconstruction.

2 Methods

PlantScan Lite uses a four step process of 3D scanning and reconstruction (Fig. 1):

Image acquisition. Multiple-view images of a calibration target, a plant of interest and background are captured using a turntable and tilting cameras. 2. 2.

Camera calibration estimates camera parameters for all captured images. 3. 3.

Image processing corrects for optical distortion and extract plant silhouettes (background removal). 4. 4.

3D reconstruction. This paper focuses on the visual hull reconstruction method which uses silhouette images and corresponding camera parameters and poses to create a 3D mesh model. The model is then processed to extract plant geometry information.

2.1 Image acquisition system

PlantScan Lite’s image acquisition system consists of (Fig. 2):

• •

Two Canon DSLR cameras (1, 2) at an angle of approximately 40 degrees. The cameras are powered by AC adapter kits and connect to a computer through USB cables for tethering capture.

• •

Two aluminum frames, one tilted (3) and one fixed to the ground (4). The tilted frame is to mount the cameras and move them up/down. The frames join together with two hinges where the axis of rotation crosses that of the turntable near the middle of a plant to be scanned. Two guiding bars (5) are attached to the lower frame (4) both to keep the tilted frame moving on a vertical plane.

• •

Two Phidget bipolar stepper motors (6, 7); one drives a turntable, the other moves the upper frame and the two cameras via a threaded rod (8). Both stepper motors have a gearbox to increase the rotation resolution (0.018 degree/step for the turntable) and torque (18 kg$\cdot$cm for the threaded rod). The turntable controller is synchronised with the cameras via software based on Phidget Python library [19] and Piggyphoto[20].

2.2 Camera calibration for turntable

Fig. 3 shows the schematic description of the system with a calibration target and a plant (notation follows [21]). A chessboard target has a local coordinate system $\{T\}$ with relative pose ${}^{T}\xi_{W}$ to the world $\{W\}$ which is fixed to the turntable. A camera at tilt position j has a local coordinate system $\{K\}$ and relative pose ${}^{K}\xi_{W}$ to the world coordinates $\{W\}$. Standard mono camera calibration only provides camera pose ${}^{K}\xi_{T}$ relative to target at angles where the chessboard pattern is visible to the camera. The aim of camera calibration for the turntable is to find intrinsic parameters and poses ${}^{K}\xi_{W}$ for all cameras at all positions. Three constraints are considered: a) turntable rotation angle is known, b) images captured by the same camera have the same intrinsic parameters, and c) different extrinsic parameters are assigned to different tilt positions as if these belong to different cameras that form a vertical arc. The calibration process consists of (1) mono and stereo camera calibrations (2) plane and circle fitting (3) camera to world pose derivation, and (4) optimisation of all camera parameters. Theses step are detailed below.

2.2.1 Mono and stereo camera calibration

The OpenCV library provides mono and stereo camera calibration routines [22, 23]. The calibration procedure starts by taking images of a chessboard target at different angles and distances. Corners of the squares on chessboard target are detected. To make this detection efficient, images of the chessboard calibration target are repeatedly scaled down 50% in a pyramid fashion to approximately 1K $\times$ 1K resolution at the top level. The detected corner positions are obtained at the top level (lowest resolution) using OpenCV’s function findChessboardCorners and tracked with subpixel accuracy using function cornerSubPix on images at lower levels of the pyramid.

Fig. 4 describes transformations between different coordinate systems. Thin arrows represent coordinate vectors. A pair of orthogonal arrows denotes a coordinate system. Thick arrows show relative pose between two coordinate systems. The coordinates (in mm) of a corner p with respect to the chessboard target coordinate system $\{T\}$ are represented as ${}^{T}p=[X,Y,Z=0,1]^{T}$. The position of p relative to a camera coordinate system $\{K\}$ is ${}^{K}p=[x,y,z]^{T}$. The relationship between ${}^{T}p$ and ${}^{K}p$ is expressed as:

[TABLE]

where $R$ is a $3\times 3$ rotation matrix, $t$ a $3\times 1$ translation vector. This rotation translation matrix $[R|t]$ represents the extrinsic parameters of camera K relative to target T. The rotation matrix can be represented as a $3\times 1$ angle-axis rotation vector, so extrinsic parameters have only 6 independent components $[r_{0},r_{1},r_{2},t_{0},t_{1},t_{2}]$. The “$\cdot$” operator is a matrix-to-vector multiplication.

Point p forms an image on camera sensor at coordinates ${}^{I}p=[u,v]^{T}$. An extended pinhole camera model is used to describe this relationship between ${}^{K}p$ and ${}^{I}p$ with radial optical distortion:

[TABLE]

[TABLE]

[TABLE]

where $f$ is focal length, $[c_{u},c_{v}]$ optical center on image, and $[d_{1},d_{2}]$ radial distortion coefficients. A vector $[f,c_{y},c_{v},d_{1},d_{2}]$ represents intrinsic camera parameters.

To calibrate the camera, multiple images of the same target are captured at different pan angles $\theta_{i}$ and tilt angles $\phi_{j}$, where i = 0 to $i_{max}$ and j = 0 to $j_{max}$. Mono camera calibration, using OpenCV’s calibrateCamera (based on [24]), takes lists of ${}^{T}p$ and ${}^{I}p$ and computes intrinsic parameters $[f,c_{y},c_{v},d_{1},d_{2}]_{k}$ for the camera, and extrinsic parameters $[R_{ijk}|t_{ijk}]$ for each of the images.

Fig. 4B shows the multiple camera setup with an additional camera $K^{\prime}$ at a different tilt position. The same point p is seen by the camera $K^{\prime}$ at coordinates ${}^{K^{\prime}}p$. The transformation between the two camera coordinate systems gives:

[TABLE]

The transformation ${}^{K^{\prime}}\xi_{K}$ between cameras K and $K^{\prime}$ is equal to stereo transformation $[R_{k,k^{\prime}}|t_{k,k^{\prime}}]$ (using OpenCV’s stereoCalibCameras). If K and $K^{\prime}$ are of the same camera but at different tilting angles $\phi_{j}$ and $\phi_{j+1}$, the transformation becomes $[R_{j,j+1}|t_{j,j+1}]$. The “$\otimes$” operator denotes matrix-to-matrix multiplication. The transformation can be applied repeatedly between successive camera pairs at different tilt positions. Given stereo transformations between successive cameras and the pose of the first camera, the poses of other cameras are also obtained.

2.2.2 Estimation of axis and centre of rotation

Location of the world coordinate system fixed to the rotation axis is found in two steps (Fig. 5): (1) compute the normal vector n of the plane containing the rotation orbit of the target’s corners, and (2) fit a circle to find the centre of rotation o. The normal vector and centre of rotation define the world coordinates.

First, rotation axis is estimated from the rotational motion of the chessboard target. From extrinsic parameters $[R_{ijk}|t_{ijk}]$, the same chessboard corner ${}^{T}p$ is seen as ${}^{K}p=[R_{ijk}|t_{ijk}]{}^{T}p$ moving on a circular orbit as shown in Fig. 5A. Note that the chessboard pattern and its positions can only be detected in 1/4 to 1/3 of the number of target images around the circular orbit. The equation of the orbit plane is expressed as:

[TABLE]

where $[x_{m},y_{m},x_{m}]$ is the centroid (or the mean) of all positions ${}^{K}p$ on the same orbit, and $[a,b,c]$ is the normal vector of the plane.

Equation (6) can be rewritten in matrix form $\mathbf{Bn}=0$ where $\mathbf{n}=[a,b,c]$ and $\mathbf{B}$ is a matrix containing chessboard positions relative to the centroid of all the positions. Vector $\mathbf{n}$ is an eigenvector corresponding to the smallest eigenvalue obtained from Singular Value Decomposition of matrix $\mathbf{B}$.

Second, to find the centre of rotation of the calibration target, a different coordinate system $\{L\}$ is used (Fig. 5B). $\{L\}$ is in fact equivalent to $\{K\}$ with a rotation transformation ${}^{L}\xi_{K}$ such that the y-axis is parallel to $\mathbf{n}$. In $\{L\}$, a 2D circle can be fitted onto the target point orbit and the centre of rotation can be obtained. The relationship between $\{L\}$ and $\{K\}$ with respect to $\{W\}$ is:

[TABLE]

The transformation ${}^{L}\xi_{W}$ has the form $[R_{\omega}|0]$, where $R_{\omega}$ is a rotation matrix whose angle-axis rotation vector $\omega\mathbf{w}=[\omega w_{x},\omega w_{y},\omega w_{z}]^{T}$ can be obtained as:

[TABLE]

[TABLE]

where $\omega$ is the rotation angle and $\mathbf{w}$ is a vector around which the rotation is applied to turn $\mathbf{n}$ to y-axis $\mathbf{y}$ of $\{W\}$. The bar denotes vector normalization.

$R_{\omega}$ is computed from $\omega\mathbf{w}$ by Rodrigues’ formula:

[TABLE]

After applying the rotation transformation $R_{\omega}$ to target positions, a circle can be fitted to $[z,x]^{T}$ coordinates by a Linear Least-Squares algorithm [25]. This fitting gives the centre of the orbit $[x_{0},y_{0},z_{0}]^{T}$, with $y_{0}$ to be the averaged y-component of the target point positions in $\{L\}$ coordinate system.

Now the world coordinate system is set at centre of rotation $\mathbf{o}$ and with its axes parallel to those of $\{L\}$, the transformation from $\{L\}$ to $\{W\}$ is:

[TABLE]

where $t_{0}=[x_{0},y_{0},z_{0}]^{T}$. As a result, pose of camera K relative to $\{W\}$ can be expressed as:

[TABLE]

2.2.3 Estimation of camera poses relative to the world coordinate system fixed to the turntable axis

Since the pose ${}^{K}\xi_{W}$ of the camera $K$ at tilt position $\phi_{j}$ is obtained, the pose of any additional camera $K^{\prime}$ relative to the world coordinate system can be obtained from a given stereo transformation ${}^{K}\xi_{K}$ as shown in Fig. 6A:

[TABLE]

Similarly, the pose of the target relative to the world coordinate system is expressed as:

[TABLE]

We are interested in the reverse transformation to obtain target point position in the world coordinate system:

[TABLE]

Since there are multiple estimates of ${}^{W}\xi_{T}$ for different rotation angle $\theta_{i}$, ${}^{W}\xi_{T0}$ for zero rotation angle is obtained by applying an inverse of the rotation to the corresponding pose:

[TABLE]

where $R_{\omega}$ is the matrix obtained from angle-axis rotation vector (equation (10)), and $R^{-1}_{\theta_{i}}=R^{T}_{\theta_{i}}$.

Since the world coordinate systems and the target are fixed, the camera coordinate system needs to move in a circle around y axis of $\{W\}$ to represent the correct relative motion seen by the camera, as shown in Fig. 6B. The pose of camera K for rotation angle $\theta_{i}$ is:

[TABLE]

The pose of camera $K^{\prime}$ for rotation angle is :

[TABLE]

2.2.4 Optimisation to refine camera parameters

A nonlinear least-square optimisation is applied to refine estimates of camera intrinsic parameters $[f,c_{u},c_{v},d_{1},d_{2}]$ and angle-rotation vector and translation vector of camera pose ${}^{K}\xi_{W}$ and ${}^{K^{\prime}}\xi_{W}$ at different tilt angles, and target inverse pose ${}^{W}\xi_{T0}$. Optimisation seeks to minimize is the pixel distance between projected target corners to camera and the corners on the actual images.

Estimated position of chessboard corner relative to camera K at rotation $\theta_{i}$:

[TABLE]

A pinhole camera projection with radial distortion is applied with equations (2), (3) and (4) to obtain corresponding image point ${}^{I\theta_{i}}p_{estimated}$. This is applied to other tilt positions and the second camera. The squared distance between image projection ${}^{I\theta_{i}}p_{estimated}$ of estimated corners and their detected image positions is minimised.

2.3 Background subtraction

Plants with large leaves can cast strong shadows, so a simple image threshold will not completely remove the background. We found that a shadow removal algorithm based on static background proposed in [15] performs background removal for thin leaves more accurately and with less computation than other techniques [26, 27, 28]. Here, we extend the technique of [15] for LAB color space [29], further improving background removal accuracy.

Suppose $L(u,v)$, $A(u,v)$ and $B(u,v)$ are the luminance and two color channels of the current image and $L^{\prime}(u,v)$, $A^{\prime}(u,v)$ and $B^{\prime}(u,v)$ the corresponding image channels of the background image. Three error functions applied to each pixel position $[u,v]$ are defined as follows:

[TABLE]

[TABLE]

[TABLE]

These error functions $\Delta(u,v)$, $\Theta(u,v)$ and $\Psi(u,v)$ represent the differences in luminance, color and texture respectively. An overall score is computed to determine a pixel as foreground or background:

[TABLE]

where the values of $\alpha$, $\beta$, $\gamma$ are found empirically.

A threshold $t$ is applied to $\Omega(u,v)$ to separate background and foreground. For our images, $\alpha=0.1$, $\beta=0.5$, $\gamma=0.4$ and $t=5$ to 10 were found to work well. Unlike [15], our proposed technique allows for segmentation of dark objects (such as the plant pots) and this is controlled via coefficient $\alpha$. Fig. 7 shows an example of background removal using our proposed algorithm.

2.4 3D reconstruction

2.4.1 Bounding box estimation

The bounding box of the subject is obtained in two steps:

An initial 3D bounding box is estimated based on silhouettes from the most horizontal camera view. These silhouette images are overlapped/combined into a single image and 2D bounding box is computed (Fig. 8). The Y axis and the origin are projected onto this overlapped image. The crossing points of the projected axis with the 2D bounding box are mapped back to the turntable axis in 3D space to obtain $y_{min}$ and $y_{max}$. The back projection of the rectangle width to the world origin gives a single magnitude for $x_{min}$, $x_{max}$, $z_{min}$ and $z_{max}$. 2. 2.

A refined bounding box is calculated from a 3D reconstruction at a low resolution ($128^{3}$ voxels) using the initial bounding box. This takes only a few seconds to compute. Particularly, we found that no thin parts of the plant are missing when reconstructed at low resolution. As a result, the refined bounding box tightly contains the 3D space of the plant.

2.4.2 Volume reconstruction

In this work, the 3D plant was reconstructed using a visual-hull volume carving approach. This method recovers sharp and thin structures common to plants (although a major drawback is that it cannot correctly recover concave surfaces, making reconstructions of curved surfaces such as leaves thicker than they should be). There may be plant movements induced by air circulation or mechanical vibration which can be accounted for by some tuning during reconstruction. The reconstruction method consisted of 3 steps:

A 3D volume equal to the bounding box is generated and split into voxels. Each voxel is repeatedly projected into the silhouettes and its 2D signed distance to the nearest boundary of each silhouette is calculated. If the distance is negative (outside) in any of the silhouettes, the voxel is flagged as empty. To accommodate some uncertainty in the silhouettes and plant movements, a voxel is set to be removed if it is outside more than a fixed number of silhouettes (3 is chosen in this paper). The process repeats until the end where remaining voxels form a 3D hull model of plant. An octree structure is used for voxel removal from lowest resolution to the highest resolution [30], giving a 3X speedup as compared to processing all voxels of a full resolution 3D volume. 2. 2.

Removal of pot and pot carrier. Since we are only interested in the plant, the pot and pot carrier need to be removed to simplify mesh analysis as well as reduce the mesh size. One method is to subtract voxels inside a given bounding tapered cylinder of the pot and pot carrier. Fig. 9 shows the snapshot of the mesh of one of the plastic plants before and after pot subtraction. This needs to be done before the 3D meshing step to produce a clean and watertight mesh. 3. 3.

3D meshing by marching cube from the remaining voxels. A grid point is checked against 8 surrounding voxels. The value of the 8 surrounding voxels is matched with a 256-element lookup table to determine if the grid is on or close to the mesh so that a polygon can be created from this grid and nearby grid points.

2.5 Mesh segmentation and geometry

Mesh segmentation algorithms involve assigning a unique label to all the vertices of the 3D mesh that belong to the same region. This paper uses a simplified version of the “hybrid” segmentation pipeline previously presented in [3]. Primarily it is based around a constrained region-growing algorithm. In short, the curvature and normal for the 3D mesh were pre-computed. A user defined curvature threshold was provided to find large “flat” regions (e.g., broad leaves) for use as seed regions. A curvature-constrained region growing was then performed from each seed region. The geometry (area, width, length, and circumference) of each segmented leaf was then extracted using the approach outlined in [3]. The result for large plastic plant is shown in Fig. 10.

3 Results

Fig. 11 shows reconstructions of different plants with different complexity and leaf shapes. 3D meshes of plants with thin and narrow leaves are reconstructed with excellent geometric agreement, although there is a minor discrepancy at the tips of the leaves. For visual comparison the pots are included in this figure, however they are removed (as in Fig. 9) before geometry measurement.

The same plastic plants are also reconstructed twice with different numbers of input images to see how this affects the reconstruction quality. Without tilting the cameras, it took 3 minutes for two cameras to capture a total of 72 images, as compared to 30 minutes to capture 360 images where the two camera were moved to 5 tilt positions (taking half of the total scanning time). A visual comparison (not shown in paper) does not show obvious differences between the two meshes of the same plant.

The leaves of the large plastic plant (bottom of Fig. 9) were dissected and scanned to measure the length, width, perimeter and area). There are 12 leaves grouped into 3 sizes shown in Fig. 12 and Tab. I. Since the leaves mostly curve along their length, the length measurement is likely to be affected. For validation of the reconstruction accuracy, the width of the leaves is chosen as this is less affected by the curving.

A quantitative comparison between the two cases is shown in Fig. 13A and B. The ground truth obtained from the 2D scans of the dissected leaves (Fig. 12) was graphed against the measurement obtained from the 3D meshes of the plant. To fit into the plots, the values of perimeter are scaled down half and the values of area are root-squared. It can be seen that the measurements agree quite well with the ground truth. The average relative error $\epsilon=\frac{1}{N}\sum\limits_{i=0}^{N-1}\frac{|truth_{i}-measurement_{i}|}{truth_{i}}$ of the measurement is 4.0% for 72 images and 3.3% for 360 images:

4 Conclusion and discussion

We have presented a complete system for automatic high-resolution 3D plant phenotyping. Several technical solutions in camera calibration, image processing and 3D reconstruction have been proposed for high accuracy of 3D mesh models. Notably, we proposed a camera calibration procedure that uses a standard chessboard calibration target that is easy to make and use in production environment. We also proposed an extension of foreground segmentation to LAB color space for improved segmentation accuracy for plants with thin leaves commonly found in major crop plants.

The system captures high quality images with accurate camera poses for image-based 3D reconstruction algorithm. The quantitative measurements using 3D visual hull algorithm provided an estimate of the accuracy of the whole system in general. We showed that useful metrics such as leaf width, length and area can be obtained with high accuracy from the 3D mesh models. Fast scanning only takes 3 minutes (72 images) per plant and still produces a reasonable measurement (4% error). More images (360 images) per plant is required for better accuracy (3.3% error) especially for complex plant structure, but requires 5 to 10 times more time to scan.

Future works include a calibration using both pan and tilt axes so that camera pose can be obtained for an arbitrary pair of pan-tilt rotation angles. This would enable a more flexible scanning trajectory other than circular rotation with fixed number of images per tilt angle.

Acknowledgments

Chuong Nguyen acknowledges the support by ARC DP120103896 and CE140100016 through ARC Centre of Excellence for Robotics Vision (http://www.roboticvision.org), CSIRO OCE Postdoctoral Scheme and the National Collaborative Infrastructure Strategy (NCRIS) project, ”Australian Plant Phenomics Centre”. Thanks to Dr. Geoff Bull of the High Resolution Plant Phenomics Centre for his valuable feedback to the manuscript.