% cviceni 5, ukol 1a - simulace ultrazvukoveho snimku sondou pracujici v
% linearnim modu.
%
% Martin Barva, 23.4.2006, Lyon

% parametry systemu
c = 1540; % rychlst zvuku v prostredi [m/s]
f0 = 3e6; % centralni frekvence sondy [Hz]
fs = 50e6; % vzorkovaci frekvence systemu [Hz]

% parametry sondy
nElements = 64; % pocet piez. elementu
kerf = 0.1 / 1000; % mezera mezi elementy [m]
width = 1 / 1000; % sirka elementu [m]
height = 5 / 1000; % vyska elementu [m]
focalPoint = [0 0 55] / 1000; % defaultni pozice ohniska [m]
widthSubdivision = 1; % deleni elementu (kvuli presnosti vypoctu)
heightSubdivision = 1; % deleni elementu (kvuli presnosti vypoctu)

% nastav vzorkovaci frekvenci systemu
set_sampling(fs);

% vytvor linearni pole elementu; prijimac / vysilac maji stejne parametry
hndEmitter = xdc_linear_array(nElements, width, height, kerf, widthSubdivision, heightSubdivision, focalPoint);
hndReceiver = xdc_linear_array(nElements, width, height, kerf, widthSubdivision, heightSubdivision, focalPoint);

% tvar impulzni odezvy systemu a budiciho signalu
impulseResponse = sin(2 * pi * f0 * (0 : 1 / fs : 2 / f0));
impulseResponse = impulseResponse .* hanning(length(impulseResponse))';
excitationSignal = sin(2 * pi * f0 * (0 : 1 / fs : 2 / f0));

% nastav impulzni odezvu (vysilace / prijimace) a budici signal (vysilace)
xdc_impulse(hndEmitter, impulseResponse);
xdc_impulse(hndReceiver, impulseResponse);
xdc_excitation(hndEmitter, excitationSignal);

% pole souradnic difuzeru a jejich koeficientu odrazu. Souradnice jsou
% ulozeny v matici typu Nx3, kde N je pocet difuzeru; v 1,2, resp. 3 tim
% sloupci jsou x-ove, y-ove, resp. z-ove souradnice difuzeru. Koeficienty
% odrazu jsou ulozeny ve sloupcovem vektoru Nx1; hodnoty amplitud odrazu
% jsou v int. <0,1>; 0 - zadny odraz, 1 - maximalni odraz. 
%   V nasem pripade je phantom tvoreny jedinym difuzerem na pozici [0 0 55] s
%   koeficientem odrazu 1.
phntPoints = [0 0 55] / 1000; % souradnice bodu v [m]
phntCoefs = [1]; % koef. odrazu v int <0,1>

% sonda s nastavenymi parametry bude skenovat oblast serii parallelnich
% paprsku ekvidistatne umistenych v rovine xz. Velikost takto skenovane
% oblasti je imageXDim, pocet paprsku je nBeams.
imageXDim = 20 / 1000; % velikost snimku v x-ovem smeru [m]
nBeams = 80; % pocet paprsku skenujici tuto oblast
xBeamCoord = -imageXDim / 2;
dxBeamCoord = imageXDim / nBeams;

% v kazde iteraci je provedena simluace prijateho echa (rf signalu) podel jednoho
% paprsku.
times = zeros(nBeams, 1);
rfSignals = cell(nBeams, 1);
h = waitbar(0,'Please wait... (calculating rf signals)');
for i = 1 : nBeams,
    % pozice paprsku se provadi prikazy xdc_center(Th, P1) a xdc_focus(Th,
    % [0], P2). Zvuk se siri podel poloprimky s pocatecnim bodem P1
    % a prochazejici bodem P2. Navic v bode P2 bude paprsek zaostren.
    % Nastaveni se provadi pro vysilac i prijimac.
    xdc_center_focus(hndEmitter, [xBeamCoord 0 0]);
    xdc_focus(hndEmitter, [0], [xBeamCoord 0 55/1000]);
    
    xdc_center_focus(hndReceiver, [xBeamCoord 0 0]);
    xdc_focus(hndReceiver, [0], [xBeamCoord 0 55/1000]);
    
    % simuluj prijaty rf signal pro nastaveneny paprsek. V promenne rfSignals{i,1} jsou
    % vzorky rf signalu, v times(i) je ulozen casovy okamzik odpovidajici prvnimu
    % vzorku.
    [rfSignals{i, 1}, times(i)] = calc_scat(hndEmitter, hndReceiver, phntPoints, phntCoefs);
    
    xBeamCoord = xBeamCoord + dxBeamCoord;
    
    waitbar(i / nBeams, h)
end
close(h)

% zpracovani jednotlivych rf signalu (detekce obalky, log komprese). 
nSamples = zeros(nBeams, 1);
h = waitbar(0,'Please wait... (processing rf signals)');
for i = 1 : nBeams,
    % detekce obalky pomoci hilbertovy transformace
    rfSignals{i,1} = abs(hilbert(rfSignals{i,1}));
    
    % zjisti pocet vzorku. Kazdy rf signal ma jiny pocet samplu.
    nSamples(i) = length(rfSignals{i,1});
    
    waitbar(i / nBeams, h)
end
close(h)

% casovy okamzik nejdrive prijateho samplu pres vsechny rf signaly
minT = min(times);

% zjisti pocet vzorku nejdelsiho rf signalu
maxSamples = max(round(times * fs) + nSamples);

% nasyp data do matice image typu maxSamples x nBeams. Kazdy rf signal
% dopln odpovidajicim poctem nul at zacina v case minT.
image = zeros(maxSamples, nBeams);
for i = 1 : nBeams,
    nPadZeros = round(times(i) * fs);
    image(1:(nPadZeros + nSamples(i)), i) = [zeros(nPadZeros, 1); rfSignals{i, 1}];
end

% dynamika hodnot ulozenych v image je prilis velka pro zobrazeni (na
% snimku by nebyly videt detaily) -> pouziva se logaritmicka komprese
image = image / max(image(:));
image = log(image + 0.01);
image = image - min(image(:));
image = image / max(image(:));

% vypocet souradni voxelu. V x-ovem smeru je pozice pixlu dana souradnici
% paprsku, na kterem vzorek lezi. V z-ovem smeru se souradnice urci, ze
% znalosti casoveho okamziku pixlu a rychlosti sireni zvuku v prostredi.
xAxis = (-(imageXDim / 2) + [0:nBeams-1] * dxBeamCoord) * 1000;
zAxis = (([0:maxSamples-1] / fs) * c) / 2 * 1000;

% vykresleni vysledku
figure; imagesc(xAxis, zAxis, image); colormap(gray); colorbar('vert');
xlabel('Lateral direction [mm]'); ylabel('Axial direction [mm]'); title('Simulated ultrasound image');


