-
Notifications
You must be signed in to change notification settings - Fork 0
/
tri_raster.jl
122 lines (107 loc) · 3.62 KB
/
tri_raster.jl
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
module TriRaster
export Triangle, rasterize, geometric_area, d, flatten, variance, average
using Random
"""
A wrapper for the 3-tuple of 2-tuples of integers.
This is probably unnecessary, but OOP habits die hard."""
struct Triangle
simplex::Tuple{Tuple{Int, Int}, Tuple{Int, Int}, Tuple{Int, Int}}
Triangle(simplex::Tuple{Tuple{Int, Int}, Tuple{Int, Int}, Tuple{Int, Int}}) = new(simplex)
Triangle((ax, ay), (bx, by), (cx, cy)) = new(((ax,ay),(bx,by),(cx,cy)))
Triangle(ax, ay, bx, by, cx, cy) = new(((ax,ay),(bx,by),(cx,cy)))
end
function Base.getindex(TRI::Triangle, i::Int64)
return TRI.simplex[i]
end
function Base.getindex(TRI::Triangle, i::Int, j::Int)
return TRI.simplex[i][j]
end
"Evaluate the trigonometric formula for the area of a triangle."
function geometric_area(TRI::Triangle)
a = d(TRI[1], TRI[2])
b = d(TRI[1], TRI[3])
γ = atan((TRI[2] .- TRI[1])...) - atan((TRI[3] .- TRI[1])...)
return 0.5 * abs(a*b*sin(γ))
end
# Utility functions
d(x1, x2) = sqrt(sum((x1 .- x2).^2))
flatten(P) = [(P...)...] # equivalently, `collect(Iterators.flatten(P))`
"Variance of `IMG` restricted to the sublattice spanned by `TRI`."
variance(IMG, TRI::Triangle) = var([IMG[i, j] for (i,j) in rasterize(TRI)])
"Average of `IMG` restricted to the sublattice spanned by `TRI`."
average(IMG, TRI::Triangle) = mean([IMG[i, j] for (i,j) in rasterize(TRI)])
"Iterate the sublattice of ℤ² intersected by `TRI`."
function rasterize(TRI)
# sort in vertically-ascending order
lo, md, hi = vertices = sort([v for v in TRI.simplex], by = v -> v[2])
x1, y1 = lo
x2, y2 = md
x3, y3 = hi
@assert y1 ≤ y2 ≤ y3
if y2 == y3
# bottom-flat triangle
return _rasterize_bf(TRI)
elseif y1 == y2
# top-flat triangle
return _rasterize_tf(TRI)
else
# decompose the triangle into top-flat and bottom-flat regions
y4 = y2
x4 = round(Int, x1 + ((y2 - y1) / (y3 - y1)) * (x3 - x1))
md2 = (x4, y4)
BTM = Triangle(lo, md, md2)
TOP = Triangle(md, md2, hi)
return union(_rasterize_bf(BTM), _rasterize_tf(TOP))
end
end
"Iterate the sublattice of ℤ² intersected by `TRI`. Case 1: `TRI` is bottom-flat."
function _rasterize_bf(TRI)
# sort in vertically-ascending order
lo, md, hi = vertices = sort([v for v in TRI.simplex], by = v -> v[2])
x1, y1 = lo
x2, y2 = md
x3, y3 = hi
@assert y1 < y2 == y3
# inverse slopes
invslope1 = (x2 - x1) / (y2 - y1)
invslope2 = (x3 - x1) / (y3 - y1)
# left, right bounds
lx = x1
rx = x1
# iterate lattice points
lines = Vector{Vector{Tuple{Int, Int}}}()
for y in y1:y2
minx = round(Int, min(lx, rx))
maxx = round(Int, max(lx, rx))
push!(lines, [(x, y) for x in minx:maxx])
lx += invslope1
rx += invslope2
end
return flatten(lines)
end
"Iterate the sublattice of ℤ² intersected by `TRI`. Case 2: `TRI` is top-flat."
function _rasterize_tf(TRI)
# sort in vertically-ascending order
lo, md, hi = vertices = sort([v for v in TRI.simplex], by = v -> v[2])
x1, y1 = lo
x2, y2 = md
x3, y3 = hi
@assert y1 == y2 < y3
# inverse slopes
invslope1 = (x3 - x1) / (y3 - y1)
invslope2 = (x3 - x2) / (y3 - y2)
# left, right bounds
lx = x3
rx = x3
# iterate lattice points
lines = Vector{Vector{Tuple{Int, Int}}}()
for y in Iterators.reverse(y1:y3)
minx = round(Int, min(lx, rx))
maxx = round(Int, max(lx, rx))
push!(lines, [(x, y) for x in minx:maxx])
lx -= invslope1
rx -= invslope2
end
return flatten(lines)
end
end