R/spherepc.R

Defines functions LPG SPC.Hauberg SPC Dist.pt Proj.Hauberg PGA GenerateCircle PrincipalCircle IntrinsicMean ExtrinsicMean Cal.recon Kernel.Gaussian Kernel.indicator Kernel.quartic Rotate.inv Rotate Crossprod Expmap Logmap Trans.sph Trans.Euclid

Documented in Cal.recon Crossprod Dist.pt Expmap ExtrinsicMean GenerateCircle IntrinsicMean Kernel.Gaussian Kernel.indicator Kernel.quartic Logmap LPG PGA PrincipalCircle Proj.Hauberg Rotate Rotate.inv SPC SPC.Hauberg Trans.Euclid Trans.sph

############################################################################
### Spherical principal curves
############################################################################

#rm(list = ls())

### install rgl, sphereplot, and geosphere packages--------------------------

#install.packages("rgl")
#install.packages("sphereplot")
#install.packages("geosphere")

#require(rgl)
#require(sphereplot)
#require(geosphere)


##########################################################################
# Define funcions
##########################################################################

  Trans.Euclid <- function(vec){               # Input: 'vec' is 2-dim spherical coordinate (longitude, latitude). Output: 3-dim Euclidean coordinate.
    vec <- as.numeric(vec)
    if (length(vec) != 2) stop("length of vector is not 2")

    return( c(cospi(vec[1] / 180) * cospi(vec[2] / 180), cospi(vec[2] / 180) * sinpi(vec[1] / 180), sinpi(vec[2] / 180)) )
    }

  
  Trans.sph <- function(vec){                  # Input: 'vec' is 3-dim Euclidean coordinate. Output: 2-dim spherical coordiante (longitude, latitude).
    vec <- as.numeric(vec)
    if (sum(vec^2) == 0) stop("Input should not be (0, 0, 0)")
    if (length(vec) != 3) stop("length of vector is not 3")
    vec <- vec/(sum(vec^2))^(1/2)              # normalize so that 'vec' is in the unit sphere.
    if ((vec[2] >= 0) && (vec[1] >= 0)){
       if (vec[1]^2 + vec[2]^2 == 0){
          if (vec[3] > 0){
            return(c(0, 90))
          }else{
            return(c(0, -90))
          }
       }else{
         return (c(atan(vec[2]/ vec[1]), asin(vec[3])) * 180/pi)
       }
    }
    if ((vec[2] >= 0) && (vec[1] < 0)){
      return (c(pi + atan(vec[2]/vec[1]), asin(vec[3])) * 180/pi)
    }
    if ((vec[2] < 0) && (vec[1] >= 0)){
      return (c(atan(vec[2]/vec[1]), asin(vec[3])) * 180/pi)
    }
    if ((vec[2] < 0) && (vec[1] < 0)){
      return (c(atan(vec[2]/vec[1]) - pi, asin(vec[3])) * 180/pi)
    }
  }

  Logmap <- function(vec) {                                                                                 # Input: 3-dim Euclidean coordinate. Output: 2-dim Euclidean coordinate.
    vec <- as.numeric(vec)
    vec <- vec/(sum(vec^2))^(1/2)                                                                           # normalize so that 'vec' is in the unit sphere.
    if (length(vec) != 3) stop("length of vector is not 3")
    if (sum(vec == c(0, 0, -1)) == 3) stop("input is the cut locus of (0, 0, 1)")                           # Input should not be (0, 0, -1).
    z <- vec[3]
    x <- vec[1]
    y <- vec[2]
    result <- c(x * asin((1 - z^2)^(1/2))/(1 - z^2)^(1/2), y * asin((1 - z^2)^(1/2))/(1 - z^2)^(1/2))
    if (is.nan(result[1]) == TRUE){
        return (c(0, 0))
      }
        return (result)
    }

  Expmap <- function(vec) {                                                                                 # Input: 2 dim-Euclidean vector. Output: 3 dim-Euclidean vector.
    vec <- as.numeric(vec)
    if (length(vec) != 2) stop("number of vector component is not 2")

    if(vec[1]^2 + vec[2]^2 == 0){
      return(c(0, 0, 1))
    }else{
      norm <- (vec[1]^2 + vec[2]^2)^(1/2)
      return( c(vec[1] * sin(norm)/norm, vec[2] * sin(norm)/norm, cos(norm)) )
    }
  }

  Crossprod <- function(vec1, vec2) {                                                                      # Cross product of 3-dim two vectors.
    if(length(vec1) != 3 | length(vec2) != 3) stop('number of vector component is not 3')
    x <- vec1[2] * vec2[3] - vec1[3] * vec2[2]
    y <- vec1[3] * vec2[1] - vec1[1] * vec2[3]
    z <- vec1[1] * vec2[2] - vec1[2] * vec2[1]
    result <- c(x, y, z)
    return (result)
   }

  Rotate <- function(pt1, pt2){                                                                            # Input pt1, pt2: spatial locations (longitude, latitude). Output: 3 dim-Euclidean coordinate.
    pt1 <- as.numeric(pt1)                                                                                 # Rotation function: Rotate 'pt2' by the extent from which rotate 'pt1' to spherical coordinate (0, 90) (= (0, 0, 1) Euclidean coordinate.)
    pt2 <- as.numeric(pt2)
    a <- Trans.Euclid(pt1)
    v <- Trans.Euclid(pt2)
    axis <- Crossprod(a, c(0, 0, 1))                                                                       # rotation axis.
    if (sum(axis^2) == 0){
      return(Trans.Euclid(pt2))
    }
    axis <- axis/(sum(axis^2))^(1/2)                                                                       # normalized unit axis.
    theta <- abs(pi/2 - pt1[2] * pi/180)                                                                   # rotation angle.
    v.rot <- v * cos(theta) + Crossprod(axis, v) * sin(theta) + axis * (sum(axis * v)) * (1 - cos(theta))  # Rodrigues' rotation formula.
    return (v.rot)
  }
 
 Rotate.inv <- function(pt1, pt2){                                                                         # Input pt1, pt2: measurements of angular form of (longitude, latitude). Output: 3-dim Euclidean vector.
    pt1 <- as.numeric(pt1)                                                                                 # inverse Rotation function:  inverse rotation of Rotate function.
    pt2 <- as.numeric(pt2)
    a <- Trans.Euclid(pt1)
    v <- Trans.Euclid(pt2)
    axis <- Crossprod(c(0, 0, 1), a)
    if (sum(axis^2) == 0){
      return(Trans.Euclid(pt2))
    }
    axis <- axis/(sum(axis^2))^(1/2)                                                                       # Normalized unit axis.
    theta <- abs(pi/2 - pt1[2] * pi/180)                                                                    # rotation angle
    v.rot <- v * cos(theta) + Crossprod(axis, v) * sin(theta) + axis * (sum(axis * v)) * (1 - cos(theta))     # Rodrigues' rotation formula.
    return (v.rot)
  }

  Kernel.quartic <- function(vec){                                    # Kernel function. If 'vec' is a vector, it returns a vector.
    result <- c()
    for (i in 1:length(vec)){
      if (abs(vec[i]) <= 1){
        result[i] <- (1 - vec[i]^2)^2
      }else{
        result[i] <- 0
        }
      }
      return(result)
    }

  Kernel.indicator <- function(vec){                                  # kernel function. If 'vec' is vector, it returns vector.
    result <- c()
      for (i in 1:length(vec)){
        if (abs(vec[i]) <= 1){
          result[i] <- 1
        }else{
          result[i] <- 0
        }
      }
    return(result)
    }

  Kernel.Gaussian <- function(vec){                                   # kernel function. If 'vec' is a vector, it returns a vector.
    result <- c()
     for (i in 1:length(vec)){
       if (abs(vec[i]) <= 1){
         result[i] <- exp(-vec[i]^2)
       }else{
         result[i] <- 0
       }
     }
    return(result)
   }

  Cal.recon <- function(data, line){                                 # calculating reconstruction error from data to line.
    r <- 6378137                                                     # earth radius(m).
    rownames(data) <- NULL                                           # data and line should be n * 2 matrix or data frame.
    if (nrow(line) == 0){
      return (0)
    }else if (nrow(line) == 1){
      return (sum((geosphere::distGeo(data, line)/r)^2))
    }else{
      proj <- geosphere::dist2Line(data, line, distfun = geosphere::distGeo)
      return (sum((proj[, 1]/r)^2))
    }
  }

 ExtrinsicMean <- function(data, weights = rep(1, nrow(data))){                      # data: longitude/latitude matrix or dataframe with two column.
    if (sum(weights^2)^(1/2) == 0) stop("weight should not be the zero vector.")     # weights: a weight vector whose length is the rownumber of points.
    weights <- weights/sum(weights)                                                  # normalize 'weights' so that its components add up to 1.
    mean.Euc <- c(0, 0, 0)
    for (i in 1:nrow(data)){
      mean.Euc <- mean.Euc + Trans.Euclid(data[i, ]) * weights[i]
    }
    if (sum(mean.Euc^2)^(1/2) == 0){
      stop("extrinsic mean is not well-defined. Data should be more localized.")
    }
    mean <- Trans.sph(mean.Euc)
    return(as.numeric(mean))
  }

 IntrinsicMean <- function(data, weights = rep(1, nrow(data)), thres = 1e-5){        # data: longitude/latitude matrix or dataframe with two column.
   if (sum(weights^2)^(1/2) == 0) stop("weight should not be the zero vector.")      # weights: n-dimensional vector.
    mu <- data[1, , drop = F]                                                        # initialize mean as a point.
    delta.mu <- c(1, 0)
    while (sum((delta.mu)^2)^(1/2) > thres){
      weights.normal <- weights/sum(weights)                                         # normalize of 'weights' so that its components add up to 1.
      summation <- c(0, 0)
        for (m in 1:length(weights)){
          rot <- Rotate(mu, data[m, ])
          summation <- summation + weights.normal[m] * Logmap(rot)
        }
      delta.mu <- summation
      exp <- Expmap(delta.mu)
      mu.Euc <- Rotate.inv(mu, Trans.sph(exp))
      mu <- Trans.sph(mu.Euc)
     }
    return(as.numeric(mu))
  }

### PrincipalCircle
 PrincipalCircle <- function(data, step.size = 1e-3, thres = 1e-5, maxit = 10000){
    # Data: a matrix or data.frame of data points represented by angular form of (longitude, latitude).
    # To convergence of this algorithm, step.size is recommended below 0.01.

    circle <- IntrinsicMean(data) * pi/180                                          # initialized by center of data to avoid being trapped at local minima. (radian measurement)
    circle[3] <- pi/2
    data <- data * pi/180                                                           # transformation from angle into radian.
    delta <- 1
    iter <- 0
    while ((delta > thres) && (iter < maxit)){
      iter <- iter + 1
      grad <- c(0, 0, 0)
      for(i in 1:nrow(data)){
        temp <- sin(circle[2]) * sin(data[i, 2]) + cos(circle[2]) * cos(data[i, 2]) * cos(circle[1] - data[i, 1])
        names(temp) <- NULL
        if (temp == 1){
         next
        }
        grad[1] <- grad[1] + (acos(temp) - circle[3]) * 1/sqrt(1-temp^2) * cos(data[i, 2]) * cos(circle[2]) * sin(circle[1] - data[i, 1])
        grad[2] <- grad[2] - (acos(temp) - circle[3]) * 1/sqrt(1-temp^2) * (cos(circle[2]) * sin(data[i, 2]) - sin(circle[2]) * cos(data[i, 2]) * cos(circle[1] - data[i, 1]))
        grad[3] <- grad[3] + circle[3] - acos(temp)
      }
      grad[1:2] <- grad[1:2] * pi/180
      circle <- circle - step.size * 2 * grad

      if (circle[3] > pi) stop("step.size should be lowered.")

      delta <- sqrt(sum((step.size * 2 * grad)^2))

    }
    circle[1:2] <- circle[1:2] * 180/pi
    return(circle)                                                                 # Output: angular form of (longitude, latitude)
  }

GenerateCircle <- function(center, radius, T = 1000){
   # center should be an angular form of (longitude, latitude).
   # The radius r is in [0, pi]
   # The number of points in circle is T.
   lon <- seq(0, 360, length.out = T)
   generate.circle <- matrix(nrow = T, ncol = 2)

   for (i in 1:T){
     generate.circle[i, ] <- Trans.sph(Rotate.inv(center, c(lon[i], (pi/2 - radius) * 180/pi)))
   }
   return(generate.circle)
 }

 PGA <- function(data, col1 = "blue", col2 = "red"){
  # data should be a matrix or data frame with 2 column (longitude / latitude).
  center <- IntrinsicMean(data)
  ant <- Trans.sph(-Trans.Euclid(center))
  cov <- matrix(c(0, 0, 0, 0), nrow = 2)
  for (i in 1:nrow(data)){
      z <- Rotate(center, data[i, ])
      y <- Logmap(z)
      cov <- cov + y %*% t(y)
  }
  eivec <- eigen(cov)$vectors[, 1]
  direc <- eivec/(sum(eivec^2))^(1/2)
  pt <- Expmap(direc)
  pt.sph <- Trans.sph(pt)
  Euc <- Rotate.inv(center, pt.sph)
  mid <- Trans.sph(Euc)
  mid.ant <- Trans.sph(-Euc)
  line <- geosphere::makeLine(rbind(center, mid, ant, mid.ant, center))
  # plot
  sphereplot::rgl.sphgrid(col.long = "black", col.lat = "black")
  sphereplot::rgl.sphpoints(data, radius = 1, col = col1, size = 12)
  sphereplot::rgl.sphpoints(line, radius = 1, col = col2, size = 6)
  fit <- list(line = line)
  return(fit)
}


 Proj.Hauberg <- function(data, line){
   # This function performs Hauberg's projections.
   # Each data is projected into the nearest points of line (not exactly).
  proj <- matrix(nrow = nrow(data), ncol = 2)
    for (i in 1:nrow(data)){
      proj[i, ] <- line[which.min(geosphere::distGeo(data[i, ], line)), ]
    }
   return(proj)
 }

 Dist.pt <- function(data){
   # This function calculates the number of distinct points.
   # Data:  matrix or data.frame of data points represented by angular form of (longitude, latitude).
   num <- nrow(data)
   if (num == 1){
     return(1)
   }
   dist.num <- 1
   for (i in 2:num){
     temp <- c()
     temp <- colSums(t(data[1:(i - 1), , drop = F]) == data[i, ]) < 2
     prod <- 1
     for (j in 1:(i - 1)){
       prod <- prod * temp[j]
     }
     dist.num <- dist.num + prod
   }
   return(dist.num)
 }

 
########################################
 ### Spherical Principal Curves (global)
 #######################################
 SPC <- function(data, q = 0.1, T = nrow(data), step.size = 1e-3, maxit = 30, type = "Intrinsic", thres = 1e-2,
                 deletePoints = FALSE, plot.proj = FALSE, kernel = "quartic", 
                 col1 = "blue", col2 = "green", col3 = "red"){
   # Data: matrix or data.frame of data points represented by angular form of (longitude, latitude).
   # q: Smoothing parameter in Expectation step.
   # T: The number of points in circle.
   # step.size: It is recommended below 0.01 owing to convergence of this algorithm.
   # maxit: The maximum number of iterations
   # deletePoints: TRUE or FALSE. If it is TRUE, then for each expecatation step delete the points which do not have adjacent data.
   # If deletePoints is FALSE, leave the points which do not have adjacent data for each expectation step.
   # col1 is color of data and col2 is that of points in the principal curves; in addition, col3 is color of principal curves line.

  r <- 6378137                                                     # earth radius(m).
  PC <- PrincipalCircle(data, step.size = step.size)
  principal.curves <- GenerateCircle(PC[1:2], PC[3], T = T)        # Initialize principal curves as principal circle.

    if (nrow(principal.curves) == 1){
       residual <- sum((geosphere::distGeo(data, principal.curves)/r)^2)
     }else{
       proj <- geosphere::dist2Line(data, principal.curves, distfun = geosphere::distGeo)
       rownames(data) <- NULL
       residual <- Cal.recon(data, principal.curves)
     }
     residual.new <- residual.old <- residual
     relative.change <- 1
     iter <- 0                                                     # Iteration of principal curve algorithm.

     while ((relative.change >= thres) && (residual.old != 0) && (iter < maxit + 1)){
       ## projection step
       if (nrow(principal.curves) == 1) {
         proj.temp <- cbind(cbind(geosphere::distGeo(data, principal.curves), rep(principal.curves[, 1], nrow(data))), rep(principal.curves[, 2], nrow(data)))
       }else{
         proj.temp <- geosphere::dist2Line(data, principal.curves, distfun = geosphere::distGeo)
       }
       ## Expectation step
       T.temp <- nrow(principal.curves)
       curve.length <- 0
       for (i in 1:(T.temp - 1)){
         curve.length <- curve.length + geosphere::distGeo(principal.curves[i, ], principal.curves[i + 1, ])/r
       }
       if (deletePoints == FALSE){
         curve.length <- curve.length + geosphere::distGeo(principal.curves[T.temp, ], principal.curves[1, ]) # length of principal curves.
       }
       for (i in 1:T.temp){
         choice <- geosphere::distGeo(proj.temp[, 2:3], principal.curves[i, ])/r < q * curve.length
         adjacent <- as.data.frame(data[choice, , drop = F])                                                  # adjacent points for each dim.redu.PC.temp point.
         if (is.na(sum(choice)) == TRUE) {warning("Errors, 'choice' have not a number or non available.")}
         if (sum(choice) == 0){                                                                               # If a point dosen't have adjacent data, there are two options: deleting and leaving the point. (deletePoints)
           if (deletePoints == TRUE){
             principal.curves[i, ] <- c(NA, NA)                                                               # delete the points which do not have adjacent data each expectation step.
           }else{
             principal.curves[i, ] <- principal.curves[i, ]                                                   # leave the points which do not have adjacent data for each expectation step.
           }
         }else{
           if ((kernel == "quartic") | (kernel == "Quartic")){

             weights <- Kernel.quartic(geosphere::distGeo(proj.temp[choice, 2:3], principal.curves[i, ])/(r * q * curve.length))
           }else if ((kernel == "Gaussian") | (kernel == "gaussian")){
             weights <- Kernel.Gaussian(geosphere::distGeo(proj.temp[choice, 2:3], principal.curves[i, ])/(r * q * curve.length))
           }else if ((kernel == "indicator") | (kernel == "Indicator")){
             weights <- Kernel.indicator(geosphere::distGeo(proj.temp[choice, 2:3], principal.curves[i, ])/(r * q * curve.length))
           }else{
             stop("Errors, kernel should be quartic, Gaussian, and indicator.")
           }

           if ((type == "Intrinsic") | (type == "intrinsic")){
             mean <- IntrinsicMean(adjacent, weights)                                                         # calculating intrinsic mean of adjacent pts
             principal.curves[i, ] <- mean
           }else if ((type == "Extrinsic") | (type == "extrinsic")){
             mean <- ExtrinsicMean(adjacent, weights)                                                         # calculating extrinsic mean of adjacent pts
             principal.curves[i, ] <- mean
           }else{
             stop("Errors, type should be Intrinsic or Extrinsic")
           }
         
        }
       }
       principal.curves <- principal.curves[!is.na(principal.curves[, 1]), , drop = F]
       residual.new <- Cal.recon(data, principal.curves)
       relative.change <- abs((residual.old - residual.new)/residual.old)
       residual.old <- residual.new
       iter <- iter + 1
     }
     
     proj <- geosphere::dist2Line(data, principal.curves)
     distnum <- Dist.pt(proj[, 2:3, drop = F])                                              # The number of distinct projections
     
     if (nrow(principal.curves) == 1){
         line <- principal.curves
     }else{ 
     line <- geosphere::makeLine(principal.curves)
     } 
     if (iter < maxit){
        converged <- TRUE
     }else{
        converged <- FALSE 
     }
     sphereplot::rgl.sphgrid(col.long = "black", col.lat = "black")                                                             # plot the data and principal curves.
     sphereplot::rgl.sphpoints(data[, 1], data[, 2], radius = 1, col = col1, size = 12)
     sphereplot::rgl.sphpoints(principal.curves[, 1], principal.curves[, 2], radius = 1, col = col2, size = 6)
     sphereplot::rgl.sphpoints(line[, 1], line[, 2], radius = 1, col = col3, size = 6)

     if (plot.proj == TRUE){
     line.proj <- matrix(ncol = 2)
       for (i in 1:nrow(proj)){                                                            # plot the projection line
         if (abs(sum((data[i, ] - proj[i, 2:3])^2)) < 10^(-10)){
           line.proj <- data[i, , drop = F]
         }else{
           line.proj <- geosphere::makeLine(rbind(data[i, ], data[i, ], proj[i, 2:3]))
         }
         sphereplot::rgl.sphpoints(line.proj[, 1], line.proj[, 2], radius = 1, col = "black")
       }
     }
     
     fit <- list(prin.curves = principal.curves, line = line, converged = converged, iteration = iter,
                 recon.error = residual.new, num.dist.pt = distnum)
     return(fit)
}

 
 SPC.Hauberg <- function(data, q = 0.1, T = nrow(data), step.size = 1e-3, maxit = 30, type = "Intrinsic",
                         thres = 1e-2, deletePoints = FALSE, plot.proj = FALSE, kernel = "quartic",
                         col1 = "blue", col2 = "green", col3 = "red"){
  # Data: matrix or data.frame of data points represented by angular form of (longitude, latitude).
  # q: Smoothing parameter in Expectation step.
  # T: The number of points in circle.
  # step.size: step size of PrincipalCircle. It is recommended below 0.01 owing to convergence of this algorithm.
  # maxit: The number of iterations
  # deletePoints: TRUE or FALSE. If it is TRUE, then for each expecatation step delete the points which do not have adjacent data.
  # If deletePoints is FALSE, leave the points which do not have adjacent data for each expectation step.
  # col1 is color of data and col2 is that of points in the principal curves; in addition, col3 is color of principal curves line.

 r <- 6378137                                                        # earth radius(m).
 PC <- PrincipalCircle(data, step.size = step.size)
 principal.curves <- GenerateCircle(PC[1:2], PC[3], T = T)           # principal circle is set to be an initialization of principal curves.
 if (nrow(principal.curves) == 1){
   residual <- sum((geosphere::distGeo(data, principal.curves)/r)^2)
 }else{
    proj <- Proj.Hauberg(data, principal.curves)
    residual <- sum((geosphere::distGeo(data, proj)/r)^2)            # Hauberg's residual.
  }
 residual.new <- residual.old <- residual
 relative.change <- 1
 iter <- 0                                                           # iteration of principal curve algorithm.

 while ((relative.change >= thres) && (residual.old != 0) && (iter < maxit + 1)){
   ## projection step
   if (nrow(principal.curves) == 1){
     proj.temp <- cbind(geosphere::distGeo(data, principal.curves), rep(principal.curves[, 1], nrow(data)), rep(principal.curves[, 2], nrow(data)))
   } else {
     proj.temp <- Proj.Hauberg(data, principal.curves)
   }
   ## Expectation step
   T.temp <- nrow(principal.curves)
   curve.length <- 0
   for (i in 1:(T.temp - 1)){
     curve.length <- curve.length + geosphere::distGeo(principal.curves[i, ], principal.curves[i + 1, ])/r
   }
   if (deletePoints == FALSE){
     curve.length <- curve.length + geosphere::distGeo(principal.curves[T.temp, ], principal.curves[1, ]) # length of principal curves.
   }
   for (i in 1:T.temp){
     choice <- geosphere::distGeo(proj.temp, principal.curves[i, ])/r < q * curve.length
     adjacent <- as.data.frame(data[choice, , drop = F])                                                  # adjacent points for each dim.redu.PC.temp point.
     if (is.na(sum(choice)) == TRUE) {warning("Errors, object 'choice' have not a number or non available.")}
     if (sum(choice) == 0){                                                                               # If a point dosen't have adjacent data, there are two options: deleting and leaving the point. (deletePoints)
       if (deletePoints == TRUE){
          principal.curves[i, ] <- c(NA, NA)                                                              # delete the points which do not have adjacent data for each expectation step.
       }else{
          principal.curves[i, ] <- principal.curves[i, ]                                                  # leave the points which do not have adjacent data for each expectation step.
       }

     }else{
       if ((kernel == "quartic") | (kernel == "Quartic")){
         weights <- Kernel.quartic(geosphere::distGeo(proj.temp[choice, , drop = F], principal.curves[i, ]) / (r * q * curve.length))
       }else if ((kernel == "Gaussian") | (kernel == "gaussian")){
         weights <- Kernel.Gaussian(geosphere::distGeo(proj.temp[choice, , drop = F], principal.curves[i, ]) / (r * q * curve.length))
       }else if ((kernel == "indicator") | (kernel == "Indicator")){
         weights <- Kernel.indicator(geosphere::distGeo(proj.temp[choice, , drop = F], principal.curves[i, ]) / (r * q * curve.length))
       }else{
         stop("Errors, kernel should be quartic, Gaussian, and indicator.")
       }
       
       if ((type == "Intrinsic") | (type == "intrinsic")){
         mean <- IntrinsicMean(adjacent, weights)                                                         # calculating intrinsic mean of adjacent pts
         principal.curves[i, ] <- mean
       }else if ((type == "Extrinsic") | (type == "extrinsic")){
         mean <- ExtrinsicMean(adjacent, weights)                                                         # calculating extrinsic mean of adjacent pts
         principal.curves[i, ] <- mean
       }else{
         stop("Errors, type should be Intrinsic or Extrinsic")
       }
       
     }
   }
   principal.curves <- principal.curves[!is.na(principal.curves[, 1]), , drop = F]                        # deleting NA rows.
   proj.temp <- Proj.Hauberg(data, principal.curves)
   residual.new <- sum((geosphere::distGeo(data, proj.temp)/r)^2)
   relative.change <- abs((residual.old - residual.new)/residual.old)                                     # relative change. (stop when it is below 'thres')
   residual.old <- residual.new
   iter <- iter + 1
 }
 proj <- Proj.Hauberg(data, principal.curves)
 distnum <- Dist.pt(proj)                                                                                 # The number of distinct projections.
 
 if (nrow(principal.curves) == 1){
   line <- principal.curves
 }else{
 line <- geosphere::makeLine(principal.curves)
 }
 if (iter < maxit){
   converged <- TRUE
 }else{
   converged <- FALSE
 }
 sphereplot::rgl.sphgrid(col.long = "black", col.lat = "black")                                                                                # Plot the data and principal curves.
 sphereplot::rgl.sphpoints(data[, 1], data[, 2], radius = 1, col = col1, size = 12)
 sphereplot::rgl.sphpoints(principal.curves[, 1], principal.curves[, 2], radius = 1, col = col2, size = 6)
 sphereplot::rgl.sphpoints(line[, 1], line[, 2], radius = 1, col = col3, size = 6)
 
 if (plot.proj == TRUE){
   line.proj <- matrix(ncol = 2)
   for (i in 1:nrow(proj)){                                                                               # Plot the projection line
     if (abs(sum((data[i, ] - proj[i, ])^2)) < 10^(-10)){
       line.proj <- data[i, , drop = F]
     }else{
       line.proj <- geosphere::makeLine(rbind(data[i, ], data[i, ], proj[i, ]))
     }
     sphereplot::rgl.sphpoints(line.proj[, 1], line.proj[, 2], radius = 1, col = "black")
   }
 }
 fit <- list(prin.curves = principal.curves, line = line, converged = converged, iteration = iter,
                recon.error = residual.new, num.dist.point = distnum)
 return(fit)
}


  ######################################
  ### Local principal geodesics (Local)
  ######################################

  LPG <- function(data, scale = 0.04, tau = scale/3, nu = 0, maxpt = 500, seed = NULL,
         kernel = "indicator", thres = 1e-4, col1 = "blue", col2 = "green", col3 = "red"){
    #  Data   a matrix or data.frame of data points represented by angular form of (longitude, latitude).
    # 'scale' is a scale parameter (LPG scaleing region).
    #  For each LPG step, the curve proceed by 'tau'. (step size of the algorithm)
    # 'nu' represents viscosity of the given data.
    # 'thres' is a threshold in IntrinsicMean function.
    set.seed(seed)                                                        # set seed number.
    nu <- nu
    h <- scale
    tau <- tau                                                            # step size of the algorithm.
    r <- 6378137                                                          # earth radius(m).
    H <- r * h                                                            # scaleing distance on earth(m).
    num.curves <- 0                                                       # the number of the resulting curves.
    # create blank matrix and list.
    dim.redu <- matrix(NA, nrow = 10000, ncol = 2)
    line <- matrix(NA, nrow = 10000, ncol = 2)
    line2 <- matrix(NA, nrow = 10000, ncol = 2)
    dim.redu.LPG <- list()
    data.raw <- data
    
    ### start LPG.
    repeat{
      if (nrow(data) == 0){
        break
      }

      # first step-------------------------------------------------------------------------------
      num.data <- nrow(data)                                              # the number of data.
      X <- Y <- matrix(nrow = 10000, ncol = 2)                            # limit nrow set 10000.
      ran.num <- sample(1:num.data, 1)                                    # choose initial number randomly.
      Y[1, ] <- X[1, ] <- as.matrix(data[ran.num, , drop = F])            # random sampling X1 from data.
      cov.local <- S <- matrix(c(0, 0, 0, 0), nrow = 2)
      num.related <- ker <- ker.sum <- 0
      direc <- direc2 <- matrix(nrow = 10000, ncol = 2)                   # limit of repeat is 10000.
      cohesive <- matrix(0, nrow = 2, ncol = 1)
      mean.initial <- c(0, 0, 0)

      adjacent <- data[geosphere::distGeo(X[1, ], data) < H, , drop = F]  # Find the center(intrinsic mean) of adjacent points.
      X[1, ] <- Y[1, ] <- IntrinsicMean(adjacent, rep(1, nrow(adjacent)), thres)

      for (i in 1:num.data){
        dist <- geosphere::distGeo(X[1, ], data[i, ])
        if (dist < H) {
          adj <- Rotate(c(X[1, 1], X[1, 2]), c(data[i, 1], data[i, 2]))
          adj.plane <- Logmap(adj)
          if ((kernel == "indicator") | (kernel == "Indicator")){
            ker <- Kernel.indicator(dist/H)
          }else if ((kernel == "quartic") | (kernel == "Quartic")){
            ker <- Kernel.quartic(dist/H)
          }else if ((kernel == "Gaussian") | (kernel == "gaussian")){
            ker <- Kernel.Gaussian(dist/H)
          }else{
            stop("kernal should be either quartic, indicator or Gaussian")
          }
          cohesive <- cohesive + adj.plane * ker                          # cohesive force.
          S <- S + adj.plane %*% t(adj.plane) * ker
          num.related <- num.related + 1
          ker.sum <- ker.sum + ker
        }
      }
      if (num.related == 0){
        dim.redu.temp <- t(as.matrix(X[1, ]))
        dim.redu  <- rbind(as.matrix(na.omit(dim.redu)), dim.redu.temp)
        line <- rbind(as.matrix(na.omit(line)), dim.redu.temp)
        data <- data[-ran.num, , drop = F]
      }else{
        cov.local <- S/ker.sum                                            # locally kernel-weighted tangent covariance at scale h.
        eivec <- eigen(cov.local)$vectors[, 1]                            # leading eigen vector of local tangent covariance.
        if (sum(cohesive^2) > 0){
          cohesive <- cohesive/(sum(cohesive^2))^(1/2)                    # making cohesive unit vector.
        }else{
          cohesive <- c(0, 0)
        }
        force <- (cohesive * nu + eivec)                                  # summation of cohesive and maximal variation and 'nu' is viscous parameter.
        direc[1, ] <- force/(sum(force^2))^(1/2)                          # normalized forward first direction.
        direc2[1, ] <- -direc[1, ]                                        # normalized backward first direction.
        forward <- tau * direc[1, ]
        backward <- tau * direc2[1, ]
        pt <- Expmap(forward)
        pt.sph <- Trans.sph(pt)
        pt.2 <- Expmap(backward)
        pt.2.sph <- Trans.sph(pt.2)
        Euc <- Rotate.inv(X[1, ], c(pt.sph[1], pt.sph[2]))                # Euclidean points(R^3) of X[2, ].
        Euc2 <- Rotate.inv(X[1, ], c(pt.2.sph[1], pt.2.sph[2]))
        X[2, ] <- Trans.sph(Euc)                                          # get spherical coordinate (longitude,latitude).
        Y[2, ] <- Trans.sph(Euc2)

        # forward diriection step---------------------------------------------------------------------------
        j <- 2                                                            # iteration of forward direction step.
        repeat{
          num.related <- ker <- ker.sum <- 0
          cov.local <- S <- matrix(c(0, 0, 0, 0), nrow = 2)
          cohesive <- sum.adj.plane <- c(0, 0)
          for (i in 1:num.data){
            dist <- geosphere::distGeo(X[j, ], data[i, ])
            adj.plane <- c(0, 0)
            if (dist < H){
              adj <- Rotate(c(X[j, 1], X[j, 2]), c(data[i, 1], data[i, 2]))                 # adjacent data of X[j, ].
              adj.plane <- Logmap(adj)                                                      # tangent plane coordinate of adjacent data.
              if ((kernel == "indicator") | (kernel == "Indicator")){
                ker <- Kernel.indicator(dist/H)
              }else if ((kernel == "quartic") | (kernel == "Quartic")){
                ker <- Kernel.quartic(dist/H)
              }else if ((kernel == "Gaussian") | (kernel == "gaussian")){
                ker <- Kernel.Gaussian(dist/H)
              }else{
                stop("kernal should be either quartic, indicator or Gaussian")
              }
              cohesive <- cohesive + adj.plane * ker
              S <- S + adj.plane %*% t(adj.plane) * ker
              ker.sum <- ker.sum + ker
              num.related <- num.related + 1
            }
            sum.adj.plane <- sum.adj.plane + adj.plane
          }
          if ((j > maxpt) | (num.related == 0)){
            break
          }else{
            loc.center <- sum.adj.plane/num.related                                         # center of adjacent data. (local center)
            cov.local <- S/ker.sum - (loc.center) %*% t(loc.center)                         # locally kernel-weighted tangent covariance at scale h.
            if (sum(cohesive^2) > 0){
              cohesive <- cohesive/(sum(cohesive^2))^(1/2)                                  # making cohesive unit vector.
            }else{
              cohesive <- c(0, 0)
            }
            eivec <- eigen(cov.local)$vectors[, 1]                                          # leading eigenvector of local tangent covariance.
            force1 <- cohesive * nu + eivec                                                # summation of cohesive and maximal variation.
            force2 <- cohesive * nu - eivec                                               
            if (sum(direc[j - 1, ] * eivec) >= 0){
              direc[j, ] <- force1/(sum(force1^2))^(1/2)                                    # making unit vector.
            }else{
              direc[j, ] <- force2/(sum(force2^2))^(1/2)
            }
            forward <- tau * direc[j, ]
            pt <- Expmap(forward)
            pt.sph <- Trans.sph(pt)
            Euc <- Rotate.inv(X[j, ], c(pt.sph[1], pt.sph[2]))                              # Euclidean point of X[j, ].
            X[j + 1, ] <- Trans.sph(Euc)                                                    # get a spatial coordinate.
            j <- j + 1
          }
        }
        X <- na.omit(X)

        # backward direction step-----------------------------------------------------------------
        j <- 2                                                                              # iteration of backward direction step.
        repeat{
          num.related <- ker <- ker.sum <- 0
          cov.local <- S <- matrix(c(0, 0, 0, 0), nrow = 2)
          cohesive <- matrix(0, nrow = 2, ncol = 1)
          sum.adj.plane <- c(0, 0)
          for (i in 1:num.data){
            dist <- geosphere::distGeo(Y[j, ], data[i, ])
            adj.plane <- c(0, 0)
            if (geosphere::distGeo(Y[j, ], data[i, ]) < H){
              adj <- Rotate(c(Y[j, 1], Y[j, 2]), c(data[i, 1], data[i, 2]))                 # adjacent data of Y[j, ].
              adj.plane <- Logmap(adj)                                                      # tangent plane coordinate of adjacent data.
               if ((kernel == "indicator") | (kernel == "Indicator")){
                 ker <- Kernel.indicator(dist/H)
               }else if ((kernel == "quartic") | (kernel == "Quartic")){
                 ker <- Kernel.quartic(dist/H)
               }else if ((kernel == "Gaussian") | (kernel == "gaussian")){
                 ker <- Kernel.Gaussian(dist/H)
               }else{
                 stop("kernal should be either quartic, indicator or Gaussian")
               }
               cohesive <- cohesive + adj.plane * ker
               S <- S + adj.plane %*% t(adj.plane) * ker
               ker.sum <- ker.sum + ker
               num.related <- num.related + 1
             }
             sum.adj.plane <- sum.adj.plane + adj.plane 
           }
          if ((j > maxpt) | (num.related == 0)){                                        # stop procedure when there is no neighborhood points.
             break
           }else{
             loc.center <- sum.adj.plane/num.related                                        # center of adjacent data. (local center)
             cov.local <- S/ker.sum - (loc.center) %*% t(loc.center)                        # local tangent covariance at scale h.
             eivec <- eigen(cov.local)$vectors[, 1]                                         # leading eigen vector of local tangent covariance.
             if (sum(cohesive^2) > 0){
               cohesive <- cohesive/(sum(cohesive^2))^(1/2)                                 # making cohesive unit vector.
             }else{
               cohesive <- c(0, 0)
             }
             force1 <- cohesive * nu + eivec                                               # summation of cohesive and maximal variation.
             force2 <- cohesive * nu - eivec                                     
             if (sum(direc2[j - 1, ] * eivec) >= 0){
               direc2[j, ] <- force1/(sum(force1^2))^(1/2)
             }else{
               direc2[j, ] <- force2/(sum(force2^2))^(1/2)
             }
             backward <- tau * direc2[j, ]
             pt <- Expmap(backward)
             pt.sph <- Trans.sph(pt)
             Euc <- Rotate.inv(Y[j, ], c(pt.sph[1], pt.sph[2]))                             # Euclidean point of Y[j, ].
             Y[j + 1, ] <- Trans.sph(Euc)                                                   # transform three-dimensional Euclidean coordinate into spherical coordinate (longitude,latitude).
             j <- j + 1
            }
          }
         Y <- na.omit(Y)                                                                    # deleting NA.
         Y.rev <- Y                                                                         # Y row reverse.
         for (i in 1:nrow(Y.rev)){
           temp <- nrow(Y.rev)
           Y.rev[i, ] <- Y[temp + 1 - i, ]
         }
         dim.redu.temp <- rbind(Y.rev, X)                                                   # connected componemt of dimension reduction.
         dim.redu <- rbind(as.matrix(na.omit(dim.redu)), dim.redu.temp)                     # sum of connected component of dimension reduction.
         line <- rbind(as.matrix(na.omit(line)), geosphere::makeLine(dim.redu.temp))        # line of dim.redu.
         proj <- geosphere::dist2Line(data, dim.redu, distfun = geosphere::distGeo)         # proj[, 2:3] are projection point.
         data <- data[proj[, 1] > H, , drop = F]                                            # remaining data set.
       }
       dim.redu.LPG[[num.curves + 1]] <- dim.redu.temp                                        # storage of dim.redu.temp.
       num.curves <- num.curves + 1
     }

     sphereplot::rgl.sphgrid(col.long = "black", col.lat = "black")                                                              # plot the result.
     sphereplot::rgl.sphpoints(data.raw, radius = 1, col = col1, size = 12)
     sphereplot::rgl.sphpoints(dim.redu, radius = 1, col = col2, size = 4)
     sphereplot::rgl.sphpoints(line, radius = 1, col = col3, size = 6)
     
     fit <- list(num.curves = num.curves, prin.curves = dim.redu.LPG, line = line)
     return(fit)
     }

Try the spherepc package in your browser

Any scripts or data that you put into this service are public.

spherepc documentation built on March 25, 2020, 5:07 p.m.