Generalized circles and inversions

[This article was first published on Saturn Elephant, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Let ι0 be the inversion with respect to the unit circle: ι0(z)=1ˉz,z0.

Consider a generalized circle GC having equation Azˉz+ˉγz+γˉz=D

which is equivalent, by dividing by zˉz, to A+ˉγ1ˉz+γ1z=Dzˉz.

Plugging z=ι10(w)=1/ˉw into this equation, we get the following equation of ι0(GC): A+ˉγw+γˉw=DwˉwDwˉwˉγwγˉw=A.

This is the equation of a generalized circle.

Now consider the inversion ι with respect to the circle of center c and radius r: ι(z)=c+r2¯zc,zc.

One has ι=f1ι0f

where f(z)=zcr, and f1(z)=rz+c.

The transformation f is a “scaling-translation”. It is clear that it maps a generalized circle to a generalized circle, and this also holds for f1. Therefore, ι maps any generalized circle to a generalized circle. This result is well-known. Here, we will give the equation of the image generalized circle.

Case when GC is a circle

Assume that GC is the circle C0 with center z0 and radius R0. Then f(C0) is the circle C1 with center z1=z0cr and radius R1=R0r. It has equation zˉz¯z1zz1ˉz=D1

with D1=R21|z1|2.

Then ι0(C1) is the generalized circle having equation D1zˉz+¯z1z+z1ˉz=1.

Case D10

If D10, the above equation of ι0(C1) is the equation of a circle C2 with center z2=z1D1 and radius R2=|z1|2D21+1D1.

Now, f1(C2)=ι(C0) is the circle with center rz2+c and radius rR2.

Case D1=0

If D1=0, our equation of ι0(C1) is ¯z1z+z1ˉz=1,

which is the equation of a line L2, with direction θ2=arg(z1) and offset d2=12|z1|2.

Then f1(L2)=ι(C0) is a line parallel to L2, that is, its direction is θ2. To get its offset, take a point w2 on L2: w2={12(z1)if (z1)=0i2(z1)if (z1)0

and take its image w3=f1(w2)=rw2+c. Then the offset is cos(θ2)(w3)+sin(θ2)(w3)=(z1)|z1|(w3)+(z1)|z1|(w3).

Case when GC is a line

Now assume that GC is the line L0 with direction θ0 and offset d0. Its image f(L0) is a line L1 parallel to L0, that is the direction of L1 is θ0. To get its offset, take a point w0 on L0: w0={d0cos(θ0)if sin(θ0)=0id0sin(θ0)if sin(θ0)0

and take its image w1=f(w0)=w0cr. Then the offset of L1 is d1=cos(θ0)(w1)+sin(θ0)(w1)
and the line L1 has equation ¯γ1z+γ1ˉz=D1
with γ1=cos(θ0)+isin(θ0) and D1=2d1.

Therefore, ι0(L1) is the generalized circle having equation D1zˉz¯γ1zγ1ˉz=0.

Case D10

If D10, then ι0(L1) is the circle C2 with center z2=γ1D1 and radius R2=|γ1||D1|. Then f1(C2)=ι(L0) is the circle with center rz2+c and radius rR2.

Case D1=0

If D1=0, then ι0(L1) is the line L2 having equation ¯γ1z+γ1ˉz=0.

Its direction is arg(γ1)=θ0 and its offset is 0.

Then f1(L2)=ι(C0) is a line L3 parallel to L2, that is, its direction is θ0. The line L2 passes through 0, therefore the point f1(0)=c belongs to L3. Consequently, the offset of L3 is cos(θ0)(c)+sin(θ0)(c).

R code

We firstly write a function iotaCircle which calculates the image of a circle.

# squared modulus
Mod2 <- function(z){
  Re(z)^2 + Im(z)^2
}
# image of a circle
iotaCircle <- function(c, r, circle, tol = sqrt(.Machine$double.eps)){
  z0 <- circle$center; R0 <- circle$radius
  z1 <- (z0-c)/r
  D1 <- (R0/r)^2 - Mod2(z1)
  if(abs(D1) > tol){
    z2 <- -z1/D1
    R2 <- sqrt(Mod2(z2) + 1/D1)
    out <- list(center = r*z2+c, radius = r*R2)
    attr(out, "type") <- "circle"
  }else{
    theta2 <- Arg(z1) %% (2*pi)
    w2 <- ifelse(Im(z1) == 0, 1/2/Re(z1), 1i/2/Im(z1))
    w3 <- r*w2 + c
    out <- list(
      theta = theta2, 
      offset = (Re(z1)*Re(w3)+Im(z1)*Im(w3))/Mod(z1)
    )
    attr(out, "type") <- "line"
  }
  out
}

Let’s check the function iotaCircle. To do so, we use the circumcircle function below, which returns the center and the radius of the circle passing by z1, z2, z3.

circumcircle <- function(z1,z2,z3){
  x1 <- Re(z1); y1 <- Im(z1); p1 <- c(x1,y1)
  x2 <- Re(z2); y2 <- Im(z2); p2 <- c(x2,y2)
  x3 <- Re(z3); y3 <- Im(z3); p3 <- c(x3,y3)
  a <- det(cbind(rbind(p1,p2,p3),1))
  q1 <- c(crossprod(p1)); q2 <- c(crossprod(p2)); q3 <- c(crossprod(p3))
  q <- c(q1,q2,q3)
  x <- c(x1,x2,x3)
  y <- c(y1,y2,y3)
  Dx <- det(cbind(q,y,1))
  Dy <- -det(cbind(q,x,1))
  c <- det(cbind(q,x,y))
  center <- c(Dx,Dy)/a/2
  r <- sqrt(c(crossprod(center-p1)))
  list(center = center[1] + 1i*center[2], radius = r)
}

Now let’s check iotaCircle.

# define iota
c <- 1 + 2i; r <- 3
iota <- function(z){
  c + r^2/Conj(z-c)
}
# define C0 - a case when iota(C0) is a circle ####
z0 <- 2 + 3i; R0 <- 2
iotaCircle(c, r, circle = list(center = z0, radius = R0))
## $center
## [1] -3.5-2.5i
## 
## $radius
## [1] 9
## 
## attr(,"type")
## [1] "circle"
circumcircle(iota(z0+R0), iota(z0-R0), iota(z0+1i*R0))
## $center
## [1] -3.5-2.5i
## 
## $radius
## [1] 9
# now, a case when iota(C0) is a line ####
R0 <- sqrt(2)
w1 <- iota(z0+R0); w2 <- iota(z0-R0); w3 <- iota(z0+1i*R0)
(w2-w1)/(w3-w1) # aligned <=> Im = 0
## [1] 3.414214+0i
(L <- iotaCircle(c, r, circle = list(center = z0, radius = R0)))
## $theta
## [1] 0.7853982
## 
## $offset
## [1] 5.303301
## 
## attr(,"type")
## [1] "line"
cos(L$theta)*Re(w1) + sin(L$theta)*Im(w1) # = L$offset
## [1] 5.303301
cos(L$theta)*Re(w2) + sin(L$theta)*Im(w2) # = L$offset
## [1] 5.303301
cos(L$theta)*Re(w3) + sin(L$theta)*Im(w3) # = L$offset
## [1] 5.303301

Now we write a function iotaLine which calculates the image of a line.

iotaLine <- function(c, r, line, tol = sqrt(.Machine$double.eps)){
  theta0 <- line$theta; d0 <- line$offset
  w0 <- ifelse(sin(theta0) == 0, d0/cos(theta0), 1i*d0/sin(theta0))
  w1 <- (w0-c)/r
  D1 <- 2 * (cos(theta0)*Re(w1) + sin(theta0)*Im(w1))
  if(abs(D1)>tol){
    gamma1 <- cos(theta0) + 1i*sin(theta0)
    z2 <- gamma1/D1
    R2 <- Mod(gamma1/D1)
    out <- list(center = r*z2+c, radius = r*R2)
    attr(out, "type") <- "circle"
  }else{
    out <- list(
      theta = theta0, 
      offset = cos(theta0)*Re(c)+sin(theta0)*Im(c)
    )
    attr(out, "type") <- "line"
  }
  out
}

Let’s check iotaLine. Firstly, a case when the image is a circle.

# define L0
theta0 <- 1; d0 <- 2
theta0 <- pi/4; d0 <- 0
# take three points on iota(L0)
x <- c(0,1,2)
y <- (d0 - cos(theta0)*x) / sin(theta0)
w1 <- iota(x[1] + 1i*y[1])
w2 <- iota(x[2] + 1i*y[2])
w3 <- iota(x[3] + 1i*y[3])
# check
circumcircle(w1,w2,w3)
## $center
## [1] -0.5+0.5i
## 
## $radius
## [1] 2.12132
iotaLine(c, r, list(theta = theta0, offset = d0))
## $center
## [1] -0.5+0.5i
## 
## $radius
## [1] 2.12132
## 
## attr(,"type")
## [1] "circle"

Now a case when the image is a line.

# define iota
c <- 1-1i 
iota <- function(z){
  c + r^2/Conj(z-c)
}
# define L0
theta0 <- pi/4; d0 <- 0
# take three points on iota(L0)
x <- c(0,2,4)
y <- (d0 - cos(theta0)*x) / sin(theta0)
w1 <- iota(x[1] + 1i*y[1])
w2 <- iota(x[2] + 1i*y[2])
w3 <- iota(x[3] + 1i*y[3])
# check
(w2-w1)/(w3-w1) # aligned <=> Im = 0
## [1] 1.5+0i
( L <- iotaLine(c, r, list(theta = theta0, offset = d0)) )
## $theta
## [1] 0.7853982
## 
## $offset
## [1] 1.110223e-16
## 
## attr(,"type")
## [1] "line"
cos(L$theta)*Re(w1) + sin(L$theta)*Im(w1) # = L$offset
## [1] -4.440892e-16
cos(L$theta)*Re(w2) + sin(L$theta)*Im(w2) # = L$offset
## [1] -4.440892e-16
cos(L$theta)*Re(w3) + sin(L$theta)*Im(w3) # = L$offset
## [1] -4.440892e-16

Finally, let’s put things together in a single function.

iotaGcircle <- function(c, r, gcircle){
  if(attr(gcircle, "type") == "circle"){
    iotaCircle(c, r, gcircle)
  }else{
    iotaLine(c, r, gcircle)
  }
}

Illustration

We will make a beautiful picture of inverted circles. This is borrowed from this circle inversion gallery.

Firstly, let’s write a function which draws a generalized circle. We use the draw.circle function of the plotrix package.

library(plotrix)
drawLine <- function(line, ...){
  theta <- line$theta; offset <- line$offset
  if(sin(theta) != 0){
    abline(a = offset/sin(theta), b = -1/tan(theta), ...)
  }else{
    abline(v = offset/cos(theta), ...)
  }
}
drawCircle <- function(circle, ...){
  draw.circle(Re(circle$center), Im(circle$center), circle$radius, ...)
}
drawGcircle <- function(gcircle, color = "black", ...){
  if(attr(gcircle, "type") == "circle"){
    drawCircle(gcircle, border = color, ...)
  }else{
    drawLine(gcircle, col = color, ...)
  }
}

We start with five circles.

# generation 0
gen0 <- sapply(c(0, pi/2, pi, 3*pi/2), function(beta){
  out <- list(center = cos(beta)+1i*sin(beta), radius = 1, gen = 0)
  attr(out, "type") <- "circle"
  out
}, simplify = FALSE)
gen0[[5]] <- list(center = 0+0i, radius = 2, gen = 0)
attr(gen0[[5]], "type") <- "circle"
# plot
par(bg = "black", mar = c(0,0,0,0))
plot(0, 0, type="n", xlim=c(-2.3,2.3), ylim=c(-2.3,2.3),   
     asp=1, axes=FALSE, xlab=NA, ylab=NA)
invisible(lapply(gen0, drawGcircle, color = "yellow", lwd = 2))

Now, we invert each of these five circles with respect to the other ones:

# generation 1
n0 <- length(gen0)
n1 <- n0*(n0-1)
gen1 <- vector("list", n1)
k <- 0
while(k < n1){
  for(j in 1:n0){
    for(i in 1:n0){
      if(i != j){
        k <- k+1
        gen1[[k]] <- 
          iotaGcircle(gen0[[i]]$center, gen0[[i]]$radius, gen0[[j]])
        gen1[[k]]$base <- i
        gen1[[k]]$gen <- 1
      }
    }
  }
}

We continue so on: we invert the obtained generalized circles with respect to the starting circles.

# generation 2
n2 <- n0*n1-n1
gen2 <- vector("list", n2)
k <- 0
while(k < n2){
  for(j in 1:n1){
    for(i in 1:n0){
      if(gen1[[j]]$base != i){
        k <- k+1
        gen2[[k]] <- 
          iotaGcircle(gen0[[i]]$center, gen0[[i]]$radius, gen1[[j]])
        gen2[[k]]$base <- i
        gen2[[k]]$gen <- 2
      }
    }
  }
}
# generation 3
n3 <- n0*n2-n2
gen3 <- vector("list", n3)
k <- 0
while(k < n3){
  for(j in 1:n2){
    for(i in 1:n0){
      if(gen2[[j]]$base != i){
        k <- k+1
        gen3[[k]] <- 
          iotaGcircle(gen0[[i]]$center, gen0[[i]]$radius, gen2[[j]])
        gen3[[k]]$gen <- 3
      }
    }
  }
}

Let’s put all the obtained generalized circles together:

gcircles <- c(gen0, gen1, gen2, gen3)
length(gcircles)
## [1] 425

There are 425 generalized circles, but some of them are duplicated. To remove the duplicates, we use the uniqueWith function below. This function takes as arguments an atomic vector or a list v and a function f of two arguments, two elements of v, and which returns TRUE or FALSE, according to whether the two elements are considered as duplicates.

uniqueWith <- function(v, f){
  size <- length(v)
  for(i in seq_len(size-1L)){
    j <- i + 1L
    while(j <= size){
      if(f(v[[j]],v[[i]])){
        v <- v[-j]
        size <- size - 1L
      }else{
        j <- j + 1L
      }
    }
  }
  v[1L:size]
}
# examples ####
v <- c(1,1,2,1,3,4,3,5,5)
uniqueWith(v, `==`)
## [1] 1 2 3 4 5
uniqueWith(v, function(x,y) (x-y) %% 3 == 0)
## [1] 1 2 3
v <- list(a="you", b="are", c="great")
uniqueWith(v, function(x,y) nchar(x) == nchar(y))
## $a
## [1] "you"
## 
## $c
## [1] "great"

So let’s define the function f which identifies two equal generalized circles.

f <- function(gcircle1, gcircle2){
  if(attr(gcircle1, "type") == attr(gcircle2, "type")){
    if(attr(gcircle1, "type") == "circle"){
      Mod(gcircle1$center-gcircle2$center) < 1e-5 && 
        abs(gcircle1$radius-gcircle2$radius) < 1e-5
    }else{
      abs(cos(gcircle1$theta)-cos(gcircle2$theta)) < 1e-5 &&
        abs(sin(gcircle1$theta)-sin(gcircle2$theta)) < 1e-5 &&
        abs(gcircle1$offset-gcircle2$offset) < 1e-5
    }
  }else{
    FALSE
  }
}

And now, let’s remove the duplicates:

gcircles <- uniqueWith(gcircles, f)
length(gcircles)
## [1] 161

Now we plot the generalized circles, with a color indicating the generation.

draw <- function(gcircle, colors=rainbow(4), ...){
  drawGcircle(gcircle, color = colors[1+gcircle$gen], ...)
}
par(mar = c(0,0,0,0), bg = "black")
plot(0, 0, type="n", xlim=c(-2.3,2.3), ylim=c(-2.3,2.3),  
     asp=1, axes=FALSE, xlab=NA, ylab=NA)
invisible(lapply(gcircles, draw, lwd=2))

To leave a comment for the author, please follow the link and comment on their blog: Saturn Elephant.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)