R+RGL: plot of spheres and segments, segments do not zoom properly

Question

Please consider the snippet below. It plots a set of spheres connected by some segments. The function to draw the smooth spheres comes from the discussion at

How to increase smoothness of spheres3d in rgl

What puzzles me is the following: when I zoom in/out the RGL plot, the spheres and the segments behave differently. In particular, if I zoom in, the segments look rather thin with respect to the spheres, whereas they look really wide when I zoom out.

Is there a way to correct this behavior, so that the proportion between the spheres and the segments is always respected regardless of the zoom level? Thanks a lot

library(rgl)
library(tidyverse)


## a function to plot good-looking spheres

sphere1.f <- function(x0 = 0, y0 = 0, z0 = 0, r = 1, n = 101, ...){
  f <- function(s,t){ 
    cbind(   r * cos(t)*cos(s) + x0,
             r *        sin(s) + y0,
             r * sin(t)*cos(s) + z0)
  }
  persp3d(f, slim = c(-pi/2,pi/2), tlim = c(0, 2*pi), n = n, add = T, ...)
}



## a set of 3D coordinates for my spheres

agg <- structure(list(X1 = c(-0.308421860438279, -1.42503395393061, 
1.10667871416591, -0.41759848570565, 0.523721760757519, 0.520653825151111, 
4.54213267745731, 2.96469370222004, 6.32495200153492, 3.78715565912871, 
5.35968114482443), X2 = c(0.183223776337368, 1.69719822686475, 
-0.992839275466541, 2.22182475540691, -0.705817674534376, -2.40358980860811, 
-0.565561169031234, -0.362260309907445, 0.326094711744554, 0.60340188259578, 
-0.00167511540165435), X3 = c(-0.712687792799106, -0.0336746884947792, 
0.0711272759107127, 1.6126544944538, -2.29999319137504, 1.36257390230441, 
-1.52942342176029, -0.316841449239697, -1.69222713171002, 1.23000775530984, 
2.30848424740017)), class = c("spec_tbl_df", "tbl_df", "tbl", 
"data.frame"), row.names = c(NA, -11L), spec = structure(list(
    cols = list(X1 = structure(list(), class = c("collector_double", 
    "collector")), X2 = structure(list(), class = c("collector_double", 
    "collector")), X3 = structure(list(), class = c("collector_double", 
    "collector"))), default = structure(list(), class = c("collector_guess", 
    "collector")), skip = 0), class = "col_spec"))



##coordinares of the segments (bonds) connecting  the spheres

bond_segments <- structure(list(X1 = c(-1.42503395393061, -0.308421860438279, 
1.10667871416591, -0.308421860438279, 0.523721760757519, -0.308421860438279, 
-0.41759848570565, -1.42503395393061, 0.520653825151111, 1.10667871416591, 
2.96469370222004, 1.10667871416591, 2.96469370222004, 4.54213267745731, 
6.32495200153492, 4.54213267745731, 3.78715565912871, 2.96469370222004, 
5.35968114482443, 3.78715565912871), X2 = c(1.69719822686475, 
0.183223776337368, -0.992839275466541, 0.183223776337368, -0.705817674534376, 
0.183223776337368, 2.22182475540691, 1.69719822686475, -2.40358980860811, 
-0.992839275466541, -0.362260309907445, -0.992839275466541, -0.362260309907445, 
-0.565561169031234, 0.326094711744554, -0.565561169031234, 0.60340188259578, 
-0.362260309907445, -0.00167511540165435, 0.60340188259578), 
    X3 = c(-0.0336746884947792, -0.712687792799106, 0.0711272759107127, 
    -0.712687792799106, -2.29999319137504, -0.712687792799106, 
    1.6126544944538, -0.0336746884947792, 1.36257390230441, 0.0711272759107127, 
    -0.316841449239697, 0.0711272759107127, -0.316841449239697, 
    -1.52942342176029, -1.69222713171002, -1.52942342176029, 
    1.23000775530984, -0.316841449239697, 2.30848424740017, 1.23000775530984
    )), class = c("spec_tbl_df", "tbl_df", "tbl", "data.frame"
), row.names = c(NA, -20L), spec = structure(list(cols = list(
    X1 = structure(list(), class = c("collector_double", "collector"
    )), X2 = structure(list(), class = c("collector_double", 
    "collector")), X3 = structure(list(), class = c("collector_double", 
    "collector"))), default = structure(list(), class = c("collector_guess", 
"collector")), skip = 0), class = "col_spec"))


open3d()
#> glX 
#>   1


##material and light effects for the spheres

material3d(ambient = "black", specular = "grey60", emission = "black", shininess = 30.0)
clear3d(type = "lights")
light3d(theta = -30, phi = 60, viewpoint.rel = TRUE, ambient = "#FFFFFF", diffuse = "#FFFFFF", specular = "#FFFFFF", x = NULL, y = NULL, z = NULL)
light3d(theta = -0, phi = 0, viewpoint.rel = TRUE,  diffuse = "gray20", specular = "gray25", ambient = "gray80", x = NULL, y = NULL, z = NULL)

## plot the spheres
agg %>%
  rowwise() %>%
  mutate(spheres = sphere1.f(X1, X2, X3, r=0.5))
#> # A tibble: 11 × 4
#> # Rowwise: 
#>        X1       X2      X3 spheres   
#>     <dbl>    <dbl>   <dbl> <rglLwlvl>
#>  1 -0.308  0.183   -0.713  15        
#>  2 -1.43   1.70    -0.0337 16        
#>  3  1.11  -0.993    0.0711 17        
#>  4 -0.418  2.22     1.61   18        
#>  5  0.524 -0.706   -2.30   19        
#>  6  0.521 -2.40     1.36   20        
#>  7  4.54  -0.566   -1.53   21        
#>  8  2.96  -0.362   -0.317  22        
#>  9  6.32   0.326   -1.69   23        
#> 10  3.79   0.603    1.23   24        
#> 11  5.36  -0.00168  2.31   25

## add the segments
segments3d(bond_segments, lwd=8, color="black")

sessionInfo()
#> R version 4.1.2 (2021-11-01)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 11 (bullseye)
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/blas/libblas.so.3.9.0
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.9.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
#>  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
#>  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] forcats_0.5.1   stringr_1.4.0   dplyr_1.0.8     purrr_0.3.4    
#>  [5] readr_2.1.1     tidyr_1.2.0     tibble_3.1.6    ggplot2_3.3.5  
#>  [9] tidyverse_1.3.1 rgl_0.108.3    
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.7        lubridate_1.8.0   assertthat_0.2.1  digest_0.6.29    
#>  [5] utf8_1.2.2        R6_2.5.1          cellranger_1.1.0  backports_1.4.1  
#>  [9] reprex_2.0.1      evaluate_0.14     httr_1.4.2        highr_0.9        
#> [13] pillar_1.6.4      rlang_1.0.1       readxl_1.3.1      R.utils_2.11.0   
#> [17] R.oo_1.24.0       rmarkdown_2.11    styler_1.6.2      htmlwidgets_1.5.4
#> [21] munsell_0.5.0     broom_0.7.10      compiler_4.1.2    modelr_0.1.8     
#> [25] xfun_0.29         pkgconfig_2.0.3   htmltools_0.5.2   tidyselect_1.1.1 
#> [29] fansi_0.5.0       crayon_1.4.2      tzdb_0.2.0        dbplyr_2.1.1     
#> [33] withr_2.4.3       R.methodsS3_1.8.1 grid_4.1.2        jsonlite_1.7.2   
#> [37] gtable_0.3.0      lifecycle_1.0.1   DBI_1.1.2         magrittr_2.0.1   
#> [41] scales_1.1.1      cli_3.1.0         stringi_1.7.6     fs_1.5.2         
#> [45] xml2_1.3.3        ellipsis_0.3.2    generics_0.1.1    vctrs_0.3.8      
#> [49] tools_4.1.2       R.cache_0.15.0    glue_1.6.0        hms_1.1.1        
#> [53] fastmap_1.1.0     yaml_2.2.1        colorspace_2.0-2  rvest_1.0.2      
#> [57] knitr_1.37        haven_2.4.3

^{Created on 2022-03-06 by the reprex package (v2.0.1)}

This is an interesting question but might be really hard. I think what you really want to be able to do is to draw *screen-aligned quads* of a specified size. *Sprites* are screen-aligned squares, but I haven't figured out yet how to manipulate their dimensions (by adjusting `userMatrix()`??) `close3d(); open3d(); sprites3d(0,0,0, lit = FALSE, alpha = 1, color = "blue", userMatrix = 4*diag(4)); spheres3d(1,1,1)` — Ben Bolker, Mar 06 '22 at 21:54
I think you should be drawing cylinders between the spheres, not lines. — user2554330, Mar 06 '22 at 23:12
Yes, resorting to cylinders may be an idea, but I need to think about how to do that. I will simplify the reprex and ask for help with that if I do not succeed. — larry77, Mar 07 '22 at 08:00

user2554330 · Accepted Answer · 2022-03-08T22:57:28.277

The previous answer https://stackoverflow.com/a/71379091/2554330 comes very close to solving the problem, but there are some minor issues:

Some of the links between the spheres are somewhat flat, because specifying the e2 argument in cylinder3d means the rotationally symmetric cross section is not perpendicular to the cylinder. Leaving it out fixes this.
You can see the facets on the cylinders (which are 6 sided by default). Since these are supposed to be interpreted as lines which resize with the scene, suppressing the lighting using the lit = FALSE material property makes them look more like fat lines.

The sphere1.f function has a noticeable seam in it where the edges of the curved surface join, because persp3d estimates the normals using interior points. Specifying the normals explicitly fixes this. They are specified with a function like f, but giving unit normals to the surface, i.e.

sphere1.f <- function(x0 = 0, y0 = 0, z0 = 0, r = 1, n = 101, ...){
  f <- function(s,t){ 
    cbind(   r * cos(t)*cos(s) + x0,
             r *        sin(s) + y0,
             r * sin(t)*cos(s) + z0)
  }
  g <- function(s,t){
    cbind( cos(t)*cos(s),
           sin(s),
           sin(t)*cos(s))
  }
  persp3d(f, slim = c(-pi/2,pi/2), tlim = c(0, 2*pi), n = n, add = T, normal = g, ...)
}

Each of the spheres drawn by sphere1.f has 101^2 vertices. rgl can handle this, but it is fairly inefficient. Since they are all identical, the sprites3d function can be used to replicate a single sphere at all the different locations. The appropriate code to do this would be
```
## plot the spheres
agg %>%
  rowwise() %>% 
  mutate(x = X1, y = X2, z = X3) %>% 
  sprites3d(shapes = sphere1.f(r = 0.5))
```

where a single sphere centred at (0, 0, 0) is redrawn at each of the computed locations. This looks the same as the original in R, but will make the output from rglwidget() much smaller. (I notice there seems to be a bug in the lighting code, so the shading looks wrong in rglwidget() with the specified lights. Commenting out the lighting code fixes it, but that shouldn't be necessary.)

This is really good! Before I mark yours as the real answer (the least I can do), I have one question: the spheres generated with sprites3d look smaller than the original ones, but if I use sprites3d(radius = 1, shapes = sphere1.f(r=0.5)) then everything appears fine. Could it be that otherwise the radius is scaled twice and it turns out as 0.25? — larry77, Mar 08 '22 at 20:17
Yes, my `radius=0.5` was multiplied by your `r=0.5`. Only one should be set. I'll delete mine. — user2554330, Mar 08 '22 at 22:56
@larry77: The bug mentioned in the last paragraph is fixed as of rgl 0.108.26, available on Github. — user2554330, Mar 12 '22 at 18:34

score 0 · Answer 2 · answered Mar 07 '22 at 09:44

Thanks for the valuable suggestions. Resorting to cylinders got the job done. For the setting up the cylinders, I really made a copy and paste of part of the discussion here

https://r-help.stat.math.ethz.narkive.com/9X5yGnh0/r-joining-two-points-in-rgl

library(rgl)
library(tidyverse)


## a function to plot good-looking spheres

sphere1.f <- function(x0 = 0, y0 = 0, z0 = 0, r = 1, n = 101, ...){
  f <- function(s,t){ 
    cbind(   r * cos(t)*cos(s) + x0,
             r *        sin(s) + y0,
             r * sin(t)*cos(s) + z0)
  }
  persp3d(f, slim = c(-pi/2,pi/2), tlim = c(0, 2*pi), n = n, add = T, ...)
}



## a set of 3D coordinates for my spheres

agg <- structure(list(X1 = c(-0.308421860438279, -1.42503395393061, 
1.10667871416591, -0.41759848570565, 0.523721760757519, 0.520653825151111, 
4.54213267745731, 2.96469370222004, 6.32495200153492, 3.78715565912871, 
5.35968114482443), X2 = c(0.183223776337368, 1.69719822686475, 
-0.992839275466541, 2.22182475540691, -0.705817674534376, -2.40358980860811, 
-0.565561169031234, -0.362260309907445, 0.326094711744554, 0.60340188259578, 
-0.00167511540165435), X3 = c(-0.712687792799106, -0.0336746884947792, 
0.0711272759107127, 1.6126544944538, -2.29999319137504, 1.36257390230441, 
-1.52942342176029, -0.316841449239697, -1.69222713171002, 1.23000775530984, 
2.30848424740017)), class = c("spec_tbl_df", "tbl_df", "tbl", 
"data.frame"), row.names = c(NA, -11L), spec = structure(list(
    cols = list(X1 = structure(list(), class = c("collector_double", 
    "collector")), X2 = structure(list(), class = c("collector_double", 
    "collector")), X3 = structure(list(), class = c("collector_double", 
    "collector"))), default = structure(list(), class = c("collector_guess", 
    "collector")), skip = 0), class = "col_spec"))



##coordinares of the segments (bonds) connecting  the spheres

bond_segments <- structure(list(X1 = c(-1.42503395393061, -0.308421860438279, 
1.10667871416591, -0.308421860438279, 0.523721760757519, -0.308421860438279, 
-0.41759848570565, -1.42503395393061, 0.520653825151111, 1.10667871416591, 
2.96469370222004, 1.10667871416591, 2.96469370222004, 4.54213267745731, 
6.32495200153492, 4.54213267745731, 3.78715565912871, 2.96469370222004, 
5.35968114482443, 3.78715565912871), X2 = c(1.69719822686475, 
0.183223776337368, -0.992839275466541, 0.183223776337368, -0.705817674534376, 
0.183223776337368, 2.22182475540691, 1.69719822686475, -2.40358980860811, 
-0.992839275466541, -0.362260309907445, -0.992839275466541, -0.362260309907445, 
-0.565561169031234, 0.326094711744554, -0.565561169031234, 0.60340188259578, 
-0.362260309907445, -0.00167511540165435, 0.60340188259578), 
    X3 = c(-0.0336746884947792, -0.712687792799106, 0.0711272759107127, 
    -0.712687792799106, -2.29999319137504, -0.712687792799106, 
    1.6126544944538, -0.0336746884947792, 1.36257390230441, 0.0711272759107127, 
    -0.316841449239697, 0.0711272759107127, -0.316841449239697, 
    -1.52942342176029, -1.69222713171002, -1.52942342176029, 
    1.23000775530984, -0.316841449239697, 2.30848424740017, 1.23000775530984
    )), class = c("spec_tbl_df", "tbl_df", "tbl", "data.frame"
), row.names = c(NA, -20L), spec = structure(list(cols = list(
    X1 = structure(list(), class = c("collector_double", "collector"
    )), X2 = structure(list(), class = c("collector_double", 
    "collector")), X3 = structure(list(), class = c("collector_double", 
    "collector"))), default = structure(list(), class = c("collector_guess", 
"collector")), skip = 0), class = "col_spec"))


lb <- nrow(bond_segments)/2

open3d()
#> glX 
#>   1

par3d(skipRedraw=TRUE)

##material and light effects for the spheres

material3d(ambient = "black", specular = "grey60", emission = "black", shininess = 30.0)
clear3d(type = "lights")
light3d(theta = -30, phi = 60, viewpoint.rel = TRUE, ambient = "#FFFFFF", diffuse = "#FFFFFF", specular = "#FFFFFF", x = NULL, y = NULL, z = NULL)
light3d(theta = -0, phi = 0, viewpoint.rel = TRUE,  diffuse = "gray20", specular = "gray25", ambient = "gray80", x = NULL, y = NULL, z = NULL)

## plot the spheres
agg %>%
  rowwise() %>%
  mutate(spheres = sphere1.f(X1, X2, X3, r=0.5))
#> # A tibble: 11 × 4
#> # Rowwise: 
#>        X1       X2      X3 spheres   
#>     <dbl>    <dbl>   <dbl> <rglLwlvl>
#>  1 -0.308  0.183   -0.713  15        
#>  2 -1.43   1.70    -0.0337 16        
#>  3  1.11  -0.993    0.0711 17        
#>  4 -0.418  2.22     1.61   18        
#>  5  0.524 -0.706   -2.30   19        
#>  6  0.521 -2.40     1.36   20        
#>  7  4.54  -0.566   -1.53   21        
#>  8  2.96  -0.362   -0.317  22        
#>  9  6.32   0.326   -1.69   23        
#> 10  3.79   0.603    1.23   24        
#> 11  5.36  -0.00168  2.31   25


### This is the bit of the code which I replaced. I now use cylinders instead
## of segments
## add the segments
## segments3d(bond_segments, lwd=8, color="black")


for (i in 1:lb){

pts <- bond_segments[(2*i-1):(2*i),]

## for the details of this, see the discussion at

## https://r-help.stat.math.ethz.narkive.com/9X5yGnh0/r-joining-two-points-in-rgl
    
shade3d(cylinder3d(pts, radius=0.1, e2=rbind(c(1,0,0),c(1,0,0))), col="black")

}


par3d(skipRedraw=F)

sessionInfo()
#> R version 4.1.2 (2021-11-01)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Debian GNU/Linux 11 (bullseye)
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.13.so
#> 
#> locale:
#>  [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB.UTF-8        LC_COLLATE=en_GB.UTF-8    
#>  [5] LC_MONETARY=en_GB.UTF-8    LC_MESSAGES=en_GB.UTF-8   
#>  [7] LC_PAPER=en_GB.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_GB.UTF-8 LC_IDENTIFICATION=C       
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] forcats_0.5.1   stringr_1.4.0   dplyr_1.0.8     purrr_0.3.4    
#>  [5] readr_2.1.0     tidyr_1.1.4     tibble_3.1.6    ggplot2_3.3.5  
#>  [9] tidyverse_1.3.1 rgl_0.107.14   
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.7        lubridate_1.8.0   assertthat_0.2.1  digest_0.6.28    
#>  [5] utf8_1.2.2        R6_2.5.1          cellranger_1.1.0  backports_1.3.0  
#>  [9] reprex_2.0.1      evaluate_0.14     httr_1.4.2        highr_0.9        
#> [13] pillar_1.6.4      rlang_1.0.1       readxl_1.3.1      R.utils_2.11.0   
#> [17] R.oo_1.24.0       rmarkdown_2.11    styler_1.6.2      htmlwidgets_1.5.4
#> [21] munsell_0.5.0     broom_0.7.10      compiler_4.1.2    modelr_0.1.8     
#> [25] xfun_0.28         pkgconfig_2.0.3   htmltools_0.5.2   tidyselect_1.1.1 
#> [29] fansi_0.5.0       crayon_1.4.2      tzdb_0.2.0        dbplyr_2.1.1     
#> [33] withr_2.4.2       R.methodsS3_1.8.1 grid_4.1.2        jsonlite_1.7.2   
#> [37] gtable_0.3.0      lifecycle_1.0.1   DBI_1.1.1         magrittr_2.0.1   
#> [41] scales_1.1.1      cli_3.1.0         stringi_1.7.5     fs_1.5.0         
#> [45] xml2_1.3.2        ellipsis_0.3.2    generics_0.1.1    vctrs_0.3.8      
#> [49] tools_4.1.2       R.cache_0.15.0    glue_1.5.0        hms_1.1.1        
#> [53] fastmap_1.1.0     yaml_2.2.1        colorspace_2.0-2  rvest_1.0.2      
#> [57] knitr_1.36        haven_2.4.3

^{Created on 2022-03-07 by the reprex package (v2.0.1)}

You might want to set `lit = FALSE` on your cylinders so that you can't see the facets. Setting `e2` is probably not a good idea; it's making some of your bonds very flat and thin. Do this using `cyl <- cylinder3d(pts, radius=0.1); cyl$material <- list(col = "black", lit = FALSE); shade3d(cyl)`. — user2554330, Mar 07 '22 at 10:34
One other comment: specifying the normals for your spheres would get rid of the seams you can see around the equators. Using 3D sprites for them would be more efficient, especially if you are planning to export to WebGL. — user2554330, Mar 07 '22 at 10:40
Thanks for improving the cylinders! Much appreciated. Could you please (if it is not too much to ask) provide some code to improve the spheres? — larry77, Mar 08 '22 at 12:15

R+RGL: plot of spheres and segments, segments do not zoom properly

2 Answers2