% Here we simulate X-ray tomographic measurement avoiding inverse crime. 
% This is done by simulating data first at higher resolution using Matlab's
% Radon.m routine and then interpolating the data to the desired (lower)
% resolution. This must be done carefully as radon.m chooses the origin of
% coordinates in a specific way, and moving between different resolutions
% involves a coordinate change.
%
% Samuli Siltanen February 2017

% For some reason we need to clear all at this point
clear all

%% Choose the relative amplitude of simulated measurement noise
noiselevel = 0.05;

%% Define the measurement setup at lower resolution
% Choose the size of the unknown. The image has size (2^j)x(2^j);
% let's denote n = (2^j)*(2^j).
j = 5;

% Choose the angles 
Nang = 33;

% Load precomputed results of the routine tomo01_RadonMatrix_comp.m
% with the particular values of j and Nang above
loadcommand = ['load data/RadonMatrix_', num2str(j), '_', num2str(Nang), ' M angles'];
eval(loadcommand)

% Construct target
%f = phantom('Modified Shepp-Logan',M);
f = SqPhantom(M);

% For comparison, construct inverse-crime measurement at the desired
% resolution. The vector s contains the linear coordinates
% of the Radon transform in pixel coordinates.
[m,s] = radon(f,angles);


%% Construct measurement data that avoids inverse crime. 
% The trick is to use higher resolution and downsample.

% Construct phantom "f2" with higher resolution. 
M2 = 2*M;
%f2 = phantom('Modified Shepp-Logan',M2);
f2 = SqPhantom(M2);

% Construct tomographic data at higher resolution. Vector tmp contains 
% the linear coordinates of the Radon transform in pixel coordinates. 
% However, since the two phantoms have different numbers of pixels, 
% we correct for the pixel size using the ratio M/M2. For extra avoidance
% of inverse crime we add random errors to the measurement angles.
%[m2,tmp]  = radon(f2,angles);
[m2,tmp]  = radon(f2,angles+0.001*360*randn(size(angles)));
ratio     = M/M2;
s2        = tmp*ratio;

% Since Matlab's radon.m routine uses different locations of the origin at
% different resolutions (type in Matlab >> help radon), we need to correct 
% for the displacement. We calculate distance between the origins in "f" 
% and "f2" matrices by constructing matching-scale coordinates for
% both images.
orind   = floor((size(f)+1)/2); % Pixel containing the origin in f
orind2  = floor((size(f2)+1)/2); % Pixel containing the origin in f2
x       = .5 + [0:(M-1)];
[X,Y]   = meshgrid(x,x);
x2      = .5 + [0:(M2-1)];
x2      = x2*ratio;
[X2,Y2] = meshgrid(x2,x2);
orx     = X(orind(1),orind(2));
ory     = Y(orind(1),orind(2));
orx2    = X2(orind2(1),orind2(2));
ory2    = Y2(orind2(1),orind2(2));
odist   = sqrt((orx-orx2)^2+(ory-ory2)^2);
% figure(13)
% clf
% plot(X,Y,'r.','markersize',10)
% hold on
% plot(X2,Y2,'b.','markersize',10)
% axis equal

% Construct measurement "mnc" from higher resolution data by interpolation. 
% Then "mnc" contains no inverse crime. Note the correction for
% displacement of the origin due to the two resolutions used.
for iii = 1:Nang
    mnc(:,iii) = interp1(s2(:)+odist*cos(2*pi*(angles(iii)+45)/360),m2(:,iii),s(:),'spline');
end

% Correct for magnitude
mnc = mnc*ratio;

% Remove NaN (not-a-number) elements resulting from possible 
% interpolation outside the domain of definition
mnc(isnan(mnc)) = 0;


%% Check the accuracy of the interpolated data. 
% When M grows, the errors become smaller. However, the error should not be 
% zero as a small positive error simulates the inevitable modeling errors 
% of practical situations.
err_sup = max(max(abs(m-mnc)))/max(max(abs(m)));
err_squ = norm(m(:)-mnc(:))/norm(m(:));
disp(['Sup norm relative error: ', num2str(err_sup)]);
disp(['Square norm relative error: ', num2str(err_squ)]);

%% Construct noisy data, save the results and plot the situation

% Noisy data without inverse crime
mncn = mnc + noiselevel*max(abs(mnc(:)))*randn(size(mnc));

% Save the result to file (with filename containing the resolution M)
savecommand = ['save data/tomo03_NoCrime_', num2str(M),'_',num2str(Nang), ' M M2 m mnc mncn angles Nang f f2 err_sup err_squ'];
eval(savecommand)

%% View the results. 
tomo03_NoCrimeData_plot(M,Nang)





