Implémentation d'algorithmes de papier technique en C ++ ou MATLAB

Je suis un étudiant en génie électrique. J'ai lu de nombreux articles techniques sur les algorithmes de traitement du signal et de l'image (reconstruction, segmentation, filtrage, etc.). La plupart des algorithmes présentés dans ces articles sont définis sur un temps continu et une fréquence continue, et donnent souvent des solutions en termes d'équations compliquées. Comment implémenteriez-vous un document technique à partir de zéro en C ++ ou MATLAB afin de reproduire les résultats obtenus dans ledit document?

Plus précisément, je regardais l'article «Un algorithme général de reconstruction de faisceaux coniques» de Wang et al ( IEEE Trans Med Imaging. 1993; 12 (3): 486-96 ), et je me demandais comment commencer. mettre en œuvre leur algorithme? L'équation 10 vous donne la formule de l'image reconstruite à. Comment coderiez-vous cela? Auriez-vous une boucle for passant par chaque voxel et calculant la formule correspondante? Comment coderiez-vous les fonctions des fonctions dans cette formule? Comment évalueriez-vous les fonctions à des points arbitraires?

J'ai lu le livre "Digital Image Processing" de Gonzalez et Woods mais je suis toujours perdu. J'ai également lu des articles sur la série de livres sur les recettes numériques. Serait-ce la bonne façon?

Quelles sont vos expériences de programmation d'algorithmes à partir d'articles de recherche? Des conseils ou des suggestions?

Les algorithmes de traitement du signal définis en temps / espace / fréquence continus sont généralement mis en œuvre en échantillonnant le signal sur une grille discrète et en convertissant les intégrales en sommes (et les dérivés en différences). Les filtres spatiaux sont mis en œuvre par convolution avec un noyau de convolution (c'est-à-dire la somme pondérée des voisins).

Il existe un énorme corpus de connaissances concernant le filtrage des signaux du domaine temporel échantillonnés; les filtres dans le domaine temporel sont mis en œuvre sous forme de filtres à réponse impulsionnelle finie , où l'échantillon de sortie actuel est calculé comme une somme pondérée des N échantillons d'entrée précédents; ou filtres à réponse impulsionnelle infinie, où la sortie de courant est une somme pondérée des entrées et sorties précédentes . Formellement, les filtres temporels discrets sont décrits en utilisant la transformée z , qui est l'analogue temporel discret de la transformée de Laplace . La transformation bilinéaire est mappée l'une à l'autre ( c2det d2cdans Matlab).

Comment évalueriez-vous les fonctions à des points arbitraires?

Lorsque vous avez besoin de la valeur d'un signal à un point qui ne se trouve pas directement sur votre grille d'échantillonnage, vous interpolez sa valeur à partir de points voisins. L'interpolation peut être aussi simple que de choisir l'échantillon le plus proche, de calculer une moyenne pondérée des échantillons les plus proches, ou d'adapter une fonction analytique arbitrairement compliquée aux données échantillonnées et d'évaluer cette fonction aux coordonnées nécessaires. L'interpolation sur une grille uniforme plus fine est un suréchantillonnage . Si votre signal d'origine (continu) ne contient pas de détails (c'est-à-dire des fréquences) plus fins que la moitié de la grille d'échantillonnage, alors la fonction continue peut être parfaitement reconstruite à partir de la version échantillonnée (le théorème d'échantillonnage de Nyquist-Shannon ). Pour voir un exemple d'interpolation en 2D, consultezinterpolation bilinéaire .

Dans Matlab, vous pouvez utiliser interp1ou interp2pour interpoler 1D ou des données 2D échantillonnées régulièrement (respectivement), ou griddatapour interpoler à partir de données 2D échantillonnées irrégulièrement.

Auriez-vous une boucle for passant par chaque voxel et calculant la formule correspondante?

Oui, exactement.

Matlab vous évite d'avoir à le faire via des boucles for explicites car il est conçu pour fonctionner sur des matrices et des vecteurs (c'est-à-dire des tableaux multidimensionnels). Dans Matlab, cela s'appelle la "vectorisation". Intégrales peuvent être approchées Definite avec sum, cumsum, trapz, cumtrapz, etc.

J'ai lu le livre "Digital Image Processing" de Gonzalez et Woods mais je suis toujours perdu. J'ai également lu des articles sur la série de livres sur les recettes numériques. Serait-ce la bonne façon?

Oui, les recettes numériques seraient un bon début. Il est très pratique et couvre la plupart des méthodes numériques dont vous aurez besoin. (Vous constaterez que Matlab implémente déjà tout ce dont vous avez besoin, mais les recettes numériques fourniront un excellent arrière-plan.)

J'ai pris un cours "algorithmes et structures de données", mais je ne vois pas la relation entre le matériel qui y est présenté et la mise en œuvre d'algorithmes scientifiques.

Le matériel traité dans les cours "Algorithmes et structures de données" a tendance à se concentrer sur des structures telles que des listes, des tableaux, des arbres et des graphiques contenant des nombres entiers ou des chaînes et des opérations comme le tri et la sélection: problèmes pour lesquels il n'y a généralement qu'un seul résultat correct. En ce qui concerne les algorithmes scientifiques, ce n'est que la moitié de l'histoire. L'autre moitié concerne les méthodes d'estimation des nombres réels et des fonctions analytiques. Vous le trouverez dans un cours sur les "Méthodes numériques" (ou "Analyse numérique"; comme celui-ci- faites défiler vers le bas pour les diapositives): comment estimer des fonctions spéciales, comment estimer des intégrales et des dérivées, etc. estimer jusqu'à ce qu'il soit suffisamment précis. (Vous pourriez vous demander comment Matlab sait comment faire quelque chose d'aussi simple que d'estimer une valeur de sin(x)pour certains x.)

À titre d'exemple simple, voici un petit script qui calcule une transformation radon d'une image dans Matlab. La transformation du radon prend des projections d'une image sur un ensemble d'angles de projection. Au lieu d'essayer de calculer la projection selon un angle arbitraire, je fais plutôt pivoter l'image entière à l'aide de imrotatesorte que la prise de projection soit toujours verticale. Ensuite, nous pouvons prendre la projection simplement en utilisant sum, car la summatrice renvoie un vecteur contenant la somme sur chaque colonne.

Vous pouvez écrire le vôtre imrotatesi vous préférez, en utilisant 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]');

Ce qui était autrefois une intégrale de la densité le long d'un rayon est maintenant une somme sur une colonne d'une image échantillonnée discrètement, qui a été trouvée à son tour en interpolant l'image originale sur un système de coordonnées transformé.

Pour ajouter à l' excellente explication de nibot , juste quelques points de plus.

• Les environnements informatiques numériques tels que MATLAB, Octave ou SciPy / NumPy vous feront économiser beaucoup d'efforts par rapport à tout faire par vous-même dans un langage de programmation générique tel que C ++. Jongler avec des doubletableaux et des boucles n'est tout simplement pas comparable à avoir des types de données tels que des nombres complexes et des opérations telles que des intégrales à portée de main. (C'est faisable à coup sûr, et un bon code C ++ peut être un ordre de grandeur plus rapide, avec de bonnes abstractions et modèles de bibliothèque, il peut même être raisonnablement propre et clair, mais il est certainement plus facile de commencer par exemple MATLAB.)

• MATLAB propose également des "boîtes à outils" pour, par exemple, le traitement d'image et le traitement numérique du signal , qui peuvent être très utiles, selon ce que vous faites.

• Le traitement numérique du signal de Mitra est un bon livre pour apprendre (dans MATLAB!) Les bases du temps discret, des filtres, des transformations, etc., ce qui est à peu près une connaissance obligatoire pour faire des algorithmes techniques décents.

Méthodes numériques. Il s'agit généralement d'un cours universitaire supérieur et d'un manuel.

Le DSP est généralement proche de l'intersection des méthodes numériques et d'une mise en œuvre efficace. Si vous ignorez l'efficacité, alors vous pourriez chercher une méthode d'approximation numérique qui pourrait produire un résultat «suffisamment précis» pour les équations du document technique qui vous intéressent. Parfois, il peut s'agir de données échantillonnées, où les théorèmes d'échantillonnage imposent des limites à la fois à la méthode d'acquisition des données (préfiltrage) et à la gamme ou à la qualité des résultats que vous pouvez obtenir avec ces données.

Parfois, Matlab, les recettes numériques ou diverses bibliothèques de traitement d'image / signal auront des algorithmes ou du code efficaces pour la solution numérique souhaitée. Mais parfois, vous devrez peut-être lancer le vôtre, il est donc utile de connaître les mathématiques derrière diverses méthodes de résolution numérique. Et c'est un gros sujet à part entière.