3D Calibration using OpenCV

The calibration was performed using the Python bindings of the OpenCV [Bradski, 2000] library. This has been made easier to use through the module fluidimage.calibration.calib_cv. We shall use this module on a set of 5 calibration images of a target which has a circle grid. The white dots on this particular target is evenly spaced at 3 mm. For the same camera position, the target coordinate z varies as [-6, -3, 0, 3, 6] mm.

We shall proceed as follows:

1. We compose a function find_origin which automatically detects the origin in pixel coordinates of the calibration target. This is achieved by using an erosion operation to fill the faint rectangle in the origin, and then using OpenCV to detect the location of this blob (origin) of minimum area 18.

2. After this we detect the image points, i.e. smaller circles in a 7x7 grid surrounding the origin and store them in an array. We repeat this operation for every calibration image.

3. We construct he object points, i.e. assign the circle grid the expected values in the world coordinate system (x, y, z) and store them as arrays using the input given to us that the circles on the target are evenly spaced by a distance equal to 3 mm.

4. Finally we calibrate the camera.

OpenCV employs a camera model based on the algorithm following Zhang [2000].

Let us start by loading a set of calibration images.

from fluidimage import get_path_image_samples

path = get_path_image_samples() / "TomoPIV" / "calibration" / "cam0"
calib_files = sorted(path.glob("*.tif"))
[path.name for path in calib_files]

['-3mm_cam0.tif',
 '-6mm_cam0.tif',
 '0mm_cam0.tif',
 '3mm_cam0.tif',
 '6mm_cam0.tif']

Detecting the origin

A typical calibration image looks like:

%matplotlib inline
import matplotlib.pyplot as plt
from fluidimage.util.util import imread

def imshow(image, ax=plt):
    ax.imshow(image, cmap='gray', vmax=255)

image = imread(str(calib_files[2]))  # z = 0 image
imshow(image)

The position of the origin (marked by a rectangle) needs to be detected for detecting the image points consistently.

from fluidimage.util.util import imread
import matplotlib.pyplot as plt
import numpy as np
from skimage.morphology import reconstruction
import warnings

def rescale(image, scale):
    """Rescale the image intensities"""
    # Scale 16-bit image to 8-bit image
    if scale is None:
        return image
    elif scale == "median":
        scale = np.median(image[image > 5])
    elif scale == "max":
        scale = image.max()

    # print("Rescaling with", scale)
    image = image * (256 / scale)
    return image

def imfill(filename):
    """Fill boundaries in an image. This is used to make the origin easy to detect."""
    image = imread(filename)
    image = rescale(image, "median")
    # Fill the rectangle at the center. Helps to detect the origin
    seed = np.copy(image)
    seed[1:-1, 1:-1] = image.max()
    mask = image
    with warnings.catch_warnings():
        warnings.simplefilter("ignore")
        # Fill the square in the center to mark the origin
        image = reconstruction(seed, mask, method='erosion')

    return image.astype(np.uint8)


test_calib_file = str(calib_files[2])
fig, axes = plt.subplots(1, 2, dpi=100)
imshow(imread(test_calib_file), axes[0])
imshow(imfill(test_calib_file), axes[1])

To detect the origin we use SimpleCircleGrid class. Although we intend to detect only one point, it works by tweaking the minArea parameter. This class will be described in the next section.

from fluidimage.calibration.calib_cv import SimpleCircleGrid
import os

def find_origin(filename):
    params = SimpleCircleGrid.create_default_params()
    params.minArea = 18.

    circle_grid = SimpleCircleGrid(params)

    # Pass the filled image to detect the square
    keypoints = circle_grid.detect_all(
        imfill(str(filename)))
    assert len(keypoints) == 1
    return keypoints[0].pt

for cfile in calib_files:
    print(f"Origin of {cfile.name.rjust(13)} detected at", find_origin(cfile))

Origin of -3mm_cam0.tif detected at (280.0, 250.0)
Origin of -6mm_cam0.tif detected at (309.0, 250.0)
Origin of  0mm_cam0.tif detected at (250.0, 250.0)
Origin of  3mm_cam0.tif detected at (220.0, 250.0)
Origin of  6mm_cam0.tif detected at (190.0, 250.0)

Detecting image points as a circle grid

The result is a list of blobs in pixel coordinates, centers in image coordinates.

from fluidimage.calibration.calib_cv import SimpleCircleGrid

params = SimpleCircleGrid.create_default_params()
params

<fluidimage.calibration.calib_cv.ParamContainerCV object at 0x7fb7258b25d0>

<SimpleBlobDetector blobColor="0" collectContours="False" filterByArea="True"
                    filterByCircularity="False" filterByColor="False"
                    filterByConvexity="True" filterByInertia="True"
                    maxArea="5000.0" maxCircularity="3.4028234663852886e+38"
                    maxConvexity="3.4028234663852886e+38"
                    maxInertiaRatio="3.4028234663852886e+38"
                    maxThreshold="220.0" minArea="0.10000000149011612"
                    minCircularity="0.800000011920929"
                    minConvexity="0.949999988079071" minDistBetweenBlobs="10.0"
                    minInertiaRatio="0.10000000149011612" minRepeatability="2"
                    minThreshold="50.0" thresholdStep="10.0"/>  

There are certain parameters which can be tweaked to detect the circles as needed. For this particular case the defaults are enough.

def construct_image_points(filename, debug=False):
    image = imread(str(filename))
    image = rescale(image, "max")
    origin = find_origin(filename)
    if debug:
        print("Origin =", origin)

    params = SimpleCircleGrid.create_default_params()
    circle_grid = SimpleCircleGrid(params)
    centers = circle_grid.detect_grid(
        image, origin, nx=7, ny=7, ds=50, debug=debug)

    return centers

centers = construct_image_points(calib_files[2], debug=True)

Origin = (250.0, 250.0)
Bbox(x0=50.0, y0=50.0, x1=450.0, y1=450.0)

Object Points

The calibrate function requires objectPoints (world coordinates) and imagePoints (image coordinates) of the blobs detected.

from fluidimage.calibration.calib_cv import construct_object_points

For example

construct_object_points(nx=3, ny=3, z=-1, ds=3)

array([[-3., -3., -1.],
       [ 0., -3., -1.],
       [ 3., -3., -1.],
       [-3.,  0., -1.],
       [ 0.,  0., -1.],
       [ 3.,  0., -1.],
       [-3.,  3., -1.],
       [ 0.,  3., -1.],
       [ 3.,  3., -1.]], dtype=float32)

Calibration

We now put together all the elements above to calibrate

from pathlib import Path
from tempfile import gettempdir
from fluidimage.calibration.calib_cv import CalibCV


path_output = Path(gettempdir()) / "fluidimage_opencv_calib"


def calibrate_camera(cam="cam0", debug=False):
    path_calib_h5 = path_output / (cam + ".h5")
    calib = CalibCV(path_calib_h5)

    objpoints = []
    imgpoints = []
    zs = []

    path = get_path_image_samples() / "TomoPIV" / "calibration" / cam
    files = sorted(list(path.glob("*.tif")))

    # Populate objpoints, imgpoints and zs
    for i, filename in enumerate(files):
        z = int(filename.name.split("mm_")[0])
        zs.append(z)
        objpoints.append(
            construct_object_points(nx=7, ny=7, z=z, ds=3)
        )
        centers = construct_image_points(str(filename))
        imgpoints.append(centers)

    im_shape = imread(str(filename)).shape[::-1]
    origin = find_origin(str(files[2]))
    return calib.calibrate(imgpoints, objpoints, zs, im_shape, origin, debug)


ret, mtx, dist, rvecs, tvecs = calibrate_camera("cam0", debug=True)

from pprint import pformat

def print_horizontally(vecs):
    vecs2 = ['']
    vecs2.extend([v.T for v in vecs])
    return pformat(vecs2)

print(f"""
        Avg. reprojection error = {ret}
        fx, fy = {mtx[0,0]}, {mtx[1,1]}
        cx, cy = {mtx[0,2]}, {mtx[1,2]}
        k1, k2, p1, p2, k3 = {dist.T}
        rotation vectors = {print_horizontally(rvecs)}
        translation vectors = {print_horizontally(tvecs)}
""")

        Avg. reprojection error = 10.975788795409809
        fx, fy = 1.0, 1.0
        cx, cy = 250.0, 250.0
        k1, k2, p1, p2, k3 = [[-5.15599720e-17 -1.35151238e-12  1.09663198e-26 -7.29811031e-21
   2.72357380e-17]]
        rotation vectors = ['',
 array([[-3.81871080e-11, -5.93309022e-05,  6.85117975e-09]]),
 array([[-5.33619016e-11, -1.49309655e-04, -1.15509541e-09]]),
 array([[-1.92346303e-11,  3.31239794e-05, -1.41583220e-09]]),
 array([[-4.94477664e-12,  1.29400604e-04,  3.82576448e-09]]),
 array([[ 4.03427685e-13,  2.12863876e-04, -3.82124696e-09]])]
        translation vectors = ['',
 array([[ 1.65673059e+00, -4.05216175e-08,  3.05369360e+00]]),
 array([[ 3.30335071e+00, -1.61159751e-08,  6.05367695e+00]]),
 array([[-9.88774057e-04,  2.13526518e-08,  5.35866173e-02]]),
 array([[-1.66530159e+00, -4.44864037e-08, -2.94664309e+00]]),
 array([[-3.31183376e+00,  2.39950118e-12, -5.94700736e+00]])]

Calibrate all 4 cameras

for i in range(4):
    calibrate_camera(f"cam{i}")

list(path_output.glob("*.h5"))

[PosixPath('/tmp/fluidimage_opencv_calib/cam2.h5'),
 PosixPath('/tmp/fluidimage_opencv_calib/cam3.h5'),
 PosixPath('/tmp/fluidimage_opencv_calib/cam0.h5'),
 PosixPath('/tmp/fluidimage_opencv_calib/cam1.h5')]