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

Delete ex2_k10.m

parent 4ad46778
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 266 deletions
ex2_k10.m deleted 100644 → 0
+ 0
− 266
View file @ 4ad46778
clear all
% close all
%% PARAMETERS %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
kappa = 10; % 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 = 2; % 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 = 10; % 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_k10 -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_k10 -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