Nothing
#' Opens default web browser and loads a web page for area estimation.
#'
#' Function \code{startweb.area()} opens a web browser and loads the page
#' for area estimation of object shown in a picture. Use this function if you
#' can isolate a part of the picture with uniformly distributed dots on an
#' elliptical domain with unclear borders. Sequence of actions on the web page
#' is as follows:
#' \enumerate{
#' \item Load a picture in JPG format
#' \item Click on upper left and lower right corner of a rectangle surrounding
#' observed object so the rectangle is drawn
#' \item Set data and estimation parameters
#' \item Click on Estimate
#' }
#'
#' The area estimation algorithm takes many horizontal and vertical (if
#' "horizontal + vertical" slicing is selected) or star-shaped (if "star"
#' slicing is selected) slices of the object. Length estimation procedure is
#' conducted on each slice and in that way set of edge points is obtained.
#' Lastly, ellipse or circle is fitted on that set of points by function
#' \code{\link[conicfit]{EllipseDirectFit}} or
#' \code{\link[conicfit]{CircleFitByPratt}} from the package \code{conicfit}
#' and area of that ellipse or circle is returned as the result. The area is
#' measured in pixels, as well as percentage of the whole image.
#'
#' Parameters that can be set on the web page are as follows:
#'
#' \strong{Data parameters}
#' \describe{
#' \item{Levels of grey}{Number of colors (shades of grey) used in
#' analysis.}
#' \item{Box size}{The algorithm takes each pixel of a picture and maps it
#' to box_size * box_size matrix. It is done in a way that the brightness
#' of the observed pixel dictates the quantity of dots in mentioned
#' matrix. Distribution of dots in matrix is uniform. Ultimately, length
#' estimation is done on the set of the resulting matrices.}
#' \item{Line thickness}{Width of the slice, i.e. the maximum length between
#' surrounding pixel and the drawn line so that pixel is to be taken into
#' account for length estimation. All surrounding pixels are orthogonally
#' projected on the central line.}
#' \item{Number of slices}{Number of slices after cutting in one direction.
#' Defaults to 10. Slices are equally thick in both directions. Smaller
#' number of cuts will be automatically applied for smaller dimension if
#' the chosen rectangle is not a square.}
#' \item{Slicing}{Sets slicing method for the edge point estimation. Can be
#' "horizontal + vertical" or "star".}
#' \item{Parallelization}{Sets whether to distribute area estimation on
#' multiple CPU cores. If set to On, total number of cores - 1 are used.}
#' \item{Object brightness}{Sets whether observed object is bright or
#' dark.}
#' \item{Represent object as}{Represent estimated object as an ellipse or as
#' a circle.}
#' }
#' \strong{Estimation parameters}
#' \describe{
#' \item{Error distribution}{Type of the error distribution. Can be Gauss,
#' Laplace, T1, T2, T3, T4 or T5 (Student).}
#' \item{Error standard deviation}{Estimation method for the error
#' standard deviation. Can be Maximum Likelihood (ML) or the Method of
#' Moments. If one does not want to estimate the deviation but to
#' explicitly enter it, he should choose "Enter value" and enter the
#' deviation in the lower field.}
#' }
#'
#' @section Note:
#' In order to have quadratic pixels on the screen, please use proportional
#' screen resolution. In the case of modern LCD (LED) displays, these are
#' usually native screen resolutions. If your display has aspect ratio
#' width:height = 16:9, these resolutions are 1280x720, 1600x900, 1920x1080,
#' etc. In the case od 16:10 display, use 1280x800, 1440x900, 1920x1200, etc.
#' If you use nonproportional screen resolution, pixels on the screen will not
#' be quadratic, so estimated values measured in pixels may not be correct.
#'
#' @examples
#' # open the web page for area estimation of an object shown in the picture
#' \donttest{startweb.area()}
#'
#' @export
startweb.area <- function() {
opencpu::ocpu_start_server(preload = "LeArEst", on_startup = function(server){
utils::browseURL(paste0(server, "/library/LeArEst/www/index_area.html"))
})
}
# Utility function for web interface call
# Not intended for user calls!
#
#' @importFrom stats runif
#' @importFrom stats kmeans
load.image.area <- function(picture, dataBright = "bright", representation,
error, var, varEst, confLevel, levelsOfGray, boxSize, thickness, parallel,
slicing, nrSlices, startX, startY, endX, endY, verbose = TRUE) {
# 0 < intensity < levelsOfGray-1 must apply: boxSize^2 >= levelsOfGray-1 !!!
var <- as.numeric(var)
conf.level <- as.numeric(confLevel)
levelsOfGray <- as.numeric(levelsOfGray)
boxSize <- as.numeric(boxSize)
thickness <- as.numeric(thickness)
nrSlices <- as.numeric(nrSlices)
parallel <- as.logical(as.numeric(parallel))
img_temp <- jpeg::readJPEG(picture)
# if the image is color, convert it to greyscale first
if (!is.na(dim(img_temp)[3])) {
img <- img_temp[, , 1] * 0.33 + img_temp[, , 2] * 0.33 +
img_temp[, , 3] * 0.33
} else {
img <- img_temp
}
rm(img_temp)
# if we watch the dark area, make the negative of the image
if (dataBright == "dark") {
img <- 1 - img
}
img <- floor(levelsOfGray * img)
height <- dim(img)[1]
width <- dim(img)[2]
startPoint <- c(as.numeric(startX), as.numeric(startY))
endPoint <- c(as.numeric(endX), as.numeric(endY))
# if the endPoint is 'more left' then the startPoint - switch them
if (endPoint[1] < startPoint[1]) {
temp <- startPoint
startPoint <- endPoint
endPoint <- temp
}
# correction of the points, if they are too close to the picture borders
if (startPoint[1] < thickness/2)
startPoint[1] <- ceiling(thickness/2)
if (startPoint[2] < thickness/2)
startPoint[2] <- ceiling(thickness/2)
if (startPoint[1] > height - thickness/2)
startPoint[1] <- height - ceiling(thickness/2)
if (startPoint[2] > width - thickness/2)
startPoint[2] <- width - ceiling(thickness/2)
if (endPoint[1] < thickness/2)
endPoint[1] <- ceiling(thickness/2)
if (endPoint[2] < thickness/2)
endPoint[2] <- ceiling(thickness/2)
if (endPoint[1] > height - thickness/2)
endPoint[1] <- height - ceiling(thickness/2)
if (endPoint[2] > width - thickness/2)
endPoint[2] <- width - ceiling(thickness/2)
# begin contrast optimization
minBrightness <- min( img[startPoint[1]:endPoint[1], startPoint[2]:endPoint[2]] )
maxBrightness <- max( img[startPoint[1]:endPoint[1], startPoint[2]:endPoint[2]] )
img <- levelsOfGray * (img - minBrightness) / (maxBrightness - minBrightness)
img[img<0] <- 0
# end contrast optimization
# prepare dataMatrix for the whole image
dataMatrix = matrix(FALSE, nrow = height * boxSize, ncol = width * boxSize)
set.seed(123)
# for (j in 1:width) {
# for (i in 1:height) {
for (j in startPoint[2]:endPoint[2]) {
for (i in startPoint[1]:endPoint[1]) {
x <- floor(runif(img[i, j], max = boxSize))
y <- floor(runif(img[i, j], max = boxSize))
dataMatrix[(i - 1) * boxSize + x, (j - 1) * boxSize + y] <- TRUE
}
} # TODO: za manje zauzece memorije smanjiti dimenziju dataMatrixa
# divide selected square into two clusters to obtain limits for chopping
kMeans <- kmeans(c(img[startPoint[1]:endPoint[1], startPoint[2]:endPoint[2]]),
2, nstart = 100)
clusterMatrix <- matrix(kMeans$cluster, nrow = endPoint[1] - startPoint[1] + 1,
ncol = endPoint[2] - startPoint[2] + 1)
# remove the outliers from clusterMatrix 'consecutive' tells how many
# consecutive rows or columns should be 'an object' to trigger clusterStartPoint
consecutive <- 20
# suppose there is no object on the selected square's borders.
# Taking that into account, adjust object to be in cluster 2.
borderClusterPoints <- c(clusterMatrix[1, ], clusterMatrix[, dim(clusterMatrix)[2]],
clusterMatrix[dim(clusterMatrix)[1], ], clusterMatrix[, 1])
if (sum(borderClusterPoints == 2) > sum(borderClusterPoints == 1)) {
clusterMatrix <- clusterMatrix + 1
clusterMatrix[clusterMatrix == 3] <- 1
}
# image(clusterMatrix)
maxrows <- apply(clusterMatrix, 1, max)
maxcols <- apply(clusterMatrix, 2, max)
for (i in 1:(length(maxrows) - consecutive + 1)) {
if (sum(maxrows[i:(i + consecutive - 1)] == 2) == consecutive) {
clusterMatrix[1:(i - 1), ] <- 1
break
}
}
for (i in (length(maxrows) - consecutive + 1):1) {
if (sum(maxrows[i:(i + consecutive - 1)] == 2) == consecutive) {
clusterMatrix[(i + consecutive - 1):length(maxrows), ] <- 1
break
}
}
for (i in 1:(length(maxcols) - consecutive + 1)) {
if (sum(maxcols[i:(i + consecutive - 1)] == 2) == consecutive) {
clusterMatrix[, 1:(i - 1)] <- 1
break
}
}
for (i in (length(maxcols) - consecutive + 1):1) {
if (sum(maxcols[i:(i + consecutive - 1)] == 2) == consecutive) {
clusterMatrix[, (i + consecutive - 1):length(maxcols)] <- 1
break
}
}
watchedPoints <- which(clusterMatrix == 2, arr.ind = TRUE)
clusterStartPoint <- c(min(watchedPoints[, 1]), min(watchedPoints[, 2])) +
startPoint # -1
clusterEndPoint <- c(max(watchedPoints[, 1]), max(watchedPoints[, 2])) +
startPoint # - 1
ellipsePoints <- matrix(NA, nrow = 1000, ncol = 2)
pointNr <- 1
biggerDim <- max((clusterEndPoint[1] - clusterStartPoint[1]),
(clusterEndPoint[2] - clusterStartPoint[2]))
lineDistance <- biggerDim %/% nrSlices
rectHeight <- endPoint[1] - startPoint[1]
rectWidth <- abs(endPoint[2] - startPoint[2])
currentX <- floor((clusterStartPoint[1] + lineDistance/2) * boxSize)
start.time <- Sys.time()
if (slicing == "hv") {
# horizontal chopping
if (parallel) {
no_cores <- parallel::detectCores() - 1
cl <- parallel::makeCluster(no_cores)
doParallel::registerDoParallel(cl)
ellipsePointsH <- foreach::foreach (currentX = seq(
from = floor((clusterStartPoint[1] + lineDistance/2) * boxSize),
to = clusterEndPoint[1] * boxSize,
by = lineDistance * boxSize)) %dopar% {
distancesVector = vector(mode = "numeric",
length = thickness * boxSize * rectWidth)
distFromStart <- 0
cnt <- 1
for (currentY in min(startPoint[2] * boxSize, endPoint[2] * boxSize):
max(startPoint[2] * boxSize, endPoint[2] * boxSize)) {
leftestPoint <- currentX - floor(thickness * boxSize / 2)
rightestPoint <- currentX + floor(thickness * boxSize / 2)
sumThisY <- sum(dataMatrix[leftestPoint:rightestPoint, currentY])
if (sumThisY > 0) {
distancesVector[cnt:(cnt + sumThisY - 1)] <- distFromStart
}
distFromStart <- distFromStart + 1
cnt <- cnt + sumThisY
}
data <- distancesVector[distancesVector != 0]
dataCenter <- mean(data)
data <- (data - dataCenter) / boxSize
if (varEst == "value") {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = varEst,
conf.level = conf.level)
}
centerY <- (min(startPoint[2], endPoint[2]) + dataCenter / boxSize)
leftPoint <- floor(c(currentX / boxSize,
centerY - unname(estResult$radius)))
rightPoint <- floor(c(currentX / boxSize,
centerY + unname(estResult$radius)))
# rm(distancesVector, data, dataCenter, data, estResult, centerY)
c(leftPoint, rightPoint)
}
# razvrti tu listu u matricu:
ellipsePoints <- matrix(unlist(ellipsePointsH), ncol = 2, byrow = TRUE)
} else {
while (currentX < clusterEndPoint[1] * boxSize) {
distancesVector = vector(mode = "numeric",
length = thickness * boxSize * rectWidth)
distFromStart <- 0
cnt <- 1
for (currentY in min(startPoint[2] * boxSize, endPoint[2] * boxSize):
max(startPoint[2] * boxSize, endPoint[2] * boxSize)) {
leftestPoint <- currentX - floor(thickness * boxSize / 2)
rightestPoint <- currentX + floor(thickness * boxSize / 2)
sumThisY <- sum(dataMatrix[leftestPoint:rightestPoint, currentY])
if (sumThisY > 0) {
distancesVector[cnt:(cnt + sumThisY - 1)] <- distFromStart
}
distFromStart <- distFromStart + 1
cnt <- cnt + sumThisY
}
data <- distancesVector[distancesVector != 0]
dataCenter <- mean(data)
data <- (data - dataCenter) / boxSize
if (varEst == "value") {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = varEst,
conf.level = conf.level)
}
centerY <- (min(startPoint[2], endPoint[2]) + dataCenter / boxSize)
leftPoint <- floor(c(currentX / boxSize, centerY - estResult$radius))
rightPoint <- floor(c(currentX / boxSize, centerY + estResult$radius))
ellipsePoints[(2 * pointNr - 1), ] <- leftPoint
ellipsePoints[(2 * pointNr), ] <- rightPoint
currentX <- currentX + lineDistance * boxSize
pointNr <- pointNr + 1
}
}
currentY <- floor((clusterStartPoint[2] + lineDistance / 2) * boxSize)
# vertical chopping
if (parallel) {
ellipsePointsV <- foreach::foreach (currentY = seq(
from = floor((clusterStartPoint[2] + lineDistance / 2) * boxSize),
to = clusterEndPoint[2] * boxSize,
by = lineDistance * boxSize)) %dopar% {
distancesVector = vector(mode = "numeric",
length = thickness * boxSize * rectHeight)
distFromStart <- 0
cnt <- 1
for (currentX in min(startPoint[1] * boxSize, endPoint[1] * boxSize):
max(startPoint[1] * boxSize, endPoint[1] * boxSize)) {
leftestPoint <- currentY - floor(thickness * boxSize / 2)
rightestPoint <- currentY + floor(thickness * boxSize / 2)
sumThisX <- sum(dataMatrix[currentX, leftestPoint:rightestPoint])
if (sumThisX > 0) {
distancesVector[cnt:(cnt + sumThisX - 1)] <- distFromStart
}
distFromStart <- distFromStart + 1
cnt <- cnt + sumThisX
}
data <- distancesVector[distancesVector != 0]
dataCenter <- mean(data)
data <- (data - dataCenter) / boxSize
if (varEst == "value") {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = varEst,
conf.level = conf.level)
}
centerX <- (min(startPoint[1], endPoint[1]) + dataCenter / boxSize)
upPoint <- floor(c(centerX - estResult$radius, currentY / boxSize))
downPoint <- floor(c(centerX + estResult$radius, currentY / boxSize))
# rm(distancesVector, data, dataCenter, data, estResult, centerY)
c(upPoint, downPoint)
}
ellipsePointsV <- matrix(unlist(ellipsePointsV), ncol = 2, byrow = TRUE)
ellipsePoints <- as.matrix(merge(ellipsePoints, ellipsePointsV,
all = TRUE))
} else {
while (currentY < clusterEndPoint[2] * boxSize) {
distancesVector = vector(mode = "numeric",
length = thickness * boxSize * rectHeight)
distFromStart <- 0
cnt <- 1
for (currentX in min(startPoint[1] * boxSize, endPoint[1] * boxSize):
max(startPoint[1] * boxSize, endPoint[1] * boxSize)) {
leftestPoint <- currentY - floor(thickness * boxSize / 2)
rightestPoint <- currentY + floor(thickness * boxSize / 2)
sumThisX <- sum(dataMatrix[currentX, leftestPoint:rightestPoint])
if (sumThisX > 0) {
distancesVector[cnt:(cnt + sumThisX - 1)] <- distFromStart
}
distFromStart <- distFromStart + 1
cnt <- cnt + sumThisX
}
data <- distancesVector[distancesVector != 0]
dataCenter <- mean(data)
data <- (data - dataCenter) / boxSize
if (varEst == "value") {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = varEst,
conf.level = conf.level)
}
centerX <- (min(startPoint[1], endPoint[1]) + dataCenter / boxSize)
upPoint <- floor(c(centerX - estResult$radius, currentY / boxSize))
downPoint <- floor(c(centerX + estResult$radius, currentY / boxSize))
ellipsePoints[(2 * pointNr - 1), ] <- upPoint
ellipsePoints[(2 * pointNr), ] <- downPoint
currentY <- currentY + lineDistance * boxSize
pointNr <- pointNr + 1
}
}
} else if (slicing == "star") {
# common part (parallel and serial):
# create matrix of "TRUE" locations in dataMatrix, filter it to data
# cluster and translate in a way that clusterStartPoint * boxSize is the
# origin of the coordinate system
mydatafull <- as.data.frame(which(dataMatrix == TRUE, arr.ind=T))
colnames(mydatafull) <- c("x", "y")
mydata <- mydatafull[mydatafull$x > clusterStartPoint[1] * boxSize &
mydatafull$y > clusterStartPoint[2] * boxSize &
mydatafull$x < clusterEndPoint[1] * boxSize &
mydatafull$y < clusterEndPoint[2] * boxSize, ]
mydata$x <- mydata$x - (clusterStartPoint[1] * boxSize)
mydata$y <- mydata$y - (clusterStartPoint[2] * boxSize)
rm(mydatafull)
# translate input points so the origin is data mean
# centerx <- mean (mydata$x)
# centery <- mean (mydata$y)
# translate input points so the origin is center point of the data cluster
centerx <- abs((clusterEndPoint[1] - clusterStartPoint[1]) * boxSize %/% 2)
centery <- abs((clusterEndPoint[2] - clusterStartPoint[2]) * boxSize %/% 2)
mydata <- mydata - c(centerx, centery)
if (parallel) {
no_cores <- parallel::detectCores() - 1
cl <- parallel::makeCluster(no_cores)
doParallel::registerDoParallel(cl)
ellipsePoints <- foreach::foreach (currAngle = seq(
from = 0, to = pi - 0.001, by = pi / nrSlices)) %dopar% {
mydataRot <- mydata
if (currAngle != 0) {
mydataRot$x <- mydata$x * cos(currAngle) - mydata$y * sin(currAngle)
mydataRot$y <- mydata$y * cos(currAngle) + mydata$x * sin(currAngle)
}
minRoty <- min(mydataRot$y)
linedata <- mydataRot[mydataRot$x > (-thickness * boxSize / 2) &
mydataRot$x < (thickness * boxSize / 2), ]
distancesVector <- sapply(linedata$y, '-', minRoty)
dataCenter <- mean(distancesVector)
data <- distancesVector - dataCenter
if (varEst == "value") {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = varEst,
conf.level = conf.level)
}
downPointRot <- c(0, dataCenter + minRoty - unname(estResult$radius))
upPointRot <- c(0, dataCenter + minRoty + unname(estResult$radius))
# rotate calculated points back and undo translate:
downPointX <- downPointRot[1] * cos(currAngle) +
downPointRot[2] * sin(currAngle) +
centerx + clusterStartPoint[1] * boxSize
downPointY <- downPointRot[2] * cos(currAngle) -
downPointRot[1] * sin(currAngle) +
centery + clusterStartPoint[2] * boxSize
upPointX <- upPointRot[1] * cos(currAngle) +
upPointRot[2] * sin(currAngle) +
centerx + clusterStartPoint[1] * boxSize
upPointY <- upPointRot[2] * cos(currAngle) -
upPointRot[1] * sin(currAngle) +
centery + clusterStartPoint[2] * boxSize
upPoint <- c(upPointX, upPointY) %/% boxSize
downPoint <- c(downPointX, downPointY) %/% boxSize
c(upPoint, downPoint)
}
# cobvert that list to a matrix:
ellipsePoints <- matrix(unlist(ellipsePoints), ncol = 2, byrow = TRUE)
} else {
for (currAngle in seq(0, pi-0.001, by = pi / nrSlices)) {
# rotate the points by currAngle
mydataRot <- mydata
if (currAngle != 0) {
mydataRot$x <- mydata$x * cos(currAngle) - mydata$y * sin(currAngle)
mydataRot$y <- mydata$y * cos(currAngle) + mydata$x * sin(currAngle)
}
minRoty <- min(mydataRot$y)
linedata <- mydataRot[mydataRot$x > (-thickness * boxSize / 2) &
mydataRot$x < (thickness * boxSize / 2), ]
distancesVector <- sapply(linedata$y, '-', minRoty)
dataCenter <- mean(distancesVector)
data <- distancesVector - dataCenter
if (varEst == "value") {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = varEst,
conf.level = conf.level)
}
downPointRot <- c(0, dataCenter + minRoty - unname(estResult$radius))
upPointRot <- c(0, dataCenter + minRoty + unname(estResult$radius))
# rotate calculated points back and undo translate:
downPointX <- downPointRot[1] * cos(currAngle) +
downPointRot[2] * sin(currAngle) +
centerx + clusterStartPoint[1] * boxSize
downPointY <- downPointRot[2] * cos(currAngle) -
downPointRot[1] * sin(currAngle) +
centery + clusterStartPoint[2] * boxSize
upPointX <- upPointRot[1] * cos(currAngle) +
upPointRot[2] * sin(currAngle) +
centerx + clusterStartPoint[1] * boxSize
upPointY <- upPointRot[2] * cos(currAngle) -
upPointRot[1] * sin(currAngle) +
centery + clusterStartPoint[2] * boxSize
upPoint <- c(upPointX, upPointY) %/% boxSize
downPoint <- c(downPointX, downPointY) %/% boxSize
ellipsePoints[(2 * pointNr - 1), ] <- upPoint
ellipsePoints[(2 * pointNr), ] <- downPoint
pointNr <- pointNr + 1
}
}
}
end.time <- Sys.time()
evalTime <- round(end.time - start.time, 4)
ellipsePointsClean <- ellipsePoints[rowSums(is.na(ellipsePoints)) != 2, ]
resultPoints <- paste(toString(t(ellipsePointsClean), sep = ", "))
if (representation == 'ellipse') {
ellipseDirect <- conicfit::EllipseDirectFit(ellipsePointsClean)
ellipseDirectG <- conicfit::AtoG(ellipseDirect)$ParG
ellipseArea <- round(ellipseDirectG[3] * ellipseDirectG[4] * pi, 2)
# center of the ellipse and two points at the axes cross
C = c(ellipseDirectG[1], ellipseDirectG[2])
A = c(C[1] - ellipseDirectG[3] * cos(ellipseDirectG[5]),
C[2] - ellipseDirectG[3] * sin(ellipseDirectG[5]))
B = c(C[1] + ellipseDirectG[4] * sin(ellipseDirectG[5]),
C[2] - ellipseDirectG[4] * cos(ellipseDirectG[5]))
percArea <- ellipseArea / (height * width ) * 100
output <- paste("Levels of grey: ", levelsOfGray,
", Box size: ", boxSize,
", Line thickness: ", thickness,
", Error distribution: ", error,
"<br><strong>Area</strong>: ", ellipseArea, " pixels (",
round(percArea, 2), "% of the image area)<br>",
"Ellipse center (", round(ellipseDirectG[1]), ", ",
round(ellipseDirectG[2]), "), ",
"semiaxes: (", round(ellipseDirectG[3], 2), ", ",
round(ellipseDirectG[4], 2) , "), ",
"angle: ", round(ellipseDirectG[5], 2),
"<br>A = (", round(A[1]), ", ", round(A[2]), "), ",
"B = (", round(B[1]), ", ", round(B[2]), ")",
"<br>Parallel: ", parallel, " Slicing: ", slicing, " Time: ", evalTime,
sep = "")
} else if (representation == 'circle') {
circle <- conicfit::CircleFitByPratt(ellipsePointsClean)
circleArea <- circle[3]**2 * pi;
percArea <- circleArea / (height * width ) * 100
if (verbose) {
output <- paste("Levels of grey: ", levelsOfGray,
", Box size: ", boxSize,
", Line thickness: ", thickness,
", Error distribution: ", error,
"<br><strong>Area</strong>: ", round(circleArea, 2), " pixels (",
round(percArea, 2), "% of the image area)<br>",
"Circle center (", round(circle[1]), ", ", round(circle[2]), "), ",
"radius: ", round(circle[3]),
"<br>Parallel: ", parallel, " Slicing: ", slicing, " Time: ", evalTime,
sep = "")
} else {
output <- paste(round(circle[1]), " ", round(circle[2]), " ",
round(circle[1]+circle[3]), " ", round(circle[2]), " ",
round(circle[1]), " ", round(circle[2]+circle[3]), sep = " ")
}
# show circle as ellipse for drawing purposes
ellipseArea <- round(circleArea, 2)
ellipseDirectG <- c(circle[1], circle[2], circle[3], circle[3], 0)
}
# return(c(output, ellipseArea, ellipseDirectG, substring(resultPoints, 3)))
return(c(output, ellipseArea, ellipseDirectG, resultPoints))
}
#' Performs area estimation of the numerically described object in plane.
#'
#' Use this function if you have a data set of uniformly distributed points on
#' an elliptical domain in the plane but captured with additive errors. The
#' estimation algorithm takes many horizontal and vertical, or star-shaped
#' slices of the object. Length estimation procedure is conducted on each slice
#' and in that way the set of edge points is obtained. An ellipse or a circle
#' is fitted to these edge points by function
#' \code{\link[conicfit]{EllipseDirectFit}} or
#' \code{\link[conicfit]{CircleFitByPratt}} from the package \code{conicfit}
#' and its semi-axes and area are returned as a result. Function optionally
#' plots input points, calculated edge points and the resulting ellipse or
#' circle.
#'
#' @param data Two-column data matrix containing the points that describe
#' observed object. First column represents \emph{x} coordinate of the point,
#' while second column represents \emph{y} coordinate.
#' @param nrSlices Number of slices applied for plain data cutting. Defaults to
#' 10.
#' @param error A character string specifying the error distribution. Must be
#' one of "laplace", "gauss" or "student". Can be abbreviated.
#' @param var.est A character string specifying the method of error variance
#' estimation. Must be given if \code{var} is not given. Can be "MM" (Method
#' of Moments) or "ML" (Maximum Likelihood).
#' @param var Explicit error variance. Needs to be given if \code{var.est} is
#' not given.
#' @param plot Logical parameter (TRUE or FALSE) that determines whether to
#' plot given object, calculated edge points and the resulting ellipse.
#' Defaults to FALSE.
#' @param representation A character string specifying the shape of an observed
#' object. Can be "ellipse" or "circle". Can be abbreviated.
#' @param parallel Logical parameter (TRUE or FALSE) that determines whether to
#' perform estimation procedure in a parallel manner. Can shorten
#' estimation time if many border points need to be calculated. Defaults to
#' FALSE.
#' @param slicing A character string specifying the method of slicing. Can be
#' "hv" (horizontal and vertical slicing) or "star" (star-shaped slicing). Can
#' be abbreviated.
#'
#' @return List containing:
#' \itemize{
#' \item area: Estimated area of the object,
#' \item points: Set of calculated object's edge points,
#' \item semiaxes: Resulting ellipse's semi-axes or circle radius.
#' }
#'
#' @examples
#' # load a data set representing the ellipse with additive Gaussian error,
#' # run area estimation on it, and plot the results
#' inputfile <- system.file("extdata", "ellipse_3_4_0.1_gauss.txt", package = "LeArEst")
#' inputdata <- read.table(inputfile)
#' area <- areaest(inputdata, error = "gauss", var.est = "ML", plot = TRUE,
#' slicing = "hv", representation = "ellipse")
#'
#' # load a data set representing the ellipse with additive Laplacian error,
#' # run area estimation on it, and plot the results
#' inputfile <- system.file("extdata", "ellipse_3_4_0.1_laplace.txt", package = "LeArEst")
#' inputdata <- read.table(inputfile)
#' area <- areaest(inputdata, error = "laplace", var = 0.1, nrSlices = 5, plot = TRUE,
#' slicing = "star", representation = "ellipse")
#'
#' @importFrom utils read.csv
#' @importFrom graphics points
#' @importFrom foreach %dopar%
#' @importFrom foreach %do%
#'
#' @export
areaest <- function(data,
nrSlices = 10,
error = c("laplace", "gauss", "student"),
var.est = c("MM", "ML"),
var = NULL,
plot = FALSE,
parallel = FALSE,
slicing = c("hv", "star"),
representation = c("ellipse", "circle")) {
# convert matrix to dataframe
mydata <- data.frame(data)
colnames(mydata) <- c("x", "y")
minx <- min(mydata$x)
maxx <- max(mydata$x)
miny <- min(mydata$y)
maxy <- max(mydata$y)
ellipsePoints <- matrix(NA, nrow = nrSlices * 4, ncol = 2)
pointNr <- 1
biggerDim <- max((maxx - minx), (maxy - miny))
lineDistance <- biggerDim / nrSlices
currentX <- minx + lineDistance / 2
currentY <- miny + lineDistance / 2
conf.level <- 0.95
start.time <- Sys.time()
# horizontal chopping
if (slicing == "hv") {
if (parallel) {
no_cores <- parallel::detectCores() - 1
cl <- parallel::makeCluster(no_cores)
doParallel::registerDoParallel(cl)
ellipsePointsH <- foreach::foreach (currentX = seq(from = minx + lineDistance / 2,
to = maxx, by = lineDistance)) %dopar% {
linedata <- mydata[mydata$x > (currentX - lineDistance / 2) &
mydata$x < (currentX + lineDistance / 2), ]
distancesVector <- sapply(linedata$y, '-', miny)
# if (length(distancesVector) < 10) {
# c(NA, NA, NA, NA)
# }
dataCenter <- mean(distancesVector)
data <- distancesVector - dataCenter
if (!is.null(var)) {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = var.est,
conf.level = conf.level)
}
leftPoint <- c(currentX, dataCenter + miny - unname(estResult$radius))
rightPoint <- c(currentX, dataCenter + miny + unname(estResult$radius))
# rm(linedata, distancesVector, dataCenter, data, estResult)
c(leftPoint, rightPoint)
}
# razvrti tu listu u matricu:
ellipsePoints <- matrix(unlist(ellipsePointsH), ncol = 2, byrow = TRUE)
} else {
while (currentX < maxx) {
linedata <- mydata[mydata$x > (currentX - lineDistance / 2) &
mydata$x < (currentX + lineDistance / 2), ]
distancesVector <- sapply(linedata$y, '-', miny)
# if (length(distancesVector) < 10) {
# currentX <- currentX + lineDistance
# next
# }
dataCenter <- mean(distancesVector)
data <- distancesVector - dataCenter
if (!is.null(var)) {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = var.est,
conf.level = conf.level)
}
leftPoint <- c(currentX, dataCenter + miny - estResult$radius)
rightPoint <- c(currentX, dataCenter + miny + estResult$radius)
ellipsePoints[(2 * pointNr - 1), ] <- leftPoint
ellipsePoints[(2 * pointNr), ] <- rightPoint
currentX <- currentX + lineDistance
pointNr <- pointNr + 1
}
}
# vertical chopping
if (parallel) {
ellipsePointsV <- foreach::foreach(currentY = seq(from = miny + lineDistance / 2,
to = maxy, by = lineDistance)) %dopar% {
linedata <- mydata[mydata$y > (currentY - lineDistance / 2) &
mydata$y < (currentY + lineDistance / 2), ]
distancesVector <- sapply(linedata$x, '-', minx)
# if (length(distancesVector) < 10) {
# c(NA, NA, NA, NA)
# }
dataCenter <- mean(distancesVector)
data <- distancesVector - dataCenter
if (!is.null(var)) {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = var.est,
conf.level = conf.level)
}
leftPoint <- c(dataCenter + minx - unname(estResult$radius), currentY)
rightPoint <- c(dataCenter + minx + unname(estResult$radius), currentY)
# rm(linedata, distancesVector, dataCenter, data, estResult)
c(leftPoint, rightPoint)
}
ellipsePointsV <- matrix(unlist(ellipsePointsV), ncol = 2, byrow = TRUE)
ellipsePoints <- as.matrix(merge(ellipsePoints, ellipsePointsV,
all = TRUE))
} else {
while (currentY < maxy) {
linedata <- mydata[mydata$y > (currentY - lineDistance / 2) &
mydata$y < (currentY + lineDistance / 2), ]
distancesVector <- sapply(linedata$x, '-', minx)
# if (length(distancesVector) < 10) {
# currentY <- currentY + lineDistance
# next
# }
dataCenter <- mean(distancesVector)
data <- distancesVector - dataCenter
if (!is.null(var)) {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = var.est,
conf.level = conf.level)
}
leftPoint <- c(dataCenter + minx - estResult$radius, currentY)
rightPoint <- c(dataCenter + minx + estResult$radius, currentY)
ellipsePoints[(2 * pointNr - 1), ] <- leftPoint
ellipsePoints[(2 * pointNr), ] <- rightPoint
currentY <- currentY + lineDistance
pointNr <- pointNr + 1
}
}
} else if (slicing == "star") {
# common part:
# translate input points so the origin is data mean
centerx <- mean (mydata$x)
centery <- mean (mydata$y)
mydata$x <- mydata$x - centerx
mydata$y <- mydata$y - centery
lineDistance <- biggerDim / 10
if (parallel) {
no_cores <- parallel::detectCores() - 1
cl <- parallel::makeCluster(no_cores)
doParallel::registerDoParallel(cl)
ellipsePoints <- foreach::foreach (currAngle = seq(
from = 0, to = pi - 0.001, by = pi / nrSlices)) %dopar% {
# rotate the points by currAngle
mydataRot <- mydata
if (currAngle != 0) {
mydataRot$x <- mydata$x * cos(currAngle) - mydata$y * sin(currAngle)
mydataRot$y <- mydata$y * cos(currAngle) + mydata$x * sin(currAngle)
}
minRoty <- min(mydataRot$y)
linedata <- mydataRot[mydataRot$x > (-lineDistance / 2) &
mydataRot$x < (lineDistance / 2), ]
distancesVector <- sapply(linedata$y, '-', minRoty)
dataCenter <- mean(distancesVector)
data <- distancesVector - dataCenter
if (!is.null(var)) {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = var.est,
conf.level = conf.level)
}
downPointRot <- c(0, dataCenter + minRoty - unname(estResult$radius))
upPointRot <- c(0, dataCenter + minRoty + unname(estResult$radius))
# rotate calculated points back and undo translate:
downPointX <- downPointRot[1] * cos(currAngle) +
downPointRot[2] * sin(currAngle) + centerx
downPointY <- downPointRot[2] * cos(currAngle) -
downPointRot[1] * sin(currAngle) + centery
upPointX <- upPointRot[1] * cos(currAngle) +
upPointRot[2] * sin(currAngle) + centerx
upPointY <- upPointRot[2] * cos(currAngle) -
upPointRot[1] * sin(currAngle) + centery
upPoint <- c(upPointX, upPointY)
downPoint <- c(downPointX, downPointY)
c(upPoint, downPoint)
}
# razvrti tu listu u matricu:
ellipsePoints <- matrix(unlist(ellipsePoints), ncol = 2, byrow = TRUE)
} else {
for (currAngle in seq(0, pi - 0.001, by = pi / nrSlices)) {
# rotate the points by currAngle
mydataRot <- mydata
if (currAngle != 0) {
mydataRot$x <- mydata$x * cos(currAngle) - mydata$y * sin(currAngle)
mydataRot$y <- mydata$y * cos(currAngle) + mydata$x * sin(currAngle)
}
minRoty <- min(mydataRot$y)
linedata <- mydataRot[mydataRot$x > (-lineDistance / 2) &
mydataRot$x < (lineDistance / 2), ]
distancesVector <- sapply(linedata$y, '-', minRoty)
dataCenter <- mean(distancesVector)
data <- distancesVector - dataCenter
if (!is.null(var)) {
estResult <- lengthest(data, error = error, sd = sqrt(var),
conf.level = conf.level)
} else {
estResult <- lengthest(data, error = error, sd.est = var.est,
conf.level = conf.level)
}
downPointRot <- c(0, dataCenter + minRoty - unname(estResult$radius))
upPointRot <- c(0, dataCenter + minRoty + unname(estResult$radius))
# rotate calculated points back and undo translate:
downPointX <- downPointRot[1] * cos(currAngle) +
downPointRot[2] * sin(currAngle) + centerx
downPointY <- downPointRot[2] * cos(currAngle) -
downPointRot[1] * sin(currAngle) + centery
upPointX <- upPointRot[1] * cos(currAngle) +
upPointRot[2] * sin(currAngle) + centerx
upPointY <- upPointRot[2] * cos(currAngle) -
upPointRot[1] * sin(currAngle) + centery
ellipsePoints[(2 * pointNr - 1), ] <- c(downPointX, downPointY)
ellipsePoints[(2 * pointNr), ] <- c(upPointX, upPointY)
pointNr <- pointNr + 1
}
}
# undo translate for original points
mydata$x <- mydata$x + centerx
mydata$y <- mydata$y + centery
}
end.time <- Sys.time()
time.taken <- end.time - start.time
cat ("Time for chopping:", time.taken)
ellipsePointsClean <- ellipsePoints[rowSums(is.na(ellipsePoints)) != 2, ]
if (representation == "ellipse") {
ellipseDirect <- conicfit::EllipseDirectFit(ellipsePointsClean)
ellipseDirectG <- conicfit::AtoG(ellipseDirect)$ParG
ellipseArea <- ellipseDirectG[3] * ellipseDirectG[4] * pi
} else if (representation == "circle") {
circle <- conicfit::CircleFitByPratt(ellipsePointsClean)
circleArea <- circle[3]**2 * pi;
# show circle as ellipse for drawing purposes
ellipseArea <- circleArea
ellipseDirectG <- c(circle[1], circle[2], circle[3], circle[3], 0)
}
if (plot) {
xLim <- range(min(minx, ellipsePointsClean[, 1]),
max(maxx, ellipsePointsClean[, 1]))
yLim <- range(min(miny, ellipsePointsClean[, 2]),
max(maxy, ellipsePointsClean[, 2]))
xyDirect<-conicfit::calculateEllipse(ellipseDirectG[1], ellipseDirectG[2],
ellipseDirectG[3], ellipseDirectG[4], 180 / pi * ellipseDirectG[5])
plot(xyDirect[, 1], xyDirect[, 2], type = 'l', col = 'cyan',
xlab = "x", ylab = "y", lwd = 3, xlim = xLim, ylim = yLim)
points(mydata)
points(ellipsePointsClean, col = "red", lwd = 7)
}
return (list(
"area" = round(ellipseArea, 2),
"points" = ellipsePointsClean,
"semiaxes" = c(ellipseDirectG[3], ellipseDirectG[4])
))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.