Generowanie powierzchni z wykorzystaniem algorytmu Marching Cube w VTK

Kolejnym sposobem - po wizualizacji objętościowej (ang. volume rendering) - renderowania danych przestrzennych jest technika generowania powierzchni (ang. surface rendering) z wykorzystaniem algorytmu Marching Cubes. Pozwala ona na wyświetlanie konturu wybranych obiektów i tworzenie dla nich dyskretnych modeli 3D. Obiekty te są reprezentowane za pomocą siatek wypełnionych trójkątami. Otrzymujemy w ten sposób model powierzchniowy obiektu, dla którego możemy dodatkowo zdefiniować stopień przeźroczystości i kolor.

Przykład
Użycie algorytmu Marching Cubes do generowania konturu w VTK jest niezwykle proste. Po odczytaniu danych z pliku wystarczy przekazać je do obiektu vtkMarchingCubes i zdefiniować wartość progową dla obiektu.

package require vtk

set inputFileName  [ lindex $argv 0]
set isoLevel [expr [ lindex $argv 1] ]
set opacity  [expr [ lindex $argv 2] ]
set color "0.9804 0.9216 0.8431"

vtkStructuredPointsReader reader
  reader SetFileName $inputFileName
  reader Update

vtkRenderer ren1
vtkRenderWindow renWin
  renWin AddRenderer ren1
vtkRenderWindowInteractor iren
  iren SetRenderWindow renWin

vtkMarchingCubes iso
  iso SetInput [reader GetOutput]
  iso SetValue 0 $isoLevel

vtkPolyDataMapper isoMapper
  isoMapper SetInput [iso GetOutput]
  isoMapper ScalarVisibilityOff

vtkActor isoActor
  isoActor SetMapper isoMapper
  eval [isoActor GetProperty] SetColor $color
  eval [isoActor GetProperty] SetOpacity $opacity

ren1 AddActor isoActor
ren1 SetBackground 0.5 0.6 0.8
ren1 SetGradientBackground 4
renWin SetSize 500 500

iren Initialize
wm withdraw .

Powyższy skrypt zawiera obsługę parametrów wejściowych. Dla przykładowych danych (tooth.vtk) wywołanie może wyglądać w następujący sposób:

tclsh vtk-iso_contour.tcl data/tooth.vtk 128 0.5

Analizując przykład, polecam zrobić eksperymenty na własnych danych zmieniając wartość konturu, kolor i przeźroczystość. Rezultaty jakie w ten sposób możemy uzyskać mogą być naprawdę różne, co widać już po samej zmianie stopnia przeźroczystości:


Możliwości VTK w zakresie generowania konturu nie kończą się na wyświetlaniu tylko jednej powierzchni, ale możemy zdefiniować sobie ich dowolną ilość. Może to być szczególnie przydatne np. w zastosowaniach klinicznych, gdzie każdy organ może zostać osobno wyświetlony.

Wadą techniki surface rendering, o której trzeba na koniec wspomnieć, jest jej wysokie zapotrzebowanie na moc obliczeniową i pamięć operacyjną. Może się okazać, że przy bardziej złożonych obiektach czas oczekiwania będzie na tyle długi, że warto będzie rozważyć skorzystanie z wizualizacji objętościowej.

Wizualizacja objętościowa w VTK

Dłuższy czas temu (rok 2007) napisałem na tym blogu wpis związany z wizualizacją medyczną. Dzisiaj nie będę powtarzał tych samych informacji, tylko skupię się na tym, jak można wyrenderować serię obrazów (np. z badania tomograficznego) w VTK za pomocą techniki Volume Rendering. Co najważniejsze, technika ta nie dość, że daje bardzo dobre rezultaty - co zobaczymy później na przykładzie - to działa nadzwyczaj sprawnie. Implementacja dostępna w VTK potrafi skorzystać z kilku rdzeni dostępnych w naszym procesorze, a podczas obrotów lub skalowania automatycznie zmniejsza rozdzielczość.

Najważniejsze klasy:
Przed rozpoczęciem analizowania przykładu warto najpierw zapoznać się z najważniejszymi klasami, z których będziemy korzystać:

• vtkPiecewiseFunction - klasa umożliwia definiowanie funkcji przeźroczystości dla danych,
  
• vtkColorTransferFunction - pozwala na definiowanie skali kolorów dla odpowiednich wartości,
  
• vtkVolumeProperty - zbiera wszystkie informacje dotyczące parametrów wyświetlania danych (funkcje koloru i przeźroczystości, sposób cieniowania i interpolacji, i wiele innych),
  
• vtkVolume - reprezentuje dane i właściwości ich wyświetlania, odpowiada za sterowanie pozycją i orientacją na scenie renderowania,
• vtkVolumeRayCastCompositeFunction i vtkVolumeRayCastMapper zawierają algorytmy, które na podstawie informacji zawartych w vtkVolumeProperty wiedzą jak zamienić obiekt 3D na obraz płaski, który w rzeczywistości widzimy na ekranie.

Przykład:
Poniżej znajduje się prosty skrypt TCL, który pokazuje jak wyrenderować serię zdjęć przedstawiających zęba za pomocą techniki Volume Rendering.
Przykładowe dane zostały zaczerpnięte z pakietu SciRun i przekonwertowane do formatu VTK Structured Points. Można je pobrać z następującej lokalizacji: link.

package require vtk
package require vtkinteraction

# Create the reader for the data
vtkStructuredPointsReader reader
  reader SetFileName "data/tooth.vtk"
  reader Update

vtkImageCast cast
  cast SetInput [reader GetOutput]
  cast SetOutputScalarTypeToUnsignedShort
  cast ClampOverflowOn

# Create transfer mapping scalar value to opacity
vtkPiecewiseFunction opacityTransferFunction
  opacityTransferFunction AddPoint    0.0  0.0
  opacityTransferFunction AddPoint   90.0  0.0
  opacityTransferFunction AddPoint  137.9  0.119
  opacityTransferFunction AddPoint  255.0  0.2

# Create transfer mapping scalar value to color
vtkColorTransferFunction colorTransferFunction
  colorTransferFunction AddRGBPoint      0.0 0.0 0.0 1.0
  colorTransferFunction AddRGBPoint    120.0 1.0 1.0 1.0
  colorTransferFunction AddRGBPoint    160.0 1.0 1.0 0.0
  colorTransferFunction AddRGBPoint    200.0 1.0 0.0 0.0
  colorTransferFunction AddRGBPoint    255.0 0.0 1.0 1.0

# Create transfer mapping scalar value to opacity based on gradient magnitude
vtkPiecewiseFunction gradientTransferFunction
  gradientTransferFunction AddPoint     0.0 0
  gradientTransferFunction AddPoint     2.5 0
  gradientTransferFunction AddPoint    12.7 1
  gradientTransferFunction AddPoint   255.0 1

# The property describes how the data will look
vtkVolumeProperty volumeProperty
  volumeProperty SetColor colorTransferFunction
  volumeProperty SetScalarOpacity opacityTransferFunction
  volumeProperty SetGradientOpacity gradientTransferFunction
  #volumeProperty ShadeOn
  volumeProperty SetInterpolationTypeToLinear

# The mapper / ray cast function know how to render the data
vtkVolumeRayCastCompositeFunction  compositeFunction
vtkVolumeRayCastMapper volumeMapper
  volumeMapper SetVolumeRayCastFunction compositeFunction
  volumeMapper SetInputConnection [cast GetOutputPort]

# The volume holds the mapper and the property and
# can be used to position/orient the volume
vtkVolume volume
  volume SetMapper volumeMapper
  volume SetProperty volumeProperty

vtkRenderer ren1
  ren1 SetBackground 0.5 0.6 0.8
  ren1 SetGradientBackground 4
vtkRenderWindow renWin
  renWin AddRenderer ren1
  renWin SetSize 600 600
vtkRenderWindowInteractor iren
  iren SetRenderWindow renWin
  
ren1 AddVolume volume
renWin Render

iren AddObserver UserEvent {wm deiconify .vtkInteract}
iren Initialize

wm withdraw .

Przykład ten został uruchomiony dwa razy - z wyłączoną i włączoną opcją przeźroczystości generowanej na podstawie gradientu (linijka: volumeProperty SetGradientOpacity gradientTransferFunction). W pierwszym przypadku widzimy struktury odsłonięte za pomocą funkcji opacityTransferFunction, natomiast w drugim przypadku możemy również zajrzeć do wnętrza tych struktur.



Przedstawiony przykład korzysta z bardzo prostego modelu, co wcale nie oznacza, że technika Volume Rendering nadaje się tylko to tego typu zastosowań. Jest ona często wykorzystywana (razem z samym VTK) w zastosowaniach klinicznych, gdzie mamy zdecydowanie więcej danych i często potrzebujemy skupić się na konkretnych tkankach lub organach. Ten sposób oceny obrazu daje zdecydowanie lepszy pogląd analizowanego obiektu badań, aniżeli przeglądanie dziesiątek, a nawet setek skanów 2D.

Jak skompilować InsightToolkit ze źródeł i podłączyć do TCLa i Pythona?

Podobnie jak w przypadku VTK instalacja InsightToolkit w Ubuntu 9.04 sprowadza się do wpisania w terminalu następującego polecenia:

sudo apt-get install insighttoolkit3-examples libinsighttoolkit3-dev libinsighttoolkit3.10 tcl8.4-insighttoolkit3 python-insighttoolkit3

Podobnie również tutaj, występuje problem z dostępem do najnowszej wersji biblioteki (wersja w repozytorium na dzień dzisiejszy - 6 paź '09 - to 3.10, a na stronie projektu jest już dostępna wersja 3.16). Dlatego jeśli zależy nam na najświeższej wersji, możemy ją sami przygotować kompilując bibliotekę ze źródeł.

Ze strony ITK, z działu Download pobieramy następujące pliki: InsightToolkit-3.16.0.tar.gz, InsightApplications-3.16.0.tar.gz, CableSwig-ITK-3.16.0.tar.gz, lub ich nowsze wersje. Po ich rozpakowaniu do wcześniej przygotowanego katalogu, najpierw skompilujemy CableSwig - pakiet używany do opakowania (ang. wrapping) kodu C++ do innych języków programowania. Nie wymaga on ustawiania żadnych niestandardowych opcji kompilacji, dlatego po przejściu do katalogu CableSwig-ITK-3.16.0 wystarczy wpisać następującą sekwencję poleceń:

mkdir build && cd build && cmake ../ && make

Następnie przechodzimy do katalogu InsightToolkit-3.16.0, zakładamy katalog build i po przejściu do niego uruchamiamy graficzną nakładkę na CMake:

mkdir build && cd build && cmake-gui ../

W kolejnym kroku ustawiamy odpowiednie opcje kompilacji. Poniższa propozycja umożliwi połączenie ITK z językami skryptowymi TCL i Python.

BUILD_EXAMPLES                 OFF
BUILD_SHARED_LIBS              ON
BUILD_TESTING                  OFF
CMAKE_BUILD_TYPE               Release
CMAKE_INSTALL_PREFIX           /opt
CableSwig_DIR                  /home/rafal/install/itk/CableSwig-ITK-3.16.0/build
ITK_USE_REVIEW                 ON
ITK_USE_PATENTED               ON
ITK_USE_REVIEW_STATISTICS      OFF
USE_WRAP_ITK                   ON
INSTALL_WRAP_ITK_COMPATIBILITY ON
WRAP_ITK_JAVA                  OFF
WRAP_ITK_PYTHON                ON
WRAP_ITK_TCL                   ON

Uruchamiamy kompilację i instalujemy bibliotekę:

make
sudo make install

Uwagi odnośnie kompilacji:

• Liczba przykładów i testów w ITK jest ogromna. Warto zatem je wyłączyć (opcje: BUILD_EXAMPLES i BUILD_TESTING) przy pierwszej kompilacji, bo w przeciwnym wypadku kompilacja może potrwać nawet kilka godzin. Jeśli proces kompilacji zakończy się poprawnie, kompilację przykładów polecam włączyć w wolnym czasie, np. w nocy.
  
• Opcje BUILD_SHARED_LIBS i ITK_USE_REVIEW są wymagane po włączeniu opcji USE_WRAP_ITK.
  
• W wersji ITK 3.16, którą testowałem, włączenie opcji ITK_USE_REVIEW_STATISTICS i USE_WRAP_ITK spowodowało błędy w kompilacji. Dlatego jeśli nam bardziej zależy na połączeniu z językami skryptowymi, polecam wyłączenie testowego modułu statystyk.

Jak połączyć ITK i Python?

Ustawiamy zmienną środowiskową PYTHONPATH na odpowiedni katalog:

export PYTHONPATH=/opt/lib/python2.6/site-packages

Jak podłączyć ITK i TCL?

Ustawiamy następujące zmienne systemowe:

export TCLLIBPATH=/opt/lib/InsightToolkit/WrapITK/Tcl
export LD_LIBRARY_PATH=/opt/lib/InsightToolkit

Testujemy połączenie z TCL:

tclsh
%puts $auto_path
%package require InsightToolkit
3.16.0

Automatyczne podłączenie TCL i Python przy starcie systemu
Aby nie ustawiać ręcznie zmiennych systemowych po każdym włączeniu komputera, możemy je dopisać do pliku ~/.bashrc, a system sam zrobi to za nas.
Poniżej przykład wpisu uwzględniający również VTK:

nano ~/.bashrc

export TCLLIBPATH="/opt/lib/vtk-5.4 /opt/lib/InsightToolkit/WrapITK/Tcl"
export LD_LIBRARY_PATH=/opt/lib/vtk-5.4:/opt/lib/InsightToolkit
export PYTHONPATH=/opt/lib/python2.6/site-packages

Instalacja ITK Application

Pakiet InsightToolkit sam w sobie zawiera tylko proste przykłady pokazujące jak pracować z biblioteką. Zbiór bardziej zaawansowanych przykładów jest dostępny w pakiecie InsightApplication (opis tych programów jest dostępny na stronie projektu, w dziale Applications).

Zdecydowana większość aplikacji korzysta z biblioteki FLTK do obsługi GUI. Możemy ją zainstalować następującym poleceniem:

sudo apt-get install libfltk1.1 libfltk1.1-dev fluid fltk1.1-doc

Następnie uruchamiamy CMake i ustawiamy, które aplikacje powinny zostać przygotowane. Ja wybrałem zdecydowaną większość z nich, z wyłączeniem wtyczek do Volview (pobierając Volview wtyczki są już skompilowane).

BUILD_SHARED_LIBS        OFF
BUILD_TESTING            OFF
CMAKE_BUILD_TYPE         Release
CMAKE_INSTALL_PREFIX     /opt
ITK_DIR                  /opt/lib/InsightToolkit
USE_FLTK                 ON
USE_VTK                  ON
VTK_DIR                  /opt/lib/vtk-5.4
USE_VolviewPlugins       OFF
USE_AUXILIARY            ON

Instalacja biblioteki InsightToolkit ze źródeł daje nam pewność, że korzystamy z jej najnowszej wersji, a możliwość wyboru parametrów kompilacji pozwala na dopasowanie jej do aktualnych potrzeb.

Jak skompilować VTK ze źródeł i podłączyć do TCLa i Pythona?

W Ubuntu 9.04, bo tego systemu będzie dotyczył ten wpis, VTK można łatwo zainstalować z repozytorium wydając następujące polecenie w konsoli:

sudo apt-get install libvtk5 libvtk5-dev libvtk5-qt4 libvtk5-qt4-dev tk8.4-dev vtkdata vtk-doc vtk-examples vtk-tcl python-vtk

Dostaniemy w ten sposób nie tylko biblioteki i pliki nagłówkowe VTK, ale również całą dokumentację, przykłady oraz dane testowe. Dodatkowo zostanie skonfigurowane połączenie z językami skryptowymi TCL i Python. Ten sposób instalacji w zdecydowanej większości przypadków powinien być wystarczający do pracy z tą biblioteką. Jednak co w przypadku jeśli potrzebujemy najnowszej wersji biblioteki (w chwili pisania tego wpisu w repozytorium była wersja 5.0.4, a na stronie oficjalnej wersja 5.4.2) lub niestandardowych parametrów pracy. Wtedy najlepiej samemu skompilować VTK ze źródeł.

Źródła VTK pobieramy ze strony oficjalnej projektu, z działu Download. Pobieramy również przykładowe dane. Po rozpakowaniu obu archiwów w wcześniej przygotowanym katalogu, przechodzimy do katalogu ze źródłami, tworzymy tam katalog build i przechodzimy do niego. Następnie uruchamiamy graficzny interfejs CMake i wybieramy opcje kompilacji (moja propozycja znajduje się poniżej):

mkdir build && cd build && cmake-gui ../

BUILD_EXAMPLES           ON
BUILD_SHARED_LIBS        ON
BUILD_TESTING            ON
CMAKE_BUILD_TYPE         Release
CMAKE_INSTALL_PREFIX     /opt
VTK_DATA_ROOT            /home/rafal/install/vtk/VTKData5.4.2/Data
VTK_USE_PARALLEL         ON
VTK_USE_GUISUPPORT       ON
VTK_USE_QVTK             ON
DESIRED_QT_VERSION       4
VTK_WRAP_JAVA            OFF
VTK_WRAP_PYTHON          ON
VTK_WRAP_TCL             ON

Po wygenerowaniu plików kompilacji w konsoli uruchamiamy program make i możemy spokojnie pójść na kawę:) - kompilacja może potrwać nawet pół godziny.

make
sudo make install

Komenda make install zainstaluje nam niezbędne pliki w katalogu /opt.

Uwagi po kompilacji:

• Dodatkowe biblioteki znalazły się również w katalogu /opt/lib/python2.6/site-packages/vtk/. Można je przekopiować do katalogu /opt/lib/vtk-5.4/. Operacja ta może w przyszłości ułatwić kompilację innych pakietów zależnych od VTK, np. InsightApplication.

  sudo cp /opt/lib/python2.6/site-packages/vtk/lib* /opt/lib/vtk-5.4/

Jak połączyć VTK i Python?
Wystarczy ustawić zmienną systemową PYTHONPATH na odpowiedni katalog:

export PYTHONPATH=/opt/lib/python2.6/site-packages

Jak podłączyć VTK i TCL?
Tutaj podobnie, należy ustawić odpowiednie zmienne systemowe:

export TCLLIBPATH=/opt/lib/vtk-5.4
export LD_LIBRARY_PATH=/opt/lib/vtk-5.4

Testujemy połączenie z TCL:

tclsh
%puts $auto_path
%package require vtk
5.4

Automatyczne podłączenie TCL i Python przy starcie systemu
Dopisujemy do pliku ~/.bashrc wcześniej opisywane zmienne:

nano ~/.bashrc

export PYTHONPATH=/opt/lib/python2.6/site-packages
export TCLLIBPATH=/opt/lib/vtk-5.4
export LD_LIBRARY_PATH=/opt/lib/vtk-5.4

Jak podłączyć się do przykładowych danych (VTK Data)?
Jeśli korzystamy ze ściągniętych danych przykładowych (VTK Data) i rozpakowaliśmy je do własnego katalogu, skrypty TCL mogą ich nie widzieć nawet jeśli ustawimy zmienne systemowe w bashu. Musimy w takim wypadku dopisać odpowiednią zmienną systemową w pliku /etc/environment:

sudo nano /etc/environment

VTK_DATA_ROOT="/home/rafal/install/vtk/VTKData5.4.2"

Jak przetestować czy wszystko działa poprawnie?

Najlepiej skorzystać z przykładów dołączonych do VTK. W katalogu Examples/Tutorial znajdują się proste przykłady zaimplementowane w kilku językach programowania.