diff --git a/src/GeometryBasics.jl b/src/GeometryBasics.jl
index 6e69d786..dbb6c1a8 100644
--- a/src/GeometryBasics.jl
+++ b/src/GeometryBasics.jl
@@ -5,9 +5,14 @@ using EarCut_jll
 
 using Base: @propagate_inbounds
 
-include("fixed_arrays.jl")
+import Base: +, -
+
+# basic concepts
 include("vectors.jl")
 include("matrices.jl")
+include("points.jl")
+
+include("fixed_arrays.jl")
 include("offsetintegers.jl")
 include("basic_types.jl")
 
@@ -26,12 +31,15 @@ include("lines.jl")
 include("boundingboxes.jl")
 
 # points
-export AbstractPoint, Point, PointMeta, PointWithUV
+export Point, Point2, Point3, Point2f, Point3f
 
 # vectors
 export Vec, Vec2, Vec3, Vec2f, Vec3f
 export vunit, vfill
 
+# TODO: review these
+export AbstractPoint, PointMeta, PointWithUV
+
 # geometries
 export AbstractGeometry, GeometryPrimitive
 export LineFace, Polytope, Line, NgonFace, convert_simplex
diff --git a/src/basic_types.jl b/src/basic_types.jl
index aac0d153..7eeb7814 100644
--- a/src/basic_types.jl
+++ b/src/basic_types.jl
@@ -13,7 +13,6 @@ Note That `Polytope{N} where N == 3` denotes a Triangle both as a Simplex or Ngo
 abstract type Polytope{Dim,T} <: AbstractGeometry{Dim,T} end
 abstract type AbstractPolygon{Dim,T} <: Polytope{Dim,T} end
 
-abstract type AbstractPoint{Dim,T} <: StaticVector{Dim,T} end
 abstract type AbstractFace{N,T} <: StaticVector{N,T} end
 abstract type AbstractSimplexFace{N,T} <: AbstractFace{N,T} end
 abstract type AbstractNgonFace{N,T} <: AbstractFace{N,T} end
@@ -57,7 +56,6 @@ Fixed Size Polygon, e.g.
 - ...
 """
 struct Ngon{Dim,T<:Real,N,Point<:AbstractPoint{Dim,T}} <: AbstractPolygon{Dim,T}
-
     points::SVector{N,Point}
 end
 
diff --git a/src/boundingboxes.jl b/src/boundingboxes.jl
index 2cd4b478..4fe736c2 100644
--- a/src/boundingboxes.jl
+++ b/src/boundingboxes.jl
@@ -8,10 +8,10 @@ boundingbox(geom) = boundingbox(coordinates(geom))
 # fallback implementation treats geometry as
 # a set of points (i.e. coordinates)
 function boundingbox(geometry::AbstractArray{<:AbstractPoint{N,T}}) where {N,T}
-    vmin = Point{N,T}(typemax(T))
-    vmax = Point{N,T}(typemin(T))
+    vmin = vfill(Vec{N,T}, typemax(T))
+    vmax = vfill(Vec{N,T}, typemin(T))
     for p in geometry
-        vmin, vmax = minmax(p, vmin, vmax)
+        vmin, vmax = minmax(coordinates(p), vmin, vmax)
     end
     Rect{N,T}(vmin, vmax - vmin)
 end
diff --git a/src/fixed_arrays.jl b/src/fixed_arrays.jl
index 632301b7..ef4c5062 100644
--- a/src/fixed_arrays.jl
+++ b/src/fixed_arrays.jl
@@ -106,23 +106,3 @@ macro fixed_vector(name, parent)
     end
     return esc(expr)
 end
-
-abstract type AbstractPoint{Dim,T} <: StaticVector{Dim,T} end
-@fixed_vector Point AbstractPoint
-
-const Pointf0{N} = Point{N,Float32}
-
-for i in 1:4
-    for T in [:Point]
-        name = Symbol("$T$i")
-        namef0 = Symbol("$T$(i)f0")
-        @eval begin
-            const $name = $T{$i}
-            const $namef0 = $T{$i,Float32}
-            export $name
-            export $namef0
-        end
-    end
-end
-
-export Pointf0
diff --git a/src/geometry_primitives.jl b/src/geometry_primitives.jl
index b4b1248f..310bde79 100644
--- a/src/geometry_primitives.jl
+++ b/src/geometry_primitives.jl
@@ -55,37 +55,18 @@ Extract all line segments in a Face.
     return v
 end
 
+# TODO: review these
 to_pointn(::Type{T}, x) where {T<:Point} = convert_simplex(T, x)[1]
 
-# disambiguation method overlords
-convert_simplex(::Type{Point}, x::Point) = (x,)
-convert_simplex(::Type{Point{N,T}}, p::Point{N,T}) where {N,T} = (p,)
-function convert_simplex(::Type{Point{N,T}}, x) where {N,T}
-    N2 = length(x)
-    return (Point{N,T}(ntuple(i -> i <= N2 ? T(x[i]) : T(0), N)),)
+# TODO: why increase the dimension of the point?
+function convert_simplex(::Type{Point{N,T}}, p::Point{M,V}) where {N,T,M,V}
+    x = coordinates(p)
+    return (Point(ntuple(i -> i <= M ? T(x[i]) : T(0), N)),)
 end
 
-function convert_simplex(::Type{Vec{N,T}}, x) where {N,T}
-    N2 = length(x)
-    return (Vec{N,T}(ntuple(i -> i <= N2 ? T(x[i]) : T(0), N)),)
-end
-
-collect_with_eltype(::Type{T}, vec::Vector{T}) where {T} = vec
-collect_with_eltype(::Type{T}, vec::AbstractVector{T}) where {T} = collect(vec)
-
-function collect_with_eltype(::Type{T}, iter) where {T}
-    # TODO we could be super smart about allocating the right length
-    # but its kinda annoying, since e.g. T == Triangle and first(iter) isa Quad
-    # will need double the length etc - but could all be figured out ;)
-    result = T[]
-    for element in iter
-        # convert_simplex always returns a tuple,
-        # so that e.g. convert(Triangle, quad) can return 2 elements
-        for telement in convert_simplex(T, element)
-            push!(result, telement)
-        end
-    end
-    return result
+# TODO: review these
+function convert_simplex(::Type{Vec{N,T}}, v::Vec{M,V}) where {N,T,M,V}
+    return (Vec(ntuple(i -> i <= M ? T(v[i]) : T(0), N)),)
 end
 
 """
diff --git a/src/interfaces.jl b/src/interfaces.jl
index a12bbc4d..2bf89235 100644
--- a/src/interfaces.jl
+++ b/src/interfaces.jl
@@ -1,8 +1,7 @@
+# TODO: review these
 """
     coordinates(geometry)
-Returns the edges/vertices/coordinates of a geometry. Is allowed to return lazy iterators!
-Use `decompose(ConcretePointType, geometry)` to get `Vector{ConcretePointType}` with
-`ConcretePointType` to be something like `Point{3, Float32}`.
+Returns the vertices of a geometry.
 """
 function coordinates(points::AbstractVector{<:AbstractPoint})
     return points
@@ -43,7 +42,7 @@ To transport this information to the various decompose methods, you can wrap it
 in the Tesselation object e.g. like this:
 
 ```julia
-sphere = Sphere(Point3f0(0), 1)
+sphere = Sphere(Point3f(0,0,0), 1.0f0)
 m1 = mesh(sphere) # uses a default value for tesselation
 m2 = mesh(Tesselation(sphere, 64)) # uses 64 for tesselation
 length(coordinates(m1)) != length(coordinates(m2))
@@ -80,46 +79,56 @@ function texturecoordinates(tesselation::Tesselation)
     return texturecoordinates(tesselation.primitive, nvertices(tesselation))
 end
 
-## Decompose methods
-# Dispatch type to make `decompose(UV{Vec2f}, primitive)` work
-# and to pass through tesselation information
-
 # Types that can be converted to a mesh via the functions below
-const Meshable{Dim,T} = Union{Tesselation{Dim,T},Mesh{Dim,T},AbstractPolygon{Dim,T},
-                              GeometryPrimitive{Dim,T},
-                              AbstractVector{<:AbstractPoint{Dim,T}}}
+const Meshable{Dim,T} = Union{Mesh{Dim,T},
+                              Tesselation{Dim,T},
+                              AbstractPolygon{Dim,T},
+                              GeometryPrimitive{Dim,T}}
 
-struct UV{T} end
-UV(::Type{T}) where {T} = UV{T}()
-UV() = UV(Vec2f)
-struct UVW{T} end
-UVW(::Type{T}) where {T} = UVW{T}()
-UVW() = UVW(Vec3f)
-struct Normal{T} end
-Normal(::Type{T}) where {T} = Normal{T}()
-Normal() = Normal(Vec3f)
+"""
+    decompose(T, meshable)
 
-function decompose(::Type{F}, primitive) where {F<:AbstractFace}
-    f = faces(primitive)
-    f === nothing && return nothing
-    return collect_with_eltype(F, f)
+Decompose a `meshable` object (e.g. Polygon) into elements of type `T`
+
+## Example
+
+```julia
+decompose(Point3, Rect3D())
+```
+"""
+function decompose(::Type{T}, primitive) where {T}
+    return collect_with_eltype(T, primitive)
 end
 
+# Specializations
+
 function decompose(::Type{P}, primitive) where {P<:AbstractPoint}
-    return collect_with_eltype(P, metafree(coordinates(primitive)))
+    convert.(P, metafree(coordinates(primitive)))
 end
 
-function decompose(::Type{Point}, primitive::Meshable{Dim,T}) where {Dim,T}
-    return collect_with_eltype(Point{Dim,T}, metafree(coordinates(primitive)))
+function decompose(::Type{F}, primitive) where {F<:AbstractFace}
+    f = faces(primitive)
+    f === nothing && return nothing
+    return collect_with_eltype(F, f)
 end
 
+# TODO: review these
 function decompose(::Type{Point}, primitive::LineString{Dim,T}) where {Dim,T}
     return collect_with_eltype(Point{Dim,T}, metafree(coordinates(primitive)))
 end
 
-function decompose(::Type{T}, primitive) where {T}
-    return collect_with_eltype(T, primitive)
-end
+# TODO: review these
+struct UV{T} end
+UV(::Type{T}) where {T} = UV{T}()
+UV() = UV(Vec2f)
+
+struct UVW{T} end
+UVW(::Type{T}) where {T} = UVW{T}()
+UVW() = UVW(Vec3f)
+
+struct Normal{T} end
+Normal(::Type{T}) where {T} = Normal{T}()
+Normal() = Normal(Vec3f)
 
 decompose_uv(primitive) = decompose(UV(), primitive)
 decompose_uvw(primitive) = decompose(UVW(), primitive)
@@ -133,14 +142,14 @@ function decompose(NT::Normal{T}, primitive) where {T}
     return collect_with_eltype(T, n)
 end
 
-function decompose(UVT::Union{UV{T},UVW{T}}, primitive) where {T}
+function decompose(UVT::Union{UV{T},UVW{T}}, primitive::Meshable{Dim,V}) where {Dim,T,V}
     # This is the fallback for texture coordinates if a primitive doesn't overload them
     # We just take the positions and normalize them
     uv = texturecoordinates(primitive)
     if uv === nothing
         # If the primitive doesn't even have coordinates, we're out of options and return
         # nothing, indicating that texturecoordinates aren't implemented
-        positions = decompose(Point, primitive)
+        positions = decompose(Point{Dim,V}, primitive)
         positions === nothing && return nothing
         # Let this overlord do the work
         return decompose(UVT, positions)
@@ -148,17 +157,20 @@ function decompose(UVT::Union{UV{T},UVW{T}}, primitive) where {T}
     return collect_with_eltype(T, uv)
 end
 
-function decompose(UVT::Union{UV{T},UVW{T}},
-                   positions::AbstractVector{<:Point}) where {T}
+function decompose(UVT::Union{UV{T},UVW{T}}, positions::AbstractVector{<:AbstractPoint}) where {T}
     N = length(T)
-    bb = boundingbox(positions) # Make sure we get this as points
+    bb = boundingbox(positions)
     return map(positions) do p
-        return (p .- minimum(bb)) ./ widths(bb)
+        return (coordinates(p) - minimum(bb)) ./ widths(bb)
     end
 end
 
-# Stay backward compatible:
-
-function decompose(::Type{T}, primitive::Meshable, nvertices) where {T}
-    return decompose(T, Tesselation(primitive, nvertices))
+function collect_with_eltype(::Type{T}, iter) where {T}
+    result = T[]
+    for element in iter
+        for telement in convert_simplex(T, element)
+            push!(result, telement)
+        end
+    end
+    return result
 end
diff --git a/src/lines.jl b/src/lines.jl
index 260826c9..d06a6748 100644
--- a/src/lines.jl
+++ b/src/lines.jl
@@ -7,10 +7,10 @@ Returns intersection_found::Bool, intersection_point
 """
 function intersects(a::Line{2,T1}, b::Line{2,T2}) where {T1,T2}
     T = promote_type(T1, T2)
-    v1, v2 = a
-    v3, v4 = b
+    v1, v2 = coordinates.(a)
+    v3, v4 = coordinates.(b)
     MT = Mat{2,2,T,4}
-    p0 = zero(Point2{T})
+    p0 = Point{2,T}(0, 0)
 
     verticalA = v1[1] == v2[1]
     verticalB = v3[1] == v4[1]
@@ -55,10 +55,10 @@ function intersects(a::Line{2,T1}, b::Line{2,T2}) where {T1,T2}
     (y < prevfloat(min(v3[2], v4[2])) || y > nextfloat(max(v3[2], v4[2]))) &&
         return false, p0
 
-    point = Point2{T}(x, y)
+    point = Point{2,T}(x, y)
     # don't forget to rotate the answer back
     if dorotation
-        point = transpose(rotation) * point
+        point = Point(transpose(rotation) * coordinates(point))
     end
 
     return true, point
diff --git a/src/meshes.jl b/src/meshes.jl
index 828a0cb1..b30577db 100644
--- a/src/meshes.jl
+++ b/src/meshes.jl
@@ -87,8 +87,8 @@ const GLNormalUVWMesh2D = GLNormalUVWMesh{2}
 const GLNormalUVWMesh3D = GLNormalUVWMesh{3}
 
 """
-    mesh(primitive::GeometryPrimitive;
-         pointtype=Point, facetype=GLTriangle,
+    mesh(primitive::Meshable{N,T};
+         pointtype=Point{N,T}, facetype=GLTriangle,
          uvtype=nothing, normaltype=nothing)
 
 Creates a mesh from `primitive`.
@@ -98,8 +98,8 @@ Note, that this can be an `Int` or `Tuple{Int, Int}``, when the primitive is gri
 It also only losely correlates to the number of vertices, depending on the algorithm used.
 #TODO: find a better number here!
 """
-function mesh(primitive::Meshable; pointtype=Point, facetype=GLTriangleFace, uv=nothing,
-              normaltype=nothing)
+function mesh(primitive::Meshable{N,T}; pointtype=Point{N,T}, facetype=GLTriangleFace,
+              uv=nothing, normaltype=nothing) where {N,T}
 
     positions = decompose(pointtype, primitive)
     faces = decompose(facetype, primitive)
@@ -124,7 +124,7 @@ function mesh(primitive::Meshable; pointtype=Point, facetype=GLTriangleFace, uv=
     if normaltype !== nothing
         primitive_normals = normals(primitive)
         if primitive_normals !== nothing
-            attrs[:normals] = decompose(normaltype, primitive_normals)
+            attrs[:normals] = convert.(normaltype, primitive_normals)
         else
             # Normals not implemented for primitive, so we calculate them!
             n = normals(positions, faces; normaltype=normaltype)
@@ -161,26 +161,29 @@ function mesh(polygon::AbstractPolygon{Dim,T}; pointtype=Point{Dim,T},
     return Mesh(positions, faces)
 end
 
-function triangle_mesh(primitive::Meshable{N}; nvertices=nothing) where {N}
+function triangle_mesh(primitive::Meshable{N,T}; nvertices=nothing) where {N,T}
     if nvertices !== nothing
         @warn("nvertices argument deprecated. Wrap primitive in `Tesselation(primitive, nvertices)`")
         primitive = Tesselation(primitive, nvertices)
     end
-    return mesh(primitive; pointtype=Point{N,Float32}, facetype=GLTriangleFace)
+    return mesh(primitive; pointtype=Point{N,T}, facetype=GLTriangleFace)
+end
+
+function triangle_mesh(points::AbstractVector{P}; nvertices=nothing) where {P<:AbstractPoint}
+    triangle_mesh(Polygon(points), nvertices=nvertices)
 end
 
 function uv_mesh(primitive::Meshable{N,T}) where {N,T}
-    return mesh(primitive; pointtype=Point{N,Float32}, uv=Vec2f, facetype=GLTriangleFace)
+    mesh(primitive; pointtype=Point{N,T}, uv=Vec{2,T}, facetype=GLTriangleFace)
 end
 
-function uv_normal_mesh(primitive::Meshable{N}) where {N}
-    return mesh(primitive; pointtype=Point{N,Float32}, uv=Vec2f, normaltype=Vec3f,
-                facetype=GLTriangleFace)
+function uv_normal_mesh(primitive::Meshable{N,T}) where {N,T}
+    mesh(primitive; pointtype=Point{N,T}, uv=Vec{2,T}, normaltype=Vec{3,T}, facetype=GLTriangleFace)
 end
 
 function normal_mesh(points::AbstractVector{<:AbstractPoint},
                      faces::AbstractVector{<:AbstractFace})
-    _points = decompose(Point3f0, points)
+    _points = decompose(Point3f, points)
     _faces = decompose(GLTriangleFace, faces)
     return Mesh(meta(_points; normals=normals(_points, _faces)), _faces)
 end
@@ -190,8 +193,7 @@ function normal_mesh(primitive::Meshable{N}; nvertices=nothing) where {N}
         @warn("nvertices argument deprecated. Wrap primitive in `Tesselation(primitive, nvertices)`")
         primitive = Tesselation(primitive, nvertices)
     end
-    return mesh(primitive; pointtype=Point{N,Float32}, normaltype=Vec3f,
-                facetype=GLTriangleFace)
+    return mesh(primitive; pointtype=Point{N,Float32}, normaltype=Vec3f, facetype=GLTriangleFace)
 end
 
 """
@@ -201,7 +203,7 @@ Calculate the signed volume of one tetrahedron. Be sure the orientation of your
 surface is right.
 """
 function volume(triangle::Triangle) where {VT,FT}
-    v1, v2, v3 = triangle
+    v1, v2, v3 = coordinates.(triangle)
     sig = sign(orthogonal_vector(v1, v2, v3) ⋅ v1)
     return sig * abs(v1 ⋅ (v2 × v3)) / 6
 end
diff --git a/src/points.jl b/src/points.jl
new file mode 100644
index 00000000..40e421d1
--- /dev/null
+++ b/src/points.jl
@@ -0,0 +1,68 @@
+abstract type AbstractPoint{N,T} end
+
+"""
+    Point{N,T}
+
+A point in `N`-dimensional space with coordinates of type `T`.
+The coordinates of the point provided upon construction are with
+respect to the canonical Euclidean basis. See [`vunit`](@ref).
+
+## Example
+
+```julia
+O = Point(0.0, 0.0) # origin of 2D Euclidean space
+```
+
+### Notes
+
+- Type aliases are `Point2`, `Point3`, `Point2f`, `Point3f`
+"""
+struct Point{N,T} <: AbstractPoint{N,T}
+    coords::SVector{N,T}
+end
+
+# convenience constructors
+Point(coords::NTuple{N,T}) where {N,T} = Point{N,T}(SVector(coords))
+Point(coords::Vararg{T,N}) where {N,T} = Point{N,T}(SVector(coords))
+
+# coordinate type conversions
+Point{N,T}(coords::NTuple{N,V}) where {N,T,V} = Point(T.(coords))
+Point{N,T}(coords::Vararg{V,N}) where {N,T,V} = Point(T.(coords))
+Base.convert(::Type{Point{N,T}}, coords) where {N,T} = Point{N,T}(coords)
+Base.convert(::Type{Point{N,T}}, p::Point) where {N,T} = Point{N,T}(p.coords)
+
+# type aliases for convenience
+const Point2  = Point{2,Float64}
+const Point3  = Point{3,Float64}
+const Point2f = Point{2,Float32}
+const Point3f = Point{3,Float32}
+
+"""
+    coordinates(A::Point)
+
+Return the coordinates of the point with respect to the
+canonical Euclidean basis. See [`vunit`](@ref).
+"""
+coordinates(A::Point) = A.coords
+
+"""
+    -(A::Point, B::Point)
+
+Return the [`Vec`](@ref) associated with the direction
+from point `A` to point `B`.
+"""
+-(A::Point, B::Point) = Vec(A.coords - B.coords)
+
+"""
+    +(A::Point, v::Vec)
+    +(v::Vec, A::Point)
+
+Return the point at the end of the vector `v` placed
+at a reference (or start) point `A`.
+"""
++(A::Point, v::Vec) = Point(A.coords + v)
++(v::Vec, A::Point) = A + v
+
+# TODO: implement rand properly with RNG, etc.
+Base.rand(::Type{Point{N,T}}) where {N,T} = Point(rand(SVector{N,T}))
+Base.rand(::Type{Point{N,T}}, n::Integer) where {N,T} = Point.(rand(SVector{N,T}, n))
diff --git a/src/primitives/cylinders.jl b/src/primitives/cylinders.jl
index 9804cf53..262d560a 100644
--- a/src/primitives/cylinders.jl
+++ b/src/primitives/cylinders.jl
@@ -1,13 +1,13 @@
 """
     Cylinder{N, T}
 
-A `Cylinder` is a 2D rectangle or a 3D cylinder defined by its origin point,
-its extremity and a radius. `origin`, `extremity` and `r`, must be specified.
+A `Cylinder` is a 2D rectangle or a 3D cylinder defined by
+its `origin` point, its `extremity` and a `radius`.
 """
 struct Cylinder{N,T} <: GeometryPrimitive{N,T}
     origin::Point{N,T}
     extremity::Point{N,T}
-    r::T
+    radius::T
 end
 
 """
@@ -15,16 +15,16 @@ end
     Cylinder3{T}
 
 A `Cylinder2` or `Cylinder3` is a 2D/3D cylinder defined by its origin point,
-its extremity and a radius. `origin`, `extremity` and `r`, must be specified.
+its extremity and a radius. `origin`, `extremity` and `radius`, must be specified.
 """
 const Cylinder2{T} = Cylinder{2,T}
 const Cylinder3{T} = Cylinder{3,T}
 
 origin(c::Cylinder{N,T}) where {N,T} = c.origin
 extremity(c::Cylinder{N,T}) where {N,T} = c.extremity
-radius(c::Cylinder{N,T}) where {N,T} = c.r
+radius(c::Cylinder{N,T}) where {N,T} = c.radius
 height(c::Cylinder{N,T}) where {N,T} = norm(c.extremity - c.origin)
-direction(c::Cylinder{N,T}) where {N,T} = (c.extremity .- c.origin) ./ height(c)
+direction(c::Cylinder{N,T}) where {N,T} = (c.extremity - c.origin) ./ height(c)
 
 function rotation(c::Cylinder{2,T}) where {T}
     d2 = direction(c)
@@ -49,11 +49,14 @@ function rotation(c::Cylinder{3,T}) where {T}
 end
 
 function coordinates(c::Cylinder{2,T}, nvertices=(2, 2)) where {T}
-    r = Rect(c.origin[1] - c.r / 2, c.origin[2], c.r, height(c))
+    o = coordinates(c.origin)
+    r = c.radius
+    h = height(c)
+    rect = Rect(o[1] - r / 2, o[2], r, h)
     M = rotation(c)
-    points = coordinates(r, nvertices)
-    vo = to_pointn(Point3{T}, origin(c))
-    return (M * (to_pointn(Point3{T}, point) .- vo) .+ vo for point in points)
+    points = Point.(coordinates(rect, nvertices))
+    vo = to_pointn(Point{3,T}, origin(c))
+    return (M * (to_pointn(Point{3,T}, point) - vo) + coordinates(vo) for point in points)
 end
 
 function faces(sphere::Cylinder{2}, nvertices=(2, 2))
@@ -72,13 +75,15 @@ function coordinates(c::Cylinder{3,T}, nvertices=30) where {T}
     range = 1:(2 * nbv + 2)
     function inner(i)
         return if i == length(range)
-            c.extremity
+            coordinates(c.extremity)
         elseif i == length(range) - 1
-            origin(c)
+            coordinates(c.origin)
         else
             phi = T((2π * (((i + 1) ÷ 2) - 1)) / nbv)
-            up = ifelse(isodd(i), 0, h)
-            (M * Point(c.r * cos(phi), c.r * sin(phi), up)) .+ c.origin
+            up = ifelse(isodd(i), T(0), h)
+            o  = coordinates(c.origin)
+            r  = c.radius
+            (M * Vec(r*cos(phi), r*sin(phi), up)) + o
         end
     end
 
diff --git a/src/primitives/rectangles.jl b/src/primitives/rectangles.jl
index 619beba0..f8f6d509 100644
--- a/src/primitives/rectangles.jl
+++ b/src/primitives/rectangles.jl
@@ -7,8 +7,8 @@ Formally it is the cartesian product of intervals, which is represented by the
 `origin` and `width` fields, whose indices correspond to each of the `N` axes.
 """
 struct HyperRectangle{N,T} <: GeometryPrimitive{N,T}
-    origin::Vec{N,T}
-    widths::Vec{N,T}
+    origin::Point{N,T}
+    widths::SVector{N,T}
 end
 
 ##
@@ -155,14 +155,15 @@ function FRect3D(x::Tuple{Tuple{<:Number,<:Number,<:Number},
     return FRect3D(Vec3f0(x[1]...), Vec3f0(x[2]...))
 end
 
-origin(prim::Rect) = prim.origin
-Base.maximum(prim::Rect) = origin(prim) + widths(prim)
-Base.minimum(prim::Rect) = origin(prim)
-Base.length(prim::Rect{N,T}) where {T,N} = N
-widths(prim::Rect) = prim.widths
+# TODO: review these
+origin(r::Rect) = r.origin
+Base.minimum(r::Rect) = coordinates(r.origin)
+Base.maximum(r::Rect) = coordinates(r.origin) + r.widths
+Base.length(::Rect{N,T}) where {T,N} = N
+widths(r::Rect) = r.widths
 
-width(prim::Rect) = prim.widths[1]
-height(prim::Rect) = prim.widths[2]
+width(r::Rect)  = r.widths[1]
+height(r::Rect) = r.widths[2]
 
 """
     split(rectangle, axis, value)
@@ -487,8 +488,9 @@ Check if a point is contained in a Rect. This will return true if
 the point is on a face of the Rect.
 """
 function Base.in(pt::Point, b1::Rect{N,T}) where {T,N}
+    cs = coordinates(pt)
     for i in 1:N
-        pt[i] <= maximum(b1)[i] && pt[i] >= minimum(b1)[i] || return false
+        cs[i] <= maximum(b1)[i] && cs[i] >= minimum(b1)[i] || return false
     end
     return true
 end
@@ -517,29 +519,29 @@ end
 function coordinates(rect::Rect2D, nvertices=(2, 2))
     mini, maxi = extrema(rect)
     xrange, yrange = LinRange.(mini, maxi, nvertices)
-    return ivec(((x, y) for x in xrange, y in yrange))
+    return ivec(Vec(x, y) for x in xrange, y in yrange)
 end
 
 function texturecoordinates(rect::Rect2D, nvertices=(2, 2))
     xrange, yrange = LinRange.((0, 1), (1, 0), nvertices)
-    return ivec(((x, y) for x in xrange, y in yrange))
+    return ivec(Vec(x, y) for x in xrange, y in yrange)
 end
 
-function normals(rect::Rect2D, nvertices=(2, 2))
-    return Iterators.repeated((0, 0, 1), prod(nvertices))
+function normals(::Rect2D, nvertices=(2, 2))
+    return Iterators.repeated(Vec(0, 0, 1), prod(nvertices))
 end
 
 ##
 # Rect3D decomposition
 function coordinates(rect::Rect3D)
     # TODO use n
-    w = widths(rect)
-    o = origin(rect)
-    points = Point{3,Int}[(0, 0, 0), (0, 0, 1), (0, 1, 1), (0, 1, 0), (0, 0, 0), (1, 0, 0),
-                          (1, 0, 1), (0, 0, 1), (0, 0, 0), (0, 1, 0), (1, 1, 0), (1, 0, 0),
-                          (1, 1, 1), (0, 1, 1), (0, 0, 1), (1, 0, 1), (1, 1, 1), (1, 0, 1),
-                          (1, 0, 0), (1, 1, 0), (1, 1, 1), (1, 1, 0), (0, 1, 0), (0, 1, 1)]
-    return ((x .* w .+ o) for x in points)
+    w  = widths(rect)
+    o  = coordinates(origin(rect))
+    xs = Vec{3,Int}[(0, 0, 0), (0, 0, 1), (0, 1, 1), (0, 1, 0), (0, 0, 0), (1, 0, 0),
+                    (1, 0, 1), (0, 0, 1), (0, 0, 0), (0, 1, 0), (1, 1, 0), (1, 0, 0),
+                    (1, 1, 1), (0, 1, 1), (0, 0, 1), (1, 0, 1), (1, 1, 1), (1, 0, 1),
+                    (1, 0, 0), (1, 1, 0), (1, 1, 1), (1, 1, 0), (0, 1, 0), (0, 1, 1)]
+    return ((x .* w .+ o) for x in xs)
 end
 
 function texturecoordinates(rect::Rect3D)
diff --git a/src/primitives/spheres.jl b/src/primitives/spheres.jl
index 7e3df377..fad544d6 100644
--- a/src/primitives/spheres.jl
+++ b/src/primitives/spheres.jl
@@ -2,12 +2,21 @@
     HyperSphere{N, T}
 
 A `HyperSphere` is a generalization of a sphere into N-dimensions.
-A `center` and radius, `r`, must be specified.
+A `center` and `radius` must be specified.
 """
 struct HyperSphere{N,T} <: GeometryPrimitive{N,T}
     center::Point{N,T}
-    r::T
+    radius::T
 end
+
+origin(s::HyperSphere) = s.center
+radius(s::HyperSphere) = s.radius
+
+# TODO: review these
+widths(s::HyperSphere{N,T}) where {N,T} = vfill(Vec{N,T}, 2s.radius)
+Base.minimum(s::HyperSphere) = coordinates(s.center) .- s.radius
+Base.maximum(s::HyperSphere) = coordinates(s.center) .+ s.radius
+
 """
     Circle{T}
 
@@ -22,53 +31,53 @@ An alias for a HyperSphere of dimension 3. (i.e. `HyperSphere{3, T}`)
 """
 const Sphere{T} = HyperSphere{3,T}
 
-HyperSphere{N}(p::Point{N,T}, number) where {N,T} = HyperSphere{N,T}(p, convert(T, number))
-
-widths(c::HyperSphere{N,T}) where {N,T} = Vec{N,T}(radius(c) * 2)
-radius(c::HyperSphere) = c.r
-origin(c::HyperSphere) = c.center
-
-Base.minimum(c::HyperSphere{N,T}) where {N,T} = Vec{N,T}(origin(c)) - Vec{N,T}(radius(c))
-Base.maximum(c::HyperSphere{N,T}) where {N,T} = Vec{N,T}(origin(c)) + Vec{N,T}(radius(c))
+function Base.in(p::AbstractPoint, s::HyperSphere)
+    x = coordinates(p)
+    c = coordinates(s.center)
+    r = s.radius
+    sum(abs2, x - c) ≤ r^2
+end
 
-function Base.in(x::AbstractPoint{2}, c::Circle)
-    @inbounds ox, oy = origin(c)
-    xD = abs(ox - x)
-    yD = abs(oy - y)
-    return xD <= c.r && yD <= c.r
+function centered(S::Type{HyperSphere{N,T}}) where {N,T}
+    center = Point(ntuple(i->zero(T),N))
+    radius = T(0.5)
+    S(center, radius)
 end
 
-centered(S::Type{HyperSphere{N,T}}) where {N,T} = S(Vec{N,T}(0), T(0.5))
 function centered(::Type{T}) where {T<:HyperSphere}
     return centered(HyperSphere{ndims_or(T, 3),eltype_or(T, Float32)})
 end
 
 function coordinates(s::Circle, nvertices=64)
-    rad = radius(s)
-    inner(fi) = Point(rad * sin(fi + pi), rad * cos(fi + pi)) .+ origin(s)
-    return (inner(fi) for fi in LinRange(0, 2pi, nvertices))
+    o = coordinates(s.center)
+    r = s.radius
+    φ = LinRange(0, 2pi, nvertices)
+    inner(φ) = Vec(r*sin(φ+pi), r*cos(φ+pi)) + o
+    return (inner(φ) for φ in φ)
 end
 
-function texturecoordinates(s::Circle, nvertices=64)
-    return coordinates(Circle(Point2f0(0.5), 0.5f0), nvertices)
+function texturecoordinates(::Circle, nvertices=64)
+    return coordinates(HyperSphere(Point2f(0.5,0.5), 0.5f0), nvertices)
 end
 
 function coordinates(s::Sphere, nvertices=24)
+    o = coordinates(s.center)
+    r = s.radius
     θ = LinRange(0, pi, nvertices)
     φ = LinRange(0, 2pi, nvertices)
-    inner(θ, φ) = Point(cos(φ) * sin(θ), sin(φ) * sin(θ), cos(θ)) .* s.r .+ s.center
+    inner(θ, φ) = Vec(r*cos(φ)*sin(θ), r*sin(φ)*sin(θ), r*cos(θ)) + o
     return ivec((inner(θ, φ) for θ in θ, φ in φ))
 end
 
-function texturecoordinates(s::Sphere, nvertices=24)
+function texturecoordinates(::Sphere, nvertices=24)
     ux = LinRange(0, 1, nvertices)
     return ivec(((φ, θ) for θ in reverse(ux), φ in ux))
 end
 
-function faces(sphere::Sphere, nvertices=24)
+function faces(::Sphere, nvertices=24)
     return faces(Rect(0, 0, 1, 1), (nvertices, nvertices))
 end
 
-function normals(s::Sphere{T}, nvertices=24) where {T}
-    return coordinates(Sphere(Point{3,T}(0), 1), nvertices)
+function normals(::Sphere{T}, nvertices=24) where {T}
+    return coordinates(HyperSphere(Point{3,T}(0,0,0), T(1)), nvertices)
 end
diff --git a/src/triangulation.jl b/src/triangulation.jl
index 6aced8f7..3e9d3f43 100644
--- a/src/triangulation.jl
+++ b/src/triangulation.jl
@@ -40,7 +40,9 @@ function area(contour::AbstractVector{<:AbstractPoint{N,T}}) where {N,T}
     p = lastindex(contour)
     q = firstindex(contour)
     while q <= n
-        A += cross(contour[p], contour[q])
+        pc = coordinates(contour[p])
+        qc = coordinates(contour[q])
+        A += cross(pc, qc)
         p = q
         q += 1
     end
@@ -64,28 +66,29 @@ end
 """
 function Base.in(P::T, triangle::Triangle) where {T<:AbstractPoint}
     A, B, C = coordinates(triangle)
-    a = C .- B
-    b = A .- C
-    c = B .- A
+    a = C - B
+    b = A - C
+    c = B - A
 
-    ap = P .- A
-    bp = P .- B
-    cp = P .- C
+    ap = P - A
+    bp = P - B
+    cp = P - C
 
     a_bp = a[1] * bp[2] - a[2] * bp[1]
     c_ap = c[1] * ap[2] - c[2] * ap[1]
     b_cp = b[1] * cp[2] - b[2] * cp[1]
 
-    t0 = zero(eltype(T))
-
-    return ((a_bp >= t0) && (b_cp >= t0) && (c_ap >= t0))
+    return ((a_bp ≥ 0) && (b_cp ≥ 0) && (c_ap ≥ 0))
 end
 
 function snip(contour::AbstractVector{<:AbstractPoint{N,T}}, u, v, w, n, V) where {N,T}
     A = contour[V[u]]
     B = contour[V[v]]
     C = contour[V[w]]
-    x = (((B[1] - A[1]) * (C[2] - A[2])) - ((B[2] - A[2]) * (C[1] - A[1])))
+    a = coordinates(A)
+    b = coordinates(B)
+    c = coordinates(C)
+    x = (((b[1] - a[1]) * (c[2] - a[2])) - ((b[2] - a[2]) * (c[1] - a[1])))
     if 0.0000000001f0 > x
         return false
     end
@@ -112,7 +115,7 @@ function decompose(::Type{FaceType},
     result = FaceType[]
 
     # the algorithm doesn't like closed contours
-    contour = if isapprox(last(points), first(points))
+    contour = if isapprox(coordinates(last(points)), coordinates(first(points)))
         @view points[1:(end - 1)]
     else
         @view points[1:end]
diff --git a/src/viewtypes.jl b/src/viewtypes.jl
index 3232950b..f4878d4f 100644
--- a/src/viewtypes.jl
+++ b/src/viewtypes.jl
@@ -77,7 +77,7 @@ x == [Line(Point(1, 2), Point(3, 4)), Line(Point(5, 6), Point(7, 8))]
     return reinterpret(Polytope(P, Point), TupleView{length(P),skip}(points))
 end
 
-@inline function connect(points::AbstractVector{T}, P::Type{<:Point{N}},
+@inline function connect(points::AbstractVector{T}, ::Type{<:AbstractPoint{N}},
                          skip::Int=N) where {T <: Real,N}
     return reinterpret(Point{N,T}, TupleView{N,skip}(points))
 end
@@ -109,7 +109,7 @@ FaceView enables to link one array of points via a face array, to generate one
 abstract array of elements.
 E.g., this becomes possible:
 ```
-x = FaceView(rand(Point3f0, 10), TriangleFace[(1, 2, 3), (2, 4, 5), ...])
+x = FaceView(rand(Point3f, 10), TriangleFace[(1, 2, 3), (2, 4, 5), ...])
 x[1] isa Triangle == true
 x isa AbstractVector{<: Triangle} == true
 # This means we can use it as a mesh:
diff --git a/test/geometrytypes.jl b/test/geometrytypes.jl
index 56c696cb..dc4513ce 100644
--- a/test/geometrytypes.jl
+++ b/test/geometrytypes.jl
@@ -3,7 +3,7 @@ using Test, GeometryBasics
 @testset "algorithms.jl" begin
     cube = FRect(Vec(-0.5,-0.5,-0.5), Vec(1,1,1))
     cube_faces = decompose(TriangleFace{Int}, faces(cube))
-    cube_vertices = decompose(Point{3,Float32}, cube)
+    cube_vertices = decompose(Point3f, cube)
     @test area(cube_vertices, cube_faces) == 6
     mesh = Mesh(cube_vertices, cube_faces)
     @test GeometryBasics.volume(mesh) ≈ 1
@@ -11,20 +11,18 @@ end
 
 @testset "Cylinder" begin
     @testset "constructors" begin
-        o, extr, r = Point2f0(1, 2), Point2f0(3, 4), 5.0f0
+        o, extr, r = Point2f(1, 2), Point2f(3, 4), 5.0f0
         s = Cylinder(o, extr, r)
         @test typeof(s) == Cylinder{2,Float32}
         @test typeof(s) == Cylinder2{Float32}
         @test origin(s) == o
         @test extremity(s) == extr
         @test radius(s) == r
-        #@test abs(height(s)- norm([1,2]-[3,4]))<1e-5
         h = norm(o - extr)
         @test isapprox(height(s), h)
-        #@test norm(direction(s) - Point{2,Float32}([2,2]./norm([1,2]-[3,4])))<1e-5
-        @test isapprox(direction(s), Point2f0(2, 2) ./ h)
-        v1 = rand(Point{3,Float64})
-        v2 = rand(Point{3,Float64})
+        @test isapprox(direction(s), Vec2f(2, 2) ./ h)
+        v1 = Point(rand(Vec3))
+        v2 = Point(rand(Vec3))
         R = rand()
         s = Cylinder(v1, v2, R)
         @test typeof(s) == Cylinder{3,Float64}
@@ -33,44 +31,43 @@ end
         @test extremity(s) == v2
         @test radius(s) == R
         @test height(s) == norm(v2 - v1)
-        #@test norm(direction(s) - Point{3,Float64}((v2-v1)./norm(v2-v1)))<1e-10
-        @test isapprox(direction(s), (v2 - v1) ./ norm(v2 .- v1))
+        @test isapprox(direction(s), (v2 - v1) ./ norm(v2 - v1))
     end
 
     @testset "decompose" begin
 
-        o, extr, r = Point2f0(1, 2), Point2f0(3, 4), 5.0f0
+        o, extr, r = Point2f(1, 2), Point2f(3, 4), 5.0f0
         s = Cylinder(o, extr, r)
-        positions = Point{3,Float32}[(-0.7677671, 3.767767, 0.0),
-                                     (2.767767, 0.23223293, 0.0),
-                                     (0.23223293, 4.767767, 0.0),
-                                     (3.767767, 1.2322329, 0.0), (1.2322329, 5.767767, 0.0),
-                                     (4.767767, 2.232233, 0.0)]
-        @test decompose(Point3f0, Tesselation(s, (2, 3))) ≈ positions
+        target = Point3f[(-0.7677671, 3.767767, 0.0),
+                         (2.767767, 0.23223293, 0.0),
+                         (0.23223293, 4.767767, 0.0),
+                         (3.767767, 1.2322329, 0.0),
+                         (1.2322329, 5.767767, 0.0),
+                         (4.767767, 2.232233, 0.0)]
+        points = decompose(Point3f, Tesselation(s, (2, 3)))
+        @test coordinates.(points) ≈ coordinates.(target)
 
         FT = TriangleFace{Int}
         faces = FT[(1, 2, 4), (1, 4, 3), (3, 4, 6), (3, 6, 5)]
         @test faces == decompose(FT, Tesselation(s, (2, 3)))
 
-        v1 = Point{3,Float64}(1, 2, 3)
-        v2 = Point{3,Float64}(4, 5, 6)
+        v1 = Point3(1, 2, 3)
+        v2 = Point3(4, 5, 6)
         R = 5.0
         s = Cylinder(v1, v2, R)
-        positions = Point{3,Float64}[(4.535533905932738, -1.5355339059327373, 3.0),
-                                     (7.535533905932738, 1.4644660940672627, 6.0),
-                                     (3.0412414523193148, 4.041241452319315,
-                                      -1.0824829046386295),
-                                     (6.041241452319315, 7.041241452319315,
-                                      1.9175170953613705),
-                                     (-2.535533905932737, 5.535533905932738,
-                                      2.9999999999999996),
-                                     (0.46446609406726314, 8.535533905932738, 6.0),
-                                     (-1.0412414523193152, -0.04124145231931431,
-                                      7.0824829046386295),
-                                     (1.9587585476806848, 2.9587585476806857,
-                                      10.08248290463863), (1, 2, 3), (4, 5, 6)]
-
-        @test decompose(Point3{Float64}, Tesselation(s, 8)) ≈ positions
+        target = Point3[(4.535533905932738, -1.5355339059327373, 3.0),
+                        (7.535533905932738, 1.4644660940672627, 6.0),
+                        (3.0412414523193148, 4.041241452319315, -1.0824829046386295),
+                        (6.041241452319315, 7.041241452319315, 1.9175170953613705),
+                        (-2.535533905932737, 5.535533905932738, 2.9999999999999996),
+                        (0.46446609406726314, 8.535533905932738, 6.0),
+                        (-1.0412414523193152, -0.04124145231931431, 7.0824829046386295),
+                        (1.9587585476806848, 2.9587585476806857, 10.08248290463863),
+                        (1, 2, 3),
+                        (4, 5, 6)]
+
+        points = decompose(Point3, Tesselation(s, 8))
+        @test coordinates.(points) ≈ coordinates.(target)
 
         faces = TriangleFace{Int}[(3, 2, 1), (4, 2, 3), (5, 4, 3), (6, 4, 5), (7, 6, 5),
                                   (8, 6, 7), (1, 8, 7), (2, 8, 1), (3, 1, 9), (2, 4, 10),
@@ -81,8 +78,9 @@ end
         m = triangle_mesh(Tesselation(s, 8))
 
         @test GeometryBasics.faces(m) == faces
-        @test GeometryBasics.coordinates(m) ≈ positions
-        m = normal_mesh(s)# just test that it works without explicit resolution parameter
+        points = metafree(coordinates(m))
+        @test coordinates.(points) ≈ coordinates.(target)
+        m = normal_mesh(s)
         @test m isa GLNormalMesh
 
         muv = uv_mesh(s)
@@ -95,14 +93,14 @@ end
     pt_expa = Point{2,Int}[(0, 0), (1, 0), (0, 1), (1, 1)]
     @test decompose(Point{2,Int}, a) == pt_expa
     mesh = normal_mesh(a)
-    @test decompose(Point2f0, mesh) == pt_expa
+    @test decompose(Point2f, mesh) == convert.(Point2f, pt_expa)
 
     b = Rect(Vec(1, 1, 1), Vec(1, 1, 1))
-    pt_expb = Point{3,Int64}[[1, 1, 1], [1, 1, 2], [1, 2, 2], [1, 2, 1], [1, 1, 1],
-                             [2, 1, 1], [2, 1, 2], [1, 1, 2], [1, 1, 1], [1, 2, 1],
-                             [2, 2, 1], [2, 1, 1], [2, 2, 2], [1, 2, 2], [1, 1, 2],
-                             [2, 1, 2], [2, 2, 2], [2, 1, 2], [2, 1, 1], [2, 2, 1],
-                             [2, 2, 2], [2, 2, 1], [1, 2, 1], [1, 2, 2]]
+    pt_expb = Point{3,Int}[(1, 1, 1), (1, 1, 2), (1, 2, 2), (1, 2, 1), (1, 1, 1),
+                           (2, 1, 1), (2, 1, 2), (1, 1, 2), (1, 1, 1), (1, 2, 1),
+                           (2, 2, 1), (2, 1, 1), (2, 2, 2), (1, 2, 2), (1, 1, 2),
+                           (2, 1, 2), (2, 2, 2), (2, 1, 2), (2, 1, 1), (2, 2, 1),
+                           (2, 2, 2), (2, 2, 1), (1, 2, 1), (1, 2, 2)]
     @test decompose(Point{3,Int}, b) == pt_expb
     mesh = normal_mesh(b)
 end
@@ -136,10 +134,10 @@ end
                (0.0, 0.0, -1.0), (0.0, 0.0, 1.0), (0.0, 0.0, 1.0),
                (0.0, 0.0, 1.0), (0.0, 0.0, 1.0), (-1.0, 0.0, 0.0),
                (-1.0, 0.0, 0.0), (-1.0, 0.0, 0.0), (-1.0, 0.0, 0.0),
-               (1.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.0, 0.0),
-               (0.0, 1.0, 0.0), (0.0, 1.0, 0.0), (0.0, 1.0, 0.0), (0.0, 1.0, 0.0),
-               (0.0, -1.0, 0.0), (0.0, -1.0, 0.0), (0.0, -1.0, 0.0),
-               (0.0, -1.0, 0.0),]
+               (1.0, 0.0, 0.0), (1.0, 0.0, 0.0), (1.0, 0.0, 0.0),
+               (1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, 1.0, 0.0),
+               (0.0, 1.0, 0.0), (0.0, 1.0, 0.0), (0.0, -1.0, 0.0),
+               (0.0, -1.0, 0.0), (0.0, -1.0, 0.0), (0.0, -1.0, 0.0)]
     n32 = map(Vec3f, n64)
     r = triangle_mesh(centered(Rect3D))
     # @test normals(coordinates(r), GeometryBasics.faces(r)) == n32
@@ -147,27 +145,32 @@ end
 end
 
 @testset "HyperSphere" begin
-    sphere = Sphere{Float32}(Point3f0(0), 1.0f0)
-
-    points = decompose(Point, Tesselation(sphere, 3))
-    point_target = Point{3,Float32}[[0.0, 0.0, 1.0], [1.0, 0.0, 6.12323e-17],
-                                    [1.22465e-16, 0.0, -1.0], [-0.0, 0.0, 1.0],
-                                    [-1.0, 1.22465e-16, 6.12323e-17],
-                                    [-1.22465e-16, 1.49976e-32, -1.0], [0.0, -0.0, 1.0],
-                                    [1.0, -2.44929e-16, 6.12323e-17],
-                                    [1.22465e-16, -2.99952e-32, -1.0]]
-    @test points ≈ point_target
+    sphere = Sphere{Float32}(Point3f(0,0,0), 1.0f0)
+
+    points = decompose(Point3f, Tesselation(sphere, 3))
+    target = Point3f[(0.0, 0.0, 1.0),
+                     (1.0, 0.0, 6.12323e-17),
+                     (1.22465e-16, 0.0, -1.0),
+                     (-0.0, 0.0, 1.0),
+                     (-1.0, 1.22465e-16, 6.12323e-17),
+                     (-1.22465e-16, 1.49976e-32, -1.0),
+                     (0.0, -0.0, 1.0),
+                     (1.0, -2.44929e-16, 6.12323e-17),
+                     (1.22465e-16, -2.99952e-32, -1.0)]
+    @test coordinates.(points) ≈ coordinates.(target)
 
     f = decompose(TriangleFace{Int}, Tesselation(sphere, 3))
     face_target = TriangleFace{Int}[[1, 2, 5], [1, 5, 4], [2, 3, 6], [2, 6, 5], [4, 5, 8],
                                     [4, 8, 7], [5, 6, 9], [5, 9, 8]]
     @test f == face_target
-    circle = Circle(Point2f0(0), 1.0f0)
-    points = decompose(Point2f0, Tesselation(circle, 20))
+    circle = HyperSphere(Point2f(0, 0), 1.0f0)
+    points = decompose(Point2f, Tesselation(circle, 20))
     @test length(points) == 20
-    tess_circle = Tesselation(circle, 32)
-    mesh = triangle_mesh(tess_circle)
-    @test decompose(Point2f0, mesh) ≈ decompose(Point2f0, tess_circle)
+    tess = Tesselation(circle, 32)
+    mesh = triangle_mesh(tess)
+    mpoints = decompose(Point2f, mesh)
+    tpoints = decompose(Point2f, tess)
+    @test coordinates.(mpoints) ≈ coordinates.(tpoints)
 end
 
 @testset "Rectangles" begin
@@ -204,8 +207,8 @@ end
 
     split1, split2 = GeometryBasics.split(rect1, 2, 1)
     @test widths(split1) == widths(split2)
-    @test origin(split1) == Vec(0, 0)
-    @test origin(split2) == Vec(0, 1)
+    @test origin(split1) == Point(0.0, 0.0)
+    @test origin(split2) == Point(0.0, 1.0)
     @test in(split1, rect1)
     @test !in(rect1, split1)
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 71ce4dea..9fde2118 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,7 +9,7 @@ using GeometryBasics: attributes
     @testset "Meshes" begin
 
         @testset "per vertex attributes" begin
-            points = rand(Point{3, Float64}, 8)
+            points = rand(Point3, 8)
             tfaces = TetrahedronFace{Int}[(1, 2, 3, 4), (5, 6, 7, 8)]
             normals = rand(Vec3, 8)
             stress = LinRange(0, 1, 8)
@@ -58,7 +58,7 @@ using GeometryBasics: attributes
     end
 
     @testset "polygon with metadata" begin
-        polys = [Polygon(rand(Point{2, Float32}, 20)) for i in 1:10]
+        polys = [Polygon(rand(Point2f, 20)) for i in 1:10]
         pnames = [randstring(4) for i in 1:10]
         numbers = LinRange(0.0, 1.0, 10)
         bin = rand(Bool, 10)
@@ -83,11 +83,7 @@ using GeometryBasics: attributes
     end
     @testset "point with metadata" begin
         p = Point(1.1, 2.2)
-        @test p isa AbstractVector{Float64}
         pm = GeometryBasics.PointMeta(1.1, 2.2; a=1, b=2)
-        p1 = Point(2.2, 3.6)
-        p2 = [p, p1]
-        @test coordinates(p2) == p2
         @test meta(pm) === (a=1, b=2)
         @test metafree(pm) === p
         @test propertynames(pm) == (:position, :a, :b)
@@ -123,7 +119,7 @@ using GeometryBasics: attributes
     end
 
     @testset "Mesh with metadata" begin
-       m = triangle_mesh(Sphere(Point3f0(0), 1))
+       m = triangle_mesh(HyperSphere(Point3f(0,0,0), 1.0f0))
        m_meta = MeshMeta(m; boundingbox=Rect(1.0, 1.0, 2.0, 2.0))
        @test meta(m_meta) === (boundingbox = Rect(1.0, 1.0, 2.0, 2.0),)
        @test metafree(m_meta) === m
@@ -133,7 +129,7 @@ end
 
 @testset "embedding MetaT" begin
     @testset "MetaT{Polygon}" begin
-        polys = [Polygon(rand(Point{2, Float32}, 20)) for i in 1:10]
+        polys = [Polygon(rand(Point2f, 20)) for i in 1:10]
         multipol = MultiPolygon(polys)
         pnames = [randstring(4) for i in 1:10]
         numbers = LinRange(0.0, 1.0, 10)
@@ -160,11 +156,7 @@ end
 
     @testset "MetaT{Point}" begin
         p = Point(1.1, 2.2)
-        @test p isa AbstractVector{Float64}
         pm = MetaT(Point(1.1, 2.2); a=1, b=2)
-        p1 = Point(2.2, 3.6)
-        p2 = [p, p1]
-        @test coordinates(p2) == p2
         @test pm.meta === (a=1, b=2)
         @test pm.main === p
         @test propertynames(pm) == (:main, :a, :b)
@@ -173,14 +165,13 @@ end
     end
 
     @testset "MetaT{MultiPoint}" begin
-        p = collect(Point{2, Float64}(x, x+1) for x in 1:5)
-        @test p isa AbstractVector
-        mpm = MetaT(MultiPoint(p); a=1, b=2)
-        @test coordinates(mpm.main) == Point{2, Float64}[(x, x+1) for x in 1:5]
+        ps = [Point2(x, x+1) for x in 1:5]
+        @test ps isa AbstractVector
+        mpm = MetaT(MultiPoint(ps); a=1, b=2)
+        @test mpm.main == ps
         @test mpm.meta === (a=1, b=2)
-        @test mpm.main == p
         @test propertynames(mpm) == (:main, :a, :b)
-        @test GeometryBasics.metafree(mpm) == p
+        @test GeometryBasics.metafree(mpm) == ps
         @test GeometryBasics.meta(mpm) == (a = 1, b = 2)
     end
 
@@ -213,7 +204,7 @@ end
 
    @testset "MetaT{Mesh}" begin
         @testset "per vertex attributes" begin
-            points = rand(Point{3, Float64}, 8)
+            points = rand(Point3, 8)
             tfaces = TetrahedronFace{Int}[(1, 2, 3, 4), (5, 6, 7, 8)]
             normals = rand(Vec3, 8)
             stress = LinRange(0, 1, 8)
@@ -252,16 +243,16 @@ end
 
     @testset "connected views" begin
         numbers = [1, 2, 3, 4, 5, 6]
-        x = connect(numbers, Point{2})
+        x = connect(numbers, Point2)
 
-        @test x == Point[(1, 2), (3, 4), (5, 6)]
+        @test x == Point{2,Int}[(1, 2), (3, 4), (5, 6)]
 
         line = connect(x, Line, 1)
         @test line == [Line(Point(1, 2), Point(3, 4)), Line(Point(3, 4), Point(5, 6))]
 
         triangles = connect(x, Triangle)
         @test triangles == [Triangle(Point(1, 2), Point(3, 4), Point(5, 6))]
-        x = connect([1, 2, 3, 4, 5, 6, 7, 8], Point{2})
+        x = connect([1, 2, 3, 4, 5, 6, 7, 8], Point2)
         tetrahedra = connect(x, NSimplex{4})
         @test tetrahedra == [Tetrahedron(x[1], x[2], x[3], x[4])]
 
@@ -283,10 +274,10 @@ end
         faces = connect([1, 2, 3], TriangleFace)
         triangles = connect(points, faces)
         @test triangles == [Triangle(Point(1, 2), Point(3, 4), Point(5, 6))]
-        x = Point{3}(1.0)
+        x = Point3(1,1,1)
         triangles = connect([x], [TriangleFace(1, 1, 1)])
         @test triangles == [Triangle(x, x, x)]
-        points = connect([1, 2, 3, 4, 5, 6, 7, 8], Point{2})
+        points = connect([1, 2, 3, 4, 5, 6, 7, 8], Point2)
         faces = connect([1, 2, 3, 4], SimplexFace{4})
         triangles = connect(points, faces)
         @test triangles == [Tetrahedron(points...)]
@@ -309,7 +300,7 @@ end
         linestring = LineString(points)
         @test linestring == [Line(points[1], points[2]), Line(points[2], points[3])]
 
-        points = rand(Point{2, Float64}, 4)
+        points = rand(Point2, 4)
         linestring = LineString(points, 2)
         @test linestring == [Line(points[1], points[2]), Line(points[3], points[4])]
 
@@ -332,11 +323,11 @@ end
 
     @testset "Polygon" begin
 
-        points = connect([1, 2, 3, 4, 5, 6], Point{2})
+        points = connect([1, 2, 3, 4, 5, 6], Point2)
         polygon = Polygon(points)
         @test polygon == Polygon(LineString(points))
 
-        points = rand(Point{2, Float64}, 4)
+        points = rand(Point2, 4)
         linestring = LineString(points, 2)
         @test Polygon(points, 2) == Polygon(linestring)
 
@@ -358,37 +349,36 @@ end
     @testset "Mesh" begin
 
         numbers = [1, 2, 3, 4, 5, 6]
-        points = connect(numbers, Point{2})
-
+        points = connect(numbers, Point2)
         mesh = Mesh(points, [1,2,3])
         @test mesh == [Triangle(points...)]
 
-        x = Point{3}(1.0)
+        x = Point3(1,1,1)
         mesh = Mesh([x], [TriangleFace(1, 1, 1)])
         @test mesh == [Triangle(x, x, x)]
 
-        points = connect([1, 2, 3, 4, 5, 6, 7, 8], Point{2})
-        faces = connect([1, 2, 3, 4], SimplexFace{4})
-        mesh = Mesh(points, faces)
+        points = connect([1, 2, 3, 4, 5, 6, 7, 8], Point2)
+        sfaces = connect([1, 2, 3, 4], SimplexFace{4})
+        mesh = Mesh(points, sfaces)
         @test mesh == [Tetrahedron(points...)]
 
-        points = rand(Point3f0, 8)
+        points = rand(Point3f, 8)
         tfaces = [GLTriangleFace(1, 2, 3), GLTriangleFace(5, 6, 7)]
-        normals = rand(Vec3f, 8)
+        normal = rand(Vec3f, 8)
         uv = rand(Vec2f, 8)
         mesh = Mesh(points, tfaces)
         meshuv = Mesh(meta(points; uv=uv), tfaces)
-        meshuvnormal = Mesh(meta(points; normals=normals, uv=uv), tfaces)
+        meshuvnormal = Mesh(meta(points; normals=normal, uv=uv), tfaces)
 
         @test mesh isa GLPlainMesh
         @test meshuv isa GLUVMesh3D
         @test meshuvnormal isa GLNormalUVMesh3D
 
         t = Tesselation(FRect2D(0, 0, 2, 2), (30, 30))
-        m = GeometryBasics.mesh(t, pointtype=Point3f0, facetype=QuadFace{Int})
-        m2 = GeometryBasics.mesh(m, facetype=QuadFace{GLIndex})
+        m = GeometryBasics.mesh(t, pointtype=Point2f, facetype=QuadFace{Int})
+        m2 = GeometryBasics.mesh(m, pointtype=Point2f, facetype=QuadFace{GLIndex})
         @test GeometryBasics.faces(m2) isa Vector{QuadFace{GLIndex}}
-        @test GeometryBasics.coordinates(m2) isa Vector{Point3f0}
+        @test GeometryBasics.coordinates(m2) isa Vector{Point2f}
 
     end
 
@@ -425,22 +415,22 @@ end
 end
 
 @testset "decompose/triangulation" begin
-    primitive = Sphere(Point3f0(0), 1)
+    primitive = HyperSphere(Point3f(0,0,0), 1.0f0)
     @test ndims(primitive) === 3
     mesh = triangle_mesh(primitive)
-    @test decompose(Point, mesh) isa Vector{Point3f0}
-    @test decompose(Point, primitive) isa Vector{Point3f0}
+    @test decompose(Point3f, mesh) isa Vector{Point3f}
+    @test decompose(Point3f, primitive) isa Vector{Point3f}
 
     primitive = Rect2D(0, 0, 1, 1)
     mesh = triangle_mesh(primitive)
 
-    @test decompose(Point, mesh) isa Vector{Point2f0}
-    @test decompose(Point, primitive) isa Vector{Point2{Int}}
+    @test decompose(Point2f, mesh) isa Vector{Point2f}
+    @test decompose(Point{2,Int}, primitive) isa Vector{Point{2,Int}}
 
     primitive = Rect3D(0, 0, 0, 1, 1, 1)
     triangle_mesh(primitive)
 
-    primitive = Sphere(Point3f0(0), 1)
+    primitive = HyperSphere(Point3f(0,0,0), 1.0f0)
     m_normal = normal_mesh(primitive)
     @test normals(m_normal) isa Vector{Vec3f}
     primitive = Rect2D(0, 0, 1, 1)
@@ -450,32 +440,32 @@ end
     m_normal = normal_mesh(primitive)
     @test normals(m_normal) isa Vector{Vec3f}
 
-    points = decompose(Point2f0, Circle(Point2f0(0), 1))
+    points = decompose(Point2f, HyperSphere(Point2f(0, 0), 1.0f0))
     tmesh = triangle_mesh(points)
     @test normals(tmesh) == nothing
 
-    m = GeometryBasics.mesh(Sphere(Point3f0(0), 1))
-    @test normals(m) == nothing
+    m = GeometryBasics.mesh(HyperSphere(Point3f(0,0,0), 1.0f0))
+    @test normals(m) === nothing
     m_normals = pointmeta(m, Normal())
     @test normals(m_normals) isa Vector{Vec3f}
 
     @test texturecoordinates(m) == nothing
     r2 = Rect2D(0.0, 0.0, 1.0, 1.0)
-    @test iterate(texturecoordinates(r2)) == ((0.0, 1.0), ((0.0, 2), (1.0, 2)))
+    @test iterate(texturecoordinates(r2)) == ([0.0, 1.0], ((0.0, 2), (1.0, 2)))
     r3 = Rect3D(0.0, 0.0, 1.0, 1.0, 2.0, 2.0)
     @test iterate(texturecoordinates(r3)) == ([0, 0, 0], 2)
     uv = decompose_uv(m)
     @test boundingbox(Point.(uv)) == Rect(0, 0, 0, 1, 1, 1)
 
-    points = decompose(Point2f0, Circle(Point2f0(0), 1))
+    points = decompose(Point2f, HyperSphere(Point2f(0, 0), 1.0f0))
     m = GeometryBasics.mesh(points)
-    @test coordinates(m) === points
+    @test coordinates(m) == points
 
     linestring = LineString(Point{2, Int}[(10, 10), (20, 20), (10, 40)])
     pts = Point{2, Int}[(10, 10), (20, 20), (10, 40)]
     linestring = LineString(pts)
     pts_decomp = decompose(Point{2, Int}, linestring)
-    @test pts === pts_decomp
+    @test pts == pts_decomp
 
     pts_ext = Point{2, Int}[(5, 1), (3, 3), (4, 8), (1, 2), (5, 1)]
     ls_ext = LineString(pts_ext)
@@ -490,7 +480,7 @@ end
 end
 
 @testset "convert mesh + meta" begin
-    m = uv_normal_mesh(Circle(Point2f0(0), 1f0))
+    m = uv_normal_mesh(HyperSphere(Point2f(0, 0), 1.0f0))
     # for 2D primitives we dont actually calculate normals
     @test !hasproperty(m, :normals)
 end
@@ -498,16 +488,14 @@ end
 @testset "convert mesh + meta" begin
     m = uv_normal_mesh(FRect3D(Vec(-1,-1,-1), Vec(1, 2, 3)))
     m_normal = normal_mesh(m)
-    # make sure we don't loose the uv
     @test hasproperty(m_normal, :uv)
     @test m == m_normal
-    # Make sure we don't create any copies
-    @test m.position === m_normal.position
-    @test m.normals === m_normal.normals
-    @test m.uv === m_normal.uv
+    @test m.position == m_normal.position
+    @test m.normals == m_normal.normals
+    @test m.uv == m_normal.uv
 
     m = GeometryBasics.mesh(FRect3D(Vec(-1,-1,-1), Vec(1, 2, 3));
-                            uv=Vec2, normaltype=Vec3, pointtype=Point3{Float64})
+                            uv=Vec2, normaltype=Vec3, pointtype=Point3)
     m_normal = normal_mesh(m)
     @test hasproperty(m_normal, :uv)
     @test m.position !== m_normal.position
@@ -518,7 +506,7 @@ end
 
 @testset "modifying meta" begin
     xx = rand(10)
-    points = rand(Point3f0, 10)
+    points = rand(Point3f, 10)
     m = GeometryBasics.Mesh(meta(points, xx=xx), GLTriangleFace[(1,2,3), (3,4,5)])
     color = rand(10)
     m = pointmeta(m; color=color)
@@ -540,37 +528,44 @@ end
     @test xxpopt === xx
 
     @testset "creating meta" begin
-        x = Point3f0[(1,3,4)]
+        x = Point3f[(1,3,4)]
         # no meta gets added, so should stay the same
         @test meta(x) === x
         @test meta(x, value=[1]).position === x
     end
-    pos = Point2f0[(10, 2)]
+    pos = Point2f[(10, 2)]
     m = Mesh(meta(pos, uv=[Vec2f(1, 1)]), [GLTriangleFace(1, 1, 1)])
     @test m.position === pos
 end
 
 @testset "mesh conversion" begin
-    s = Sphere(Point3(0.0), 1.0)
-    m = GeometryBasics.mesh(s)
-    @test m isa Mesh{3, Float64}
-    @test coordinates(m) isa Vector{Point{3, Float64}}
+    s = HyperSphere(Point3(0,0,0), 1.0)
+    m = GeometryBasics.mesh(s, pointtype=Point3)
+    @test m isa Mesh{3,Float64}
+    @test coordinates(m) isa Vector{Point3}
     @test GeometryBasics.faces(m) isa Vector{GLTriangleFace}
-    # Check, that decompose isn't making a copy for matching eltype
-    @test coordinates(m) === decompose(Point{3, Float64}, m)
+    points1 = coordinates(m)
+    points2 = decompose(Point3, m)
+    @test coordinates.(points1) ≈ coordinates.(points2)
 
+    m = GeometryBasics.mesh(s, pointtype=Point3f)
     tmesh = triangle_mesh(m)
     @test tmesh isa GLPlainMesh
-    @test coordinates(tmesh) === decompose(Point3f0, tmesh)
+    points1 = coordinates(tmesh)
+    points2 = decompose(Point3f, tmesh)
+    @test coordinates.(points1) ≈ coordinates.(points2)
 
     nmesh = normal_mesh(m)
     @test nmesh isa GLNormalMesh
-    @test metafree(coordinates(nmesh)) === decompose(Point3f0, nmesh)
-    @test normals(nmesh) === decompose_normals(nmesh)
-
-    m = GeometryBasics.mesh(s, pointtype=Point3f0)
-    @test m isa Mesh{3, Float32}
-    @test coordinates(m) isa Vector{Point3f0}
+    points1 = metafree(coordinates(nmesh))
+    points2 = decompose(Point3f, nmesh)
+    @test coordinates.(points1) ≈ coordinates.(points2)
+    normals1 = normals(nmesh)
+    normals2 = decompose_normals(nmesh)
+    @test isequal(normals1, normals2)
+
+    @test m isa Mesh{3,Float32}
+    @test coordinates(m) isa Vector{Point3f}
     @test GeometryBasics.faces(m) isa Vector{GLTriangleFace}
 end
 
@@ -586,7 +581,7 @@ end
     @test intersects(a, c) === (false, Point(0.0, 0.0))
     @test intersects(d, d) === (false, Point(0.0, 0.0))
     found, point = intersects(d, e)
-    @test found && point ≈ Point(0.0, 4.0)
+    @test found && coordinates(point) ≈ [0.0, 4.0]
     @test intersects(a, f) === (false, Point(0.0, 0.0))
 end
 
@@ -622,8 +617,8 @@ end
 end
 
 @testset "MetaT and heterogeneous data" begin
-    ls = [LineString([Point(i, (i+1)^2/6), Point(i*0.86,i+5), Point(i/3, i/7)]) for i in 1:10]
-    mls = MultiLineString([LineString([Point(i+1, (i)^2/6), Point(i*0.75,i+8), Point(i/2.5, i/6.79)]) for i in 5:10])
+    ls = [LineString([Point(i*1.0, (i+1)^2/6), Point(i*0.86,i+5.0), Point(i/3, i/7)]) for i in 1:10]
+    mls = MultiLineString([LineString([Point(i+1.0, (i)^2/6), Point(i*0.75,i+8.0), Point(i/2.5, i/6.79)]) for i in 5:10])
     poly = Polygon(Point{2, Int}[(40, 40), (20, 45), (45, 30), (40, 40)])
     geom = [ls..., mls, poly]
     prop = Any[(country_states = "India$(i)", rainfall = (i*9)/2) for i in 1:11]
@@ -642,8 +637,7 @@ end
     (LineString{2,Float64,Point{2,Float64},Base.ReinterpretArray{GeometryBasics.Ngon{2,Float64,2,Point{2,Float64}},1,Tuple{Point{2,Float64},Point{2,Float64}},TupleView{Tuple{Point{2,Float64},Point{2,Float64}},2,1,Array{Point{2,Float64},1}}}},
     (:country_states, :rainfall), Tuple{String,Float64})
 
-
-    @test StructArrays.createinstance(typeof(feat[1]), LineString([Point(1, (2)^2/6), Point(1*0.86,6), Point(1/3, 1/7)]), "Mumbai", 100) isa typeof(feat[1])
+    @test StructArrays.createinstance(typeof(feat[1]), LineString([Point(1.0, (2)^2/6), Point(1*0.86,6.0), Point(1/3, 1/7)]), "Mumbai", 100) isa typeof(feat[1])
 
     @test Base.getindex(feat[1], 1) isa Line
     @test Base.size(feat[1]) == (2,)