Nothing
##### Core functions for polynomials approaches on open outlines
#' Calculate orthogonal polynomial fits on open outlines
#'
#' Calculates orthogonal polynomial coefficients,
#' through a linear model fit (see \link{lm}), from a matrix of (x; y) coordinates or
#' a \link{Opn} object
#'
#' @param x a matrix (or a list) of (x; y) coordinates
#' @param degree polynomial degree for the fit (the Intercept is also returned)
#' @param baseline1 numeric the \eqn{(x; y)} coordinates of the first baseline
#' by default \eqn{(x= -0.5; y=0)}
#' @param baseline2 numeric the \eqn{(x; y)} coordinates of the second baseline
#' by default \eqn{(x= 0.5; y=0)}
#' @param nb.pts number of points to sample and on which to calculate polynomials
#' @param ... useless here
#' @return a list with components when applied on a single shape:
#' \itemize{
#' \item \code{coeff} the coefficients (including the intercept)
#' \item \code{ortho} whether orthogonal or natural polynomials were fitted
#' \item \code{degree} degree of the fit (could be retrieved through \code{coeff} though)
#' \item \code{baseline1} the first baseline point (so far the first point)
#' \item \code{baseline2} the second baseline point (so far the last point)
#' \item \code{r2} the r2 from the fit
#' \item \code{mod} the raw lm model
#' }
#' otherwise an \link{OpnCoe} object.
#' @note Orthogonal polynomials are sometimes called Legendre's polynomials. They are
#' preferred over natural polynomials since adding a degree do not change lower orders coefficients.
#' @family polynomials
#' @examples
#' data(olea)
#' o <- olea[1]
#' op <- opoly(o, degree=4)
#' op
#' # shape reconstruction
#' opi <- opoly_i(op)
#' coo_plot(o)
#' coo_draw(opi)
#' lines(opi, col='red')
#' # R2 for degree 1 to 10
#' r <- numeric()
#' for (i in 1:10) { r[i] <- opoly(o, degree=i)$r2 }
#' plot(2:10, r[2:10], type='b', pch=20, col='red', main='R2 / degree')
#' @rdname opoly
#' @export
opoly <- function(x,
...){UseMethod("opoly")}
#' @rdname opoly
#' @export
opoly.default <- function(x,
degree,
...) {
coo <- x
coo <- coo_check(coo)
if (missing(degree)) {
degree <- 5
message("'degree' not provided and set to ", degree)
}
x <- poly(coo[, 1], degree = degree, raw = FALSE)
mod <- lm(coo[, 2] ~ x)
r2 <- summary(mod)$r.squared
return(list(coeff = mod$coefficients, ortho = TRUE, degree = degree,
baseline1 = coo[1, ], baseline2 = coo[nrow(coo), ], r2 = r2,
mod = mod))
}
#' @rdname opoly
#' @export
opoly.Opn <- function(x,
degree,
baseline1 = c(-0.5, 0),
baseline2 = c(0.5, 0),
nb.pts = 120,
...) {
Opn <- x
# verify
Opn %<>% verify()
# we check a bit
min.pts <- min(sapply(Opn$coo, nrow))
if (nb.pts > min.pts) {
if (missing(nb.pts)) {
nb.pts <- min.pts
message("'nb.pts' missing and set to ", nb.pts)
} else {
nb.pts <- min.pts
message("at least one outline has less coordinates than 'nb.pts' ", nb.pts, "points")
}
}
if (missing(degree)) {
degree <- 5
message("'degree' missing and set to ", degree)
}
# we normalize
Opn <- coo_sample(Opn, nb.pts)
coo <- Opn$coo
coo <- lapply(coo, coo_baseline, ldk1 = 1, ldk2 = nb.pts,
t1 = baseline1, t2 = baseline2)
# we prepare the coe matrix
rn <- names(coo)
cn <- paste0("x", 1:degree)
cn <- c("Intercept", cn)
coe <- matrix(NA, nrow = length(Opn), ncol = degree + 1,
dimnames = list(rn, cn))
r2 <- numeric(length(Opn))
mod <- list()
# the loop
for (i in seq(along = coo)) {
mod <- opoly(coo[[i]], degree = degree)
# mod[[i]] <- pol
coe[i, ] <- mod$coeff
r2[i] <- mod$r2
}
# mod$coefficients <- rep(NA, length(mod$coefficients))
method <- "opoly"
res <- OpnCoe(coe = coe, fac = Opn$fac, method = method,
baseline1 = baseline1, baseline2 = baseline2, r2 = r2,
mod = mod)
res$cuts <- ncol(res$coe)
return(res)
}
#' @rdname opoly
#' @export
opoly.list <- function(x, ...){
lapply(x, opoly, ...)
}
#' Calculate natural polynomial fits on open outlines
#'
#' Calculates natural polynomial coefficients,
#' through a linear model fit (see \link{lm}), from a matrix of (x; y) coordinates
#' or an \link{Opn} object
#'
#' @param x a matrix (or a list) of (x; y) coordinates or an \link{Opn} object
#' @param degree polynomial degree for the fit (the Intercept is also returned)
#' @param baseline1 numeric the \eqn{(x; y)} coordinates of the first baseline
#' by default \eqn{(x= -0.5; y=0)}
#' @param baseline2 numeric the \eqn{(x; y)} coordinates of the second baseline
#' by default \eqn{(x= 0.5; y=0)}
#' @param nb.pts number of points to sample and on which to calculate polynomials
#' @param ... useless here
#' @return when applied on a single shape, a list with components:
#' \itemize{
#' \item \code{coeff} the coefficients (including the intercept)
#' \item \code{ortho} whether orthogonal or natural polynomials were fitted
#' \item \code{degree} degree of the fit (could be retrieved through \code{coeff} though)
#' \item \code{baseline1} the first baseline point (so far the first point)
#' \item \code{baseline2} the second baseline point (so far the last point)
#' \item \code{r2} the r2 from the fit
#' \item \code{mod} the raw lm model
#' }
#'
#' otherwise, an \link{OpnCoe} object.
#' @family polynomials
#' @examples
#' data(olea)
#' o <- olea[1]
#' op <- opoly(o, degree=4)
#' op
#' # shape reconstruction
#' opi <- opoly_i(op)
#' coo_plot(o)
#' coo_draw(opi, border="red")
#' # R2 for degree 1 to 10
#' r <- numeric()
#' for (i in 1:10) { r[i] <- npoly(o, degree=i)$r2 }
#' plot(2:10, r[2:10], type='b', pch=20, col='red', main='R2 / degree')
#' @rdname npoly
#' @export
npoly <- function(x, ...){UseMethod("npoly")}
#' @rdname npoly
#' @export
npoly.default <- function(x, degree, ...) {
coo <- x
coo <- coo_check(coo)
if (missing(degree)) {
degree <- 5
message("'degree' not provided and set to: ", degree)
}
x <- poly(coo[, 1], degree = degree, raw = TRUE)
mod <- lm(coo[, 2] ~ x)
r2 <- summary(mod)$r.squared
return(list(coeff = mod$coefficients, ortho = FALSE, degree = degree,
baseline1 = coo[1, ], baseline2 = coo[nrow(coo), ], r2 = r2,
mod = mod))
}
#' @rdname npoly
#' @export
npoly.Opn <- function(x,
degree,
baseline1 = c(-0.5, 0),
baseline2 = c(0.5, 0),
nb.pts = 120, ...) {
Opn <- x
# verify
Opn %<>% verify()
# we check a bit
min.pts <- min(sapply(Opn$coo, nrow))
if (nb.pts > min.pts) {
if (missing(nb.pts)) {
nb.pts <- min.pts
message("'nb.pts' missing and set to: ", nb.pts)
} else {
nb.pts <- min.pts
message("at least one outline has less coordinates than 'nb.pts': ", nb.pts)
}
}
if (missing(degree)) {
degree <- 5
message("'degree' missing and set to: ", degree)
}
# we normalize
Opn <- coo_sample(Opn, nb.pts)
coo <- Opn$coo
coo <- lapply(coo, coo_baseline, ldk1 = 1, ldk2 = nb.pts,
t1 = baseline1, t2 = baseline2)
# we prepare the coe matrix
rn <- names(coo)
cn <- paste0("x", 1:degree)
cn <- c("Intercept", cn)
coe <- matrix(NA, nrow = length(Opn), ncol = degree + 1,
dimnames = list(rn, cn))
r2 <- numeric(length(Opn))
mod <- list()
# the loop
for (i in seq(along = coo)) {
mod <- npoly(coo[[i]], degree = degree)
# mod[[i]] <- pol
coe[i, ] <- mod$coeff
r2[i] <- mod$r2
}
# mod$coefficients <- rep(NA, length(mod$coefficients))
method <- "npoly"
res <- OpnCoe(coe = coe, fac = Opn$fac, method = method,
baseline1 = baseline1, baseline2 = baseline2, r2 = r2,
mod = mod)
res$cuts <- ncol(res$coe)
return(res)
}
#' @rdname npoly
#' @export
npoly.list <- function(x, ...){
lapply(x, npoly, ...)
}
#' Calculates shape from a polynomial model
#'
#' Returns a matrix of (x; y) coordinates when passed with a list obtained with
#' \link{opoly} or \link{npoly}.
#' @param pol a pol list such as created by \link{npoly} or \link{opoly}
#' @param nb.pts the number of points to predict. By default (and cannot be higher)
#' the number of points in the original shape.
#' @param reregister logical whether to reregister the shape with the original baseline.
#' @return a matrix of (x; y) coordinates.
#' @family polynomials
#' @examples
#' data(olea)
#' o <- olea[5]
#' coo_plot(o)
#' for (i in 2:7){
#' x <- opoly_i(opoly(o, i))
#' coo_draw(x, border=col_summer(7)[i], points=FALSE) }
#' @rdname poly_i
#' @export
opoly_i <- function(pol, nb.pts = 120, reregister = TRUE) {
.mprod <- function (m, s) {
res <- m
for (i in 1:ncol(m)) {
res[, i] <- m[, i] * s[i]
}
return(res)
}
x.new <- seq(pol$baseline1[1], pol$baseline2[1], length = nb.pts)
degree <- length(pol$coeff) - 1
x.poly <- poly(x.new, degree = degree)
y.pred <- predict(x.poly, x.new)
y.new <- pol$coeff[1] + apply(.mprod(m = y.pred, s = pol$coeff[-1]),
1, sum)
coo <- cbind(x.new, y.new)
if (reregister) {
coo <- coo_baseline(coo, 1, nrow(coo), t1 = pol$baseline1, t2 = pol$baseline2)
}
colnames(coo) <- c("x", "y")
return(coo)
}
#' @rdname poly_i
#' @export
npoly_i <- function(pol, nb.pts = 120, reregister = TRUE) {
x.new <- seq(pol$baseline1[1], pol$baseline2[1], length = nb.pts)
degree <- length(pol$coeff) - 1
y.pred <- numeric(nb.pts)
for (i in 1:degree) {
y.pred <- y.pred + (x.new^i * pol$coeff[i + 1])
}
y.new <- y.pred + pol$coeff[1]
coo <- cbind(x.new, y.new)
if (reregister) {
coo <- coo_baseline(coo, 1, nrow(coo), t1 = pol$baseline1, t2 = pol$baseline2)
}
colnames(coo) <- c("x", "y")
return(coo)
}
#opoly_shape #TODO
#npoly_shape #TODO
##### end Polynomials
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.