Skip to content

Commit 6fc8e45

Browse files
authored
Merge branch 'main' into optimize_com
2 parents f591550 + 02a35cd commit 6fc8e45

9 files changed

Lines changed: 173 additions & 76 deletions

File tree

Project.toml

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -6,20 +6,16 @@ version = "0.2.0"
66
projects = ["test", "docs", "benchmark"]
77

88
[deps]
9-
Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
109
ConsensusFitting = "463569c3-d6e9-43b4-bb3a-bd02a13deb89"
1110
CoordinateTransformations = "150eb455-5306-5404-9cee-2592286d6298"
12-
Distances = "b4f34e82-e78d-54a5-968a-f98e89d6e8f7"
1311
ImageTransformations = "02fcd773-0e25-5acc-982a-7f6622650795"
1412
NearestNeighbors = "b8a86587-4115-5ab1-83bc-aa920d37bbce"
1513
Photometry = "af68cb61-81ac-52ed-8703-edc140936be4"
1614
TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"
1715

1816
[compat]
19-
Combinatorics = "1"
2017
ConsensusFitting = "1"
2118
CoordinateTransformations = "0.6"
22-
Distances = "0.10"
2319
ImageTransformations = "0.10"
2420
NearestNeighbors = "0.4"
2521
Photometry = "0.9"

benchmark/benchmarks.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@ const SUITE = BenchmarkGroup()
77
SUITE["core"] = BenchmarkGroup()
88

99
SUITE["core"]["_triangle_invariants"] = @benchmarkable _triangle_invariants(phot) setup=(phot = Table(xcenter = rand(Float64, 100), ycenter = rand(Float64, 100)))
10-
SUITE["core"]["_build_correspondences"] = @benchmarkable _build_correspondences(C_from, ℳ_from, C_to, ℳ_to) setup=begin
10+
SUITE["core"]["_build_correspondences"] = @benchmarkable _build_correspondences(C_from, ℳ_from, phot_from, C_to, ℳ_to, phot_to) setup=begin
1111
phot = Table(xcenter = rand(Float64, 100), ycenter = rand(Float64, 100))
1212
C_from, ℳ_from = _triangle_invariants(phot)
1313
C_to, ℳ_to = _triangle_invariants(phot)

docs/src/walkthrough.md

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -345,7 +345,7 @@ end
345345
We next build our list of candidate correspondences in this invariant space via a nearest neighbors search:
346346

347347
```@example walkthrough
348-
correspondences = Astroalign._build_correspondences(C_from, ℳ_from, C_to, ℳ_to)
348+
correspondences = Astroalign._build_correspondences(C_from, ℳ_from, phot_from, C_to, ℳ_to, phot_to)
349349
350350
println("Candidate triangle matches: $(size(correspondences, 4))")
351351
```

src/Astroalign.jl

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,9 +1,7 @@
11
module Astroalign
22

3-
using Combinatorics: combinations
43
using ConsensusFitting: ransac
54
using CoordinateTransformations: kabsch, AffineMap
6-
using Distances: euclidean
75
using ImageTransformations: warp
86
using NearestNeighbors: nn, KDTree
97
using Photometry:

src/register.jl

Lines changed: 150 additions & 61 deletions
Original file line numberDiff line numberDiff line change
@@ -4,61 +4,148 @@
44
Returns all combinations (``C``) of three candidate point sources from the table of sources `phot` returned by [`Photometry.Aperture.photometry`](@extref), and the computed invariant ``\\mathscr M`` for each according to Eq. 3 from [_Beroiz, M., Cabral, J. B., & Sanchez, B. (2020)_](https://ui.adsabs.harvard.edu/abs/2020A%26C....3200384B/abstract).
55
"""
66
function _triangle_invariants(phot)
7-
C = combinations(phot, 3)
8-
= map(C) do (pa, pb, pc)
9-
a, b, c = (
10-
(pa.ycenter, pa.xcenter),
11-
(pb.ycenter, pb.xcenter),
12-
(pc.ycenter, pc.xcenter),
13-
)
14-
Ls = sort!([euclidean(a, b), euclidean(b, c), euclidean(a, c)])
15-
(Ls[3] / Ls[2], Ls[2] / Ls[1])
16-
end |> stack
7+
@inline function sort3(a, b, c)
8+
a, b = minmax(a, b)
9+
b, c = minmax(b, c)
10+
a, b = minmax(a, b)
11+
return a, b, c # sorted ascending: small, mid, large
12+
end
13+
xs = phot.xcenter
14+
ys = phot.ycenter
15+
n = length(phot)
16+
ntriangles = binomial(n, 3)
17+
= Matrix{Float64}(undef, 2, ntriangles)
18+
C = Vector{NTuple{3,Int}}(undef, ntriangles)
19+
20+
# Enumerate all combinations(1:n,3) triangles via nested loops over strictly increasing index triples
21+
# (i < j < l), storing results in pre-allocated C and ℳ indexed by the flat counter k.
22+
k = 0
23+
for i in 1:n-2
24+
xi, yi = xs[i], ys[i]
25+
for j in i+1:n-1
26+
xj, yj = xs[j], ys[j]
27+
dx_ji, dy_ji = xj - xi, yj - yi
28+
d2_ji = dx_ji^2 + dy_ji^2 # squared distance between sources i and j
29+
for l in j+1:n
30+
k += 1
31+
xl, yl = xs[l], ys[l]
32+
dx_lj, dy_lj = xl - xj, yl - yj
33+
dx_li, dy_li = xl - xi, yl - yi
34+
d2_lj = dx_lj^2 + dy_lj^2 # squared distance between sources j and l
35+
d2_li = dx_li^2 + dy_li^2 # squared distance between sources i and l
36+
37+
L1, L2, L3 = sort3(d2_ji, d2_lj, d2_li)
38+
ℳ[1, k] = sqrt(L3 / L2)
39+
ℳ[2, k] = sqrt(L2 / L1)
40+
C[k] = _canonical_vertex_order(i, j, l, d2_ji, d2_lj, d2_li, xs, ys)
41+
end
42+
end
43+
end
44+
1745
return C, ℳ
1846
end
1947

20-
"""
21-
_canonical_vertex_order(pa, pb, pc)
48+
@doc raw"""
49+
_canonical_vertex_order(i, j, l, d2_ji, d2_lj, d2_li, xs, ys)
2250
23-
Re-order the three vertices of a triangle so that:
51+
Re-order the three vertices of a triangle with `x` and `y` coordinates
52+
`x = (xs[i], xs[j], xs[l])` and `y = (ys[i], ys[j], ys[l])`
53+
into a canonical form so that:
2454
2555
1. The apex (vertex opposite the longest edge) is last.
26-
1. The two base vertices are ordered counter-clockwise (positive cross product).
56+
2. The two base vertices are ordered counter-clockwise (positive cross product).
2757
28-
This canonical form is preserved under rotation and translation, so corresponding triangles in two images receive the same vertex permutation and produce geometrically consistent point correspondences.
29-
"""
30-
function _canonical_vertex_order(pa, pb, pc)
31-
xa, ya = pa.xcenter, pa.ycenter
32-
xb, yb = pb.xcenter, pb.ycenter
33-
xc, yc = pc.xcenter, pc.ycenter
34-
35-
d2_ab = (xb - xa)^2 + (yb - ya)^2
36-
d2_bc = (xc - xb)^2 + (yc - yb)^2
37-
d2_ac = (xc - xa)^2 + (yc - ya)^2
38-
39-
# Identify the two base vertices (endpoints of longest edge) and the apex
40-
if d2_ab >= d2_bc && d2_ab >= d2_ac
41-
v1, v2, apex = pa, pb, pc
42-
elseif d2_bc >= d2_ab && d2_bc >= d2_ac
43-
v1, v2, apex = pb, pc, pa
44-
else
45-
v1, v2, apex = pa, pc, pb
46-
end
58+
This canonical form is preserved under rotation and translation, so corresponding triangles
59+
in two images receive the same vertex permutation and produce geometrically consistent point
60+
correspondences.
61+
62+
# Examples
63+
64+
Consider the triangle with vertices
65+
66+
```text
67+
(1) = (0,4)
68+
|\
69+
| \
70+
4 | \ 5
71+
| \
72+
| \
73+
(0,0) = (2)----(3) = (3,0)
74+
3
75+
```
76+
77+
The longest edge is (1)-(3), so the apex is vertex (2).
78+
The initial canonical ordering is therefore (1,3,2).
79+
80+
However, this ordering is clockwise in image coordinates, so the
81+
base vertices are swapped to enforce counter-clockwise winding,
82+
yielding the final canonical ordering (3,1,2) for these vertices.
4783
48-
# Enforce CCW winding so that a rotation does not change the order
49-
cross = (v2.xcenter - v1.xcenter) * (apex.ycenter - v1.ycenter) -
50-
(v2.ycenter - v1.ycenter) * (apex.xcenter - v1.xcenter)
51-
cross < 0 && ((v1, v2) = (v2, v1))
84+
```jldoctest
85+
julia> using Astroalign: _canonical_vertex_order
5286
53-
return (v1, v2, apex)
87+
julia> i, x_i, y_i = 1, 0.0, 4.0;
88+
89+
julia> j, x_j, y_j = 2, 0.0, 0.0;
90+
91+
julia> l, x_l, y_l = 3, 3.0, 0.0;
92+
93+
julia> d2_ji = (x_j - x_i)^2 + (y_j - y_i)^2
94+
16.0
95+
96+
julia> d2_lj = (x_l - x_j)^2 + (y_l - y_j)^2
97+
9.0
98+
99+
julia> d2_li = (x_l - x_i)^2 + (y_l - y_i)^2
100+
25.0
101+
102+
julia> _canonical_vertex_order(i, j, l, d2_ji, d2_lj, d2_li, [x_i, x_j, x_l], [y_i, y_j, y_l])
103+
(3, 1, 2)
104+
```
105+
106+
# Notes
107+
108+
The squared edge lengths `d2_ji`, `d2_lj`, `d2_li` are accepted as arguments rather than
109+
recomputed, since they are already available at the call site in [`_triangle_invariants`](@ref).
110+
"""
111+
@inline function _canonical_vertex_order(i, j, l, d2_ji, d2_lj, d2_li, xs, ys)
112+
# There are three possible longest edges:
113+
# 1. ji is longest (c1 == true)
114+
# 2. lj is longest (c2 == true)
115+
# 3. li is longest (neither c1 nor c2)
116+
ji_longest = (d2_ji >= d2_lj) & (d2_ji >= d2_li)
117+
lj_longest = (d2_lj >= d2_ji) & (d2_lj >= d2_li)
118+
119+
# The nested ifelse expressions are equivalent to:
120+
#
121+
# if c1
122+
# v1, v2, apex = i, j, l
123+
# elseif c2
124+
# v1, v2, apex = j, l, i
125+
# else
126+
# v1, v2, apex = i, l, j
127+
# end
128+
# Using ifelse keeps this branchless
129+
v1 = ifelse(ji_longest, i, ifelse(lj_longest, j, i))
130+
v2 = ifelse(ji_longest, j, ifelse(lj_longest, l, l))
131+
apex = ifelse(ji_longest, l, ifelse(lj_longest, i, j))
132+
133+
# Enforce CCW winding: swap base vertices if cross product is negative
134+
@inbounds cross = (xs[v2] - xs[v1]) * (ys[apex] - ys[v1]) - (ys[v2] - ys[v1]) * (xs[apex] - xs[v1])
135+
swap = cross < 0
136+
vv1 = ifelse(swap, v2, v1)
137+
vv2 = ifelse(swap, v1, v2)
138+
return vv1, vv2, apex
54139
end
55140

56141
"""
57-
_build_correspondences(C_from, ℳ_from, C_to, ℳ_to)
142+
_build_correspondences(C_from, ℳ_from, phot_from, C_to, ℳ_to, phot_to)
58143
59144
Build a `2 × 3 × 2 × N` array of candidate triangle-level correspondences
60145
between the `from` and `to` frames. The `C` and `ℳ` are the combinations of three
61-
points and their invariants as returned by [`Astroalign._triangle_invariants`](@ref).
146+
points and their invariants as returned by [`_triangle_invariants`](@ref), which
147+
guarantees that both `C_from` and `C_to` are already in canonical vertex order
148+
(base vertices CCW, apex last).
62149
63150
The axes are `[coord, vertex, frame, match]`:
64151
@@ -69,29 +156,31 @@ The axes are `[coord, vertex, frame, match]`:
69156
70157
So `out[:, v, 1, n]` is the `(x, y)` position of vertex `v` in the `from`
71158
frame for match `n`, and `out[:, v, 2, n]` is the corresponding position in
72-
the `to` frame. Vertices are ordered canonically via
73-
[`_canonical_vertex_order`](@ref), so corresponding triangles receive the
74-
same geometric vertex assignment.
159+
the `to` frame. Because both frames share the same canonical vertex ordering,
160+
corresponding vertices across frames are geometrically consistent and suitable
161+
for direct use in transform estimation.
75162
"""
76-
function _build_correspondences(C_from, ℳ_from, C_to, ℳ_to)
77-
C_from_list = collect(C_from)
78-
C_to_list = collect(C_to)
79-
80-
(isempty(C_from_list) || isempty(C_to_list)) && return zeros(2, 3, 2, 0)
81-
82-
# Get most similar triangle in to-frame for each from-frame triangle, using the invariants as features
83-
idxs, _ = nn(KDTree(ℳ_to), ℳ_from)
84-
85-
out = Array{Float64}(undef, 2, 3, 2, length(C_from_list))
86-
87-
for i in eachindex(C_from_list)
88-
canon_from = _canonical_vertex_order(C_from_list[i]...)
89-
canon_to = _canonical_vertex_order(C_to_list[idxs[i]]...)
163+
function _build_correspondences(C_from, ℳ_from, phot_from, C_to, ℳ_to, phot_to)
164+
(isempty(C_from) || isempty(C_to)) && return zeros(2, 3, 2, 0)
165+
166+
# Build the KD-tree once; query per-triangle inside the loop to avoid
167+
# materialising the full idxs vector
168+
tree = KDTree(ℳ_to)
169+
170+
out = Array{Float64}(undef, 2, 3, 2, length(C_from))
171+
xs_from, ys_from = phot_from.xcenter, phot_from.ycenter
172+
xs_to, ys_to = phot_to.xcenter, phot_to.ycenter
173+
174+
for i in eachindex(C_from)
175+
# Find nearest neighbor of ℳ_from[:, i] in ℳ_to, returning the index of the corresponding triangle in the to-frame
176+
idx, _ = nn(tree, view(ℳ_from, :, i))
177+
C_from_i = C_from[i]
178+
C_to_i = C_to[idx]
90179
for v in 1:3
91-
out[1, v, 1, i] = canon_from[v].xcenter
92-
out[2, v, 1, i] = canon_from[v].ycenter
93-
out[1, v, 2, i] = canon_to[v].xcenter
94-
out[2, v, 2, i] = canon_to[v].ycenter
180+
out[1, v, 1, i] = xs_from[C_from_i[v]]
181+
out[2, v, 1, i] = ys_from[C_from_i[v]]
182+
out[1, v, 2, i] = xs_to[C_to_i[v]]
183+
out[2, v, 2, i] = ys_to[C_to_i[v]]
95184
end
96185
end
97186

src/warp.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -165,7 +165,7 @@ function find_transform(img_from, img_to;
165165
C_to, ℳ_to = _triangle_invariants(phot_to)
166166

167167
# Step 3: Build candidate correspondence pool via nearest neighbors triangle matching
168-
correspondences = _build_correspondences(C_from, ℳ_from, C_to, ℳ_to)
168+
correspondences = _build_correspondences(C_from, ℳ_from, phot_from, C_to, ℳ_to, phot_to)
169169

170170
size(correspondences, 4) < 1 &&
171171
error("align_frames: not enough candidate correspondences ($(size(correspondences, 4))); " *

test/Project.toml

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,3 +8,4 @@ ParallelTestRunner = "d3525ed8-44d0-4b2c-a655-542cee43accc"
88
Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
99
StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
1010
Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
11+
TypedTables = "9d95f2ec-7b3d-5a63-8d20-e2491e220bb9"

test/runtests.jl

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -47,11 +47,11 @@ const init_code = quote
4747
),
4848

4949
combinations_to = [
50-
[(xcenter = 0, ycenter = 0), (xcenter = 0, ycenter = 4), (xcenter = 3, ycenter = 0)],
50+
(1, 2, 3),
5151
],
5252

5353
combinations_from = [
54-
[(xcenter = 0, ycenter = -3), (xcenter = 0, ycenter = 1), (xcenter = 3, ycenter = -3)],
54+
(1, 2, 3),
5555
],
5656

5757
invariants = [

test/test-functions.jl

Lines changed: 17 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -1,14 +1,27 @@
11
@testset "triangle_invariants" begin
2-
using Astroalign: _triangle_invariants
2+
using Astroalign: _triangle_invariants, _canonical_vertex_order
3+
using TypedTables: Table
4+
5+
# This signature computes the square distances
6+
function _canonical_vertex_order_test(i, j, l, xs, ys)
7+
d2_ji = (xs[j] - xs[i])^2 + (ys[j] - ys[i])^2
8+
d2_lj = (xs[l] - xs[j])^2 + (ys[l] - ys[j])^2
9+
d2_li = (xs[l] - xs[i])^2 + (ys[l] - ys[i])^2
10+
return _canonical_vertex_order(i, j, l, d2_ji, d2_lj, d2_li, xs, ys)
11+
end
312

413
points = Data.points_to
514
combinations = Data.combinations_to
615
invariants = Data.invariants
716

8-
C, ℳ = _triangle_invariants(points)
9-
17+
# We always call _triangle_invariants with a Table, not a Vector{NamedTuple},
18+
# so we test that the output is correct for the Table input
19+
points_table = Table(points)
20+
C, ℳ = _triangle_invariants(points_table)
1021
@test length(C) == 1
11-
@test collect(C) == combinations
22+
# _canonical_vertex_order is called within _triangle_invariants, so we test that the
23+
# output C is consistent with the expected canonical vertex order for the given points and combinations
24+
@test C == [_canonical_vertex_order_test(c..., points_table.xcenter, points_table.ycenter) for c in combinations]
1225
@test== invariants
1326
end
1427

0 commit comments

Comments
 (0)