2525};
2626"""
2727
28- _3d_image_1nn_str = r"""
28+
29+ def _3d_image_1nn_str_func (masked : bool ) -> str :
30+ return """
2931 ull z = i / (height * width);
3032 ull y = (i % (height * width)) / width;
3133 ull x = i % width;
3436
3537 #pragma unroll
3638 for (int j = 0; j < d_size; ++j)
37- {
39+ {{
3840 ull nz = z + dz[j];
3941 ull ny = y + dy[j];
4042 ull nx = x + dx[j];
4143
4244 // zero check is not needed since we are using unsigned integers
4345 if (nz < depth && ny < height && nx < width)
44- {
46+ {{
4547 ull nidx = IDX(nz, ny, nx, height, width);
46- if (mask[nidx]) {
48+ {if_masked} { {
4749 ull nval = VAL_IDX(image[nidx], nidx);
48- if (nval < min_value) {
50+ if (nval < min_value) {{
4951 min_value = nval;
50- }
51- }
52- }
53- }
52+ }}
53+ }}
54+ }}
55+ }}
5456
5557 roots[i] = min_value;
56- """
58+ """ .format (if_masked = "if (mask[nidx])" if masked else "" )
59+
5760
5861_3d_image_1nn = cp .ElementwiseKernel (
59- r"raw uint16 image, raw bool mask, int64 depth, int64 height, int64 width" ,
62+ r"raw uint16 image, int64 depth, int64 height, int64 width" ,
6063 r"raw uint64 roots" ,
61- f"if (mask[i]) {{ \n { _3d_image_1nn_str } \n }} \n " ,
64+ _3d_image_1nn_str_func ( masked = False ) ,
6265 r"_3d_image_1nn" ,
6366 preamble = preamble ,
6467)
6568
69+ _3d_image_1nn_masked = cp .ElementwiseKernel (
70+ r"raw uint16 image, raw bool mask, int64 depth, int64 height, int64 width" ,
71+ r"raw uint64 roots" ,
72+ f"if (mask[i]) {{\n { _3d_image_1nn_str_func (masked = True )} \n }}\n " ,
73+ r"_3d_image_1nn_masked" ,
74+ preamble = preamble ,
75+ )
76+
6677_3d_image_1nn_sparse = cp .ElementwiseKernel (
6778 r"raw int64 indices, raw uint16 image, raw bool mask, int64 depth, int64 height, int64 width" ,
6879 r"raw uint64 roots" ,
69- f"i = indices[i];\n { _3d_image_1nn_str } \n " ,
80+ f"i = indices[i];\n { _3d_image_1nn_str_func ( masked = True ) } \n " ,
7081 r"_3d_image_1nn_sparse" ,
7182 preamble = preamble ,
7283)
8091"""
8192
8293_assign_root = cp .ElementwiseKernel (
94+ r"" ,
95+ r"raw uint64 roots" ,
96+ _assign_root_str ,
97+ r"_assign_root" ,
98+ preamble = preamble ,
99+ )
100+
101+ _assign_root_masked = cp .ElementwiseKernel (
83102 r"raw bool mask" ,
84103 r"raw uint64 roots" ,
85104 f"if (mask[i]) {{\n { _assign_root_str } \n }}\n " ,
86- r"_assign_root " ,
105+ r"_assign_root_masked " ,
87106 preamble = preamble ,
88107)
89108
95114 preamble = preamble ,
96115)
97116
98- _merge_flat_zones_str = r"""
117+
118+ def _merge_flat_zones_str_func (masked : bool ) -> str :
119+ return """
99120 ull r = roots[i];
100121 ull root_idx = r & IDX_MASK;
101122 ull i_value = image[i];
110131
111132 #pragma unroll
112133 for (int j = 0; j < d_size; ++j)
113- {
134+ {{
114135 ull nz = z + dz[j];
115136 ull ny = y + dy[j];
116137 ull nx = x + dx[j];
117138
118139 if (nz < depth && ny < height && nx < width)
119- {
140+ {{
120141 ull nidx = IDX(nz, ny, nx, height, width);
121- if (mask[nidx]) {
142+ {if_masked} { {
122143 ull nr = roots[nidx];
123- if (nr < r) { // deeper root
144+ if (nr < r) {{ // deeper root
124145 r = nr;
125- }
126- }
127- }
128- }
129- if (r != roots[i]) {
130- if (root_idx != i) {
146+ }}
147+ }}
148+ }}
149+ }}
150+ if (r != roots[i]) {{
151+ if (root_idx != i) {{
131152 roots[i] = r; // using atomicMin only when necessary
132- }
153+ }}
133154 atomicMin(&roots[root_idx], r);
134155 atomicAdd(&changed[0], 1);
135- }
136- """
156+ }}
157+ """ .format (if_masked = "if (mask[nidx])" if masked else "" )
158+
137159
138160_merge_flat_zones = cp .ElementwiseKernel (
139- r"raw uint16 image, raw bool mask, int64 h, int64 depth, int64 height, int64 width" ,
161+ r"raw uint16 image, int64 h, int64 depth, int64 height, int64 width" ,
140162 r"raw uint64 roots, raw uint64 changed" ,
141- f"if (mask[i]) {{ \n { _merge_flat_zones_str } \n }} \n " ,
163+ _merge_flat_zones_str_func ( masked = False ) ,
142164 r"_merge_flat_zones" ,
143165 preamble = preamble ,
144166)
145167
168+ _merge_flat_zones_masked = cp .ElementwiseKernel (
169+ r"raw uint16 image, raw bool mask, int64 h, int64 depth, int64 height, int64 width" ,
170+ r"raw uint64 roots, raw uint64 changed" ,
171+ f"if (mask[i]) {{\n { _merge_flat_zones_str_func (masked = True )} \n }}\n " ,
172+ r"_merge_flat_zones_masked" ,
173+ preamble = preamble ,
174+ )
175+
146176_merge_flat_zones_sparse = cp .ElementwiseKernel (
147177 r"raw int64 indices, raw uint16 image, raw bool mask, int64 h, int64 depth, int64 height, int64 width" ,
148178 r"raw uint64 roots, raw uint64 changed" ,
149- f"i = indices[i];\n { _merge_flat_zones_str } \n " ,
179+ f"i = indices[i];\n { _merge_flat_zones_str_func ( masked = True ) } \n " ,
150180 r"_merge_flat_zones_sparse" ,
151181 preamble = preamble ,
152182)
164194
165195def watershed_from_minima (
166196 image : cp .ndarray ,
167- mask : cp .ndarray ,
197+ mask : cp .ndarray | None = None ,
168198 h : int = 0 ,
169199 sparse : bool = True ,
170200) -> cp .ndarray :
@@ -200,9 +230,13 @@ def watershed_from_minima(
200230
201231 if image .ndim == 2 :
202232 image = image [None , ...]
203- mask = mask [None , ...]
233+ if mask is not None :
234+ mask = mask [None , ...]
235+
236+ if sparse and mask is None :
237+ raise ValueError ("'mask' is required for 'sparse' mode" )
204238
205- if image .shape != mask .shape :
239+ if mask is not None and image .shape != mask .shape :
206240 raise ValueError (f"Image and mask must have the same shape: { image .shape } != { mask .shape } " )
207241
208242 if image .dtype != np .uint16 :
@@ -214,8 +248,9 @@ def watershed_from_minima(
214248
215249 roots = cp .arange (size , dtype = cp .uint64 )
216250
217- flat_mask = mask .ravel ()
218251 flat_image = image .ravel ()
252+ if mask is not None :
253+ flat_mask = mask .ravel ()
219254
220255 # TODO:
221256 # - reduce branching of kernels by having different branch for border and non-border pixels
@@ -256,33 +291,60 @@ def watershed_from_minima(
256291 n_iters += 1
257292
258293 else :
259- _3d_image_1nn (
260- flat_image ,
261- flat_mask ,
262- int (image .shape [0 ]),
263- int (image .shape [1 ]),
264- int (image .shape [2 ]),
265- roots ,
266- size = size ,
267- )
268- while changed [0 ]:
269- _assign_root (flat_mask , roots , size = size )
270- changed [0 ] = 0
271- _merge_flat_zones (
294+ if mask is not None :
295+ _3d_image_1nn_masked (
272296 flat_image ,
273297 flat_mask ,
274- h ,
275298 int (image .shape [0 ]),
276299 int (image .shape [1 ]),
277300 int (image .shape [2 ]),
278301 roots ,
279- changed ,
280302 size = size ,
281303 )
304+ else :
305+ _3d_image_1nn (
306+ flat_image ,
307+ int (image .shape [0 ]),
308+ int (image .shape [1 ]),
309+ int (image .shape [2 ]),
310+ roots ,
311+ size = size ,
312+ )
313+ while changed [0 ]:
314+ if mask is not None :
315+ _assign_root_masked (flat_mask , roots , size = size )
316+ else :
317+ _assign_root (roots , size = size )
318+
319+ changed [0 ] = 0
320+ if mask is not None :
321+ _merge_flat_zones_masked (
322+ flat_image ,
323+ flat_mask ,
324+ h ,
325+ int (image .shape [0 ]),
326+ int (image .shape [1 ]),
327+ int (image .shape [2 ]),
328+ roots ,
329+ changed ,
330+ size = size ,
331+ )
332+ else :
333+ _merge_flat_zones (
334+ flat_image ,
335+ h ,
336+ int (image .shape [0 ]),
337+ int (image .shape [1 ]),
338+ int (image .shape [2 ]),
339+ roots ,
340+ changed ,
341+ size = size ,
342+ )
282343 n_iters += 1
283344
284345 LOG .info ("Performed %d minima merging iterations" , n_iters )
285346
286- _relabel_inplace (flat_mask , roots )
347+ if mask is not None :
348+ _relabel_inplace (flat_mask , roots )
287349
288350 return roots .reshape (orig_shape )
0 commit comments