Szumy

Zdjęcia są zazwyczaj w pewnym stopniu zaszumione. Szumem na zdjęciu będziemy nazywać wszystko, co zniekształca / degraduje idealny obraz. Szumy mogą mieć poważny wpływ na proces przetwarzania obrazu. Aby sprawnie przetwarzać obrazy musimy być w stanie usuwać (bądź istotnie redukować) zakłócenia.

Najczęściej używaną miarą zaszumienia jest stosunek sygnału do szumu (the signal to noise ratio). Dla obrazu $f(i,j)$ stosunek sygnału do szumu definiuje się jako:

\[ \frac{S}{N}ration = \frac{\sum_{(i,j)} f^2(i,j) }{ \sum_{(i,j)} v^2(i,j) } \]

gdzie $v(i,j)$ jest szumem.

Zacznijmy od omówienie podstawowych szumów oraz nauczmy się je generować. Będzie to przydatne w dalszej części tego rozdziału, aby móc symulować zakłócenia i móc weryfikować poprawność metod usuwania szumów.

Typy szumów

Najczęściej spotykanymi szumami są:

  • szum Gaussowski,
  • szum Jednostajny,
  • szum typu pieprz i sól.

Szum Gaussowski 

Szum Gaussowski dobrze symuluje zakłócenia, które pojawiają się w rzeczywistości.

Szum Gausowski $v(i,j)$ (w punkcie $(i,j)$) jest modelowany za pomocą rozkładu Gaussa o średniej $\mu$ (która zazwyczaj jest równa 0) oraz odchyleniu standardowym $\sigma$.  

Przypomnijmy sobie wzór na jednowymiarowy rozkład normalny Rozkład normalny:

\[ N(\mu, \sigma)(x) = \frac{1}{\sigma\sqrt{2\pi}} \exp \left(-\frac{\left(x-\mu\right)^2}{2\sigma^2} \right) \]

Na poniższym rysunku przedstawiamy wykres gęstości rozkładu normalnego o różnych parametrach:

Przypomnijmy jeszcze zasadę trzech sigm, która mówi, że około $68,3%$ pola pod wykresem krzywej znajduje się w odległości jednego odchylenia standardowego od średniej, około $95,5%$ w odległości dwóch odchyleń standardowych i około $99,7%$ w odległości trzech.

Zanim przejdziemy do szumów na obrazach zacznijmy od prostego przykładu jednowymiarowego. Zacznijmy od wygenerowania sygnału (wygenerujmy falę). Będziemy używać Pakietu R.

wave = 2*sin(2*pi*(1:200)/50 + .6*pi)
plot.ts(wave, ylim=c(-3,3), main="wave")

 

 

Teraz dodajmy do niego szum typu Gausowskiego:

wave = 2*sin(2*pi*(1:200)/50 + .6*pi)
norm = rnorm(200,0,1)  
plot.ts(wave+norm, ylim=c(-3,3), main="wave")

Jak widzimy informacja o naszym sygnale została prawie całkowicie utracona. 

 

Dodajmy teraz do obrazu Leny szum typu Gaussowskiego. Obraz poniżej został wygenerowany za pomocą kodu: OpenCV (C++)

    

Ilustracja: Klasyczny obraz Leny oraz obraz z dodanym szumem: $N(0,10)$, $N(0,100)$. 

Szum jednostajny

Szum jednostajny jest rzadziej spotykany niż pozostałe dwa (Gaussowski oraz Sól i Pieprz). Niemniej jednak wspomnieć należy o nim właśnie tutaj, ponieważ różni się od szumu Gaussowskiego tylko rozkładem zmiennej losowej opisującej zakłócenia. Jak łatwo się domyślić używa się rozkładu jednostajnego na odcinku (Rozkład jednostajny). Przypomnijmy wzór na gęstość rozkładu jednostajnego na odcinku $[a,b]$:

\[ U[a,b](x) = \begin{cases} \frac{1}{b - a} & \text{for } x \in [a,b] \\ 0 & \text{otherwise} \end{cases} \]

Podobnie jak wcześniej zacznijmy od danych jednowymiarowych. Dodajmy do naszego sygnału szum typu jednostajnego:

wave = 2*sin(2*pi*(1:200)/50 + .6*pi)
unif = runif(200, min=0, max=1)
plot.ts(wave+unif, ylim=c(-3,3), main="wave + Gaussian noise")

 

Jak widzimy rozkład jednostajny powoduje, że szum, który zakłóca oryginalny sygnał jest jednostajny. Dodanie szumu typu jednostajnego do obrazu (OpenCV (C++)).

    
Ilustracja: Klasyczny obraz Leny oraz obraz z dodanym szumem: $U(0,50)$, $U(0,100)$.

Szum typu pieprz i sól

Szum typu pieprz i sól polega na dodaniu wartości ekstremalnych (np -1,1) które w przypadku obrazów w odcieniach szarości przyjmuje barwę czarną i białą. Podobnie jak w poprzednich przypadkach zacznijmy od danych jednowymiarowych:

wave = 2*sin(2*pi*(1:200)/50 + .6*pi)
saltPaper = sample(-1:1,200,rep=TRUE,prob=c(.05,.9,.05))
plot.ts(wave+saltPaper, ylim=c(-3,3), main=" wave + salt and paper noise")

Szum typu pieprz i sól polega na dodaniu czarnych i białych pikseli w przypadku zdjęć w odcieniach szarości. W przypadku zdjęć kolorowych dodajemy skrajne wartości na odpowiednich kanałach. Oczywiście z jakimś prawdopodobieństwem. Poniższe obrazy zostały uzyskane za pomocą kodu: OpenCV (C++)

  

  

Modelowanie szumów

Aby dobrze modelować szum na obrazach musimy dodać zaburzenie do oryginalnego obrazu. Możemy użyć wiele różnych sposobów aby uzyskać taki efekt. Najpopularniejsze z nich to:

  1. szum addytywny - gdy dodany szum nie zależy od obrazu:
    \[ f (i, j) = g(i, j) + v(i, j)  \] gdzie $g$ to oryginalny obraz a $v$ jest szumem.
  2. szum multiplikatywny - występuje gdy szum jest zależny od obrazu:
    \[ f (i, j) = g(i, j) + g(i, j) \cdot v(i, j) \] gdzie $g$ to oryginalny obraz a $v$ jest szumem.

 Miary zaszumienia

Aby zbadać poziom zaszumienia obrazu można użyć następujących miar:

  1. Mean Quadratic Difference  \[ \sum \sum \left( g(i,j) - f(i,j) \right)^2 \]
  2. Mean Absolute Difference  \[ \sum \sum | g(i,j) - f(i,j) | \]

Możemy również wyodrębnić sam szum \[ v(i,j) = f(i,j) - g(i,j) \] 

W naturalny sposób powyższe miary mają jedną bardo dużą wadę. A mianowicie musimy znać zarówno oryginalny obraz $g$ jak i zaszumiony $f$. W rzeczywistości jest to bardzo trudne. Dzięki temu, że umiemy generować szum mamy oba powyższe obiekty.

Filtry

Usuwanie (a raczej redukowanie) szumów to jedno z najbardziej podstawowych zagadnień analizy obrazów. Istnieje wiele metod, które umożliwiają oczyszczenie obrazów z różnego rodzaju zaburzeń. Dobór techniki zależy od rodzaju szumu oraz typu obrazów.

Najczęściej stosowane podejścia do redukcji szumów opierają się na liniowych transformacjach czyli takich, w których obliczenia mogą być wyrażone jako sumy liniowe. Usuwanie szumów za pomocą liniowych przekształceń zazwyczaj powoduje zacieranie ostrych krawędzi. 

Do filtrowania używa się również przekształceń nieliniowych czyli takie, które nie mogą być wyrażone jako prosta sumy. Pozwalają one usuwać szumy bez rozmywania (przynajmniej z mniejszym efektem rozmycia) krawędzi. 

Filtr średniej ruchomej (średniej kroczącej)

Jedną z metod usuwania szumu jest filtr średniej ruchomej (średniej kroczącej z ang. Moving average). Zacznijmy od danych jednowymiarowych. Filtr średniej ruchomej (średniej kroczącej) polega na przesuwaniu okna o określonej szerokości (na poniższym rysunku przedstawiony jest filtr średniej ruchomej z oknem o szerokości 3).  


Ilustracja. Przykład obrazujący efekt filtru średniej ruchomej z szerokością okna 3.

Używa się okien o nieparzystej szerokości. W każdym punkcie danych przykłada się okno o danej szerokości i element jest zamieniany na średnią z elementów w oknie.

Niech będzie dany punkt $x_k$ oraz jego otoczenie o szerokości $2N+1$ gdzie $N \in \mathrm{N}$ 

\[ (x_{k-N}, x_{k-N+1}, \ldots ,x_{k}, \ldots , x_{k+N-1}, x_{k+N} ) \]

 wtedy po przejściu filtru średniej ruchomej zostanie on zamieniony na: 

\[ y_{k}=mean( x_{k-N}, x_{k-N+1}, \ldots ,x_{k}, \ldots , x_{k+N-1}, x_{k+N} ) = \frac{ x_{k-N} + x_{k-N+1} + \ldots + x_{k} + \ldots + x_{k+N-1} + x_{k+N} }{2N+1}\]

Przetestujmy jak działa filtr średniej ruchomej na danych jednowymiarowych.

Szum Gaussowski

wave = 2*sin(2*pi*(1:200)/50 + .6*pi)
norm = rnorm(200,0,1)  
par(mfrow=c(3,1))
plot.ts(wave, ylim=c(-3,3), main="wave", ylab="")
ts.plot(ts(wave+norm), wave, ylim=c(-3,3), col=1:2, main="wave + Gaussian noise", ylab="")
ts.plot(filter(wave+norm,sides=2, rep(1/3,3)), wave, col=1:2,  ylim=c(-3,3), main="moving average", ylab="")

Szum jednostajny

wave = 2*sin(2*pi*(1:200)/50 + .6*pi)
unif = runif(200, min=0, max=1)
par(mfrow=c(3,1))
plot.ts(wave, ylim=c(-3,3), main="wave", ylab="")
ts.plot(ts(wave+unif), wave, ylim=c(-3,3), col=1:2, main="wave + Uniform noise")
ts.plot(filter(wave+unif,sides=2, rep(1/3,3)), wave, col=1:2,  ylim=c(-3,3), main="moving average", ylab="")

Szum typu sól i pieprz

wave = 2*sin(2*pi*(1:200)/50 + .6*pi)
saltPaper = sample(-1:1,200,rep=TRUE,prob=c(.05,.9,.05))
par(mfrow=c(3,1))
plot.ts(wave, ylim=c(-3,3), main="wave", ylab="")
ts.plot(ts(wave+saltPaper), wave, ylim=c(-3,3), col=1:2, main=" wave + salt and paper noise")
ts.plot(filter(wave+saltPaper,sides=2, rep(1/3,3)), wave, col=1:2,  ylim=c(-3,3), main="moving average", ylab="")

Jeżeli filtrujemy obraz w odcieniach szarości, to mamy do dyspozycji dane dwuwymiarowe (tablicę/ macierz). Każdy punkt zamienimy na średnią z pewnego otoczenia. Najmniejsze otoczenie punktu w takim przypadku zawiera 9 wartości.

Ilustracja. Przykład działania filtru średniej kroczącej w przypadku zdjęcia.

Jeżeli piksel jest zastępowany za pomocą średniej z danych znajdujących się w jego otoczeniu, to taki filtr średniej ruchomej nazywamy lokalnym rozmyciem (z ang. local averaging). Oczywiście zamiast rozważać średnie możemy wziąć pod uwagę jakąś inną operację.

Ilustracja. Schemat działania filtrów lokalnych.

Jeżeli zostanie użyta średnia ważona z wagami odpowiadającymi wartością gęstości dwuwymiarowego rozkładu normalnego to otrzymamy filtr Gaussowski, natomiast jeżeli będziemy rozważać medianę to dostaniemy filtr medianowy.

Efekt filtru średniej ruchomej w przypadku zdjęcia Leny prezentuje poniższa ilustracja (OpenCV (C++)):

    

    

Wynik filtru średniej ruchomej z różnymi parametrami można prześledzić za pomocą programu OpenCV (C++).

Porównajmy jeszcze jak działa taki filtr w przypadku szumów, które omówiliśmy na początku tego rozdziału (OpenCV (C++)):

  

Ilustracja. Obraz z szumem Gaussowski $N(0,100)$ oraz obraz po rozmyciu filtrem średniej ruchomej z oknem $5 \times 5$

  

Ilustracja. Obraz z szumem Jednostajnym $U(0,100)$ oraz obraz po rozmyciu filtrem średniej ruchomej z oknem $5 \times 5$

  

Ilustracja. Obraz z szumem typu pieprz i sól oraz obraz po rozmyciu filtrem średniej ruchomej z oknem $5 \times 5$

Filtr medianowy 

Filtr medianowy jest bardzo podobny do filtru średniej ruchomej. Tak naprawdę zamieniamy średnią w oknie na medianę Mediana z danych w oknie.

Niech będzie dany punkt $x_k$ oraz jego otoczenie o szerokości $2N+1$ gdzie $N \in \mathrm{N}$ 

\[ (x_{k-N}, x_{k-N+1}, \ldots ,x_{k}, \ldots , x_{k+N-1}, x_{k+N} ) \]

 wtedy po przejściu filtru medianowego zostanie on zamieniony na:

\[ y_{k}=med( x_{k-N}, x_{k-N+1}, \ldots ,x_{k}, \ldots , x_{k+N-1}, x_{k+N} ) \].

Podobnie jak w przypadku filtru średniej biegnącej rozważana jest ona o nieparzystej ilości elementów. W takim przypadku łatwo wyliczyć medianę. Musimy tylko posortować próbkę i wyciągnąć środkowy element.

Prześledźmy jak działa filtr medianowy w przypadku danych jednowymiarowych.

Szum Gaussowski

wave = 2*sin(2*pi*(1:200)/50 + .6*pi)
norm = rnorm(200,0,1)  
par(mfrow=c(3,1))
plot.ts(wave, ylim=c(-3,3), main="wave", ylab="")
ts.plot(ts(wave+norm), wave, ylim=c(-3,3), col=1:2, main="wave + Gaussian noise", ylab="")
ts.plot(ts(runmed(wave+norm, k=3)), wave, col=1:2,  ylim=c(-3,3), main="median filter", ylab="")

Szum Jednostajny

wave = 2*sin(2*pi*(1:200)/50 + .6*pi)
unif = runif(200, min=0, max=1)
par(mfrow=c(3,1))
plot.ts(wave, ylim=c(-3,3), main="wave", ylab="")
ts.plot(ts(wave+unif), wave, ylim=c(-3,3), col=1:2, main="wave + Uniform noise")
ts.plot(ts(runmed(wave+unif, k=3)), wave, col=1:2,  ylim=c(-3,3), main="median filter", ylab="")

 

Szum typu pieprz i sól

wave = 2*sin(2*pi*(1:200)/50 + .6*pi)
saltPaper = sample(-1:1,200,rep=TRUE,prob=c(.05,.9,.05))
par(mfrow=c(3,1))
plot.ts(wave, ylim=c(-3,3), main="wave", ylab="")
ts.plot(ts(wave+saltPaper), wave, ylim=c(-3,3), col=1:2, main=" wave + salt and paper noise")
ts.plot(ts(runmed(wave+saltPaper, k=3)), wave, col=1:2,  ylim=c(-3,3), main="median filter", ylab="")

 

Jak widzimy na powyższych rysunkach, filtr medianowy dobrze sobie radzi w przypadku szumu typu pieprz i sól jak i szumów Gaussowskich. 

Podobnie jak w przypadku średniej kroczącej gdy mamy do czynienia z obrazami (tablicami / macierzami dwuwymiarowymi), będziemy rozważać dwuwymiarowe otoczenia punktu. 

Przykład:

Zauważmy, że otoczenie punktu $3 \times 3$ składa się z liczb 25, 21, 23, 25, 18, 255, 30, 13, 22. Po posortowaniu otrzymujemy 13, 18, 22, 21, 23, 25, 25, 30, 255.
Mediana to środkowy element 23. Natomiast średnia wynosi 48.

Efekt filtru medianowego w przypadku zdjęcia Leny prezentuje poniższa ilustracja (OpenCV (C++)):

   

   

Jak łatwo zauważyć filtr medianowy w dużo mniejszym stopniu rozmywa krawędzie. Jednym z powodów jest fakt, że kolor każdego piksela jest zamieniany na inny kolor z obrazu (medianę z danych).

Wynik filtru medianowego z różnymi parametrami można prześledzić za pomocą programu OpenCV (C++).

Porównajmy jeszcze jak działa taki filtr w przypadku szumów, które omówilismy na początku tego rozdziału (OpenCV (C++)):

  

Ilustracja. Obraz z szumem Gaussowski $N(0,100)$ oraz obraz po rozmyciu filtrem medianowym z oknem $5 \times 5$

  

Ilustracja. Obraz z szumem Jednostajnym $U(0,100)$ oraz obraz po rozmyciu filtrem medianowym z oknem $5 \times 5$

  

Ilustracja. Obraz z szumem typu pieprz i sól oraz obraz po rozmyciu filtrem medianowym z oknem $5 \times 5$

Jak widzimy filtr medianowy działa bardzo dobrze. Ponadto uzyskujemy mniejszy stopień rozmycia krawędzi. Należy jednak zauważyć, że filtr medianowy usuwa cienkie linie oraz rozmywa rogi. Problem ten można rozwiązać używając bardziej skomplikowanych masek.

Drugim problemem związanym z filtrem medianowym jest jego złożoność obliczeniowa. W standardowej wersji złożoność wynosi $O(k^2 \log(k))$. 

W pracy S. Perreault, P. Hebert, Median Filtering in Constant Time, IEEE Transactions on Image Processing, 2389--2394, 2007 zaproponowany został algorytm, który redukuje tą procedurę do $O(1)$ (w czasie stałym). 

Aby zrozumieć ideę, która kryje się za powyższą metodą trzeba zdać sobie sprawę z dwóch bardzo ważnych faktów. 

Pierwszy mówi, że jeżeli mamy dwa histogramy zbudowane na dwóch rozłącznych zbiorach $A$, $B$ $(A \cap B = \emptyset)$, to wspólny histogram możemy uzyskać w czasie stałym (w stosunku do ilości dostępnych wartości, w przypadku obrazów [0,255]).

\[ H( A \cup B ) = H(A) + H(B) \]

Powyższą operację można rozumieć jako dodawanie dwóch wektorów po współrzędnych.

Drugi fakt to zwykłe spostrzeżenie, że z histogramu możemy łatwo uzyskać medianę.

Powyższy algorytm dzieli okno (które będzie się poruszać po obrazie na kolumny i w każdej z nich liczy histogram). Podczas poruszania się maski mogą nastąpić dwie sytuacje:

  • jeżeli maska porusza się poziomo (prawy obrazek na poniższej ilustracji) - musimy wtedy usunąć jedną kolumnę (jeden histogram) z lewej strony i dodać jedną kolumnę z prawej,
  • jeżeli maska porusza się w pionie (prawy obrazek na poniższej ilustracji) - wtedy w każdym z histogramów musimy usunąć jeden element znajdujący się na samej górze i dodać jedną wartość piksela z dołu.

Obie powyższe operacje są proste (z numerycznego punktu widzenia).

Pseudokod algorytmu:

Input: Image $X$ of size $m ,n$, kernel radius $r$
Output: Image $Y$ of the same size
Initialize kernel histogram $H$ and column histograms $h_{1...n}$
for $i = 1$ to $m$ do
  for $j = 1$ to $n$ do
    Remove $X_{i-r-1,j+r}$ from $h_{j+r}$
    Add $X_{i+r, j+r}$  to $h_{j+r}$
    $H \leftarrow  H + h_{j+r} - h_{j-r-1}$
    $Y_{i,j} \leftarrow  median(H)$
  end for
end for

Rozmycie Gaussa

Trzecią i ostatnią metodą, którą omówimy w kontekście usuwania szumów jest rozmycie Gaussa. Na początku tego rozdziału wprowadziliśmy wzór na gęstość jednowymiarowego rozkładu normalnego. Teraz będzie nam potrzebna wersja dwuwymiarowa:

\[ N( \mu, \Sigma)(x)= \frac {1}{(2\pi) \det(\Sigma)^{1/2}} \exp\left( -\frac{1}{2} (X - \mu)^T \Sigma^{-1} (X -\mu)\right) \],

gdzie $\mu = [\mu_1,\mu_2]^T$ oraz $\Sigma = \begin{bmatrix} \sigma_{11} & \sigma_{12} \\ \sigma_{21} & \sigma_{22} \end{bmatrix}$.

Wykresem funkcji gęstości rozkładu normalnego (2D) jest kapelusz  Gaussa.

Ilustracja. Gęstość rozkładu normalnego z różnymi parametrami.

W rozmyciu Gaussowskim, podobnie jak w filtrze średniej kroczącej będziemy liczyć średnią z wartości kolorów pikseli w otoczeniu danego punktu. Tym razem użyjemy średniej ważonej. Jako wagi będziemy używać wartości rozkładu normalnego w punkcie. W takim przypadku punkt z obrazu o współrzędnych $(i,j)$ (o kolorze $f(i,j)$) zmieni odcień na $f_{new}(x,y)$:

\[ f_{new}(x,y) = \sum_{i=-r}^{r} \sum_{j=-r}^{r} f(x+i,y+j) N(\mu,\Sigma)( i, j )\].

  

Ilustracja. Wyznaczanie wag w filtrze Gaussowskim.

Jak widzimy w przypadku filtru Gaussowskiego musimy podać dwa parametry szerokość okna oraz $\Sigma$, ponieważ najczęściej używa się sferycznych rozkłdów normalnych czyli takich, że:  

 \[ \Sigma = \begin{bmatrix} \sigma & 0 \\ 0 & \sigma \end{bmatrix} \]

oraz średniej równej $\mu = [0,0]^T$ więc nasza gęstość pobiera tylko jeden parametr $\sigma$ odnośnie rozkładu normalnego. 

Efekt filtru medianowego w przypadku zdjęcia Leny prezentuje poniższa ilustracja OpenCV (C++)

   

   

Wynik filtru Gaussowskiego z różnymi parametrami można prześledzić za pomocą programu OpenCV (C++).

Porównajmy jeszcze jak działa taki filtr w przypadku szumów, które omówiliśmy na początku tego rozdziału (OpenCV (C++)):

  

Ilustracja. Obraz z szumem Gaussowski $N(0,100)$ oraz obraz po rozmyciu filtrem Gaussowskim z oknem $5 \times 5$ oraz $\sigma = 1,5$.

  

Ilustracja. Obraz z szumem Jednostajnym $U(0,100)$ oraz obraz po rozmyciu filtrem Gaussowskim z oknem $5 \times 5$ oraz $\sigma = 1,5$.

  

Ilustracja. Obraz z szumem typu pieprz i sól oraz obraz po rozmyciu filtrem Gaussowskim z oknem $5 \times 5$ oraz $\sigma = 1,5$.

Na zakończenie tej części wspomnijmy jeszcze, że można uogólnić nasz filtr Gaussowski dla dowolnych wag. Możemy używać otoczenia oraz jądra lub maski - czyli macierzy odpowiadającej naszemu otoczeniu, które definiuje wagi. Dokonując operacji brania średniej ważonej tak naprawdę wykonujemy operację splotu (będziemy jeszcze o tym mówić - tutaj tylko zaznaczamy że coś takiego istnieje).

Możemy na przykład użyć maski:

\[ \begin{bmatrix} 0 & -1 & 0 \\ -1 & 5 & -1 \\ 0 & -1 & 0 \end{bmatrix} \] 

W takim przypadku otrzymamy efekt wyostrzenia (dokładnie operacja ta zostanie opisana w późniejszych rozdziałach)

  

Powyższy efekt można uzyskać za pomocą programów OpenCV (C++).