Autor Wątek: Dekonwolucja 2D  (Przeczytany 3421 razy)

Offline komorra

  • Użytkownik
    • Blog naszego teamu (o grze Voxelfield)

# Maj 26, 2015, 11:03:56
Czy jest jakiś sposób aby przeprowadzić dekonwolucję 2D na tablicy float[,] w sposób nie wymagający użycia np. FFT? Chciałbym to osiągnąć używając jedynie pętli i prostych działań arytmetycznych. Wiem, że dla 1D da się bez problemu. Dodam, że kernel jest mi znany.

Offline Mr. Spam

  • Miłośnik przetworów mięsnych

Offline Krzysiek K.

  • Redaktor
    • DevKK.net

# Maj 26, 2015, 14:13:49
Cytuj
Czy jest jakiś sposób aby przeprowadzić dekonwolucję 2D na tablicy float[,] w sposób nie wymagający użycia np. FFT? Chciałbym to osiągnąć używając jedynie pętli i prostych działań arytmetycznych.
FTT to przecież też są wyłącznie pętle i proste działania matematyczne. :)

Cytuj
Dodam, że kernel jest mi znany.
W takim razie potraktuj całość jako układ równań liniowych aM = b ( a - dane sprzed splotu, M - macierz wykonująca splot, b - wynik splotu ) i rozwiąż. Sprowadzi się to do odwrócenia macierzy M.

Offline komorra

  • Użytkownik
    • Blog naszego teamu (o grze Voxelfield)

# Maj 26, 2015, 14:43:23
Szybka (i może głupia) wypowiedź/pytanie: czyli mogę przyjąć, że splot to to samo co mnożenie macierzowe? (zaraz to sprawdzę ale wydaje mi się coś nie tak)

Offline komorra

  • Użytkownik
    • Blog naszego teamu (o grze Voxelfield)

# Maj 26, 2015, 15:08:58
No wygląda na to, że to dwa różne pojęcia i nie mogę tutaj użyć odwrotności albo mnożenia macierzowego.

https://en.wikipedia.org/wiki/Convolution

https://en.wikipedia.org/wiki/Matrix_multiplication

Może pójdę w stronę FFT, ale chciałbym bardziej mieć coś nisko-kosztowego jeżeli chodzi o cepeużerność.

Offline CheshireCat

  • Użytkownik

# Maj 26, 2015, 15:11:12
Może pójdę w stronę FFT, ale chciałbym bardziej mieć coś nisko-kosztowego jeżeli chodzi o cepeużerność.

To może zrobić to na GPU, jeżeli platforma pozwala?
Ludzie od renderowania oceanów robią często FFT na GPU właśnie.

Offline komorra

  • Użytkownik
    • Blog naszego teamu (o grze Voxelfield)

# Maj 26, 2015, 15:12:34
Odpada GPU - platforma targetowa to ios, android, wp, ogólnie mobilne i to powinno chodzić dobrze na starych modelach (z naciskiem na nie właściwie).

Offline Krzysiek K.

  • Redaktor
    • DevKK.net

# Maj 26, 2015, 15:35:50
Cytuj
No wygląda na to, że to dwa różne pojęcia i nie mogę tutaj użyć odwrotności albo mnożenia macierzowego.
To dwa różne pojęcia, ale oba dotyczą operacji liniowych. A że masz sygnał dyskretny, a kernel jest stały, to splot z nim staje się operacją liniową na skończonej liczbie wartości - więc można przedstawić to jako operację macierzową.

Przykładowo masz sygnał [0 1 2 3 4] i kernel [1 1]. Splot można zapisać wtedy macierzowo jako:
              [1 1 0 0 0]
              [0 1 1 0 0]
[0 1 2 3 4] * [0 0 1 1 0]
              [0 0 0 1 1]
              [1 0 0 0 1]

I wychodzi z tego coś w wyniku. Jeśli teraz bym miał sam wynik, a nie miał sygnału wejściowego, to jak najbardziej mogę go wyliczyć, bo macierz robiąca splot jest odwracalna.

W przypadku sygnału 2D o wymiarach NxM macierz będzie miała N*M wierszy oraz kolumn (jedna kolumna wyznacza pojedynczy piksel), więc robi się tego bardzo dużo i dla obrazów sensownych wielkości trzeba już pokombinować trochę teoretycznie, ale jest to zawsze jakiś punkt wyjścia.

Offline mkk12259

  • Użytkownik

# Maj 26, 2015, 17:47:35
Raz: Jeżeli dobrze rozumiem propozycję wyżej, to dostaniesz układ równań o wymiarze N, gdzie N to liczba pikseli. Przy czym macierz A, w równaniu Ax=b, jest rzadka. Do rozwiązywania układów równań z macierzami rzadkimi masz mnóstwo metod: iteracyjne (sparse iterative solvers), albo przez bezpośrednie odwracanie macierzy rzadkich (sparse direct solvers - superLU, UMFPACK, PARDISO, etc). Te drugie są stabilne i szybkie ale wymagają raczej dużo pamięci operacyjnej. Więc dla mobilnych platform zostają metody iteracyjne.

Dwa: trzeba uważać bo nie wszystkie rodzaje konwolucji są odwracalne :)

Offline komorra

  • Użytkownik
    • Blog naszego teamu (o grze Voxelfield)

# Maj 27, 2015, 20:50:17
Dobra, poszedłem w stronę FFT ale ni kij nie wiem co jest nie tak, próbuję próbuję i nic.



Lewy, górny - obraz wejściowy
Prawy, górny - kernel (suma pól = 1)
Lewy, dolny - konwolucja wejścia i kernela
Prawy, dolny - niby ? dekonwolucja (a powinien być identyczny z lewym górnym).

Kod leci mniej więcej tak:

Kod: (csharp) [Zaznacz]
            FFT fker = new FFT(ckernel); //obiekt transformaty kernela 16x16
            ckernel = fker.FFT2D(ckernel, 16, 16, 1); //transformata FFT w przód
            ckernel = fker.FFT2D(ckernel, 16, 16, 1);//transformata FFT w przód powtórnie

            FFT tcon = new FFT(cconv); //obiekt transformaty konwolucji 16x16
            cconv = tcon.FFT2D(cconv, 16, 16, 1); //transformata FFT w przód jeden raz

            //Poniżej stosuję wzór FFT(konwolucja)/FFT(FFT(kernel))
            var result = new COMPLEX[conv.W(),conv.H()];
            for (int i = 0; i < conv.W(); i++)
            {
                for (int j = 0; j < conv.H(); j++)
                {
                    result[i, j] = new COMPLEX(cconv[i,j].real/ckernel[i,j].real, 0);           
                }
            }

            //Następnie odwracam transformatę dla tablicy z powyższymi ilorazami
            var rf = new FFT(result);
            result = rf.FFT2D(result, 16, 16, -1);           

            //Kopiuję wszystko do tablicy floatów
            var res = new float[conv.W(), conv.H()];
            for (int i = 0; i < res.W(); i++)
            {
                for (int j = 0; j < res.H(); j++)
                {
                    res[i, j] = (float) result[i, j].real;
                }
            }

            res.Balance(); //Na koniec balans - czyli taka "normalizacja" rozciągająca od 0 do 1 w skali
            control2DDeconv.Data = res; //podstawienie tablicy do kontrolki prawej, dolnej

Bo zgodnie z http://research.stowers-institute.org/efg/Report/FourierAnalysis.pdf (slajd około 45)

powinienem podzielić FFT konwolucji (wyniku) przez funkcję FFT(PSF) (a więc FFT od FFT kernela) i wynik przetransformować inverse FFT.

Rozumiem, że taką operację mam wykonać dla każdej komórki tablicy (ten iloraz), a końcową tablicę przetransformować inverse FFT i powinno dać obraz wejściowy. (A jak widać nie dzieje się tak)

Będę wdzięczny za jakieś wskazówki, jakby ktoś miał czas i ochotę :)

Offline smajler

  • Użytkownik

# Maj 28, 2015, 09:50:12
Kiedys chcialem zrobic dekonwolucje ale bez znajomosci kernela, suwakiem sie ja budowalo metoda prob i bledow. Uzywalem wtedy do tego ImageMagick. Kupa bibliotek do majstrowania obrazem z linii komend jak i zrodla C, C++ to wbudowania w programik. Jest tam gotowy przyklad konwolucji i dekonwolucji ze skryptem. O dokumentacje ciezko, bardzo ciezko mi wtedy pomagal jakis obcokrajowiec i nawet konwolucja i dekonwolucja w C++ wyszla.
Tutaj masz przyklady i skrypt jakby to dzialalo z linii polecen. Jak odpowiada to daj znac poszukam kod C++ chyba ze nie znajde :)
http://www.imagemagick.org/Usage/fourier/fft_math/

Offline komorra

  • Użytkownik
    • Blog naszego teamu (o grze Voxelfield)

# Maj 28, 2015, 10:05:28
@smajler - dzięki :) Znalazłem też wrapper na ImageMagick-a: https://magick.codeplex.com/
W domu zobacze co uda się zdziałać.

Offline mkk12259

  • Użytkownik

  • +1
# Maj 28, 2015, 10:07:50
Zobacz sobie jeszcze na to: http://cimg.sourceforge.net/
Biblioteka składa się z jednego pliku nagłówkowego, i ma m.in. FFT. Do każdej funkcji jest przykład.

Offline komorra

  • Użytkownik
    • Blog naszego teamu (o grze Voxelfield)

# Maj 28, 2015, 10:28:31
O to też się nada, dzięki mkk12259 oraz wszystkim za pomoc :) Myślę, że powinienem teraz dać radę :)

Offline smajler

  • Użytkownik

# Maj 28, 2015, 20:25:56
Wraperry masz na glownej stronie imagemagick. Z tego co pamietam IM.NET byl bardzo ubogi i konwolucji tam nie bylo. Sciagnij sobie pod C++ i jak musisz miec w NET to wykorzystaj DLL napisana w C++, ja tak robilem.