TI/Programowanie dla Fizyków Medycznych/Manipulacja obrazem

Maniplacja obrazem

Ćwiczenia z mainuplacji obrazem rozpocznijmy od wczytania pliku który będziemy przetwarzać. Dla fizyków medycznych naturalnie będzie to plik DICOM.

import dicom
import numpy as np
import pylab as py

plik=dicom.read_file('I00001.dcm')
pixel=plik.pixel_array
print pixel.shape
>>>(2964, 2364)

Stworzona w ten sposób tablica ma dosyć duże rozmiary. O ile wyświetlenie takiego pliku nie jest problemem dla znajdujących się w pracowni komputerów, to bardziej złożona modyfikacja takiej grafiki mogła by z czasem obliczeń znacznie wykraczać poza czas przewidziany na zajęcia. Wyłącznie dla celów dydaktycznych, by ułatwić i przyspieszyć pracę zmniejszymy rozmiar badanego obrazu.

pixel=pixel[::5,::5]
print min(pixel.flatten()),max(pixel.flatten())
>>> 0 3907

Wartości zapisane w tablicy odpowiadającej zmniejszonemu obrazowi mają wartości w przedziale [0,3907]. Metoda wyświetlając imshow najwyższej wartości czyli 3908 przypisze kolor biały, wartości 0 kolor czarny. Aby ułatwić sobie manipulacje kolorami przypiszmy im wartości z zakresu [0,1].

pixel=pixel*1.0/max(pixel.flatten())

Tak powstały obraz możemy dowolnie modyfikować. Aby lepiej widoczne były skutki manipulacji obrazem dodajmy u góry obrazka þlynne przejście od czerni do bieli

def dodajPasek(tablica, n):
    tablica[:n,:]=np.linspace(0,1,tablica.shape[1])
dodajPasek(pixel,30)
py.imshow(pixel,cmap=py.cm.gray,interpolation='nearest')
py.show()

Możemy teraz zdefiniować dowolną funkcję przyjmującą za arguemnty wartości z przedziału [0,1] i zwracającą wartości z tego samego przedziału. Działając taką funkcją na każdy pixel obrazu dokonamy jego transformacji. Najprostrze funkcje jakie mogą przyjść od razu do głowy posiadają już swoje tradycyjne nazwy związane z wpływem jaki mają na obraz.

Jasność

Najprostszą funkcję od której na Wydziale Fizyki zaczynamy zawsze jest funkcja liniowa. Musimy jedynie zapewnić aby przy wspólczynniku nachylenia większym od 1.0 wynikowe wartości nie przekraczały jedności. Zdefiniujmy zatem jednoparametrową funkcje medyfikującą obraz.

def jasnosc(tablica, a):
    def f(x,a):
        return min(a*x,1.0)
    wynik=tablica.copy()
    for x,y in np.ndindex(tablica.shape):
        wynik[x,y]=f(tablica[x,y],a)
    dodajPasek(wynik,15)
    return wynik

Aby móc porównać wynik w pierwotnym obrazem znów wstawiamy pasek szarości ale już o połowę węższy. dzieki temu możemy porównać pasek po operacji przekształcenia i przed nią. Dla a>1 obraz powinien się rozjaśnić, natomiast dla a<1 powinien być ciemniejszy. Zaobaczmy teraz na dwóch przykładach w jaki sposób funkcja jasnosc modyfikuje obraz

py.imshow(jasnosc(pixel,2),cmap=py.cm.gray,interpolation='nearest')
py.show()

py.imshow(jasnosc(pixel,0.5),cmap=py.cm.gray,interpolation='nearest')
py.show()

Gamma

Kolejną narzucającą się funkcją jest podnoszenie wartości zapisanej w tablicy do potęgi. Dla wykładników większych od jeden jasne kolory staną się jeszcze jaśniejsze a przejście między bardzo podobnymi ciemniejszymi odcieniami szarości staną się bardziej widoczne. Dla wykładników mniejszych od jeden, ciemne kolory staną się jeszcze ciemniejsze, natomiast bardzo biskie siebie jasne odcienie staną sie lepiej rozróżnialne. Zdefiniujmy zatem funkcję gamma.

def gamma(tablica, a):
    def f(x,a):
        return min(x**a,1.0)
    wynik=tablica.copy()
    for x,y in np.ndindex(tablica.shape):
        wynik[x,y]=f(tablica[x,y],a)
    dodajPasek(wynik,15)
    return wynik

Możemy teraz na dwóch przykładach obejrzeć efekt działania przekształcenia gamma.

        
py.imshow(gamma(pixel,2),cmap=py.cm.gray,interpolation='nearest')
py.show()

py.imshow(gamma(pixel,0.5),cmap=py.cm.gray,interpolation='nearest')
py.show()

Próg

def prog(tablica, a):
    def f(x,a):
        return 0 if x<a else 1
    wynik=tablica.copy()
    for x,y in np.ndindex(tablica.shape):
        wynik[x,y]=f(tablica[x,y],a)
    dodajPasek(wynik,15)
    return wynik

py.imshow(prog(pixel,0.5),cmap=py.cm.gray,interpolation='nearest')
py.show()

Schodek

def schodek(tablica, a,b):
    def f(x,a,b):
        return 1 if ((x<a) or (x>b)) else 0
    wynik=tablica.copy()
    for x,y in np.ndindex(tablica.shape):
        wynik[x,y]=f(tablica[x,y],a,b)
    dodajPasek(wynik,15)
    return wynik

py.imshow(schodek(pixel,0.2,0.6),cmap=py.cm.gray,interpolation='nearest')
py.show()

Rozmycie i wyostrzenie

from scipy import ndimage
imag=ndimage.gaussian_filter(pixel,sigma=4)
print imag

WYOSTRZ

def wyostrz(tablica, sigma, a):
    rozmyty=np.array(ndimage.gaussian_filter(tablica,sigma=sigma))
    roznica=tablica-rozmyty
    wynik=rozmyty+a*roznica
    wynik/=max(wynik.flatten())
    dodajPasek(wynik,15)
    return wynik

py.imshow(wyostrz(pixel,4,0),cmap=py.cm.gray,interpolation='nearest')
py.show()

py.imshow(wyostrz(pixel,4,3),cmap=py.cm.gray,interpolation='nearest')
py.show()

"Programowanie dla Fizyków Medycznych"