title: R + plotly: solid of revolution
body: <p>I have a function <code>r(x)</code> that I want to rotate around the <code>x</code> axis to get a <a href="https://en.wikipedia.org/wiki/Solid_of_revolution" rel="nofollow noreferrer">solid of revolution</a> that I want to add to an existing <code>plot_ly</code> plot using <code>add_surface</code> (colored by <code>x</code>).</p>
  <p>Here is an example:</p>
  <pre class="lang-r prettyprint-override"><code>library(dplyr)
  library(plotly)
  
  # radius depends on x
  r &lt;- function(x) x^2
  
  # interval of interest
  int &lt;- c(1, 3)
  
  # number of points along the x-axis
  nx &lt;- 20
  
  # number of points along the rotation
  ntheta &lt;- 36
  
  # set x points and get corresponding radii
  coords &lt;- data_frame(x = seq(int[1], int[2], length.out = nx), r = r(x))
  
  # for each x: rotate r to get y and z coordinates
  # edit: ensure 0 and pi are both amongst the angles used
  coords %&lt;&gt;%
    rowwise() %&gt;%
    do(data_frame(x = .$x, r = .$r,
                  theta = seq(0, pi, length.out = ntheta / 2 + 1) %&gt;%
                  c(pi + .[-c(1, length(.))]))) %&gt;%
  
    ungroup %&gt;%
    mutate(y = r * cos(theta), z = r * sin(theta))
  
  # plot points to make sure the coordinates define the desired shape
  coords %&gt;%
    plot_ly(x = ~x, y = ~y, z = ~z, color = ~x) %&gt;%
    add_markers()
  </code></pre>
  <p><a href="https://i.stack.imgur.com/BMu6j.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/BMu6j.png" alt="3D scatter plot" /></a></p>
  <h1>How can I generate the shape indicated by the points above as a <code>plotly</code> surface (ideally open on both ends)?</h1>
  <hr />
  <p><strong>edit (1)</strong>:</p>
  <p>Here is my best attempt so far:</p>
  <pre class="lang-r prettyprint-override"><code># get all x &amp; y values used (sort to connect halves on the side)
  xs &lt;-
    unique(coords$x) %&gt;%
    sort
  ys &lt;-
    unique(coords$y) %&gt;%
    sort
  
  # for each possible x/y pair: get z^2 value
  coords &lt;-
    expand.grid(x = xs, y = ys) %&gt;%
    as_data_frame %&gt;%
    mutate(r = r(x), z2 = r^2 - y^2)
  
  # format z coordinates above x/y plane as matrix where columns
  # represent x and rows y
  zs &lt;- matrix(sqrt(coords$z2), ncol = length(xs), byrow = TRUE)
  
  # format x coordiantes as matrix as above (for color gradient)
  gradient &lt;-
    rep(xs, length(ys)) %&gt;%
    matrix(ncol = length(xs), byrow = TRUE)
    
  # plot upper half of shape as surface
  p &lt;- plot_ly(x = xs, y = ys, z = zs, surfacecolor = gradient,
               type = &quot;surface&quot;, colorbar = list(title = 'x'))
  
  # plot lower have of shape as second surface
  p %&gt;%
    add_surface(z = -zs, showscale = FALSE)
  </code></pre>
  <p><a href="https://i.stack.imgur.com/j7t11.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/j7t11.png" alt="two 3D surfaces attempt" /></a></p>
  <p>While this gives the desired shape,</p>
  <ol>
  <li>It has 'razor teeth' close to the <code>x</code>/<code>y</code> plane.</li>
  <li><s>The halves parts don't touch.</s> (<em>resolved by including <code>0</code> and <code>pi</code> in the <code>theta</code> vectors</em>)</li>
  <li><s>I didn't figure out how to color it by <code>x</code> instead of <code>z</code> (though I didn't look much into this so far).</s> (<em>resolved by <code>gradient</code> matrix</em>)</li>
  </ol>
  <hr />
  <p><strong>edit (2)</strong>:</p>
  <p>Here is an attempt using a single surface:</p>
  <pre class="lang-r prettyprint-override"><code># close circle in y-direction
  ys &lt;- c(ys, rev(ys), ys[1])
  
  # get corresponding z-values
  zs &lt;- rbind(zs, -zs[nrow(zs):1, ], zs[1, ])
  
  # as above, but for color gradient
  gradient &lt;-
    rbind(gradient, gradient[nrow(gradient):1, ], gradient[1, ])
  
  # plot single surface
  plot_ly(x = xs, y = ys, z = zs, surfacecolor = gradient,
          type = &quot;surface&quot;, colorbar = list(title = 'x'))
  </code></pre>
  <p>Surprisingly, while this should connect the two halves orthogonal to the <code>x</code> / <code>y</code> plane to create the full shape,
  it still suffers from the same 'razor teeth' effect as the above solution:</p>
  <p><a href="https://i.stack.imgur.com/CJ0RH.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/CJ0RH.png" alt="single 3D surface attempt" /></a></p>
  <hr />
  <p><strong>edit (3)</strong>:</p>
  <p>It turns out the missing parts result from <code>z</code>-values being <code>NaN</code> when close to 0:</p>
  <pre class="lang-r prettyprint-override"><code># color points 'outside' the solid purple
  gradient[is.nan(zs)] &lt;- -1
  
  # show those previously hidden points
  zs[is.nan(zs)] &lt;- 0
  
  # plot exactly as before
  plot_ly(x = xs, y = ys, z = zs, surfacecolor = gradient,
          type = &quot;surface&quot;, colorbar = list(title = 'x'))
  </code></pre>
  <p><a href="https://i.stack.imgur.com/3pDTO.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/3pDTO.png" alt="NaN-override attempt" /></a></p>
  <p><s>This could be caused by numerical instability of the substraction when <code>r^2</code> and <code>y</code> get too close, resulting in negative input for <code>sqrt</code> where the actual input is still non-negative.</s></p>
  <p>This seams unrelated to numerical issues as even when considering +-4 'close' to zero, the 'razor teeth' effect can not be avoided completely:</p>
  <pre class="lang-r prettyprint-override"><code># re-calculate z-values rounding to zero if 'close'
  eps &lt;- 4
  zs &lt;- with(coords, ifelse(abs(z2) &lt; eps, 0, sqrt(z2))) %&gt;%
        matrix(ncol = length(xs), byrow = TRUE) %&gt;%
        rbind(-.[nrow(.):1, ], .[1, ])
  
  # plot exactly as before
  plot_ly(x = xs, y = ys, z = zs, surfacecolor = gradient,
          type = &quot;surface&quot;, colorbar = list(title = 'x'))
  </code></pre>
  <p><a href="https://i.stack.imgur.com/xn3Rj.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/xn3Rj.png" alt="eps-attempt" /></a></p>