Simulando la formación de imágenes#
Show code cell content
%load_ext autoreload
%autoreload 2
# Default imports
import sys
sys.path.append('../../../')
from helpers import *
from matplotlib import pyplot as plt
from myst_nb import glue
import numpy as np
from scipy import ndimage
Introducción#
La macro Diferencia de Gaussianas desarrollada en ImageJ: escritura de macros fue útil, pero bastante simple. Esta sección contiene una práctica mas amplia, cuyo objetivo es desarrollar una macro algo más sofisticada que tome una imagen «ideal» y luego simule cómo se vería después de ser registrada por un microscopio de fluorescencia. Puede utilizarse no sólo para comprender mejor el proceso de formación de imágenes, sino también para generar datos de prueba para algoritmos de análisis. Al crear simulaciones con diferentes configuraciones, podemos investigar cómo nuestros resultados podrían verse afectados por cambios en la adquisición y la calidad de la imagen.
Resumen de formación de imágenes#
El siguiente es un resumen de los aspectos de la formación de imágenes discutidos en este libro:
Las imágenes se componen de píxeles, cada uno de los cuales tiene un valor numérico único (¡no un color!).
El valor de un píxel en microscopía de fluorescencia se relaciona con una cantidad de fotones detectados o, más técnicamente, la carga de los electrones producidos por los fotones que chocan contra un detector.
Las imágenes pueden tener muchas dimensiones. El número de dimensiones es esencialmente el número de cosas que necesitas saber para identificar cada píxel (por ejemplo, punto temporal, número de canal, coordenada x, coordenada y, segmento z).
Los dos factores principales que limitan la calidad de la imagen son el desenfoque y el ruido. Ambos son inevitables y ninguno puede superarse por completo.
El desenfoque se caracteriza por la función de dispersión de puntos (PSF) del microscopio, que es el volumen 3D que sería el resultado de obtener imágenes de un único punto emisor de luz. Actúa como una convolución.
En el plano focal, la imagen de un punto es un patrón Airy. La mayor parte de la luz está contenida dentro de una región central, el disco Airy.
La resolución espacial es una medida de la separación que debe existir entre estructuras antes de que puedan distinguirse adecuadamente como separadas, y se relaciona con el tamaño del PSF (o disco de Airy en 2D).
Los dos tipos principales de ruido son el ruido de fotones y el ruido de lectura. El primero es causado por la aleatoriedad de la emisión de fotones, mientras que el segundo surge de imprecisiones en la cuantificación del número de fotones detectados.
Detectar más fotones ayuda a superar los problemas causados por ambos tipos de ruido.
Los distintos tipos de microscopio tienen diferentes ventajas y costes en términos de información espacial, resolución temporal y ruido.
Los PMT se utilizan para detectar fotones de un solo píxel, mientras que los CCD y EMCCD se utilizan para detectar fotones de muchos píxeles en paralelo.
La macro de este capítulo funcionará para imágenes 2D y simulará los tres componentes principales:
el desenfoque del PSF
ruido de fotones
ruido de lectura
Además, la macro finalmente se escribirá de tal manera que nos permita investigar los efectos de cambiar algunos parámetros adicionales:
el tamaño del PSF (relacionado con la apertura numérica de la lente del objetivo y el tipo de microscopio)
la cantidad de fluorescencia que se emite desde la región más brillante
la cantidad de fondo (de luz desviada y otras fuentes)
el tiempo de exposición (y por lo tanto el número de fotones detectados)
compensación del detector
la ganancia del detector
El ruido de lectura del detector
binning (o agrupación) de cámara
profundidad de bits
Grabar los pasos principales.#
Realmente no importa qué imagen uses para esto, pero recomiendo una imagen 2D de un solo canal que comience sin ningún ruido obvio (por ejemplo, la imagen Happy cell.tif). Después de iniciar la grabadora de macros, completa los siguientes pasos para crear la estructura principal de la macro:
Figura 163 Ejemplo de imagen de entrada (input)#
Figura 164 Ejemplo de imagen de salida (output)#
Asegúrate de que la imagen inicial sea de 32 bits.
Ejecuta usando un valor sigma de 2 para simular la convolución con el PSF.
Aquí asumiremos que hay algunos fotones de fondo de otras fuentes, pero aproximadamente el mismo número en cada píxel de la imagen. Así que simplemente podemos agregar una constante a esta imagen usando . El valor debería ser pequeño, quizás 10.
La imagen ahora contiene las «tasas promedio de emisión de fotones» que normalmente nos gustaría tener para un tiempo de exposición particular (es decir, está libre de ruido). Si cambiamos el tiempo de exposición, deberíamos cambiar los valores de píxeles de manera similar para que las tasas sigan siendo las mismas. Debido a que ajustar la exposición funciona como una multiplicación, podemos usar el comando . Configúralo en un valor “predeterminado” de 1 por ahora.
Para convertir las tasas de emisión de fotones en recuentos de fotones reales que potencialmente podríamos detectar, necesitamos simular el ruido de los fotones reemplazando cada píxel por un valor aleatorio de una distribución de Poisson que tenga el mismo \(\lambda\) que la tasa misma. Al momento de escribir este artículo, no hay ningún comando incorporado en ImageJ o Fiji para simular el ruido de Poisson, pero podemos instalar RandomJ siguiendo las instrucciones en https://imagescience.org/meijering/software/randomj/. Luego agrega ruido llamando a , asegurándote de configurar el valor de Insertion: a Modulatory. El valor Mean se ignorará, por lo que su configuración no importa.
Observa que todos los píxeles ahora deberían tener valores enteros: no puedes detectar partes de fotones.
La ganancia del detector aumenta la cantidad de electrones producidos por los fotones detectados. Hazle espacio incluyendo otra multiplicación, aunque por ahora establece el valor en 1 (es decir, sin ganancia adicional).
Simula la compensación del detector agregando otra constante, p.e. 100.
Agrega ruido de lectura con , estableciendo la desviación estándar en 5.
Recorta cualquier valor negativo ejecutando y definiendo el valor en 0.
Recorta cualquier valor positivo que exceda la profundidad de bits ejecutando . Para asumir una imagen de 8 bits, define el valor en 255.
Ahora es un buen momento para limpiar el código eliminando líneas innecesarias, agregando comentarios adecuados y colocando cualquier variable interesante en la parte superior de la macro para que pueda modificarse fácilmente más adelante (como en Diferencia de Gaussianos). También debemos duplicar la imagen para evitar modificar accidentalmente el original. El resultado final debería verse así:
// Variables to change
psfSigma = 2;
backgroundPhotons = 10;
exposureTime = 1;
readStdDev = 5;
detectorGain = 1;
detectorOffset = 100;
maxVal = 255;
// Duplicate the image, to avoid changing the original
// Including " " means the default name will be used &
// we don't need to specify any other
run("Duplicate...", " ");
// Retain a reference so we can close the duplicate at the end
// (since RandomJ Poisson creates another new window)
idTemp = getImageID();
// Ensure image is 32-bit
run("32-bit");
// Simulate PSF blurring
run("Gaussian Blur...", "sigma="+psfSigma);
// Add background photons
run("Add...", "value="+backgroundPhotons);
// Multiply by the exposure time
run("Multiply...", "value="+exposureTime);
// Simulate photon noise
run("RandomJ Poisson", "mean=1.0 insertion=modulatory");
// Simulate the detector gain
run("Multiply...", "value="+detectorGain);
// Simulate the detector offset
run("Add...", "value="+detectorOffset);
// Simulate read noise
run("Add Specified Noise...", "standard="+readStdDev);
// Clip any negative values
run("Min...", "value=0");
// Clip the maximum values based on the bit-depth
un("Max...", "value="+maxVal);
Tendrías una macro perfectamente respetable si te detuvieras ahora, pero la siguiente sección contiene algunas formas en las que se puede mejorar.
Haciendo mejoras#
Normalizando la imagen#
Los resultados que obtengas al ejecutar la macro anterior cambiarán según el rango original de la imagen que utilices: es decir, una imagen que comienza con píxeles de alto valor terminará teniendo mucho menos ruido. Para compensar esto un poco, primero podemos normalizar la imagen para que todos los píxeles estén en el rango 0–1. Para hacer esto, necesitamos determinar el rango actual de valores de píxeles, que se puede encontrar usando la función macro:
getStatistics(area, mean, min, max);
Después de ejecutar esto, se crean cuatro variables que proporcionan los valores de píxeles «media», «mínimo» y «máximo» en la imagen, junto con el «área» total de la imagen. La normalización ahora es posible usando los comandos y y ajustando sus valores. Al final esto nos da
getStatistics(area, mean, min, max);
run("Subtract...", "value="+min);
divisor = max - min;
run("Divide...", "value="+divisor);
Variación de la tasa de emisión de fluorescencia.#
El nuevo problema que tendremos tras la normalización es que habrá una tasa máxima de emisión de fotones de 1 en la parte más brillante de la imagen, lo que nos dará una imagen dominada completamente por el ruido. Esto lo podemos cambiar multiplicando de nuevo los píxeles, y así definir cuál queremos que sea la tasa de emisión en la parte más brillante de la imagen. Sugiero crear una variable para esto y definir su valor en 10. Luego agrega la siguiente línea inmediatamente después de la normalización:
run("Multiply...", "value="+maxPhotonEmission);
Modificar este valor te permite cambiar entre simular muestras que tienen una fluorescencia más o menos intensa. Para una muestra menos brillante, lo más probable es que necesites aumentar el tiempo de exposición para obtener una cantidad similar de señal, pero ten en cuenta que aumentar el tiempo de exposición también implica recolectar más fotones de fondo inútiles, por lo que no es tan bueno como tener una muestra donde las partes importantes son intrínsecamente más brillantes.
Simulando agrupación o binning#
La idea principal de binning es que los electrones de múltiples píxeles se suman antes de la lectura, de modo que la cantidad de electrones que se cuantifican sea mayor en relación con el ruido de lectura. Para la agrupación 2×2, esto implica dividir la imagen en distintos bloques de 2×2 píxeles y crear otra imagen en la que el valor de cada píxel sea la suma de los valores dentro del bloque correspondiente.
Esto se podría hacer usando el comando , con un shrink factor de 2 y el método bin Sum. La grabadora de macros se puede utilizar nuevamente para obtener el código principal que se necesita. Después de algunas modificaciones, esto se convierte en
if (doBinning) {
run("Bin...", "x=2 y=2 bin=Sum");
}
Al encerrar la línea dentro de un bloque de código (limitado por las llaves) y comenzar el bloque con if (doBinning), es fácil controlar si se aplica la agrupación o no. Simplemente agrega una variable adicional a tu lista al inicio de la macro
doBinning = true;
para activar la agrupación, o
doBinning = false;
para apagarlo. Las líneas de código que realizan la agrupación deben insertarse antes de agregar ruido de lectura.
Profundidad de bits variables#
Variar las profundidades de bits simuladas cambiando el valor máximo permitido en la imagen requiere un poco de trabajo: necesitas saber que el valor máximo en una imagen de 8 bits es 255, mientras que para una imagen de 12 bits es 4095 y así sucesivamente. Es más intuitivo simplemente cambiar la profundidad de bits de la imagen y dejar que la macro haga el cálculo por ti. Para hacer esto, puedes reemplazar la variable maxVal = 255; al inicio de la macro con nBits = 8; y luego actualizar el código de recorte posterior para que se convierta en
maxVal = pow(2, nBits) - 1;
run("Max...", "value="+maxVal);
Aquí, pow(2, nBits) es una función que le da el valor de 2nBits. Ahora es más fácil explorar la diferencia entre imágenes de 8, 12 y 14 bits (que son las principales profundidades de bits normalmente asociadas con los detectores de microscopio, incluso si la imagen resultante se almacena como 16 bits).
Redondeo a valores enteros#
La macro ya ha recortado la imagen a una profundidad de bits específica, pero todavía contiene datos de 32 bits y, por lo tanto, potencialmente tiene valores no enteros que no se pueden almacenar en las imágenes de 8 o 16 bits que un microscopio normalmente proporciona como output. Por lo tanto, queda redondear los valores al número entero más cercano.
Hay algunas formas de hacer esto: podemos convertir la imagen usando los comandos , aunque luego debemos tener cuidado sobre si aplicaremos alguna escala. Sin embargo, podemos evitar pensar en esto si aplicamos el redondeo nosotros mismos. Para hacerlo necesitamos visitar cada píxel, extraer su valor, redondear el valor al número entero más cercano y volver a colocarlo en la imagen. Esto requiere el uso de bucles. El código, que debe agregarse al final de la macro, tiene este aspecto:
// Get the image dimensions
width = getWidth();
height = getHeight();
// Loop through all the rows of pixels
for (y = 0; y < height; y++) { // Loop through all the columns of pixels
for (x = 0; x < width; x++) { // Extract the pixel value at coordinate (x, y)
value = getPixel(x, y); // Round the pixel value to the nearest integer
value = round(value); // Replace the pixel value in the image
setPixel(x, y, value);
}
}
Esto crea dos variables, x e y, que se utilizan para almacenar las coordenadas horizontales y verticales de un píxel. Cada uno comienza establecido en 0 (por lo que comenzamos con el píxel en 0,0, es decir, en la parte superior izquierda de la imagen). El código en el medio se ejecuta para establecer el valor del primer píxel, luego la variable x se incrementa hasta convertirse en 1, porque x++ significa sumar 1 a x. Este proceso se repite siempre que «x» sea menor que el ancho de la imagen, «x <ancho». Cuando x finalmente se vuelve igual al ancho, significa que todos los valores de píxeles en la primera fila de la imagen se han redondeado. Luego, se incrementa «y» y «x» se restablece a cero, antes de que el proceso se repita y la siguiente fila también se redondee. Esto continúa hasta que y sea igual a la altura de la imagen, momento en el que se completa el procesamiento [^fn_loops].
Código final#
El código final de mi versión de la macro se proporciona a continuación:
// Variables to change
psfSigma = 2;
backgroundPhotons = 10;
exposureTime = 10;
readStdDev = 5;
detectorGain = 1;
detectorOffset = 100;
nBits = 8;
maxPhotonEmission = 10;
doBinning = false;
// Duplicate the image, to avoid changing the original
// Including " " means the default name will be used &
// we don't need to specify any other
run("Duplicate...", " ");
// Retain a reference so we can close the duplicate at the end
// (since RandomJ Poisson creates another new window)
idTemp = getImageID();
// Ensure image is 32-bit
run("32-bit");
// Normalize the image to the range 0-1
getStatistics(area, mean, min, max);
run("Subtract...", "value="+min);
divisor = max - min;
run("Divide...", "value="+divisor);
// Define the photon emission at the brightest point
run("Multiply...", "value="+maxPhotonEmission);
// Simulate PSF blurring
run("Gaussian Blur...", "sigma="+psfSigma);
// Add background photons
run("Add...", "value="+backgroundPhotons);
// Multiply by the exposure time
run("Multiply...", "value="+exposureTime);
// Simulate photon noise
run("RandomJ Poisson", "mean=1.0 insertion=modulatory");
// Simulate the detector gain
// (note this should really add Poisson noise too!)
run("Multiply...", "value="+detectorGain);
// Simulate binning (optional)
if (doBinning) {
run("Bin...", "x=2 y=2 bin=Sum");
}
// Simulate the detector offset
run("Add...", "value="+detectorOffset);
// Simulate read noise
run("Add Specified Noise...", "standard="+readStdDev);
// Clip any negative values
run("Min...", "value=0");
// Clip the maximum values based on the bit-depth
maxVal = pow(2, nBits) - 1;
run("Max...", "value="+maxVal);
// Get the image dimensions
width = getWidth();
height = getHeight();
// Round the pixels to integer values
for (y = 0; y < height; y++) {
// Loop through all the columns of pixels
for (x = 0; x < width; x++) {
// Extract the pixel value at coordinate (x, y)
value = getPixel(x, y);
// Round the pixel value to the nearest integer
value = round(value);
// Replace the pixel value in the image
setPixel(x, y, value);
}
}
// Reset the display range (i.e. image contrast)
resetMinAndMax();
// Close the temporary image
selectImage(idTemp);
close();
Usos y limitaciones#
Por supuesto, la macro anterior se basa en algunas suposiciones y simplificaciones. Por ejemplo, trata la ganancia como una simple multiplicación del número de fotones, pero el proceso de amplificación de la ganancia también implica cierta aleatoriedad, lo que introduce ruido adicional. Debido a que este ruido se comporta estadísticamente como el ruido de los fotones, se puede considerar que el efecto disminuye la cantidad de fotones que se detectaron. Además, hemos tratado el fondo como una constante que es la misma en todas partes de una imagen. En la práctica, el fondo normalmente consiste principalmente en luz desenfocada de otros planos de la imagen, por lo que realmente debería cambiar en diferentes partes de la imagen, particularmente en el caso de campo amplio.
Sin embargo, se han tenido en cuenta bastantes factores. Al explorar diferentes combinaciones de configuraciones, puedes tener una idea de cómo afectan la calidad general de la imagen. Por ejemplo, podrías intentar:
Aumentar el fondo, manteniendo igual la emisión máxima de fotones
Eliminar la compensación del detector o establecerlo en un valor negativo
Comparar los efectos de la agrupación para imágenes con recuentos de fotones altos y bajos
Crear varias imágenes a partir de los mismos datos de origen y luego promediarlas juntas para ver cómo cambia el ruido.
Cuando se planea implementar alguna estrategia de análisis, especialmente si se realizan mediciones de intensidad de fluorescencia, también puede resultar útil probar su eficacia utilizando esta macro. Para hacerlo, necesitarías
Crear de alguna manera una imagen de ejemplo «perfecta», sin ruido ni desenfoque, ya sea manualmente o desconvolucionando una imagen de muestra similar
Aplica tu algoritmo a esta imagen perfecta para saber qué detecta y qué conclusiones podrías sacar
Aplica exactamente el mismo algoritmo a una versión de la imagen que pasó por el simulador y ve cuán diferentes serían tus mediciones y conclusiones.
Idealmente, los resultados deberían ser los mismos en ambos casos: esto sugeriría que tu método de análisis es sólido y puede manejar imágenes con diferentes niveles de calidad. Sin embargo, es más probable que los resultados sean diferentes. En ese caso, la comparación te da una idea de cuán afectadas están tus mediciones por el proceso de obtención de imágenes y, por lo tanto, de que tan confiable es la relación con la muestra subyacente «real».