Skip to content

Optimize triangles - #55

Merged
icweaver merged 19 commits into
mainfrom
optimize-triangles
May 21, 2026
Merged

Optimize triangles#55
icweaver merged 19 commits into
mainfrom
optimize-triangles

Conversation

@cgarling

@cgarling cgarling commented May 16, 2026

Copy link
Copy Markdown
Member

Optimizes the building of the correspondences and triangles to minimize the impact of that step in the case that it is repeated in workflows when fitting multiple images (#52). Also added a benchmark script for you to compare if you'd like. Note that comparisons to main are not totally fair as canonical_vertex_order is now called internally to _triangle_invariants whereas before it was an additional post-processing step (so there's actually more work being done in _triangle_invariants than before. Even so, _triangle_invariants for 100 stars now takes ~1.5 ms on my machine compared to 10.9 ms on main.

  • Replaced Combinatorics.combinations iterator with explicit triple-nested for loops over strictly increasing index triples (i < j < l), eliminating iterator allocation.
  • Replaced Distances.euclidean with inlined squared-distance calculation, avoiding sqrt until the final ratio computation (the sqrt may not even be strictly necessary, we could do matching in distance^2 space, but I digress)
  • Replaced sort!([...]) per triangle with a branchless 3-element sorting network
  • Replaced map + stack result construction with a single pre-allocated Matrix{Float64} written in-place, eliminating intermediate allocations and a second pass over the data
  • Integrated _canonical_vertex_order directly into _triangle_invariants, so C is returned in canonical vertex order (apex last, base vertices CCW) without a separate post-processing pass; the squared distances already computed in the loop are reused directly
  • Removed the Combinatorics, Distances dependencies

@codecov

codecov Bot commented May 16, 2026

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.
✅ Project coverage is 99.45%. Comparing base (8841f99) to head (b2e9551).

Additional details and impacted files
@@            Coverage Diff             @@
##             main      #55      +/-   ##
==========================================
+ Coverage   99.37%   99.45%   +0.08%     
==========================================
  Files           5        5              
  Lines         161      185      +24     
==========================================
+ Hits          160      184      +24     
  Misses          1        1              
Flag Coverage Δ
unittests 99.45% <100.00%> (+0.08%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

@cgarling
cgarling force-pushed the optimize-triangles branch from e289665 to 7f8ece8 Compare May 16, 2026 22:54
@cgarling
cgarling requested a review from icweaver May 17, 2026 00:59
@cgarling

Copy link
Copy Markdown
Member Author

Ready for review

Just a note, the _build_correspondences is still somewhat slow (~50 ms for 100 stars) because it is building all correspondences between all triangles. In particular making the triangles is cheap but the KDTree build and nn calls in _build_correspondences are more expensive the more triangles there are. astroalign's solution to this is to trim the triangle list --

To make triangles for one image,

  • walk down list of stars
  • only construct triangles for n nearest stars
  • at the end, deduplicate the triangles

This logic will be easier to implement with the explicit loop architecture I use in _triangle_invariants here.
Changing that logic will be significant speedup but I think we should merge this first before we do that.

@icweaver

Copy link
Copy Markdown
Member

This is brilliant, thanks Chris! I am also going to use this as an opportunity for me to finally learn about AirspeedVeclocity.jl and implement it in this repo. Once ConsensusFitting.jl is registered, I think the PR here would be a great first test case, so I will just leave it open until then.

btw, I am seeing similar speed-ups on my end too 👍🏾

┌────────────────────────┬─────────────┬────────────┬─────────┐
│ Benchmark              │ Median Time │     Memory │  Allocs │
├────────────────────────┼─────────────┼────────────┼─────────┤
│ _build_correspondences │  159.065 ms │ 143.37 MiB │ 2910662 │
│ _triangle_invariants   │   15.031 ms │  39.48 MiB │  808508 │
└────────────────────────┴─────────────┴────────────┴─────────┘

┌────────────────────────┬─────────────┬───────────┬────────┐
│ Benchmark              │ Median Time │    Memory │ Allocs │
├────────────────────────┼─────────────┼───────────┼────────┤
│ _build_correspondences │   87.688 ms │ 34.80 MiB │ 323439 │
│ _triangle_invariants   │    1.651 ms │  6.17 MiB │      6 │
└────────────────────────┴─────────────┴───────────┴────────┘

@cgarling

Copy link
Copy Markdown
Member Author

Sounds good, FYI I think you will need to add an extra-pkgs argument to the airspeedvelocity action, see https://github.com/cgarling/TRGBDistances.jl/pull/10/changes

@cgarling

Copy link
Copy Markdown
Member Author

You could even get into stuff like

@inline function _canonical_vertex_order(i, j, l, d2_ji, d2_lj, d2_li, xs, ys)
    # There are three possible longest edges:
    #   1. ji is longest  (c1 == true)
    #   2. lj is longest  (c2 == true)
    #   3. li is longest  (neither c1 nor c2)
    ji_longest = (d2_ji >= d2_lj) & (d2_ji >= d2_li)
    lj_longest = (d2_lj >= d2_ji) & (d2_lj >= d2_li)

    # The nested ifelse expressions are equivalent to:
    #
    #   if c1
    #       v1, v2, apex = i, j, l
    #   elseif c2
    #       v1, v2, apex = j, l, i
    #   else
    #       v1, v2, apex = i, l, j
    #   end
    # Using ifelse keeps this branchless
    v1 = ifelse(ji_longest, i, ifelse(lj_longest, j, i)) 
    v2 = ifelse(ji_longest, j, ifelse(lj_longest, l, l))
    apex = ifelse(ji_longest, l, ifelse(lj_longest, i, j))

    # Enforce CCW winding: swap base vertices if cross product is negative
    @inbounds cross = (xs[v2] - xs[v1]) * (ys[apex] - ys[v1]) - (ys[v2] - ys[v1]) * (xs[apex] - xs[v1])
    swap = cross < 0
    vv1 = ifelse(swap, v2, v1)
    vv2 = ifelse(swap, v1, v2)
    return vv1, vv2, apex
end

This emits better LLVM and results in speedup for _triangle_invariants for 100 stars from 1.4 ms -> 1.0 ms on my machine. I haven't pushed it because it's less readable IMO so there's some tradeoff between keeping the code readable and optimizing for performance. I can push it if you'd like, or not, either way is fine with me

Comment thread src/register.jl
Comment thread docs/src/walkthrough.md
Comment thread src/register.jl Outdated
@icweaver

Copy link
Copy Markdown
Member

Thanks Chris, just left my very light review responses above.

Sounds good, FYI I think you will need to add an extra-pkgs argument to the airspeedvelocity action, see https://github.com/cgarling/TRGBDistances.jl/pull/10/changes

Thanks for the heads up! I'm also getting up to speed on what gaps might still exist with [workspace] so that stuff like this will just work™

This emits better LLVM and results in speedup for _triangle_invariants for 100 stars from 1.4 ms -> 1.0 ms on my machine. I haven't pushed it because it's less readable IMO so there's some tradeoff between keeping the code readable and optimizing for performance. I can push it if you'd like, or not, either way is fine with me

That's a good concern. I think just having good docs and/or docstrings for things like this would probably be enough if you think it would be worth it for the performance improvement

@icweaver
icweaver merged commit 02a35cd into main May 21, 2026
13 of 14 checks passed
@icweaver
icweaver deleted the optimize-triangles branch May 21, 2026 05:04
icweaver referenced this pull request Jun 4, 2026
If you copy-paste the revised function definitions from this PR into
your REPL you can run the following comparison which will show the new
implementation is ~3x faster and allocation free **EDIT: If you have to
collect the input (as we do because of how Photometry.jl passes us the
cutout) then it will be 2 allocations and adds ~100 ns to the base case
speed**. I chose to do `com_psf(T::Type{<:AbstractFloat}, img_ap,
rel_thresh)` because Float64 is not really any slower on my testing so
this way it makes it easier for us to switch to Float64 in the future if
we every wanted to.

```julia
using PSFModels: gaussian
using BenchmarkTools
import Astroalign

const x = 1:20
const y = 1:20
const T = Float32
model(x, y, amp) = gaussian(T, 4, 4; x, y, fwhm=3) * amp
data = model.(x, y', 10) .+ T(0.1) * randn(T, length(x), length(y))
result_orig = Astroalign.com_psf(data; rel_thresh=0.1f0)
result_revised = com_psf(data; rel_thresh=0.1f0)
println("Original COM: ", (result_orig.psf_params.x, result_orig.psf_params.y))
println("Revised COM: ", (result_revised.psf_params.x, result_revised.psf_params.y))
println("Original FWHM: ", result_orig.psf_params.fwhm)
println("Revised FWHM: ", result_revised.psf_params.fwhm)
println("Original benchmark: ")
display(@benchmark Astroalign.com_psf($data; rel_thresh=$0.1f0))
println("Revised benchmark: ")
display(@benchmark com_psf($data; rel_thresh=$0.1f0))
```

My result:

**EDIT:** Because Photometry.jl will pass us an object from
Transducers.jl we have to call `collect`, giving us 2 allocations and
runtime +100 ns over not collecting (if we had a pure matrix input).

```julia
julia> println("Original COM: ", (result_orig.psf_params.x, result_orig.psf_params.y))
Original COM: (3.9974918f0, 3.9858212f0)

julia> println("Revised COM: ", (result_revised.psf_params.x, result_revised.psf_params.y))
Revised COM: (3.9974945f0, 3.9858336f0)

julia> println("Original FWHM: ", result_orig.psf_params.fwhm)
Original FWHM: (2.31581660320762, 2.3766457470426725)

julia> println("Revised FWHM: ", result_revised.psf_params.fwhm)
Revised FWHM: (2.3164616f0, 2.3772223f0)

julia> display(@benchmark Astroalign.com_psf($data; rel_thresh=$0.1f0))
BenchmarkTools.Trial: 10000 samples with 10 evaluations per sample.
 Range (min … max):  1.408 μs …  1.050 ms  ┊ GC (min … max):  0.00% … 99.19%
 Time  (median):     1.544 μs              ┊ GC (median):     0.00%
 Time  (mean ± σ):   2.372 μs ± 15.628 μs  ┊ GC (mean ± σ):  21.39% ±  4.62%

  ▆█▇▅▃▁▁▁                ▃▄▃▂                               ▂
  █████████▇▇▆▇▇▇▇▆▇▆▆▆▇▆██████▆▆▆▅▄▅▄▅▅▃▅▆███▇▅▄▄▄▁▁▄▁▆▅▇██ █
  1.41 μs      Histogram: log(frequency) by time     5.56 μs <

 Memory estimate: 10.41 KiB, allocs estimate: 20.

julia> println("Revised benchmark: ")
Revised benchmark:

julia> @benchmark com_psf($data; rel_thresh=$0.1f0)
BenchmarkTools.Trial: 10000 samples with 183 evaluations per sample.
 Range (min … max):  571.716 ns …   9.490 μs  ┊ GC (min … max): 0.00% … 88.05%
 Time  (median):     600.563 ns               ┊ GC (median):    0.00%
 Time  (mean ± σ):   678.928 ns ± 532.693 ns  ┊ GC (mean ± σ):  8.48% ±  9.61%

  █▅▂▁▁                                                         ▁
  ██████▇▇▆▆▄▅▁▄▃▁▃▁▁▁▁▁▃▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▃▁▃▄▃▄▄▄▄▄▄▄▃▄▅ █
  572 ns        Histogram: log(frequency) by time        4.5 μs <

 Memory estimate: 1.64 KiB, allocs estimate: 2.
```

---------

Co-authored-by: Ian Weaver <weaveric@gmail.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants