[Python] Skeletonized microglial cell branches not reconnecting after binarization, and...

I'm working on a Python script for processing immunohistochemistry microscopy images of microglial cells. The goal is to:

1. Extract the green channel.
2. Apply CLAHE for contrast enhancement.
3. Perform Otsu thresholding to binarize the image.
4. Skeletonize the binary mask.
5. Reconnect fragmented branches (that were cut during binarization) based on their proximity and direction.
6. Remove background structures that are not connected to a microglial soma (i.e., eliminate isolated or floating branches).

Despite using a function to reconnect endpoints (rellenar_segmentos_faltantes), the branches are still not reconnecting properly, even when visually aligned and close in the original image.

Also, the background removal step does not seem to eliminate unconnected or floating structures, especially long branches not connected to any cell body.

My code is the following:

import os
import cv2
import numpy as np
import matplotlib.pyplot as plt
from skimage import io, img_as_ubyte, img_as_float, morphology, measure
from skimage.morphology import disk, skeletonize, remove_small_objects
from skimage.filters import threshold_otsu
from skimage.io import imsave
from scipy.ndimage import convolve
from skimage.draw import line

# === CONFIGURACIÓN ===
input_dir = "xxxxxx"
output_dir = "resultado_celulas_filtradas"
os.makedirs(output_dir, exist_ok=True)

def show_images(title, images, titles, cmap='gray'):
fig, axs = plt.subplots(1, len(images), figsize=(4*len(images), 4))
fig.suptitle(title, fontsize=16)
for ax, img, t in zip(axs, images, titles):
ax.imshow(img, cmap=cmap)
ax.set_title(t)
ax.axis('off')
plt.tight_layout()
plt.show()

def find_endpoints(skel):
kernel = np.array([[1,1,1], [1,10,1], [1,1,1]])
filtered = convolve(skel.astype(np.uint8), kernel, mode='constant', cval=0)
return np.logical_and(skel, filtered == 11)

def get_direction_vector(skel, y, x):
neighbors = [(-1,0), (1,0), (0,-1), (0,1), (-1,-1), (-1,1), (1,-1), (1,1)]
for dy, dx in neighbors:
ny, nx = y + dy, x + dx
if 0 <= ny < skel.shape[0] and 0 <= nx < skel.shape[1]:
if skel[ny, nx]:
return np.array([dy, dx])
return None

def rellenar_segmentos_faltantes(skeleton, max_dist=25, max_angle=30):
endpoints = np.argwhere(find_endpoints(skeleton))
conectadas = np.zeros_like(skeleton, dtype=bool)
for i in range(len(endpoints)):
y1, x1 = endpoints
v1 = get_direction_vector(skeleton, y1, x1)
if v1 is None:
continue
for j in range(i+1, len(endpoints)):
y2, x2 = endpoints[j]
dist_euc = np.hypot(y2 - y1, x2 - x1)
if dist_euc < max_dist:
v2 = np.array([y2 - y1, x2 - x1])
cos_angle = np.dot(v1, v2) / (np.linalg.norm(v1) * np.linalg.norm(v2) + 1e-5)
angle = np.degrees(np.arccos(np.clip(cos_angle, -1.0, 1.0)))
if angle < max_angle:
rr, cc = line(y1, x1, y2, x2)
conectadas[rr, cc] = True
return np.logical_or(skeleton, conectadas)

# === PROCESAMIENTO ===
for filename in os.listdir(input_dir):
if filename.lower().endswith(('.jpg', '.jpeg', '.tif', '.tiff', '.png')):
path = os.path.join(input_dir, filename)
print(f"\nProcesando: {filename}")
img_rgb = io.imread(path)

# 1. Usar canal verde
canal_verde = img_rgb[:, :, 1]
img_gray = img_as_ubyte(img_as_float(canal_verde))

# 2. CLAHE
clahe = cv2.createCLAHE(clipLimit=2.0, tileGridSize=(8, 8))
enhanced = clahe.apply(img_gray)

# 3. Umbral de Otsu
otsu_val = threshold_otsu(enhanced)
binarizada = (enhanced > otsu_val).astype(np.uint8)

# 4. Cierre morfológico
cerrada = morphology.binary_closing(binarizada, disk(2))

# 5. Esqueletizar y rellenar segmentos faltantes con dirección
esqueleto = skeletonize(cerrada)
esqueleto_conectado = rellenar_segmentos_faltantes(esqueleto, max_dist=25, max_angle=30)

# 6. Limpiar el fondo eliminando estructuras muy cortas
limpio = remove_small_objects(esqueleto_conectado.astype(bool), min_size=15)

# 7. Combinar con máscara original cerrada
final = np.logical_or(cerrada, limpio).astype(np.uint8)

# 8. Visualización
show_images("Reconstrucción conectada con limpieza de fragmentos",
[img_gray, enhanced, binarizada, final],
["Canal verde", "CLAHE", "Binarizada", "Final conectada"])

# 9. Guardar
nombre_base = os.path.splitext(filename)[0]
output_path = os.path.join(output_dir, f"filtrada_{nombre_base}.tif")
imsave(output_path, img_as_ubyte(final))
print(f"Guardado: {output_path}")

1. How can I improve the reconnection of skeletonized microglial branches that appear cut after binarization, even though they are visually aligned?
2. What would be a better approach to remove background filaments not connected to a microglial soma (likely the cell body)? I want to retain only those structures attached to a central region.
3. Is there a more robust way to identify and keep structures attached to a large binary region (the soma) and discard the rest?

I've attached example images showing the results of each processing step — including the original green channel, the CLAHE-enhanced version, the binarized image and the final output — to illustrate the issue more clearly.

Any suggestions or insights would be greatly appreciated! Thanks in advance.

enter image description here

Continue reading...