Cvičení 3 : Registrace

Na třetím cvičení se budeme zabývat registrací obrazů. Ta je potřebná ve chvíli, kdy máme dva obrázky I, J (např. rozdílné modality, rozdílné časy vytvoření snímku, apod.) a chceme informace z těchto dvou vstupů spojit dohromady. Cílem je tedy nalézt transformaci T, která převede obrázek J na obrázek I s co nejlepším možným překryvem.

Na cvičení budeme pracovat s matlabem a vyzkoušíme si jak manuální, tak i automatickou registraci různých obrázků.

Zadání

  • Landmark registration Pracujte s dvojicí obrázků histology.HEstain a histology.PanCytokeratin, které zobrazují skoro identický řez ve dvou různých barveních (Hematoxylin+Eosin a Cytokeratin). Určete dostatečný počet kontrolních bodů a nejděte co možno nejlepší transformaci. Vyzkoušejte alespoň dvě transformace, ukažte výsledek a diskutujte volbu transformace, rozmístění kontrolních bodů ve zprávě [1b]
  • Automatic registration Máme k dispozici dva MRI snímky (stejného pacienta) - T1.nii a FLAIR.nii, a dále segmentaci lézí na snímku FLAIR - obrázek LesionSeg.nii. Spočítejte použitím automatické registrace transformaci snímku FLAIR na snímek T1 a touto pak přetransformujte snímek segmentací. [2b] Spočtěte histogram intenzit snímku T1 přes segmentované oblasti. [BONUS +1b]
  • Deformable registration Zde máme dva MRI snímky kolene od různých pacientů - 58_MRI.nii a 64_MRI.nii. Spočítejte použitím automatické registrace transformaci snímku 64 na snímek 58. Protože se jedná o dvě rozdílné anatomie, musíme použít metodu pružné a konzistentní deformace, abychom dosáhli rozumného zarovnání. [2b]

Data, funkce

Vstupní data a pomocné funkce pro načítání a zobrazování jsou ke stažení jako ZIP Archiv. V archivu je adresář registration_images a cv3_matlabfun. Po rozbalení je potřeba přidat adresář cv3_matlabfun v Matlabu (add path), pak jsou funkce dostupné.

  • Načítání obrázků
    fixed = imread('registration_images/brain.t1.png');
  • Načítání 3D obrázků funguje za pomocí metody load_nii
    image1 = load_nii('58_MRI.nii');
    fixed = image1.img;
    Funkce vytvoří struct image1, pod image1.img se pak nachází obrazová data.
  • Zobrazení dvou obrázků imshowpair
    figure, imshowpair( m_poly, fixed );
    pro 3D obrázky si lze buď zobrazit n-tou vrstvu pomocí image(:,:,n) nebo použijte metodu imshow3D z přiloženého balíčku pomocných funkcí.
    imshow3D(moving, fixed, 'montage', 1);
    Parametry funkce imshow3D jsou dva obrázky, které chceme zobrazit, dále metoda překryvu montage|falsecolor|blend a jako poslední dimenze, kterou chceme zobrazit.

Landmark registration

Registrace pomocí kontrolních bodů je jednou ze základních poloautomatických metod. Začneme tím, že ručně určíme páry bodů  (x_i, y_i), i=1, \ldots N pro významné pozice, které lze identifikovat v obou obrázcích a pak hledáme transformaci T, která minimalizuje vzdálenost  x_i a  T(y_i) .

Postup

  1. Spustíme interaktivní rozhraní cpselect pro definici kontrolních bodů
    h = cpselect(moving, fixed);
  2. Před zavřením dialogového okna je potřeba uložit páry bodů do souborů, např fixedPoints a movingPoints
  3. S těmito množinami bodů pak pracujeme dále a hledáme optimální transformaci pomocí fitgeotrans
    t_sim = fitgeotrans(movingPoints, fixedPoints, 'similarity');
    Typ transformace (v příkladu similarity) je závislý na vstupních datech, možné další třídy transformací jsou v uvedené v dokumentaci
  4. Nakonec transformujeme moving image spočtenou transformací
    m_sim = imwarp(moving, t_sim, 'OutputView', imref2d(size(fixed)));
    Poslední argument imref2d specifikuje, že se vstupní obrázek má transformovat na mřížku o velikosti vstupního obrázku fixed.

Automatic registration

Příklad s MRI snímky hlavy je typickým problémem při vyhodnocování medicínským obrazů. K jednomu pacientovi máme dva různé snímky – MRI sekvence T1 a FLAIR. Na druhém jsou dobře vidět léze jako světlé oblasti a zajímají nás hodnoty druhém obrázku ve stejných oblastech. Kvůli pohnutí hlavy mezi snímáním jednotlivých sekvencí nejsou obrázky zarovnány a tudíž momentálně nesouhlasí segmentace a T1 snímek, jak je zřejmé z levého obrázku následující ilustrace.

 T1-Image FLAIR Image

Metoda

Automatická registrace v MATLABu je dostupná přes funkci imregister (vrací přetransformovaný obrázek) nebo imregtform (vrací transformaci). Abychom ji mohli spočíst, musíme specifikovat postup optimalizace přes parametr optimizer a dále pak vhodnou metriku, která měří chybu zarovnání – parameter metric. Jedná se o snímky s rozdílnými kontrasty, proto použijeme metriku MutualInformation

  registration.metric.MattesMutualInformation
 
  Properties:
    NumberOfSpatialSamples: 500
     NumberOfHistogramBins: 50
              UseAllPixels: 1
Nastavení parametrů uvedených pod Properties speficikují, jakým způsobem se metrika vyhodnocuje a má výrazný vliv na průběh a výsledek registrace. Výpočet lze urychlit, pokud nebudeme používat všechny pixely (UseAllPixels = 0), pak ale musíme počítat metriku přes dostatečné množství vzorků (NumberOfSpatialSamples = 50000 a více). Posledním parametrem před spuštěním registrace je nastavení třídy transformací. Jedná se o stejného pacienta, tudíž bychom mohli použít rigidní transformaci, lepší volbou je affiní transformace, která spolu s pohybem hlavy do určité míry kompenzuje rozdílná zkreslení v sekvencí i T1 a FLAIR.

Objekt optimizer získáme následovně

opt_gd = registration.optimizer.RegularStepGradientDescent
opt_oneplus = registration.optimizer.OnePlusOneEvolutionary

Postup

  1. Načteme soubory brain_T1.nii a brain_FLAIR.nii pomocí funkce load_nii obsažené v balíčku s pomocnými funkcemi, od verze R2017b lze použít příkaz niftyread
  2. Nastavíme objekty metric a optimizer a zavoláme funkci pro registraci obrázků
    moving_transformed = imregister(moving, fixed, 'rigid', optimizer, metric);
    % pokud počítáme transformaci, pak použijeme
    moving_t = imregtform(moving, fixed, 'rigid', optimizer, metric);
    moving_transformed = imwarp(moving, moving_t, imref3d(size(fixed)));
  3. Výslednou transformaci aplikujeme na FLAIR image, abychom zjistili kvalitu registrace, pokud je dobrá, pak můžeme transformaci aplikovat i na snímek se segmentací
  4. Histogram spočteme pomocí funkce imhist, nejprve ale musíme aplikovat transformovatnou masku

Deformable registration

Při registraci obrázků s rozdílnou anatomií není affiní transformace dostačující a je potřeba registrovaný obrázek deformovat. Tato metoda se používá ve chvíli, kdy chceme z množiny snímků vytvořit atlas (tj. mít obrázek, jak v průměru daná anatomická struktura vypadá), nebo když nás zajímá, jak se tkáň mezi dvěma snímky získaných s delším časovým odstupem změnila.

Abychom nemuseli počítat příliš velké deformace, spočteme v prvním kroku affiní transformaci (viz. předchozí úloha) a deformaci budeme určovat až mezi výsledkem affiní registrace a cílovým obrázkem

Postup

  1. Spočteme moving_aff pomocí funkce imregister. Snímky mají stejný rozsah intenzit, použijeme tedy metriku registration.metric.MeanSquare, optimalizovat budeme pomocí registration.optimizer.RegularStepGradientDescent
  2. Spočteme deformaci mezi moving_aff a fixed pomocí metody imregdemons
    [D, moving_def] = imregdemons( moving_aff, fixed, [400 150 60], 'AccumulatedFieldSmoothing', 1.8);
    Program pracuje s tzv. pyramidou - nejprve počítá deformaci na vstupech s redukovabným rozlišením, takto spočtenou deformaci postupně zjemňuje na úrovních pyramidy se rostoucím rozlišením, poslední úroveň je pak vstupní obrázek. Parameter [400 150 60] určuje, kolik iterací se počítá na jednotlivých úrovních pyramidy. Druhý parameter - AccumulatedFieldSmoothing regularizuje výslednou deformaci. Menší hodnoty umožňí detailnější deformace, vyšší hodnoty zvyšují hladkost výsledného deformačního pole.

Vizualizaci výsledné deformace můžeme získat aplikací deformačního pole na obrázek se 3D-šachovnicí s velikostí odpovídající vstupním obrázkům.

slice = checkerboard( 32, 4, 4 );
 
for i = 1:1:60
if mod( floor(i / 10), 2) > 0
    check(i, :, :) = uint8(255 * slice(:, :));
else
    check(i, :, :) = uint8( 255 * (1 - slice(:, :))); 
end
end
 
warped_checkerboard = imwarp( check, D );
 
% Zobrazíme deformovaný vstup a vizualizaci pole vedle sebe
imshow3d( moving_def, warped_checkerboard, 'montage', 1);

courses/a6m33zsl/lab03_registrace.txt · Last modified: 2018/03/05 17:34 by herinjan