Draw a Free Hand Circle

Abstruse:

In 2007 Alexander Overwijk went viral with his 'Perfect Circumvolve' video. The same year a World Freehand Circumvolve Drawing Championship was organized, which he won. In this mail we testify how a mobile camera, R and the imager packet can be used to develop an image analysis based method to guess future instances of the championship.

This work is licensed nether a Creative Eatables Attribution-ShareAlike 4.0 International License. The markdown+Rknitr source lawmaking of this blog is available under a GNU Full general Public License (GPL v3) license from github.

Introduction

A few years back I watched with awe the 2007 video of Alexander Overwijk's freehand drawing a 1m diameter circle:

Note: Depending on your browser you might need to click "Scout on Youtube" to see the video.

Ever since watching that video I have wondered how one would go about to judge the winner of such an alleged Globe Freehand Circle Drawing Title (WFHCDC). While researching for this post I finally figured it out. On his web page Alexander in the story backside the video "reveals":

They have a light amplification by stimulated emission of radiation machine called the circleometer that creates the perfect circumvolve closest to the one you lot drew. The circleometer then calculates the divergence in surface area between the laser circle and the circumvolve that yous drew. The machine and then calibrates the area deviation as if you had fatigued a circle with radius one meter. The person with the smallest area difference is alleged the world freehand circle drawing champion.

Aha! Imaginary circleometers are expensive and my dean of study well-nigh likely isn't open for investing in perfect circle measurement equipment… So here is a cheaper solution involving a mobile device camera, R and the imager package past Simon Barthelmé et al. Altogether a combination of modern data science tools, which my dean of study is nearly probable to approve! Nosotros'll utilise a screenshot from the perfect circle video every bit motivating example to guide through the iii phases of the method:

1. Epitome Rectification
2. Freehand circumvolve identification and perfect circle estimation
3. Quantifying deviation from the perfect circle

We beginning by loading the screenshot into R using imager:

library("imager") file <- "circle2.png" img <- imager::load.image(file.path(fullFigPath, file))

Image Rectification

The image clearly suffers from perspective distortions acquired by the photographic camera being positioned to the right of the circle and, hence, not being orthogonal to the blackboard plane. Furthermore, small lense distortions are as well visible - for case 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, nosotros will ignore lense distortions in this postal service. If such information is bachelor one can utilize a programme such as RawTherapee (available under a GNU GPL v3 license) to read the EXIF information in the meta information of the paradigm and automatically right for lens baloney.

To rectify the image we judge the parameters of the 2D projection based on iv ground control points (GPC). We use R's locator role to determine the pixel location of the iv corner points of the blackboard in the image, but could but as well utilize any image analysis program such every bit Gimp. Furthermore, nosotros need the truthful object coordinates of these GPC. Unfortunately, these are only approximately available to due lack of noesis of the size of the blackboard in the classroom. As a issue a guesstimate of the horizontal length is used.

plot(img) p <- locator(4) p <- circular(cbind(p$x, p$y)) dump(list=c("p"), "")

These points are at present used to rectify the epitome past a Straight Linear Transformation (DLT) based on exactly 4 control points (Hartley and Zisserman 2004, Affiliate 4) ¹. That is the parameters of the 3x3 transformation matrix \(H\) in homogeneous coordinates are estimated such that \(p' = H p\), meet the lawmaking on github for details.

We can implement the rectifying transformation using the imager::warp function:

##Transform image coordinates (x',y') to (10,y), i.eastward. 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(ten,y,1)   list(x=out_image[i,]/out_image[3,], y=out_image[two,]/out_image[3,]) } ##Pad dx_blackboard pixels to the right to make infinite for blackboard img_padded <- pad(img, nPix=dx_blackboard, axes="x", pos=i) ##Warp prototype warp <- imwarp(img_padded, map=function(ten,y) map.persp.inv(ten,y,solve(H)),coordinates="absolute", management="backward")

The effect looks equally follows: Please notice the dissimilar x-axes of the two images when comparing them. For faster computation and meliorate visualization in the remainder of this post, we crop the x-axis of the image to the relevant parts of the circle.

warp <- imsub(warp, ten %inr% c(dx_blackboard, nrow(warp)))

Freehand circumvolve identification

As described in the imager edge detection tutorial we use length of the gradient to decide the edges in the image. This tin be done by applying filters to the image.

##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 algorithm with a priority map inverse proportional to gradient magnitude. This includes a few points outside the circumvolve and a few points inside the circumvolve. Note: a perfect circumvolve would have no edge, but when drawing a circle with a piece of chalk it's destined to have a thin edge line.

Nosotros can now extract the circumvolve by:

##Merely the circle freehandCircle <- (warp * (mask==two) > 0) %>% grayscale ##Total expanse covered past the circle freehandDisc <- label(freehandCircle, high_connectivity=TRUE) > 0 dilatedDisc <- freehandDisc %>% dilate_rect(sx=3,sy=iii) freehandCircleThinBorder <- (freehandDisc - dilatedDisc) != 0

Perfect circle interpretation

Once the freehand circumvolve path in the prototype has been identified, we need to notice the best fitting perfect circle matching this path. This trouble is elegantly solved past Coope (1993), who formulates the trouble every bit finding center and radius of the circle minimizing the squared Euclidean distance to \(m\) information points \(a_j\), \(j=1,\ldots,m\). Cogent by \(c\) the center of the circle and by \(r>0\) the radius we want to detect the solution of

\[ \min_{c\in \mathbb{R^two}, r>0} \sum_{j=1}^m F_j(c,r)^ii, \quad\text{where}\quad F_j(c,r) = \left|r - ||c-a_j||_2\right|, \]

and \(||ten||_2\) denotes Euclidean distance. Because the curve fitting minimizes the distance betwixt an observed point \(a_j\) and its closest point on the circle and thus involves both the \(x\) and the \(y\) management , this is a so called full least squares problem. The problem is non-linear and can merely exist solved by iterative numerical methods. Yet, the dimension of the parameter infinite can be reduced by one, because given the centre \(c\) we can determine that \(r(c)=\frac{one}{m} \sum_{j=1}^k ||c-a_j||_2\).

##Compute radius given center radius_given_center <- function(center, dist=NULL) {   if (is.nil(dist)) {     a <- equally.matrix(where(freehandCircleThinBorder > 0))     dist <- sqrt((a[,1] - center[1])^2 + (a[,2] - center[2])^two)   }   return(mean(dist)) }  ##Target functin of the full to the lowest degree squares criterion of Coope (1993) target_tls <- function(theta) {   ##Extract parameters   middle <- exp(theta[1:two])    ##Total least squares criterion from Coope (1993)   a <- as.matrix(where(freehandCircleThinBorder > 0))   dist <- sqrt((a[,1] - heart[i])^2 + (a[,2] - centre[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,one], y=groundwork[1,2])), fn=target_tls) middle <- exp(res_tls$par) fit_tls <- c(center,radius=radius_given_center(center)) fit_tls ##        ten        y   radius  ## 894.3885 707.3191 518.4119

We illustrate the freehand circle (in black) and the fitted circle (magenta) on summit of each other using the alpha channel. You have to study the image carefully to observe differences between the two curves!

Quantifying the circularness of the freehand circumvolve

We quantify the circularness of the freehand circle past contrasting the expanse covered by it with the surface area of the fitted perfect circumvolve. The closer this ratio is to 1 the more perfect is the freehand circle.

##Area of the freehand drawn disc areaFreehandDisc <- sum(freehandDisc)  ##Area of the disc corresponding to the idealized circumvolve fitted ##to the freehand circle areaIdealDisc <- pi * fit_tls["radius"]^2  ##Ratio between the two areas ratio_area <- as.numeric(areaFreehandDisc  / areaIdealDisc) ratio_area ## [1] 0.9971778

Yup, it's a pretty perfect circumvolve! Note also that the scale is to a circle with an area of 1 pixel unit and not a circumvolve with diameter of 1m as described in the above text about the "circleometer". Since the fitted circumvolve already takes the desired shape into account, my intuition is that this ratio is a pretty skillful way to quantify circularness. However, to avert measurehacks, nosotros use every bit fill-in measure the circleometer arroyo: for each point on the freehand circle we mensurate its distance to the freehand circle and integrate/sum this upward over the path of the freehand circle. Nosotros can guess this integration using image pixels as follows.

##Create a pixel based circle in an image of the same size as the ##freehandCircle img. For visibility we apply a border of 'border' pixels ##s.t. circumvolve goes [radius - border/ii, radius + edge/two]. Circle <- function(middle, radius, edge) {   equally.cimg(role(x,y) {     lhs <- (x-eye[1])^2 + (y-center[2])^2     return( (lhs >= (radius-border/2)^2) & (lhs <= (radius+border/2)^ii))   }, dim=dim(freehandCircle)) }  ##Build pixel circle based on the fitted parameters C_tls <- Circle(fit_tls[one:2], fit_tls[3], border=1) ##Calculate Euclidean altitude to circle for each pixel in the image dist <- distance_transform(C_tls, value=one, metric=2) ##Distance between outer border of freehand circle and perfect circle area_difference <- sum(dist[freehandCircleThinBorder>0])  ##Compute area difference and scaled it past the area of the fitted disc ratio_areadifference <- as.numeric(area_difference / areaIdealDisc)

The image beneath illustrates this by overlaying the result on top of the distance map. For meliorate 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 surface area between the two paths is summed upwardly over the entire path of the freehand circumvolve:

Nosotros obtain ratio_areadifference= 0.01735. Thus also this measure out tells us: it'due south a pretty perfect circumvolve! To summarise: The output on the brandish of the guess'southward Circle-O-Meter App (bachelor under a GPL v3 license) at the World Freehand Circumvolve Cartoon Title would be equally follows: ?

Discussion

Nosotros took elements of computer vision, paradigm analysis and full least squares to segment a chalk-fatigued 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 lilliputian guesstimation was necessary, notwithstanding, the results show that it was a pretty circular freehand circle!

With the machinery in identify 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 the github WFHCDC repository. You lot tin spend the anxious waiting time practicing your freehand 1m diameter circles - information technology's a good way to loosen upward long & unproductive meetings!

Appendix

If we instead of the full sum of squares criterion involving \(F_j(c,r)\) mentioned in the text solve the related benchmark \[ \sum_{j=1}^1000 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. Coope (1993) explains that this alternative criterion geometrically corresponds to minimizing the product

\[ \text{(distance to the closest point on the circumvolve)}\times \text{(distance to the furthest away bespeak signal on the circle)} \]

over the measurement point. In social club 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 \\ \stop{matrix} \right] \quad \text{and let} \quad b_j = \left[ \begin{matrix} a_{j1} \\ a_{j2} \\ 1 \terminate{matrix} \right]. \] The minimization problem then becomes \[ \min_{y \in \mathbb{R}^3} \sum_{j=one}^m \left\{ a_j^T a_j - b_j^T y \right\}, \] which can exist written as a linear least square (LLS) expression \[ \min_y ||By - d||_2^2, \] where \(B\) is a \(3\times g\) matrix with the \(b_j\)-vectors as columns and \(d=||a_j||_2^2\). This expression is then easily solved using the standard to the lowest degree squares machinery.

##Fast linear to the lowest degree squares problem as described in Coope (1993) fitCircle_lls <- function(freehandCircle) {   a <- equally.matrix(where(freehandCircle > 0))   b <- cbind(a,1)   B <- b   d <- a[,one]^2 + a[,ii]^2   y <- solve(t(B) %*% B) %*% t(B) %*% d   x <- 1/2*y[i:2]   r <- as.numeric(sqrt(y[iii] + t(ten) %*% x))   return(c(x=ten[1], y=ten[2], radius=r)) }  ##Fit using linear least squares process of Coole (1993) fit_lls <- fitCircle_lls(freehandCircleThinBorder)  ##Compare TLS and LLS fit rbind(lls=fit_lls,tls=fit_tls) ##            x        y   radius ## lls 894.3666 707.3901 518.4295 ## tls 894.3885 707.3191 518.4119

In other words: the results are nearly identical.

Literature

Coope, I.D. 1993. "Circumvolve Fitting by Linear and Nonlinear Least Squares." Journal of Optimization Theory and Applications 76 (2): 381–88. doi:10.1007/BF00939613.

---

1. Alternatively, run into slide 18 and onward in https://ags.cs.uni-kl.de/fileadmin/inf_ags/3dcv-ws11-12/3DCV_WS11-12_lec04.pdf↩

thackertreff1945.blogspot.com

Source: https://www.r-bloggers.com/2018/07/judging-freehand-circle-drawing-competitions/

Share this post