Troubles using watershed algorithm

[DiegoPerezDones]

Hello everyone,

I am having trouble using watershed function. I have followed the tutorials but still haven't been able to get an acceptable solution. I would be very greatful if anyone could give me some insights in what ccoul be happening here.

a = TiffImages.load("C://Users//i7//Desktop//imagenes_mario\\24_1.tif");

# convert to Gray to discard useless green and blue channels (all info is in red)
# and reshape since the color info is intercalated in the z dimension
b = reshape(Gray.(a), size(a)[[1,2]]..., 3, :)

#Select just the channels with cell nuclei
nuclei = b[:,:,3,:]
img = nuclei

#Obtain binary image based on a filters pipeline, it is not important how the pipeline actually works.
binary = processing_nuclei_multiThreading(nuclei,6)

#Invert binary image
binary_inv = binary .< 1

#For ease, we will aply our code just in one slice of the stack but bear in mind I do the same for each slice in a for loop
i = 17
dist = distance_transform(feature_transform(binary_inv[:,:,i]))
if sum(dist) != 0 #Here we asses if the slice if empty of objects.
        peaks = zeros(size(dist))
        peaks[findlocalmaxima(dist, window = (6,6))] .= 1
        markers = label_components(peaks)
        segments = watershed(dist, Int.(markers), mask = binary[:,:,i] .> 0, compactness = 1)
        labels[:,:,i] = p = labels_map(segments)
        IndirectArray(p .+1, distinguishable_colors(maximum(p)+1)) #This is just for print a coloured version of the labels.
    else
        labels[:,:,i] = img[:,:,i]
    end

So after all this pipeline, the final result of the segmentation is the next:

As you may see, the result is not terrible but could be improved. You should also take into account that, in this slice in particular, this is working quite fine in comparison with, for example the first slice. Any advices to improve the result?

Thank you all in advance! @tlnagy

[tlnagy]

@DiegoPerezDones Do you mind tweaking this to be a minimal working example? I think at the minimum we would need access to the TIF and I suspect processing_nuclei_multiThreading is a custom function of yours that we don't have access to. So I would just upload the finished binary mask output of that function.

Now looking at the output, my best guess is that the seeding of the watershed algorithm is the problem, but without being able to run your code it's hard to say for sure.

[DiegoPerezDones]

Sorry for the late reply.
Without problem I can upload both things, binary mask after processing and the function also.

You are absolutely right about the seeds. We are trying to find the maxima after distance transform but it doesnt seem to look okey. On the other hand, apparently when watershed algorithm does the region growing of the seeds, although these are rigth it divides the objects in a weird fashion.

function processing_nuclei_multiThreading(img,radius)
  output = Array{Gray{Bool}}(undef, size(img)...)
  Threads.@threads for t in 1:size(output)[3]
      img2 = img[:, :, t]
      img2 = mapwindow(median, img2, Int.((radius*2+1,radius*2+1)))
      duplicate2 = copy(img2)
      duplicate2 = imfilter(duplicate2, Kernel.gaussian((radius*2,radius*2)))
      duplicate2 = Gray{N0f16}.(duplicate2)
      duplicate2 = mapwindow(minimum, duplicate2, Int.((radius*4+1,radius*4+1)))
      duplicate2 = mapwindow(maximum, duplicate2, Int.((radius*4+1,radius*4+1)))
      result = img2 .- duplicate2
      binary = Gray.(result .> 0)
      binary = mapwindow(minimum, binary, Int.((ceil((radius/2)+3.5),ceil((radius/2)+3.5))))
      output[:, :, t] = mapwindow(maximum, binary, Int.((ceil((radius/2)+3.5),ceil((radius/2)+3.5))))
  end
  return output
end

This is the function although its more like a pipeline of filters.

https://drive.google.com/file/d/1jFnO8WgeobxeqiunM3PPpI1_OuFwEVqh/view?usp=share_link
And here is the link to the binary after processing, I cannot directly upload tif files.

[DiegoPerezDones]

@tlnagy @jsundram @kmsquire @ViralBShah I would be really grateful if any of you could help us.

Thanks in advance!

[DiegoPerezDones]

Sorry for reopening this issue. We are right now trying to implement a whole 3D segementation algorithm in Julia and this step is crucial for us. If any of you could help us I would be very grateful.