Skip to content
Snippets Groups Projects
Commit f5353d97 authored by Roland Griesmaier's avatar Roland Griesmaier
Browse files

Upload New File

parent bce73bc5
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 266 additions and 0 deletions
ex2.m 0 → 100644
+ 266
− 0
View file @ f5353d97
clear all
close all
%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kappa = 1; % wave number
noiselevel = 0.0; % noise level of forward data (relative error)
nobs = kappa*64;
subsampling = 2; % increase the number of discretization points for the Nystrom method if nobs is too low
nobs_factor = 4; % size of extended sampling grid for the far field pattern
nobs_extended = nobs_factor*nobs;
model_error = 1e-3; % parameter in get_poles_via_aaa ( relative error of large argument asymptotics )
eta = 1e-6; % parameter in get_N_from_data
rebinning = 'gaussian';
epsilon = 1.3; % parameter in rebinning algorithm
% rebinning = 'nearest';
% epsilon = 0.8; % parameter in rebinning algorithm
axis_length = 16; % size of region of interest
b = 3*epsilon; % band-width of the ramp filter in filtered backprojection
RadonDeltaS = pi/b;
RadonDeltaOmega = 360 / (ceil(360 / (sqrt(2)/axis_length*pi/b)));
plot_flag = 0;
%% INITIALIZE GRID FOR VIRTUAL ORIGIN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
dist_Sensors = 25:1:29; % circular sampling grid
othetagrid = (0:2*pi/180:2*pi);
othetagrid = othetagrid(1:end-1);
origins = [];
for l = 1:length(dist_Sensors)
origins = [origins , dist_Sensors(l)*exp(1i*othetagrid)];
end
clear othetagrid;
nr_of_origins = length(origins);
%% SIMULATE SCATTERING DATA (OBSTACLES) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% [F,geometry] = computefarfield_obstacle('one_tiny_obstacle',kappa,subsampling*nobs,0);
[F,geometry] = computefarfield_obstacle('two_tiny_obstacles',kappa,subsampling*nobs,0);
% [F,geometry] = computefarfield_obstacle('three_tiny_obstacles',kappa,subsampling*nobs,0);
%% SAME FOR OBSTACLES AND INHOMOGENEOUS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F = F(1:subsampling:end,1:subsampling:end);
nr_of_scatterers = size(geometry,2)/2;
farfield = F(:,1); % far field data for plane wave with direction of propagation d=[1,0]
clear F Fcomponents;
farfield = addnoise(farfield, noiselevel); % add noise
%% EVALUATE EXTENDED FAR FIELD %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fkoeffs_extended = fft(farfield)/length(farfield); % Fourier coefficients of far field pattern
Nstar = get_N_from_data(fkoeffs_extended, eta);
Rmax = 2/exp(1)*Nstar/kappa;
fkoeffs_extended(Nstar+2:end-Nstar) = 0; % cut off tails
farfield_extended = ifft(fkoeffs_extended)*length(farfield);
norm(farfield - farfield_extended) / norm(farfield)
fkoeffs_extended = [fkoeffs_extended(1:length(farfield_extended)/2);zeros((nobs_factor-1)*length(farfield_extended),1);fkoeffs_extended(length(farfield_extended)/2+1:end)];
farfield = ifft(fkoeffs_extended)*length(fkoeffs_extended);
theta_farfield = (0:nobs_extended-1)/nobs_extended*2*pi;
xhat_farfield = exp(1i*theta_farfield); % observation points of far field pattern on S1
clear Nstar fkoeffs_extended farfield_extended
%% VISUALIZE GEOMETRY AND VIRTUAL ORIGINS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
clf
hold on
for iter = 1 : nr_of_scatterers
plot(geometry(:,2*iter-1), geometry(:,2*iter),'b-','LineWidth',2);
end
axis([-32 32 -32 32])
axis square
grid on
for iter=1:nr_of_origins
plot(origins(iter)+eps*1i,'r+','LineWidth',0.5);
end
plot([-axis_length axis_length axis_length -axis_length -axis_length],[axis_length axis_length -axis_length -axis_length, axis_length],'k--','LineWidth',2);
hold off
set(gca,'Fontsize', 15)
set(gca,'LooseInset',get(gca,'TightInset'));
print plots/ex2_geom -depsc
%% INITIALIZE SAMPLING GRID FOR RADON DATA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
RadonRotAngles = (0:RadonDeltaOmega:360);
RadonRotAngles = RadonRotAngles(1:end-1);
SensorPos_max = ceil(1.1*sqrt(2)*axis_length/RadonDeltaS)*RadonDeltaS;
RadonSensorPos = (-SensorPos_max:RadonDeltaS:SensorPos_max);
RadonData = zeros(length(RadonSensorPos),length(RadonRotAngles));
%% COMPUTE SINOGRAM %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
l = 1; % Laufindex fuer origin=0
min_poles = 0;
while (l <= nr_of_origins)
if (plot_flag==1) && (l>1)
for j=1:length(directions)
set(plot_directions(j),'color','k');
end
end
%% SHIFT ORIGIN %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
origin = origins(l);
farfieldshifted = exp(1i*kappa*real(xhat_farfield'*origin)) .* farfield;
fkoeffs_shifted = fft(farfieldshifted)/length(farfieldshifted); % Fourier coefficients of far field pattern
norm_fkoeffs = norm(fkoeffs_shifted);
%% ESTIMATE RADIUS OF SMALLES BALL CONTAINING ALL SCATTERERS %%%%%%%%%%
Nstar_shifted = get_N_from_data(fkoeffs_shifted, eta);
R_shifted = 2/exp(1)*Nstar_shifted/kappa;
if plot_flag == 1
figure(19)
hold on
plot(origin+eps*1i,'r*');
if l == 1
circle_plot = plot(origin+R_shifted*exp(1i*(0:1000)/1000*2*pi),'k:');
origin_plot = plot(origin,'ob');
else
set(circle_plot,'XData',real(origin + R_shifted*exp(1i*(0:1000)/1000*2*pi)),...
'YData',imag(origin + R_shifted*exp(1i*(0:1000)/1000*2*pi)))
set(origin_plot,'XData',real(origin),'YData',imag(origin))
end
end
M = floor(sqrt(kappa*(abs(origin)-Rmax)));
a = [fkoeffs_shifted(end-M+1:end) ; fkoeffs_shifted(1:M+1)];
delta = sqrt(pi*M^2/2) * model_error*norm(a);
[poles,res] = get_poles_via_aaa(a,min_poles,norm_fkoeffs,delta);
if plot_flag == 1
if l == 1
figure(19)
pole_plot = plot(origin+poles','r+','Markersize',10);
unit_circle_plot = plot(origin+exp(1i*(0:1000)/1000*2*pi),'k:');
else
figure(19)
set(pole_plot,'XData',real(origin+poles'),'YData',imag(origin+poles'))
set(unit_circle_plot,'XData',real(origin+exp(1i*(0:1000)/1000*2*pi)),...
'YData',imag(origin+exp(1i*(0:1000)/1000*2*pi)))
end
end
pole_count = length(poles);
directions = conj(poles);
for j=1:length(directions)
if plot_flag == 1
plot_directions(j) = plot(origin+[-directions(j),directions(j)]*5*R_shifted,'r');
end
dir_perp = 1i*directions(j)/abs(directions(j)); % Einheitsnormale auf directions(j)
dir_angle = angle(dir_perp)*180/pi;
if dir_angle < 0
dir_angle = dir_angle + 360;
end
dir_dist = real(origin.*conj(dir_perp));
[val_angle,idx_angle] = min(abs(dir_angle-RadonRotAngles));
switch rebinning
case 'gaussian'
RadonData_help = epsilon/sqrt(2*pi)*exp(-1/2*(RadonSensorPos'-dir_dist).^2*epsilon^2); % smooth with Gaussian in dist-direction
case 'nearest'
[val_dist,idx_dist] = min(abs(dir_dist-RadonSensorPos));
RadonData_help = zeros(length(RadonSensorPos),1);
RadonData_help(idx_dist,1) = 1;
end
% RadonData_help = abs(res(j)) * RadonData_help;
RadonData(:,idx_angle) = RadonData(:,idx_angle) + RadonData_help;
end
l = l+1;
end
if max(isnan(RadonData)) == 1
error('NaN in RadonData!')
end
figure(31)
imagesc(RadonRotAngles,RadonSensorPos,RadonData)
set(gca,'YDir','normal')
% colorbar
axis square
set(gca,'Fontsize', 15)
colormap(1-sqrt(gray))
print plots/ex2_sino -depsc
% RadonData_res = flipud( [RadonData(:,1),(RadonData(:,181:end)),RadonData(:,2:180)] );
% RadonData_res = (RadonData + RadonData_res);
%% FILTERED BACKPROJECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[I,xI] = fbp(RadonData,RadonRotAngles*pi/180,RadonSensorPos,axis_length,b);
% [I,xI] = fbp(RadonData_res,RadonRotAngles*pi/180,RadonSensorPos,axis_length,b);
I(I<0)=0;
figure(44)
imagesc(xI,xI,I)
set(gca,'YDir','normal')
axis square
set(gca,'Fontsize', 15)
% colorbar
colormap(flipud(hot))
print plots/ex2_recon -depsc
figure(45)
imagesc(xI,xI,I)
set(gca,'YDir','normal')
axis square
% colorbar
nr_of_obstacles = size(geometry,2) / 2;
hold on
for iter = 1:nr_of_obstacles
plot(geometry(:,2*iter-1), geometry(:,2*iter), '-w', 'linewidth', 2)
end
hold off
set(gca,'Fontsize', 15)
colormap(flipud(hot))
print plots/ex2_recon+geom -depsc
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Help | Imprint | Privacy policy | Accessibility | Contact