From ca0fa3c4f1bf3eeee520ecf033f0c33f63223fa3 Mon Sep 17 00:00:00 2001 From: Aymeric Stamm Date: Thu, 8 Feb 2024 10:46:05 +0100 Subject: [PATCH 1/3] Fix bugs in find_rotation_seed_coord. --- R/curve_functions.R | 148 +++++++++++++++++++++++--------------------- 1 file changed, 78 insertions(+), 70 deletions(-) diff --git a/R/curve_functions.R b/R/curve_functions.R index 2582bf0..0fbac59 100644 --- a/R/curve_functions.R +++ b/R/curve_functions.R @@ -155,84 +155,92 @@ shift_f <- function(f, tau){ return(fn) } - -find_rotation_seed_coord <- function(beta1, beta2, mode="O", rotation=TRUE, - scale=TRUE){ - n = nrow(beta1) - T1 = ncol(beta1) - q1 = curve_to_q(beta1, scale)$q - - scl = 4 - minE = 1000 - if (mode=="C"){ - end_idx = floor(T1/scl) +find_rotation_seed_coord <- function(beta1, beta2, + mode = "O", + rotation = TRUE, + scale = TRUE) { + n <- nrow(beta1) + T1 <- ncol(beta1) + out <- curve_to_q(beta1, scale = scale) + q1 <- out$q + len1 <- out$len + + scl <- 4 + minE <- 1000 + if (mode == "C") + end_idx <- floor(T1 / scl) + else + end_idx <- 0 + + for (ctr in 0:end_idx) { + if (mode == "C") { + if ((scl*ctr) <= end_idx) + beta2n <- shift_f(beta2, scl * ctr) + else + break + } else + beta2n <- beta2 + + if (rotation) { + out <- find_best_rotation(beta1, beta2n) + beta2n <- out$q2new + Rout <- out$R + } else + Rout <- diag(nrow(beta2n)) + + out <- curve_to_q(beta2n, scale = scale) + q2n <- out$q + len2 <- out$len + + if (norm(q1 - q2n, 'F') > 0.0001) { + q1i <- q1 / sqrt(innerprod_q2(q1, q1)) + q2ni <- q2n / sqrt(innerprod_q2(q2n, q2n)) + dim(q1i) <- c(T1 * n) + dim(q2ni) <- c(T1 * n) + gam0 <- .Call('DPQ', PACKAGE = 'fdasrvf', q1i, q2ni, n, T1, 0, 1, 0, rep(0,T1)) + gamI <- invertGamma(gam0) + gam <- (gamI - gamI[1]) / (gamI[length(gamI)] - gamI[1]) + beta2new <- group_action_by_gamma_coord(beta2n, gam) + out <- curve_to_q(beta2new, scale = scale) + q2new <- out$q + len2 <- out$len + if (mode == "C") + q2new <- project_curve(q2new) } else { - end_idx = 0 + q2new <- q2n + beta2new <- beta2n + gam <- seq(0, 1, length.out = T1) } - for (ctr in 0:end_idx){ - if (mode=="C"){ - if ((scl*ctr) <= end_idx){ - beta2n = shift_f(beta2, scl*ctr) - } else { - break - } - } else { - beta2n = beta2 - } - - if (rotation){ - out = find_best_rotation(beta1, beta2n) - beta2n = out$q2new - q2n = curve_to_q(beta2n, scale)$q - Rout = out$R - } else { - q2n = curve_to_q(beta2n, scale)$q - Rout = diag(nrow(beta2n)) - } - + dist <- innerprod_q2(q1, q2new) + if (dist < -1) dist <- -1 + if (dist > 1) dist <- 1 - if (norm(q1-q2n,'F') > 0.0001){ - q1i = q1/sqrt(innerprod_q2(q1, q1)) - q2ni = q2n/sqrt(innerprod_q2(q2n, q2n)) - dim(q1i) = c(T1*n) - dim(q2ni) = c(T1*n) - gam0 = .Call('DPQ', PACKAGE = 'fdasrvf', q1i, q2ni, n, T1, 0, 1, 0, rep(0,T1)) - gamI = invertGamma(gam0) - gam = (gamI-gamI[1])/(gamI[length(gamI)]-gamI[1]) - beta2n = q_to_curve(q2n) - beta2new = group_action_by_gamma_coord(beta2n, gam) - q2new = curve_to_q(beta2new)$q - if (mode=="C"){ - q2new = project_curve(q2new) - } - } else{ - q2new = q2n - beta2new = beta2n - gam = seq(0,1,length.out=T1) - } - dist = innerprod_q2(q1,q2new) - if (dist < -1){ - dist = -1 - } - if (dist > 1){ - dist = 1 - } - Ec = acos(dist) - if (Ec < minE){ - Rbest = Rout - beta2best = beta2new - q2best = q2new - gambest = gam - minE = Ec - tau = scl*ctr - } + Ec <- acos(dist) + if (Ec < minE) { + Rbest <- Rout + beta2best <- beta2new + q2best <- q2new + gambest <- gam + minE <- Ec + tau <- scl * ctr } + } - return(list(beta2best=beta2best,q2best=q2best,Rbest=Rbest,gambest=gambest,tau=tau)) + ratio <- 1 + if (scale) + ratio <- len1 / len2 + + list( + beta2best = beta2best * ratio, + q2best = q2best, + Rbest = Rbest, + gambest = gambest, + tau = tau, + scale = ratio + ) } - find_rotation_seed_unqiue <- function(q1, q2, mode="O", rotation=TRUE, scale=TRUE, lam=0.0){ n1 = nrow(q1) From 9bbd9e1df68951c2be8b33f9744be675d599d683 Mon Sep 17 00:00:00 2001 From: Aymeric Stamm Date: Thu, 8 Feb 2024 10:47:07 +0100 Subject: [PATCH 2/3] Use correctly transformed beta2best from output of find_rotation_seed_coord() instead of duplicated the rotation and curve_to_q() operations. --- R/curve_pair_align.R | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/R/curve_pair_align.R b/R/curve_pair_align.R index 9b3d5c0..aba3541 100644 --- a/R/curve_pair_align.R +++ b/R/curve_pair_align.R @@ -37,9 +37,7 @@ curve_pair_align <- function(beta1, beta2, mode="O", rotation=TRUE, scale=TRUE){ out = find_rotation_seed_coord(beta1, beta2, mode, rotation, scale) gam = out$gambest q2n = out$q2best - beta2n = out$Rbest %*% shift_f(beta2, out$tau) - beta2n = group_action_by_gamma_coord(beta2n, gam) - q2n = curve_to_q(beta2n, scale)$q + beta2n = out$beta2best return(list(beta2n=beta2n, q2n=q2n, gam=gam, q1=q1, beta1=beta1, beta2=beta2, R=out$Rbest, tau=out$tau)) From a819c67d9ae7e7a2b8aeeb063827f56aa022bc57 Mon Sep 17 00:00:00 2001 From: Aymeric Stamm Date: Thu, 8 Feb 2024 10:48:19 +0100 Subject: [PATCH 3/3] Some reformatting, improved documentation and more output information for calc_shape_dist(). It should now do the same thing as in curve_pair_align(). --- R/calc_shape_dist.R | 147 ++++++++++++++++++++++++++--------------- man/calc_shape_dist.Rd | 58 +++++++++++----- 2 files changed, 133 insertions(+), 72 deletions(-) diff --git a/R/calc_shape_dist.R b/R/calc_shape_dist.R index 862d990..7d30575 100644 --- a/R/calc_shape_dist.R +++ b/R/calc_shape_dist.R @@ -1,73 +1,112 @@ #' Elastic Shape Distance #' -#' Calculate elastic shape distance between two curves beta1 and beta2. If the -#' curves beta1 and beta2 are describing multidimensional functional data, then -#' `rotation == FALSE` and `mode == 'O'` +#' Calculates the elastic shape distance between two curves `beta1` and `beta2`. +#' If the input curves are describing multidimensional functional data, then +#' most of the time the user should set `rotation == FALSE`, `scale == FALSE` +#' and `mode == "O"`. +#' +#' @param beta1 A numeric matrix of shape \eqn{L \times M} specifying an +#' \eqn{L}-dimensional curve evaluated on \eqn{M} sample points. +#' @param beta2 A numeric matrix of shape \eqn{L \times M} specifying an +#' \eqn{L}-dimensional curve evaluated on \eqn{M} sample points. This curve +#' will be aligned to `beta1`. +#' @param mode A character string specifying whether the input curves should be +#' considered open (`mode == "O"`) or closed (`mode == "C"`). Defaults to +#' `"O"`. +#' @param rotation A boolean specifying whether the distance should be made +#' invariant by rotation. Defaults to `TRUE`. +#' @param scale A boolean specifying whether the distance should be made +#' invariant by scaling. This is effectively achieved by making SRVFs having +#' unit length and switching to an appropriate metric on the sphere between +#' normalized SRVFs. Defaults to `TRUE`. +#' @param include.length A boolean specifying whether to include information +#' about the actual length of the SRVFs in the metric when using normalized +#' SRVFs to achieve scale invariance. This only applies if `scale == TRUE`. +#' Defaults to `FALSE`. +#' +#' @return A list with the following components: +#' - `d`: the amplitude (geodesic) distance; +#' - `dx`: the phase distance; +#' - `q1`: the SRVF of `beta1`; +#' - `q2n`: the SRVF of `beta2` after alignment and possible optimal rotation +#' and scaling; +#' - `beta1`: the input curve `beta1`; +#' - `beta2n`: the input curve `beta2` after alignment and possible optimal +#' rotation and scaling. +#' - `R`: the optimal rotation matrix that has been applied to the second curve; +#' - `betascale`: the optimal scaling factor that has been applied to the second +#' curve; +#' - `gam`: the optimal warping function that has been applied to the second +#' curve. #' -#' @param beta1 curve1, provided as a matrix of sizes \eqn{n\times T} for -#' \eqn{n}-dimensional curve on \eqn{T} sample points -#' @param beta2 curve 2, provided as a matrix of sizes \eqn{n\times T} for -#' \eqn{n}-dimensional curve on \eqn{T} sample points -#' @param mode Open (`"O"`) or Closed (`"C"`) curves -#' @param rotation Include rotation (default = `TRUE`) -#' @param scale scale curves to unit length (default = `TRUE`) -#' @param include.length include length in distance calculation (default = `FALSE`) -#' this only applies if `scale=TRUE` -#' @return Returns a list containing \item{d}{geodesic distance} -#' \item{dx}{phase distance} -#' \item{q1}{srvf of curve 1} -#' \item{q2n}{srvf of aligned curve 2} #' @keywords distances +#' #' @references Srivastava, A., Klassen, E., Joshi, S., Jermyn, I., (2011). Shape -#' analysis of elastic curves in euclidean spaces. Pattern Analysis and Machine -#' Intelligence, IEEE Transactions on 33 (7), 1415-1428. +#' analysis of elastic curves in euclidean spaces. Pattern Analysis and +#' Machine Intelligence, IEEE Transactions on 33 (7), 1415-1428. #' @references Kurtek, S., Srivastava, A., Klassen, E., and Ding, Z. (2012), -#' “Statistical Modeling of Curves Using Shapes and Related Features,” Journal -#' of the American Statistical Association, 107, 1152–1165. +#' “Statistical Modeling of Curves Using Shapes and Related Features,” Journal +#' of the American Statistical Association, 107, 1152–1165. +#' #' @export #' @examples #' out <- calc_shape_dist(beta[, , 1, 1], beta[, , 1, 4]) -calc_shape_dist <- function(beta1, beta2, mode="O", rotation=TRUE, - scale=TRUE, include.length=FALSE){ - T1 = ncol(beta1) - centroid1 = calculatecentroid(beta1) - dim(centroid1) = c(length(centroid1),1) - beta1 = beta1 - repmat(centroid1,1,T1) - out1 = curve_to_q(beta1, scale) - q1 = out1$q - lenq1 = out1$lenq - lenq2 = curve_to_q(beta2, scale)$lenq +calc_shape_dist <- function(beta1, beta2, + mode = "O", + rotation = TRUE, + scale = TRUE, + include.length = FALSE) { + T1 <- ncol(beta1) + centroid1 <- calculatecentroid(beta1) + dim(centroid1) <- c(length(centroid1), 1) + beta1 <- beta1 - repmat(centroid1, 1, T1) + out1 <- curve_to_q(beta1, scale = scale) + q1 <- out1$q + lenq1 <- out1$lenq + lenq2 <- curve_to_q(beta2, scale = scale)$lenq - centroid1 = calculatecentroid(beta2) - dim(centroid1) = c(length(centroid1),1) - beta2 = beta2 - repmat(centroid1,1,T1) + centroid1 <- calculatecentroid(beta2) + dim(centroid1) <- c(length(centroid1), 1) + beta2 <- beta2 - repmat(centroid1, 1, T1) - out = find_rotation_seed_coord(beta1, beta2, mode, rotation, scale) - q1dotq2 = innerprod_q2(q1, out$q2best) + out <- find_rotation_seed_coord( + beta1, beta2, + mode = mode, + rotation = rotation, + scale = scale + ) - if (q1dotq2 > 1){ - q1dotq2 = 1 - } else if(q1dotq2 < -1){ - q1dotq2 = -1 - } + q1dotq2 <- innerprod_q2(q1, out$q2best) + + if (q1dotq2 > 1) q1dotq2 <- 1 + if (q1dotq2 < -1) q1dotq2 <- -1 - if (scale){ - if (include.length){ - d = sqrt(acos(q1dotq2)^2+log(lenq1/lenq2)^2) - } else{ - d = acos(q1dotq2) - } + if (scale) { + if (include.length) + d <- sqrt(acos(q1dotq2) ^ 2 + log(lenq1 / lenq2) ^ 2) + else + d <- acos(q1dotq2) } else { - v = q1-out$q2best - d = sqrt(innerprod_q2(v, v)) + v <- q1 - out$q2best + d <- sqrt(innerprod_q2(v, v)) } - gam = out$gambest - time1 <- seq(0,1,length.out=T1) + gam <- out$gambest + time1 <- seq(0, 1, length.out = T1) binsize <- mean(diff(time1)) - psi <- sqrt(gradient(gam,binsize)) - v <- inv_exp_map(rep(1,length(gam)), psi) - dx <- sqrt(trapz(time1, v^2)) + psi <- sqrt(gradient(gam, binsize)) + v <- inv_exp_map(rep(1, length(gam)), psi) + dx <- sqrt(trapz(time1, v ^ 2)) - return(list(d=d,dx=dx,q1=q1,q2n=out$q2best)) + list( + d = d, + dx = dx, + q1 = q1, + q2n = out$q2best, + beta1 = beta1, + beta2n = out$beta2best, + R = out$Rbest, + betascale = out$scale, + gam = out$gambest + ) } diff --git a/man/calc_shape_dist.Rd b/man/calc_shape_dist.Rd index 05c8fbe..86c4946 100644 --- a/man/calc_shape_dist.Rd +++ b/man/calc_shape_dist.Rd @@ -14,39 +14,61 @@ calc_shape_dist( ) } \arguments{ -\item{beta1}{curve1, provided as a matrix of sizes \eqn{n\times T} for -\eqn{n}-dimensional curve on \eqn{T} sample points} +\item{beta1}{A numeric matrix of shape \eqn{L \times M} specifying an +\eqn{L}-dimensional curve evaluated on \eqn{M} sample points.} -\item{beta2}{curve 2, provided as a matrix of sizes \eqn{n\times T} for -\eqn{n}-dimensional curve on \eqn{T} sample points} +\item{beta2}{A numeric matrix of shape \eqn{L \times M} specifying an +\eqn{L}-dimensional curve evaluated on \eqn{M} sample points. This curve +will be aligned to \code{beta1}.} -\item{mode}{Open (\code{"O"}) or Closed (\code{"C"}) curves} +\item{mode}{A character string specifying whether the input curves should be +considered open (\code{mode == "O"}) or closed (\code{mode == "C"}). Defaults to +\code{"O"}.} -\item{rotation}{Include rotation (default = \code{TRUE})} +\item{rotation}{A boolean specifying whether the distance should be made +invariant by rotation. Defaults to \code{TRUE}.} -\item{scale}{scale curves to unit length (default = \code{TRUE})} +\item{scale}{A boolean specifying whether the distance should be made +invariant by scaling. This is effectively achieved by making SRVFs having +unit length and switching to an appropriate metric on the sphere between +normalized SRVFs. Defaults to \code{TRUE}.} -\item{include.length}{include length in distance calculation (default = \code{FALSE}) -this only applies if \code{scale=TRUE}} +\item{include.length}{A boolean specifying whether to include information +about the actual length of the SRVFs in the metric when using normalized +SRVFs to achieve scale invariance. This only applies if \code{scale == TRUE}. +Defaults to \code{FALSE}.} } \value{ -Returns a list containing \item{d}{geodesic distance} -\item{dx}{phase distance} -\item{q1}{srvf of curve 1} -\item{q2n}{srvf of aligned curve 2} +A list with the following components: +\itemize{ +\item \code{d}: the amplitude (geodesic) distance; +\item \code{dx}: the phase distance; +\item \code{q1}: the SRVF of \code{beta1}; +\item \code{q2n}: the SRVF of \code{beta2} after alignment and possible optimal rotation +and scaling; +\item \code{beta1}: the input curve \code{beta1}; +\item \code{beta2n}: the input curve \code{beta2} after alignment and possible optimal +rotation and scaling. +\item \code{R}: the optimal rotation matrix that has been applied to the second curve; +\item \code{betascale}: the optimal scaling factor that has been applied to the second +curve; +\item \code{gam}: the optimal warping function that has been applied to the second +curve. +} } \description{ -Calculate elastic shape distance between two curves beta1 and beta2. If the -curves beta1 and beta2 are describing multidimensional functional data, then -\code{rotation == FALSE} and \code{mode == 'O'} +Calculates the elastic shape distance between two curves \code{beta1} and \code{beta2}. +If the input curves are describing multidimensional functional data, then +most of the time the user should set \code{rotation == FALSE}, \code{scale == FALSE} +and \code{mode == "O"}. } \examples{ out <- calc_shape_dist(beta[, , 1, 1], beta[, , 1, 4]) } \references{ Srivastava, A., Klassen, E., Joshi, S., Jermyn, I., (2011). Shape -analysis of elastic curves in euclidean spaces. Pattern Analysis and Machine -Intelligence, IEEE Transactions on 33 (7), 1415-1428. +analysis of elastic curves in euclidean spaces. Pattern Analysis and +Machine Intelligence, IEEE Transactions on 33 (7), 1415-1428. Kurtek, S., Srivastava, A., Klassen, E., and Ding, Z. (2012), “Statistical Modeling of Curves Using Shapes and Related Features,” Journal