1. No category

Ultraljudslabb
TSBB31 Medicinska Bilder
Ansvarig lärare: Mats Andersson
mats.x.andersson@liu.liu.se
Inehåll
•
•
•
•
•
•
•
Uppgiften
Läsa in RF-data
Enveloppdetektion
Subsampling
Histogram transformation
Skannkonvertering
Egen enveloppdetektion
Uppgiften
Uppgiften i denna labb är att skapa en Brightness-mode (B-mode) ultraljudsbild
från Radio frequency (RF)-data. RF-data är de ekon som detekteras i ultraljudsproben direkt efter att ett pulspaket har skickats ut och kristallerna i proben
har kopplats om till att lyssna efter ekon.
Exempel på en B-mode ultraljudsbild av ett hjärta.
RF-signalen har samma grundfrekvens som den utskickade pulsen och informationen som vi vill återskapa i bilden återfinns i hur signalens amplitud ändras
med tiden när ekon från olika djup i vävnaden detekteras i proben. Amplituden varierar mycket långsammare än RF-signalen men vi måste sampla signalen
med minst dubbla RF-frekvensen för att undvika vikning (samplingsteoremet).
Första steget är att detektera enveloppen på signalen. Vi kommer i denna labb
att använda oss av kvadraturfilter för att detta. Eftersom enveloppen varierar
mycket långsammare än grundtonen i RF-signalen kan vi nu subsampla utan
att förlora någon information.
Om vi tittar på bilden som den ser ut nu kommer vi att märka att vi tydligt
kan se ställen där vi har stor skillnad i akustisk impedans mellan intilliggande
1
vävnad men där skillnaden är liten är det mycket svårt att uppfatta gränsen
mellan olika vävnad. För att göra bilden lättare att tolka behöver vi förstärka
små skillnader i signalen. Vi gör en kompression av signalen. Denna del kallas
också för histogram transformation.
Sista steget mot en färdig bild är att göra en skannkonvertering, bilden är lagrad
i en matris där varje kolumn fysiskt motsvarar bildinnehållet i ett vinkelsegment.
För att få en geometriskt korrekt bild så måste vi göra en omsampling från
en rektangelformad bildarea till en solfjäderformad. Vi kommer att använda
matlabs interp2 för detta och vi vill göra omsamplingen så att vi vet exakt
vilken storlek i mm en pixel har så att vi enkelt kan mäta i bilden.
Sista uppgiften är att fundera ut ett eget sätt att göra enveloppdetektionen och
jämföra med de kvadraturfilter vi använde tidigare.
Läsa in RF-data
RF-data ligger i filen rfdata.mat som fås av lektionsassistenen eller laddas ner
från hemsidan.
load rfdata;
Varje kolumn i matrisen rf innehåller informationen från en skannstråle enligt
figuren nedan.
Ultraljudsprob
Lagrad datas
format
Stråle
Bildområde
Sampel
Utskickade
ultraljudspulser
Skannstrålarna är representerade som kolumner i matrisen rf.
I matlab kör
info
info.help
för att för reda på mer information om RF-data. Alla avstånd anges i meter och
vinklar i radianer. Vinklar och koordinat system finns också i figuren nedan.
Vinkel till den första skannstrålen är α0 och var och en av kolumnerna motsvarar
ett vinkelinkrement om deltaalpha
2
α0
X
r0
Y
R
Vinklar, avstånd och koordinatsystem för RF-data
Fråga: Vi ser att bilden börjar ett visst avstånd r0 från proben. Varför kan vi
inte avbilda delar av patienten som befinner sig godtyckligt nära proben?
Titta på rf som en gråbild med
imagesc(rf), title(’RF-data’)
colormap gray
RF−data
500
1000
1500
2000
2500
10
20
30
40
50
60
70
80
90
notera att imagesc automatiskt skalar om bilden, det mest negativa värdet blir
svart och det största vitt. rf har medelvärde 0 och detta motsvarar alltså grått i
bilden. Notera också storleken. Om vi visar rf med kvadratiska pixlar blir höjd
bredd förhållandet över 30ggr.
3
Envelopp detektion
Vi börjar med att titta på en enskild skannstråle (en kolumn i rf-matrisen)
col = 20;
h = rf(:,col); % h blir en vektor från
kolumn ’col’ i rf
Beräkna Fouriertransformen H av h mha fft. Tänkt på att använda fftshift
eftersom fft utgår från att origo liggre i början på vectorn.
H = .....
H = abs(H);
% amplitud spektrum
H = H/max(H); % normalisera map maxvärde
Plotta skannstrålen och dess amplitudspektrum
figure
subplot(211),plot(h),grid on
title(’RF-scanline amplitude and spectra’);
subplot(212)
plot(H), grid on
4
5
RF−scanline amplitude and spectra
x 10
0
−5
0
500
1000
1500
2000
2500
3000
0
500
1000
1500
2000
2500
3000
1
0.8
0.6
0.4
0.2
0
Använd horizontal zoom för att studera plotten av h.
I amplitud spektrat ovan har vi med både positiva och negativa frekvenser, origo
ligger alltså i mitten och Nyquistfrekvensen (som vi kallar π) i det samplade
systemet finns längst ut.
Fråga: Vilken fysisk frekvens motsvarar π? Hur stor är RF-frekvensen relativt
denna? Stämmer det med amplitud spektrat?
Fråga: I amblitudspektrat ser vi att det finns en liten topp på ungefär dubbla
RF-frekvensen. Vad beror det på?
Vi vill bara behålla den del av spektrat som motsvarar toppen runt RFfrekvensen, övriga delar kommer bara att bidra till ett försämrat signal brus
förhållande (SNR). Signalen behöver bandpass filtreras med ett bandpass filter
där center frekvensen u0 och bandbredden B är anpassad till signalens spektrum.
Vi kommer att använda kvadraturfilter både till att välja ut rätt frekvensband
och för att gör envelopp detektering. Först koncentrerar vi oss på frekvensbandet.
4
En mängd olika 1D kvadraturfilter finns i filen quadraturefilters.mat. Ett
kvadraturfilter med namn qfilt_pi_4_B_20 har en centerfrekvens u0 = π/4 och
en bandbredd B = 2.0 oktaver. En oktav innebär en fördubbling av frekvensen.
Bandbredden mäts från 6dB gränsen (dvs när filtrets amplitud har sjunkit till
hälften). För ett kvadraturfilter med bandbredd B = 2 innebär det att den
övre gränsfrekvensen, uh , är 4ggr så stor som den undre, ul . Vidare gäller att
√
u0 = ul uh och B = log 2(uh /ul ).
Läs in kvadraturfiltren och välj ut ett filter. Du kan räkna ut vilken u0 som böra
vara rätt annars får du prova dig fram.
load quadraturefilters.mat
f = qfilt_pi_...
f =
0.0001
0.0001
-0.0007
-0.0037
-0.0115
.
.
.
+
+
+
+
+
0.0005i
0.0026i
0.0070i
0.0140i
0.0225i
.
.
.
Som du ser är filtret komplexvärt. Plotta realdel och imaginärdel
figure
plot(real(f),’g’),hold on
plot(imag(f),’r’),grid on, hold off
title(’complex valued conv kernel f’)
complex valued conv kernel f
0.3
0.25
0.2
0.15
0.1
0.05
0
−0.05
−0.1
−0.15
−0.2
0
5
10
15
20
25
Fråga: Vad kan man säga om symmetrin för filtrets real och imaginärdel? Om
du tittar på flera filter upptäcker du att filter med lägre u0 är lite större. Vad
beror det på?
Plotta nu Fouriertransformen av filtret. Vi vill ha samma storlek på Fouriertransformen av filtret som som längden på rf för att kunna jämföra med signalens spektrum H. Vi paddar på nollor före och efter f så att vi får en vektor som
5
är lika lång som rf. Räkna ut det tal n som gör detta. Eftersom filtret alltid
har en udda storlek och storleken på rf är jämn så behöver vi en nolla mer på
den första delen
n = ...
fpad = [zeros(n,1); f; zeros(n-1,1)];
Beräkna Fouriertransformen av den reella delen av filtret först och plotta den
tillsammans med amplitud spektrum.
Fre = fftshift( fft( fftshift( real(fpad))));
Fre = real(Fre);
Fre borde vara reell men mycket små imaginära värden kan finnas. Kolla gärna!
Plotta Fre och amplitud spectrum i samma plot.
figure
plot(H,’g’),hold on
plot(Fre,’k’), hold off, grid on
title(’FT of real(f) and amplitude spectrum’)
Hitta en centerfrekvens och bandbredd som passar bra till amplitud spektrum.
Prova sen att titta på Fouriertransformen av imaginär delen i*imag(f) och hela
f och plotta resultatet;
Fim = ...
Fim = real(Fim);
F = ...
F = real(F);
figure
plot(Fre,’g’),hold on
plot(Fim,’r’)
plot(F,’k’), hold off, grid on
title(’FT of real(f), imag(f) and f’)
FT of real(f), imag(f) and f
2
1.5
1
0.5
0
−0.5
−1
0
500
1000
1500
2000
2500
3000
Fråga:Vad kan man säga om frekvensgången och symmetrin för den reella och
imaginära delen av f i Fourierdomänen, jämför med Fouriertransformen av sinus
och cosinus.
6
Nu har vi ett filter f som består av en jämn del som är reell i tidsdomänen
och en udda (och imaginär del) som har samma amplitud spektrum men släcker
ut varandra i ena halvan av Fourierdomänen och samverkar i den andra. Detta
är ett säkert kännetecken på ett kvadraturfilter! Nu ska vi titta på vad detta
innebär i tidsplanet. Filtrera (Falta) rf med f. Vi använder matlabs imfilter
för filtreringen.
q = imfilter(rf, f, ’replicate’);
Titta på resultatet för den reella delen och den imaginära delen i samma figur
för en skannstråle
figure
plot(real(q(:,col)),’g’); hold on
plot(imag(q(:,col)),’r’); hold off, grid on
title(’Real and imaginary part of f on rf’)
4
2
Real and imaginary part of f on rf
x 10
1.5
1
0.5
0
−0.5
−1
−1.5
−2
250
300
350
400
450
500
Zooma in i figuren, verkar det som om de två komponenterna är fasförskjutna
π/2?
Plotta slutligen rf och abs(q) i samma plot
q = abs(q);
figure
plot(rf(:,col)); hold on
plot(q(:,col),’r’);
plot(-q(:,col),’r’); hold off
title(’rf and abs(q)’)
4
2
rf and abs(q)
x 10
1.5
1
0.5
0
−0.5
−1
−1.5
−2
250
300
350
400
450
500
Om du har gjort rätt ska abs(q) nu följa amplituden på rf.
7
Titta nu på q som en gråbild, jämför med bilden av rf
figure
imagesc(q), title(’envelopp av rf’)
colormap gray
envelopp av rf
500
1000
1500
2000
2500
10
20
30
40
50
60
70
80
90
Subsampling
Som vi ser i plotten ovan så varierar enveloppen q mycket långsammare än
rf-signalen. Vi kan nu subsampla utan att förlora någon information. För att
avgöra hur mycket vi kan subsampla tittar vi på spektrum. Vi drar först bort
medelvärdet från en skann stråle i q. Medelvärdet innehåller ingen information
och dominerar spektrumet.
q0 = q(:,col) -mean(q(:,col)); % dra bort medelvärdet
Q = fftshift(abs(fft(q0)));
% amplitud spektrum av q0
Q = Q/max(Q);
figure
plot(Q),grid on
Vi bedömmer att när spektrum har sjunkit till 1/100 (40dB) kan vi subsampla
utan risk. Avgör hur mycket vi kan subsampla och använd resample Observera
att vi vill bara subsampla q utmed dess kolumner. Gör help resample för att
ta reda på hur det går till.
q = resample(q, ... , ...);
Titta på resultatet som en gråbild, ser du nån skillnad mot tidigare?
figure
imagesc(q), title(’subsamplad envelopp av rf’)
colormap gray
8
Histogram transformation
om vi tittar på bilden av q och dess histogram
figure
subplot(1,2,1),imagesc(q)
colormap gray,
title(’q’)
subplot(1,2,2)
hist(q(:),256),
title(’hist(q)’)
q
hist(q)
6000
100
5000
200
300
4000
400
3000
500
600
2000
700
800
1000
900
20
40
60
0
80
0
2
4
6
4
x 10
så ser vi att väldigt många pixlar är nästan svarta och det är bara i områden där
vi har stora skillnader i akustisk impedans som vi kan se något av strukturen.
Vi behöver förstärka små värden i bilden relativt de större för att få en mer
lätttolkad bild. Idealet är ett klockformat histogram liknande det i bilden nedan.
Att skriva q(:) är ett matlab trick som drar ut värdena i matrisen q till en
vektor. På så sätt får vi ett histogram för hela bilden, annars får vi ett histogram
för varje kolumn.
Ideal histogram
1000
900
800
700
600
500
400
300
200
100
0
2
3
4
5
6
7
8
9
10
11
För att åstakomma detta måste vi mappa q genom en funktion som förstärker de små värdena relativt de stora. Prova några olika t.ex qe=sqrt(q) eller
qe=log(q). Titta på bilden och på histogrammet.
När du är nöjd med bilden/histogrammet så justerar vi bilden så att värdena
ligger mellan 0 och 1 för att lätt kunna spara ner bilden på ett standardformat
senare.
qe = qe -min(qe(:));
qe = qe/max(qe(:));
figure
imagesc(qe),colormap gray
title(’Histogram equalized image gray values in range[0 , 1]’)
9
Skannkonvertering
Nu har vi en bild som ser bra ut utom i ett avseende. Bilden är insamlad i polära
koordinater, (vinkel, radie), och lagrad i en (kartesisk) matris. När vi tittar på
bilden av qe så har hjärtat fel geometri och storleken på objekt i bilden varierar
med positionen i bilden. Vi måste sampla om bilden till en solfjäder form. Vi
använder matlabs interp2 för detta. Vi kommer att behöva ange fem inargument som alla är matriser: I = interp2(qalpha, qrad, qe, Ialpha, Irad)
De tre första argumenten har med indata att göra och qe är vår histogramtranformerade bild. Matriserna qalpha och qrad har samma storlek som qe men
anger vilken vinkel resp. radie som varje position motsvarar. Titta gärna igen
på figuren i början på labpm som beskriver vinklar och koordinatsystem samt
samt gör info och info.help. Alla vinklar är i radiander och avstånd i meter.
Vi gör först två vektorer som beskriver hur alpha och r varierar i qe.
alpha0 = info.alpha0;
deltaalpha = info.deltaalpha;
nbeams = info.nbeams;
alphamax = alpha0 + deltaalpha*(nbeams-1);
alphavec = linspace(alpha0, alphamax, nbeams);
Matlabs linspace genererar en vektor som börjar i alphamax och slutar i
alpha0 med längden nbeams. Gör en motsvarande vektor för radien.
r0 = ...
R = ...
rvec = linspace(... , ..., length(qe));
Nu kan vi använda matlabs meshgrid för att generera matriserna qalpha och
qrad från vektorerna alphavec och rvec
[qalpha, qrad] = meshgrid(alphavec,rvec);
kontrollera att matriserna har samma storlek som qe och att det är regelbundna avstånd mellan rader/kolumner. Det är en förutsättning för att använda
interp2
De två sista argumenten i interp2 relaterar till de koordinater vi vill ha de nya
sampelpunkterna i. Vi börjar med en övning i geometri genom att räkna ut min
och max för X och Y koordinaterna (se figur i början på labpm).
Ymax
Ymin
Xmax
Xmin
=
=
=
=
R;
...
...
...
Bestäm storleken (i meter) på en pixel i den geometriskt korrekt bilden. Välj en
storlek på 0.5mm eller mindre.
dpix = ....
Gör en X och en Y vektor
10
Xvec = [Xmin:dpix:Xmax];
Yvec = [Ymin:dpix:Ymax];
Använd sedan meshgrid igen för att få motsvarande matriser med X och Y
koordinater i interpolerade bilden.
[X,Y] = meshgrid(Xvec,Yvec);
Matriserna X och Y innehåller de kartesiska koordinaterna för den färdiga bilden
I men för att vi ska kunna interpolera i originalbilden behöver vi räkna om X
och Y till polära koordinater.
Irad
=
Ialpha =
sqrt(X.^2 + Y.^2);
...
om du har gjort rätt ska Ialpha ha värden mellan π och 2π och vara störst på
höger sidan. Titta gärna på Irad och Ialpha som bilder om du är osäker. När
du är säker på att det är rätt så interpoleras den slutgiltiga bilden fram med
interp2.
I = interp2(qalpha, qrad, qe, Ialpha, Irad);
En del koordinater i I hamnar utanför qe. Matlab sätter dessa värden till NaN
(Not a Number). Vi sätter dessa värden till 0. Och tittar på resultatet.
I(isnan(I))=0;
figure
imagesc(I), colormap gray;
title(’Min färdiga B-mode bild’)
Min färdiga B−mode bild
50
100
150
200
250
300
350
400
450
50
100
150
200
250
300
350
400
450
Är du nöjd med resultatet? Gå gärna tillbaka och testa andra filter, mappningar,
pixelstorlekar mm.
Om du vill spara din bild i något standard format som t.ex PNG kan du använda
imwrite
imwrite(I,’myBmodeim.png’,’PNG’)
11
Egen enveloppdetektion
Sista delen av laborationen är att göra en egen envelopp detektion. Det räcker att
titta på en skannstråle h från rf. Ni kan t.ex prova med att hel eller halvvågs
likrikta h följt av ett lågpass filter eller kanske prova att demodulera? Några
matlab funktioner som kan vara till hjälp h = abs(h); (helvågslikriktning)
h = max(0,h); (halvvågs likriktning)
För att göra ett 1D LP-filter kan man använda fspecial
flp = fspecial(’gaussian’,[11,1],2)
flp =
0.0088
0.0271
0.0651
0.1216
...
...
...
som ger ett Gaussiskt LP-filter i y-led med storlek 11 och sigma=2. Lycka till!
12