Add Cone primitive (#257)

* add Cone primitive

* add tests

* fix docs ci

* fix normals

* update test

* clean up some code duplication

* add comment on nvertices choice

* add docs

* fix typos

Author: ffreyer

.github/workflows/ci.yml

@@ -43,7 +43,7 @@ jobs:
           token: ${{ secrets.CODECOV_TOKEN }}
   docs:
     name: Documentation
-    runs-on: ubuntu-20.04
+    runs-on: ubuntu-latest
     env:
       JULIA_PKG_SERVER: ""
     steps:

docs/src/primitives.md

@@ -9,7 +9,7 @@ GeometryBasics comes with a few predefined primitives:
 
 #### HyperRectangle
 
-A `Rect{D, T} = HyperRectangle{D, T}` is a D-dimensional axis-aligned 
+A `Rect{D, T} = HyperRectangle{D, T}` is a D-dimensional axis-aligned
 hyperrectangle defined by an origin and a size.
 
 ```@repl rects
@@ -33,7 +33,7 @@ Shorthands:
 
 #### Sphere and Circle
 
-`Circle` and `Sphere` are the 2 and 3 dimensional variants of `HyperSphere`. 
+`Circle` and `Sphere` are the 2 and 3 dimensional variants of `HyperSphere`.
 They are defined by an origin and a radius.
 While you can technically create a HyperSphere of any dimension, decomposition
 is only defined in 2D and 3D.
@@ -54,23 +54,34 @@ The coordinates of Circle are defined in anti-clockwise order.
 
 A `Cylinder` is a 3D shape defined by two points and a radius.
 
-
 ```@setup cylinder
 using GeometryBasics
 ```
 ```@repl cylinder
 c = Cylinder(Point3f(-1, 0, 0), Point3f(0, 0, 1), 0.3f0) # start point, end point, radius
 ```
 
-Cylinder supports normals an Tessellation, but currently no texture coordinates.
+Cylinder supports normals and Tessellation, but currently no texture coordinates.
+
+#### Cone
+
+A `Cone` is also defined by two points and a radius, but the radius decreases to 0 from the start point to the tip.
+
+```@setup cone
+using GeometryBasics
+```
+```@repl cone
+c = Cone(Point3f(-1, 0, 0), Point3f(0, 0, 1), 0.3f0) # start point, tip point, radius
+```
+
+Cone supports normals and Tessellation, but currently no texture coordinates.
 
 #### Pyramid
 
 `Pyramid` corresponds to a pyramid shape with a square base and four triangles
 coming together into a sharp point.
 It is defined by by the center point of the base, its height and its width.
 
-
 ```@setup pyramid
 using GeometryBasics
 ```
@@ -132,7 +143,7 @@ end
 ```
 
 To connect these points into a mesh, we need to generate a set of faces.
-The faces of a parallelepiped are parallelograms, which we can describe with `QuadFace`. 
+The faces of a parallelepiped are parallelograms, which we can describe with `QuadFace`.
 Here we should be conscious of the winding direction of faces.
 They are often used to determine the front vs the backside of a (2D) face.
 For example GeometryBasics normal generation and OpenGL's backface culling assume a counter-clockwise winding direction to correspond to a front-facing face.
@@ -187,7 +198,7 @@ function GeometryBasics.texturecoordinates(::Parallelepiped{T}) where {T}
     uvs = [Vec2f(x, y) for x in range(0, 1, length=4) for y in range(0, 1, 3)]
     fs = QuadFace{Int}[
         (1, 2, 5, 4),   (2, 3, 6, 5),
-        (4, 5, 8, 7),   (5, 6, 9, 8), 
+        (4, 5, 8, 7),   (5, 6, 9, 8),
         (7, 8, 11, 10), (8, 9, 12, 11)
     ]
     return FaceView(uvs, fs)

src/GeometryBasics.jl

@@ -17,6 +17,7 @@ include("primitives/spheres.jl")
 include("primitives/cylinders.jl")
 include("primitives/pyramids.jl")
 include("primitives/particles.jl")
+include("primitives/Cone.jl")
 
 include("interfaces.jl")
 include("viewtypes.jl")
@@ -57,7 +58,7 @@ export triangle_mesh, triangle_mesh, uv_mesh
 export uv_mesh, normal_mesh, uv_normal_mesh
 
 export height, origin, radius, width, widths
-export HyperSphere, Circle, Sphere
+export HyperSphere, Circle, Sphere, Cone
 export Cylinder, Pyramid, extremity
 export HyperRectangle, Rect, Rect2, Rect3, Recti, Rect2i, Rect3i, Rectf, Rect2f, Rect3f, Rectd, Rect2d, Rect3d, RectT
 export before, during, meets, overlaps, intersects, finishes

src/primitives/Cone.jl

@@ -0,0 +1,105 @@
+"""
+    Cone{T}(origin::Point3, tip::Point3, radius)
+
+A Cone is a cylinder where one end has a radius of 0. It is defined by an
+`origin` with a finite `radius` which linearly decreases to 0 at the `tip`.
+"""
+struct Cone{T} <: GeometryPrimitive{3, T}
+    origin::Point3{T}
+    tip::Point3{T}
+    radius::T
+end
+
+function Cone(origin::Point3{T1}, tip::Point3{T2}, radius::T3) where {T1, T2, T3}
+    T = promote_type(T1, T2, T3)
+    return Cone{T}(origin, tip, radius)
+end
+
+origin(c::Cone) = c.origin
+extremity(c::Cone) = c.tip
+radius(c::Cone) = c.radius
+height(c::Cone) = norm(c.tip - c.origin)
+direction(c::Cone) = (c.tip .- c.origin) ./ height(c)
+
+# Note:
+# nvertices is matched with Cylinder, where each end has half the vertices. That
+# results in less than nvertices for Cone, but allows a Cylinder and a Cone to
+# be seamless matched with the same `nvertices`
+
+function coordinates(c::Cone{T}, nvertices=30) where {T}
+    nvertices += isodd(nvertices)
+    nhalf = div(nvertices, 2)
+
+    R = cylinder_rotation_matrix(direction(c))
+    step = 2pi / nhalf
+
+    ps = Vector{Point3{T}}(undef, nhalf + 2)
+    for i in 1:nhalf
+        phi = (i-1) * step
+        ps[i] = R * Point3{T}(c.radius * cos(phi), c.radius * sin(phi), 0) + c.origin
+    end
+    ps[end-1] = c.tip
+    ps[end] = c.origin
+
+    return ps
+end
+
+function normals(c::Cone, nvertices = 30)
+    nvertices += isodd(nvertices)
+    nhalf = div(nvertices, 2)
+
+    R = cylinder_rotation_matrix(direction(c))
+    step = 2pi / nhalf
+
+    ns = Vector{Vec3f}(undef, nhalf + 2)
+    # shell at origin
+    # normals are angled in z direction due to change in radius (from radius to 0)
+    # This can be calculated from triangles
+    z = radius(c) / height(c)
+    norm = 1.0 / sqrt(1 + z*z)
+    for i in 1:nhalf
+        phi = (i-1) * step
+        ns[i] = R * (norm * Vec3f(cos(phi), sin(phi), z))
+    end
+
+    # tip - this is undefined / should be all ring angles at once
+    # for rendering it is useful to define this as Vec3f(0), because tip normal
+    # has no useful value to contribute to the interpolated fragment normal
+    ns[end-1] = Vec3f(0)
+
+    # cap
+    ns[end] = Vec3f(normalize(c.origin - c.tip))
+
+    faces = Vector{GLTriangleFace}(undef, nvertices)
+
+    # shell
+    for i in 1:nhalf
+        faces[i] = GLTriangleFace(i, mod1(i+1, nhalf), nhalf+1)
+    end
+
+    # cap
+    for i in 1:nhalf
+        faces[i+nhalf] = GLTriangleFace(nhalf + 2)
+    end
+
+    return FaceView(ns, faces)
+end
+
+function faces(::Cone, facets=30)
+    nvertices = facets + isodd(facets)
+    nhalf = div(nvertices, 2)
+
+    faces = Vector{GLTriangleFace}(undef, nvertices)
+
+    # shell
+    for i in 1:nhalf
+        faces[i] = GLTriangleFace(i, mod1(i+1, nhalf), nhalf+1)
+    end
+
+    # cap
+    for i in 1:nhalf
+        faces[i+nhalf] = GLTriangleFace(i, mod1(i+1, nhalf), nhalf+2)
+    end
+
+    return faces
+end

src/primitives/cylinders.jl

@@ -21,8 +21,15 @@ radius(c::Cylinder) = c.r
 height(c::Cylinder) = norm(c.extremity - c.origin)
 direction(c::Cylinder) = (c.extremity .- c.origin) ./ height(c)
 
-function rotation(c::Cylinder{T}) where {T}
-    d3 = direction(c)
+"""
+    cylinder_rotation_matrix(direction::VecTypes{3})
+
+Creates a basis transformation matrix `R` that maps the third dimension to the
+given `direction` and the first and second to orthogonal directions. This allows
+you to encode a rotation around `direction` in the first two components and
+transform it with `R * rotated_point`.
+"""
+function cylinder_rotation_matrix(d3::VecTypes{3, T}) where {T}
     u = Vec{3, T}(d3[1], d3[2], d3[3])
     if abs(u[1]) > 0 || abs(u[2]) > 0
         v = Vec{3, T}(u[2], -u[1], T(0))
@@ -39,9 +46,9 @@ function coordinates(c::Cylinder{T}, nvertices=30) where {T}
     nvertices += isodd(nvertices)
     nhalf = div(nvertices, 2)
 
-    R = rotation(c)
+    R = cylinder_rotation_matrix(direction(c))
     step = 2pi / nhalf
-    
+
     ps = Vector{Point3{T}}(undef, nvertices + 2)
     for i in 1:nhalf
         phi = (i-1) * step
@@ -61,9 +68,9 @@ function normals(c::Cylinder, nvertices = 30)
     nvertices += isodd(nvertices)
     nhalf = div(nvertices, 2)
 
-    R = rotation(c)
+    R = cylinder_rotation_matrix(direction(c))
     step = 2pi / nhalf
-    
+
     ns = Vector{Vec3f}(undef, nhalf + 2)
     for i in 1:nhalf
         phi = (i-1) * step

test/geometrytypes.jl

@@ -696,4 +696,70 @@ end
     @test all(getindex.(Ref(mp), 1:10) .== ps1)
     @test size(mp) == (10, ) # TODO: Does this make sense?
     @test length(mp) == 10
+end
+
+@testset "Cone" begin
+    @testset "constructors" begin
+        v1 = rand(Point{3,Float64})
+        v2 = rand(Point{3,Float64})
+        R = rand()
+        s = Cone(v1, v2, R)
+        @test typeof(s) == Cone{Float64}
+        @test origin(s) == v1
+        @test extremity(s) == v2
+        @test radius(s) == R
+        @test height(s) == norm(v2 - v1)
+        @test isapprox(direction(s), (v2 - v1) ./ norm(v2 .- v1))
+    end
+
+    @testset "decompose" begin
+        v1 = Point{3,Float64}(1, 2, 3)
+        v2 = Point{3,Float64}(4, 5, 6)
+        R = 5.0
+        s = Cone(v1, v2, R)
+        positions = Point{3,Float64}[
+            (4.535533905932738, -1.5355339059327373, 3.0),
+            (3.0412414523193148, 4.041241452319315, -1.0824829046386295),
+            (-2.535533905932737, 5.535533905932738, 2.9999999999999996),
+            (-1.0412414523193152, -0.04124145231931431, 7.0824829046386295),
+            (4, 5, 6),
+            (1, 2, 3)
+        ]
+
+        @test decompose(Point3{Float64}, Tessellation(s, 8)) ≈ positions
+
+        _faces = TriangleFace[
+            (1,2,5), (2,3,5), (3,4,5), (4,1,5),
+            (1,2,6), (2,3,6), (3,4,6), (4,1,6)]
+
+        @test _faces == decompose(TriangleFace{Int}, Tessellation(s, 8))
+
+        m = triangle_mesh(Tessellation(s, 8))
+        @test m === triangle_mesh(m)
+        @test GeometryBasics.faces(m) == decompose(GLTriangleFace, _faces)
+        @test GeometryBasics.coordinates(m) ≈ positions
+
+        m = normal_mesh(s) # just test that it works without explicit resolution parameter
+        @test hasproperty(m, :position)
+        @test hasproperty(m, :normal)
+        @test faces(m) isa AbstractVector{GLTriangleFace}
+
+        ns = Vec{3, Float32}[
+            (0.90984505, -0.10920427, 0.40032038),
+            (0.6944946, 0.6944946, -0.18802801),
+            (-0.10920427, 0.90984505, 0.40032038),
+            (0.106146194, 0.106146194, 0.9886688),
+            (0.0, 0.0, 0.0),
+            (-0.57735026, -0.57735026, -0.57735026),
+        ]
+        fs = [
+            GLTriangleFace(1, 2, 5), GLTriangleFace(2, 3, 5), GLTriangleFace(3, 4, 5), GLTriangleFace(4, 1, 5),
+            GLTriangleFace(6, 6, 6), GLTriangleFace(6, 6, 6), GLTriangleFace(6, 6, 6), GLTriangleFace(6, 6, 6)
+        ]
+
+        @test FaceView(ns, fs) == decompose_normals(Tessellation(s, 8))
+
+        muv = uv_mesh(s)
+        @test !hasproperty(muv, :uv) # not defined yet
+    end
 end