77Usage:
88 python diagnose_junctions.py --tif test.tif --out /tmp/junc
99"""
10- import argparse , os
10+ import argparse
11+ import os
1112import numpy as np
1213from PIL import Image
1314from skimage import morphology
@@ -26,9 +27,12 @@ def load_tif_2d(path):
2627 except ImportError :
2728 arr = np .array (Image .open (path ).convert ('L' ))
2829 if arr .ndim == 3 :
29- if arr .shape [0 ] <= 4 : arr = arr [0 ]
30- elif arr .shape [2 ] <= 4 : arr = (0.299 * arr [:,:,0 ]+ 0.587 * arr [:,:,1 ]+ 0.114 * arr [:,:,2 ]).astype (arr .dtype )
31- else : arr = arr .max (axis = 0 )
30+ if arr .shape [0 ] <= 4 :
31+ arr = arr [0 ]
32+ elif arr .shape [2 ] <= 4 :
33+ arr = (0.299 * arr [:,:,0 ]+ 0.587 * arr [:,:,1 ]+ 0.114 * arr [:,:,2 ]).astype (arr .dtype )
34+ else :
35+ arr = arr .max (axis = 0 )
3236 return arr
3337
3438
@@ -38,20 +42,23 @@ def skeleton_stats(binary, dilation_r=0):
3842 if dilation_r > 0 :
3943 mask = morphology .dilation (mask , morphology .disk (dilation_r ))
4044 skel = morphology .skeletonize (mask )
41- kernel = np .ones ((3 ,3 ), dtype = np .uint8 ); kernel [1 ,1 ] = 0
45+ kernel = np .ones ((3 ,3 ), dtype = np .uint8 )
46+ kernel [1 ,1 ] = 0
4247 nc = convolve (skel .astype (np .uint8 ), kernel , mode = 'constant' , cval = 0 )
4348 branch_mask = skel & (nc >= 3 )
4449 by , bx = np .where (branch_mask )
4550 return skel , branch_mask , nc , np .column_stack ((bx .astype (float ), by .astype (float ))) if len (by ) else np .zeros ((0 ,2 ))
4651
4752
4853def cluster (coords , r = 4.0 ):
49- if len (coords ) == 0 : return coords
54+ if len (coords ) == 0 :
55+ return coords
5056 tree = cKDTree (coords )
5157 visited = np .zeros (len (coords ), dtype = bool )
5258 out = []
5359 for i in range (len (coords )):
54- if visited [i ]: continue
60+ if visited [i ]:
61+ continue
5562 nb = tree .query_ball_point (coords [i ], r )
5663 out .append (coords [nb ].mean (axis = 0 ))
5764 visited [nb ] = True
@@ -93,7 +100,8 @@ def main():
93100 '(red=branch, green=edge, blue=endpoint, yellow=clustered vertex)' , fontsize = 11 )
94101 plt .tight_layout ()
95102 p = os .path .join (args .out , 'dilation_comparison.png' )
96- plt .savefig (p , dpi = 130 , bbox_inches = 'tight' ); plt .close ()
103+ plt .savefig (p , dpi = 130 , bbox_inches = 'tight' )
104+ plt .close ()
97105 print (f"Saved: { p } " )
98106
99107 # ── Figure 2: zoom into missed junctions (low-branch regions) ────────
@@ -118,7 +126,7 @@ def main():
118126 top_idx = np .argsort (scores )[::- 1 ][:9 ]
119127 top_junctions = ep_clustered [top_idx ]
120128
121- print (f "\n Top missed junction candidates (by nearby endpoint count):" )
129+ print ("\n Top missed junction candidates (by nearby endpoint count):" )
122130 for i , (cx , cy ) in enumerate (top_junctions ):
123131 print (f" [{ i } ] center=({ cx :.0f} ,{ cy :.0f} ) nearby_endpoints={ scores [top_idx [i ]]} " )
124132
@@ -139,21 +147,26 @@ def main():
139147
140148 # Overlay dilation=0 skeleton (green)
141149 crop0 = skel0 [y0 :y1 , x0 :x1 ]
142- ys0 , xs0 = np .where (crop0 ); ax .scatter (xs0 , ys0 , c = 'lime' , s = 2 , alpha = 0.7 )
150+ ys0 , xs0 = np .where (crop0 )
151+ ax .scatter (xs0 , ys0 , c = 'lime' , s = 2 , alpha = 0.7 )
143152
144153 # Overlay dilation=1 skeleton (magenta) for comparison
145154 crop1 = skel1 [y0 :y1 , x0 :x1 ]
146- ys1 , xs1 = np .where (crop1 ); ax .scatter (xs1 , ys1 , c = 'magenta' , s = 2 , alpha = 0.5 )
155+ ys1 , xs1 = np .where (crop1 )
156+ ax .scatter (xs1 , ys1 , c = 'magenta' , s = 2 , alpha = 0.5 )
147157
148158 # Branch pixels (dilation=0): red; dilation=1: yellow
149- br0c = bm0 [y0 :y1 , x0 :x1 ]; yb0 ,xb0 = np .where (br0c )
150- br1c = bm1 [y0 :y1 , x0 :x1 ]; yb1 ,xb1 = np .where (br1c )
159+ br0c = bm0 [y0 :y1 , x0 :x1 ]
160+ yb0 ,xb0 = np .where (br0c )
161+ br1c = bm1 [y0 :y1 , x0 :x1 ]
162+ yb1 ,xb1 = np .where (br1c )
151163 ax .scatter (xb0 , yb0 , c = 'red' , s = 20 , zorder = 10 )
152164 ax .scatter (xb1 , yb1 , c = 'yellow' , s = 20 , zorder = 10 , marker = '*' )
153165
154166 ax .set_title (f'Junction { i } ({ cx } ,{ cy } )\n '
155167 f'red=branch(d=0) yellow★=branch(d=1)' , fontsize = 8 )
156- ax .set_xlim (0 , x1 - x0 ); ax .set_ylim (y1 - y0 , 0 )
168+ ax .set_xlim (0 , x1 - x0 )
169+ ax .set_ylim (y1 - y0 , 0 )
157170 ax .axis ('off' )
158171
159172 for j in range (n_show , len (axes )):
@@ -168,7 +181,8 @@ def main():
168181 'yellow stars = branches recovered by disk(1) dilation' , fontsize = 11 )
169182 plt .tight_layout ()
170183 p2 = os .path .join (args .out , 'missed_junctions_zoom.png' )
171- plt .savefig (p2 , dpi = 130 , bbox_inches = 'tight' ); plt .close ()
184+ plt .savefig (p2 , dpi = 130 , bbox_inches = 'tight' )
185+ plt .close ()
172186 print (f"Saved: { p2 } " )
173187
174188 # ── Print recommendation ──────────────────────────────────────────────
@@ -179,7 +193,7 @@ def main():
179193 print (f"\n Vertex counts by dilation: { dict (zip (dils , counts ))} " )
180194 best = dils [int (np .argmax (counts ))]
181195 print (f"→ Recommended dilation: disk({ best } ) (maximizes vertex count)" )
182- print (f "\n In topology.py _labels_to_boundary or extract_topology, add:" )
196+ print ("\n In topology.py _labels_to_boundary or extract_topology, add:" )
183197 print (f" boundary = morphology.dilation(boundary, morphology.disk({ best } ))" )
184198
185199
0 commit comments