% 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 2015
% For some reason we need to clear all at this point
clear all
% Construct target image at desired resolution
N = 32;
target = phantom('Modified Shepp-Logan',N);
% Choose relative noise level in simulated noisy data
noiselevel = 0.01;
% Choose the angles between 0 and 180 degrees
T = 20;
th = ([0:(T-1)]/T)*180;
% 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(target,th);
% Construct phantom "target2" with higher resolution.
N2 = 2*N;
target2 = phantom('Modified Shepp-Logan',N2);
% 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 N/N2. For extra avoidance
% of inverse crime we add random errors to the measurement angles.
%[m2,tmp] = radon(target2,th);
[m2,tmp] = radon(target2,th+0.001*180*randn(size(th)));
ratio = N/N2;
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 "target"
% and "target2" matrices by constructing matching-scale coordinates for
% both images.
orind = floor((size(target)+1)/2); % Pixel containing the origin in target
orind2 = floor((size(target2)+1)/2); % Pixel containing the origin in target2
x = .5 + [0:(N-1)];
[X,Y] = meshgrid(x,x);
x2 = .5 + [0:(N2-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:T
mnc(:,iii) = interp1(s2(:)+odist*cos(2*pi*(th(iii)+45)/360),m2(:,iii),s(:),'spline');
end
% Correct for magnitude
mnc = mnc*ratio;
% Remove possible NaN (not-a-number) elements resulting from possible
% interpolation outside the domain of definition
mnc(isnan(mnc)) = 0;
% Monitor the accuracy of the interpolated data. When N 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)]);
% Reset the random number generator to ensure repeatable results
reset(RandStream.getGlobalStream);
% Construct noisy data
mncn = mnc + noiselevel*max(abs(mnc(:)))*randn(size(mnc));
% Save the result to file (with filename containing the resolution N)
savecommand = ['save data/XRMC_NoCrime_', num2str(N),'_',num2str(T), ' N N2 m mnc mncn th target target2 err_sup err_squ'];
eval(savecommand)
% View the results.
XR04_NoCrimeData_plot(N,T)