diff --git a/src/Meshes.jl b/src/Meshes.jl
index 0dc6a8ff8..69f1efc0c 100644
--- a/src/Meshes.jl
+++ b/src/Meshes.jl
@@ -41,7 +41,7 @@ import Distances: evaluate
 import NearestNeighbors: MinkowskiMetric
 
 # Transforms API
-import TransformsBase: Transform, →
+import TransformsBase: Transform, →, SequentialTransform, Identity
 import TransformsBase: isrevertible, isinvertible
 import TransformsBase: apply, revert, reapply, inverse
 import TransformsBase: parameters, preprocess
@@ -78,6 +78,9 @@ include("domains.jl")
 # utilities
 include("utils.jl")
 
+# transforms
+include("transforms.jl")
+
 # domain indices
 include("indices.jl")
 
@@ -120,9 +123,6 @@ include("discretization.jl")
 include("refinement.jl")
 include("coarsening.jl")
 
-# transforms
-include("transforms.jl")
-
 # miscellaneous
 include("rand.jl")
 include("distances.jl")
diff --git a/src/indices.jl b/src/indices.jl
index 722379bd7..ae5bae4bf 100644
--- a/src/indices.jl
+++ b/src/indices.jl
@@ -7,71 +7,95 @@
 # ---------------
 
 """
-    indices(domain, geometry)
+		indices(domain, geometry)
 
 Return the indices of the elements of the `domain` that intersect with the `geometry`.
 """
 indices(domain::Domain, geometry::Geometry) = findall(intersects(geometry), domain)
 
 function indices(grid::OrthoRegularGrid, point::Point)
-  point ∉ grid && return Int[]
+	point ∉ grid && return Int[]
 
-  # grid properties
-  orig = minimum(grid)
-  spac = spacing(grid)
-  dims = size(grid)
+	# grid properties
+	orig = minimum(grid)
+	spac = spacing(grid)
+	dims = size(grid)
 
-  # integer coordinates
-  coords = ceil.(Int, (point - orig) ./ spac)
+	# integer coordinates
+	coords = ceil.(Int, (point - orig) ./ spac)
 
-  # fix coordinates that are on the grid border
-  coords = clamp.(coords, 1, dims)
+	# fix coordinates that are on the grid border
+	coords = clamp.(coords, 1, dims)
 
-  # convert to linear index
-  [LinearIndices(dims)[coords...]]
+	# convert to linear index
+	[LinearIndices(dims)[coords...]]
 end
 
 function indices(grid::OrthoRegularGrid, chain::Chain)
-  dims = size(grid)
-  mask = falses(dims)
+	dims = size(grid)
+	mask = falses(dims)
 
-  for segment in segments(chain)
-    p₁, p₂ = vertices(segment)
-    _bresenham!(mask, grid, true, p₁, p₂)
-  end
+	for segment in segments(chain)
+		p₁, p₂ = vertices(segment)
+		_bresenham!(mask, grid, true, p₁, p₂)
+	end
 
-  LinearIndices(dims)[mask]
+	LinearIndices(dims)[mask]
 end
 
 function indices(grid::OrthoRegularGrid, poly::Polygon)
-  dims = size(grid)
-  mask = zeros(Int, dims)
-  cpoly = poly ∩ boundingbox(grid)
-  isnothing(cpoly) && return Int[]
+	dims = size(grid)
+	mask = zeros(Int, dims)
+	cpoly = poly ∩ boundingbox(grid)
+	isnothing(cpoly) && return Int[]
 
-  for (i, triangle) in enumerate(simplexify(cpoly))
-    _fill!(mask, grid, i, triangle)
-  end
+	for (i, triangle) in enumerate(simplexify(cpoly))
+		_fill!(mask, grid, i, triangle)
+	end
 
-  LinearIndices(dims)[mask .> 0]
+	LinearIndices(dims)[mask .> 0]
 end
 
 function indices(grid::OrthoRegularGrid, box::Box)
-  # cartesian range
-  range = cartesianrange(grid, box)
+	# cartesian range
+	range = cartesianrange(grid, box)
 
-  # convert to linear indices
-  LinearIndices(size(grid))[range] |> vec
+	# convert to linear indices
+	LinearIndices(size(grid))[range] |> vec
 end
 
 indices(grid::OrthoRegularGrid, multi::Multi) = mapreduce(geom -> indices(grid, geom), vcat, parent(multi)) |> unique
 
 function indices(grid::OrthoRectilinearGrid, box::Box)
-  # cartesian range
-  range = cartesianrange(grid, box)
+	# cartesian range
+	range = cartesianrange(grid, box)
 
-  # convert to linear indices
-  LinearIndices(size(grid))[range] |> vec
+	# convert to linear indices
+	LinearIndices(size(grid))[range] |> vec
+end
+
+# Fallback: inverse if possible, otherwise error
+function _reversetransform(t::Transform)
+	isrevertible(t) ? inverse(t) : throw(ArgumentError("Irreversible transform: $(typeof(t))"))
+end
+
+# Treat specific irreversible transforms as Identity (no-op)
+_reversetransform(::Morphological) = Identity()
+_reversetransform(::Repair) = Identity()
+_reversetransform(::Bridge) = Identity()
+
+function _reversetransform(st::SequentialTransform)
+	# Recursively reverse steps, replacing exceptions with Identity
+	newsteps = map(_reversetransform, st.transforms)
+	reduce(→, reverse(newsteps))
+end
+
+function indices(grid::TransformedGrid, geometry::Geometry)
+	# construct reverse transform from revertible steps
+	revtrans = _reversetransform(transform(grid))
+
+	# find indices in non-transformed space
+	indices(parent(grid), revtrans(geometry))
 end
 
 # ----------------
@@ -79,7 +103,7 @@ end
 # ----------------
 
 """
-    cartesianrange(grid, box)
+		cartesianrange(grid, box)
 
 Return the Cartesian range of the elements of the `grid` that intersect with the `box`.
 """
@@ -90,97 +114,97 @@ _manifoldrange(::Type{<:𝔼}, grid::Grid, box::Box) = _euclideanrange(grid, box
 _manifoldrange(::Type{<:🌐}, grid::Grid, box::Box) = _geodesicrange(grid, box)
 
 function _euclideanrange(grid::OrthoRegularGrid, box::Box)
-  # grid properties
-  or = minimum(grid)
-  sp = spacing(grid)
-  sz = size(grid)
+	# grid properties
+	or = minimum(grid)
+	sp = spacing(grid)
+	sz = size(grid)
 
-  # intersection of boxes
-  lo, up = extrema(boundingbox(grid) ∩ box)
+	# intersection of boxes
+	lo, up = extrema(boundingbox(grid) ∩ box)
 
-  # Cartesian indices of new corners
-  ijkₛ = max.(ceil.(Int, (lo - or) ./ sp), 1)
-  ijkₑ = min.(floor.(Int, (up - or) ./ sp) .+ 1, sz)
+	# Cartesian indices of new corners
+	ijkₛ = max.(ceil.(Int, (lo - or) ./ sp), 1)
+	ijkₑ = min.(floor.(Int, (up - or) ./ sp) .+ 1, sz)
 
-  # Cartesian range from corner to corner
-  CartesianIndex(Tuple(ijkₛ)):CartesianIndex(Tuple(ijkₑ))
+	# Cartesian range from corner to corner
+	CartesianIndex(Tuple(ijkₛ)):CartesianIndex(Tuple(ijkₑ))
 end
 
 function _euclideanrange(grid::OrthoRectilinearGrid, box::Box)
-  # grid properties
-  nd = paramdim(grid)
+	# grid properties
+	nd = paramdim(grid)
 
-  # intersection of boxes
-  lo, up = to.(extrema(boundingbox(grid) ∩ box))
+	# intersection of boxes
+	lo, up = to.(extrema(boundingbox(grid) ∩ box))
 
-  # integer coordinates of lower point
-  ijkₛ = ntuple(nd) do i
-    findlast(x -> x ≤ lo[i], xyz(grid)[i])
-  end
+	# integer coordinates of lower point
+	ijkₛ = ntuple(nd) do i
+		findlast(x -> x ≤ lo[i], xyz(grid)[i])
+	end
 
-  # integer coordinates of upper point
-  ijkₑ = ntuple(nd) do i
-    findfirst(x -> x ≥ up[i], xyz(grid)[i])
-  end
+	# integer coordinates of upper point
+	ijkₑ = ntuple(nd) do i
+		findfirst(x -> x ≥ up[i], xyz(grid)[i])
+	end
 
-  # integer coordinates of elements
-  CartesianIndex(ijkₛ):CartesianIndex(ijkₑ .- 1)
+	# integer coordinates of elements
+	CartesianIndex(ijkₛ):CartesianIndex(ijkₑ .- 1)
 end
 
 function _geodesicrange(grid::Grid, box::Box)
-  nlon, nlat = vsize(grid)
-
-  boxmin = convert(LatLon, coords(minimum(box)))
-  boxmax = convert(LatLon, coords(maximum(box)))
-
-  a = convert(LatLon, coords(vertex(grid, (1, 1))))
-  b = convert(LatLon, coords(vertex(grid, (nlon, 1))))
-  c = convert(LatLon, coords(vertex(grid, (1, nlat))))
-
-  swaplon = a.lon > b.lon
-  swaplat = a.lat > c.lat
-
-  loninds = swaplon ? (nlon:-1:1) : (1:1:nlon)
-  latinds = swaplat ? (nlat:-1:1) : (1:1:nlat)
-
-  gridlonₛ, gridlonₑ = swaplon ? (b.lon, a.lon) : (a.lon, b.lon)
-  gridlatₛ, gridlatₑ = swaplat ? (c.lat, a.lat) : (a.lat, c.lat)
-
-  lonmin = max(boxmin.lon, gridlonₛ)
-  latmin = max(boxmin.lat, gridlatₛ)
-  lonmax = min(boxmax.lon, gridlonₑ)
-  latmax = min(boxmax.lat, gridlatₑ)
-
-  iₛ = findlast(loninds) do i
-    p = vertex(grid, (i, 1))
-    c = convert(LatLon, coords(p))
-    c.lon ≤ lonmin
-  end
-  iₑ = findfirst(loninds) do i
-    p = vertex(grid, (i, 1))
-    c = convert(LatLon, coords(p))
-    c.lon ≥ lonmax
-  end
-
-  jₛ = findlast(latinds) do i
-    p = vertex(grid, (1, i))
-    c = convert(LatLon, coords(p))
-    c.lat ≤ latmin
-  end
-  jₑ = findfirst(latinds) do i
-    p = vertex(grid, (1, i))
-    c = convert(LatLon, coords(p))
-    c.lat ≥ latmax
-  end
-
-  if iₛ == iₑ || jₛ == jₑ
-    throw(ArgumentError("the passed limits are not valid for the grid"))
-  end
-
-  iₛ, iₑ = swaplon ? (iₑ, iₛ) : (iₛ, iₑ)
-  jₛ, jₑ = swaplat ? (jₑ, jₛ) : (jₛ, jₑ)
-
-  CartesianIndex(loninds[iₛ], latinds[jₛ]):CartesianIndex(loninds[iₑ] - 1, latinds[jₑ] - 1)
+	nlon, nlat = vsize(grid)
+
+	boxmin = convert(LatLon, coords(minimum(box)))
+	boxmax = convert(LatLon, coords(maximum(box)))
+
+	a = convert(LatLon, coords(vertex(grid, (1, 1))))
+	b = convert(LatLon, coords(vertex(grid, (nlon, 1))))
+	c = convert(LatLon, coords(vertex(grid, (1, nlat))))
+
+	swaplon = a.lon > b.lon
+	swaplat = a.lat > c.lat
+
+	loninds = swaplon ? (nlon:-1:1) : (1:1:nlon)
+	latinds = swaplat ? (nlat:-1:1) : (1:1:nlat)
+
+	gridlonₛ, gridlonₑ = swaplon ? (b.lon, a.lon) : (a.lon, b.lon)
+	gridlatₛ, gridlatₑ = swaplat ? (c.lat, a.lat) : (a.lat, c.lat)
+
+	lonmin = max(boxmin.lon, gridlonₛ)
+	latmin = max(boxmin.lat, gridlatₛ)
+	lonmax = min(boxmax.lon, gridlonₑ)
+	latmax = min(boxmax.lat, gridlatₑ)
+
+	iₛ = findlast(loninds) do i
+		p = vertex(grid, (i, 1))
+		c = convert(LatLon, coords(p))
+		c.lon ≤ lonmin
+	end
+	iₑ = findfirst(loninds) do i
+		p = vertex(grid, (i, 1))
+		c = convert(LatLon, coords(p))
+		c.lon ≥ lonmax
+	end
+
+	jₛ = findlast(latinds) do i
+		p = vertex(grid, (1, i))
+		c = convert(LatLon, coords(p))
+		c.lat ≤ latmin
+	end
+	jₑ = findfirst(latinds) do i
+		p = vertex(grid, (1, i))
+		c = convert(LatLon, coords(p))
+		c.lat ≥ latmax
+	end
+
+	if iₛ == iₑ || jₛ == jₑ
+		throw(ArgumentError("the passed limits are not valid for the grid"))
+	end
+
+	iₛ, iₑ = swaplon ? (iₑ, iₛ) : (iₛ, iₑ)
+	jₛ, jₑ = swaplat ? (jₑ, jₛ) : (jₛ, jₑ)
+
+	CartesianIndex(loninds[iₛ], latinds[jₛ]):CartesianIndex(loninds[iₑ]-1, latinds[jₑ]-1)
 end
 
 # -----------------
@@ -188,98 +212,98 @@ end
 # -----------------
 
 function _fill!(mask, grid, val, triangle)
-  v = vertices(triangle)
-
-  # fill edges of triangle
-  _bresenham!(mask, grid, val, v[1], v[2])
-  _bresenham!(mask, grid, val, v[2], v[3])
-  _bresenham!(mask, grid, val, v[3], v[1])
-
-  # fill interior of triangle
-  j₁ = findfirst(==(val), mask).I[2]
-  j₂ = findlast(==(val), mask).I[2]
-  for j in j₁:j₂
-    i₁ = findfirst(==(val), @view(mask[:, j]))
-    i₂ = findlast(==(val), @view(mask[:, j]))
-    mask[i₁:i₂, j] .= val
-  end
+	v = vertices(triangle)
+
+	# fill edges of triangle
+	_bresenham!(mask, grid, val, v[1], v[2])
+	_bresenham!(mask, grid, val, v[2], v[3])
+	_bresenham!(mask, grid, val, v[3], v[1])
+
+	# fill interior of triangle
+	j₁ = findfirst(==(val), mask).I[2]
+	j₂ = findlast(==(val), mask).I[2]
+	for j in j₁:j₂
+		i₁ = findfirst(==(val), @view(mask[:, j]))
+		i₂ = findlast(==(val), @view(mask[:, j]))
+		mask[i₁:i₂, j] .= val
+	end
 end
 
 # Bresenham's line algorithm: https://en.wikipedia.org/wiki/Bresenham's_line_algorithm
 function _bresenham!(mask, grid, val, p₁, p₂)
-  o = minimum(grid)
-  s = spacing(grid)
-
-  # integer coordinates
-  x₁, y₁ = ceil.(Int, (p₁ - o) ./ s)
-  x₂, y₂ = ceil.(Int, (p₂ - o) ./ s)
-
-  # fix coordinates of points that are on the grid border
-  xmax, ymax = size(grid)
-  x₁ = clamp(x₁, 1, xmax)
-  y₁ = clamp(y₁, 1, ymax)
-  x₂ = clamp(x₂, 1, xmax)
-  y₂ = clamp(y₂, 1, ymax)
-
-  if abs(y₂ - y₁) < abs(x₂ - x₁)
-    if x₁ > x₂
-      _bresenhamlow!(mask, val, x₂, y₂, x₁, y₁)
-    else
-      _bresenhamlow!(mask, val, x₁, y₁, x₂, y₂)
-    end
-  else
-    if y₁ > y₂
-      _bresenhamhigh!(mask, val, x₂, y₂, x₁, y₁)
-    else
-      _bresenhamhigh!(mask, val, x₁, y₁, x₂, y₂)
-    end
-  end
+	o = minimum(grid)
+	s = spacing(grid)
+
+	# integer coordinates
+	x₁, y₁ = ceil.(Int, (p₁ - o) ./ s)
+	x₂, y₂ = ceil.(Int, (p₂ - o) ./ s)
+
+	# fix coordinates of points that are on the grid border
+	xmax, ymax = size(grid)
+	x₁ = clamp(x₁, 1, xmax)
+	y₁ = clamp(y₁, 1, ymax)
+	x₂ = clamp(x₂, 1, xmax)
+	y₂ = clamp(y₂, 1, ymax)
+
+	if abs(y₂ - y₁) < abs(x₂ - x₁)
+		if x₁ > x₂
+			_bresenhamlow!(mask, val, x₂, y₂, x₁, y₁)
+		else
+			_bresenhamlow!(mask, val, x₁, y₁, x₂, y₂)
+		end
+	else
+		if y₁ > y₂
+			_bresenhamhigh!(mask, val, x₂, y₂, x₁, y₁)
+		else
+			_bresenhamhigh!(mask, val, x₁, y₁, x₂, y₂)
+		end
+	end
 end
 
 function _bresenhamlow!(mask, val, x₁, y₁, x₂, y₂)
-  dx = x₂ - x₁
-  dy = y₂ - y₁
-  yi = 1
-  if dy < 0
-    yi = -1
-    dy = -dy
-  end
-
-  D = 2dy - dx
-  y = y₁
-
-  for x in x₁:x₂
-    mask[x, y] = val
-
-    if D > 0
-      y = y + yi
-      D = D + 2dy - 2dx
-    else
-      D = D + 2dy
-    end
-  end
+	dx = x₂ - x₁
+	dy = y₂ - y₁
+	yi = 1
+	if dy < 0
+		yi = -1
+		dy = -dy
+	end
+
+	D = 2dy - dx
+	y = y₁
+
+	for x in x₁:x₂
+		mask[x, y] = val
+
+		if D > 0
+			y = y + yi
+			D = D + 2dy - 2dx
+		else
+			D = D + 2dy
+		end
+	end
 end
 
 function _bresenhamhigh!(mask, val, x₁, y₁, x₂, y₂)
-  dx = x₂ - x₁
-  dy = y₂ - y₁
-  xi = 1
-  if dx < 0
-    xi = -1
-    dx = -dx
-  end
-
-  D = 2dx - dy
-  x = x₁
-
-  for y in y₁:y₂
-    mask[x, y] = val
-
-    if D > 0
-      x = x + xi
-      D = D + 2dx - 2dy
-    else
-      D = D + 2dx
-    end
-  end
+	dx = x₂ - x₁
+	dy = y₂ - y₁
+	xi = 1
+	if dx < 0
+		xi = -1
+		dx = -dx
+	end
+
+	D = 2dx - dy
+	x = x₁
+
+	for y in y₁:y₂
+		mask[x, y] = val
+
+		if D > 0
+			x = x + xi
+			D = D + 2dx - 2dy
+		else
+			D = D + 2dx
+		end
+	end
 end
diff --git a/test/indices.jl b/test/indices.jl
index 6d5f2ccc4..ba4e6fad0 100644
--- a/test/indices.jl
+++ b/test/indices.jl
@@ -1,205 +1,206 @@
 @testitem "Indices" setup = [Setup] begin
-  g = cartgrid(10, 10)
-  b = Box(cart(1, 1), cart(5, 5))
-  v = view(g, b)
-  @test v == CartesianGrid(cart(0, 0), cart(6, 6), dims=(6, 6))
-
-  p = PointSet(vertices(g))
-  v = view(p, b)
-  @test centroid(v, 1) == cart(1, 1)
-  @test centroid(v, nelements(v)) == cart(5, 5)
-
-  # boxes
-  g = cartgrid(10, 10)
-  p = PointSet(vertices(g))
-  b = Ball(cart(0, 0), T(2))
-  v = view(g, b)
-  @test nelements(v) == 4
-  @test v[1] == g[1]
-  v = view(p, b)
-  @test nelements(v) == 6
-  @test to.(v) == vector.([(0, 0), (1, 0), (2, 0), (0, 1), (1, 1), (0, 2)])
-
-  b1 = Box(cart(0.6, 0.7), cart(1.0, 1.0))
-  b2 = Box(cart(0.6, 0.7), cart(1.2, 1.2))
-  b3 = Box(cart(0.0, 0.0), cart(0.6, 0.3))
-  b4 = Box(cart(-0.2, -0.2), cart(0.6, 0.3))
-  b5 = Box(cart(0.25, 0.15), cart(0.95, 0.85))
-  b6 = Box(cart(0.35, 0.25), cart(0.85, 0.75))
-  b7 = Box(cart(0.05, 0.05), cart(0.65, 0.35))
-  b8 = Box(cart(0.55, 0.65), cart(0.95, 0.95))
-  x = range(zero(T), stop=one(T), length=6)
-  y = T[0.0, 0.1, 0.3, 0.7, 0.9, 1.0]
-  g = RectilinearGrid(x, y)
-  linds = LinearIndices(size(g))
-  @test issetequal(indices(g, b1), [linds[4, 4], linds[5, 4], linds[4, 5], linds[5, 5]])
-  @test issetequal(indices(g, b2), [linds[4, 4], linds[5, 4], linds[4, 5], linds[5, 5]])
-  @test issetequal(indices(g, b3), [linds[1, 1], linds[2, 1], linds[3, 1], linds[1, 2], linds[2, 2], linds[3, 2]])
-  @test issetequal(indices(g, b4), [linds[1, 1], linds[2, 1], linds[3, 1], linds[1, 2], linds[2, 2], linds[3, 2]])
-  @test issetequal(
-    indices(g, b5),
-    [
-      linds[2, 2],
-      linds[3, 2],
-      linds[4, 2],
-      linds[5, 2],
-      linds[2, 3],
-      linds[3, 3],
-      linds[4, 3],
-      linds[5, 3],
-      linds[2, 4],
-      linds[3, 4],
-      linds[4, 4],
-      linds[5, 4]
-    ]
-  )
-  @test issetequal(
-    indices(g, b6),
-    [
-      linds[2, 2],
-      linds[3, 2],
-      linds[4, 2],
-      linds[5, 2],
-      linds[2, 3],
-      linds[3, 3],
-      linds[4, 3],
-      linds[5, 3],
-      linds[2, 4],
-      linds[3, 4],
-      linds[4, 4],
-      linds[5, 4]
-    ]
-  )
-  @test issetequal(
-    indices(g, b7),
-    [
-      linds[1, 1],
-      linds[2, 1],
-      linds[3, 1],
-      linds[4, 1],
-      linds[1, 2],
-      linds[2, 2],
-      linds[3, 2],
-      linds[4, 2],
-      linds[1, 3],
-      linds[2, 3],
-      linds[3, 3],
-      linds[4, 3]
-    ]
-  )
-  @test issetequal(
-    indices(g, b8),
-    [
-      linds[3, 3],
-      linds[4, 3],
-      linds[5, 3],
-      linds[3, 4],
-      linds[4, 4],
-      linds[5, 4],
-      linds[3, 5],
-      linds[4, 5],
-      linds[5, 5]
-    ]
-  )
-
-  # convex polygons
-  tri = Triangle(cart(5, 7), cart(10, 12), cart(15, 7))
-  pent = Pentagon(cart(6, 1), cart(2, 10), cart(10, 16), cart(18, 10), cart(14, 1))
-
-  grid = cartgrid(20, 20)
-  linds = LinearIndices(size(grid))
-  @test linds[10, 10] ∈ indices(grid, tri)
-  @test linds[10, 6] ∈ indices(grid, pent)
-
-  grid = CartesianGrid(cart(-2, -2), cart(20, 20), T.((0.5, 1.5)))
-  linds = LinearIndices(size(grid))
-  @test linds[21, 7] ∈ indices(grid, tri)
-  @test linds[21, 4] ∈ indices(grid, pent)
-
-  grid = CartesianGrid(cart(-100, -100), cart(20, 20), T.((2, 2)))
-  linds = LinearIndices(size(grid))
-  @test linds[57, 54] ∈ indices(grid, tri)
-  @test linds[55, 53] ∈ indices(grid, pent)
-
-  # non-convex polygons
-  poly1 = PolyArea(cart.([(3, 3), (9, 9), (3, 15), (17, 15), (17, 3)]))
-  poly2 = PolyArea([pointify(pent), pointify(tri)])
-
-  grid = cartgrid(20, 20)
-  linds = LinearIndices(size(grid))
-  @test linds[12, 6] ∈ indices(grid, poly1)
-  @test linds[10, 3] ∈ indices(grid, poly2)
-
-  grid = CartesianGrid(cart(-2, -2), cart(20, 20), T.((0.5, 1.5)))
-  linds = LinearIndices(size(grid))
-  @test linds[22, 6] ∈ indices(grid, poly1)
-  @test linds[17, 4] ∈ indices(grid, poly2)
-
-  grid = CartesianGrid(cart(-100, -100), cart(20, 20), T.((2, 2)))
-  linds = LinearIndices(size(grid))
-  @test linds[57, 54] ∈ indices(grid, poly1)
-  @test linds[55, 53] ∈ indices(grid, poly2)
-
-  # rotate
-  poly1 = poly1 |> Rotate(Angle2d(T(π / 2)))
-  poly2 = poly2 |> Rotate(Angle2d(T(π / 2)))
-
-  grid = CartesianGrid(cart(-20, 0), cart(0, 20), T.((1, 1)))
-  linds = LinearIndices(size(grid))
-  @test linds[12, 12] ∈ indices(grid, poly1)
-  @test linds[16, 11] ∈ indices(grid, poly2)
-
-  grid = CartesianGrid(cart(-22, -2), cart(0, 20), T.((0.5, 1.5)))
-  linds = LinearIndices(size(grid))
-  @test linds[26, 8] ∈ indices(grid, poly1)
-  @test linds[36, 9] ∈ indices(grid, poly2)
-
-  grid = CartesianGrid(cart(-100, -100), cart(20, 20), T.((2, 2)))
-  linds = LinearIndices(size(grid))
-  @test linds[46, 57] ∈ indices(grid, poly1)
-  @test linds[48, 55] ∈ indices(grid, poly2)
-
-  # multi
-  multi = Multi([tri, pent])
-  grid = cartgrid(20, 20)
-  linds = LinearIndices(size(grid))
-  @test linds[10, 10] ∈ indices(grid, multi)
-  @test linds[10, 6] ∈ indices(grid, multi)
-
-  # clipping
-  tri = Triangle(cart(-4, 10), cart(5, 19), cart(5, 1))
-  grid = cartgrid(20, 20)
-  linds = LinearIndices(size(grid))
-  @test linds[3, 10] ∈ indices(grid, tri)
-
-  # out of grid
-  tri = Triangle(cart(-12, 8), cart(-8, 14), cart(-4, 8))
-  grid = cartgrid(20, 20)
-  @test isempty(indices(grid, tri))
-
-  # chain
-  seg = Segment(cart(2, 12), cart(16, 18))
-  rope = Rope(cart(8, 1), cart(5, 9), cart(9, 13), cart(17, 10))
-  ring = Ring(cart(8, 1), cart(5, 9), cart(9, 13), cart(17, 10))
-  grid = cartgrid(20, 20)
-  linds = LinearIndices(size(grid))
-  @test linds[9, 15] ∈ indices(grid, seg)
-  @test linds[7, 11] ∈ indices(grid, rope)
-  @test linds[12, 5] ∈ indices(grid, ring)
-
-  # points
-  p1 = cart(0, 0)
-  p2 = cart(0.5, 0.5)
-  p3 = cart(1, 1)
-  p4 = cart(2, 2)
-  p5 = cart(10, 10)
-  p6 = cart(11, 11)
-  grid = cartgrid(10, 10)
-  linds = LinearIndices(size(grid))
-  @test linds[1, 1] == only(indices(grid, p1))
-  @test linds[1, 1] == only(indices(grid, p2))
-  @test linds[1, 1] == only(indices(grid, p3))
-  @test linds[2, 2] == only(indices(grid, p4))
-  @test linds[10, 10] == only(indices(grid, p5))
-  @test isempty(indices(grid, p6))
+	g = cartgrid(10, 10)
+	b = Box(cart(1, 1), cart(5, 5))
+	v = view(g, b)
+	@test v == CartesianGrid(cart(0, 0), cart(6, 6), dims=(6, 6))
+
+	p = PointSet(vertices(g))
+	v = view(p, b)
+	@test centroid(v, 1) == cart(1, 1)
+	@test centroid(v, nelements(v)) == cart(5, 5)
+
+	# boxes
+	g = cartgrid(10, 10)
+	p = PointSet(vertices(g))
+	b = Ball(cart(0, 0), T(2))
+	v = view(g, b)
+	@test nelements(v) == 4
+	@test v[1] == g[1]
+	v = view(p, b)
+	@test nelements(v) == 6
+	@test to.(v) == vector.([(0, 0), (1, 0), (2, 0), (0, 1), (1, 1), (0, 2)])
+
+	b1 = Box(cart(0.6, 0.7), cart(1.0, 1.0))
+	b2 = Box(cart(0.6, 0.7), cart(1.2, 1.2))
+	b3 = Box(cart(0.0, 0.0), cart(0.6, 0.3))
+	b4 = Box(cart(-0.2, -0.2), cart(0.6, 0.3))
+	b5 = Box(cart(0.25, 0.15), cart(0.95, 0.85))
+	b6 = Box(cart(0.35, 0.25), cart(0.85, 0.75))
+	b7 = Box(cart(0.05, 0.05), cart(0.65, 0.35))
+	b8 = Box(cart(0.55, 0.65), cart(0.95, 0.95))
+	x = range(zero(T), stop=one(T), length=6)
+	y = T[0.0, 0.1, 0.3, 0.7, 0.9, 1.0]
+	g = RectilinearGrid(x, y)
+	linds = LinearIndices(size(g))
+	@test issetequal(indices(g, b1), [linds[4, 4], linds[5, 4], linds[4, 5], linds[5, 5]])
+	@test issetequal(indices(g, b2), [linds[4, 4], linds[5, 4], linds[4, 5], linds[5, 5]])
+	@test issetequal(indices(g, b3), [linds[1, 1], linds[2, 1], linds[3, 1], linds[1, 2], linds[2, 2], linds[3, 2]])
+	@test issetequal(indices(g, b4), [linds[1, 1], linds[2, 1], linds[3, 1], linds[1, 2], linds[2, 2], linds[3, 2]])
+	@test issetequal(
+		indices(g, b5),
+		[
+			linds[2, 2],
+			linds[3, 2],
+			linds[4, 2],
+			linds[5, 2],
+			linds[2, 3],
+			linds[3, 3],
+			linds[4, 3],
+			linds[5, 3],
+			linds[2, 4],
+			linds[3, 4],
+			linds[4, 4],
+			linds[5, 4],
+		],
+	)
+	@test issetequal(
+		indices(g, b6),
+		[
+			linds[2, 2],
+			linds[3, 2],
+			linds[4, 2],
+			linds[5, 2],
+			linds[2, 3],
+			linds[3, 3],
+			linds[4, 3],
+			linds[5, 3],
+			linds[2, 4],
+			linds[3, 4],
+			linds[4, 4],
+			linds[5, 4],
+		],
+	)
+	@test issetequal(
+		indices(g, b7),
+		[
+			linds[1, 1],
+			linds[2, 1],
+			linds[3, 1],
+			linds[4, 1],
+			linds[1, 2],
+			linds[2, 2],
+			linds[3, 2],
+			linds[4, 2],
+			linds[1, 3],
+			linds[2, 3],
+			linds[3, 3],
+			linds[4, 3],
+		],
+	)
+	@test issetequal(
+		indices(g, b8),
+		[
+			linds[3, 3],
+			linds[4, 3],
+			linds[5, 3],
+			linds[3, 4],
+			linds[4, 4],
+			linds[5, 4],
+			linds[3, 5],
+			linds[4, 5],
+			linds[5, 5],
+		],
+	)
+
+	# convex polygons
+	tri = Triangle(cart(5, 7), cart(10, 12), cart(15, 7))
+	pent = Pentagon(cart(6, 1), cart(2, 10), cart(10, 16), cart(18, 10), cart(14, 1))
+
+	grid = cartgrid(20, 20)
+	linds = LinearIndices(size(grid))
+	@test linds[10, 10] ∈ indices(grid, tri)
+	@test linds[10, 6] ∈ indices(grid, pent)
+
+	grid = CartesianGrid(cart(-2, -2), cart(20, 20), T.((0.5, 1.5)))
+	linds = LinearIndices(size(grid))
+	@test linds[21, 7] ∈ indices(grid, tri)
+	@test linds[21, 4] ∈ indices(grid, pent)
+
+	grid = CartesianGrid(cart(-100, -100), cart(20, 20), T.((2, 2)))
+	linds = LinearIndices(size(grid))
+	@test linds[57, 54] ∈ indices(grid, tri)
+	@test linds[55, 53] ∈ indices(grid, pent)
+
+	# non-convex polygons
+	poly1 = PolyArea(cart.([(3, 3), (9, 9), (3, 15), (17, 15), (17, 3)]))
+	poly2 = PolyArea([pointify(pent), pointify(tri)])
+
+	grid = cartgrid(20, 20)
+	linds = LinearIndices(size(grid))
+	@test linds[12, 6] ∈ indices(grid, poly1)
+	@test linds[10, 3] ∈ indices(grid, poly2)
+
+	grid = CartesianGrid(cart(-2, -2), cart(20, 20), T.((0.5, 1.5)))
+	linds = LinearIndices(size(grid))
+	@test linds[22, 6] ∈ indices(grid, poly1)
+	@test linds[17, 4] ∈ indices(grid, poly2)
+
+	grid = CartesianGrid(cart(-100, -100), cart(20, 20), T.((2, 2)))
+	linds = LinearIndices(size(grid))
+	@test linds[57, 54] ∈ indices(grid, poly1)
+	@test linds[55, 53] ∈ indices(grid, poly2)
+
+	# rotate
+	poly1 = poly1 |> Rotate(Angle2d(T(π / 2)))
+	poly2 = poly2 |> Rotate(Angle2d(T(π / 2)))
+
+	grid = CartesianGrid(cart(-20, 0), cart(0, 20), T.((1, 1)))
+	linds = LinearIndices(size(grid))
+	@test linds[12, 12] ∈ indices(grid, poly1)
+	@test linds[16, 11] ∈ indices(grid, poly2)
+
+	grid = CartesianGrid(cart(-22, -2), cart(0, 20), T.((0.5, 1.5)))
+	linds = LinearIndices(size(grid))
+	@test linds[26, 8] ∈ indices(grid, poly1)
+	@test linds[36, 9] ∈ indices(grid, poly2)
+
+	grid = CartesianGrid(cart(-100, -100), cart(20, 20), T.((2, 2)))
+	linds = LinearIndices(size(grid))
+	@test linds[46, 57] ∈ indices(grid, poly1)
+	@test linds[48, 55] ∈ indices(grid, poly2)
+
+	# multi
+	multi = Multi([tri, pent])
+	grid = cartgrid(20, 20)
+	linds = LinearIndices(size(grid))
+	@test linds[10, 10] ∈ indices(grid, multi)
+	@test linds[10, 6] ∈ indices(grid, multi)
+
+	# clipping
+	tri = Triangle(cart(-4, 10), cart(5, 19), cart(5, 1))
+	grid = cartgrid(20, 20)
+	linds = LinearIndices(size(grid))
+	@test linds[3, 10] ∈ indices(grid, tri)
+
+	# out of grid
+	tri = Triangle(cart(-12, 8), cart(-8, 14), cart(-4, 8))
+	grid = cartgrid(20, 20)
+	@test isempty(indices(grid, tri))
+
+	# chain
+	seg = Segment(cart(2, 12), cart(16, 18))
+	rope = Rope(cart(8, 1), cart(5, 9), cart(9, 13), cart(17, 10))
+	ring = Ring(cart(8, 1), cart(5, 9), cart(9, 13), cart(17, 10))
+	grid = cartgrid(20, 20)
+	linds = LinearIndices(size(grid))
+	@test linds[9, 15] ∈ indices(grid, seg)
+	@test linds[7, 11] ∈ indices(grid, rope)
+	@test linds[12, 5] ∈ indices(grid, ring)
+
+	# points
+	p1 = cart(0, 0)
+	p2 = cart(0.5, 0.5)
+	p3 = cart(1, 1)
+	p4 = cart(2, 2)
+	p5 = cart(10, 10)
+	p6 = cart(11, 11)
+	grid = cartgrid(10, 10)
+	linds = LinearIndices(size(grid))
+	@test linds[1, 1] == only(indices(grid, p1))
+	@test linds[1, 1] == only(indices(grid, p2))
+	@test linds[1, 1] == only(indices(grid, p3))
+	@test linds[2, 2] == only(indices(grid, p4))
+	@test linds[10, 10] == only(indices(grid, p5))
+	@test isempty(indices(grid, p6))
+
 end