Implementacja technicznych algorytmów papierowych w C ++ lub MATLAB

Jestem studentem inżynierii elektrycznej. Czytałem wiele artykułów technicznych na temat algorytmów przetwarzania sygnałów i obrazów (rekonstrukcja, segmentacja, filtrowanie itp.). Większość algorytmów przedstawionych w tych dokumentach jest definiowana w czasie ciągłym i częstotliwości ciągłej i często daje rozwiązania w postaci skomplikowanych równań. Jak zaimplementowałbyś od podstaw dokument techniczny w C ++ lub MATLAB, aby powtórzyć wyniki uzyskane we wspomnianym artykule?

Mówiąc dokładniej, patrzyłem na artykuł „Ogólny algorytm rekonstrukcji wiązki stożkowej” autorstwa Wanga i wsp. ( IEEE Trans Med Imaging. 1993; 12 (3): 486-96 ) i zastanawiałem się, jak mam zacząć wdrażanie ich algorytmu? Równanie 10 podaje formułę zrekonstruowanego obrazu przy. Jak byś to kodował? Czy miałbyś pętlę for przechodzącą przez każdy woksel i obliczającą odpowiednią formułę? Jak kodowałbyś funkcje funkcji w tej formule? Jak oceniasz funkcje w dowolnych punktach?

Przeczytałem książkę „Digital Image Processing” Gonzalez i Woods, ale wciąż jestem zagubiony. Czytałem również o serii książek z numerycznymi przepisami. Czy to byłby właściwy sposób?

Jakie są twoje doświadczenia w programowaniu algorytmów z prac naukowych? Wszelkie wskazówki lub sugestie?

Algorytmy przetwarzania sygnałów zdefiniowane w ciągłym czasie / przestrzeni / częstotliwości są zwykle implementowane przez próbkowanie sygnału na dyskretnej siatce i przekształcanie całek w sumy (i pochodne w różnice). Filtry przestrzenne są implementowane poprzez splot z jądrem splotu (tj. Ważoną sumą sąsiadów).

Istnieje ogromna wiedza na temat filtrowania próbkowanych sygnałów w dziedzinie czasu; filtry w dziedzinie czasu są implementowane jako filtry o skończonej odpowiedzi impulsowej , w których bieżąca próbka wyjściowa jest obliczana jako ważona suma poprzednich N próbek wejściowych; lub nieskończone filtry odpowiedzi impulsowej, w których prąd wyjściowy jest ważoną sumą poprzednich sygnałów wejściowych i poprzednich sygnałów wyjściowych . Formalnie, filtry czasu dyskretnego są opisywane przy użyciu transformaty z , która jest analogowym czasem dyskretnym do transformaty Laplace'a . Metoda tustina odwzorowuje jeden do drugiego ( c2di d2cw Matlab).

Jak oceniasz funkcje w dowolnych punktach?

Gdy potrzebujesz wartości sygnału w punkcie, który nie leży bezpośrednio na siatce próbkowania, interpolujesz jego wartość z pobliskich punktów. Interpolacja może być tak prosta, jak wybranie najbliższej próbki, obliczenie średniej ważonej najbliższych próbek lub dopasowanie arbitralnie skomplikowanej funkcji analitycznej do próbkowanych danych i ocena tej funkcji przy wymaganych współrzędnych. Interpolacja na jednolitą drobniejszą siatkę to próbkowanie w górę . Jeśli twój oryginalny (ciągły) sygnał nie zawiera szczegółów (tj. Częstotliwości) drobniejszych niż połowa siatki próbkowania, to funkcję ciągłą można doskonale zrekonstruować na podstawie próbkowanej wersji ( twierdzenie Nyquista-Shannona o próbkowaniu ). Aby zobaczyć przykład interpolacji w 2D, zobaczinterpolacja dwuliniowa .

W Matlab możesz używać interp1lub interp2interpolować dane 2D lub regularnie próbkowane dane 2D (odpowiednio) lub griddatainterpolować dane 2D nieregularnie próbkowane.

Czy miałbyś pętlę for przechodzącą przez każdy woksel i obliczającą odpowiednią formułę?

Tak, dokładnie.

Matlab ratuje Cię przed koniecznością robienia tego za pomocą jawnych pętli for, ponieważ jest zaprojektowany do pracy na macierzach i wektorach (tj. Tablicach wielowymiarowych). W Matlabie nazywa się to „wektoryzacją”. Całka mogą być przybliżone sum, cumsum, trapz, cumtrapz, itd.

Przeczytałem książkę „Digital Image Processing” Gonzalez i Woods, ale wciąż jestem zagubiony. Czytałem również o serii książek z numerycznymi przepisami. Czy to byłby właściwy sposób?

Tak, przepisy numeryczne byłyby świetnym początkiem. Jest to bardzo praktyczne i obejmuje większość metod numerycznych, których będziesz potrzebować. (Przekonasz się, że Matlab już implementuje wszystko, czego potrzebujesz, ale Przepisy numeryczne zapewnią doskonałe tło).

Wziąłem klasę „algorytmy i struktury danych”, ale nie widzę związku między prezentowanym tam materiałem a wdrażaniem algorytmów naukowych.

Materiał omawiany na kursach „Algorytmy i struktury danych” ma tendencję do koncentrowania się na strukturach takich jak listy, tablice, drzewa i wykresy zawierające liczby całkowite lub łańcuchy oraz operacje, takie jak sortowanie i wybieranie: problemy, dla których zwykle występuje jeden poprawny wynik. Jeśli chodzi o algorytmy naukowe, to tylko połowa historii. Druga połowa dotyczy metod szacowania liczb rzeczywistych i funkcji analitycznych. Znajdziesz to na kursie „Metody numeryczne” (lub „Analiza numeryczna”; taki jak ten- przewiń w dół po slajdy): jak oszacować funkcje specjalne, jak oszacować całki i pochodne itp. Tutaj jednym z głównych zadań jest oszacowanie dokładności wyniku, a jednym wspólnym wzorem jest iteracja rutyny, która poprawia oszacuj, aż będzie wystarczająco dokładny. (Możesz zadać sobie pytanie, skąd Matlab wie, jak zrobić coś tak prostego, jak oszacowanie wartości sin(x)dla niektórych x).

Jako prosty przykład, oto krótki skrypt, który oblicza transformatę radonową obrazu w Matlabie. Transformacja radonu przenosi projekcje obrazu na zbiór kątów projekcji. Zamiast próbować obliczyć projekcję pod dowolnym kątem, zamiast tego obracam cały obraz za pomocą imrotate, aby projekcja była zawsze pionowa. Następnie możemy po prostu użyć rzutu sum, ponieważ summacierz zwraca wektor zawierający sumę z każdej kolumny.

Możesz napisać własny, imrotatejeśli wolisz, używając 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]');

To, co kiedyś było całką gęstości wzdłuż promienia, jest teraz sumą nad kolumną dyskretnie próbkowanego obrazu, który z kolei został znaleziony przez interpolację oryginalnego obrazu nad przekształconym układem współrzędnych.

Dodając do doskonałego wyjaśnienia nibota , jeszcze kilka punktów.

• Środowiska komputerowe numeryczne, takie jak MATLAB, Octave lub SciPy / NumPy, pozwolą Ci zaoszczędzić dużo wysiłku w porównaniu do robienia tego samemu w ogólnym języku programowania, takim jak C ++. Żonglowanie doubletablicami i pętlami po prostu nie jest porównywalne z posiadaniem typów danych, takich jak liczby zespolone i operacje, takie jak całki na wyciągnięcie ręki. (Z pewnością jest to wykonalne, a dobry kod C ++ może być o rząd wielkości szybszy, z dobrymi abstrakcjami bibliotek i szablonami może być nawet dość czysty i przejrzysty, ale zdecydowanie łatwiej jest zacząć np. Od MATLAB.)

• MATLAB ma również „zestawy narzędzi” do np. Przetwarzania obrazu i cyfrowego przetwarzania sygnałów , które mogą bardzo pomóc, w zależności od tego, co robisz.

• Cyfrowe przetwarzanie sygnałów Mitry jest dobrą książką do nauki (w MATLAB!) Podstaw dyskretnego czasu, filtrów, transformacji itp., Która jest prawie obowiązkową wiedzą do wykonania wszelkich przyzwoitych algorytmów technicznych.

Metody numeryczne. Zwykle jest to kurs uniwersytecki wyższej klasy i podręcznik.

DSP zwykle znajduje się w pobliżu skrzyżowania metod numerycznych i wydajnej implementacji. Jeśli zignorujesz efektywność, to możesz szukać jakiejkolwiek metody przybliżenia numerycznego, która może dać „wystarczająco dokładny” wynik dla równania (ów) technicznych. Czasami można mieć do czynienia z próbkowanymi danymi, gdzie twierdzenia o próbkowaniu narzucą pewne ograniczenia zarówno na metodę akwizycji danych (wstępne filtrowanie), jak i na zakres lub jakość wyników, jakie można uzyskać, biorąc pod uwagę te dane.

Czasami Matlab, receptury numeryczne lub różne biblioteki przetwarzania obrazu / sygnału będą miały wydajne algorytmy lub kod dla pożądanego rozwiązania numerycznego. Ale czasami może być konieczne samodzielne wykonanie, więc pomaga poznać matematykę związaną z różnymi metodami rozwiązywania problemów numerycznych. I to sam w sobie duży temat.