Search
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 <latex>I, J</latex> (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 <latex>T</latex>, která převede obrázek <latex>J</latex> 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ů.
histology.HEstain
histology.PanCytokeratin
T1.nii
FLAIR.nii
LesionSeg.nii
58_MRI.nii
64_MRI.nii
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é.
registration_images
cv3_matlabfun
fixed = imread('registration_images/brain.t1.png');
load_nii
image1 = load_nii('58_MRI.nii'); fixed = image1.img;
image1.img
figure, imshowpair( m_poly, fixed );
image(:,:,n)
imshow3D
imshow3D(moving, fixed, 'montage', 1);
montage|falsecolor|blend
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ů <latex> (x_i, y_i), i=1, \ldots N </latex> pro významné pozice, které lze identifikovat v obou obrázcích a pak hledáme transformaci <latex>T</latex>, která minimalizuje vzdálenost <latex> x_i </latex> a <latex> T(y_i) </latex>.
h = cpselect(moving, fixed);
fixedPoints
movingPoints
t_sim = fitgeotrans(movingPoints, fixedPoints, 'similarity');
similarity
moving
m_sim = imwarp(moving, t_sim, 'OutputView', imref2d(size(fixed)));
imref2d
fixed
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.
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
optimizer
metric
MutualInformation
registration.metric.MattesMutualInformation Properties: NumberOfSpatialSamples: 500 NumberOfHistogramBins: 50 UseAllPixels: 1
Properties
UseAllPixels = 0
NumberOfSpatialSamples = 50000
Objekt optimizer získáme následovně
opt_gd = registration.optimizer.RegularStepGradientDescent opt_oneplus = registration.optimizer.OnePlusOneEvolutionary
brain_T1.nii
brain_FLAIR.nii
niftyread
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)));
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
moving_aff
imregister
registration.metric.MeanSquare
registration.optimizer.RegularStepGradientDescent
imregdemons
[D, moving_def] = imregdemons( moving_aff, fixed, [400 150 60], 'AccumulatedFieldSmoothing', 1.8);
[400 150 60]
AccumulatedFieldSmoothing
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);