Implementación de algoritmos técnicos en papel en C ++ o MATLAB

Soy un estudiante de Ingeniería Eléctrica. He leído muchos documentos técnicos sobre algoritmos de procesamiento de señales e imágenes (reconstrucción, segmentación, filtrado, etc.). La mayoría de los algoritmos que se muestran en estos documentos se definen sobre el tiempo continuo y la frecuencia continua, y a menudo dan las soluciones en términos de ecuaciones complicadas. ¿Cómo implementaría un documento técnico desde cero en C ++ o MATLAB para replicar los resultados obtenidos en dicho documento?

Más específicamente, estaba mirando el documento "Un algoritmo general de reconstrucción de haz de cono" de Wang et al ( IEEE Trans Med Imaging. 1993; 12 (3): 486-96 ), y me preguntaba, ¿cómo puedo comenzar? implementando su algoritmo? La ecuación 10 le da la fórmula de la imagen reconstruida en. ¿Cómo codificarías eso? ¿Tendría un bucle for pasando por cada vóxel y calculando la fórmula correspondiente? ¿Cómo codificaría funciones de funciones en esa fórmula? ¿Cómo evaluaría las funciones en puntos arbitrarios?

He leído el libro "Procesamiento digital de imágenes" de González y Woods, pero todavía estoy perdido. También he leído sobre la serie de libros Recetas Numéricas. ¿Sería esa la forma correcta?

¿Cuáles son sus algoritmos de programación de experiencias de trabajos de investigación? ¿Algún consejo o sugerencia?

Los algoritmos de procesamiento de señales definidos en tiempo / espacio / frecuencia continuos se implementan típicamente muestreando la señal en una cuadrícula discreta y convirtiendo integrales en sumas (y derivadas en diferencias). Los filtros espaciales se implementan mediante convolución con un núcleo de convolución (es decir, suma ponderada de vecinos).

Existe un gran conocimiento sobre el filtrado de señales de dominio de tiempo muestreadas; los filtros en el dominio del tiempo se implementan como filtros de respuesta de impulso finito , donde la muestra de salida actual se calcula como una suma ponderada de las N muestras de entrada anteriores; o filtros de respuesta de impulso infinito, donde la salida de corriente es una suma ponderada de las entradas y salidas anteriores . Formalmente, los filtros de tiempo discreto se describen utilizando la transformación z , que es el análogo de tiempo discreto de la transformada de Laplace . La transformación bilineal se mapea una a la otra ( c2dy d2cen Matlab).

¿Cómo evaluaría las funciones en puntos arbitrarios?

Cuando necesita el valor de una señal en un punto que no se encuentra directamente en su cuadrícula de muestreo, interpola su valor desde puntos cercanos. La interpolación puede ser tan simple como elegir la muestra más cercana, calcular un promedio ponderado de las muestras más cercanas o ajustar una función analítica arbitrariamente complicada a los datos muestreados y evaluar esta función en las coordenadas necesarias. Interpolar en una cuadrícula más fina y uniforme es el muestreo ascendente . Si su señal original (continua) no contiene detalles (es decir, frecuencias) más finos que la mitad de la cuadrícula de muestreo, entonces la función continua puede reconstruirse perfectamente a partir de la versión muestreada (el teorema de muestreo Nyquist-Shannon ). Para ver un ejemplo de cómo puede interpolar en 2D, veainterpolación bilineal .

En Matlab puede usar interp1o interp2para interpolar datos 2D 1D o muestreados regularmente (respectivamente), o griddatapara interpolar datos 2D muestreados irregularmente.

¿Tendría un bucle for pasando por cada vóxel y calculando la fórmula correspondiente?

Sí exactamente.

Matlab le evita tener que hacer esto a través de bucles for explícitos porque está diseñado para operar en matrices y vectores (es decir, matrices multidimensionales). En Matlab esto se llama "vectorización". Integrales definidas se pueden aproximar con sum, cumsum, trapz, cumtrapz, etc.

He leído el libro "Procesamiento digital de imágenes" de González y Woods, pero todavía estoy perdido. También he leído sobre la serie de libros Recetas Numéricas. ¿Sería esa la forma correcta?

Sí, las recetas numéricas serían un gran comienzo. Es muy práctico y cubre la mayoría de los métodos numéricos que terminarás necesitando. (Encontrará que Matlab ya implementa todo lo que necesita, pero las Recetas Numéricas le proporcionarán una excelente experiencia).

He tomado una clase de "algoritmos y estructuras de datos", pero no veo la relación entre el material presentado allí y la implementación de algoritmos científicos.

El material tratado en cursos de "Algoritmos y estructuras de datos" tiende a concentrarse en estructuras como listas, matrices, árboles y gráficos que contienen enteros o cadenas y operaciones como ordenar y seleccionar: problemas para los que generalmente hay un único resultado correcto. Cuando se trata de algoritmos científicos, esto es solo la mitad de la historia. La otra mitad se refiere a métodos para estimar números reales y funciones analíticas. Encontrará esto en un curso sobre "Métodos numéricos" (o "Análisis numérico"; como este- desplácese hacia abajo para las diapositivas): cómo estimar funciones especiales, cómo estimar integrales y derivadas, etc. Aquí una de las tareas principales es estimar la precisión de su resultado, y un patrón común es iterar una rutina que mejore estimar hasta que sea lo suficientemente preciso. (Puede preguntarse cómo sabe Matlab cómo hacer algo tan simple como estimar un valor sin(x)para algunos x).

Como un simple ejemplo, aquí hay un script corto que calcula una transformación de radón de una imagen en Matlab. La transformación de radón toma proyecciones de una imagen sobre un conjunto de ángulos de proyección. En lugar de tratar de calcular la proyección a lo largo de un ángulo arbitrario, en lugar de eso giro toda la imagen imrotate, de modo que la toma de proyección sea siempre vertical. Entonces podemos tomar la proyección simplemente usando sum, ya que el sumde una matriz devuelve un vector que contiene la suma sobre cada columna.

Puedes escribir el tuyo imrotatesi lo prefieres, usando interp2.

%%# Home-made Radon Tranform

%# load a density map (image).  
A = phantom;

n_pixels = size(A, 1);  %# image width (assume square)

%# At what rotation angles do we want to take projections?
n_thetas = 101;
thetas = linspace(0, 180, n_thetas);

result = zeros(n_thetas, n_pixels);

%# Loop over angles
for ii=1:length(thetas)
    theta = thetas(ii);
    rotated_image = imrotate(A, theta, 'crop');
    result(ii, :) = sum(rotated_image);
end

%# display the result
imagesc(thetas, 1:n_pixels, result.');
xlabel('projection angle [degrees]');

Lo que una vez fue una integral de la densidad a lo largo de un rayo ahora es una suma sobre una columna de una imagen muestreada discretamente, que a su vez se encontró al interpolar la imagen original sobre un sistema de coordenadas transformado.

Agregando a la excelente explicación de nibot , solo un par de puntos más.

• Los entornos informáticos numéricos como MATLAB, Octave o SciPy / NumPy le ahorrarán mucho esfuerzo en comparación con hacerlo solo en un lenguaje de programación genérico como C ++. Hacer malabarismos con doublematrices y bucles simplemente no se compara con tener tipos de datos como números complejos y operaciones como integrales a su alcance. (Es seguro, y un buen código C ++ puede ser un orden de magnitud más rápido, con buenas abstracciones de biblioteca y plantillas, incluso puede ser razonablemente limpio y claro, pero definitivamente es más fácil comenzar, por ejemplo, MATLAB).

• MATLAB también tiene "juegos de herramientas" para, por ejemplo, el procesamiento de imágenes y el procesamiento de señales digitales , que pueden ser de gran ayuda, dependiendo de lo que haga.

• El procesamiento digital de señales de Mitra es un buen libro para aprender (¡en MATLAB!) Los conceptos básicos de tiempo discreto, filtros, transformaciones, etc., que es un conocimiento bastante obligatorio para hacer cualquier algoritmo técnico decente.

Métodos numéricos. Suele ser un curso universitario de división superior y un libro de texto.

El DSP generalmente está cerca de la intersección de los métodos numéricos y la implementación eficiente. Si ignora la eficiencia, entonces lo que podría estar buscando es cualquier método de aproximación numérica que pueda producir un resultado "lo suficientemente preciso" para la (s) ecuación (es) de interés del documento técnico. A veces, uno podría estar tratando con datos muestreados, donde los teoremas de muestreo pondrán algunos límites tanto en el método de adquisición de datos (prefiltrado) como en el rango o la calidad de los resultados que puede obtener con esos datos.

A veces, Matlab, recetas numéricas o varias bibliotecas de procesamiento de imágenes / señales tendrán algoritmos o códigos eficientes para la solución numérica deseada. Pero a veces puede que tenga que rodar el suyo, por lo que es útil conocer las matemáticas detrás de varios métodos de solución numérica. Y ese es un gran tema por derecho propio.