Obrazowanie:Obrazowanie Medyczne/ImageJ cwiczenia: Różnice pomiędzy wersjami
(Nie pokazano 20 pośrednich wersji utworzonych przez tego samego użytkownika) | |||
Linia 1: | Linia 1: | ||
+ | powrót do strony głównej wykładu [[Obrazowanie Medyczne]] | ||
+ | |||
=Ćwiczenia z Obrazowania medycznego z wykorzystaniem programu ImageJ= | =Ćwiczenia z Obrazowania medycznego z wykorzystaniem programu ImageJ= | ||
Linia 7: | Linia 9: | ||
Wykonaj fantom do analizy rozdzielczości. Na przykład taki: | Wykonaj fantom do analizy rozdzielczości. Na przykład taki: | ||
[[Plik:Fantom1.png]] | [[Plik:Fantom1.png]] | ||
− | W tym celu utwórz nowy obraz: File->New->Image i podaj jego rozmiary 128 x 128 punktów. | + | |
+ | W tym celu | ||
+ | # utwórz nowy obraz: File->New->Image i podaj jego rozmiary 128 x 128 punktów. | ||
+ | # wejdź do okna Process->Math->Macro i wpisz w linijce: '''v=128*sin(2*PI/32*x)+128'''. Jak widzisz takie funkcje jak sinus lub taka wielkość jak PI są w programie predefiniowane.(np. potęgę trzecią osiąga się funkcją pow(x,3), a pierwiastek sqrt(x) itp.) | ||
+ | # zapisz obraz we własnym katalogu | ||
+ | |||
+ | Zastanów się jak | ||
+ | # zmienić częstość przestrzenną fantomu. Stwórz i zapisz dwa inne fantomy o większej częstości. | ||
+ | # utworzyć fantom o symetrii kołowej, na przykład taki: [[Plik:Fantom2.png]] | ||
+ | # utworzyć fantom o zmiennej częstości, na przykład taki:[[Plik:Fantom3.png]] | ||
+ | |||
+ | Aby utworzyć fantom zawierający funkcję schodkową (a nie sinusoidalną) proponuję wykonać własne makro. | ||
+ | # wybierz Plugins->New->Macro | ||
+ | # skopiuj poniższy skrypt | ||
+ | # uruchom makro w okienku edytora: Macros->Run macro (ctrl+R) | ||
+ | |||
+ | <source lang='c'> | ||
+ | //to jest komentarz. Poniżej znajduje się zawartość makra tworzacego fantom z funkcją schodkową | ||
+ | imax=512; | ||
+ | jmax=512; | ||
+ | T=32; | ||
+ | newImage ("rozdzielczosc", "8-bit black", imax, jmax, 1); | ||
+ | for (i=0; i<imax; i++) | ||
+ | for (j=0; j<jmax; j++) { | ||
+ | if (i%T<T/2) | ||
+ | setPixel(i, j, 255); | ||
+ | else | ||
+ | setPixel(i, j, 0); | ||
+ | } | ||
+ | </source> | ||
+ | |||
+ | Jeśli chcesz obejrzeć profil jakiejś części fantomu, możesz zaznaczyć obszar i wybrać Analyze->Plot profile (ctr+K). Jest możliwe obejrzeć też profil wzdłuż linii, jeśli taką się zaznaczy na obrazie. | ||
+ | |||
+ | # Zaobserwuj jak wpłynie na rozdzielczość dodanie szumu do każdego z tych fantomów: Process->Noise->Add noise. | ||
+ | # Wykonaj analizę FFT fantomów: Process->FFT->FFT. Pomnóż uzyskany obraz przez funkcję, która usunie niskie częstości lub wysokie częstości i wykonaj odwrotną transformację Fouriera: Process->FFT->Inverse FFT. | ||
+ | Np. pomnożenie przez dwuwymiarową funkcję Gaussa można wykonać tak: | ||
+ | '''v=v*exp(-(pow(64-x,2)+pow(64-y,2))/100)''' | ||
+ | |||
+ | Filtrowanie za pomocą splotu | ||
+ | # utwórz nowy obraz o wymiarach 15 x 15 | ||
+ | # wypełnij go funkcją Gaussa, np. '''v=exp(-pow(d,2)/10)*255''' | ||
+ | # zamień obraz na liczby: Image->Transform->Image to results | ||
+ | # skopiuj całą tabelę (zaznacz w opcjach okienka wyników Results->Options, by nie kopiować nagłówków kolumn i wierszy) | ||
+ | # wróć do obrazka z fantomem | ||
+ | # wybierz Process->Filters->Convolve i wklej tam zawartość swojego Gaussa | ||
+ | |||
+ | ==2. Sinogram i algorytm wstecznej projekcji== | ||
+ | |||
+ | Poniżej zawartość dwóch makr. Pierwsze tworzy sinogram z istniejącego obrazu. Drugie tworzy rekonstrukcję z sinogramu. W międzyczasie można dokonać filtrowania sinogramu. | ||
+ | |||
+ | Filtrowanie polega na usunięciu z sinogramu niskich częstości poprzez: | ||
+ | # wykonanie szybkiej transformaty Fouriera (FFT) sinogramu | ||
+ | # pomnożenie transformaty przez funkcję <math> \frac{|u|}{\pi} </math>, gdzie <math>u</math> jest poziomą częstością przestrzenną z przedziału <math> 0 \div \pi</math> | ||
+ | # wykonanie odwrotnej transformaty Fouriera (Inverse FFT). | ||
+ | |||
+ | Zadaniem będzie popracować z paroma różnymi fantomami (jednym z nich będzie [https://en.wikipedia.org/wiki/Shepp%E2%80%93Logan_phantom fantom Sheppa-Logana]). Najlepiej zacząć pracę z niewielkim fantomem, dla którego szybko pójdą obliczenia. Zakłada się, że każdy fantom jest obrazem w skali szarości. | ||
+ | |||
+ | Pożądanym wynikiem powinno być automatyczne zapisanie 10 rekonstrukcji dla różnej ilości projekcji (powinien to zrobić komputer, nie chodzi o robienie każdej rekonstrukcji osobno i ręczne jej zapisanie). Jeśli to się komuś uda, znaczy, że zrozumiał jak pisać makra w ImageJ. Pod tym adresem można znaleźć [https://imagej.nih.gov/ij/developer/macro/macros.html help]. | ||
+ | |||
+ | <source lang='c'> | ||
+ | macro "sinogram" { | ||
+ | imax=getWidth(); | ||
+ | jmax=getHeight(); | ||
+ | kmax=360; | ||
+ | a=floor(sqrt(imax*imax+jmax*jmax)+1); | ||
+ | d_fi=2*PI/kmax; | ||
+ | |||
+ | sinogram=newArray(kmax*a); | ||
+ | for (i=0; i<kmax*a; i++) | ||
+ | sinogram[i]=0; | ||
+ | |||
+ | for (k=0; k<kmax; k++) { | ||
+ | fi=k*d_fi; | ||
+ | x0=-a/2*cos(fi); | ||
+ | y0=a/2*sin(fi); | ||
+ | |||
+ | for (i=0; i<imax; i++) | ||
+ | for (j=0; j<jmax; j++) { | ||
+ | l=(i-imax/2-x0)*cos(fi)-(j-jmax/2-y0)*sin(fi); | ||
+ | sinogram[k*a+floor(l)]+=getPixel(i,j); | ||
+ | } | ||
+ | |||
+ | } | ||
+ | |||
+ | maks=0; | ||
+ | for (i=0; i<kmax*a; i++) | ||
+ | maks=maxOf(sinogram[i],maks); | ||
+ | razy=(256*256-1)/maks; | ||
+ | |||
+ | newImage("Sinogram","16-bit",a,kmax,0); | ||
+ | for (k=0; k<kmax; k++) | ||
+ | for (i=0; i<a; i++) { | ||
+ | setPixel(i, k, sinogram[k*a+i]*razy); | ||
+ | } | ||
+ | |||
+ | } | ||
+ | </source> | ||
+ | <source lang='c'> | ||
+ | macro "rekonstrukcja" { | ||
+ | a=getWidth(); | ||
+ | kmax=getHeight(); | ||
+ | imax=floor(a*sqrt(2)/2); | ||
+ | jmax=imax; | ||
+ | a=floor(sqrt(imax*imax+jmax*jmax)+1); | ||
+ | d_fi=2*PI/kmax; | ||
+ | rekon=newArray(imax*jmax); | ||
+ | |||
+ | for (i=0; i<imax*jmax; i++) | ||
+ | rekon[i]=0; | ||
+ | for (k=0; k<kmax; k++) { | ||
+ | fi=k*d_fi; | ||
+ | x0=-a/2*cos(fi); | ||
+ | y0=a/2*sin(fi); | ||
+ | |||
+ | for (i=0; i<imax; i++) | ||
+ | for (j=0; j<jmax; j++) { | ||
+ | l=(i-imax/2-x0)*cos(fi)-(j-jmax/2-y0)*sin(fi); | ||
+ | rekon[i+j*imax]+=getPixel(floor(l),k); | ||
+ | } | ||
+ | } | ||
+ | |||
+ | maks2=0; | ||
+ | for (i=0; i<imax*jmax; i++) | ||
+ | maks2=maxOf(rekon[i],maks2); | ||
+ | razy2=(256*256-1)/maks2; | ||
+ | |||
+ | newImage("Rekonstrukcja "+kmax+" projekcji","16-bit",imax,jmax,0); | ||
+ | for (i=0; i<imax; i++) | ||
+ | for (j=0; j<jmax; j++) { | ||
+ | setPixel(i,j,rekon[i+j*imax]*razy2); | ||
+ | } | ||
+ | } | ||
+ | </source> |
Aktualna wersja na dzień 10:59, 22 lis 2017
powrót do strony głównej wykładu Obrazowanie Medyczne
Ćwiczenia z Obrazowania medycznego z wykorzystaniem programu ImageJ
Otwórz program ImageJ.
1. Rozdzielczość
Wykonaj fantom do analizy rozdzielczości. Na przykład taki:
W tym celu
- utwórz nowy obraz: File->New->Image i podaj jego rozmiary 128 x 128 punktów.
- wejdź do okna Process->Math->Macro i wpisz w linijce: v=128*sin(2*PI/32*x)+128. Jak widzisz takie funkcje jak sinus lub taka wielkość jak PI są w programie predefiniowane.(np. potęgę trzecią osiąga się funkcją pow(x,3), a pierwiastek sqrt(x) itp.)
- zapisz obraz we własnym katalogu
Zastanów się jak
- zmienić częstość przestrzenną fantomu. Stwórz i zapisz dwa inne fantomy o większej częstości.
- utworzyć fantom o symetrii kołowej, na przykład taki:
- utworzyć fantom o zmiennej częstości, na przykład taki:
Aby utworzyć fantom zawierający funkcję schodkową (a nie sinusoidalną) proponuję wykonać własne makro.
- wybierz Plugins->New->Macro
- skopiuj poniższy skrypt
- uruchom makro w okienku edytora: Macros->Run macro (ctrl+R)
//to jest komentarz. Poniżej znajduje się zawartość makra tworzacego fantom z funkcją schodkową
imax=512;
jmax=512;
T=32;
newImage ("rozdzielczosc", "8-bit black", imax, jmax, 1);
for (i=0; i<imax; i++)
for (j=0; j<jmax; j++) {
if (i%T<T/2)
setPixel(i, j, 255);
else
setPixel(i, j, 0);
}
Jeśli chcesz obejrzeć profil jakiejś części fantomu, możesz zaznaczyć obszar i wybrać Analyze->Plot profile (ctr+K). Jest możliwe obejrzeć też profil wzdłuż linii, jeśli taką się zaznaczy na obrazie.
- Zaobserwuj jak wpłynie na rozdzielczość dodanie szumu do każdego z tych fantomów: Process->Noise->Add noise.
- Wykonaj analizę FFT fantomów: Process->FFT->FFT. Pomnóż uzyskany obraz przez funkcję, która usunie niskie częstości lub wysokie częstości i wykonaj odwrotną transformację Fouriera: Process->FFT->Inverse FFT.
Np. pomnożenie przez dwuwymiarową funkcję Gaussa można wykonać tak: v=v*exp(-(pow(64-x,2)+pow(64-y,2))/100)
Filtrowanie za pomocą splotu
- utwórz nowy obraz o wymiarach 15 x 15
- wypełnij go funkcją Gaussa, np. v=exp(-pow(d,2)/10)*255
- zamień obraz na liczby: Image->Transform->Image to results
- skopiuj całą tabelę (zaznacz w opcjach okienka wyników Results->Options, by nie kopiować nagłówków kolumn i wierszy)
- wróć do obrazka z fantomem
- wybierz Process->Filters->Convolve i wklej tam zawartość swojego Gaussa
2. Sinogram i algorytm wstecznej projekcji
Poniżej zawartość dwóch makr. Pierwsze tworzy sinogram z istniejącego obrazu. Drugie tworzy rekonstrukcję z sinogramu. W międzyczasie można dokonać filtrowania sinogramu.
Filtrowanie polega na usunięciu z sinogramu niskich częstości poprzez:
- wykonanie szybkiej transformaty Fouriera (FFT) sinogramu
- pomnożenie transformaty przez funkcję [math] \frac{|u|}{\pi} [/math], gdzie [math]u[/math] jest poziomą częstością przestrzenną z przedziału [math] 0 \div \pi[/math]
- wykonanie odwrotnej transformaty Fouriera (Inverse FFT).
Zadaniem będzie popracować z paroma różnymi fantomami (jednym z nich będzie fantom Sheppa-Logana). Najlepiej zacząć pracę z niewielkim fantomem, dla którego szybko pójdą obliczenia. Zakłada się, że każdy fantom jest obrazem w skali szarości.
Pożądanym wynikiem powinno być automatyczne zapisanie 10 rekonstrukcji dla różnej ilości projekcji (powinien to zrobić komputer, nie chodzi o robienie każdej rekonstrukcji osobno i ręczne jej zapisanie). Jeśli to się komuś uda, znaczy, że zrozumiał jak pisać makra w ImageJ. Pod tym adresem można znaleźć help.
macro "sinogram" {
imax=getWidth();
jmax=getHeight();
kmax=360;
a=floor(sqrt(imax*imax+jmax*jmax)+1);
d_fi=2*PI/kmax;
sinogram=newArray(kmax*a);
for (i=0; i<kmax*a; i++)
sinogram[i]=0;
for (k=0; k<kmax; k++) {
fi=k*d_fi;
x0=-a/2*cos(fi);
y0=a/2*sin(fi);
for (i=0; i<imax; i++)
for (j=0; j<jmax; j++) {
l=(i-imax/2-x0)*cos(fi)-(j-jmax/2-y0)*sin(fi);
sinogram[k*a+floor(l)]+=getPixel(i,j);
}
}
maks=0;
for (i=0; i<kmax*a; i++)
maks=maxOf(sinogram[i],maks);
razy=(256*256-1)/maks;
newImage("Sinogram","16-bit",a,kmax,0);
for (k=0; k<kmax; k++)
for (i=0; i<a; i++) {
setPixel(i, k, sinogram[k*a+i]*razy);
}
}
macro "rekonstrukcja" {
a=getWidth();
kmax=getHeight();
imax=floor(a*sqrt(2)/2);
jmax=imax;
a=floor(sqrt(imax*imax+jmax*jmax)+1);
d_fi=2*PI/kmax;
rekon=newArray(imax*jmax);
for (i=0; i<imax*jmax; i++)
rekon[i]=0;
for (k=0; k<kmax; k++) {
fi=k*d_fi;
x0=-a/2*cos(fi);
y0=a/2*sin(fi);
for (i=0; i<imax; i++)
for (j=0; j<jmax; j++) {
l=(i-imax/2-x0)*cos(fi)-(j-jmax/2-y0)*sin(fi);
rekon[i+j*imax]+=getPixel(floor(l),k);
}
}
maks2=0;
for (i=0; i<imax*jmax; i++)
maks2=maxOf(rekon[i],maks2);
razy2=(256*256-1)/maks2;
newImage("Rekonstrukcja "+kmax+" projekcji","16-bit",imax,jmax,0);
for (i=0; i<imax; i++)
for (j=0; j<jmax; j++) {
setPixel(i,j,rekon[i+j*imax]*razy2);
}
}