PL196564B1 - Sposób i urządzenie do rekonstruowania rentgenowskich obrazów tomograficznych - Google Patents

Abstract

1. Sposób rekonstruowania rentgenowskich obrazów tomograficznych zawieraj acych dane z orbity ko lowej i dane z orbity liniowej, na podstawie strefy cienia okre slonej w przestrzeni Radona, utworzonej przez rzutowanie danych orbity ko lowej, przy czym stosuje si e metod e rekonstrukcji obrazu z zastosowaniem transformacji Fouriera (FT) oraz odwrotnej transformacji Fouriera, znamienny tym, ze przedstawia si e obraz uzyskany z orbity ko lowej w ko lowej obj eto sci bry ly w kszta lcie rogu w przestrzeni Radona oraz przedstawia si e obraz uzyskany z orbity liniowej w odpo- wiedniej liniowej obj eto sci bry ly w kszta lcie rogu w prze- strzeni Radona, po czym przekszta lca si e (FT) obrazy dwuwymiarowe (2D) w p laszczyznach pionowych ka zdej ko lowej obj eto sci bry ly w kszta lcie rogu w przestrzeni Ra- dona i ka zdej liniowej obj eto sci bry ly w kszta lcie rogu ……. 13. Urz adzenie do rekonstruowania rentgenowskich ob- razów tomograficznych zawieraj acych dane z orbity ko lowej i dane z orbity liniowej, na podstawie strefy cienia w prze- strzeni Radona utworzonej przez rzutowanie danych orbity ko lowej, zawieraj ace zród lo promieniowania rentgenow- skiego oraz procesor i po laczone z nim podsystemy obs lugi i przetwarzania obrazu, obejmuj ace pami ec i wy swietlacz, znamienne tym, ze zawiera po laczone z centralnym pro- cesorem (CPU) (16), podsystem do przedstawiania obrazu uzyskanego z orbity ko lowej w ko lowej obj eto sci bry ly w kszta lcie rogu, ………………………………………………. PL PL PL PL

Description

(21) Numer zgłoszenia: 338536 ^{(51) Int.Cl.}

G06T 11/00 (2006.01) A61B 6/03 (2006.01)

Urząd Patentowy Rzeczypospolitej Polskiej (22) Data zgłoszenia: 21.02.2000

Opis patentowy przedrukowano ze względu na zauważ one błędy (54) Sposób i urządzenie do rekonstruowania rentgenowskich obrazów tomograficznych

(30) Pierwszeństwo: 22.02.1999,US,09/253,706	(73) Uprawniony z patentu: General Electric Company,Schenectady,US
(43) Zgłoszenie ogłoszono: 28.08.2000 BUP 18/00	(72) Twórca(y) wynalazku: Wen-Tai Lin,Schenectady,US Mehmet Yavuz,Clifton Park,US
(45) O udzieleniu patentu ogłoszono: 31.01.2008 WUP 01/08	(74) Pełnomocnik: Balińska Ewa, POLSERVICE, Kancelaria Rzeczników Patentowych Sp. z o.o.

(57) 1. Sposób rekonstruowania rentgenowskich obrazów tomograficznych zawierających dane z orbity kołowej i dane z orbity liniowej, na podstawie strefy cienia określonej w przestrzeni Radona, utworzonej przez rzutowanie danych orbity kołowej, przy czym stosuje się metodę rekonstrukcji obrazu z zastosowaniem transformacji Fouriera (FT) oraz odwrotnej transformacji Fouriera, znamienny tym, że przedstawia się obraz uzyskany z orbity kołowej w kołowej objętości bryły w kształcie rogu w przestrzeni Radona oraz przedstawia się obraz uzyskany z orbity liniowej w odpowiedniej liniowej objętości bryły w kształcie rogu w przestrzeni Radona, po czym przekształca się (FT) obrazy dwuwymiarowe (2D) w płaszczyznach pionowych każdej kołowej objętości bryły w kształcie rogu w przestrzeni Radona i każdej liniowej objętości bryły w kształcie rogu.......

13. Urządzenie do rekonstruowania rentgenowskich obrazów tomograficznych zawierających dane z orbity kołowej i dane z orbity liniowej, na podstawie strefy cienia w przestrzeni Radona utworzonej przez rzutowanie danych orbity kołowej, zawierające źródło promieniowania rentgenowskiego oraz procesor i połączone z nim podsystemy obsługi i przetwarzania obrazu, obejmujące pamięć i wyświetlacz, znamienne tym, że zawiera połączone z centralnym procesorem (CPU) (16), podsystem do przedstawiania obrazu uzyskanego z orbity kołowej w kołowej objętości bryły w kształcie rogu,.......................................................

UKŁAD PRÓBKĄ	^CY
32	34

FFT

U

1Q

FIG. 1

PL 196 564 B1

Opis wynalazku

Przedmiotem wynalazku jest sposób i urządzenie do rekonstruowania rentgenowskich obrazów tomograficznych.

W tomografii komputerowej wykorzystują cej technologię rzutowania wiązki stożkowej stosuje się zwykle źródło promieniowania X. Orbita skanowania źródła promieniowania X to krzywa, wzdłuż której wierzchołek wiązki stożkowej przesuwa się podczas skanowania. W ustalonej odległości od źródła wiązki stożkowej rozmieszczony jest zbiór detektorów. Znane są sposoby ustalania zależności pomiędzy rzutami wiązki stożkowej i pierwszą pochodną trójwymiarowych (3D) transformat Radona takich rzutów. Teoretycznie każde położenie źródła „S” może dostarczyć dane Radona rozmieszczone na sferze, powłokę Radona. Odpowiednie sfery na kołowej orbicie wycinają przestrzeń Radona obiektu, określaną jako torus. Obszary wewnątrz najmniejszej sfery zawierającej obiekt, ale zewnętrzne względem torusa, są zwane strefą cienia. Strefa cienia charakteryzuje obszar, w którym brakuje danych Radona z orbity kołowej.

Przy użyciu orbity kołowej i liniowej stosuje się zwykle dane orbity liniowej do wypełnienia strefy cienia. Mimo że strefa cienia jest bardzo mała, w porównaniu do wielkości łącznych poprawnych danych Radona dostarczonych przez orbitę kołową, obecne sposoby rzutowania wiązki stożkowej są złożone obliczeniowo. Odpowiednio do tego potrzebne jest uproszczenie rekonstrukcji obrazu oparte na orbicie kołowej i liniowej.

W amerykańskim opisie patentowym nr US 5,444,792 ujawniono sposób i urz ądzenie do rekonstruowania trójwymiarowych obrazów obiektu z zastosowaniem dwóch orbit kołowych (tj. kołowych geometrycznych trajektorii dostępu), które to rozwiązanie ma zastosowanie do obrazowania obiektów w dziedzinie medycyny lub nieniszcz ą cych metodach inspekcji obiektów. Przy rekonstruowaniu obrazów trójwymiarowych obiektu, w tym rozwiązaniu przemieszcza się co najmniej jedną dwuwymiarową tablicę czujników, czułych na promieniowanie przechodzące przez obiekt lub pochodzących z obiektu, w co najmniej dwóch kołowych trajektoriach, uzyskują c co najmniej jeden zestaw danych pomiarowych. Stosuje się cyfrową metodę umożliwiającą syntezę uzyskanych informacji (zmierzonych danych), a w szczególności konwersję uzyskanych danych pomiarowych w odpowiednich ramkach (granicach), powiązanych z trajektoriami w obrębie odpowiedniej dyskretnej ramki. Dwuwymiarową tablicę czujników stosuje się do uzyskania co najmniej dwóch zestawów danych pomiarowych w funkcji wybranych punktów obiektu określonych przez współrzędne w ramce odniesienia obiektu, która to funkcja reprezentuje obiekt i wynika z własności promieniowania odbieranego przez czujnik, przy czym co najmniej jedna tablica przemieszcza się wzdłuż co najmniej dwóch kołowych trajektorii wokół obiektu. Czujniki każdej, co najmniej jednej tablicy czujników są zogniskowane na odpowiednim wspólnym punkcie ogniskowania usytuowanym poza obiektem, w stałym położeniu w stosunku do tych czujników tak, aby każdy punkt ogniskowania przesuwał się wzdłuż odpowiedniej kołowej trajektorii a czujniki „widziały” stożkową przestrzeń i jedną ramkę trajektorii odniesienia oraz jeden zestaw danych pomiarowych, związany z każdą z trajektorii ogniskowych. Rekonstrukcja obrazu z zestawów danych pomiarowych z tablic czujników obejmuje sumowanie komputerowe tych zestawów danych pomiarowych na płaszczyznach odpowiedniego zestawu w dziedzinie dwuwymiarowej, a odpowiedni zespół płaszczyzn jest określony przez współrzędne trajektorii ramki odniesienia, łączenie uzyskanych danych sumowania na płaszczyznach w ramce odniesienia obrazu i przekształcanie odwrotne połączonych danych sumowania w ramce odniesienia obrazu w celu uzyskania zrekonstruowanego obrazu obiektu.

Sposób rekonstruowania rentgenowskich obrazów tomograficznych zawierających dane z orbity kołowej i dane z orbity liniowej, na podstawie strefy cienia określonej w przestrzeni Radona, utworzonej przez rzutowanie danych orbity kołowej, przy czym stosuje się metodę rekonstrukcji obrazu z zastosowaniem transformacji Fouriera FT oraz odwrotnej transformacji Fouriera, według wynalazku charakteryzuje się tym, że przedstawia się obraz uzyskany z orbity kołowej w kołowej objętości bryły w kształcie rogu w przestrzeni Radona oraz przedstawia się obraz uzyskany z orbity liniowej w odpowiedniej liniowej objętości bryły w kształcie rogu w przestrzeni Radona, po czym przekształca się (FT) obrazy dwuwymiarowe (2D) w płaszczyznach pionowych każdej kołowej objętości bryły w kształcie rogu w przestrzeni Radona i każdej liniowej objętości bryły w kształcie rogu w przestrzeni Radona na zbiór dwuwymiarowych płaszczyzn (2D) Fouriera, usuwa się rzuty trójwymiarowego (3D) stożka cienia w każ dej z dwuwymiarowych płaszczyzn (2D) Fouriera z danych orbity kołowej i zastępuje się te rzuty odpowiednimi rzutami uzyskanymi z danych orbity liniowej oraz na tej podstawie generuje się zbiór

PL 196 564 B1 poprawionych dwuwymiarowych płaszczyzn (2D) Fouriera, po czym przeprowadza się odwrotną dwuwymiarową transformację Fouriera na zbiorze poprawionych dwuwymiarowych płaszczyzn (2D) Fouriera, uzyskując zbiór dwuwymiarowych płaszczyzn (2D) obrazu reprezentujący otrzymaną objętość bryły w kształcie rogu w przestrzeni Radona oraz przekształca się tę otrzymaną objętość bryły w kształ cie rogu z powrotem na rentgenowski obraz tomograficzny.

Korzystnie, w trakcie przedstawiania obrazu uzyskanego z orbity kołowej wybiera się ponadto stożek cienia mający przeciwległy kąt 2θ₀, przy czym kąt θ₀ określa się jako kąt, który wyznacza granicę strefy cienia w przestrzeni Radona.

Korzystnie, w trakcie przedstawiania obrazu uzyskanego z orbity liniowej generuje się zbiór pionowych wycinków obrazu, w płaszczyznach pionowych, z danych orbity liniowej, tworząc liniową objętość bryły w kształcie rogu.

Korzystnie, w trakcie generowania zbioru pionowych wycinków obrazu, w płaszczyznach pionowych, rekonstruuje się zbiór równoległych wiązek wygenerowanych z rzutu danych orbity liniowej.

Korzystnie, do zrekonstruowania zbioru pionowych wycinków obrazu, w płaszczyznach pionowych, stosuje się metodę projekcji wstecznej z filtrowaniem (FBP).

Korzystnie, do zrekonstruowania zbioru pionowych wycinków obrazu, w płaszczyznach pionowych, używa się bezpośredniej metody Fouriera (DFM).

Korzystnie, w trakcie użycia bezpośredniej metody Fouriera (DFM) stosuje się odwrotną interpolację dwuliniową do przekształcenia danych z siatki biegunowej na siatkę prostokątną.

Korzystnie, wykonuje się odwrotną interpolację dwuliniową na każdej wartości „x” danych Fouriera orbity liniowej względem czterech najbliższych miejsc na siatce gi,j, gi+1,j, gi,j+1 i gi+1,j+1.

Korzystnie, w trakcie wykonywania odwrotnej interpolacji dwuliniowej na każdej wartości „x” danych orbity liniowej wyznacza się znormalizowany podobszar, aoo, a10, a01, a11, pomiędzy położeniem każdej odpowiedniej wartości „x” danych orbity liniowej i każdym z czterech najbliższych miejsc na siatce gi,j, gi+1,j, gi,j+1 i gi+1,j+1.

Korzystnie, dla każdego punktu siatki gi,i sąsiadującego z odpowiednią wartością „x” danych orbity liniowej pobiera się i przypisuje wagę wij.

Korzystnie, w trakcie wyznaczania znormalizowanego podobszaru aoo, a10, a01, a11 aktualizuje się każde miejsce na siatce gi,j, gi+1,j, gi,j+1 i gi+1,j+1 zgodnie z następującymi równościami:

g'i,j = gi,j + a11x g'i,1,j = gi+1,j + a01x ^g'i,j+1 ^{= g}i,j+1 ^{+ a}10^{x g'}i+1,j+1 ^{= g}i+1,j+1 ^{+ a}oo^x

Korzystnie, normalizuje się każdą daną siatki gi,j dzieląc odpowiednie gi,j przez wij.

Urządzenie do rekonstruowania rentgenowskich obrazów tomograficznych zawierających dane z orbity kołowej i dane z orbity liniowej, na podstawie strefy cienia w przestrzeni Radona utworzonej przez rzutowanie danych orbity kołowej, zawierające źródło promieniowania rentgenowskiego oraz procesor i połączone z nim podsystemy obsługi i przetwarzania obrazu, obejmujące pamięć i wyświetlacz, według wynalazku charakteryzuje się tym, że zawiera połączone z centralnym procesorem CPU, podsystem do przedstawiania obrazu uzyskanego z orbity kołowej w kołowej objętości bryły w kształcie rogu, podsystem do przedstawiania obrazu uzyskanego z orbity liniowej w odpowiedniej liniowej objętości bryły w kształcie rogu, podsystem do przekształcania (FT) obrazów dwuwymiarowych (2D) w płaszczyznach pionowych każdej kołowej objętości bryły w kształcie rogu i każdej liniowej objętości bryły w kształcie rogu na zbiór dwuwymiarowych płaszczyzn (2D) Fouriera zawierający procesor szybkiej transformaty Fouriera (FFT), urządzenie próbkujące oraz tablicę detektorów obrazu, podsystem do usuwania trójwymiarowych rzutów (3D) stożka cienia z każdej z dwuwymiarowych płaszczyzn (2D) Fouriera z danych orbity kołowej i do zastępowania tych rzutów odpowiednimi rzutami uzyskanymi z danych orbity liniowej oraz generowania zbioru poprawionych dwuwymiarowych płaszczyzn (2D) Fouriera, obejmujący interpolator, a także podsystem do przeprowadzania odwrotnej dwuwymiarowej transformacji Fouriera na zbiorze poprawionych dwuwymiarowych płaszczyzn (2D) Fouriera dla uzyskania zbioru dwuwymiarowych płaszczyzn (2D) obrazu reprezentującego otrzymaną objętość bryły w kształcie rogu, a także podsystem do przekształcania tej otrzymanej objętości w kształcie rogu z powrotem na rentgenowski obraz tomograficzny.

Korzystnie, podsystem do przedstawiania obrazu uzyskanego z orbity kołowej zawiera układ do wybierania stożka cienia mającego przeciwległy kąt 20_o, przy czym kąt θ₀ jest określony jako kąt, który wyznacza granicę strefy cienia w przestrzeni Radona.

PL 196 564 B1

Korzystnie, podsystem do przedstawiania obrazu uzyskanego z orbity liniowej zawiera zespół do generowania zbioru pionowych wycinków obrazu z danych orbity liniowej.

Korzystnie, zespół do generowania zbioru pionowych wycinków obrazu zawiera układ do rekonstruowania zbioru równoległych wiązek, wygenerowanych z rzutu danych orbity liniowej.

Korzystnie, urządzenie zawiera układ do stosowania metody projekcji wstecznej z filtrowaniem FBP do zrekonstruowania zbioru pionowych wycinków obrazu.

Korzystnie, urządzenie zawiera układ do stosowania bezpośredniej metody Fouriera (DFM) do zrekonstruowania zbioru pionowych wycinków obrazu.

Korzystnie, układ do stosowania bezpośredniej metody Fouriera DFM zawiera układ do stosowania odwrotnej interpolacji dwuliniowej do przekształcenia danych z siatki biegunowej na siatkę prostokątną.

Zaleta rozwiązania według wynalazku polega na zapewnieniu systemu rekonstrukcji rentgenowskich obrazów tomograficznych w oparciu o orbitę kołową i liniową. System zawiera etapy, w których łączy się obrazy odtworzone oddzielnie z danych orbity kołowej oraz danych orbity liniowej w dziedzinie częstotliwości. W szczególności przekształca się strefę cienia na stożek cienia w dziedzinie częstotliwości, odwzorowuje się stożek cienia w trójwymiarowej przestrzeni Fouriera (3D) na zbiór dwuwymiarowych płaszczyzn (2D) Fouriera, usuwa się dane leżące w obszarze oznaczonym jako rzut stożka cienia z transformaty Fouriera (FT) danych orbity kołowej i wstawia się w ich miejsce dane z transformaty Fouriera (FT) danych orbity liniowej, przekształ ca się każ d ą dwuwymiarow ą pł aszczyznę (2D) Fouriera na odpowiednią płaszczyznę obrazu dwuwymiarowego, oraz przekształca się objętość w kształcie rogu z powrotem na objętość siatkową. Stosuje się także technikę interpolacji do rekonstruowania danych orbity liniowej przy użyciu bezpośredniej metody Fouriera (DFM).

Przedmiot wynalazku przedstawiono w przykładzie wykonania na rysunku, na którym fig. 1 przedstawia schemat blokowy urządzenia według niniejszego wynalazku, fig. 2 - przekrój poprzeczny orbity kołowej utworzonej przez urządzenie z fig. 1, fig. 3 - strefę cienia z fig. 2 i odpowiadający jej stożek cienia, fig. 4a i 4b przedstawiają odpowiednio przekrój poprzeczny i podłużny stożka cienia z fig. 2 rzutowany na pł aszczyznę , fig. 5 przedstawia etapy przetwarzania przyrostowego wykonywane w kierunku „z” przez urządzenie z fig. 1, fig. 6 - zbiór danych rzutu równoległego pod kątem θ_παχ uzyskany z M-tego wiersza tablicy detektorów przedstawionej na fig. 5, fig. 7 - zbiór danych rzutu równoległego pod kątem 0_min uzyskany z pierwszego wiersza tablicy detektorów przedstawionej na fig. 5, a fig. 8 - odwrotną interpolację dwuliniową z urządzenia z fig. 1.

Figura 1 przedstawia schemat blokowy urządzenia CT 10 do rekonstruowania rentgenowskich obrazów tomograficznych zgodny z przedstawionym przykładem wykonania wynalazku. Urządzenie CT 10 zawiera procesor centralny (CPU) 16, urządzenie próbkujące 12 oraz inne podsystemy obsługi i przetwarzania obrazu 14, 18, 20, 22. Urządzenie próbkujące 12 zawiera źródło promieniowania X 34 i tablicę detekcji obrazu 32 wraz z przyrządami uruchamiającymi. Urządzenia obsługi obrazu 14, 18, 22 zawierają procesor 14 szybkiej transformaty Fouriera (FFT), interpolator 18, pamięć 20 i wyświetlacz 22.

W dwuwymiarowej (2D) rekonstrukcji tomograficznej alternatywą dla projekcji wstecznej z filtrowaniem (FBP) jest projekcja wsteczna według bezpośredniej metody Fouriera (DFM). Projekcja wsteczna według bezpośredniej metody Fouriera (DFM) opiera się na twierdzeniu o wycinku rzutu, które mówi, że jednowymiarowa (ID) transformata Fouriera (FT) rzutu pod określonym kątem odpowiada przekrojowi poprzecznemu dwuwymiarowej FT obiektu pod tym samym kątem. DFM będzie stosowany jako alternatywa dla zwykłej rekonstrukcji z projekcją wsteczną z użyciem orbity liniowej.

W urządzeniu CT 10 wykorzystującym dane orbity kołowej i alternatywnie d_rane orbity kołowej i liniowej, które należy zrekonstruować, wyrażone jako funkcja określona przez f(r), można wyrazić następują_rcym rów_rnaniem_r: _r f(_rr) = fc0(r) + fc1(r)+ f1(r) gdzie (r)- jest wektorem fc0(_rr) , fc1(r)- oznaczają funkcję danych z orbit kołowych c0 i c1 f1(r)- funkcję danych z orbity_rliniowej_r

Pierwsze dwa elementy, fc0 (r) i f_c1(r), oblicza się przy użyciu danych i techniki skanowania wzdłuż orbity kołowej. Drugi element jest potrzebny tylko wtedy, gdy trzeci składnik oblicza się w przePL 196 564 B1 strzeni Radona. Ponieważ będziemy obliczać trzeci składnik w przestrzeni Fouriera, drugi składnik zostanie zignorowany. _r

Trzeci element f1(r) jest obliczany z danych orbity liniowej. Stosuje się go zwykle do wypełniania strefy cienia 28 określonej przez torus 24, która jest wycinana przez orbitę kołową wokół osi obrotu, oraz najmniejszą sferę zawierającą obiekt 40 (figura 2). Mimo że ta strefa cienia 28, która jest określona w przestrzeni Radona, jest mała w porównaniu do torusa 24, obecnie st_rosowana metoda projekcji wstecznej jest złożona obliczeniowo, co jest oczywiste ze złożoności f1(r) w równaniu 1 podanym poniżej.

f1(r)=1

4π (d+r x')

π

J <^^ΘΗ(^zo, ^Θ,¹⁾1=_Υο sin_&+Z₀cos Θ'

Θ=0 (1) gdzie

H(Ζο, Θ, 1) = cos (-)w(Ζο, Θ, 1) d² + 1² δ''Σζο(1, Θ) _d2 δ²1 δΣ ζο (1, Θ) Ί d² δ1 J

Σ_ζο (1, Θ) = JJΡ_ζο (Y, Z)δ(Y sin Θ + Z cos Θ - 1)dYdZ, w(z₀,e>,1) = {1 - gdy 2lz₀ cos Θ +z²0 cos² Θ - d² sin² Θ > 0 w przeciwnym przypadku oraz gdzie H(z_o^,1) jest ważoną drugą pochodną danych Radona obliczonych z przekroju wykonanego przy z=zo. Dla każdego punktu w przestrzeni obiektu, r, punkt jest związany z punktem (Yo,Zo) na płaszczyźnie detektora. Jeżeli narysuje się okrąg na płaszczyźnie detektora przy użyciu początku (0,0) i (Y_o,Z_o) jako dwóch końców średnicy, wszystkie poprawne dane Radona, H(z_o^,1), przecinane przez okrąg, spełniają warunek 1=Y_osinθ+z_ocosθ. Dla każdego punktu (Y_o,Z_o) istnieje odpowiadający mu zbiór danych Radona, H(z_o^,1) do zsumowania i następnie wstecznej projekcji na przestrzeń obiektu. Złożoność obliczeniowa tego procesu jest więc bliska O(N⁵).

Według wynalazku dokonuje się efektywnej obliczeniowo rekonstrukcji obrazów tomograficznych przy użyciu danych uzyskanych ze skanowania orbity liniowej obiektu 40 (przedstawionego na fig. 2). Najpierw powiększa się strefę cienia 28 w przestrzeni Radona i przekształca na stożek cienia 30 w przestrzeni częstotliwości (w dziedzinie częstotliwości figura 3). Po drugie proces uzupełniania wykonywany przez procesor CPU 16 upraszcza się od początkowej trójwymiarowej przestrzeni częstotliwości do wielu dwuwymiarowych wycinków przestrzeni częstotliwości (przekrojów).

Pierwszy etap opiera się na fakcie, że transformata Fouriera linii radialnej w trójwymiarowej przestrzeni Radona (3D) jest równoważna tej samej linii radialnej w trójwymiarowej przestrzeni Fouriera. Jeżeli linia radialna w trójwymiarowej przestrzeni Radona przecina strefę cienia 28, jej odpowiednik w trójwymiarowej przestrzeni Fouriera jest „zanieczyszczony” ponieważ brakuje pewnych danych Radona w linii radialnej w przestrzeni Radona. Ponieważ wszystkie linie radialne przecinające strefę cienia 28 tworzą stożek cienia 30, zbiór ich odpowiedników w trójwymiarowej przestrzeni Fouriera również tworzy stożek niewystarczający co do częstotliwości. Stożek niewystarczający co do częstotliwości (albo stożek zanieczyszczony) pierwszego składnika równania (1), który jest określony w przestrzeni Fouriera, będzie zastąpiony trzecim składnikiem równania (1).

Drugi etap polega na pokazaniu, że orbita liniowa odpowiada dostatecznej liczbie danych o częstotliwości, aby zastąpić zanieczyszczony stożek (wynikający z orbity kołowej), po czym następuje procedura zastępowania zanieczyszczonego stożka w dwuwymiarowej przestrzeni Fouriera.

Figura 5 pokazuje, że orbita skanowania liniowego może być zrealizowana przez utrzymywanie stacjonarnego źródła promieniowania X 34 i płyty detektora 32, przy przesuwaniu obiektu 40 wzdłuż osi „z”, która jest również osią obrotową orbity kołowej. Obiekt 40 jest pokazany jako poruszający się w „kierunku z” w nieciągłych przyrostach określonych jako przyrosty „kroku z” 36. Każda pionowa kolumna detektorów 32 wraz ze źródłem promieniowania X 34 określa płaszczyznę przechodzącą przez obiekt 40. W każdej płaszczyźnie pionowej każde położenie detektora wyznacza określony kąt, a dane zebrane przez ten detektor odpowiadają równoległemu rzutowi obiektu, spróbkowanego w przedziałach „kroku z” 36. Figura 6 pokazuje pewną liczbę równoległych rzutów obiektu 40 pod kątem θm=θma_x.

PL 196 564 B1 a figura 7 pokazuje równoległe rzuty obiektu 40 pod kątem θ_π=θπί_η. Te równoległe rzuty stosuje się do wypełnienia odpowiedniego brakującego segmentu kątowego w płaszczyźnie dwuwymiarowej FT. θίπ=θπίη to kąty wybrane tak, aby ilość danych dla stożka cienia 30 (fig. 3) stosowanego do zastępowania danych strefy cienia 28 obejmowała dane strefy cienia. Warunek ten jest spełniony gdy θ₀=θ_π, jak przedstawiono na fig. 2, gdzie θ₀ = sin^-1 przy czym „R” jest promieniem najmniejszej sfery zawierającej obiekt 40, a „d” jest promieniem orbity kołowej.

Stożek cienia 30 (fig. 3) w trójwymiarowej przestrzeni Fouriera (3D) może być rzutowany jako dwuwymiarowy stożek cienia na zbiór dwuwymiarowych płaszczyzn Fouriera (2D). Odwrotna transformata każdej z tych dwuwymiarowych płaszczyzn Fouriera (2D) odpowiada płaszczyźnie obrazu określonej przez źródło promieniowania „X” i kolumnę „N” detektorów 32 (fig. 5).

Najpierw niech F(u,v) będzie transformatą Fouriera (FT) funkcji dwuwymiarowej f(x,y). Jednowymiarowe dane Fouriera na x=xo wyznacza się następująco:

F_x0 (v) = JF(u, v)e^i2nvdu = J{JJf(x y)e^-j2n(xu+^yv’dxdy}e^i2nx°^Udu = JJf(x, y)e-^Jπ {J e^{--2nx x)U}du}dxdy = JJ f (x, y )δ(x - x₀ )e-^Jnvdxdy = J f (x₀, y )e^-jnvdy.

(2)

Podobnie dla transformaty Fouriera (FT) funkcji trójwymiarowej f(x,y,z) dwuwymiarową transformatę Fouriera wycinka obiektu dla x=xo wyznacza się jak następuje:

F_{x 0} (v, w ) = J F (u, v, w)e^{i 2nX0V}du = J { JJ f (x, y, z )e ^{-} 2π( xu}+'^zw) dxdydz}e^{i 2nX0U}du = JJ f (x, y, z)e-^12π(yv+^zw){Je~^}2π(x-X0)^udu}dxdydz = JJ f (x, y, ζ)δ(x - x₀ )e-^j2π(yv+^zw)dxdydz = J f (x₀, y, z)e^-J2π(yv+^zw)dydz.

(3)

W przypadku gdy żądana n-wymiarowa przestrzeń Fouriera nie jest ortogonalna do wymiaru n+1 powiększonej n+1-wymiarowej przestrzeni Fouriera, na przykład jeżeli żądana jednowymiarowa transformata Fouriera leży na linii określonej przez równanie xcosθ+ysinθ=p=0, właściwy jest obrót początkowych współrzędnych dwuwymiarowych, (x, y) na nowe współrzędne (x', y') o kąt θ. Niech F(u', v') będzie transformatą Fouriera f(x', y'). Wtedy F(u', v') także uzyskuje się przez obrót F(u,v) o kąt θ. Transformatę Fouriera (FT) wzdłuż linii xcosθ+ysinθ=p=0 wyznacza się następująco:

F_x._ _p (v') = J F (u',v ' )e^{j 2nPU}'dU (4)

Podobnie dwuwymiarowa transformata Fouriera wzdłuż płaszczyzny określonej przez xcos<p+ysin<p=p=0 jest dana przez

F_x-p (v', w) = J F (U, v', w )e^{j :xpu}du (5)

I znowu F(u',v',w) uzyskuje się przez obrót współrzędnych (u,v) o kąt φ, jak przedstawiono na figurze 4A, pozostawiając współrzędną „w” bez zmian. Obrócony układ współrzędnych jest przedstawiony graficznie na figurach 4a i 4b, gdzie nowe współrzędne, (u', v', w) są tak dobrane, że ,,u'” jest równoległe do normalnej rozważanej płaszczyzny pokazanej jako linia SO'. Tak więc dysk cienia 26, który jest przekrojem poprzecznym stożka cienia 30 (fig. 4a), przecinający płaszczyznę poziomą 27

PL 196 564 B1 (fig. 4b) jest rzutowany na płaszczyźnie SO' jako segment liniowy ze środkiem „O'” i mający punkty końcowe „f” i „e”. Długość odcinka liniowego „ fe ” jest równa średnicy wycinka przestrzeni stożka cienia 30 ze środkiem „O”.

Podsumowując, stożek cienia 30 (fig. 4a), określający obszar w przestrzeni Radona w trójwymiarowej przestrzeni Fouriera, jest rzutowany na dwuwymiarową podprzestrzeń Fouriera przez zastosowanie odwrotnej jednowymiarowej transformaty Fouriera (FT). Innymi słowy, zamiast zastępować stożek cienia w trójwymiarowej przestrzeni Fouriera (3D), skanowany obiekt 40 jest reprezentowany w szeregu dwuwymiarowych płaszczyzn Fouriera (2D) i zastępuje rzutowany dwuwymiarowy dysk cienia 26 na poszczególnych płaszczyznach.

Przy stosowaniu opisanej wyżej zasady przyjmuje się następujące parametry tablicy i skanera. Tablica detektorów 32 (fig. 5) będzie określona wyrażeniem D(MxN), przy czym tablica detektorów 32 ma „M” wierszy i „N” kolumn. Dane rzutu orbity kołowej będą określone wyrażeniem Pc(MxNxVc), gdzie „Vc” jest łączną liczbą rzutów. Dane rzutu orbity liniowej będą określone wyrażeniem P1(MxNxV1), gdzie „V1” jest także łączną liczbą rzutów. Obraz zrekonstruowany z danych orbity kołowej będzie miał postać Ic (KxKxK) albo Ics (KxNxK); pierwsza objętość to objętość trójwymiarowej bryły w jednostkach sześciennych, natomiast druga to objętość bryły „w kształcie rogu” określona przez źródło promieniowania X i płytę detektora. Obraz zrekonstruowany z danych orbity liniowej będzie miał postać I1(KxKxK) albo I1s(KxNxK), podobnie pierwsza objętość to objętość trójwymiarowej bryły w jednostkach sześciennych, a druga to obję tość bryły „w kształcie rogu”.

Obraz rekonstruuje się z danych orbity liniowej. W pierwszym kroku modyfikuje się dane rzutu orbity liniowej, aby uwzględnić efekt rozbieżności gdy zbiór promieni przechodzi od źródła promieniowania X 34 do tablicy detektora 32, co przedstawiono graficznie na fig. 5, 6 i 7. Rozbieżny promień wykryty na płycie detektora modyfikuje się mnożąc wartości wskazań detektora przez funkcję wagową, jak następuje:

W(m,n) _D_ yjD² + (m - m0)² a¹ + (n - n0)² b² gdzie „D” w tym przypadku jest odległością od źródła do detektora (fig. 5), a ISO jest rzutowane na (mo, no) płyty detektora 32, gdzie ISO oznacza „środek skanera” (isocenter), to jest punkt leżący wzdłuż osi obrotu, który ma tą samą współrzędną „z” co źródło. Wartości „a” i „b” to odstępy detektora, odpowiednio w wymiarze wierszy i kolumn. Odstępy detektora to standardowe rozdzielczości określające odległość pomiędzy środkami dwóch sąsiednich detektorów.

Następnie rekonstruuje się n-tą kątową płaszczyznę pionową. Aby zrekonstruować n-tą kątową płaszczyznę pionową, dane rzutu P1n(MxV1) pobiera się z n-tej kolumny tablicy P1, gdzie m-ty wiersz odpowiada stałemu kątowi wiązki 0_m.

Elementy wierszą Pm przesuwa się tak, aby środek obiektu tego pionowego wycinka został ^v1 przesunięty do v = — . Niech P_1n· będzie tablicą rzutu przesuniętego.

Aby powrócić do przestrzeni obiektu, można użyć każdego z dwóch sposobów. W pierwszej opcji stosuje się jednowymiarową FFT dla każdego wiersza P1n'. Umieszcza się każdy zbiór danych FFT na linii radialnej pod kątem 0_m dla m=1...M. Stosuje się odwrotną transformatę Fouriera FFT^-1 do przekształcenia zbioru danych Fouriera na przestrzeń obiektu.

W drugiej opcji filtrowany wiersz P1n' dla wszystkich „N” kolumn rzutuje się wstecznie (stosuje się projekcję wsteczną), aby uzyskać objętość bryły w kształcie rogu, zawierającą „N” pionowych płaszczyzn określonych przez źródło promieniowania X 34 i „N” kolumn detektora. Niech fn(x,y) będzie obrazem zrekonstruowanym na podstawie n-tej kolumny, a Qn tablicą rzutu filtrowanego, wtedy fn(x,y) \^θΜ -θΐ\

M ^Qnm^{(x cos} dm + y ^sin 3m).

(6) m=1

Wyniki zastosowania powyższego równania dla przesuniętych elementów wiersza są wyznaczane jako tablica danych objętości Ils mającej kształt rogu, gdzie

I1s - obraz zrekonstruowany z danych z rzutu orbity kołowej dla c = 1.

PL 196 564 B1

Końcowy obraz można zrekonstruować jednym z dwóch możliwych sposobów. Po pierwsze, stożek cienia 30 można zastąpić danymi uzyskanymi ze skanowania orbity liniowej w trójwymiarowej przestrzeni Fouriera (jak objaśniono poniżej). Alternatywnie, stożek cienia 30 można zastąpić odpowiednimi danymi w dwuwymiarowych przestrzeniach Fouriera i następnie interpolować z objętości bryły o kształcie rogu z powrotem do objętości trójwymiarowej bryły we współrzędnych sześciennych.

Najpierw przedstawia się zastąpienie stożka cienia 30 w przestrzeni trójwymiarowej Fouriera. W pierwszym kroku stosuje się algorytm Feldkampa do zrekonstruowania obrazu Ic związanego z torusem 24, z danych orbity kołowej. Następnie uzyskuje się z Ic trójwymiarową transformatę (FT) Fouriera F(Ic).

Bryłę o objętości w kształcie rogu I1S interpoluje się, aby uzyskać odpowiadającą jej objętość bryły trójwymiarowej w jednostkach sześciennych I1. Z objętości I1 uzyskuje się trójwymiarową transformatę Fouriera F(I1).

Część stożka cienia F(lc) zastępuje się jej odpowiednikiem z F(Ic), aby uzyskać poprawioną tablicę danych Fouriera. Następnie stosuje się odwrotną transformatę Fouriera dla poprawionej tablicy danych Fouriera, aby powrócić do przestrzeni trójwymiarowej we współrzędnych sześciennych.

Następnie przedstawiony jest alternatywny sposób, w którym następuje zastępowanie stożka cienia w dwuwymiarowych przestrzeniach Fouriera (2D). W pierwszym kroku znowu stosuje się algorytm Feldkampa, aby zrekonstruować obraz Ic z danych orbity kołowej. Następnie interpoluje się Ic do Ics, gdzie Is określa zbiór kątowych płaszczyzn pionowych (patrz fig. 5). IIs uzyskuje się z danych orbity liniowej, jak opisano wyżej.

Transformatę Fouriera F(lcs) i F(lIs) uzyskuje się z Ics i IIs. Za stożek cienia podstawia się w każdej płaszczyźnie dwuwymiarowej F(Ics) jego odpowiednik określony w F(lIs). Stosuje się odwrotną transformatę Fouriera dla każdego z poprawionych dwuwymiarowych zbiorów danych Fouriera. Na końcu objętość bryły w kształcie rogu interpoluje się z powrotem do objętości bryły w przestrzeni we współrzędnych sześciennych.

Opisana wyżej metoda DFM wymaga interpolacji ze współrzędnych biegunowych na współrzędne prostokątne.

Opisana wyżej rekonstrukcja wiązki stożkowej jest oparta na orbitach kołowych i liniowych, które wymagają, aby brakująca strefa cienia 28 danych orbity kołowej w dziedzinie Fouriera została zastąpiona przez dane orbity liniowej. Dla orientacji orbity liniowej pokazanej na fig. 5, obiekt 40 jest pokazany jako poruszający się w „kierunku z” w nieciągłych przyrostach określonych jako przyrosty „kroku z” 36. Każda pionowa kolumna detektorów 32 wraz ze źródłem promieniowania „X” 34 określa płaszczyznę przechodzącą przez obiekt 40. W każdej płaszczyźnie pionowej każde położenie detektora wyznacza określony kąt, a dane zebrane przez ten detektor odpowiadają równoległemu rzutowi obiektu, spróbkowanemu w przedziałach „kroku z” (fig. 5, 6 i 7). Figura 6 pokazuje pewną liczbę równoległych rzutów obiektu 40 pod kątem 0_m ⁼emax, a fig. 7 pokazuje równoległe rzuty obiektu 40 pod kątem θ^θ^η. Te równoległe rzuty stosuje się do wypełnienia odpowiedniego brakującego segmentu kątowego w płaszczyźnie dwuwymiarowej FT. θ^ i ^_in to kąty wybrane tak, aby ilość danych dla stożka cienia 30 (fig. 3) stosowanych do zastępowania danych strefy cienia 28 obejmowała dane strefy cienia. Warunek ten jest spełniony gdy %=)_m, jak przedstawiono na fig. 2.

Pierwsze podejście polega na rekonstrukcji filtrowanego obrazu projekcji wstecznej (FBP) dla każdej płaszczyzny kątowej przy użyciu równoległych rzutów ze zbioru ograniczonych kątowych rzutów i następnie obliczeniu transformaty 2D FT i jej zastosowaniu do zastąpienia brakującego segmentu kątowego. Inną alternatywą jest przekształcenie danych siatki biegunowej z ID FT równoległych rzutów na siatkę prostokątną bezpośrednio w dziedzinie 2D FT. Sposób ten jest podobny do interpolacji wymaganej w projekcji wstecznej według bezpośredniej metody Fouriera (DFM).

W orbicie liniowej dla każdej płaszczyzny uzyskuje się pewną liczbę rzutów (równą „M” detektorom w każdej kolumnie pionowej), dla zakresu kątowego od ^_in do θ_Π3Χ. Techniki interpolacji dostarczają metodologię rozwijania opisanych poniżej technik interpolacji, które są użyteczne w obliczeniach dla sposobów rekonstrukcji wiązki stożkowej. Najpierw niech μφ(χ, y) będzie obrazem w ustalonym układzie współrzędnych x-y, a μ_φ(X, y) niech reprezentuje obiekt w systemie współrzędnych X - y będącym obrotem x - y o kąt φ. Rzutem obiektu pod katem φ jest lDFTp^x) jest dany przez

PL 196 564 B1

P_t(u) = J p_t( e^-2nxudu (8)

Jeżeli Μφ(ρ,φ) oznacza transformatę 2D FT μ_φ(χ, y) we współrzędnych biegunowych, w oparciu o twierdzenie o wycinku rzutu,

Μ(ρ,φ) = Ρφ(ρ), gdzie ρ > 0, 0 < φ < 2 π (9)

Zakładając, że Μ (ρ,φ) jest znane, te dane transformaty FT we współrzędnych biegunowych mogą zostać przekształcone na współrzędne kartezjańskie przez wyrażenie M(ρcosφ, ρβΐπφ^Μφ/φ). Obiekt jest rekonstruowany przez odwrotną transformatę 2D FT przy użyciu wyrażenia:

μ( x, y) = J M_c (u, v )e^{- j 2(ux}+^vy )dudy (10)

Jednak Μ(ρ,φ) nie jest znane we wszystkich położeniach, a tylko w skończonym zbiorze dyskretnych punktów (ρ,,φιΟ. Problemem staje się wtedy interpolacja ze znanych wartości punktów biegunowych na wartości wymagane na prostokątnej siatce kartezjańskiej.

Do rozwiązywania tego problemu stosuje się istniejące twierdzenie. Najpierw niech x(t) będzie funkcją okresową zmiennej czasowej „t” w przestrzeni czasowej o okresie „T”. Jeżeli x(t) jest kątowo \ 1 ograniczone w paśmie do wartości I c— I, gdzie „c” jest stałą, a „K” jest liczbą całkowitą odnoszącą się do liczby próbek x(t), to x(t) można zapisać jako:

2K +1 x(t) = Σ x(Ts)

1=0

sin	T (t - 1Ts)
(2K + 1)sin	n (t - 1Ts)

gdzie T_s =^s 2K +1 (11)

Dowód wynika z twierdzenia Shannona o próbkowaniu, jak następuje:

x(t) = Σ x(nTs )sin c n=-»

N-1

T (t - nTs) ^Ts (12) = ... + Σ x((1 - N )Τχ )sin c

1=0

T-(t - (1 - N)Tx) ^Τ

N-1 + Σ x(1Ts )sin c

1=0

N-1

T- (t - T) ^T + Σ x((1 + N T )sin c

1=0

N-1 = Σ x(1Ts )sin c

1=0

T- (t - (1+N T) ^T

A (t - (1Ts + kT)) + ...

(13) gdzie liczba kolumn tablicy danych N = 2K + 1.

W interpolacji kątowej stosuje się funkcje Μ(ρ,φ), gdzie zmienna φ ma okres 2π. Ponieważ Μ(ρ,φ) jest kątowo ograniczone w paśmie do „K” (czyli współczynniki szeregu Fouriera spełniają cn=0 dla |n|>K), to

2K-1 Λ

M (ρφ) =Σ M \ρ k=0 V kn K +1 φkn

K +1) (14) gdzie

PL 196 564 B1 (15) σ(φ) = \ 2K +1 ?

^sinl—φ (2K +1) sin( 2 φ)

Następnie opisana jest interpolacja radialna. Jeżeli obiekt jest ograniczony przestrzennie do średnicy 2A, gdzie „A” jest równoważne z promieniem „R, określonym poprzednio jako promień najmniejszej sfery zawierającej obiekt 40, jego rzuty są także ograniczone do średnicy 2A i dlatego transformatę FT rzutu P_ł(u) można zrekonstruować jako ^ρΦ^(υ) = Σ I^{sin c} (16) n=-o>

przy użyciu twierdzenia Shannona o próbkowaniu. Po połączeniu interpolacji kątowej i radialnej:

M (ρ,φ) ! 2K +1 ,

Σ ΣΜ n=-co k=0 n kn 2A ’ K +1 sin c I 2A I ρ-l l 2A σ φkn K +1 (17)

Mc (u1, vm) oblicza się na przestrzeni kartezjańskiej jako:

Mc (U1,Vm ) = mJu

-1 +^vm^{,cos 1} u² + v² “m ¹ * m (18)

Aby zmniejszyć błąd interpolacji, można podjąć pewną liczbę kroków, jak opisano poniżej. Najpierw zostanie rozważony błąd nieodpowiedniego próbkowania w przestrzeni.

W orbicie liniowej dane rzutu równoległego są próbkowane w przedziałach o określonym „kroku z” 36. Ponieważ rzuty są ograniczone w przestrzeni, nie będą ograniczone w paśmie (co do szerokości pasma). Próbkowanie rzutów na określonym oknie „kroku z” (na przykład 2A) efektywnie zmniejszy szerokość pasma i wprowadzi błędy utożsamiania do dziedziny 2D FT. W konsekwencji „krok z” 36 jest dobrany jako wystarczająco mały, tak aby efektywna szerokość pasma,

B =-¹-(punktów na cykl) w zasadzie obejmowała częstotliwościową zawartość rzutów.

* krok _ z

Innym źródłem błędu jest nieodpowiednie próbkowanie w częstotliwości radialnej. Jeżeli obiekt jest ograniczony przestrzennie do średnicy-, jego rzuty także będą ograniczone w przestrzeni

2*A do 2A, tak więc FT jest próbkowana jednorodnie w przedziałach przynajmniej - [punktów na cykl],

2A aby uniknąć błędu utożsamiania.

Dla danych liniowych efektywna szerokość pasma jest określona jako „B”. Liczba radialnych *B 1 1 próbek „N” jest tak dobrana, że -=-<-. Również odstęp siatki w dziedzinie

N * krok _ z 2 * A

2D FT wynosi najwyżej-[punktów na cykl].

2A

Innym źródłem błędu jest niedostateczna liczba rzutów. Obiekty nieizotropiczne i niewycentrowane nie będą kątowo ograniczone w paśmie. Tak więc liczba rzutów może być bardzo duża. Jednak w praktyce liczba rzutów jest tak dobrana, że w dziedzinie 2D FT próbkowanie pomiędzy dwoma są1 siednimi liniami radialnymi będzie w przybliżeniu równe maksymalnemu odstępowi siatki - [punk2A tów na cykl].

PL 196 564 B1

Następnie opisana jest interpolacja w dziedzinie 2D FT. Dla danych orbity liniowej stosuje się odstęp siatki prostokątnej-[punktów na cykl] w dziedzinie 2D FT.

2A

Zwykle w realizacji komputerowej stosuje się szybką transformację Fouriera (FFT) do obliczania transformaty 1D FT każdego rzutu. Aby ulepszyć próbkowanie radialne wykonuje się uzupełnianie danych rzutu zerami przed zastosowaniem transformacji FFT. W tym opisie uzupełnianie zerami polega na dodawaniu zer do zbioru danych, aby zminimalizować zakłócenia po wykonaniu transformacji FFT. Po uzyskaniu próbkowania radialnego, stosuje się interpolację najbliższego sąsiada albo interpolację liniową, aby obliczyć wartości pomiędzy nowymi punktami próbkowania.

Jak wspomniano uprzednio, dla orbity liniowej wymaganych jest wiele rzutów dla małego segmentu kątowego (czyli pobiera się więcej próbek w kierunku kątowym w porównaniu do próbkowania radialnego). Jednym z problemów jest efektywne wykorzystanie tych danych. Prosty schemat interpolacji taki jak z najbliższym sąsiadem nie wykorzystuje wszystkich dostępnych danych. Alternatywnym podejściem jest ograniczenie liczby rzutów przez pogrupowanie wierszy i następnie zsumowanie wierszy z każdej grupy, aby zwiększyć stosunek sygnału do szumu (SNR) przy zastosowaniu techniki odwrotnej interpolacji dwuliniowej, jak opisano poniżej.

W rozwiązaniu według wynalazku stosuje się technikę odwrotnej interpolacji dwuliniowej dla danych w formacie siatki biegunowej. Dla opcjonalnej interpolacji w kierunku azymutalnym σ(φ) można zastosować funkcję określoną w równaniu (13). Jednak ta funkcja jest niepraktyczna obliczeniowo. Dla realizacji komputerowej trzeba obcinać wartości w funkcji interpolacyjnej, co powoduje zakłócenia.

Z obliczeniowego punktu widzenia najprostszym algorytmem interpolacji do zastosowania jest algorytm najbliższego sąsiada (NN), w którym każdemu pikselowi nadaje się wartość próbki, która jest najbliższa piksela. NN jest równoważne splataniu spróbkowanych danych przez funkcję prostokątną. Splot z funkcją prostokątną w dziedzinie przestrzennej jest równoważny pomnożeniu sygnału w dziedzinie częstotliwości przez funkcję sinc(.), która jest złym filtrem dolnoprzepustowym, ponieważ ma wyraźne występy boczne.

Interpolacja liniowa oznacza splatanie spróbkowanych danych przez okno trójkątne. Funkcja ta odpowiada w dziedzinie częstotliwości funkcji sinc²( . ), która ma lepszą właściwość filtrowania dolnoprzepustowego niż NN.

Metoda interpolacji NN opiera się na pojedynczym punkcie. Algorytm interpolacji liniowej interpoluje na podstawie dwóch najbliższych punktów. Zastosowanie trzech punktów do interpolacji daje dwa punkty po jednej stronie i tylko jeden po drugiej stronie. Sposób z funkcjami B-sklejanymi umożliwia interpolację na dwóch najbliższych punktach w każdym kierunku. Metoda sześciennych funkcji B-sklejanych obejmuje cztery sploty funkcji prostokątnej z metody NN. Ponieważ sześcienne funkcje B-sklejane są symetryczne, trzeba je zdefiniować tylko w przedziale (0,2). Matematycznie funkcje B-sklejane można zapisać jako:

f (x) = 2x³ - x² + 2 χ _e (0,1) _{- x} ³ ₄ f (x) = —+ x² - 2x + 2 x _e (1,2)

Równania sześciennych funkcji sklejanych określone powyżej można dostosować tak, aby pasowały do gładkiego wielomianu rzędu trzeciego w dostępnym zbiorze danych stożka cienia 30.

Inna funkcja interpolacyjna zawiera obcięcie σ(φ). Jak wspomniano wyżej, pożądane jest znalezienie idealnej funkcji interpolacyjnej w kątowym kierunku funkcji σ(φ). Jednak funkcja σ(φ) ma znaczną energię, więc nie zanika szybko ze wzrostem odległości, dlatego nie może być łatwo zrealizowana, tj. nie nadaje się do splotu w dziedzinie przestrzennej. Typowe jest obcięcie funkcji ma małym odcinku, ale przy obcięciu ignoruje się dużo energii wspomnianej wyżej. Ponadto obcięcie w dziedzinie przestrzennej spowoduje zakłócenia pierścieniowe gdy obiekt jest przekształcany z dziedziny przestrzennej na dziedzinę częstotliwości.

Interpolacja dwuliniowa jest metodą interpolacji dwuwymiarowej. W tej metodzie interpolacji oblicza się f(x, y) jako kombinację liniową f(n₁,n₂) w czterech najbliższych punktach dla n₁T < x < (n₁+1)T i n₂T < y < (n₂ + 1)T. Wartość interpolowana f(x, y) w metodzie dwuliniowej to:

PL 196 564 B1

F(x,y) = (1 - A_x) (1 - A_y) f (n₁, n₂) + (1 - A_x) A_yf (n₁, n₂ + 1) + A_x (1 - A_y) f (n₁ + 1, n₂) + A_xA_yf (n₁ + 1, n₂ + 1) (x -n₁T) (y -n₂T)) gdzie A _x =-—_ ¹ oraz A_y =—W odwrotnej interpolacji dwuliniowej przekształca się dane siatki biegunowej na dane siatki prostokątnej poprzez odwrotne ważenie dwuliniowe. Tak więc, każdy element siatki biegunowej rozkłada się ponownie na cztery najbliższe punkty siatki prostokątnej przy użyciu czterech podobszarów jako znormalizowanych wag. Metoda jest niezawodna w tym sensie, że początkowe dane nie muszą być rozkładane regularnie, o ile dane mają przynajmniej jeden element w każdym kwadracie. Tak więc niepożądane elementy danych spowodowane przez defekty albo rowki w tablicy detektorów 32 (fig. 5) można zignorować. Inną zaletą tego procesu jest to, że stosuje się wszystkie poprawne dane, co w konsekwencji daje lepszy stosunek sygnału do szumu (SNR).

Wartość danej na orbicie liniowej „x” i związane z nią punkty siatki stosowane w odwrotnej interpolacji dwuliniowej są pokazane graficznie na fig. 8. Wartość danej „x” jest rozkładana na cztery najbliższe punkty siatki jak następuje:

g'i,j = gi,j + a11x g'i+1,j = gi+1,j + a01x ^g'i,j+1 ^{= g}i,j+1 ^{+ a}10^{x g'}i+1,j+1 ^{= g}i+1,j+1 ^{+ a}oo^x

Każdy punkt siatki ma wartość skumulowaną, gi,j, która określa skumulowane ważone wartości „x” z czterech sąsiednich kwadratów g'i,j oznacza wartość zaktualizowaną). Aby znormalizować wartość skumulowaną w każdym punkcie siatki, trzeba utrzymywać inną zmienną, Wij, która wyznacza skumulowane wagi (aij) związane z dwuliniowym rozkładem piksela (i,j). Jest to faktycznie krok normalizacyjny na końcu procedury interpolacji, ponieważ każda wartość siatkowa jest dzielona przez skumulowaną wagę wij.

Przeprowadzono sześć zestawów symulacji przy użyciu opisanych metod interpolacji dla dziedziny 2D FT i dziedziny obrazu. Wyniki symulacji ilustrują odwrotną interpolację dwuliniową w porównaniu do metod interpolacji najbliższego sąsiada (NN), liniowej, z sześciennymi funkcjami B-sklejanymi i projekcją wsteczną z filtrowaniem (FBP). Wyniki są przestawione w tabelach I do VI.

Dla pierwszych trzech symulacji zastosowano interpolację liniową w kierunku radialnym. W kierunku azymutalnym zastosowano techniki interpolacji najbliższego sąsiada, interpolacji liniowej i z sześciennymi funkcjami B-sklejanymi. Rekonstrukcję segmentu kątowego w dziedzinie 2D FT symulowano przy użyciu danych orbity liniowej. Rozważanym wycinkiem był centralny wycinek przechodzący przez obiekt na fig. 5. Obiekt składał się z trzech poziomych dysków, tak więc miał ostre krawędzie w kierunku „z” (pionowym).

Metodę interpolacji z obcięciem σ(φ) zastosowano także w kierunku azymutalnym (z 8 sąsiednimi punktami). W kierunku radialnym zastosowano interpolację liniową. Interpolacja najbliższego sąsiada i interpolacja liniowa wygenerowały lepsze wyniki robocze niż interpolacja z sześciennymi funkcjami B-sklejanymi. Jednak dla sześciennych funkcji B-sklejanych stosowano tylko 4 punkty, a nie 8 punktów jak w interpolacji z obcięciem σ(φ). Tak więc sześcienna funkcja B-sklejana jest bardziej do przyjęcia obliczeniowo niż interpolacja z obcięciem σ(φ) dla 8 punktów, więc interpolacji z obcięciem a σ(φ) nie uwzględniono w porównaniu.

Oceniono także interpolację dwuliniową opartą na punktach biegunowych. Najpierw obliczono transformatę 2D FT rozważanego wycinka i zastosowano jako odniesienie podstawowe dla porównywania czterech metod podanych wyżej. W porównaniu zawarto także transformację 2D FT rekonstrukcji FBP. Rekonstrukcja FBP wymaga rekonstrukcji obrazu 2D FBP i następnie obliczenia transformaty 2D FT. Jednak biegunowa interpolacja dwuliniowa nie wymaga konwersji wstecznej na dziedzinę obrazu. Daje to w rezultacie korzyść obliczeniową w porównaniu do metody FBP.

Aby ocenić działanie każdej metody, obliczono interpolowane wartości siatki prostokątnej przy użyciu powyższych metod dla różnych liczb rzutów i próbek radialnych. Błąd obliczono za pomocą rzeczywistej transformaty 2D FT jak następuje: Norm1 (łączna różnica bezwzględna, tabela I); Norm2 (pierwiastek kwadratowy łącznej odległości kwadratowej, tabela II); oraz Norm/ (maksymalna różnica bezwzględna, tabela 3). Należy zauważyć, że określenia Norm1, Norm2 i Norm/ to standardowe określenia stosowane w analizie statystycznej.

PL 196 564 B1

T a b e l a I

	BŁĄD %NORM1 w dziedzinie 2D FT
liczba rzutów	256	256	64	16
liczba próbek radialnych	1280	256	256	256
najbliższy sąsiad	58,4	58,0	60,9	75,2
liniowa	49,7	49,5	51,6	67,1
sześcienna funkcja B-sklejana	45,7	45,6	47,9	64,9
odwrotna dwuliniowa	41,2	40,7	43,3	60,9
FBP	77,6	77,6	77,6	77,6

T a b e l a II

	BŁĄD %NORM2 w dziedzinie 2D FT
liczba rzutów	256	256	64	16
liczba próbek radialnych	1280	256	256	256
najbliższy sąsiad	31,6	31,6	33,0	38,9
liniowa	26,5	26,5	26,3	33,8
sześcienna funkcja B-sklejana	24,4	24,5	24,8	32,9
odwrotna dwuliniowa	21,8	21,3	22,4	31,9
FBP	54,5	54,5	54,5	54,5

T a b e l a III

	BŁĄD %NORM/ w dziedzinie 2D FT
liczba rzutów	256	256	64	16
liczba próbek radialnych	1280	256	256	256
najbliższy sąsiad	8,88	8,88	8,88	8,88
liniowa	8,87	8,87	7,82	6,78
sześcienna funkcja B-sklejana	8,61	8,61	7,70	6,77
odwrotna dwuliniowa	7,58	6,08	6,08	7,26
FBP	80,25	80,25	80,25	80,25

Obliczono również odpowiednie obrazy ograniczonego rzutu (tabele IV-VI) przy użyciu odwrotnej transformaty 2D FT z interpolowanych siatek prostokątnych. Porównano błędy w obrazach zrekonstruowanych z obrazem rzeczywistym. W porównaniu zawarte są także obrazy FBP. Należy zauważyć, że w rzeczywistej realizacji dla orbity liniowej nie trzeba rekonstruować obrazów przed uzupełnieniem danych 2D FT w danych 2D FT Feldkampa.

PL 196 564 B1

T a b e l a IV

	BŁĄD %NORM1 w dziedzinie obrazu
liczba rzutów	256	256	64	16
liczba próbek radialnych	1280	256	256	256
najbliższy sąsiad	24,1	24,0	24,2	27,4
liniowa	19,1	19,1	18,4	21,2
sześcienna funkcja B-sklejana	16,1	16,1	15,7	18,9
odwrotna dwuliniowa	12,5	11,8	12,5	22,3
FBP	38,5	38,5	38,5	38,5

T a b e l a V

	BŁĄD %NORM2 w dziedzinie obrazu
liczba rzutów	256	256	64	16
liczba próbek radialnych	1280	256	256	256
najbliższy sąsiad	93,3	93,1	87,5	108,4
liniowa	52,8	52,6	50,1	77,6
sześcienna funkcja B-sklejana	50,8	50,7	49,4	75,7
odwrotna dwuliniowa	40,0	40,3	45,8	89,8
FBP	84,5	84,5	84,5	84,5

T a b e l a VI

	BŁĄD %NORMx w dziedzinie obrazu
liczba rzutów	256	256	64	16
liczba próbek radialnych	1280	256	256	256
najbliższy sąsiad	82,0	81,7	93,8	133,5
liniowa	67,8	67,5	78,4	121,5
sześcienna funkcja B-sklejana	66,5	66,2	77,7	117,0
odwrotna dwuliniowa	61,8	61,8	68,0	89,3
FBP	97,8	97,8	97,8	97,8

Można zauważyć, że zwiększenie próbkowania radialnego nie zmniejsza błędu w dziedzinie obrazu. Jednak zmniejszenie próbkowania radialnego wprowadza błędy, ponieważ odstępy siatki prostokątnej nie są już tego samego rzędu co w próbkowaniu radialnym. Jednak wartość „kroku z” 36 jest stała z powodu ograniczeń szerokości pasma. Dostępnych jest więc „N” próbek dla każdego rzutu, tak że „N” pomnożone przez „krok z” = 2A.

W każdym przypadku odwrotna interpolacja dwuliniowa ma mniejszy błąd zarówno w dziedzinie 2D FT, jak i w dziedzinie obrazu w porównaniu do innych podanych metod interpolacji. Liczbę rzutów w dziedzinie 2D FT i w dziedzinie obrazu można także zmniejszyć z 256, nie zwiększając znacząco błędu obrazu, co można na przykład zauważyć z danych w tabeli VI.

PL 196 564 B1

Claims (19)

1. Zastrzeżenia patentowe

   1. Sposób rekonstruowania rentgenowskich obrazów tomograficznych zawierających dane z orbity kołowej i dane z orbity liniowej, na podstawie strefy cienia określonej w przestrzeni Radona, utworzonej przez rzutowanie danych orbity kołowej, przy czym stosuje się metodę rekonstrukcji obrazu z zastosowaniem transformacji Fouriera (FT) oraz odwrotnej transformacji Fouriera, znamienny tym, że przedstawia się obraz uzyskany z orbity kołowej w kołowej objętości bryły w kształcie rogu w przestrzeni Radona oraz przedstawia się obraz uzyskany z orbity liniowej w odpowiedniej liniowej objętości bryły w kształcie rogu w przestrzeni Radona, po czym przekształca się (FT) obrazy dwuwymiarowe (2D) w płaszczyznach pionowych każdej kołowej objętości bryły w kształcie rogu w przestrzeni Radona i każdej liniowej objętości bryły w kształcie rogu w przestrzeni Radona na zbiór dwuwymiarowych płaszczyzn (2D) Fouriera, usuwa się rzuty trójwymiarowego (3D) stożka cienia (30) w każdej z dwuwymiarowych płaszczyzn (2D) Fouriera z danych orbity kołowej i zastępuje się te rzuty odpowiednimi rzutami uzyskanymi z danych orbity liniowej oraz na tej podstawie generuje się zbiór poprawionych dwuwymiarowych płaszczyzn (2D) Fouriera, po czym przeprowadza się odwrotną dwuwymiarową transformację Fouriera na zbiorze poprawionych dwuwymiarowych płaszczyzn (2D) Fouriera, uzyskując zbiór dwuwymiarowych płaszczyzn (2D) obrazu reprezentujący otrzymaną objętość bryły w kształcie rogu w przestrzeni Radona oraz przekształca się tę otrzymaną objętość bryły w kształcie rogu z powrotem na rentgenowski obraz tomograficzny.
2. 2. Sposób według zastrz. 1, znamienny tym, że w trakcie przedstawiania obrazu uzyskanego z orbity kołowej wybiera się ponadto stożek cienia (30) mający przeciwległy kąt 2%, przy czym kąt % określa się jako kąt, który wyznacza granicę strefy cienia (28) w przestrzeni Radona.
3. 3. Sposób według zastrz. 1, znamienny tym, że w trakcie przedstawiania obrazu uzyskanego z orbity liniowej generuje się zbiór pionowych wycinków obrazu, w płaszczyznach pionowych, z danych orbity liniowej, tworząc liniową objętość bryły w kształcie rogu.
4. 4. Sposób według zastrz. 3, znamienny tym, że w trakcie generowania zbioru pionowych wycinków obrazu, w płaszczyznach pionowych, rekonstruuje się zbiór równoległych wiązek wygenerowanych z rzutu danych orbity liniowej.
5. 5. Sposób według zastrz. 4, znamienny tym, że do zrekonstruowania zbioru pionowych wycinków obrazu, w płaszczyznach pionowych, stosuje się metodę projekcji wstecznej z filtrowaniem (FBP).
6. 6. Sposób według zastrz. 4, znamienny tym, że do zrekonstruowania zbioru pionowych wycinków obrazu, w płaszczyznach pionowych, używa się bezpośredniej metody Fouriera (DFM).
7. 7. Sposób według zastrz. 6, znamienny tym, że w trakcie użycia bezpośredniej metody Fouriera (DFM) stosuje się odwrotną interpolację dwuliniową do przekształcenia danych z siatki biegunowej na siatkę prostokątną.
8. 8. Sposób według zastrz. 7, znamienny tym, że wykonuje się odwrotną interpolację dwuliniową na każdej wartości „x” danych Fouriera orbity liniowej względem czterech najbliższych miejsc na siatce gi,j, gi+1,j, gi,j+1 i gi+1,j+1.
9. 9. Sposób według zastrz. 8, znamienny tym, że w trakcie wykonywania odwrotnej interpolacji dwuliniowej na każdej wartości „x” danych orbity liniowej wyznacza się znormalizowany podobszar, aoo, a10, a01, a11, pomiędzy położeniem każdej odpowiedniej wartości „x” danych orbity liniowej i każdym z czterech najbliższych miejsc na siatce gi,j, gi+1,j, gi,j+1 i gi+1,j+1.
10. 10. Sposób według zastrz. 9, znamienny tym, że dla każdego punktu siatki gi,i sąsiadującego z odpowiednią wartością „x” danych orbity liniowej pobiera się i przypisuje wagę wij.
11. 11. Sposób według zastrz. 10, znamienny tym, że trakcie wyznaczania znormalizowanego podobszaru aoo, a10, a01, a11 aktualizuje się każde miejsce na siatce gi,j, gi+1,j, gi,j+1 i gi+1,j+1 zgodnie z następującymi równościami:

    g'i,j = gi,j + a11x g'i,1,j = gi+1,j + a01x ^g'i,j+1 ^{= g}i,j+1 ^{+ a}10^{x g'}i+1,j+1 ^{= g}i+1,j+1 ^{+ a}oo^x
12. 12. Sposób według zastrz. 11, znamienny tym, że normalizuje się każdą daną siatki gi,j dzieląc odpowiednie gi,j przez wi,j.
13. 13. Urządzenie do rekonstruowania rentgenowskich obrazów tomograficznych zawierających dane z orbity kołowej i dane z orbity liniowej, na podstawie strefy cienia w przestrzeni Radona utworzonej przez rzutowanie danych orbity kołowej, zawierające źródło promieniowania rentgenowskiego

    PL 196 564 B1 oraz procesor i połączone z nim podsystemy obsługi i przetwarzania obrazu, obejmujące pamięć i wyświetlacz, znamienne tym, że zawiera połączone z centralnym procesorem (CPU) (16), podsystem do przedstawiania obrazu uzyskanego z orbity kołowej w kołowej objętości bryły w kształcie rogu, podsystem do przedstawiania obrazu uzyskanego z orbity liniowej w odpowiedniej liniowej objętości bryły w kształcie rogu, podsystem do przekształcania (FT) obrazów dwuwymiarowych (2D) w płaszczyznach pionowych każdej kołowej objętości bryły w kształcie rogu i każdej liniowej objętości bryły w kształcie rogu na zbiór dwuwymiarowych płaszczyzn (2D) Fouriera zawierający procesor (14) szybkiej transformaty Fouriera (FFT), urządzenie próbkujące (12) oraz tablicę detektorów obrazu (32), podsystem do usuwania trójwymiarowych rzutów (3D) stożka cienia (30) z każdej z dwuwymiarowych płaszczyzn (2D) Fouriera z danych orbity kołowej i do zastępowania tych rzutów odpowiednimi rzutami uzyskanymi z danych orbity liniowej oraz generowania zbioru poprawionych dwuwymiarowych płaszczyzn (2D) Fouriera, obejmujący interpolator (18), a także podsystem do przeprowadzania odwrotnej dwuwymiarowej transformacji Fouriera na zbiorze poprawionych dwuwymiarowych płaszczyzn (2D) Fouriera dla uzyskania zbioru dwuwymiarowych płaszczyzn (2D) obrazu reprezentującego otrzymaną objętość bryły w kształcie rogu, a także podsystem do przekształcania tej otrzymanej objętości w kształcie rogu z powrotem na rentgenowski obraz tomograficzny.
14. 14. Urządzenie według zastrz. 13, znamienne tym, że podsystem do przedstawiania obrazu uzyskanego z orbity kołowej zawiera układ do wybierania stożka cienia (30) mającego przeciwległy kąt 2θ0, przy czym kąt θ0 jest określony jako kąt, który wyznacza granicę strefy cienia (28) w przestrzeni Radona.
15. 15. Urządzenie według zastrz. 13, znamienne tym, że podsystem do przedstawiania obrazu uzyskanego z orbity liniowej zawiera zespół do generowania zbioru pionowych wycinków obrazu z danych orbity liniowej.
16. 16. Urządzenie według zastrz. 15, znamienne tym, że zespół do generowania zbioru pionowych wycinków obrazu zawiera układ do rekonstruowania zbioru równoległych wiązek, wygenerowanych z rzutu danych orbity liniowej.
17. 17. Urządzenie według zastrz. 16, znamienne tym, że zawiera układ do stosowania metody projekcji wstecznej z filtrowaniem (FBP) do zrekonstruowania zbioru pionowych wycinków obrazu.
18. 18. Urządzenie według zastrz. 16, znamienne tym, że zawiera układ do stosowania bezpośredniej metody Fouriera (DFM) do zrekonstruowania zbioru pionowych wycinków obrazu.
19. 19. Urządzenie według zastrz. 18, znamienne tym, że układ do stosowania bezpośredniej metody Fouriera (DFM) zawiera układ do stosowania odwrotnej interpolacji dwuliniowej do przekształcenia danych z siatki biegunowej na siatkę prostokątną.