|
| 1 | +#set page( |
| 2 | + width: 28.8cm, |
| 3 | + height: 16.2cm, |
| 4 | + margin: (top: 1.4cm, bottom: 1cm, left: 1.8cm, right: 1.8cm), |
| 5 | + fill: white, |
| 6 | + numbering: "1", |
| 7 | + number-align: right, |
| 8 | +) |
| 9 | + |
| 10 | +#set text(font: "New Computer Modern", size: 16pt) |
| 11 | +#set par(justify: true, leading: 0.65em) |
| 12 | + |
| 13 | +#let accent = rgb("#1d4ed8") |
| 14 | +#let accent-light = rgb("#eff6ff") |
| 15 | +#let accent-mid = rgb("#93c5fd") |
| 16 | +#let muted = rgb("#6b7280") |
| 17 | +#let dark = rgb("#111827") |
| 18 | +#let green-bg = rgb("#f0fdf4") |
| 19 | +#let green-border = rgb("#bbf7d0") |
| 20 | +#let green-text = rgb("#166534") |
| 21 | + |
| 22 | +#let slide(title, body) = { |
| 23 | + pagebreak() |
| 24 | + block( |
| 25 | + width: 100%, |
| 26 | + inset: (top: 8pt, bottom: 10pt, left: 0pt, right: 0pt), |
| 27 | + stroke: (bottom: 2.5pt + accent), |
| 28 | + )[ |
| 29 | + #text(size: 24pt, fill: accent, weight: "bold")[#title] |
| 30 | + ] |
| 31 | + v(0.5cm) |
| 32 | + body |
| 33 | +} |
| 34 | + |
| 35 | +#let highlight(body) = { |
| 36 | + block( |
| 37 | + width: 100%, |
| 38 | + fill: accent-light, |
| 39 | + inset: 14pt, |
| 40 | + radius: 6pt, |
| 41 | + stroke: 1pt + accent-mid, |
| 42 | + )[ |
| 43 | + #body |
| 44 | + ] |
| 45 | +} |
| 46 | + |
| 47 | +#let proofbox(body) = { |
| 48 | + block( |
| 49 | + width: 100%, |
| 50 | + fill: green-bg, |
| 51 | + inset: 14pt, |
| 52 | + radius: 6pt, |
| 53 | + stroke: 1pt + green-border, |
| 54 | + )[ |
| 55 | + #body |
| 56 | + ] |
| 57 | +} |
| 58 | + |
| 59 | +#align(center + horizon)[ |
| 60 | + #block(width: 80%)[ |
| 61 | + #align(center)[ |
| 62 | + #text(size: 34pt, fill: accent, weight: "bold")[ |
| 63 | + VesSkel\ |
| 64 | + ] |
| 65 | + #v(1.0cm) |
| 66 | + #text(size: 20pt, fill: dark)[ |
| 67 | + Vessel Skeletonization and Graph-Based\ |
| 68 | + Phenotype Analysis in Retinal Fundus Images |
| 69 | + ] |
| 70 | + #v(1.2cm) |
| 71 | + #line(length: 40%, stroke: 1.5pt + accent) |
| 72 | + #v(0.8cm) |
| 73 | + #text(size: 16pt, fill: muted)[ |
| 74 | + Simon Wittmann |
| 75 | + |
| 76 | + Supervisor: Anna Möller |
| 77 | + ] |
| 78 | + #v(0.3cm) |
| 79 | + #text(size: 14pt, fill: muted)[ |
| 80 | + 21. April 2026 |
| 81 | + ] |
| 82 | + ] |
| 83 | + ] |
| 84 | +] |
| 85 | + |
| 86 | +#slide("Project Overview")[ |
| 87 | + #columns(2, gutter: 1.5em)[ |
| 88 | + *Objectives* |
| 89 | + + Efficient numba-parallelized Lee94 skeletonization (2D + 3D) |
| 90 | + + Application to HRF retinal vessel masks |
| 91 | + + Graph-based feature extraction |
| 92 | + + Phenotype differentiation |
| 93 | + |
| 94 | + #colbreak() |
| 95 | + |
| 96 | + *Key Deliverables* |
| 97 | + + `vesskel` Python package |
| 98 | + + Optimized 2D & 3D thinning |
| 99 | + + (kinda slow) Graph construction + basic feature extraction |
| 100 | + + Napari plugin for interactive use |
| 101 | + + Regression test suite (45 HRF + 3D brain) |
| 102 | + ] |
| 103 | +] |
| 104 | + |
| 105 | +#slide("Implementation Progress")[ |
| 106 | + #v(-0.3cm) |
| 107 | + #table( |
| 108 | + columns: (auto, 1fr), |
| 109 | + stroke: 0.5pt + rgb("#e5e7eb"), |
| 110 | + inset: 7pt, |
| 111 | + fill: (x, y) => if y == 0 { accent-light } else if calc.odd(y) { rgb("#f9fafb") }, |
| 112 | + [*Date*], [*Milestone*], |
| 113 | + [Apr 01], [Project setup, HRF dataset integration], |
| 114 | + [Apr 03], [Initial 3D Lee94 adapted for 2D images], |
| 115 | + [Apr 03], [First regression test (skeletonization)], |
| 116 | + [Apr 04], [Proof: Euler invariant is redundant in 2D], |
| 117 | + [Apr 04], [Parallelized 2D implementation (without Euler)], |
| 118 | + [Apr 09], [Numba parallelization for 2D thinning], |
| 119 | + [Apr 14], [Graph construction + feature extraction], |
| 120 | + [Apr 16], [Feature regression tests], |
| 121 | + [Apr 18], [3D implementation (skimage port)], |
| 122 | + [Apr 18], [3D Parallelization (candidate marking + adj LUTs)], |
| 123 | + [Apr 18], [Napari Plugin], |
| 124 | + [Apr 18], [CI and PyPi Release], |
| 125 | + ) |
| 126 | +] |
| 127 | + |
| 128 | +#slide("Phase 1: 3D -> 2D Adaptation")[ |
| 129 | + The initial implementation embedded 2D images into a 3D volume and ran the full Lee94 algorithm unchanged: |
| 130 | + |
| 131 | + $ V_(p,r,c) = cases(I_(r-1,c-1) "if" p = 1, 0 "otherwise") $ |
| 132 | + |
| 133 | + #v(0.3cm) |
| 134 | + |
| 135 | + Observations: |
| 136 | + + The 26-neighborhood collapses to an effective 8-neighborhood in the $p=1$ plane |
| 137 | + + All neighbors in the $p=0$ and $p=2$ planes are identically zero |
| 138 | + + Every foreground pixel is automatically a "border point" since the $p$-direction neighbors are always background. The 2 extra border check directions provides no constraining power in 2D. |
| 139 | + + Topology is preserved -- but at the cost of carrying redundant 3D machinery |
| 140 | +] |
| 141 | + |
| 142 | +#slide("Key Insight: Euler Check is Redundant in 2D")[ |
| 143 | + Lee94 preserves topology through four conditions: |
| 144 | + |
| 145 | + #v(0.3cm) |
| 146 | + #table( |
| 147 | + columns: (auto, 1fr, auto), |
| 148 | + stroke: 0.5pt + rgb("#e5e7eb"), |
| 149 | + inset: 7pt, |
| 150 | + align: (center, left, center), |
| 151 | + fill: (x, y) => if y == 0 { accent-light }, |
| 152 | + [*Nr.*], [*Condition*], [*Role*], |
| 153 | + [1], [Border point], [Restricts deletion to boundary], |
| 154 | + [2], [Endpoint], [Preserves line endings], |
| 155 | + [3], [Euler invariant ($delta chi = 0$)], [Prevents hole/tunnel changes], |
| 156 | + [4], [Simple point ($O = 1$)], [Prevents disconnection], |
| 157 | + ) |
| 158 | + |
| 159 | + #v(0.4cm) |
| 160 | + In 3D, conditions 3 and 4 are *complementary*:\ |
| 161 | + Euler detects holes/tunnels, Simple detects disconnection. |
| 162 | + |
| 163 | + #v(0.3cm) |
| 164 | + #highlight[ |
| 165 | + *In 2D, the Euler check is provably redundant.*\ |
| 166 | + The Simple-Point-Check alone suffices for topology preservation. |
| 167 | + ] |
| 168 | +] |
| 169 | + |
| 170 | +#slide("Proof: Impossibility of Hole Creation")[ |
| 171 | + #proofbox[ |
| 172 | + *Lemma 1:* _In 2D Lee94 thinning, no border point removal can create a hole ($delta H > 0$)._ |
| 173 | + ] |
| 174 | + |
| 175 | + #v(0.3cm) |
| 176 | + *Proof.* For a hole to be created by removing pixel $x$, all 4-neighbors must be foreground: |
| 177 | + |
| 178 | + $ forall y in N_4(x) : V_y = 1 quad arrow quad x "fails border point condition" $ |
| 179 | + |
| 180 | + If any 4-neighbor is background, it provides a 4-connected path from the newly-background $x$ to the exterior, so no new isolated background component (hole) is formed: |
| 181 | + |
| 182 | + $ exists y in N_4(x) : V_y = 0 quad arrow quad "no hole created" $ |
| 183 | + |
| 184 | + Since border point candidates must have at least one background 4-neighbor, hole creation is impossible. $square$ |
| 185 | +] |
| 186 | + |
| 187 | +#slide("Proof: Impossibility of Hole Elimination & Theorem")[ |
| 188 | + #proofbox[ |
| 189 | + *Lemma 2:* _No simple point removal can eliminate a hole ($delta H < 0$)._ |
| 190 | + ] |
| 191 | + |
| 192 | + #v(0.2cm) |
| 193 | + If removing $x$ would eliminate a hole, foreground pixels in $N(x)$ must form $>= 2$ disconnected components wrapping around separate background regions. But the Simple-Point-Check requires $O(S inter N(x)) = 1$. Contradiction. $square$ |
| 194 | + |
| 195 | + #v(0.4cm) |
| 196 | + |
| 197 | + #block( |
| 198 | + width: 100%, |
| 199 | + fill: rgb("#fef3c7"), |
| 200 | + inset: 14pt, |
| 201 | + radius: 6pt, |
| 202 | + stroke: 1pt + rgb("#fde68a"), |
| 203 | + )[ |
| 204 | + *Theorem:* _For 2D Lee94 thinning, the Euler invariant check is redundant._ |
| 205 | + |
| 206 | + #v(0.2cm) |
| 207 | + Since $delta H > 0$ is impossible (Lemma 1) and $delta H < 0$ is impossible (Lemma 2): |
| 208 | + $ delta H = 0 quad arrow quad delta chi = delta O - delta H = delta O $ |
| 209 | + The Simple-Point-Check ensures $delta O = 0$, so $delta chi = 0$ is automatically satisfied. $square$ |
| 210 | + ] |
| 211 | +] |
| 212 | + |
| 213 | +#slide("Practical Implications of the Proof")[ |
| 214 | + The proof enables significant simplifications in the 2D implementation: |
| 215 | + |
| 216 | + #v(0.3cm) |
| 217 | + + *No 3D embedding required* -- operate directly on padded 2D image $I in {0,1}^((H+2) times (W+2))$ |
| 218 | + + *Border detection simplifies* from 6-neighbor 3D checks to 4-neighbor 2D checks (N, S, E, W) |
| 219 | + + *Euler-Invariant-Check omitted entirely* -- no octree construction, no Euler LUT |
| 220 | + + *Simple-Point-Check via direct flood fill* -- count 8-connected foreground components with DFS (no octrees) |
| 221 | + |
| 222 | + #v(0.4cm) |
| 223 | + The thinning loop iterates over 4 directional sub-iterations, each: |
| 224 | + 1. Collecting border candidates along one direction in parallel |
| 225 | + 2. Sequential recheck for simple-point preservation |
| 226 | +] |
| 227 | + |
| 228 | +#slide("Phase 2: Native 2D Implementation")[ |
| 229 | + Refactored from 3D volume to direct 2D processing (`thin_2d.py`): |
| 230 | + |
| 231 | + #v(0.3cm) |
| 232 | + #table( |
| 233 | + columns: (1fr, 1fr), |
| 234 | + stroke: 0.5pt + rgb("#e5e7eb"), |
| 235 | + inset: 8pt, |
| 236 | + fill: (x, y) => if y == 0 { accent-light }, |
| 237 | + [*3D Embedded (before)*], [*2D Native (after)*], |
| 238 | + [3D volume $3 times (H+2) times (W+2)$], [Padded 2D $(H+2) times (W+2)$], |
| 239 | + [6 border directions], [4 border directions (N, S, E, W)], |
| 240 | + [26-neighborhood lookup], [8-neighborhood lookup], |
| 241 | + [Euler check + Simple check], [Simple check *only*], |
| 242 | + [Octree-based connectivity], [Direct flood-fill DFS], |
| 243 | + ) |
| 244 | + |
| 245 | + #v(0.4cm) |
| 246 | + #highlight[ |
| 247 | + *Performance:* Thinning all 45 HRF Images takes $approx 9 "seconds" (approx 0.197 s "per image")$. \ |
| 248 | + ] |
| 249 | +] |
| 250 | + |
| 251 | +#slide("Phase 3: Numba Parallelization")[ |
| 252 | + The thinning loop has an inherent data dependency -- deletions in one iteration affect the next. But *candidate marking* is data-parallel: |
| 253 | + |
| 254 | + #v(0.3cm) |
| 255 | + |
| 256 | + #highlight[ |
| 257 | + *Strategy: parallel candidate marking, sequential deletion* |
| 258 | + ] |
| 259 | + |
| 260 | + #v(0.3cm) |
| 261 | + + *Parallel phase* (`prange` over rows): each row independently marks candidates that pass border + endpoint + simple-point checks |
| 262 | + + *Merge phase*: collect per-row results into a flat candidate array |
| 263 | + + *Sequential phase*: iterate candidates, recheck simple-point condition, delete if still valid |
| 264 | +] |
| 265 | + |
| 266 | +#slide("Phase 4: Full 3D Implementation")[ |
| 267 | + Ported the scikit-image Cython implementation to pure Python + Numba (`thin_3d.py`): |
| 268 | + |
| 269 | + #v(0.3cm) |
| 270 | + + Faithful port of `skimage.morphology._skeletonize_lee_cy` |
| 271 | + + Euler LUT + octant index table for Euler invariant computation |
| 272 | + + 26-neighborhood adjacency via pre-computed lookup tables |
| 273 | + + All checks active: border, endpoint, Euler invariant, simple point |
| 274 | + |
| 275 | + #v(0.3cm) |
| 276 | + #highlight[ |
| 277 | + *Validated:* The output is *bit-identical* to `skimage.morphology.skeletonize` on the |
| 278 | + scikit-image brain volume test case. |
| 279 | + ] |
| 280 | + |
| 281 | + Further Refinement: `_mark_removable_candidates()` marks deletable pixels in parallel; sequential `_apply_removals()` rechecks and deletes |
| 282 | +] |
| 283 | + |
| 284 | +#slide("Phase 5: Graph Construction & Feature Extraction")[ |
| 285 | + Skeletons are transformed into graph representations via `skan.Skeleton` |
| 286 | + |
| 287 | + #columns(2, gutter: 1.5em)[ |
| 288 | + *Topology* |
| 289 | + + Number of nodes |
| 290 | + + Number of edges |
| 291 | + + Number of endpoints |
| 292 | + + Number of bifurcations |
| 293 | + + Number of connected components |
| 294 | + + Mean node degree |
| 295 | + + Maximum node degree |
| 296 | + |
| 297 | + #colbreak() |
| 298 | + |
| 299 | + *Geometry* |
| 300 | + + Total vessel length |
| 301 | + + Mean segment length |
| 302 | + + Std segment length |
| 303 | + + Max segment length |
| 304 | + + Min segment length |
| 305 | + + Mean tortuosity |
| 306 | + + Std tortuosity |
| 307 | + ] |
| 308 | + |
| 309 | + #v(0.4cm) |
| 310 | + *Example:* Sample `01_h` from HRF: |
| 311 | + #table( |
| 312 | + columns: (auto, auto, auto, auto, auto), |
| 313 | + stroke: 0.5pt + rgb("#e5e7eb"), |
| 314 | + inset: 6pt, |
| 315 | + fill: (x, y) => if y == 0 { accent-light }, |
| 316 | + [*Nodes*], [*Edges*], [*Bifurcations*], [*Mean tortuosity*], [*Components*], |
| 317 | + [720], [765], [404], [1.086], [14], |
| 318 | + ) |
| 319 | +] |
| 320 | + |
| 321 | +#slide("Napari Plugin: Interactive Skeletonization")[ |
| 322 | + `vesskel` now also ships a Napari plugin for visual exploration: |
| 323 | + |
| 324 | + #v(0.3cm) |
| 325 | + + Registered via `napari.yaml` manifest + `napari.manifest` entry point |
| 326 | + + *Lee94 Thinning* widget: select an image layer, run thinning, see the skeleton as a new labels layer |
| 327 | + |
| 328 | + |
| 329 | + #v(0.3cm) |
| 330 | + #highlight[ |
| 331 | + The Dispatcher (`thin.py`) auto-selects 2D or 3D based on `img.ndim`, so it works for both. |
| 332 | + ] |
| 333 | +] |
| 334 | + |
| 335 | +#slide("Some Testing")[ |
| 336 | + Three-tier regression test suite with baseline comparison: |
| 337 | + |
| 338 | + #v(0.3cm) |
| 339 | + #table( |
| 340 | + columns: (auto, 1fr, auto), |
| 341 | + stroke: 0.5pt + rgb("#e5e7eb"), |
| 342 | + inset: 7pt, |
| 343 | + fill: (x, y) => if y == 0 { accent-light }, |
| 344 | + [*Test*], [*What it checks*], [*Scope*], |
| 345 | + [2D Regression], [Thinning + features on HRF vs. saved baselines], [45 samples], |
| 346 | + [3D Regression], [Thinning + features on brain volume vs. saved baselines], [1 volume], |
| 347 | + [3D Comparison], [vesskel vs. `skimage.morphology.skeletonize` -- bit-identical output], [1 volume], |
| 348 | + ) |
| 349 | +] |
| 350 | + |
| 351 | +#slide("Package Structure")[ |
| 352 | + ```text |
| 353 | + vesskel |
| 354 | + ├── __init__.py # public API (currently only exposes thin) |
| 355 | + ├── _napari.py # napari widget |
| 356 | + ├── benchmark # some benchmark files |
| 357 | + ├── features.py # graph construction & feature extraction |
| 358 | + ├── hrf.py # dataset loader |
| 359 | + ├── napari.yaml # plugin manifest |
| 360 | + ├── thin.py # dispatcher with lazy imports based on ndim |
| 361 | + ├── thin_2d.py # 2D Lee |
| 362 | + └── thin_3d.py # 3D Lee |
| 363 | + ``` |
| 364 | +] |
| 365 | + |
| 366 | +#slide("Summary & Next Steps")[ |
| 367 | + |
| 368 | + #columns(2, gutter: 1.5em)[ |
| 369 | + *What I did:* |
| 370 | + + Native 2D thinning without Euler invariant check (with proof) |
| 371 | + + Full 3D thinning -- bit-identical to scikit-image |
| 372 | + + Numba parallelization for both 2D and 3D |
| 373 | + + Rudimentary graph-based feature extraction |
| 374 | + + Rudimentary Napari plugin |
| 375 | + + Regression test suite |
| 376 | + + CI/CD with PyPi releases |
| 377 | + |
| 378 | + #colbreak() |
| 379 | + |
| 380 | + *Next steps:* |
| 381 | + + More Graph Features (e.g. fractal dimension) |
| 382 | + + DIY Graph Assembly? |
| 383 | + + Statistical analysis of phenotype differentiation |
| 384 | + + 3D Benchmarking with huge volumes: `vesskel` vs. `vesselvio` |
| 385 | + + Parallelize rechecking phase? (maybe too much overhead) |
| 386 | + + Toggleable Feature extraction in Napari |
| 387 | + + Configurable/Savable/Loadable Batch jobs in Napari |
| 388 | + ] |
| 389 | +] |
| 390 | + |
0 commit comments