diff --git a/src/BiodiversityObservationNetworks.jl b/src/BiodiversityObservationNetworks.jl
index ef45d16..cfb0cfe 100644
--- a/src/BiodiversityObservationNetworks.jl
+++ b/src/BiodiversityObservationNetworks.jl
@@ -37,6 +37,9 @@ export AdaptiveSpatial
 include("cubesampling.jl")
 export CubeSampling
 
+include("fractaltriad.jl")
+export FractalTriad
+
 include("uniqueness.jl")
 export Uniqueness
 
diff --git a/src/adaptivespatial.jl b/src/adaptivespatial.jl
index aee82ac..3482712 100644
--- a/src/adaptivespatial.jl
+++ b/src/adaptivespatial.jl
@@ -5,9 +5,9 @@
 
 **numsites**, an Integer (def. 50), specifying the number of points to use.
 """
-Base.@kwdef mutable struct AdaptiveSpatial{T <: Integer, F<: AbstractFloat} <: BONRefiner
+Base.@kwdef mutable struct AdaptiveSpatial{T <: Integer, F <: AbstractFloat} <: BONRefiner
     numsites::T = 30
-    uncertainty::Array{F,2} = rand(50,50)
+    uncertainty::Array{F, 2} = rand(50, 50)
     function AdaptiveSpatial(numsites, uncertainty)
         as = new{typeof(numsites), typeof(uncertainty[begin])}(numsites, uncertainty)
         check_arguments(as)
@@ -15,18 +15,18 @@ Base.@kwdef mutable struct AdaptiveSpatial{T <: Integer, F<: AbstractFloat} <: B
     end
 end
 
+maxsites(as::AdaptiveSpatial) = prod(size(as.uncertainty))
 function check_arguments(as::AdaptiveSpatial)
     check(TooFewSites, as)
-    
-    max_num_sites = prod(size(as.uncertainty))
-    check(TooManySites, as, max_num_sites)
+    check(TooManySites, as)
+    return nothing
 end
 
 function _generate!(
     coords::Vector{CartesianIndex},
     pool::Vector{CartesianIndex},
     sampler::AdaptiveSpatial,
-) 
+)
     # Distance matrix (inlined)
     d = zeros(Float64, Int((sampler.numsites * (sampler.numsites - 1)) / 2))
 
@@ -79,7 +79,6 @@ function _D(a1::T, a2::T) where {T <: CartesianIndex{2}}
     return sqrt((x1 - x2)^2.0 + (y1 - y2)^2.0)
 end
 
-
 # ====================================================
 #
 #   Tests
@@ -99,4 +98,4 @@ end
     @test_throws TooFewSites AdaptiveSpatial(numsites = 1)
     @test_throws TooFewSites AdaptiveSpatial(numsites = 0)
     @test_throws TooFewSites AdaptiveSpatial(numsites = -1)
-end
\ No newline at end of file
+end
diff --git a/src/balancedacceptance.jl b/src/balancedacceptance.jl
index 6f13abb..3845159 100644
--- a/src/balancedacceptance.jl
+++ b/src/balancedacceptance.jl
@@ -18,7 +18,8 @@ maxsites(bas::BalancedAcceptance) = prod(bas.dims)
 
 function check_arguments(bas::BalancedAcceptance)
     check(TooFewSites, bas)
-    check(TooManySites, bas, maxsites(bas))
+    check(TooManySites, bas)
+    return nothing
 end
 
 function _generate!(
diff --git a/src/cubesampling.jl b/src/cubesampling.jl
index a2afa90..001b0cb 100644
--- a/src/cubesampling.jl
+++ b/src/cubesampling.jl
@@ -22,7 +22,6 @@ Base.@kwdef struct CubeSampling{I <: Integer, M <: Matrix, V <: Vector} <: BONRe
     function CubeSampling(numsites, fast, x, πₖ)
         cs = new{typeof(numsites), typeof(x), typeof(πₖ)}(numsites, fast, x, πₖ)
         _check_arguments(cs)
-
         return cs
     end
 end
@@ -32,22 +31,23 @@ maxsites(cubesampling::CubeSampling) = size(cubesampling.x, 2)
 
 function check_arguments(cubesampling::CubeSampling)
     check(TooFewSites, cubesampling)
-    check(TooManySites, maxsites(cubesampling))
+    check(TooManySites, cubesampling)
 
-    if numsites > length(πₖ)
+    if numsites > length(cubesampling.πₖ)
         throw(
             ArgumentError(
                 "You cannot select more points than the number of candidate points.",
             ),
         )
     end
-    if length(πₖ) != size(x, 2)
+    if length(cubesampling.πₖ) != size(cubesampling.x, 2)
         throw(
             DimensionMismatch(
                 "The number of inclusion probabilites does not match the dimensions of the auxillary variable matrix.",
             ),
         )
     end
+    return
 end
 
 function check_conditions(coords, pool, sampler)
diff --git a/src/exceptions.jl b/src/exceptions.jl
index b2ea8d9..c012d22 100644
--- a/src/exceptions.jl
+++ b/src/exceptions.jl
@@ -8,19 +8,22 @@ Base.showerror(io::IO, e::E) where {E <: BONException} =
     )
 
 function _check_arguments(sampler::S) where {S <: Union{BONSeeder, BONRefiner}}
-    return sampler.numsites > 1 || throw(TooFewSites(sampler.numsites))
+    sampler.numsites > 1 || throw(TooFewSites())
+    return nothing
 end
 
 @kwdef struct TooFewSites <: BONException
     message = "Number of sites to select must be at least two."
 end
-function check(TooFewSites, sampler)
-    return sampler.numsites > 1 || throw(TooFewSites())
+function check(::Type{TooFewSites}, sampler)
+    sampler.numsites > 1 || throw(TooFewSites())
+    return nothing
 end
 
 @kwdef struct TooManySites <: BONException
     message = "Cannot select more sites than there are candidates."
 end
-function check(TooManySites, sampler, max_sites)
-    return sampler.numsites <= max_sites || throw(TooManySites())
+function check(::Type{TooManySites}, sampler)
+    sampler.numsites <= maxsites(sampler) || throw(TooManySites())
+    return nothing
 end
diff --git a/src/fractaltriad.jl b/src/fractaltriad.jl
index 79ad3c6..c0b3b23 100644
--- a/src/fractaltriad.jl
+++ b/src/fractaltriad.jl
@@ -1,12 +1,184 @@
-Base.@kwdef struct FractalTriad{I <: Integer, F <: AbstractFloat} <: SpatialSampler
-    numsites::I = 50
+"""
+    FractalTriad
+
+A `BONSeeder` that generates `FractalTriad` designs
+"""
+Base.@kwdef struct FractalTriad{I <: Integer, F <: AbstractFloat} <: BONSeeder
+    numsites::I = 81
     horizontal_padding::F = 0.1
-    vetical_padding::F = 0.1
-    dims::Tuple{I, I}
+    vertical_padding::F = 0.1
+    dims::Tuple{I, I} = (100, 100)
+    function FractalTriad(numsites, horizontal_padding, vertical_padding, dims)
+        ft = new{typeof(numsites), typeof(horizontal_padding)}(
+            numsites,
+            horizontal_padding,
+            vertical_padding,
+            dims,
+        )
+        check_arguments(ft)
+        return ft
+    end
+end
+
+maxsites(ft::FractalTriad) = maximum(ft.dims) * 10  # gets numerically unstable for large values because float coords belong to the the same cell in the raster, and arctan breaks  
+function check_arguments(ft::FractalTriad)
+    check(TooManySites, ft)
+    ft.numsites > 9 || throw(TooFewSites("FractalTriad requires more than 9 or more sites"))
+    ft.numsites % 9 == 0 ||
+        throw(ArgumentError("FractalTriad requires number of sites to be a multiple of 9"))
+    return
+end
+
+"""
+    _outer_triangle
+
+Takes a FractalTriad generator `ft` and returns the outermost triangle based on the size of the raster and the padding on the horizontal and vertical.
+"""
+function _outer_triangle(ft)
+    x, y = ft.dims
+    Δx, Δy =
+        Int.([floor(0.5 * ft.horizontal_padding * x), floor(0.5 * ft.vertical_padding * y)])
+
+    return CartesianIndex.([(Δx, Δy), (x ÷ 2, y - Δy), (x - Δx, Δy)])
 end
 
-function _generate!(ft::FractalTriad)
-    response = zeros(ft.numsites, 2)
+"""
+    _hexagon(A, B, C)
+
+Takes a set of `CartesianIndex`'s that form the vertices of a triangle, and returns the `CartisianIndex`'s for each of the points that form the internal hexagon points
+
+For input vertices A, B, C, `_hexagon` returns the points on the edges of the triangle below in the order `[d,e,f,g,h,i]`
+
+                        B 
+
+                  e           f
+
+
+             d                      g
+
+        A          i          h            C
+
+--
+
+After running `vcat(triangle, hex)`, the resulting indices form the 2-level triad with indices corresponding to points in the below manner:
+
+                         2 
+
+
+                  5           6
+
+
+             4                      7
+
+        1         9            8         3
+
+  - 
+γ = |AB|
+χ = |BC|
+
+θ = ⦤ BAC
+α = ⦤ BCA
+
+TODO: this always assumes |AC| is horizontal. This could be changed later.
+"""
+function _hexagon(A, B, C)
+    γ = sqrt((B[1] - A[1])^2 + (B[2] - A[2])^2) # left side length
+    χ = sqrt((B[1] - C[1])^2 + (B[2] - C[2])^2) # right side length
+
+    θ = atan((B[2] - A[2]) / (B[1] - A[1]))
+    α = atan((B[2] - C[2]) / (C[1] - B[1]))
+
+    d = (A[1] + (γ * cos(θ) / 3), A[2] + γ * sin(θ) / 3)
+    e = (A[1] + 2γ * cos(θ) / 3, A[2] + 2γ * sin(θ) / 3)
+    f = (A[1] + (C[1] - A[1]) - (2χ * cos(α) / 3), A[2] + 2χ * sin(α) / 3)
+    g = (A[1] + (C[1] - A[1] - (χ * cos(α) / 3)), A[2] + χ * sin(α) / 3)
+    h = (A[1] + 2(C[1] - A[1]) / 3, A[2])
+    i = (A[1] + (C[1] - A[1]) / 3, A[2])
+    return [CartesianIndex(Int.([floor(x[1]), floor(x[2])])...) for x in [d, e, f, g, h, i]]
+end
+
+"""
+    _fill_triangle!(coords, traingle, count)
+
+Takes a set of vertices of a triangle `triangle`, and fills the internal hexagon for those points.
+"""
+function _fill_triangle!(coords, triangle, count)
+    start = count
+    hex = _hexagon(triangle...)
+    for i in eachindex(hex)
+        coords[count] = CartesianIndex(hex[i])
+        count += 1
+        if count > length(coords)
+            return coords[start:(count - 1)], count
+        end
+    end
+
+    return coords[start:(count - 1)], count
+end
+
+"""
+    _generate!(coords::Vector{CartesianIndex}, ft::FractalTriad)
+
+Fills `coords` with a set of points generated using the `FractalTriad` generator `ft`.
+"""
+function _generate!(
+    coords::Vector{CartesianIndex},
+    ft::FractalTriad,
+)
+    base_triangle = _outer_triangle(ft)
+    coords[1:3] .= base_triangle
+    count = 4
+
+    triangle = coords[1:3]
+    hex, count = _fill_triangle!(coords, triangle, count)
+    pack = vcat(triangle, hex)
+    vert_idxs = [[5, 2, 6], [1, 4, 9], [8, 7, 3]]
+
+    pack_stack = [pack]
+
+    while length(pack_stack) > 0
+        pack = popat!(pack_stack, 1)
+        for idx in vert_idxs
+            triangle = pack[idx]
+            hex, count = _fill_triangle!(coords, triangle, count)
+            if count > ft.numsites
+                return coords
+            end
+            push!(pack_stack, vcat(triangle, hex))
+        end
+    end
+    return coords
+end
+
+# ====================================================
+#
+#   Tests
+#
+# =====================================================
+
+@testitem "FractalTriad is correct subtype" begin
+    @test FractalTriad <: BONSeeder
+    @test FractalTriad <: BONSampler
+end
+
+@testitem "FractalTriad default constructor works" begin
+    @test typeof(FractalTriad()) <: FractalTriad
+end
+
+@testitem "FractalTriad can change number of sites with keyword argument" begin
+    ft = FractalTriad(; numsites = 18)
+    @test typeof(ft) <: FractalTriad
+    @test ft.numsites == 18
+
+    ft = FractalTriad(; numsites = 27)
+    @test typeof(ft) <: FractalTriad
+    @test ft.numsites == 27
+
+    @test_throws ArgumentError FractalTriad(numsites = 20)
+end
 
-    return response
+@testitem "FractalTriad throws error when too few points as passed as an argument" begin
+    @test_throws TooFewSites FractalTriad(; numsites = 9)
+    @test_throws TooFewSites FractalTriad(; numsites = -1)
+    @test_throws TooFewSites FractalTriad(; numsites = 0)
 end
diff --git a/src/simplerandom.jl b/src/simplerandom.jl
index 82f34f4..444c38a 100644
--- a/src/simplerandom.jl
+++ b/src/simplerandom.jl
@@ -12,11 +12,12 @@ Base.@kwdef struct SimpleRandom{I <: Integer} <: BONSeeder
         return srs
     end
 end
+maxsites(srs::SimpleRandom) = prod(srs.dims)
 
-function check_arguments(srs::SimpleRandom)
+function check_arguments(srs::SRS) where {SRS <: SimpleRandom}
     check(TooFewSites, srs)
-    max_num_sites = prod(srs.dims)
-    check(TooManySites, srs, max_num_sites)
+    check(TooManySites, srs)
+    return nothing
 end
 
 function _generate!(
diff --git a/src/spatialstratified.jl b/src/spatialstratified.jl
index 2e23a4a..e03a34b 100644
--- a/src/spatialstratified.jl
+++ b/src/spatialstratified.jl
@@ -20,8 +20,7 @@ maxsites(ss::SpatiallyStratified) = prod(size(ss.strata))
 
 function check_arguments(ss::SpatiallyStratified)
     check(TooFewSites, ss)
-    check(TooManySites, ss, maxsites(ss))
-
+    check(TooManySites, ss)
 
     length(unique(ss.strata)) == length(ss.inclusion_probability_by_stratum) || throw(
         ArgumentError(
@@ -29,8 +28,9 @@ function check_arguments(ss::SpatiallyStratified)
         ),
     )
 
-    return sum(ss.inclusion_probability_by_stratum) ≈ 1.0 ||
-           throw(ArgumentError("Inclusion probabilities for each strata do not sum to 1."))
+    sum(ss.inclusion_probability_by_stratum) ≈ 1.0 ||
+        throw(ArgumentError("Inclusion probabilities for each strata do not sum to 1."))
+    return nothing
 end
 
 function _default_strata(sz)
diff --git a/src/uniqueness.jl b/src/uniqueness.jl
index f486a53..60b7425 100644
--- a/src/uniqueness.jl
+++ b/src/uniqueness.jl
@@ -17,14 +17,16 @@ maxsites(uniq::Uniqueness) = prod(size(uniq.layers)[1:2])
 
 function check_arguments(uniq::Uniqueness)
     check(TooFewSites, uniq)
-    check(TooManySites, uniq, maxsites(uniq))
+    check(TooManySites, uniq)
 
-    return length(size(uniq.layers)) == 3 ||
-           throw(
-        ArgumentError(
-            "You cannot have a Uniqueness sampler without layers passed as a 3-dimenisional array, where the first two axes are latitude and longitude.",
-        ),
-    )
+    length(size(uniq.layers)) == 3 ||
+        throw(
+            ArgumentError(
+                "You cannot have a Uniqueness sampler without layers passed as a 3-dimenisional array, where the first two axes are latitude and longitude.",
+            ),
+        )
+
+    return nothing
 end
 
 function _generate!(
diff --git a/src/weightedbas.jl b/src/weightedbas.jl
index 12ef085..83b8e89 100644
--- a/src/weightedbas.jl
+++ b/src/weightedbas.jl
@@ -14,19 +14,16 @@ Base.@kwdef struct WeightedBalancedAcceptance{I <: Integer, F <: Real} <: BONSee
     end
 end
 
+maxsites(wbas::WeightedBalancedAcceptance) = prod(size(wbas.uncertainty))
+
 function check_arguments(wbas::WeightedBalancedAcceptance)
     check(TooFewSites, wbas)
-
-    max_num_sites = prod(size(wbas.uncertainty))
-    max_num_sites >= wbas.numsites || throw(
-        TooManySites(
-            "Number of sites to select $(wbas.numsites) is greater than number of possible sites $(max_num_sites)",
-        ),
-    )
-    return wbas.α > 0 ||
-           throw(
-        ArgumentError("WeightedBalancedAcceptance requires α to be greater than 0 "),
-    )
+    check(TooManySites, wbas)
+    wbas.α > 0 ||
+        throw(
+            ArgumentError("WeightedBalancedAcceptance requires α to be greater than 0 "),
+        )
+    return nothing
 end
 
 function _generate!(
@@ -124,13 +121,13 @@ end
     @test numpts == length(coords)
 end
 
-@testitem "BalancedAcceptance can take bias parameter α as keyword argument" begin
+@testitem "WeightedBalancedAcceptance can take bias parameter α as keyword argument" begin
     α = 3.14159
     wbas = WeightedBalancedAcceptance(; α = α)
     @test wbas.α == α
 end
 
-@testitem "BalancedAcceptance can take number of points as keyword argument" begin
+@testitem "WeightedBalancedAcceptance can take number of points as keyword argument" begin
     N = 40
     wbas = WeightedBalancedAcceptance(; numsites = N)
     @test wbas.numsites == N