Three-dimensional reconstruction and measurements of zebrafish larvae from high-throughput axial-view in vivo imaging

: High-throughput imaging is applied to provide observations for accurate statements on phenomena in biology and this has been successfully applied in the domain of cells, i.e. cytomics. In the domain of whole organisms, we need to take the hurdles to ensure that the imaging can be accomplished with a su ﬃ cient throughput and reproducibility. For vertebrate biology, zebraﬁsh is a popular model system for high-throughput applications. The development of the Vertebrate Automated Screening Technology (VAST BioImager), a microscope mounted system, enables the application of zebraﬁsh high-throughput screening. The VAST BioImager contains a capillary that holds a zebraﬁsh for imaging. Through the rotation of the capillary, multiple axial-views of a specimen can be acquired. For the VAST BioImager, ﬂuorescence and / or confocal microscopes are used. Quantitation of a speciﬁc signal as derived from a label in one ﬂuorescent channel requires insight in the zebraﬁsh volume to be able to normalize quantitation to volume units. However, from the setup of the VAST BioImager, a specimen volume cannot be straightforwardly derived. We present a high-throughput axial-view imaging architecture based on the VAST BioImager. We propose proﬁle-based 3D reconstruction to produce 3D volumetric representations for zebraﬁsh larvae using the axial-views. Volume and surface area can then be derived from the 3D reconstruction to obtain the shape characteristics in high-throughput measurements. In addition, we develop a calibration and a validation of our methodology. From our measurements we show that with a


Introduction
The application of high-throughput imaging is wide-spread in modern molecular genetics based biology. Imaging takes an important position in a variety of high-throughput applications and the image based applications are found on the cellular, tissue-culture and organismal level. High-throughput is a demanding process and as throughput is required, the optics setup needs to cover for sufficient image quality in further processing of the samples.
Zebrafish is a popular model system to study a broad range of biological phenomena. A lot of these phenomena are rooted in molecular genetics, i.e. developmental biology, toxicology, cancer, infectious diseases and drug targeting. In an experimental setting zebrafish are easy to breed and per crossing about 200 eggs can be obtained. The early embryonic stages are easily studied through microscopy as well as easy to manipulate. In the past decade, through genetic engineering, a large amount of transgenic zebrafish lines incorporating Green Fluorescent Protein (GFP, and the like) as a reporter gene have become available. These lines, and the genes they represent, support the research in that the gene-expression from the transgenic line helps indicating specific events in space and time. The GFP-like reporter genes are studied through fluorescent and confocal microscopy and the fluorescence facilitates measurability of the expression so that a numerical representation of experiments can be obtained.
Zebrafish, augmented with the large amount of reporter lines, are very suitable for highthroughput experiments given the number of embryos that can be obtained for experimentation. Furthermore, the genomics of zebrafish is close to human. Therefore, zebrafish are now massively used as a model in disease studies which are often performed in a high-throughput setting. Initially, these experiments were designed such that the read-out was realized using 2D images and subsequent processing of these images. In this manner a kind of mid-throughput could be accomplished. The space and time features of these images could, de facto, be used efficiently. In order to scale experiments to true high-throughput, other solutions need be probed [1]. Moreover, spatial acuity is required to enable the 3D image analysis. Apparently, confocal laser scanning microscopy (CLSM) seems the best option. So, in addition to the microscope and automated acquisition protocols, devices are required to make it possible to do accurate and fast imaging, both 2D and 3D, on life specimens and, if possible, in a time-lapse fashion. A solution has been presented through the adaptation of screening hardware from the invertebrate field and making it suitable for the vertebrate field, more specifically for zebrafish; i.e. the Vertebrate Automated Screening Technology (VAST BioImager) [2] (http://www.unionbio.com/vast/). The VAST BioImager extends a microscope, equipped with (epi-)fluorescence and/or confocal imaging options. Through a capillary system, zebrafish embryos/larvae are delivered, transported and kept in a fixed position for imaging. The capillary is mounted such that it can rotate while the specimen is kept in its fixed position. In this manner, images from different angles, i.e. axial-view images, can be acquired and analysed.
The VAST BioImager can operate with bright-field, fluorescence and confocal microscopy; a schematic of the setup is shown in Fig. 1. A system of stepper-motors is used to fix and rotate the capillary; the positioning of the specimen in the field of view of the microscope is accomplished with an extra camera which is part of the VAST BioImager system. We refer to this camera as the VAST-camera which has a mounted lens that observes the capillary via a prism. In a synchronized acquisition protocol, images from different views, e.g lateral, dorsal and ventral, are acquired with the microscope mounted camera or photomultiplier (in case of CLSM) in a full revolution of the capillary. In addition, the VAST-camera can also acquire images of the specimen, at a lower resolution, but always with a full view of the specimen. In comparison, the microscope view renders more detail, due to the characteristics of its lenses, but does not necessarily produce a complete view of the specimen. Some typical examples of the axial views obtained by the VAST-camera are depicted in Fig. 2(A).
In order to obtain successful image analysis of zebrafish, accurate and precise shape measure-

Reservoir for zebrafish
Positioning capillary Stepper motor Prism VAST-camera #1

Image acquisition management
Microscope; Microscope camera #2 Fig. 1. A schematic illustration of the axial-view imaging system based on the VAST BioImager. Zebrafish larvae are loaded from a reservoir. The VAST system detects the position and direction of the object in the capillary. The stepper motors manipulate the position and view of the specimen. Images of a corresponding axial-view can be acquired by the associated VAST-camera #1. The VAST-camera is, in particular, suitable for obtaining overview images of the specimen. Originally, this camera is concerned with the control of the system but, by us, it is adapted for axial image acquisition. In addition, the whole system is mounted on a microscope equipped with high resolution camera #2. The microscope camera enables both of the organ-and cellular-level imaging. The blue dashed curves represent the communication between the imaging components and the control platform.
ments are important, especially for providing a reference of the overall size and shape. In fluorescence, only the specific shape characteristics from specific labels can be analysed. Brightfield images, on the other hand, will give all information on shape and size if they present an overview of the specimen. Therefore, we investigate how the axial-view images can be employed to obtain precise and accurate measurements on size and shape [3]. We will use the VAST-system to produce a full revolution of the specimen and acquire, with the VAST-camera, a set of axialview images. Subsequently, we will investigate how these images can be used to reconstruct a volume and obtain the size and shape measurements. These measurements then serve as a reference of size and shape of more specific features. A general overview of the proposed method is shown in Fig. 2. Given the microscope setup with the VAST BioImager, axial-view profiles of the zebrafish as shown in Fig. 2(B) are extracted from the original images as shown in Fig.2(A) and, henceforth, a 3D volume can be produced through combination of the profile images. The underlying idea of the 3D reconstruction from axial-views is similar to the problem of multi-view based 3D vision. To that end, a camera calibration model needs to be developed as framed in Fig. 2(D). The 3D reconstruction can be accomplished once the segmentation and camera calibration are in place. An example of the 3D reconstructed zebrafish model is depicted in Fig. 2(E). We need, however, to make sure that the result is a good representation of the shape; i.e. the measurements that we derive from the model can be considered correct. So, validation experiments and evaluation of the methodology on zebrafish larvae will be performed. In order to clarify our method, the following issues are addressed in the next sections of this paper. In section 2, a number of 3D reconstruction methods are reviewed and the proposed method is motivated by the 3D reconstruction based on profiles, i.e. the binary masks obtained from segmentation of the object. In order to apply the method, an accurate estimation of the camera configuration is required. Therefore, we elaborate this for microscope optics and include some necessary optimization steps for efficient computation in section 3. In section 4, we continue with a validation of the proposed method by introducing calibration beads that are used in an initial experiment testing if our setup can produce accurate reconstructions. In addition, we test if our measurement algorithms are sound. Subsequently, we have performed experiments on a large set of zebrafish larvae and shown successful application of the proposed method for zebrafish larvae. The experiment deals with a range of developmental stages in zebrafish and the outcome of the measurements is interpreted. One of the results is a graph of the distribution for volume and surface area of the zebrafish larvae. We explain how this distribution is derived from the experimental results. We present the major conclusions from this research and extrapolate to future work in section 5.

Background and method motivation
In multi-view 3D reconstruction methods, two major categories can be discerned. One category is based on epipolar geometry and aims estimating the depth information by matching the corresponding points from correlated images of one identical scene using geometrical clues [4,5]. These methods use texture mapping to define point similarity or disparity, but require pixel-wise correspondence. The other category concerns the multi-view 3D reconstruction by producing a volumetric representation of the object from a series of binary masks of the shape, which sometimes is referred as the silhouette-based 3D reconstruction.
For the 3D measurements that we wish to extract from the zebrafish, the 3D shape plays an important role. Our method should, therefore, have an innate shape preserving ability. Due to the partial transparency of the zebrafish embryos, it is difficult to determine point correspondences from the object surface within the multi-view images. Therefore, the silhouette-based 3D reconstruction method is more suitable for this type of data. In this work, a binary mask is derived from an axial-view image which is preferably defined as a profile. As of now, we therefore use, in the remainder of this paper, the term profile-based 3D reconstruction to indicate our reconstruction method.
In the silhouette-based method, the initial idea takes into account the parallel projections derived from each individual profile [6] to generate a 3D intersection. Various 3D volumetric representations are then proposed. Hierarchical resolution voxels are used for an efficient construction of the 3D object [7]. The voxels that are populating the object surface are represented with a high resolution grid, while the voxels inside the object are represented by much coarser grid. Such data structure is commonly referred as an Octree. Due to the characteristics of camera imaging system, cone-shaped projections have been derived from the camera configurations [8,9]. The silhouette-based 3D reconstruction methods have become more feasible through the introduction of the concept of the visual hull [10][11][12]. The space carving algorithm [13,14], as conceptualized in the space-carving theory, is an important implementation of this concept. In this method, based on the pinhole camera projection model, all the silhouettes are back-projected to the 3D space represented as a voxel-space. According to space carving theory [14], only the voxels visible to each of the silhouettes are preserved and used to reconstruct the final 3D model.
Other methods are based on combining the information from texture mappings of the original images; one of which includes the surface reflectance to generate more natural 3D scene [15]. In recent approaches [16][17][18][19], the 3D reconstruction is accomplished by searching for an optimal surface. The 3D surface is embedded in a higher dimensional space. Each voxel candidate included or excluded by the enclosed surface is assigned with a probability which indicates the possibility of the voxel candidate belonging to the object or the background. The probability is modelled by a multi-variate Gaussian kernel which is applied on the texture mappings interactively indicated by users. An energy function is subsequently defined by incorporating the total probability and the total surface area. Finally, the variational framework is taken to solve the optimization problem. Recently, innovations [30, 31] were reported, which are carefully formulated to preserve and refine the details across the projective rays on the object surface by introducing various constraints, and achieve elaborate results on public dataset.
Our aim in this research is to accomplish 3D measurements from our high-throughput imaging system. From the related work, we have taken inspiration to further develop the multi-view based 3D reconstruction and introduce new ideas for the calibration of the optical setup in microscopy. Our point of departure is an existing opto-electronic configuration for high-throughput imaging to accomplish a reconstruction of zebrafish larvae. As a consequence, in order to be able to produce reconstructions from multi-view imaging, a calibration of the microscope imaging system needs be performed. To that end, a solution for the calibration problem is proposed using a revisited objective function in the optimization. This function is defined as the voxel residual volume (VRV) maximization problem, which is a simplified formulation of the area coherence [20]. In terms of implementation, VRV is more efficient compared to the silhouette coherence [21]. In order to validate the proposed method, a large dataset of zebrafish larvae of different developmental stages is produced. For each of the instances in this dataset, a multi-view image set is acquired. From our method, 3D reconstructions are produced including the generation of dense surface points and triangulated meshes [24]. From the volumetric reconstructions and the meshes, volume and surface area can be computed. In addition, the method is applied to calibrated spheres in order to obtain metrical references. With a statistical analysis of the dataset, we have obtained a metrical reference for the most frequently used developmental stages of zebrafish in high-throughput imaging. This analysis will provide a baseline evaluation for future research in high-throughput imaging of zebrafish as well as a baseline for volume and surface estimations in zebrafish development providing reference values for a range of other features.

Materials and methods
In this section, we elaborate on the imaging system, the profile-based 3D reconstruction and parameterizations that are required to realize the proposed method. Our starting point is images that are acquired with the VAST BioImager. From these images, the profiles (silhouettes) need be extracted. For the reconstruction, we need a sound camera model. In this section we will introduce and optimize this camera model.

Imaging method
We first describe the imaging system and procedures. Fig. 1 shows the imaging architecture based on the VAST BioImager as designed for zebrafish high-throughput imaging. One at a time, the specimen is loaded into a dedicated capillary which is held by a pair of stepper motors. The motors will control the capillary motion in full revolution so that any axial-view of the specimen can be acquired by the VAST-camera; a standard color CCD camera (Allied Vision Systems, Pro Silica GE 1050). Originally, the VAST-camera is used for visual detection of the location and orientation of the specimen in the capillary so that the system will be able to manipulate the position and perspective of the specimen and initialize the imaging. Due to the characteristics of VAST BioImager, the imaging can produce the an overview image of the specimen; in fact, it can do so very efficiently for a large amount of specimen. This is very suitable for our application of 3D reconstruction and measurements of the specimen. Therefore, we adopted the system to be the basis of our 3D reconstructions. It should be noted that the VAST BioImager is mounted on a microscope, of which the high resolution camera will be able to produce more detailed images of the specimen. In combination with the VAST-camera, the further analysis of cellular-, organand overview-level of the specimen will become possible. In our method, we require the binary representations of the object segmented from the color images acquired from the VAST BioImager. In order to obtain a good and solid reconstruction, these binary representations should reflect the whole subject; in the case of transparent zebrafish, we have to ensure that the translucent parts are also included in the shape representation. We concluded that standard segmentation methods are not sufficient and therefore we have implemented a hybrid method that incorporates the mean shift algorithm (MS) [23] and an improved level-set method (ILS) [32]. From Fig. 3, one can observe that the the MS preserves a whole shape representation for the object but fails in edge sensitivity for the original clear boundaries. The ILS can obtain a more compact 2D shape but is hampered by the problem of edge leakage. So, we combine the two segmentation methods followed by a refinement to obtain  Fig. 4. A visualization of the parameterization of 3D transformation from camera center to object center. We assume that the object is positioned in the focal plane and a good quality image is acquired, so the projection line from camera center to the object center is defined as the focal length. The ∠α and ∠γ represent the 3D rotation angles of camera center around Y and Z axis. The ∠ϕ is defined as the "translation" angle from the object center to the image center. It models the 3D translations of camera center over Y and Z axis. accurate segmentations of the object, i.e. the zebrafish. For the rest of the paper, the results from this segmentation method will be the basis for the 3D reconstructions. The segmentation of a set of axial views will be referred to as the profiles. At http://bio-imaging.liacs.nl/ galleries/VAST-3Dimg/ a link is provided with a further elaboration on this segmentation method.

Camera model parameterization
Key to the profile-based 3D reconstruction is a feasible and sufficiently accurate parameterization of the camera imaging system. For our particular application, the camera position is static while the object rotates through the revolution of the capillary that holds the object, over a given profile axis. In order to make parameter parameterization and visualization feasible, the camera is represented as a range of fixed cameras around the object in a circular path as shown in Fig. 2(D). Now, for each of these cameras, the pinhole imaging principle can be used to interpret the imaging procedure. We start from a straightforward case of the camera model in which case the camera center is defined as the world origin. Let X = (X, Y, Z, 1) T be a point in 3D space and x = (x, y, 1) T be the corresponding point in a 2D image plane, represented in homogeneous coordinates. Then, the mapping between the point and the corresponding image can be geometrically represented as a proportional projection of x = ( f X )/Z and y = ( f Y )/Z, where f denotes the focal length of the camera lens. The scaling factors k x and k y , image center (u x , u y ) T and the skew factor s are, de facto, the intrinsic properties of the camera. Taking these intrinsic properties into consideration, the matrix representation of the projection can then be formulated as: wherex = xZ. To simplify the formulation, an implicit transform of Eq. 1 is denoted asx = PX, where P = K (I | 0), and P is defined as the camera projection matrix and K is defined as the intrinsic camera matrix. In practice, the more generic case for the camera model is to define the object center as the world origin, which is more convenient for computation and parameterization. The decomposed parameterization for the camera model allows convenient incorporation of the extrinsic camera parameters which are the 3D transformation of the camera center with respect to the world center. Suppose that for one object we have N views, the intrinsic camera matrix K is shared by all the cameras, however, a corresponding 3D transformation for each individual view should be formulated. This is illustrated in Fig. 2(C) and Fig. 4. For the i th camera, the 3D translation of the camera center from the object center is parametrized as the total 3D rotation, the revised formulation of the camera projection matrix is then represented by The "translation" angle of ϕ and the rotation angle of γ around Z axis are specified by each of the cameras. In this manner we intend solving the problem that the rotation axis is not exactly parallel to the object center. In practice, with the VAST BioImager one can accurately align the object, so these two angles can be shared by all N cameras. Finally, a camera parameter vector is generated through the concatenation of all parameters.
Given a specific camera parameter vector ψ, a series of correlated camera configurations can be generated, so that the multi-view 3D projections P can be modelled.

Profile-based zebrafish 3D reconstruction
We aim to recover the volumetric representation of the object from a range of 2D axial-view images. The profile-based 3D reconstruction method shows to be suitable for shape preservation in this application. Given a series of (object) profiles extracted from the axial-view images, combined with the corresponding camera projection matrices, the 3D representation can be generated by back-projecting the profiles to the 3D world frame so as to obtain the intersections.
Let the 3D point X j ∈ X ( j ∈ [1, M]) be the index of the j th voxel candidate center included by a collection of voxels, where M is the total number of the voxels; if the collection of voxels is spatially constrained, then M will determine the 3D volumetric resolution. The larger the M is, the higher the resolution will be. Next, one can find the corresponding pixel of the j th voxel in the i th image plane by x j i = P i X j . According to the space-carving theorem, if the voxel X j is visible to the view of i th camera, its corresponding image should be covered by the i th profile S i . In mathematical terms this can be formulated as: where X i denotes the volumetric representation for the object generated from the i th camera. When involving all the views, we can locate the pixel location in each image plane for each voxel.
In general, only the voxels which are visible to all of the N profiles are preserved. Thus, where X * denotes the optimal 3D volumetric representation of the object. Eq. 4 can be efficiently implemented by iteratively discarding the invisible voxels from the original voxel collection. However, taking the eventual imperfection in the profiles segmentation into consideration, some voxels which should be included in the 3D object are excluded during this space carving operation. Instead, we assign a score to each of the voxels indicating its visibility. This score is formalized as where 1[·] denotes an indicator function which takes the value of 1 if the image of the j th voxel lies in the i th profile and 0 otherwise. The score function assigns a value between 0 and N for each individual voxel. If, in addition, a constrained error tolerance rate is defined, the problem of imperfection in the profile segmentation can be more or less solved by thresholding the membership function of a voxel. The 3D reconstruction model of Eq. 4 is, therefore, reparameterized accordingly as: For our particular application in zebrafish imaging, we have empirically established that the parameter takes the value of 0.05, meaning the voxels visible to 95% of the profiles will contribute to the 3D reconstruction. With the reconstructed 3D volumetric representation, an initial 3D surface model can be obtained by the marching cubes algorithm; this surface model is represented as dense surface points and triangulated meshes, which can be further optimized [24].

Camera system optimization
From the previous sections it becomes clear that the camera projection matrix, i.e. the camera configuration, plays an important role in the profile-based 3D reconstruction. A good estimation of the camera configuration will ensure accuracy in the 3D projective geometry. If a set of specific camera parameters as denoted in Eq. 2 are given, the camera projection matrix P can be correspondingly computed. The camera parameters, more specifically, the extrinsic camera properties, are usually unknown or not sufficiently accurate; e.g. through a drift of the image center from the object center. Therefore, an automated camera configuration estimation, i.e. a camera calibration, should be performed. The importance of a camera calibration is demonstrated in Fig. 5 by the results of applying the method to 3D reconstruction of zebrafish shapes. Standard camera calibration methods [28, 29] for our imaging environment are not available for the small scales of microscope imaging, i.e. the VAST-camera. Therefore we have proposed [3] a method similar to ideas on texture registration and stitching [20]. We define the area coherence as the area of overlap between the image projected from the 3D reconstructed model and the groundtruth profile. Accordingly, we can define an energy function as: where the camera matrix is reparameterized as P i (ψ) meaning that P is a function of ψ. By evaluating the camera configuration space, the optimal suggestion of ψ can be obtained by maximizing Eq. 7 Visualizations corresponding with this approach are depicted in the lower rows of Fig. 5(A) and (B). We argue, however, that the expression of Eq. 7 is equivalent to the, so called, voxel residual volume (VRV) which is defined as Equation 9 denotes the total number of the voxel candidates of the 3D volumetric representation. Similarly, the reparameterization of X * (ψ) indicates the functional relation between X * and ψ. In principle, keeping the maximal total area of overlap between the projected images and the groundtruth profiles is actually equivalent with keeping the maximal VRV of the reconstructed object. The VRV is, however, computationally less complex while the same performance can be obtained. This is because the computation of the area coherence are unnecessary. For such an unconstrained optimization problem, multiple optimizers are available, for example, the Nelder-Mead simplex [25] and the Evolution Strategy [26,27]. One should be aware that, although small, there always exists the risk of a local minimum of the optimizers which may lead to an inaccurate camera configuration estimation. However, feasible initializations are possible so as to provide accelerations and improvements to the optimization. The intrinsic camera parameters including the focal length and the scaling factors are provided by the camera and lens specifications. The extrinsic parameters including the 3D rotations and translations can be assigned through reasonable estimations; e.g. the parameter of ω i is initialized as ω i = (i − 1) * 2π/N. With these sensible initializations, more efficient and accurate camera configurations become possible.

Experimental results
In this section we explore the method through validation experiments. These experiments are performed on calibration particles as well as on zebrafish larvae. This will result in numerical and statistical validity of the method so that it can be successfully applied in further research.

Sampling of axial-views and volume for the experiments
With respect to the sampling, the physical limitations imposed by the hardware and the imaging scale are confined. There are two sampling modes that co-operate in the reconstruction process. The production of axial-view images is depending on the properties of the hardware, which is referred to as the Axial-view Sampling Density (ASD) and determined by the step size of the stepper motors that operate the capillary. For the experiments we will use a range of steps in order to find a good operational ASD from sparse to dense. For the physical reconstruction, a volumetric representation is required. From the hardware and the spatial constraints of specimens, a confined volume of 5000 × 1000 × 1000 µm 3 (width × height × depth) is determined. The sampling involves the variation of the size of the isotropic elements (voxels) in this volume, which can be inferred as the total number of the elements. The smallest sampling size will be equal to the sampling elements on the CCD sensor in which case the number of sampling elements is 30 × 10 6 ; this represents the largest number of sampling elements. The smallest number of elements is 1 × 10 6 which gives rise to a considerably different performance. The given sampling volume is used in all experiments described in this paper.

Validation of the proposed method
In order to validate the accuracy of the proposed method, particles, i.e. beads, with a known size distribution are used (GP 500 µm Control Particles, High Fluorescence, lot 110701, Union Biometrica). The particles are used to calibrate the VAST BioImager system. Some typical examples are depicted in Fig. 6(A), from which we can observe some variation in the diameter of the particles. We assume that this variation in size is according to a Gaussian distribution as given by the manufacturer, i.e. N (500, 25) measured in µm. We have acquired image sets from the VAST BioImager for 25 calibration particles. Subsequently, the proposed 3D reconstruction method is applied to these image sets.
If the size measured from the reconstructions corresponds with the distribution of the size of the original calibration particles, it means that the proposed method can accurately recover the 3D object in terms of size and shape. To measure the diameter of the reconstructed calibration particles, a 3D sphere is estimated to fit the dense surface points by minimizing a least squares function. The result of the profile-based 3D reconstruction method can be affected by the Axialview Sampling Density (ASD) and the voxel size (sampling density). As far as the ASD is concerned, more views may result in a more accurate and natural shape, similarly a higher sampling density of the voxels may result in a more accurate representation of the shape. Therefore, the experiments are designed to verify the accuracy according to these two parameters as well as establishing the correctness of the measurement itself by comparing to analytical representations or empirical observations. The results are presented in Table 1 -3; here the sampling density is denoted as SD given the size of the isotropic sampling element (in µm). For all tables the horizontal direction, from left to right, indicates an increasing sampling density from 17.10 µm per voxel to 5.85 µm per voxel. The vertical direction indicates an increase of the ASD from 4 to 84 views. We have configured our imaging system to take 84 images, corresponding with a step-size of 4.3 • per axial-view, in a full revolution of the shape. Initial experiments have shown that a strategy of even axial-view sampling in a full revolution accomplishes the most stable performance. By resampling of the 84 views of 4.3 • , a range of different ASD's is obtained for the experiments.
Evaluations for diameter of the calibration particles We interpret the measurements of the diameter in terms of the mean and standard deviation over the 25 calibration particles and present the results in Table 1. The results show that from the proposed method accurate estimations are obtained which are in range of the distribution (500 ± 25 µm) provided by the manufacturer. The results vary with different ASD and SD. More specifically, an increase in the sampling density results in an increase of the estimated diameter which asymptotically stabilizes at the higher sampling densities. A high sampling density generates a smoother 3D surface, while for lower resolutions the volumetric representation tends to be a little bit expanded. An increase of the ASD results in an decrease of the estimated diameter while for higher ASD it stabilizes. Using more views to reconstruct the 3D shapes causes less voxels to be preserved and the 3D shapes seem to be more compact compared to the ones reconstructed by smaller ASD. We tested whether our measurements were sufficiently accurate. In Table 1, in each second row, the T-test scores are given per ASD. We state in the null hypothesis that the measured diameter from the 3D reconstructed models equals to the values (mean, sigma) given by the manufacturer, i.e. 500 µm. For a significance level of 0.01 (two-sided T-test), a value of 2.787 is given from the table of selected values. From Table 1 one can appreciate that most T-test scores are smaller than the selected value; this means that our measurements for the diameter of the calibration particles from the 3D reconstructions are accurate. Especially, the values from the ASD of 21 combined an SD of 6.93 or smaller come out very well.
Evaluations for volume and surface area of the calibration particles The calibration particles have been provided with an indication of the size distribution, i.e. given diameter. However, the real distributions of the volume and surface area are unknown. In Table 2 and 3 the volume and surface area of the calibration particles are presented as computed from the reconstructed shapes in a range of parameters for ASD and SD. Volumes are computed directly from the voxels in the shape, i.e. integrating over the residual voxels after reconstruction. For each object a surface model is composed from a dense surface point cloud and a mesh triangulation from these points. The surface area can be obtained by integrating over the facets of the meshes using Heron's formula [24] A = (s(s − a)(s − b)(s − c)) (1/2) , where s = (a + b + c)/2 and a, b, c represent the edge-length of a triangle. We assume that the shape of the calibration particles is a standard sphere, consequently the diameters from Table 1 are used to obtain the analytical computation of the volume and surface area by V = (4/3) * π (d/2) 3 and S = 4π (d/2) 2 for each set of parameters. For each ASD, the values measured from the reconstructions are listed in each first row of Table 2 and 3, indicated with M. The second row for each ASD in Table 2 and 3 lists the analytical computations, indicated with an A. The trends in Table 2 (Volume) and 3 (Surface area) are similar to those observed in Table 1. While the ASD increases, the reconstructed shapes are smoother and less voxels are preserved. Moreover, the variation in the results decreases at higher SD with voxels smaller than than 6.93 µm.
Taking the reconstruction effects shown in Fig. 6(C) into account, it can be appreciated that a ASD of 21 is a reasonable choice to generate a smooth and natural shape with very few carving artifacts. Moreover, one can observe very small differences between the measured and analytical values for an ASD of 21 as of the SD of 6.93; i.e. the shape of the calibration particle is well preserved and reconstructed with such set of parameters. In comparison, a larger value of ASD requires more camera projections and the 3D reconstruction is therefore computationally more expensive; a sampling density smaller than 6.93 µm requires more memory and computation to ensure the precision. So, the results in this experiment indicate that the proposed method can obtain accurate 3D reconstructions measured in size, volume and surface area; the ASD of 21

Evaluations of performance on application with zebrafish larvae
In order to obtain 3D reconstructions and evaluate the performance of the proposed method as applied to zebrafish specimens, we first assemble a large number of image sets of zebrafish larvae with the help of the axial-view imaging architecture presented in section 3.1 and put the images in a dataset. We have chosen the most frequently used larval stages (for VAST BioImager) to be represented in the dataset; meaning 3 groups, i.e. 3, 4 and 5 days post fertilization (dpf), containing 12, 24 and 24 instances, respectively. Although the VAST BioImager is suitable for live imaging, for this experiment, we have used fixed samples of one and the same strain, i.e. ABxTL wild type. All samples were fixed in 4% PFA/PBS (PFA: paraformaldehyde; PBS: Phosphate-Buffered saline). For the imaging, the samples were kept in a solution of PBS with 0.3% Tween-80. For each specimen, the same imaging setup of the VAST BioImager is consistently used (cf. calibration experiment configuration). Specifically, in one full revolution, 84 views are captured with equal step size. Other ASDs, e.g. 42, 21 etc., are evenly sampled from the these 84 views. In the acquisition, the rotation steps between any two contiguous views are considered to be equal. However, due to the drift on the electro-mechanical parts of the VAST imager, we have to take into account that there can be a variation in the step size. If present, this variation can seriously affect the 3D reconstruction of the zebrafish specimens and therefore a calibration of the camera system is, therefore, very necessary.
In Table 4 and 5, we list the volume and surface area for the three groups of larval stages. For these tables the horizontal direction indicates an increase in ASD. In order to get a better idea of the the data, we visualized Table 4 and 5 in Fig. 7 and 8, in which the voxel residual volume (VRV) and the surface area are plotted as a function of the ASD, respectively. For all three groups the trend is similar, as the ASD increases the VRV decreases. This can be explained from the methodology as only voxels which are visible to most of the views are preserved. Consequently, less voxels will be preserved when more views are employed. The decrease of VRV rapidly changes in an asymptotic fashion of VRV for an ASD larger than 7. The asymptote indicates that we do not need an abundance of sampling views, an ASD of 42 or 84 does not result in a much better reconstruction and therefore these are unnecessary for the reconstruction of the zebrafish larvae in our application. The VRV changes with the growth of the zebrafish for the different larval stages; i.e. as the zebrafish grow, the VRV increases. Similar trends can also be found in Fig. 8 which depicts the asymptotics of surface area for the 3D zebrafish models. External conditions influence growth as it also depends on the feeding policy. At some point in development the larvae have used all their yolk and are thereafter depending on nutrients in the medium. So, if the amount of nutrients is low, a delayed growth will occur. From Fig. 7 it can be concluded that an ASD of 7 is already a good option for this application. However, integrating the results from all experiments, shows that an ASD of 21 achieves the best performance measured in terms of size, volume and area. In order to appreciate this, the reconstructed results are visualized for an identical zebrafish specimen in Fig. 9. For a ASD of 4 or 7, the 3D models exhibit lots of artificialities. The space carving effects generate some flat regions on the surface. However, with a higher ASD value, such as 12 or 14, the 3D reconstruction results improve; the surface of the shape is still rough though. As of an ASD of 21 or larger, one can appreciate that results are more or less similar; the 3D shapes all appear smooth and natural. The visualization quality does not improve much for an ASD larger than 21, while a large ASD is computationally much more expensive. The empirical ASD of 21 is therefore a motivated tradeoff over a range of considerations regarding the 3D reconstruction of zebrafish larvae and the measurement thereof.
A number of typical examples of models from 3D reconstructed zebrafish larvae are shown in Fig. 9. For each zebrafish developmental stage in our experiment two individuals have been selected; the 3D shapes are visualized with and without texture mapping. The texture mappings are just serving visualization purposes, where texture is used to add some realism to the view. The specimens, in this stage, are still partially transparent, therefore, some artificial mapping was applied to obtain a realistic visualization. For measurement the binary volume models are used (depicted in green).

3D measurements on zebrafish larvae
From the aforementioned experiments we can conclude that an accurate reconstruction can be obtained from our method and from the 3D reconstruction accurate measurements can be made. In order to assess the shape and size in a high-throughput setup we use the database of our samples to find the distributions per larval stage. From our empirically established parameters we reconstruct our images to 3D models using a ASD of 21 and a voxel sampling density of 6.93 µm.
On these models, the computations are performed, the volume and surface area expressed in µm 3 and µm 2 respectively. Over the measurements per group, the mean and standard derivation are calculated and a Gaussian distribution is applied to model the data. In order to compare the 3 groups of larval stages the density distributions are normalized. The distribution of the volume is shown in Fig. 11. The trend as seen in Fig 7 is reflected in the results in that the volume of 5 dpf zebrafish is the largest, and that of 3 dpf zebrafish is the smallest. The volume is more discriminative for the 3 dpf and 5 dpf zebrafish. The 3 dpf and 4 dpf zebrafish appear to be very similar in volume, which can be seen from the large overlap of the two distributions. Indeed, these larval stages are very similar in shape and appearance and this is also illustrated in Fig. 10.
The distribution of surface area is shown in Fig. 12 and can be interpreted in similar fashion to the volume measurement. A plot for the the joint distribution of volume and surface area is shown in Fig. 13. Although the distributions of 4 and 5 dpf zebrafish are close to each other we can clearly separate the 3 centres. This experiment provides significant quantitative information on zebrafish larvae because from the distributions, given the larval stage, the regularity of the shape can be directly assessed. The distributions provide an absolute reference baseline of volume and surface area for zebrafish in 3 larval stages. This is very important for applications in toxicology and drug targeting in which quantities of specific substances need to be related to volume and shape and phenotypical information on dose and effect. These need to be derived from measurements on the specimen. Moreover, as the readout of effects of substances is often measured with fluorescent probes, a relation to volume is even more important. In addition to these basic measurements other 3D features can be considered, e.g. moments, wavelets, shape descriptors, etc. This will further strengthen the phenotypical description and assessment of zebrafish larvae in high-throughput imaging applications. In order to make sure that we could measure the specific characteristics of the given larval stages, we used fixed specimens. In a control experiment we have used living larvae staged by specialists and imaged with the VAST BioImager under exactly the same conditions. We used 7 examples from 4 dpf and 8 examples from 5 dpf. The results of the measurements are shown in Table 6. From these results one can appreciate that they are all within the range of the fixed specimens.
So, we validated our methodology for the 3D reconstruction and measurements and at the same time showed the importance of the baseline measurements for the 3 larval stages; the major stages in compound screens with zebrafish. The reference values for volume and surface area as well as their distributions for the different larval stages can be directly used in experiments. Consequently, for phenotype analysis an assessment for relative size and shape can be made. Moreover, features that are frequently used for phenotyping zebrafish and thereby evaluating anomalies can be related to size and shape to better assess the effect under study.  Fig. 8. Surface area of zebrafish larval stages against ASD. The trend in surface area is similar to that of the volume as depicted in Fig. 7. The surface area is computed from a mesh, to suppress noise a mesh smoothing, i.e. implicit fairing, with 10 iterations was applied to all meshes [33].

Conclusions and future work
An automated system is presented in which images of zebrafish larvae are acquired and used to reconstruct to 3D models so that 3D measurements can be applied on the whole specimen.  and an ASD=14 (D) results in an improvement of the 3D models, but some carving artifacts remain. Using an ASD=21(E) to an ASD=84 (H) results in accurate and natural-shaped 3D models. For our particular problem domain, this suggests that both sparse and dense axial-view sampling are not the optimal configuration for zebrafish larva 3D reconstruction. Sparse axial-view sampling produces poor 3D models while dense axial-views sampling requires more computation time whilst the effects do not improve accordingly.
The system integrates the axial-view imaging, segmentation, camera system optimization and a profile-based 3D reconstruction method. The 3D measurements of volume and surface area, Fig. 10. Visualization of 3D models of 3 zebrafish larval stages (3 dpf, 4 dpf and 5 dpf). Each box represents a reconstructed 3D model for one specific zebrafish larvae visualized from three different viewpoints. The 3D volumetric representations are shown in green on the left side in each box. The models with texture-mapping are in the right side of each box. (A) 3D models of two selected 3 dpf zebrafish larvae. (B) 3D models of two selected 4 dpf zebrafish larvae. (C) 3D models of two selected 5 dpf zebrafish larvae. Variation in size and shape between stages and within stages (interclass and intraclass) can be appreciated from the visualizations. A remarkable intraclass discrimination originates from the size and color of the yolk. In addition, animations of the 3D zebrafish models are available at: http://bio-imaging.liacs.nl/galleries/VAST-3Dimg/.
as well as the 3D visualization become available through the proposed method on zebrafish larvae using high-throughput imaging. From the research presented in this paper, the important conclusions are that (1) accurate 3D measurements for the zebrafish larvae can be made in high-throughput; (2) as of now a metrical reference of shape descriptors for 3 important larval stages is produced and (3) the methodology for these measurements is validated. The results are directly applicable for a range of different in vivo applications, such as toxicology, drug targeting and infection studies. Further research will be directed to 3D shape alignment so that an even better match between different zebrafish instances can be given. Subsequently, other 3D features will be explored to accomplish elaborate classification of subtle differences between different treatment groups of zebrafish in experiments. The VAST-camera provides good overview images which have been employed for this study. The same methodology can be applied with the microscope camera. The higher numerical aperture and the magnification of the lenses in the microscope will then provide 3D details of the specimen in similar fashion. The tradeoff of using a confocal or wide field approach for such applications needs to be assessed from future results. Alternative to the axial-reconstruction with the VAST BioImager, Optical Projection Tomography (OPT) imaging can be probed. This imaging approach has been successfully applied in zebrafish research [34] and, recently, successful high-throughput applications in zebrafish have been reported [35]. As indicated, the measurements and their distributions provide reference values for volume and surface area (size and shape). This supports the assessment of other observations as we now  Fig. 11. Distribution of volume of zebrafish larval stages (3 -5 dpf). X axis: volume; Y axis: normalized probability density. The color-filled triangles on the Volume-axis indicate the locations of the mean for each of the 3 distributions. The numerical values of the mean and standard deviation are indicated with corresponding double-sided colored arrows. One can see that the volume of older zebrafish larvae is always larger than the younger ones. The growth rate increases more from 4 dpf to 5 dpf compared to that from 3 dpf to 4 dpf.  can relate the effect of a change in size to a norm. In like manner, the VAST system can be used to provide distributions of other morphometrical landmarks in zebrafish development; to this end the microscope can be used to acquire axial-views at higher resolution and, in that manner, provide 3D features. These features can be combined with volume/surface area features to enable  Fig. 13. Joint distribution of volume and surface area of zebrafish larval stages (3 -5 dpf). X axis: volume; Y axis: surface area. The 3 joint distributions have overlap with respect to each other, especially for the 4 dpf zebrafish larvae. Nevertheless, individual distribution centers can still be separated. The color schema is similar to that of Fig. 11 and Fig. 12 relative comparisons between experiments. Such analysis will further enhance our understanding of phenotypes in zebrafish and their exposure to experimental conditions which is a necessity for high-throughput analysis.

Funding
This research was funded by the Netherlands Organisation for Scientific Research (NWO), grant number 834.14.001. Yuanhao Guo was financially supported through the China Scholarship Council (CSC) to participate in the PhD programme of Leiden University.