dahtah / imager

R package for image processing
GNU Lesser General Public License v3.0
187 stars 43 forks source link

Recognizing an arc in an image? #132

Closed mconsidine closed 4 years ago

mconsidine commented 4 years ago

Hi, I am trying to use imager to extract an arc from an image and from that calculate the "circle of best fit" that would describe the arc (all I really need is either the radius or diameter of the circle associated with it). I am getting sort of close, but wonder if there is a better and/or different way that eludes me given the features in imager. A working example follows - thanks in advance for any help or pointers. mconsidine Example: `library(imager)

test_url <- "http://www.considine.net/test3.jpg"

im <- load.image(test_url) plot(im)

out <- threshold(im, "80%")

highlight(out)

how get the diameter of the "circle of best fit" for the red arc?

getting the diameter of a circle from a chord and a perpendicular:

https://math.stackexchange.com/questions/564058/calculate-the-radius-of-a-circle-given-the-chord-length-and-height-of-a-segment

temp <- highlight(out)

looks like x/y values in "temp" are for the highlighted outline???. If so, take averages

minxval <- (min(temp[[1]]$x)+min(temp[[2]]$x))/2 maxxval <- (max(temp[[1]]$x)+max(temp[[2]]$x))/2 minyval <- (min(temp[[1]]$y)+min(temp[[2]]$y))/2 maxyval <- (max(temp[[1]]$y)+max(temp[[2]]$y))/2 boxwidth <- abs( maxxval - minxval) boxheight <- abs( maxyval - minyval ) cornerlength <- sqrt(boxwidth^2+boxheight^2) #pixels midpt <- ceiling(length(temp[[1]]$x)/2) a <- (temp[[1]]$x[midpt]+temp[[2]]$x[midpt])/2 b <- (temp[[1]]$y[midpt]+temp[[2]]$y[midpt])/2 c <- ceiling(minxval+boxwidth/2) d <- ceiling(minyval+boxheight/2) e <- sqrt(abs(a-c)^2+abs(b-d)^2) #length of perpendicular f <- (4(e^2)+cornerlength^2)/(8(e)) #implied radius g <- 2*f #implied diameter

but this is subject to noise and is not quite what I'm looking for`

mconsidine commented 4 years ago

Looks like what I might want is the hough_circle function, though ideally without having to specify a radius ahead of time ...

mconsidine commented 4 years ago

Solved... sort of. Below is example code showing how I derived the coordinates of two chords from the image(s) using get.locations. The here creates one sample image and then takes a handful of subsets of the image to try to test how robust the approach is. I'm not quite sure why some derived values come out so widely different from the actual radius used. But since I'm taking a lot of samples and then looking at the median value, this seems like it will work for my purposes.

If someone knows of a more elegant way, I'd be interested to hear of it. (I could not use the hough_circle function because it requires you to specify a radius. And as-written I don't think it would work on arcs...)

`library(imager)

set.seed(123)

width <- 1900 height <- 1200

blankim <- cimg(array(0.2, c(width,height,1,3)))

x0 <- 900 y0 <- 600 aradius <- 431 blankim <- draw_circle(blankim, x0, y0, aradius, opacity = 1, color="WHITE",filled = TRUE) plot(blankim)

im1 <- imsub(blankim,x>300,x<500,y>300,y<500) im2 <- imsub(blankim,x>300,x<500,y<500) im3 <- imsub(blankim,x>300,x<500,y>500) im4 <- imsub(blankim,x>700,y>300,y<500) im5 <- imsub(blankim,x>700,y>300) im6 <- imsub(blankim,x>700,y<500) im7 <- imsub(blankim,x>950,y<500)

im <- im7 ########## START TEST ############### px <- cannyEdges(im, t1=0.1, t2=0.2, sigma=1) #args needed for test images plot(px) px2 <- get.locations(px,function(v) v==TRUE)

return_radius <- function(alist,i,j,k){

Method 3 from here:

https://www.qc.edu.hk/math/Advanced%20Level/circle%20given%203%20points.htm

mid1 <- array(c( (alist$x[i]+alist$x[j])/2, (alist$y[i]+alist$y[j])/2 ),2) mid2 <- array(c( (alist$x[j]+alist$x[k])/2, (alist$y[j]+alist$y[k])/2 ),2) grad1 <- (alist$y[j]-alist$y[i])/(alist$x[j]-alist$x[i]) grad2 <- (alist$y[k]-alist$y[j])/(alist$x[k]-alist$x[j])
gradL1 <- -1/grad1 gradL2 <- -1/grad2

ctrx <- ((mid2[2]-mid1[2])-(gradL1 (-1) mid1[1]-gradL2 (-1) mid2[1]))/(gradL1-gradL2) ctry <- gradL1 * (ctrx-mid1[1])+mid1[2]

ctrx, ctry = implied center

return(sqrt((alist$x[i]-ctrx)^2 + (alist$y[i]-ctry)^2)) #return radius }

numobs <-length(px2$x) numsamples <- 5000 estr <- array(NA,dim=numsamples) i<-1 while (i <= numsamples) { idx1 <- sample(1:ceiling(numobs 4/12),1) idx2 <- sample(ceiling(numobs 4/12):ceiling(numobs 8/12),1) idx3 <- sample(ceiling(numobs 8/12):numobs,1) tempval <- (return_radius(px2,idx1,idx2,idx3)) if (!is.na(tempval) && (abs(tempval)<5000)){ i <- i+1 estr[i] <- tempval } }

try to get at a good estimate of radius

hist(estr, breaks=2000) hist(estr[which(estr<2000)], breaks=2000) estr <- estr[which(estr<2000)] mu <- mean(estr) mdn <- median(estr) stdevn <- sd(estr) lowerbnd <- max(0,mu-3 stdevn) upperbnd <- min(2000,mu+3 stdevn) hist(estr[lowerbnd:upperbnd], breaks=2000) median(estr[lowerbnd:upperbnd],na.rm=TRUE)`