Ever since watching that video I have wondered how one would go about to judge the winner of such an alleged [World Freehand Circle Drawing Championship](https://www.youtube.com/watch?v=u1J5ANnq0T8) (WFHCDC). While researching for this post I finally figured it out. On his web page Alexander in the story behind the video ["reveals"](http://slamdunkmath.blogspot.com): *They have a laser machine called the **circleometer** that creates the perfect circle closest to the one you drew. The circleometer then calculates the difference in area between the laser circle and the circle that you drew. The machine then calibrates the area difference as if you had drawn a circle with radius one meter. The person with the smallest area difference is declared the world freehand circle drawing champion.* Aha! Imaginary circleometers are expensive and my dean of study most likely isn't open for investing in perfect circle measurement equipment... So here is a cheaper solution involving a mobile device camera, [R](https://www.r-project.org) and the [`imager`](https://cran.r-project.org/web/packages/imager/index.html) package by [Simon Barthelmé](https://sites.google.com/site/simonbarthelme/) et al. Altogether a combination of modern **data science** tools, which my dean of study is most likely to approve! We'll use a screenshot from the perfect circle video as motivating example to guide through the 3 phases of the method: 1. Image rectification 2. Freehand circle identification and perfect circle estimation 3. Quantifying deviation from the perfect circle We start by loading the screenshot into R using `imager`: ```{r LOADIMAGE, echo=TRUE, message=FALSE, warning=FALSE} library("imager") file <- "circle2.png" img <- imager::load.image(file.path(fullFigPath, file)) ``` ```{r PLOTSCREENSHOT} plot(img) ``` ## Image rectification The image clearly suffers from perspective distortions caused by the camera being positioned to the right of the circle and, hence, not being orthogonal to the blackboard plane. Furthermore, small lense distortions are also visible - for example the right vertical line of the blackboard arcs slightly. Since the video contains no details about what sort of lense equipment was used, for the sake of simplicity, we will ignore lense distortions in this post. If such information is available one can use a program such as [RawTherapee](http://rawtherapee.com) (available under a GNU GPL v3 license) to read the EXIF information in the meta data of the image and automatically correct for lens distortion. To rectify the image we estimate the parameters of the 2D projection based on 4 ground control points (GPC). We use R's `locator` function to determine the pixel location of the four corner points of the blackboard in the image, but could just as well use any image analysis program such as [Gimp](https://www.gimp.org). Furthermore, we need the true object coordinates of these GPC. Unfortunately, these are only approximately available to due lack of knowledge of the size of the blackboard in the classroom. As a consequence a *guesstimate* of the horizontal length is used. ```{r AQUIREGPC, echo=TRUE, eval=FALSE} plot(img) p <- locator(4) p <- round(cbind(p$x, p$y)) dump(list=c("p"), "") ``` ```{r DEFINEGPC} ## Corner coordinates of blackboard in image (found with Gimp or locator()) p <- rbind(c(0,318),c(2016,100),c(1988,1242),c(0,980)) ## True coordinates - since we don't know the exact measurements of the ## blackboard we have to guess a little. Move coordinates to the right, s.t. all ## object coordinates are positive (otherwise they will not appear in the image. ## Offset to the right here: 1000 pixels dx_blackboard <- 1200 pp <- rbind(c(-dx_blackboard ,100),c(2016,100),c(1988,1242),c(-dx_blackboard,1242)) + matrix(c(dx_blackboard, 0),4,2,byrow=TRUE) ``` These points are now used to rectify the image by a Direct Linear Transformation (DLT) based on exactly 4 control points [@hartley_zisserman2004, Chapter 4]^[Alternatively, see slide 18 and onward in https://ags.cs.uni-kl.de/fileadmin/inf_ags/3dcv-ws11-12/3DCV_WS11-12_lec04.pdf]. That is the parameters of the 3x3 transformation matrix $H$ in homogeneous coordinates are estimated such that $p' = H p$, see the [code](`r paste0("https://raw.githubusercontent.com/hoehleatsu/hoehleatsu.github.io/master/_source/",current_input())`) on github for details. ```{r DLT, results='hide'} ## Direct Linear Transformation (DLT) algorithm using inhomogeneous ## solution, i.e. H[3,3] = 1 and w_i=1. See Sect 4.1.2 of ## Hartley and Zisserman (2004) ## ## (H[1,1]*x + H[1,2]*y + H[1,3]) - x' * (H[3,1]*x + H[3,2]*y + 1) = 0 ## (H[2,1]*x + H[2,2]*y + H[2,3]) - y' * (H[3,1]*x + H[3,2]*y + 1) = 0 ## ## [[ x y 1 0 0 0 -x'*x -x'*y = x'] , ## [ 0 0 0 x y 1 -y'*x -y'*y = y']] A <- list() b <- list() ##Build RHS and LHS matrix components for each point pair for (i in 1:4) { A[[paste0(i)]] <- matrix(c( p[i,1],p[i,2],1, 0, 0, 0, -pp[i,1]*p[i,1], -pp[i,1]*p[i,2], 0, 0, 0, p[i,1], p[i,2], 1, -pp[i,2]*p[i,1], -pp[i,2]*p[i,2]), 2,8, byrow=TRUE) b[[paste0(i)]] <- matrix(c(pp[i,1], pp[i,2]),2,1) } ##Glue matrices together A_matrix <- do.call(rbind, A) b_matrix <- do.call(rbind, b) ##Solve equation system h_vec <- solve(A_matrix, b_matrix) ##Form H matrix from the solution H <- matrix(c(h_vec,1), 3,3, byrow=TRUE) H ``` We can implement the rectifying transformation using the `imager::warp` function: ```{r RECTIFY, echo=TRUE} ##Transform image coordinates (x',y') to (x,y), i.e. note we specify ##the back transformation p = H * p', so H here is the inverse. map.persp.inv <- function(x,y, H) { out_image <- H %*% rbind(x,y,1) list(x=out_image[1,]/out_image[3,], y=out_image[2,]/out_image[3,]) } ##Pad dx_blackboard pixels to the right to make space for the blackboard coming closer img_padded <- pad(img, nPix=dx_blackboard, axes="x", pos=1) ##Warp image warp <- imwarp(img_padded, map=function(x,y) map.persp.inv(x,y,solve(H)),coordinates="absolute", direction="backward") ``` The result looks as follows: ```{r PLOTRECTIFICATION, fig.height=8} layout(c(1,2)) plot(img, main="Original") lines(c(p[,1],p[1,1]), c(p[,2],p[1,2]), col="magenta",lwd=2) points(p[,1], p[,2], cex=3, pch=20, col="steelblue") plot(warp, main="Rectified") lines(c(pp[,1],pp[1,1]), c(pp[,2],pp[1,2]), col="magenta",lwd=2) points(pp[,1], pp[,2], cex=3, pch=20, col="steelblue") ``` Please notice the different x-axes of the two images when comparing them. For faster computation and better visualization in the remainder of this post, we crop the x-axis of the image to the relevant parts of the circle. ```{r, echo=TRUE} warp <- imsub(warp, x %inr% c(dx_blackboard, nrow(warp))) ``` ```{r PLOTRECTIFIED} plot(warp, main="Rectified and Cropped") ``` ## Freehand circle identification As described in the `imager` [edge detection tutorial](https://dahtah.github.io/imager/canny.html) we use length of the gradient to determine the edges in the image. This can be done by applying filters to the image. ```{r, echo=TRUE} ##Edge detection function. Sigma is the size of the blur window. detect.edges <- function(im, sigma=1) { isoblur(im,sigma) %>% imgradient("xy") %>% enorm() %>% imsplit("c") %>% add } #Edge detection filter sequence. edges <- detect.edges(warp,1) %>% sqrt ``` To detect the circle from this we specify a few seed points for a [watershed](https://en.wikipedia.org/wiki/Watershed_%28image_processing%29) algorithm with a priority map inverse proportional to gradient magnitude. This includes a few points inside and outside the circle and a few points *on* the circle. Note: a perfect circle would have no border, but when drawing a circle with a piece of chalk it's destined to have a thin border line. ```{r} #Priority map with priority inverse proportional to gradient magnitude pmap <- 1/(1+edges) ``` ```{r, eval=FALSE} layout(1) plot(pmap) background <- locator(2) background <- round(cbind(background$x, background$y)) dump("background","") foreground <- locator(2) foreground <- round(cbind(foreground$x, foreground$y)) dump("foreground","") ``` ```{r SEGMENTATION, results='hide'} ##Seed image seeds <- imfill(dim=dim(pmap)) background <- structure(c(1400, 956, 239, 507), .Dim = c(2L, 2L)) foreground <- structure(c(508, 1078, 1018, 234), .Dim = c(2L, 2L)) seeds[cbind(background,1,1)] <- 1 seeds[cbind(foreground,1,1)] <- 2 ##Number of seed points sum(seeds >0) ##Run watershed algorithm wt <- watershed(seeds,pmap) mask <- add.colour(wt) #We copy along the three colour channels layout(t(1:2)) plot(warp * (mask==1),main="Background") points(background, col="steelblue", pch=20) plot(warp * (mask==2),main="Foreground") points(foreground, col="magenta", pch=20) ``` We can now extract the circle by: ```{r, echo=TRUE} ##Just the circle freehandCircle <- (warp * (mask==2) > 0) %>% grayscale ##Total area covered by the circle freehandDisc <- label(freehandCircle, high_connectivity=TRUE) > 0 ``` ...and by morphological operations we get just the outer border as a hairline ```{r THINBORDER, echo=TRUE} dilatedDisc <- freehandDisc %>% dilate_rect(sx=3,sy=3) freehandCircleThinBorder <- (freehandDisc - dilatedDisc) != 0 ``` ```{r PLOTTHINBORDER, fig.width=8, fig.height=4} par(mar=c(2,2,2,2)) plot(1-freehandCircleThinBorder, main="Extracted freehand circle") ``` ## Perfect circle estimation Once the freehand circle path in the image has been identified, we need to find the best fitting *perfect* circle matching this path. This problem is elegantly solved by @coope1993, who formulates the problem as finding center and radius of the circle minimizing the squared Euclidean distance to $m$ data points $a_j$, $j=1,\ldots,m$. Denoting by $c$ the center of the circle and by $r>0$ the radius we want to find the solution of $$ \min_{c\in \mathbb{R^2}, r>0} \sum_{j=1}^m F_j(c,r)^2, \quad\text{where}\quad F_j(c,r) = \left|r - ||c-a_j||_2\right|, $$ and $||x||_2$ denotes Euclidean distance. Because the curve fitting minimizes the distance between an observed point $a_j$ and its closest point on the circle and thus involves both the $x$ and the $y$ direction , this is a so called **[total least squares](https://en.wikipedia.org/wiki/Total_least_squares)** problem. The problem is non-linear and can only be solved by iterative numerical methods. However, the dimension of the parameter space can be reduced by one, because given the center $c$ we can determine that $r(c)=\frac{1}{m} \sum_{j=1}^m ||c-a_j||_2$. ```{r FITCIRCLE, echo=TRUE} ##Compute radius given center radius_given_center <- function(center, dist=NULL) { if (is.null(dist)) { a <- as.matrix(where(freehandCircleThinBorder > 0)) dist <- sqrt((a[,1] - center[1])^2 + (a[,2] - center[2])^2) } return(mean(dist)) } ##Target functin of the total least squares criterion of Coope (1993) target_tls <- function(theta) { ##Extract parameters center <- exp(theta[1:2]) ##Total least squares criterion from Coope (1993) a <- as.matrix(where(freehandCircleThinBorder > 0)) dist <- sqrt((a[,1] - center[1])^2 + (a[,2] - center[2])^2) ##Compute radius given center radius <- radius_given_center(center, dist) F <- abs( radius - dist) sum(F^2) } res_tls <- optim(par=log(c(x=background[1,1], y=background[1,2])), fn=target_tls) center <- exp(res_tls$par) fit_tls <- c(center,radius=radius_given_center(center)) fit_tls ``` We illustrate the freehand circle (in black) and the fitted circle (magenta) on top of each other using the alpha channel. You have to study the image carefully to detect differences between the two curves! ```{r DRAWCIRCLE, fig.width=8, fig.height=4} par(mar=c(2,2,0,2)) nPoints <- 1000 t <- seq(0,2*pi, length=nPoints) circle <- cbind(fit_tls[3]*cos(t) + fit_tls[1], fit_tls[3]*sin(t)+fit_tls[2]) plot(1-freehandCircleThinBorder) lines(circle[,1], circle[,2], col=rgb(0.8,0.4,0.9,0.3),lwd=3) points(fit_tls[1], fit_tls[2], pch=20, cex=1, col=rgb(0.8,0.4,0.9,0.4)) ``` ## Quantifying the circularness of the freehand circle We quantify the **circularness** of the freehand circle by contrasting the area covered by it with the area of the fitted perfect circle. The closer this ratio is to 1 the more perfect is the freehand circle. ```{r, echo=TRUE} ##Area of the freehand drawn disc areaFreehandDisc <- sum(freehandDisc) ##Area of the disc corresponding to the idealized circle fitted ##to the freehand circle areaIdealDisc <- pi * fit_tls["radius"]^2 ##Ratio between the two areas ratio_area <- as.numeric(areaFreehandDisc / areaIdealDisc) ratio_area ``` Yup, it's a pretty perfect circle! Since the fitted circle already takes the desired shape into account, my intuition is that this ratio is a pretty good way to quantify circularness. However, to avoid **measurehacks**, we use as backup measure the circleometer approach: for each point on the freehand circle we measure its distance to the freehand circle and integrate/sum this up over the path of the freehand circle. We can approximate this integration using image pixels as follows. ```{r DISTANCE, echo=TRUE} ##Create a pixel based circle in an image of the same size as the ##freehandCircle img. For visibility we use a border of 'border' pixels ##s.t. circle goes [radius - border/2, radius + border/2]. Circle <- function(center, radius, border) { as.cimg(function(x,y) { lhs <- (x-center[1])^2 + (y-center[2])^2 return( (lhs >= (radius-border/2)^2) & (lhs <= (radius+border/2)^2)) }, dim=dim(freehandCircle)) } ##Build pixel circle based on the fitted parameters C_tls <- Circle(fit_tls[1:2], fit_tls[3], border=1) ##Calculate Euclidean distance to circle for each pixel in the image dist <- distance_transform(C_tls, value=1, metric=2) ##Distance between outer border of freehand circle and perfect circle area_difference <- sum(dist[freehandCircleThinBorder>0]) ##Compute area difference and scaled it by the area of the fitted disc ratio_areadifference <- as.numeric(area_difference / areaIdealDisc) ``` The image below illustrates this by overlaying the result on top of the distance map. For better visualization we zoom in on the 270-300 degree part of the circle (i.e. the bottom right). In magenta is the fitted perfect circle, in gray the freehand circle and the area between the two paths is summed up over the entire path of the freehand circle: ```{r PLOTAREADIFF} CircleDisc <- function(center, radius, border) { as.cimg(function(x,y) { lhs <- (x-center[1])^2 + (y-center[2])^2 return( lhs <= (radius+border/2)^2) }, dim=dim(freehandCircle)) } CDisc <- CircleDisc(fit_tls[1:2], fit_tls[3], border=1) xr <- c(fit_tls[1], fit_tls[1] + 0.7*fit_tls[3]) yr <- c(fit_tls[2]+0.7*fit_tls[3], fit_tls[2] + 1.1*fit_tls[3]) distsub <- imsub(dist, x%inr% xr, y %inr% yr) par(mar=c(0,0,0,0)) plot( distsub/max(distsub) + imsub(freehandDisc * (CDisc==0), x %inr% xr, y %inr% yr), axes=FALSE) highlight(imsub(C_tls>0, x %inr% xr, y %inr% yr), col="magenta") highlight(imsub(freehandCircleThinBorder, x %inr% xr, y %inr% yr), col="gray") ``` We obtain `ratio_areadifference`= `r sprintf("%.5f",ratio_areadifference)`. Thus also this measure tells us: it's a pretty perfect circle! To summarise: The output on the display of the judge's Circle-O-Meter App (available under a GPL v3 license) at the World Freehand Circle Drawing Championship would be as follows: `r emo::ji("wink")` ```{r HUD} ##Show our findings Head-up-display (HUD) style on top of image par(mar=c(0,0,3,0)) plot(img, axes=FALSE,main="Athlete: Alexander Overwijk (CAN), Rank: 1") lines(c(p[,1],p[1,1]), c(p[,2],p[1,2]), col="magenta",lwd=2) points(p[,1], p[,2], cex=3, pch=20, col="steelblue") ##Points of the circle t <- seq(0,2*pi, length=nPoints) circle <- cbind(fit_tls[3]*cos(t) + fit_tls[1], fit_tls[3]*sin(t)+fit_tls[2]) circle[,1] <- circle[,1] + dx_blackboard ##Transform and show points of circle p_circle <- map.persp.inv(circle[,1], circle[,2], H=solve(H)) lines(p_circle$x, p_circle$y, col=rgb(0.8,0.8,0.8,0.5),lwd=3, lty=2) ##Transform and show center p_center <- map.persp.inv(fit_tls[1] + dx_blackboard, fit_tls[2], H=solve(H)) points(p_center$x, p_center$y, pch=20, cex=1, col=rgb(0.8,0.8,0.8,0.8)) ##Add HUD text text(p_center$x, p_center$y - 100, substitute(r[A] == x, list(x=sprintf("%.2f%%",100*ratio_area))), col=rgb(0.8,0.8,0.8,0.8)) text(p_center$x, p_center$y + 100, substitute(r[Delta * "A"] == x, list(x=sprintf("%.2f%%",100*ratio_areadifference))), col=rgb(0.8,0.8,0.8,0.8)) ``` ## Discussion We took elements of computer vision, image analysis and total least squares to segment a chalk-drawn circle on a blackboard and provided measures of it's circularness. Since we did not have direct access to the measurements of the blackboard in object space, a little guesstimation was necessary, nevertheless, the results show that it was a pretty circular freehand circle! With the machinery in place for judging freehand circles, its time to send out the call for contributions to the **2nd World Freehand Circle Drawing Championship** (online edition). Stay tuned for the call: participants would upload their photo plus minor modifications of a general analysis R-script computing the two area ratios measures and submit their contribution by a pull request to the github [WFHCDC repository](https://github.com/hoehleatsu/worldfreehandcirclechampionship). You can spend the anxious waiting time practicing your freehand 1m diameter circles - it's a good way to loosen up long and unproductive meetings! ## Appendix If we instead of the total sum of squares criterion involving $F_j(c,r)$ mentioned in the text solve the related criterion $$ \sum_{j=1}^m f_j(c,r)^2, \quad\text{where}\quad f_j(c,r) = ||c-a_j||_2^2 - r^2, $$ then a much simpler solution emerges. @coope1993 explains that this alternative criterion geometrically corresponds to minimizing the product $$ \text{(distance to the closest point on the circle)}\times \text{(distance to the furthest away point point on the circle)} $$ over the measurement point. In order to obtain the solution write the residuals $f_j$ as $f_j(c,r) = c^T c - 2 c^T a_j + a_j^T a_j - r^2$ and perform a change of variables from $(c_1, c_2, r)'$ to $$ y = \left[ \begin{matrix} 2 c_1 \\ 2 c_2 \\ r^2 - c^T c \\ \end{matrix} \right] \quad \text{and let} \quad b_j = \left[ \begin{matrix} a_{j1} \\ a_{j2} \\ 1 \end{matrix} \right]. $$ The minimization problem then becomes $$ \min_{y \in \mathbb{R}^3} \sum_{j=1}^m \left\{ a_j^T a_j - b_j^T y \right\}, $$ which can be written as a linear least square (LLS) expression $$ \min_y ||By - d||_2^2, $$ where $B$ is a $3\times m$ matrix with the $b_j$-vectors as columns and $d=||a_j||_2^2$. This expression is then easily solved using the standard least squares machinery. ```{r, echo=TRUE} ##Fast linear least squares problem as described in Coope (1993) fitCircle_lls <- function(freehandCircle) { a <- as.matrix(where(freehandCircle > 0)) b <- cbind(a,1) B <- b d <- a[,1]^2 + a[,2]^2 y <- solve(t(B) %*% B) %*% t(B) %*% d x <- 1/2*y[1:2] r <- as.numeric(sqrt(y[3] + t(x) %*% x)) return(c(x=x[1], y=x[2], radius=r)) } ##Fit using linear least squares procedure of Coole (1993) fit_lls <- fitCircle_lls(freehandCircleThinBorder) ##Compare TLS and LLS fit rbind(lls=fit_lls,tls=fit_tls) ``` In other words: the results are nearly identical. ```{r, eval=FALSE} layout(1) C_lls<- Circle(center=fit_lls[1:2], radius=fit_lls[3],border=1) plot(colorise(1-freehandCircleThinBorder, C_tls>0, col="blue", alpha=0.5)) highlight(C_lls, col="magenta") ``` ## Literature