Calibration of multiple cameras for large-scale experiments using a freely moving calibration target

Obtaining accurate experimental data from Lagrangian tracking and tomographic velocimetry requires an accurate camera calibration consistent over multiple views. Established calibration procedures are often challenging to implement when the length scale of the measurement volume exceeds that of a typical laboratory experiment. Here, we combine tools developed in computer vision and non-linear camera mappings used in experimental fluid mechanics, to successfully calibrate a four-camera setup that is imaging inside a large tank of dimensions ∼10×25×6m3\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\sim 10 \times 25 \times 6 \; \mathrm {m}^3$$\end{document}. The calibration procedure uses a planar checkerboard that is arbitrarily positioned at unknown locations and orientations. The method can be applied to any number of cameras. The parameters of the calibration yields direct estimates of the positions and orientations of the four cameras as well as the focal lengths of the lenses. These parameters are used to assess the quality of the calibration. The calibration allows us to perform accurate and consistent linear ray-tracing, which we use to triangulate and track fish inside the large tank. An open-source implementation of the calibration in Matlab is available.


Introduction
New studies in biophysics and fluid mechanics require the quantitative imaging of large-scale field experiments. Such studies include the large-scale Lagrangian tracking of bats and bird flocks (Theriault et al. 2014;Attanasi et al. 2015), super-large-scale particle image velocimetry measurements using natural snowfall (Toloui et al. 2014) and recent advancements in tomographic-PIV (Jux et al. 2018).
Obtaining reliable two-and three-dimensional imaging data in these large-field experiments is challenging and requires a camera calibration that is accurate down to the smallest physical length scale of interest. Non-linear polynomial camera mappings (Soloff et al. 1997;Wieneke 2005Wieneke , 2008 are often used in laboratory experiments (Kühn et al. 2010), but their application at length scales beyond that of the laboratory is, in practice, limited. First, the size of the calibration target is limited, such that it only covers a small portion of the measurement. Second, the absence of conventional laboratory equipment, providing access to the measurement volume, does not support the accurate spatial positioning of the target, required in conventional calibration procedures.
In the present work, we combine the pinhole camera model (Tsai 1987) with non-linear polynomial camera mappings used in experimental fluid mechanics (Soloff et al. 1997) to perform a multiple camera calibration over a large-scale measurement volume inside the tank of the aquarium located in the Rotterdam zoo. Our method integrates the use of the pinhole camera model with a nonlinear camera mapping to correct for optical distortion across refractive interfaces (Belden 2013). Our approach uses the framework of projective geometry in computer vision (Hartley and Zisserman 2004) and apply advanced self-calibration techniques (Svoboda et al. 2005;Shen et al. 2008).
Here we apply the planar checkerboard calibration technique by Zhang (2000), see also Zhang (1998Zhang ( , 1999, Sturm and Maybank (1999), Menudet et al. (2008) and Bouguet (2015). This approach eliminates the need to accurately position the calibration target, as required in conventional calibration procedures. Instead, the checkerboard calibration target is moved to arbitrary and unknown positions and orientations, here with the help of a team of divers. Second, by sequentially acquiring multiple calibration images while freely moving the calibration target, we achieve a camera calibration that spans over length scales much larger than the calibration target itself. This approach yields an accurate calibration over the measurement volume with a characteristic length scale on the order of several tens of meters.
We process the camera calibration in steps (Heikkilä and Silvén 1997). First, we correct for optical distortions (Fryer and Brown 1986) by rectifying the curved lines of the checkerboard images (Prescott and McLean 1997;Devernay and Faugeras 2001). Second, we perform a calibration based on a single view for each camera following (Zhang 2000). Finally, we combine the single views and find the positions and orientations of the cameras over multiple views (Geiger et al. 2012) and optimize the calibration for spatial accuracy and consistency between the different views.
The camera calibration yields accurate results. To assess the validity of the camera calibration, we compare the estimated effective focal length obtained from it against the true value of the focal length of the lenses. We quantify the spatial accuracy of the camera calibration, by computing the skewness of optical rays associated with the multiple views. Our calibration allows us to use linear ray-tracing (Hedrick 2008) to track and triangulate multiple fish swimming over the entire visual depth of the tank. The method is versatile and can be implemented in field experiments over large length scales and for measurement volumes that are challenging to access experimentally.

Camera setup and calibration procedure
We image inside the large tank of the aquarium in the Rotterdam zoo in a measurement volume of dimensions ∼ 10 × 25 × 6 m 3 (Fig. 1). We use four 5.5 megapixel sCMOS cameras with wide-angle lenses of focal length f lens = 24 mm (Nikkor AF-24 mm ) that cause significant variations in magnification over the large depth-of-field of DOF ≈ 25 m (Adrian and Westerweel 2011).
The camera setup is positioned behind an acrylic window of thickness ∼ 50 cm . Optical distortions in the image plane are due to refraction at the water ( n water = 1.363 ) /acrylic ( n acryl = 1.51 ) and the acrylic/air ( n air = 1.0003 ) interfaces (Sedlazeck and Koch 2012), where n is the refractive index. The optical access is limited to the acrylic window, which constrains the spacing between the cameras to H ≈ 1 m in height and W ≈ 6 m in width, and limits the relative angles between the cameras from 5 • to 20 • (Fig. 1).
To calibrate the camera setup we image a planar checkerboard calibration target of dimensions 1.5 × 1.8 m 2 with 5 × 6 tiles (Geiger et al. 2012) that are of area A tile = 30 × 30 cm 2 . This calibration target is moved within the aquarium by a team of divers that swim with the checkerboard under arbitrary and unknown positions and orientations throughout the aquarium.

Image processing
The different images of the calibration target are processed to identify the curves and the nodes corresponding to the gridlines and intersections between the tiles of the checkerboard (Fig. 2). In our application, using a checkerboard calibration target is advantageous over using a pattern of dots. This is because the image gradient obtained from a checkerboard determines grid points more accurately and more robustly over the large depth-of-field of our experiment.
The calibration images are converted to the image gradient using a Savitsky-Golay image differentiation approach (Meer and Weiss 1992) to mark locations on the gridlines between the tiles. For each image, we then fit a set of polynomial curves i (t) to the I = 9 gridlines between the tiles of the checkerboards using the local intensity values from the image gradient. These fitting curves are written as i (t) = ∑ k k i t k−1 , where k i are two-dimensional vectors. Here the parameter t varies within the interval t ∈ [0, 1] over the checkerboard image, such that i (t = 0) and i (t = 1) correspond to the beginning and end-points of the gridline of the checkerboard image, see the lines on Fig. 2. At last, we find the J = 20 intersections between all the gridlines as a set of nodes j in the image plane of each camera. Here j = [x j y j ] T with (•) T the vector transpose, where the numbering j = 1 ⋯ J is consistent between the different camera views (Geiger et al. 2012) (Fig. 2).

Distortion correction
Following the path of an optical ray for a single camera in Fig. 2, the linear path is refracted across the air/water interface (Belden 2013). This causes optical distortions in the image plane for each camera and, therefore, we rectify the image plane by dewarping the optical distortion. The coordinates = [x y] T in the image plane are mapped to distortion-corrected image coordinates ̂ = [xŷ] T by determining a distortion mapping ̂ = m( ) . This mapping ensures that collinear points in the object domain are projected as collinear points in the dewarped image plane and, therefore, supports linear ray tracing. For moderate optical distortions, a polynomial distortion map (Soloff et al. 1997) is sufficient. In this study, we write the distortion map as The distortion map is determined by rectifying the curved gridlines in the N original calibration images. We follow Devernay and Faugeras (2001) and minimize the percentage of deflection along the gridlines. We consider the nodes n j (1) = + ∑ k k k ( ) = + 1 x 2 + 2 xy + 3 y 2 + 4 x 3 + 4 x 2 y + 5 xy 2 + 6 y 3 . Fig. 1 The four-camera setup at the Rotterdam zoo. a Schematic representation of the large measurement volume of ∼ 10 × 25 × 6 m 3 along the DOF ≈ 25 m , including the flat checkerboard of dimensions 1.5 × 1.8 m 2 that is positioned by a team of divers. b Optical access through the large window spanning 2 × 8 m 2 including the camera setup at relative spacing of H ≈ 1 m by W ≈ 6 m . c An example of calibration images acquired while positioning the checkerboard 7 Page 4 of 12 along a particular gridline n i (t) and their images ̂ n j and ̂n i (t) in the dewarped image plane through the map m( ) . We fit a straight line ̂ n i through the nodes ̂ n j and compute the point-line distance d(̂ n j ,̂ n i ) . The parameters k defining the distortion map are then determined by solving a minimization problem over all gridlines in all calibration images: This minimization problem can be solved efficiently using a steepest descend algorithm and numeric integration techniques described in Boyd and Vandenberghe (2004). This approach directly extends to larger optical distortions requiring more elaborate distortion models (see Appendix A for the air/water interface and Supplementary- Fig. 8). An example of a dewarped image can be found in Fig. 2. (2)

Camera calibration and projective geometry
We consider a physical point in the object domain of coordinates = [X Y Z] T in a world coordinate system and its projected image ̂ = [xŷ] T in the dewarped image plane of a single camera. The calibration is defined by the mapping function F, such that ̂ = F( ) . Our method uses the framework of projective geometry to express the mapping function F and implicitly assumes a pinhole camera model. In the following, we outline the main notations used in projective geometry; a more complete introduction can be found in Hartley and Zisserman (2004). We make use of augmented vectors to represent points in both the image plane and in the object domain. The coordinates in the dewarped image plane ̂ are augmented to the ray-tracing vector ̃ such that ̃ = [kx kŷ k] T , where k is a scaling parameter in the direction of the principal optical axis. The associated inverse function that projects ̃ back to ̂ is defined as the projection p(̃ ) = [xỹ] T ∕k =̂ . Similarly, the world coordinates are augmented to a homogeneous vector as ̃ = [X Y Z 1] T (Hartley and Zisserman 2004). Using augmented vectors a geometric transformation, consisting of a rotation and a translation, is simply written as a matrix multiplication [R ]̃ , where R is a rotation matrix and is a translation vector.
With these notations, the mapping function can be written in the following form that is widely used in projective geometry: (Hartley and Zisserman 2004): Here K is the 3 × 3 camera calibration matrix and [R ] is a 3 × 4 matrix, with R the 3 × 3 rotation matrix, and the 3 × 1 translation vector. The matrix K in Eq. 3 has the following form: where x = fr x and y = fr y are scale factors, with f the focal length of the lens in mm and r x and r y the pixel pitch of the sCMOS sensor in (px∕mm) . s is the pixel skew, characterizing the angle between the x and the y pixel axes, and [p x p y ] T are the coordinates of the principal point at the intersection between the optical axis and the dewarped image plane (Hartley and Zisserman 2004). The elements of K are often referred to as the intrinsic camera parameters, representing the characteristic properties of the camera itself (Hartley and Zisserman 2004), while [R ] are referred to as the extrinsic parameters representing the position of the camera with respect to the world coordinate system (Fig. 2).
Together, K, R, define the mapping function F of Eq. 3 and have to be determined for each of the cameras separately.
x y x ŷX c Fig. 2 Image processing and the geometry of the optical path. a Processed calibration image with the identified gridlines fitted as secondorder polynomial curves (green lines) and the nodes at the gridline intersections (red squares). The gridlines are numbered from i = 1 to 9 and the nodes from j = 1 to 20. b The dewarped calibration image corresponding to a in which the gridlines are rectified for minor optical distortion using the mapping of Eq. 1. Supplementary- Fig. 8 provides an example of dewarping for severe optical distortion. c The geometry of the optical path where the optical rays are refracted across the air/water interface (neglecting the acrylic window for simplicity) and the positioning of the (virtual-)camera coordinate system [xỹ k] T with the origin at the camera center c in world coordinate system [X Y Z] T . d The local coordinate system [X o Y o ] T attached to the plane of the checkerboard. The indexing i corresponds to the gridlines and j to the locations of the nodes In the following, we use the superscript c = 1 ⋯ 4 and the notations K c , R c and c , when we distinguish explicitly between the different cameras. We omit the superscript for clarity when no distinction between the cameras is needed; the details of the algorithms can be found in Appendices.

Single-camera calibration
First, the camera matrix K is determined for each of the four separate cameras by calibrating a single camera using the method developed by Zhang (2000). We consider a local to the planar checkerboard in the object domain, where Z o = 0 corresponds to the plane of the checkerboard. In this coordinate system, the J nodes at the intersections between the gridlines have known The nodes o j are mapped to their images ̂ n j following the formalism used in Eq. 3. This transformation can be written as p(K[R n n ]̃ ) o j , where the rotation matrix R n and translation vector n characterize the position of the checkerboard in the object domain. Geometrically, this corresponds to a rotation and translation of the checkerboard plane in the object domain, followed by a projection on the image plane. The camera calibration matrix K and the positioning of the checkerboard, characterized by R n and n , are determined following Zhang (2000) and are refined by minimizing the following functional:

Multiple-camera calibration
To complete the calibration over the multiple cameras, we determine the rotation R c and the translation c representing the position of each of the cameras in the world coordinate system. R c and c correspond to the extrinsic camera parameters described in Sect. 2.3 and a first estimate is deduced directly from the calibration of single cameras performed in the previous step. Selecting two different calibrated cameras, we use the Kabsch algorithm (Kabsch 1976) to estimate the relative positions of the two cameras by comparing the positions of the checkerboards, R n,c and n,c , for these two views. Considering all camera pairs, we determine the relative positions between all the views and deduce a first estimate for R c and c ; see Appendix B for detail. Next, for each of the N calibration images, we estimate position R n and n of the checkerboard in the world coordinate system. We do this by averaging the position estimates R n,c and n,c obtained from the four separate single-camera calibrations.
In the last step, we compute the final values for K c , R c and c by minimizing the reprojection errors from all cameras and calibration images. For this optimization step, we use as initial conditions, the camera matrices K c obtained from Sect. 2.4, the positions R c and c obtained from the Kabsch algorithm and the estimates of the checkerboard positions R n and n . We define the reprojection error n,c j for each of the N calibration images and in each camera view as The reprojection error n,c j is normalized by Â n,c j , which is the area in pixels of a single tile on the checkerboard projection in the dewarped image plane (see Appendix C). Hence, n,c j provides a measure of the error relative to the size of the checkerboard. This error is relatively independent of the location of the calibration target within the large depth of field and the positions of the cameras and, therefore, independent of the apparent size of the checkerboard in the image plane. The final parameters for the calibration function F for each view are obtained from the minimization of the following summation: The resulting camera calibration is shown in Fig. 3.

Assessment of the calibration method
In the following, we evaluate the performance of the calibration. First, we report the extrinsic and intrinsic camera parameters from Table 1 and discuss their physical interpretation. Second, we assess the robustness and convergence of the method as a function of the number of calibration images used. Third, we study the spatial accuracy of the calibration and identify the sources of error.

Intrinsic and extrinsic camera parameters
First, we consider the numerical values of the extrinsic camera parameters obtained from a calibration, see Table 1. As discussed in Sect. 2.3, these parameters characterize the spatial position and orientation of each camera. The reconstructed positions of the cameras are in agreement with the experimental setup, with cameras 4 and 1 positioned above cameras 3 and 2, and cameras 4 and 3 positioned on the left-hand side while cameras 1 and 2 are located on the right-hand side (as shown in Fig. 1), see Table 1. Furthermore, the relative distances between cameras are also in agreement with the experimental scene as we find the horizontal distance between cameras W ≈ 5.6 m and a vertical distance between cameras H ≈ 1.2 m , as deduced from Table 1. Likewise, the reconstructed camera orientations are consistent with the cameras on the right-hand side oriented with a positive angle , while the cameras on the left-hand side are oriented with a comparable negative angle .
In addition to reconstructing the position of the camera, the calibration procedure reconstructs the intrinsic camera properties, which we compare to the specification of the instrumentation. The coefficients of the camera calibration matrix K of Eq. 4 are provided for each camera in Table 1. We focus on the values of the focal length of the lenses and deduce an effective focal length f eff directly from the coefficients as where r is the resolution of the camera sensor in px/mm , which is known from the camera specifications, J represents a correction factor for the image expansion due to optical distortion (see Appendix D) and ñ = n air ∕n water corrects for the magnification due to refraction at the air/water interface. Our calibration yields values for the effective focal length of the four cameras of f eff = 23.73 ± 0.82 mm . This reconstruction of the focal length lies within 1-2 % of the actual focal length f lens = 24 mm of the lenses that were used. Hence, we find that both the extrinsic and intrinsic camera parameters deduced from our calibration procedure are in agreement with the dimensions and characteristics of our experimental setup.

Convergence and robustness
The camera calibration is obtained by minimizing the sum of the squared reprojection errors n,c j over all four cameras, all N calibration images and all nodes on the checkerboard, see Eq. 7. The camera calibration converges to low values for n,c j see Table 1. Here, the average normalized reprojection error is on the order of n,c j ∼ 2% of the size of a checkerboard tile, which corresponds to an error of less than 1 cm.
The camera calibration requires a minimum of two noncoplanar checkerboard images (Zhang 2000). Increasing the number of calibration images increases the sampling Fig. 3 The resulting camera calibration. a The calibrated views that include physical scales on a reference plane at k = 1 m in the depth of field for each view. b Three-dimensional reconstruction of a random selection of checkerboards used for the calibration, in yellow the reconstructed (virtual) cameras in front of the acrylic window of Fig. 1

Table 1 Numerical values of the calibration parameters
From top to bottom: the expansions factor of the distortion mapping J , the intrinsic camera properties from the matrix K, the extrinsic camera positions X, Y and Z and orientations in pitch-yaw-roll angles , and (Figs. 1, 2). The effective focal lengths f eff deduced from K with Eq. 8, the reprojection errors n,c j in percentage and equivalent pixel-dimensions (•) * , and the total number of calibration images N used for each camera of the measurement volume and, therefore, improves the reliability of the calibration. We further characterize the performance of the method as a function of the number of calibration images used, by randomly selecting different subsets of checkerboard images. We first consider the effective focal length f eff of Eq. 8 deduced from the matrix K to characterize the quality of the position in space of the checkerboards and the cameras. In Fig. 4, we select different subsets of calibration images ranging from N = 2 to N = 50 images and compute f eff from the associated calibration. Figure 4a represents the ratio between f eff and f lens as a function of N. With a low number of N = 2 to 10 calibration images, the ratio f eff ∕f lens is far from one and the focal length is under-and over-estimated by 50% , indicating an unreliable calibration. Increasing the number of calibration images to N = 15 shows that the estimated focal length f eff approaches f lens and represents a clear improvement of the calibration. Further increasing N beyond N = 15 , the ratio f eff ∕f lens does not significantly converge further (Fig. 4a).
Second, we consider the reprojection errors n,c j in Fig. 4b. For a low number of N = 2 to 10 calibration images, the average reprojection error is as high as n,c j ∼ 50 to 60 % . Increasing the number of calibration images from N = 15 to 50 shows an additional decrease of the normalized reprojection error from n,c j ∼ 5 to 2 % (Fig. 4b) and inset. This shows that 15 calibration images are sufficient to achieve a valid calibration. Further increasing the number of calibration images improves the convergence for the camera calibration while the ratio f eff ∕f lens remains at a value close to 1.

Spatial accuracy of the camera calibration
By inverting the camera matrix K, one can directly associate an optical ray in the object domain to a point in the dewarped image plane. For an ideal calibration, the four optical rays associated with the images of the same point on each of the four cameras should intersect at a unique location in the object domain. In practice, the four optical rays are skew lines and do not intersect at a single point. Here, we characterize the spatial accuracy of our camera calibration by estimating the skewness among the four optical rays.
For this, we use the nodes identified at gridline intersections on the N calibration images. We proceed by evaluating the four optical rays associated with each node of each calibration image. We then triangulate the location of each node by finding the point n j in the object domain that minimizes the sum of the squared distances from the point to the four optical rays. We report for each node the skewness s n j as the average distance from the triangulated location n j to the four optical rays (see Appendix E for more detail).
The calibration images were acquired over the entire depth of the tank and used to characterize the spatial accuracy, by reporting the skewness as a function of the position along the Z-axis of the world coordinate system (Fig. 5a). We find the skewness to be mostly uniform over the large depth of the measurement volume Z = 5 − 25 m . Our calibration yields a high spatial accuracy characterized by an average skewness of less than 1 cm . Only a slight increase in the skewness can be observed towards the back of the aquarium although the average skewness still remains below 1 cm . This small increase is due to the decrease in spatial resolution and the decrease in angles between the optical path from the different views.
The variation in skewness over the height and width of the tank are represented in Fig. 5b, c. Figure 5b is a map of the skewness in a XY-plane at the front of the aquarium, while Fig. 5c provides a similar map at the back of the aquarium. Both are qualitatively similar and the skewness remains small over the depth of the tank. For reference, Fig. 5d, e shows the sampling density, which indicates the number of checkerboard images that were used at a location to compute the calibration. It is noteworthy that the center of the tank was better sampled than the sides and the bottom. We find that the skewness, and hence the spatial accuracy, is relatively constant over a large part of the measurement volume, but increases towards the edges of the tank. This is likely a result of the lower sampling away from the camera

Application to field experiments
We demonstrate the effectiveness of the calibration method by performing three-dimensional measurements in the large aquarium at the Rotterdam Zoo. We begin by focusing on an element which is easily identifiable on each camera view and show that we can accurately triangulate the position. We end our validation by tracking the three-dimensional position of fish of various sizes that are freely swimming over the depth of the tank.
Large predatory fish in the aquarium, such as sharks, swim through the entire aquarium. They provide a good target to evaluate the robustness of the camera calibration as their sharply tipped fins provide easily recognizable and well-defined features. Figure 6 shows how accurately such features can be triangulated with our calibration method. We first identify the vertex of the right-hand side pectoral fin of a shark on each camera view. Similar to Sect. 3.3, we directly associate the vertices, identified on each of the four images, with four optical rays in the object domain. We triangulate the position of the vertex and find a skewness of only 0.35 cm , which demonstrates the accuracy of the method. This small skewness is illustrated in Fig. 6, where, for each camera view, the optical rays associated with the other camera views are projected on the image plane onto the epipolar lines (Hartley and Zisserman 2004). The epipolar lines intersect nearly perfectly at the identified vertex of the shark fin. The inset in Fig. 6 presents a closeup view from which one can infer the reprojection error from the slight skew between the optical rays. This associated reprojection error is of only 1.11 px.
Further, we demonstrate that our calibration supports the tracking of several fish simultaneously over large distances by tracking a small group of six tuna fish (Fig. 7). By triangulation, we reconstruct the three-dimensional time-resolved position of the group (Fig. 7b). Using an in-house automated tracking code, we track the group of six individual fish (tuna) The sampling density of checkerboards associated with (b, c), respectively Fig. 6 Triangulation of the vertex of a shark fin. For each camera view: the point corresponding to the vertex of the shark fin is identified with a marker, while the three lines correspond to the epipolar lines associated with the three markers of the remaining three camera views. The color coding is consistent across the multiple views, for example, the vertex on camera 1 is identified with a red marker and the epipolar lines associated with this point in cameras 2, 3, 4 are red. Likewise the vertex in cameras 2, 3, 4 and the corresponding epipolar lines are respectively represented in green, blue and magenta. The inset on camera 2 zooms on the vertex of the fin and shows that the epipolar lines intersect at pixel accuracy swimming away from the cameras over a large distance of 7 m . Together with the triangulated shark fin, this experiment demonstrates the robustness and accuracy of the calibration and its potential use in large-scale field experiments. The calibration supports accurate triangulation over a large distance along the depth of field. This makes it of interest to further application for the tracking of particles (Schanz et al. 2015), birds (Attanasi et al. 2015), insects (van der Vaart et al. 2019), fish and other animals, and the study of fluid motion using tomographic methods (Elsinga et al. 2006) for large-scale field experiments.

Conclusion
Here, we have described and characterized a versatile, accurate and robust calibration method, which supports the three-dimensional tracking and triangulation of multiple fish. Our method is of particular interest to large-scale fields experiments, when spatial access to the measurement volume is limited and laboratory equipment to precisely position the target cannot be installed. The method does not require a large calibration target to be moved with known displacements. Rather, we use a freely moving checkerboard calibration target, which is much smaller than the measurement volume itself. The calibration mapping uses the framework of projective geometry, which assumes linear ray tracing. It combines a pinhole camera model for the linear ray tracing, with a non-linear camera mappings commonly used in experimental fluid mechanics to correct for optical distortion. All the algorithms necessary for the implementation of the calibration method are described here with details provided in Appendices.
The calibration method has been implemented to obtain an accurate and consistent multiple-view camera calibration in the large-scale aquarium of the Rotterdam Zoo. Here, the calibration target was positioned arbitrarily by a team of divers. We have characterized the spatial accuracy of the calibration to be less than 2% of the size of a checkerboard tile, corresponding to 1 cm over the entire measurement volume that spans several tens of meters. When correcting for the difference in refractive index of air and seawater, we find that our method reconstructs both the camera position and the intrinsic properties of the camera such as the focal length of the lenses. Selecting different subsets of calibration images in a quality assessment, we show that in our case 15 calibration images are sufficient to achieve a valid calibration. Increasing the number of checkerboard images and better sampling the measurement volume further improves the spatial accuracy. The resulting camera calibration allows accurate imaging and threedimensional tracking over a large measurement volume, here with application to biological fluid mechanics and the tracking of fish.

A: Optical distortion across an interface
Refraction across an interface between two media of different refractive index is governed by Snell's law: where ñ = n 1 ∕n 2 with n the refractive index, 1 is the incident angle and 2 is the exit angle to the normal vector on the interface. For a flat interface, the image plane of a camera can be warped parallel to the interface by the linear image mapping H: where A is a 2 × 2 matrix that describes an affine image transformation that changes the aspect ratio and skew of the image, centers the image at the principle ray, and T describes a perspective change (Hartley and Zisserman 2004). With these notations, the distortion mapping can be written as where = (1 −ñ 2 )∕f with f the focal length in pixels dimensions.

B: Relative camera positioning from calibrated views
We determine the rigid body motion from a view c to another view c ′ , by first finding the best rotation, using the Kabsch algorithm (Kabsch 1976). We first compute the cross-covariance matrix A = ∑ n,j n,c � j ( n,c j ) T using the paired object coordinates n,c j and n,c ′ j from multiple checkerboards. We then compute the singular value decomposition of the cross-covariance matrix A = USV T and extract the best rotation matrix R c,c ′ as Knowing R c,c ′ , we use the rigid body motion n,c � j = R c,c � n,c j + c,c � to compute the translation vector c,c ′ between the views.
From the relative position between the views we find a unique extrinsic camera positioning R c , c by taking the first view at reference and solving the minimization problem for the translation c , and the rotation R c , where F is the Frobenius norm that sums over all matrix components squared and I is the identity matrix. These linear least squares problems have a direct solution using methods described in Boyd and Vandenberghe (2004) and generalize to any number of views (Fig. 8).