| Title: | Piecewise Geodesic Smoothing for Spherical Data |
|---|---|
| Description: | Fitting a smooth path to a given set of noisy spherical data observed at known time points. It implements a piecewise geodesic curve fitting method on the unit sphere based on a velocity-based penalization scheme. The proposed approach is implemented using the Riemannian block coordinate descent algorithm. To understand the method and algorithm, one can refer to Bak, K. Y., Shin, J. K., & Koo, J. Y. (2023) <doi:10.1080/02664763.2022.2054962> for the case of order 1. Additionally, this package includes various functions necessary for handling spherical data. |
| Authors: | Jae-Hwan Jhong [aut] (<https://orcid.org/0000-0003-2266-4986>, Chungbuk National University), Ja-Yong Koo [aut] (Korea University), Seyoung Lee [aut] (Sungshin Women's University), Kwan-Young Bak [aut, cre, cph] (<https://orcid.org/0000-0002-4541-160X>, Sungshin Women's University) |
| Maintainer: | Kwan-Young Bak <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.1.3 |
| Built: | 2024-12-07 10:52:19 UTC |
| Source: | CRAN |
A polar wander dataset presented in Kent and Irving (2010). The 17 Triassic/Jurassic cratonic poles from other major cratons are rotated into North American coordinates and combined with the 14 observations from North America. Our method is applied to these 31 observations ranging in age from 243 to 144 Ma (millions of years ago), which covers the late Triassic and Jurassic periods. The first column represents the time points, and the remaining two columns provides the observed spherical coordinates.
This function calculates the loss function based on the squared spherical distances between observed values and predicted values on the curve.
calculate_loss(y, gamma)calculate_loss(y, gamma)
y |
Matrix of observed values. |
gamma |
Matrix of predicted values. |
Loss value.
This function converts Cartesian coordinates to spherical coordinates.
cartesian_to_spherical(x, byrow = TRUE)cartesian_to_spherical(x, byrow = TRUE)
x |
A matrix where each row represents a point in Cartesian coordinates. |
byrow |
logical. If TRUE (the default) the matrix is filled by rows, otherwise the matrix is filled by columns. |
The Cartesian coordinates (x, y, z) are converted to spherical coordinates (theta, phi). Theta represents the inclination angle (0 to pi), and phi represents the azimuth angle (0 to 2*pi).
A matrix where each row represents a point in spherical coordinates.
#example1 cartesian_points1 <- matrix(c(1/sqrt(3), 1/sqrt(3), 1/sqrt(3),-1/sqrt(3), 1/sqrt(3), -1/sqrt(3)), ncol = 3, byrow = TRUE) cartesian_to_spherical(cartesian_points1) #example2 cartesian_points2 <- matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1),ncol = 3, byrow = TRUE) cartesian_to_spherical(cartesian_points2)#example1 cartesian_points1 <- matrix(c(1/sqrt(3), 1/sqrt(3), 1/sqrt(3),-1/sqrt(3), 1/sqrt(3), -1/sqrt(3)), ncol = 3, byrow = TRUE) cartesian_to_spherical(cartesian_points1) #example2 cartesian_points2 <- matrix(c(1, 0, 0, 0, 1, 0, 0, 0, 1),ncol = 3, byrow = TRUE) cartesian_to_spherical(cartesian_points2)
This function computes the cross product of two input vectors u and v.
cross(u, v, normalize = FALSE)cross(u, v, normalize = FALSE)
u |
Numeric vector. |
v |
Numeric vector. |
normalize |
logical. If TRUE, returns the normalized vector of the cross product result. |
Numeric vector representing the cross product of u and v.
cross(c(1,0,0), c(0,1,0))cross(c(1,0,0), c(0,1,0))
This function computes the dot product of two input vectors u and v.
dot(u, v)dot(u, v)
u |
Numeric vector. |
v |
Numeric vector. |
Numeric value representing the dot product of u and v.
dot(c(1,2,3), c(4,5,6))dot(c(1,2,3), c(4,5,6))
This function computes the equal-distance projection of a point p onto the xy plane.
edp(p)edp(p)
p |
Numeric vector representing a point in Cartesian coordinates. |
Numeric vector representing the equal-distance projection of p onto the xy plane.
This function computes the exponential map on the unit sphere given a base point x and a vector v.
exp_map(x, v)exp_map(x, v)
x |
Numeric vector representing the base point. |
v |
Numeric vector representing a point. |
Numeric vector representing the result of the exponential map.
exp_map(c(0,0,1), c(1,1,0))exp_map(c(0,0,1), c(1,1,0))
This function computes the value of the geodesic curve connecting two points p and q on the unit sphere at specified time points.
geodesic(t, p, q, a, b)geodesic(t, p, q, a, b)
t |
Numeric vector representing time points for the geodesic path. |
p |
Numeric vector representing the starting point on the sphere. |
q |
Numeric vector representing the ending point on the sphere. |
a |
Start time parameter. |
b |
End time parameter. |
Numeric matrix representing points along the geodesic path at specified time points.
geodesic(c(0.25, 0.5, 0.75), c(1,0,0), c(0,1,0), 0, 1)geodesic(c(0.25, 0.5, 0.75), c(1,0,0), c(0,1,0), 0, 1)
This function computes points along the geodesic connecting two points p and q on the unit sphere.
geodesic_lower(t, p, q, a, b)geodesic_lower(t, p, q, a, b)
t |
Time parameter for the geodesic path. |
p |
Numeric vector representing the starting point. |
q |
Numeric vector representing the ending point. |
a |
Start time parameter. |
b |
End time parameter. |
Numeric vector representing a point along the geodesic path.
geodesic_lower(0.5, c(1,0,0), c(0,1,0), 0, 1)geodesic_lower(0.5, c(1,0,0), c(0,1,0), 0, 1)
A tropical cyclone dataset provided by the Regional Specialized Meteorological Center (RSMC) Tokyo Typhoon Center. We select a cyclone called 'Goni' observed over time on August, 2015. The first column represents the time points, and the remaining two columns provides the observed spherical coordinates.
This generates a sequence of knots for a given set of time points based on the quantiles.
knots_quantile(x, dimension, tiny = 1e-05)knots_quantile(x, dimension, tiny = 1e-05)
x |
Numeric vector representing time points for the geodesic path. |
dimension |
Numeric vector the number of knots. |
tiny |
Numeric value representing a small constant that slightly expands the boundary. |
Numeric vector representing knots sequence in the time domain.
knots_quantile(seq(0, 1, length.out = 100), 10)knots_quantile(seq(0, 1, length.out = 100), 10)
This function computes the L2 norm (Euclidean norm) of the input vector u.
norm2(u)norm2(u)
u |
Numeric vector. |
Numeric value representing the L2 norm of u.
norm2(c(1,2,3))norm2(c(1,2,3))
This function normalizes the rows of the input matrix x by dividing each row by its L2 norm (Euclidean norm).
normalize(x)normalize(x)
x |
Numeric matrix. |
Numeric matrix with normalized rows.
normalize(matrix(c(1,2,3,4,5,6), nrow = 2, byrow = TRUE))normalize(matrix(c(1,2,3,4,5,6), nrow = 2, byrow = TRUE))
This function normalizes the the input vector v by dividing its L2 norm (Euclidean norm).
normalize_lower(v)normalize_lower(v)
v |
Numeric vector. |
Numeric vector with normalized.
normalize_lower(1:6)normalize_lower(1:6)
This function fits a penalized piecewise geodesic curve (linear spherical spline) to the given data.
penalized_linear_spherical_spline( t, y, initial_control_points = NULL, dimension, initial_knots, lambdas, step_size = 1, maxiter = 1000, epsilon_iter = 0.001, jump_eps = 1e-04, verbose = FALSE )penalized_linear_spherical_spline( t, y, initial_control_points = NULL, dimension, initial_knots, lambdas, step_size = 1, maxiter = 1000, epsilon_iter = 0.001, jump_eps = 1e-04, verbose = FALSE )
t |
A numeric vector representing the time or location. |
y |
A matrix where each row represents a data point on the sphere. |
initial_control_points |
An optional matrix specifying initial control points. Default is NULL. |
dimension |
An integer specifying the dimension of the spline. |
initial_knots |
An optional numeric vector specifying initial knots. Default is NULL. |
lambdas |
A numeric vector specifying the penalization parameters. |
step_size |
A numeric value specifying the step size for optimization. Default is 1. |
maxiter |
An integer specifying the maximum number of iterations. Default is 1000. |
epsilon_iter |
A numeric value specifying the convergence criterion for iterations. Default is 1e-03. |
jump_eps |
A numeric value specifying the threshold for pruning control points based on jump size. Default is 1e-04. |
verbose |
A logical value indicating whether to print progress information. Default is FALSE. |
The goal is to find the optimal piecewise geodesic curve for the given spherical data while controlling model complexity through penalty terms. This function computes the optimal control points and knots for the given data and returns the fitted result. Internally, coordinate-wise gradient descent is used to minimize the loss function, and a penalty term is added to control the complexity of the model. The BIC (Bayesian Information Criterion) value is calculated according to the model's complexity to provide information for model selection. The function constructs piecewise curves using the piecewise_geodesic function and employs penalty terms to control the complexity of the model by updating control points and knots. To see how to use the function in practical applications, refer to the README or https://github.com/kybak90/spheresmooth.
A list containing the fitted result for each complexity parameter and BIC values for model selection. One might choose the element that corresponds to the minimum BIC values as illustrated in the example.
apw_cartesian = spherical_to_cartesian(apw_spherical[, 2:3]) t = apw_spherical[, 1] dimension = 3 initial_knots = knots_quantile(t, dimension = dimension) lambda_seq = exp(seq(log(1e-04), log(1), length = 3)) fit = penalized_linear_spherical_spline(t = t, y = apw_cartesian, dimension = dimension, initial_knots = initial_knots, lambdas = lambda_seq) # choose a curve that minimizes the BIC best_index = which.min(fit$bic_list) best_index # obtained control points for the piecewise geodesic curve fit[[best_index]]$control_pointsapw_cartesian = spherical_to_cartesian(apw_spherical[, 2:3]) t = apw_spherical[, 1] dimension = 3 initial_knots = knots_quantile(t, dimension = dimension) lambda_seq = exp(seq(log(1e-04), log(1), length = 3)) fit = penalized_linear_spherical_spline(t = t, y = apw_cartesian, dimension = dimension, initial_knots = initial_knots, lambdas = lambda_seq) # choose a curve that minimizes the BIC best_index = which.min(fit$bic_list) best_index # obtained control points for the piecewise geodesic curve fit[[best_index]]$control_points
This function computes a piecewise geodesic path between control points.
piecewise_geodesic(t, control_points, knots)piecewise_geodesic(t, control_points, knots)
t |
A numeric vector representing the time or location. |
control_points |
A matrix of control points where each row represents a control point. |
knots |
A numeric vector of knot values. |
This function calculates the piecewise geodesic curve between control points based on the provided knots. The geodesic curve is computed segment by segment between adjacent control points. It interpolates the path between control points in a geodesic manner, ensuring the shortest path along the surface.
A matrix containing the piecewise geodesic path.
# `rgl` package and `sphereplot` pacakges are needed for the visualizaiton of the following example. # Define control points and knots library(rgl) library(sphereplot) control_points <- matrix(c(1, 0, 0, # Control point 1 1/sqrt(2), 1/sqrt(2), 0, # Control point 2 -1/sqrt(3), 1/sqrt(3), 1/sqrt(3), # Control point 3 0, 0, 1), # Control point 4 nrow = 4, byrow = TRUE) knots <- c(1, 2, 3, 3.5) # Knots indicating transitions # Example of generating piecewise geodesic curve t_example <- seq(0, 4, by = 0.01) gamma_example <- piecewise_geodesic(t_example, control_points, knots) # Plotting the piecewise geodesic curve rgl.sphgrid(deggap = 45, col.long = "skyblue", col.lat = "skyblue") spheres3d(x = 0, y = 0, z = 0, radius = 1, col = "grey", alpha = 0.05) pch3d(control_points, col = "blue", cex = 0.2, pch = 19) lines3d(gamma_example, col = "red", lty = 1, lwd = 2)# `rgl` package and `sphereplot` pacakges are needed for the visualizaiton of the following example. # Define control points and knots library(rgl) library(sphereplot) control_points <- matrix(c(1, 0, 0, # Control point 1 1/sqrt(2), 1/sqrt(2), 0, # Control point 2 -1/sqrt(3), 1/sqrt(3), 1/sqrt(3), # Control point 3 0, 0, 1), # Control point 4 nrow = 4, byrow = TRUE) knots <- c(1, 2, 3, 3.5) # Knots indicating transitions # Example of generating piecewise geodesic curve t_example <- seq(0, 4, by = 0.01) gamma_example <- piecewise_geodesic(t_example, control_points, knots) # Plotting the piecewise geodesic curve rgl.sphgrid(deggap = 45, col.long = "skyblue", col.lat = "skyblue") spheres3d(x = 0, y = 0, z = 0, radius = 1, col = "grey", alpha = 0.05) pch3d(control_points, col = "blue", cex = 0.2, pch = 19) lines3d(gamma_example, col = "red", lty = 1, lwd = 2)
Fitting a smooth path to a given set of noisy spherical data observed at known time points. It implements a piecewise geodesic curve fitting method on the unit sphere based on a velocity-based penalization scheme. The proposed approach is implemented using the Riemannian block coordinate descent algorithm. To understand the method and algorithm, one can refer to Bak, K. Y., Shin, J. K., & Koo, J. Y. (2023) <doi:10.1080/02664763.2022.2054962> for the case of order 1. Additionally, this package includes various functions necessary for handling spherical data.
A spheresmooth package
This function calculates the spherical distance between two vectors.
spherical_dist(x, y)spherical_dist(x, y)
x |
A numeric vector. |
y |
A numeric vector. |
The distance between vectors x and y.
x <- c(1, 0, 0) y <- c(0, 1, 0) spherical_dist(x, y)x <- c(1, 0, 0) y <- c(0, 1, 0) spherical_dist(x, y)
This function converts spherical coordinates (theta, phi) to Cartesian coordinates.
spherical_to_cartesian(theta_phi, byrow = TRUE)spherical_to_cartesian(theta_phi, byrow = TRUE)
theta_phi |
A matrix where each row contains the spherical coordinates (theta, phi) of a point. |
byrow |
logical. If TRUE (the default) the matrix is filled by rows, otherwise the matrix is filled by columns. |
A matrix where each row contains the Cartesian coordinates (x, y, z) of a point.
theta_phi <- matrix(c(pi/4, pi/3, pi/6, pi/4), ncol = 2, byrow = TRUE) spherical_to_cartesian(theta_phi)theta_phi <- matrix(c(pi/4, pi/3, pi/6, pi/4), ncol = 2, byrow = TRUE) spherical_to_cartesian(theta_phi)