Target-less and targeted multi-camera color calibration

(1)

Independent degree project - first cycle

Datateknik

Computer Engineering

Target-less and targeted multi-camera color calibration Emil Hammarlund

(2)

Examiner: Ulf Jennehag, Ulf.Jennehag@miun.se

Supervisor: Dima Elijs, Mittuniversitetet, Elijs.Dima@miun.se Author: Emil Hammarlund, emha1503@student.miun.se

Degree programme: Degree of Master of Science in Engineering: Computer Engineering, 300 credits

(3)

Abstract

Multiple camera arrays are beginning to see more widespread use in a variety of different applications, be it for research purposes or for enhancing the viewing experience in entertainment. However, when using multiple cameras the images produced are often not color consistent due to a variety of different reasons such as differences in lighting, chip-level differences e.t.c. To address this there exists a multitude of different color calibration algorithms. This paper ex- amines two different color calibration algorithms one targeted and one target- less. Both methods were implemented in Python using the libraries OpenCV, Matplotlib, and NumPy. Once the algorithms had been implemented, they were evaluated based on two metrics; color range homogeneity and color accuracy to target values. The targeted color calibration algorithm was more effective improving the color accuracy to ground truth then the target-less color calibration algorithm, but the target-less algorithm deteriorated the color range homogeneity less than the targeted color calibration algorithm. After both methods where tested, an improvement of the targeted color calibration algorithm was attempted. The resulting images were then evaluated based on the same two criteria as before, the modified version of the targeted color calibration algorithm performed better than the original targeted algorithm with respect to color range homogeneity while still maintaining a similar level of performance with respect to color accuracy to ground truth as before. Further- more, when the color range homogeneity of the modified targeted algorithm was compared with the color range homogeneity of the target-less algorithm.

The performance of the modified targeted algorithm performed similarly to the target-less algorithm. Based on these results, it was concluded that the targeted color calibration was superior to the target-less algorithm.

(4)

Foreword

I would like to begin with personally thanking my supervisor Dima Elijs, during the many meetings that we have had during the course of this project he has always provided fantastic guidance and discussions that have been some of the most fruitful learning experiences I have experienced so far. I would also like to thank Rasmus Holm for the many great brews of coffee he has provided, and Joel Nilson for the insight he has brought during our discussions surrounding the subject that this project has covered. Last but not least, I would like to thank my loving family and friends that have provided me with the moral support that has kept me motivated throughout the project.

(5)

Terminology

FLANN Short for ”Fast Library for Approximate Nearest Neighbor”. Approx- imates a search for the matching neighbor nearest neighbor for matching feature points [1]. 22–24

RANSAC Short for ”Random sample consensus”. RANSAC is a resampling technique that generates candidate solutions to a set of data while trying to use a minimum amount of observations. The RANSAC algortihm is, within the context of this project, more specifically used for finding matching points between images[2, 3]. 22–24

SIFT Short for ”Scale Invariant Feature Transform. Is used for extracting dis- tinctive features from images [4]. 3, 22, 23

(8)

1 Introduction

1.1 Background and Problem motivation

Multiple camera arrays are beginning to see more widespread use both for research purposes as well as for enhancing the viewing experience in the entertainment industry. Researchers from the university of Hong Kong have used multiple camera arrays for image based rendering of ancient chinese artifacts [5]. GoPro, an action camera manufacturer has used multiple camera arrays in order to create ”bullet time video” [6].

When using multiple camera arrays, the images produced are often not color consistent. This problem arises due to a variety of different reasons be it chip level differences, varying intensities of light from different light sources etc. This poses a problem, especially when it comes to light field rendering. To solve this issue there are a multitude of different color calibration algorithms.

Among these methods there are two main ways of approaching the problem.

Targeted color calibration and target-less color calibration. Targeted color calibration uses a color pattern board as its calibration source. Target-less color calibration extracts color relation information from the image without relying on a color pattern board and uses this information as its calibration source.[7]

There are advantages and disadvantages to both targeted and target-less color calibration. Targeted color calibration has the advantage of being able to include the entire spectrum of colors in one image due to it’s use of a color pattern board, while target-less color calibration only can calibrate the colors that are supplied by the image [7]. Yet the use of a color pattern board can be limiting in some cases, for instance if there is a light source with a high intensity the targeted color calibration will not take the area with higher intensity into account. In this case target-less color calibration will have the upper hand since it extracts its color relation information from the whole image.

1.2 Overall purpose

The purpose of this thesis is to implement and compare two already existing targeted and target-less color calibration algorithms. From this, this thesis aims to develop insight into the performance of targeted and target-less color calibration algorithms and possibly verify the advantages and disadvantages to targeted and target-less color calibration algorithms.

Both methods will be tested with data supplied from the Realistic 3D group’s camera system. The results from the color calibration shall be compared between the two methods with respect to color range homogeneity and color accuracy to target values. After the comparison between the two already existing methods is complete further improvements will be implemented and tested.

(9)

1.3 Limitations

The process of color calibrating cameras is done in two steps, the first step involves adjusting camera parameters such as gamma, brightness, gain etc. and the second step involves modifying the captured images/video. This thesis is limited to the second step in color calibration.

Only two different color calibration algorithms will be implemented, one targeted and one target-less. These will not be selected based on which avail- able algorithms have the best performance necessarily, rather they will be selected based on feasibility of implementation and how well documented the method is in earlier work.

This thesis will not look into a wide variety of test data or cameras. Test data will be limited to a smaller subset of data and a limited variety of cameras.

1.4 Concrete and verifiable goals

In order to execute this thesis successfully the following goals must be achieved.

1. Implement a targeted color calibration algorithm 2. Implement a target-less color calibration algorithm

3. Test and evaluate the algorithms with respect to ground truth

4. Test and evaluate the algorithms with respect to color range homogeneity 5. Improve further upon the targeted color calibration algorithm

6. Improve further upon the target-less color calibration algorithm

1.5 Overview

Following this section, chapter 2 discusses related work surrounding color calibration. Chapter 3 explains the methodology used in this thesis, chapter 4 will cover the construction of the two chosen methods, chapter 5 presents the results and finally chapter 6 discusses the conclusions that can be derived from the results.

(10)

2 Related work

This thesis focuses on the work from three articles: ”Color Calibration for Multi-camera Imaging Systems” written by Sabri Gubruz, Masahiro Kawakita and Hiroshi Ando[8], ”Ensuring Color Consistency across Multiple Cameras”

by Adrian Ilie and Greg Welch [9] and ”Image Registration Using Adaptive Polar Transform” written by Rittavee Matungka, Yuan F. Zheng and Robert L.

Ewing [11]. The papers ”Color Calibration for Multi-camera Imaging Systems”

and ”Ensuring Color Consistency across Multiple Cameras” explain the two color calibration methods that will be examined. ”Image Registration Using Adaptive Polar Transform” explains a sub process used in ”Color Calibration for Multi-camera Imaging Systems” that is used for verifying matching feature points. The following sub chapters summarize their work.

2.1 Color Calibration for Multi-camera Imaging Systems

The paper ”Color Calibration for Multi-camer Imaging Systems” written by Sabri Gubruz, Masahiro Kawakita and Hiroshi Ando [8] utilizes a set of stereo correspondences from a reference camera to achieve a target-less color calibration. This sub section explains their work.

The extraction of color relation information in [8] is done in four steps.

The first step involves utilizing Scale Invarient Feature Transform (SIFT) to locate initial stereo correspondences [4]. The second step further refines these correspondences using a polar transform verification algorithm by Rittavee Matungka, Yuan F. Zheng and Robert L. Ewing [11]. In the third stage a fundamental matrix F is calculated from the refined correspondences between con- secutive camera pairs. This fundamental matrix is then used in the epipolar constraint:

m_lF m^t_r= 0. (1)

To verify their their stereo geometry. In equation 1, mland mr are stereo correspondences from the left and right image. In the final stage a color discrepancy function for the remaining stereo correspondences is calculated using equation 2:

∆E =p

(Rl− Rr)²+ (Gl− Gr)²+ (Bl− Br)². (2) In which Rl, Gl and Bl are the Red, Green and Blue components of the correspondences from the left image and Rr, Gr and Br are the Red, Green and Blue components of the correspondences from the right image. The dis- tribution of ∆E is then modeled after a Gaussian probability density function N (µ, σ). Correspondences that are ∆E > kσ are discarded in order to avoid outliers (where 1 < k < 3).

For each time instance, that is each frame in the video, a 3x4 matrix for color

(11)



 R_k−1 Gk−1

Bk−1



= [M_k|t_k]





 R_k Gk

Bk

1







, (3)

where (Rk−1, Gk−1, Bk−1)and (Rk, Gk, Bk)are RGB color values from the same object points from camera number k and camera number k − 1. Mk is a 3x3 coefficient matrix multiplier. The variable tk is a 3x1 translation vector

t_RtGtB

T

that is used for compensation for the offset between camera number kand camera number k − 1. Mkand tkare both estimated using linear regression from the set of correspondences between camera k and camera k − 1.

With Mk and tkestimated, each pixel in image number k can then be calibrated using equation 4:



 R⁰_k G⁰_k B_k⁰



= [Mk|tk]





 R_k G_k B_k 1







. (4)

In which R_k⁰, G⁰_kand B⁰_kare the new RGB values for camera number k.

(12)

2.2 Adaptive Polar Transform Verification

Adaptive polar transform verification is performed with a set of matching points as it’s input and outputs a set of matching points that pass the verification process. This verification for two given matching points is performed in three steps.

1. Adaptive polar transform the two given points.

2. Calculate their radius and angle projections 3. Calculate the distance coefficient between them

Once the distance coefficient between the two images has been calculated, the mean and standard deviation of each of the projections is calculated. Those matches that have a distance coefficient εz> |µl+ σl− (µr+ σr)|are discarded.

2.2.1 Adaptive Polar transform

The paper ”Image Registration Using Adaptive Polar Transform” by Rittavee M, Yuan F. Z and Robert L. E explains adaptive polar transform [11]. This subsection explains their work.

Let nrand nθdenote the numbers of samples in the radius and the angular directions respectively. For an adaptive polar transform given a square of size 2 · Rmax× 2 · Rmax. The parameter nrshould be greater then or equal to Rmax. Let Ridenote the size of the radius in pixels for sample i in the radius direction and Ui be the sampling circumference at sample i. Therefore, Ri = ^i·r_n^max In order to sample the image effectively nθneeds to be adaptive. Since eachr

circumference Ui is approximately 2πRi pixels, the number of samples for a given nθishould take this into account. For simplicity’s sake, 2πRiis estimated to be approximately 8 · Ri. Given this, the sampling parameters nrand nθican be denoted as seen in equation 5,

nr= Rmax, nθ_i= 8 · Ri. (5) The complete algorithm for transforming a square image of shape 2·Rmax× 2 · Rmax pixels in the Cartesian to an adaptive polar images can be seen in algorithm 1.

Algorithm 1Algorithm for adaptive polar transforming an image for i ← 1 : nrdo

nθ_i← 8 · i for j ← 1 : nθ_ido

IP (i, j) ← I(Rmax+ Ricos(^2πj_n

θi), Rmax+ Risin(^2πj_n

θi)) end for

(13)

In algorithm 1, IP is the adaptive polar transformed image and I is the Cartesian image that is being adaptive polar transformed.

Algorithm 1 can be further generalized to adaptive polar transform any given point of an image with any given radius. Let nrbe equal to the desired sampling radius, x be the desired sampling point in the x direction and y be the desired sampling point int the y direction. The sampling points must all be within the bounds of the image. This generalization can be seen in algorithm 2.

Algorithm 2Generalization of algorithm 1 for i ← 1 : nrdo

n_θ_i← 8 · i for j ← 1 : nθido

IP (i, j) ← I(x + Ricos(^2πj_n

θi), y + Risin(^2πj_n

θi)) end for

end for

2.2.2 Calculation of radius and angular projection

In order to keep the advantages of scale and rotation invariance from a con- ventional log polar transform, the paper ”Image Registration Using Adaptive Polar Transform” by Rittavee M, Yuan F. Z and Robert L. E uses an innova- tive ”projection transform” method for creating 2D plots of the adaptive polar transformed image inboth the radius and the angle directions [11]. This subsection explains their work.

Given a adaptive polar transformed image IP with nr samples in which each sample has length nθ_i for i = 1, ..., nr. Let ˆnθ be the number of samples in the angular direction at Ri= Rmax, < be the radius projection and Θ be the projection in the angular direction. The mathematical expressions for calculating the radius projection and the angular projection can be seen in equations 5 and 6,

<(i) =

n_θi

X

j=1

IP (i, j), (6)

Θ(j) =

n_r

X η¹IP

i, ceil j − 1

+ η²IP

i, ceil j

. (7)

(14)

2.2.3 Calculation of distance coefficient between projections

With the radius and angular projections for each of two given points calculated.

The final step in finding if two positions are a match or not could be performed.

Namely, the distance coefficient could now be calculated. The paper ”Image Registration Using Adaptive Polar Transform” by Rittavee M, Yuan F. Z and Robert L. E proposes two different methods for finding the discrepancy ˆdfor two given projections [11]. One uses the correlation function C( ˆΘ^l, ˆΘ^r)to find dˆas can be seen in equation 8:

d = argmax(C( ˆˆ Θ^l, ˆΘ^r)), (8) in which argmax is a function that finds the maximum value of a given set of data and ˆΘ^land ˆΘ^rare angular projections of two matching points from the left and right image. The other method uses the scale parameter azbetween two projections and calculates the discrepancy ˆdas seen in equation 9:

d = log(aˆ z), (9)

in which log calculates the natural logarithm of a given number. The scale parameter azis calculated based on a property of a Fourier transform, namely given a the Fourier transform of a signal x(t) that can be expressed as X(ω) the scaled version of a Fourier transform can be expressed as seen in equations 10 and 11,

x(t) ↔ X(ω), (10)

x(a_zt) ↔ 1

|az|X(ω az

). (11)

With the property expressed in equations 10 and 11, the scale parameter between the two given projections can be found as seen in equation 12 and 13:

Θˆ^r(i) = ˆΘ^l(azi) (12)

|az| = F T^l(0)

F T^r(0) (13)

The distance coefficient between the angular projections is defined as seen in equation 14:

εz= v u u t

ˆ n_θ

X

i=1

[ˆθ^r− ˆθ^l(i − ˆd)]², (14)

l r

(15)

2.3 Ensuring Color Consistency across Multiple Cameras

The paper ”Ensuring Color Consistency across Multiple Cameras” by Adrian Ilie and Greg Welch [9] uses a color calibration method that is done in two steps. First a calibration is done by adjusting hardware register settings in a closed loop until the calibration converges or is close enough to the target values. In the second step, a software based refinement is performed on the images. Here the paper explores three different post-processing methods: ”Linear Least Squares Matching” from E. Weisstein [12], ”RGB to RGB transform” by N.

Joshi [13] and ”General Polynomial Transform” a method proposed by Adrian Ilie and Greg Welch that aims to take in account for the shortcomings of RGB to RGB transform [9]. These three methods are then compared against each other. The following sub sections explains these three post-processing color calibration methods.

2.3.1 Linear Least Squares Matching

In linear least squares matching two coefficients ac and bc are computed in order to provide the best linear transform that can map the color values from the image to the target values for the color channel c{R, G, B}. The transforms are computed by minimizing equation 15:

N S

X

s=1

( ~ICs− (acT~cs0bc))², c{R, G, B}. (15)

In equation 5 ~I_C_s is the component for color channel c of image sample ~I_s and ~Tc_sis the component for color channel c of target value image sample ~Ts.

2.3.2 RGB to RGB transform

In this method a 3x3 matrix that best transforms the samples from the camera images to the target values is found using equation 16:





 I~1

I~₂ ... I~_n







nx3





trrtrgtrb

t_grt_ggt_gb t_brt_bgt_bb





3x3

≈





 T~1

T~₂ ... T~_n







nx3

, (16)

in which each ~I_nis defined as ~I_n =I_r I_g I_b and represents the RGB values

(16)

dress the shortcomings in RGB to RGB transform [9]. This is achieved by gen- eralizing RGB to RGB transforms by adding higher degree terms to be able to account for non linearity’s and adding a bias to allow for translations. Equation 17 shows the general formula for color c{r, g, b} of sample s:

D

X

k=1

(trc_kI_rk

s + tgc_kI_gk

s + tbc_kI_bk

s) + Tc0 ≈ tc_s. (17) In which D is the degree of the polynomial. Ir_s^k, I_gk

s and Ib^k_s are the Red, Green and Blue components of image sample s raised to the power k. The variable Tcs is the value of color channel c{r, g, b} of the target image sample s. txck is the polynomial coefficient of the k^thorder that indicates how much x{r, g, b}contributes to color chanel c and tc0 is a additive term that is used for translation to the color channel c. Modifying equation 6 for n samples with degree D = 2 produces equation 18:





 I_r₁I_r²

1I_g₁I_g²

1I_b₁I_b²

11 Ir₂I_r²₂Ig₂I_g²₂Ib₂I_b²

21 . . . .

Ir_nI_r²_nIg_nI_g²_nIb_nI_b²

n1







nx7





 trc1

trcn

tgc1

tgcn

tbcn

tr0







7x1

≈





 Tc1

Tc2

... Tcn







nx3

. (18)

(17)

3 Method

Figure 1: Flow chart of the method of this work

The methodology of this thesis will follow a layered progression. Figure 1 illustrates the structure of this method. The following subsections explain each step in more detail.

3.1 Literature study

In the first stage, a literature study will be conducted. This will consist of reading a broad spectrum of different articles describing multiple camera array color calibration methods and their performance. Having built a broader un- derstanding of the subject from the first stage, one targeted and one target-less color calibration algorithm will be chosen based on feasibility of implementa-

(18)

3.2 Implementation of algorithms

With the literature study complete, it will be time to implement the chosen algorithms based on the description from the papers: ”Ensuring Color Consis- tency across Multiple Cameras” written by Adrian Llie and Greg Welch [9] (see section 2.3), and ”Color Calibration for Multi-camera Imaging Systems” written by Sabri Gurbuz, Masahiro Kawakita, and Hiroshi Ando [8] (can be found section 2.1).

These methods will be implemented using python [17] with the libraries OpenCV [18], numpy [19] and Matplotlib [20].

3.3 Testing and evaluation of algorithms

Once the implementation of the algorithms is complete the algorithms will be tested and evaluated. The algorithms will be evaluated based on two measure- ments. Color range homogeneity and color accuracy to ground truth.

Color range homogeneity:

In the paper ”Ensuring Color Consistency across Multiple Cam- eras” written by Adrian Llie and Greg Welch [9]. Color range homogeneity is measured through mean inter-sample standard deviation. That is how far apart the color samples from the camera images are with respect to each other. Equation 19 shows how this will be calculated for a given color:

interStdDev[i, 0] =r P(µbs− c[i, 0])² 10 − 1

interStdDev[i, 1] =r P(µg_s− c[i, 1])² 10 − 1

interStdDev[i, 2] =r P(µ_r_s− c[i, 2])² 10 − 1 , (19) in which interStdDev is a 24 by 3 matrix that stores the inter-sample standard deviation for each color and color channel, µb_s, µg_s, µr_sare the average values of the blue, green and red channels of each color sample and c is an iterator that iterates over all ten samples of a given target color. This metric will be used for evaluating the algorithms based on inter-sample standard deviation.

Color accuracy to ground truth:

Color accuracy with respect to ground truth will be evaluated with a root median squared (Euclidean distance) between the measured and desired color values in BGR colorspace, this will be calculated using equation 20:

(19)

In which Ir_s, Ig_s and Ib_s are the Red, Green and Blue components of image sample s and Gr_s, Gg_sand Gb_sare the Red Green and Blue components of the selected squares s.

3.4 Further improvements

Further improvements to the existing algorithms will be implemented, tested and evaluated in a manner similar to the above describes the method. The implementation stage of the targeted and the target-less color calibration methods will help in deciding what can be improved and the verification of the original algorithms will help determine which method is worth trying to improve.

(20)

4 Implementation and Test Construction

This chapter explains how the chosen algorithms where implemented beginning with subsection 4.1 that explains the implementation of the targeted color calibration. Subsection 4.2 explains the implementation of the target-less color calibration method. Subsections 4.3 and 4.4 explain the construction of the evaluation with respect to ground truth and the evaluation with respect to color range homogeneity respectively. Finally, things are rounded off with subsection 4.5 that motivates why and explains the attempted improvement of the targeted color calibration method.

4.1 Targeted color calibration

The following subsections describe the implementation of General Polynomial Transform. This method takes in a set of images as it’s input and outputs a set of calibrated images of the same size. The targeted color calibration in this thesis is performed in the following steps.

1. Detect centroids calibration squares for all images

2. Extract the BGR value from each color calibration square in each image 3. Calculate the average value of each color

4. Estimate matrix multiplier

5. Use matrix multiplier to calibrate images These steps are visualized in figure 2.

Figure 2: Flowchart illustrating the targeted color calibration

(21)

4.1.1 Detection of calibration squares

In order to perform a targeted color calibration, the calibration squares must first be extracted from the images that are to be calibrated. Naturally, this led to the first step in the implementation of the targeted method to be a way of finding the calibration squares. The flowchart in figure 3 illustrates the process of detecting these color calibration squares. The source code for detecting the color calibration squares can be seen in appendix A, lines 6-77 in the function

”findSquares”.

Figure 3: Flowchart illustrating detection of color calibrations squares Using the OpenCV library, the simplest way to detect shapes was to use contour detection, but contour detection was not effective unless some image processing was done first. First, most of the image does not contain relevant information for the detection of squares, so the region of interest selector was used to prompt the user for the location of the color calibration chart. With the region of interest selected, the image was converted to grayscale, a bilateral filter was used to reduce noise and lastly, a Gaussian blur was used to smooth out the image. Figure 4 illustrates the effect of this step.

(22)

Figure 4: Color chart before and after converting to gray scale filtering and smoothing

With the image converted to grayscale, noise reduced and smoothed the next step was to threshold the image and finally detect the contours of the image in order to locate the squares. This was performed by first using adaptive thresholding, the image this gave a fairly good result but the resulting lines had gaps in them. To alleviate this the image was then dilated and eroded using the morphologyEx function from OpenCV. Figure 5 shows the result from these two steps.

Figure 5: Thresholded image before dilation (black) and gaps that where filled after dilation (red)

Finally, the image contours where extracted using findContours from OpenCV.

But the set of contours from findContours still contained contours that were not squares as seen in figure 6.

(23)

Figure 6: Detected contours before filtering

In order to make ensure that the contours detected where only squares filtering was needed. A parameter named ”squareSize” was added to the list of function parameters and a maximum and minimum size for a contour to be valid were calculated as shown in equation 21:

maxSize = (squareSize · 1.3)²

minSize = squareSize², (21)

where the constant 1.3 was calculated through trial and error in test runs.

(24)

The set of all detected contours was iterated over, the area of each contour was calculated using ”contourArea” from OpenCV and the contours that did not have an area between ”maxSize” and ”minSize” where discarded. Each contour that was within ”minSize” and ”maxSize” was then fitted to the best fitting square and appended to the set of filtered contours. Figure 7 illustrates the set of detected squares.

Figure 7: The initial detected squares

As seen in figure 7, many squares are still not detected and some of the squares that were detected before the filtering have been discarded as a result of not being close to the valid square size. In order to solve this, several methods where attempted.

First attempt at detecting missing color calibration squares In the first attempt, each detected squares centroids height was calculated. Using this height for a given detected squares centroid a set of candidate centroids was generated by taking the x position of the current centroid and creating 5 new points that where a distance in the x-direction of xP os − i · distance from the current centroid and 5 new points that where xP os + i · distance from the current position. The points that were out of the scope of the image and the points that were too close to an already generated point from a previous generation of points where discarded. This did work to a degree, but for images that had a color chart that was not flat towards the camera the generated points where not correct.

Second attempt at detecting missing color calibration squares To address this shortcoming a slightly refined attempt of a similar method was attempted.

Instead of having each centroid of a given detected square generate points with the same height, detected squares that belonged to the same row where grouped together and a line for each row was generated using linear regression. Similarly, the distance in between each generated centroid was the current centroid being generated plus or minus the given position with the small difference that some trigonometry was used since the distance between points was now the hypotenuse of an equilateral triangle. Finally, the generated centroids that were outside of the scope of the image where discarded. This new

(25)

where discarded, images where there was significant extra space between the end of the color chart and the end of the image ended up having excess centroids.

Final attempt at detecting missing color calibration squares and selection of fallback method Finally, one last attempt was made to detect the centroids of the missing squares. This time, lines were estimated using linear regression for the rows and the columns of the color chart. The intersections of the lines were then used to estimate the position of the centroids for the color calibration chart. In theory, this seemed like a good idea, but unfortunately, the method did not work as intended. For most color charts, the number of detected squares was sufficient but not always. When there where no detected squares on any given row or column. There would still be squares that were undetected.

It was at this point time started to become an issue and in order to complete the project in a feasible amount of time. The detection of missing squares from the color calibration chart was settled to manual selection of missing squares.

The set of missing centroids of color chart squares was then emplaced into the array of automatically detected color calibration squares. The final algorithm for detecting the color calibrations squares can be seen in figure 3.

(26)

4.1.2 Estimation matrix multiplier

Now that the colors from the color calibration charts have been extracted for all ten images. The next step in the implementation of ”General Polynomial Transform” was to estimate the matrix multiplier.

Firstly, in order to avoid selecting an outlier to estimate the matrix multiplier. An average BGR value for each color was calculated. This was achieved by the algorithm shown in algorithm 1. The source code for algorithm 3 can be seen in appendix A, lines 230-244 in the function ”averageOutColors”.

Algorithm 3Algorithm for finding average BGR value of each color

1: averages ← zeros(24, 3)

2: for i ← 0 : 9 do

3: for j ← 0 : 23 do

4: if i == 0 then

5: averages[j, 0] ← calibrationSquares[i, j][0]

8: else

9: averages[j, 0] ←

q(averages[j,0]²+calibrationSquares[i,j][0]²) 2

12: end if

13: end for

14: end for

Next, now that the average value of each color calibration square has been calculated. The matrix multiplier was calculated using the average of each of the 24 colors to create the 24 by 7 matrix used for estimating the matrix multiplier. Algorithm 4 shows how the 24 by 7 matrix was calculated.

Algorithm 4Calculation of matrix of observed values

1: matM ul ← zeros(24, 7)

2: i ← 0

3: fora in averages do

4: matM ul[i, :] ←a[0] a[0]² a[1] a[1]² a[2] a[2]² 1

5: i ← i + 1

6: end for

Finally, the calibration matrix multiplier was calculated as can be seen in algorithm 5.

(27)

Algorithm 5Estimation of calibration matrix multiplier

1: calibM atM ul ← zeros(24, 3)

2: for i ← 0 : 2 do

3: calibM atM ul[:, i] ← P inv(matM ul) · groundT ruth[:, i]

4: end for

Where in algorithm 5. ”groundTruth” is a 24 by 3 matrix containing the ground truth values from the color calibration chart and Pinv is a function that calculates the pseudo-inverse of a given matrix. The source code for algorithms 2 and 3 can be seen in appendix A, lines 182-193 in the function ”calculateMu- tipliers”.

(28)

4.1.3 Calibration of images

With the calibration matrix multiplier calculated, the images could now be calibrated. To achieve this, all ten images where iterated over and every pixel in each image was accessed, the values of the blue, green and red channels of the current pixel where then used to calculate the current multiplier as shown in equation 22:

currM ultiplier =c[0] c[0]² c[1] c[1]² c[2] c[2]² 1 , (22) in which c i a 3 by 1 matrix containing the Blue, Red and Green channels of the current pixel.

The new BGR value for the current pixel was then calculated as shown in algorithm 4. The soruce code that is described in algorithm 6 can be seen in appendix A, lines 196-228 in the function ”calibrateImage”.

Algorithm 6Calculation of new BGR value

1: blue ← currM ultiplier · calibM atM ul[:, 0]

2: green ← currM ultiplier · calibM atM ul[:, 1]

3: red ← currM ultiplier · calibM atM ul[:, 2]

4: if blue > 255 then

5: currImage[i, j, 0] ← 255

6: else if blue < 0 then

8: else

9: currImage[i, j, 0] ← blue

10: end if

11: if green > 255 then

12: currImage[i, j, 1] ← 255

13: else if green < 0 then

15: else

16: currImage[i, j, 1] ← green

17: end if

18: if red > 255 then

19: currImage[i, j, 2] ← 255

20: else if red < 0 then

22: else

23: currImage[i, j, 2] ← red

24: end if

(29)

4.2 Target-Less Color Calibration

The following subsections describe the implementation of the target-less color calibration. This method takes in a set of images as it’s input and outputs a set of calibrated images of the same size. The target-less color calibration follows the following steps.

1. Detect the initial set of feature points between two given images using SIFT

2. Find the initial set of matching pairs between two given images using the Fast Library for Approximate Nearest Neighbors (FLANN) based matcher, a Lowe’s ratio test and RANSAC

3. Verify matching points using adaptive polar transform verification 4. Extract BGR values from the matching points and discard matching points

that have a ∆E > 2σ

5. Calculate the matrix multiplier 6. Calibrate the images

These steps are visualized in the flowchart that can be seen in figure 8.

Figure 8: Flowchart illustrating the complete target-less color calibration method

(30)

4.2.1 Extracting of feature points using SIFT

Given two images, the initial set of feature points of each image was extracted using SIFT. This was achieved by first converting both images to gray scale and creating a SIFT detector using the OpenCV library [14]. Figure 9 shows an example the initial set of feature points for two images.

Figure 9: The initial set of feature points for two images

4.2.2 Matching feature points

The initial set of correspondences between two given images was found by first using a FLANN based matcher from the OpenCV library and using the

”knnMatch” function [1] this was achieved by following the documentation for finding a set of matching feature points [15]. This gave a set of matching feature points. The matches where then iterated through and the matches that passed a Lowe’s ratio test where saved [4]. If the number of correspondences was larger then a given threshold, then the matching feature points where filtered even further by using RANSAC from OpenCV. If the amount of matching feature points was smaller then the threshold, then set of matching feature points were not filtered further using RANSAC. The code for finding the initial set of matching feature points between two given images can be found in appendix B on lines 2-69. Figure 10 shows an example of what the initial correspondences between the two example images look like.

(31)

Figure 10: The initial set of matching points filtered by FLANN based matcher, Lowe’s ratio test and RANSAC

(32)

4.2.3 Adaptive polar transform verification

With the initial set of correspondences found, the next step was to verify these matches with adaptive polar transform verification. Subsection 2.2 explains the related work that was followed when implementing adaptive polar transform verification. The following text explains how adaptive polar transform verification was implemented.

The first step in the implementation of adaptive polar transform verification was to reconstruct the algorithm that transforms images from cartesian coordinates to polar coordinates by sampling a given point adaptively. This was done as the following describes. First, the image that is being sampled is converted to gray scale. A new image is then initialized with all points set to zero with a height equal to the desired sampling radius and a width of 8 · radius. With the image converted to grayscale and the new image having been initialized to store the new image, the transform can now be performed. Algorithm 7 shows the pseudocode for this transform. The source code for algorithm 5 can be seen in appendix B, lines 106-133 in the function ”adaptivePolarTransform”.

Algorithm 7Algorithm for performing adaptive polar transform

1: for i ← 1 : f loor(radius) do

2: nθ_i ← 8 · i

3: for j ← 1 : nθ_ido

4: trasf ormedImage[i, j] ← sampleImage[yP os + i·cos(2·π·j)

n_θi , xP os +

i·sin(2·π·j) n_θi ]

5: end for

6: end for

Figure 11 shows an example of an image before and after an adaptive polar transform. The center point of the image was used as the starting point for sampling the image and the radius was the height of the image divided by two.

(33)

Figure 11: Example image before adaptive polar transform (above) and after (below)

Next, in order to be able to determine if two points are a match with scale and rotation invariance. The adaptively polar-transformed images are projection transformed with respect to scale (radius) and rotation. From these two projections, the scale and rotation parameters can be calculated. The pseudocode in algorithm 8 shows how the radial projection is calculated and the pseudocode in algorithm 9 shows how the rotation projection was calculated.

The source code for algorithm 8 and 9 can be found in appendix B, lines 135-189 in the function ”projectionTransform”.

Algorithm 8Calculation of radius projection

1: nˆ_θ← r_max· 8

2: radiusP rojection ← []

3: for i ← 1 : radius do

4: nθ_i ← i · 8

5: Ωi← _nⁿ^ˆ^θ

θi

6: radiusSun ← 0

7: for j ← 1 : nθ_ido

8: radiusSum ← radiusSum + polartransf ormedImage[i, j]

9: end for

10: radiusP rojection.append(Ωi· radiusSum)

11: end for

(34)

Algorithm 9Calculation of rotation projection

1: rotP roj ← []

2: for j ← 1 : ˆn_θdo

3: sum ← 0

4: for i ← 1 : radius do

5: nθ_i← i · 8

6: Ωi← _nⁿ^ˆ^θ

θi

7: n¹_ij ← Ωi· (j − 1) − f loor(Ωi· (j − 1))

8: n²_ij ← 1 − n¹_ij

9: lhs_jpos← ceil(^j−1_Ω

i )

10: rhsjpos← ceil(_Ω^j

i)

11: if lhsjpos== 0 ∧ rhsjpos== 0then

12: lef tOperand ← 0

13: rightOperand ← 0

14: else if rhsjpos == 0then

15: rightOperand ← 0

16: lef tOperand ← polartransf ormedImage[i, lhsjpos]

17: else if lhsjpos== 0then

18: lhsjpos← 0

19: rhs_jpos← polartransf ormedImage[i, rhsjpos

20: else

21: lef tOperand ← polartransf ormedImage[i, lhs_jpos]

22: rhs_jpos← polartransf ormedImage[i, rhs_jpos

23: end if

24: sum ← sum + n¹_ij· lef tOperand + n²_ij· rightOperand

25: end for

26: rotP roj.append(sum)

27: end for

Figure 12 shows the example point being projected, figure 13 shows the radius projection of the example point from figure 12 and figure 14 shows the rotation projection of the example point from figure 12.

(35)

Figure 12: The example point from both images that was projected

Figure 13: Radius projection of the example point, the left side contains the projection of the example point from the left image and the right side contains the projection of the example point from the right image

(36)

With the radius and rotation projections complete, the next step in implementing adaptive polar transform verification was to calculate the distance between the projections of matching feature points. In order to be able to calculate the distance between two given points, the scale parameter between their rotation projections had to be calculated. Projections of unequal length where discarded, for the remaining projections the distance coefficient was calculated.

Equation 23 shows how the scale parameter was calculated, θ^l(i) = θ^r(a · i)

a = F T^l(0)

F T^r(0). (23)

In which a is the scale parameter between the two rotation projections, F T^lis the Fourier transform of the left projection and F T^ris the Fourier transform of the right projection. Using the scale parameter the displacement between the two projections can then be calculated as shown in equation 24:

d = log(a),ˆ (24)

where ˆdis the disparity between the two projections.

(37)

Algorithm 10 shows the full algorithm for calculating the distance coefficient between projections. The source code for algorithm 8 can be seen in appendix B, lines 200-223 in the function ”calcDistCoefficient”.

Algorithm 10Algorithm for calculating the distance coefficient between projections

1: if len(θl)! = len(θr)then

2: ./*− inf is a reserved value to indicate that the match is invalid and will later be discarded*/

3: distance ← − inf

4: else

5: a ← _{F T}^{F T}_r^l⁽⁰⁾₍₀₎

6: if a! = N aN then

7: θ_r← θ_r· a

8: end if

9: if a == 0 then

10: d ← 0ˆ

11: else

12: d ← log(a)ˆ

13: end if

14: distance ← 0

15: for i ← 0 : len(θl)do

16: if ˆd == N aN then

17: distance ← distance + (θ_r[i] − θ_l[i])²

18: else

19: if round(i − ˆd) < len(θl) ∧ round(i − ˆd) >= 0then

20: distance ← distance + (θr[i] − θl[round(i − ˆd)])²

21: end if

22: end if

23: end for

24: distance ←√

distance

25: end if

With the distance for each matching pair calculated, the matches that have too large of a distance could now be filtered out. Algorithm 11 illustrates how this was achieved. The source code for algorithm 9 can be seen in appendix B, lines 88-104 in the function ”adaptivePolarTransformVerification”.

(38)

Algorithm 11Discarding of invalid matches based on distance between projections

1: for m : matches do

2: µθ_l← average(m.θl)

3: µθ_r ← average(m.θr)

4: if distance! = − inf then

5: stdDev_r← 0

6: stdDev_l← 0

7: for i ← 0 : len(θl)do

8: stdDev_r= stdDev_r+ (µ_θ_r− m.θ_r[i])²

9: stdDev_l= stdDev_l+ (µ_θ_l− m.θ_l[i])²

10: end for

11: stdDevr←q

stdDev_r len(m.θ_r)−1

12: stdDevl←q

stdDev_l len(m.θ_l)−1

13: if distance > |µl+ stdDevl− (µr− stdDevr)|then

14: matches.remove(m)

15: end if

16: else

18: end if

19: end for

Figure 15 shows the set of matching points after the adaptive polar transform verification and the set of points that where removed by adaptive polar transform verification.

Figure 15: The set of matching points between two example images after adaptive polar transform verification (in blue) and the set of points that where removed (in red)

(39)

4.2.4 Extraction of BGR values from matching points and delta E filtering After the set of reliable correspondences between two given images was found.

The average BGR value of each point was calculated by sampling points of each point adaptively (in a similar manner as in adaptive polar transform) and averaging out the value of each sample together into one BGR value.

Next, matching points that posed a ∆E (defined in equation 2 in chapter 2.1) larger then 2σ ± avg (defined in chapter 2.1) where discarded as seen in algorithm 12. The source code for algorithm 10 can be seen in appendix B, lines 257-276 in the function ”calcDeltaE”.

Algorithm 12Filtering out of matching points that pose a larger ∆E then 2σ

1: avg ← 0

2: for m : matchess do

3: m.∆E ← p(m.BGRl[0] − m.BGRr[0])²+ (m.BGRl[1] − m.BGRr[1])² +(m.BGRl[2] − m.BGRr[2])²

4: avg ← avg + m.∆E

5: end for

6: avg ← len(matches)^avg 7: stdDev ← 0

8: for m : matches do

9: stdDev ← stdDev + (m.∆E − avg)

10: end for

11: stdDev ←q _stdDev

len(stdDev)−1 12: for m : matches do

13: if m.∆E > average + 2 · stdDev ∨ m.∆E < average − 2 · stdDev then

15: end if

16: end for

(40)

4.2.5 Calculation of matrix multiplier

Now that the set of matching feature points has been distilled to reliable matches from the above filtering methods. The matrix multiplier that best transforms all the BGR values from the image to be calibrated to the BGR values from the model image was estimated using linear regression. Algorithm 13 illustrates this process. The source code for algorithm 11 can be seen in appendix B, lines 278-295 and lines 297-313 in the functions ”calculateMultiplier” and

”calculateAlphaAndBeta”.

Algorithm 13Calculation of matrix multiplier for calibration of image

1: sumRef ← zeros(1, 3)

2: sumT oCal ← zeros(1, 3)

3: for i ← 0 : 2 do

4: sumRef ←P ref erenceV alues[i][:]

5: sumtoCal ←P toCalibrateV alues[i][:]

6: end for

7: squareSumRef ← zeros(1, 3)

8: toCalT imesRef Sum ← zeros(1, 3)

9: for i ← 0 : 2 do

10: for j ← 0 : len(ref erenceV alues) do

11: squareSumRef [i] ← squareSumRef [i] + ref erencevalues[i][j]²

12: toCalT imesRef Sum[i] ← toCalT imesRef Sum[i] + ref erencevalues[i][j] · toCalibrateV alues[i][j]

13: end for

14: end for

15: n ← len(ref erenceV alues[i][:])

16: a ← zeros(1, 3)

17: b ← zeros(1, 3)

18: for i ← 0 : 2 do

19: a[i] ← toCalT imesRef Sum[i]−n·(^{sumRef [i]}_n ·sumT oCal[i]

n )

squareSumRef [i]−n·(^sumRef_n )² 20: b[i] ← sumT oCal[i]

n −^{sumRef [i]}_n

21: end for

22: matmul ←





a[0] 0 0 b[0]

0 a[1] 0 b[1]

0 0 a[2] b[2]





(41)

4.2.6 Calibration of images

Once the calibration matrix multiplier was found for two given images. The image to be calibrated could be calibrated with respect to the reference image.

This was performed for all ten images such that image 2 was calibrated with respect to image 1, image 3 with respect to image 2 and so on and so forth.

Figure 16 illustrates the calibration process for 10 images.

Figure 16: Color calibration process for ten images

(42)

4.3 Testing performance of calibration with respect to ground truth

In order to evaluate the accuracy towards target values, the BGR value of the calibration squares was extracted using the same method as in section 4.1.1.

The Euclidean distance for each image with respect to ground truth for each color in each image was calculated as seen in algorithm 12. The code for algorithm 14 can be seen in appendix C, lines 1-6 in the function ”calculateRMS”.

Algorithm 14Calculation of euclidean distance to ground truth

1: distances ← zeros(10, 24)

2: for i ← 0 : 9 do

3: for j ← 0 : 23 do

4: distances[i, j] ←p(Ibs[i, j] − G_b_s[j])²+ (I_g_s[i, j] − G_g_s[j])²+ (I_r_s[i, j]−

Gr_s[j])²

5: end for

6: end for

In which Ir_s, Ig_sand Ib_s are the Red, Green and Blue components of image sample s and Gr_s, Gg_s and Gb_s are the Red Green and Blue components of the ground truth sample s.

The same test was also performed for the same images before calibration in order to be able to compare color accuracy before calibration versus color accuracy after calibration.

(43)

4.4 Testing performance of calibration with respect to color range homogeneity

As before, the method for extracting color calibration squares BGR values from section 4.1.1 was used in order to extract the BGR values from each color calibration square. The inter-sample standard deviation was then calculated for each color channel and each color calibration square. Algorithm 15 shows how the inter-sample standard deviation for each color calibration square across all ten images. The source code for algorithm 15 can be seen in appendix C, lines 7-20 in the function ”calculateInterSampleStdDev”.

Algorithm 15Calculation of inter sample standard deviation

1: averages ← zeros(24, 3)

2: for c : colors do

3: for i ← 0 : 23 do

4: averages[i, :] ← averages[i, :] + c[i].BGRV al

5: end for

6: end for

7: for i ← 0 : 23 do

8: averages[i, :] ← [averages[i,0]

len(colors) ,averages[i,1]

len(colors) ,averages[i,2]

len(colors) ]

9: end for

10: stdDevs ← zeros(24, 3)

11: for c : colors do

12: for i ← 0 : 23 do

13: stdDevs[i, :] ← stdDevs[i, :] + [(c[i].BGRV al[0] − averages[i, 0])², (c[i].BGRV al[1] − averages[i, 1])², (c[i].BGRV al[2] − averages[i, 2])²]

14: end for

15: end for

16: for i ← 0 : 23 do

17: stdDevs[i, :] ← [q stdDevs[i,0]

len(colors)−1,q stdDevs[i,1]

len(colors)−1,q stdDevs[i,2]

len(colors)−1]

18: end for

(44)

4.5 Modification of General Polynomial Transform

The extracted colors used for calculating the calibration matrix multiplier uses an average of the color charts from all ten images in order to avoid selecting an outlier as it’s reference. This seemed to have an inherent flaw in avoiding outliers, namley even those colors that where outliers where a part of the average. To avoid this the Euclidean distance to ground truth for every color sample for every color sample. For each color from all ten images, an average and standard deviation of their Euclidean distances to ground truth was calculated. Those color samples that where within one standard deviation from the mean where added to the average of a given color. Algorithm 16 visualizes this process.

Target-less and targeted multi-camera color calibration

Independent degree project - first cycle

Abstract

Foreword

Table of Contents

Terminology

1 Introduction

1.1 Background and Problem motivation

1.2 Overall purpose

1.3 Limitations

1.4 Concrete and verifiable goals

1.5 Overview

2 Related work

2.1 Color Calibration for Multi-camera Imaging Systems

2.2 Adaptive Polar Transform Verification

2.3 Ensuring Color Consistency across Multiple Cameras

3 Method

3.1 Literature study

3.2 Implementation of algorithms

3.3 Testing and evaluation of algorithms

3.4 Further improvements

4 Implementation and Test Construction

4.1 Targeted color calibration

4.2 Target-Less Color Calibration

4.3 Testing performance of calibration with respect to ground truth

4.4 Testing performance of calibration with respect to color range homogeneity

4.5 Modification of General Polynomial Transform