I'm analyzing the texture of cell nuclei in fluorescence microscopy images. I’ve attached two example nuclei: one that is less bright but has strong texture — full of aggregates or “granular” patterns — and another one that is much brighter but more homogeneous.
I’ve been trying to quantify this difference using several texture measures in MATLAB, mainly based on the Gray-Level Co-occurrence Matrix (GLCM). I also varied the distance parameter since the aggregates are relatively large. However, I keep getting homogeneity and contrast values that are completely inverted compared to what I would expect (the granular nucleus gives higher homogeneity, for instance).
I also tried calculating entropy, which works slightly better, but still doesn’t consistently capture the difference across all nuclei.
The visual difference between the two images is so obvious that it’s really frustrating not to be able to translate it into quantitative values.
I’m attaching some .tif examples images (including the nuclei of the "problem" in case anyone wants to test directly with them, and below is the part of my MATLAB code that loops through each nucleus and calculates the texture parameters.
Any recommendations? I'm lost!
Thanks in advance for your help!
%----------------------------------------------------------------------------------------------------------------------------------------------
% CÁLCULO DE LOS PARÁMATROS
%--------------------------------------------------------------------------------------------------------------------------------------------
% Calcular el número de nucleos de la imagen en la que estamos trabajando en esta ronda
nnuclei=max(In(:));
% Inicializar vectores para almacenar los datos extraídos
intensities= zeros(1, nnuclei);
contrasts = zeros(1, nnuclei);
homogeneities = zeros(1, nnuclei);
% Para cada uno de los núcleos de la imagen
for n = 1:nnuclei
count=count+1;
% Crear una máscara para el núcleo actual
nuclei_n = (In == n); %figure; imshow(nuclei_n)
%------------------------------------FLUORESCENCE INTENSITY--------------------------------------------------------------------------------------------
% Extraer los valores de intensidad de la imagen para este núcleo
nuclei_n_intensities = Idapi(nuclei_n);
intensities(n) = mean(nuclei_n_intensities);
%------------------------------------TEXTURA--------------------------------------------------------------------------------------------
%La matriz de coocurrencia (GLCM) es una tabla que registra cuántas veces aparecen juntas diferentes combinaciones de tonos de gris entre píxeles vecinos.
% Por ejemplo, cuenta cuántas veces un píxel oscuro está al lado de uno claro o cuántas veces dos tonos medios coinciden.
% Al analizar esa distribución de combinaciones, se obtiene una medida objetiva de la textura de la imagen: si la mayoría de los píxeles vecinos tienen tonos similares,
% la imagen es uniforme, con alta homogeneidad y bajo contraste; en cambio, si abundan las combinaciones entre tonos muy distintos, la imagen es irregular o rugosa,
% con alto contraste y baja homogeneidad. En definitiva, la GLCM permite traducir una percepción visual —como que un núcleo parece más granular o más liso—
% en un valor numérico que describe la variación de intensidades dentro de esa región.
%Vamos a cuantizar la imagen a 32 niveles. Así la GLCM es manejable y comparable entre núcleos.
%La matriz de coocurrencia (GLCM) necesita que la imagen tenga un número limitado de niveles de gris para poder contar las combinaciones de intensidades vecinas de forma eficiente.
% Si se trabajara con todos los valores posibles (por ejemplo, 65.536 en una imagen de 16 bits), la matriz sería enorme y redundante, porque muchos tonos son muy parecidos.
% Por eso se cuantiza la imagen, es decir, se agrupan las intensidades en un número reducido de categorías o "cajones”.
% Elegir L = 32 significa dividir el rango de grises en 32 intervalos y asignar cada píxel a uno de ellos.
% Este valor ofrece un buen equilibrio: mantiene suficiente detalle para distinguir texturas diferentes (como núcleos con y sin agregados) pero evita que la matriz sea demasiado
% grande o sensible al ruido, lo que hace que las medidas de contraste y homogeneidad sean más estables y comparables entre imágenes
L = 32;
% Recortamos la imagen sólo a ese nucleo en concreto. Buscamos el rectángulo mínimo que lo contiene: la primera y última fila/columna donde hay píxeles del núcleo.
%Luego cortamos tanto la imagen original como la máscara de los núcleos
[min_row,max_row] = deal(find(any(nuclei_n, 2),1,'first'), find(any(nuclei_n, 2),1,'last'));
[min_col,max_col] = deal(find(any(nuclei_n, 1),1,'first'), find(any(nuclei_n, 1),1,'last'));
In_crop = Idapi(min_row:max_row, min_col:max_col);
nuclei_n_crop = nuclei_n(min_row:max_row, min_col:max_col);
% Cuantizar SOLO los píxeles del núcleo a L niveles
In_q = gray2ind(mat2gray(In_crop), L); figure(image);imshow(In_crop, [])
In_q(~nuclei_n_crop) = 0; % fuera del núcleo a 0 (valor inválido)
%Una GLCM depende de una dirección (qué vecino miramos). Para no depender de la orientación de la textura, usamos 4 direcciones típicas.
dists = [5 7 12 15 20]; % prueba 1–4 px (ajusta a tu tamaño de gránulo)
offsets = [];
for d = dists
offsets = [offsets; 0 d; -d d; d 0; d d]; %0°(vecino a la derecha), 45° (diagonal superior derecha), 90° (debajo), 135° (diagonal inferior derecha)
end
% glcm_sum es la matriz acumulada donde iremos sumando la GLCM de cada dirección.
glcm_sum = zeros(L,L);
%para cada una de las direcciones
for k = 1:size(offsets,1)
off = offsets(k,:); %es el desplazamiento del vecino que miramos (por ejemplo, [0 1] = vecino de la derecha).
% Píxeles válidos donde el par (p, p+off) cae dentro del núcleo
Mvalid = nuclei_n_crop & circshift(nuclei_n_crop, [-off(1), -off(2)]); %El & dice: válido solo donde los dos (píxel y vecino) están dentro del núcleo. Así evitas pares con fondo.
A = In_q; B = circshift(In_q, [-off(1), -off(2)]); % circshift(Mcrop, -off...) mueve la máscara como si alinearas cada píxel con su vecino.
A = A(Mvalid); B = B(Mvalid);
% filtrar: índices deben ser enteros en 1..L (sin ceros ni NaN)
idx = A>=1 & A<=L & B>=1 & B<=L & isfinite(A) & isfinite(B);
if ~any(idx)
continue % no hay pares válidos para este offset (núcleo pequeño)
end
A = A(idx); B = B(idx);
%Aquí hacemos literalmente una tabla de frecuencias en las filas pones el tono del píxel, en las columnas el tono del vecino,
% y cada vez que aparece esa pareja (por ejemplo, “nivel 5 con nivel 7”) sumas 1.
% Esa tabla M es la GLCM para esa dirección. Imagina que cuentas: "He visto 10 veces un píxel gris medio (nivel 5) al lado de uno claro (nivel 7)." → En la celda (5,7) de la tabla pones un 10.
M = accumarray([A(:), B(:)], 1, [L L], @sum, 0);
%La sumamos en glcm_sum para acumular las 4 direcciones.En cada vuelta del bucle se analiza una dirección distinta (derecha, diagonal, abajo, etc.)
% y se construye una tabla que cuenta cuántas veces aparecen juntas las distintas combinaciones de tonos de gris entre píxeles vecinos.
% Cada una de esas tablas se va sumando en glcm_sum, que al final reúne la información de todas las direcciones.
% Es como recorrer el núcleo con una lupa, anotando qué tonos aparecen juntos al mirar en cada orientación y, después, juntar todas esas notas en una única tabla.
% Esa tabla global describe todas las relaciones de tonos posibles dentro del núcleo y permite calcular medidas como el contraste y la homogeneidad de su textura.
% Al final en la casilla 1:1 tendré cuantas veces han aparecido los valores de intensidad de 1 y 1 juntos (en cualquier dirección) y así con el 1:2, 1:3, 1:32...2:1, 2:2, 2:32....todos
glcm_sum = glcm_sum + M;
end
% Extraer las propiedades que nos interesan
props = graycoprops(glcm_sum, {'Contrast','Homogeneity'});
contrasts(n) = props.Contrast;
homogeneities(n) = props.Homogeneity;
