Obrazowanie:Obrazowanie Medyczne/ImageJ cwiczenia: Różnice pomiędzy wersjami

Z Brain-wiki
 
(Nie pokazano 18 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 8: Linia 10:
 
[[Plik:Fantom1.png]]
 
[[Plik:Fantom1.png]]
  
W tym celu  
+
W tym celu
 
# utwórz nowy obraz: File->New->Image i podaj jego rozmiary 128 x 128 punktów.
 
# 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 taką wielkość jak PI są w programie predefiniowane.
+
# 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
 
# zapisz obraz we własnym katalogu
  
Linia 16: Linia 18:
 
# zmienić częstość przestrzenną fantomu. Stwórz i zapisz dwa inne fantomy o większej częstości.
 
# 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 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: Fantom1.png

W tym celu

  1. utwórz nowy obraz: File->New->Image i podaj jego rozmiary 128 x 128 punktów.
  2. 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.)
  3. zapisz obraz we własnym katalogu

Zastanów się jak

  1. zmienić częstość przestrzenną fantomu. Stwórz i zapisz dwa inne fantomy o większej częstości.
  2. utworzyć fantom o symetrii kołowej, na przykład taki: Fantom2.png
  3. utworzyć fantom o zmiennej częstości, na przykład taki:Fantom3.png

Aby utworzyć fantom zawierający funkcję schodkową (a nie sinusoidalną) proponuję wykonać własne makro.

  1. wybierz Plugins->New->Macro
  2. skopiuj poniższy skrypt
  3. 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.

  1. Zaobserwuj jak wpłynie na rozdzielczość dodanie szumu do każdego z tych fantomów: Process->Noise->Add noise.
  2. 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

  1. utwórz nowy obraz o wymiarach 15 x 15
  2. wypełnij go funkcją Gaussa, np. v=exp(-pow(d,2)/10)*255
  3. zamień obraz na liczby: Image->Transform->Image to results
  4. skopiuj całą tabelę (zaznacz w opcjach okienka wyników Results->Options, by nie kopiować nagłówków kolumn i wierszy)
  5. wróć do obrazka z fantomem
  6. 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:

  1. wykonanie szybkiej transformaty Fouriera (FFT) sinogramu
  2. 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]
  3. 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);
        }
}