循环浏览文件，获取theta并按R中的theta旋转曲线 [英] Loop through files, obtain theta and rotate curve by theta in R

问题描述

我问了一个问题

I asked a question here where Miff helped me figure out how to find the point on a curve which is perpendicular to the line between the endpoints at its midpoint.

To do this, the user pointed out it was necessary to rotate the curve by the gradient of the line joining the endpoints so that the line is flat, using approx, and then rotate in the opposite direction using the lava package. I can make this work on a case by case basis, establishing theta and working from there. I'm having really bad luck trying to embed this within a function, however.

I am having difficulty using dplyr to rotate each set of 42 points by a theta value within a function.

Here's a sample set of data. The real data has hundreds of curves I need to work through.

data <- structure(list(X = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, 4.9046, 
6.1424, 7.275, 8.5851, 10.0373, 11.9981, 13.7726, 15.0731, 16.0664, 
18.1945, 21.2666, 24.2093, 26.7119, 28.8037, 30.7135, 32.1351, 
33.1982, 34.2341, 35.7587, 37.2147, 38.4303, 39.625, 40.4596, 
42.0938, 42.7428, 42.7593, 43.5085, 43.7419, 43.5989, 44.0841, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, -14.845, -11.9052, 
-8.7897, -5.8034, -2.6756, 0.3316, 3.4003, 6.5281, 9.6517, 12.804, 
15.9861, 19.1769, 22.2929, 25.4089, 28.3392, 31.0054, 33.1847, 
35.081, 36.7227, 38.1544, 39.1697, 40.049, 40.9647, 41.5014, 
41.8874, 42.1778, 42.3435, 42.2681, 42.3745, 42.4619, NA, NA, 
NA, NA), Y = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, -9.9938, -7.4596, 
-4.8647, -2.2903, 0.3158, 2.9302, 5.7262, 8.7033, 11.8007, 14.9847, 
16.7225, 16.7813, 15.6921, 14.2964, 11.5579, 8.2378, 5.183, 1.5938, 
-2.0712, -5.195, -7.1447, -9.0446, -11.1269, -13.0979, -15.3295, 
-17.1898, -19.4376, -21.4781, -23.8426, -25.6343, NA, NA, NA, 
NA, NA, NA, NA, NA, NA, NA, NA, 8.0113, 9.1826, 9.838, 10.7908, 
11.175, 12.0393, 12.6813, 12.8828, 13.2281, 13.5102, 13.6637, 
13.5493, 12.8699, 12.2191, 10.9208, 9.0209, 6.2158, 3.2466, 0.2169, 
-2.7807, -6.0439, -9.1262, -11.8684, -14.7779, -17.5825, -20.2452, 
-22.807, -25.3519, -27.6105, -29.7536, NA, NA, NA, NA), fan_line = c(1L, 
2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 14L, 15L, 
16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 27L, 28L, 
29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L, 40L, 41L, 
42L, 1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 11L, 12L, 13L, 
14L, 15L, 16L, 17L, 18L, 19L, 20L, 21L, 22L, 23L, 24L, 25L, 26L, 
27L, 28L, 29L, 30L, 31L, 32L, 33L, 34L, 35L, 36L, 37L, 38L, 39L, 
40L, 41L, 42L)), class = "data.frame", row.names = c(NA, -84L
))

Currently within the function, I loop through each set of 42 XY coordinates that constitute my curves of interest, and obtain the start and end XY coordinates for each curve:

plyr, dplyr, tidyr and lava are loaded

data <- data %>% mutate(Group = rep(1:(n()/42), each = 42)) %>% dplyr::group_by(Group) %>% 
    mutate(start = min(which(!is.na(X))), end = max(which(!is.na(X))), midpoint = round((start+end)/2, digits = 0)) %>% ungroup()

for (i in 1:nrow(data)){
    if (data[i, "fan_line"] == data[i, "start"]){
      data[i, "start_val_x"] = data[i, "X"]
      data[i, "start_val_y"] = data[i, "Y"]
    }
    else{data[i, "start_val_y"] = NA
    data[i, "start_val_x"] = NA}
  }
  
  for (i in 1:nrow(data)){
    if (data[i, "fan_line"] == data[i, "end"]){
      data[i, "end_val_x"] = data[i, "X"]
      data[i, "end_val_y"] = data[i, "Y"]
    }
    else{data[i, "end_val_y"] = NA
    data[i, "end_val_x"] = NA}
  }

data <- data %>%  group_by(Group) %>% fill(c(start_val_x, start_val_y), .direction = "down") %>% fill(c(start_val_x, start_val_y), .direction = "up")
data2 <- data %>%  group_by(Group) %>% fill(c(end_val_x, end_val_y), .direction = "down") %>% fill(c(end_val_x, end_val_y), .direction = "up")

Here I find points along a line running between the start and endpoints along each curve. In this cases, it is the one third point.

data3 <- data2 %>% group_by(Group) %>% mutate(x_value_one_third_hypo = ((start_val_x/3)*2 + (end_val_x*0.33)),
                                            y_value_one_third_hypo = ((start_val_y/3)*2 + (end_val_y*0.33)),
                                            third_x = approx(X, Y, (min(X, na.rm = T)/3)*2 + (max(X, na.rm = T)/3))$x,
                                            third_y = approx(X, Y, (min(X, na.rm = T)/3)*2 + (max(X, na.rm = T)/3))$y)

I calculate theta for each curve like so:

data3 <- data3 %>% group_by(Group) %>% mutate(theta = max(atan(diff(c(start_val_y, end_val_y))/diff(c(start_val_x, end_val_x))), na.rm = T))

However, I then run into the problem of trying to rotate using this theta value - I get an error message informing me that object theta is not found.

  data3.5 <- data3 %>% bind_cols(as_tibble(rotate2(as.matrix(.)[,1:2], theta = theta)))

I'm not sure how I can rotate the coordinates of each set of xy coordinates by their corresponding theta values.

解决方案

I'm sure there's an easier way to do this, but it works. That's assuming that you're trying to rotate about the origin. If you wanted to use an arbitrary axis, this won't be right.

I did have to change the data to work through this because there is no 80 to match to end, so I changed it to 38 just to work through the code.

I thought that the lava function required a square matrix, but that didn't work, so I went back to the more traditional rotation matrix.

Using the theta you calculated, I created this matrix in R

# create the 2D rotation matrix
tMat = rbind(c(cos(unique(data3$theta)), 
               -1* sin(unique(data3$theta))),
             c(sin(unique(data3$theta)), 
               cos(unique(data3$theta))))

# an empty list to store the results
trXY <- vector("list")   # transformed xy storage

# transform points - extract x & y, put into a vector, then matrix mult
for(i in 1:nrow(data3)){
  # collect original x and y
  coordData = data3[i, 1:2] %>% unlist()
 
  # skip the NA rows
  if(is.na(data3[i, 1]) == F) {
     # calculate
     results <- tMat %*% coordData
     # store the results
     trXY[i] <- list(results)
   } # end if
   else{
     # is NA
     trXY[i] <- NA
   } # end else
 } # end for

 # check output
 trXY[38] # it's stored in rows x is row 1; y is row 2

 # extract coordinates from the list and put them in the data frame
 # x coord
 data3$trX <- apply(trXY, function(x) ifelse(length(x) == 1,
                                             NA,
                                             x[[1]])) %>% unlist()
 # y coord
 data3$trY <- apply(trXY, function(x) ifelse(length(x) == 1,
                                             NA,
                                             x[[2]])) %>% unlist()

  # check rotation
  (p <- plotly::plot_ly(data3, 
                        x = ~X, 
                        y = ~Y, 
                        name = "original data", 
                        type = "scatter", 
                        mode = "lines", 
                        color = ~Group %>% as.factor(), 
                        colors = c("#cd0c18","#1660a7")) %>% 
     plotly::add_trace(x = ~trX, 
                       y = ~trY, 
                       name = "transformed", 
                       mode = "lines", 
                       color = ~Group %>% as.factor(), 
                       line = list(dash = "dash")))

This is what it looks like:

这篇关于循环浏览文件，获取theta并按R中的theta旋转曲线的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！