seventh code

% Obtain sinograms for each slice in the projection
sinograms = zeros(size(image_stack, 1), size(image_stack, 2), num_slices);
for i = 1:num_slices
    sinograms(:, :, i) = radon(image_stack(:, :, i), angles);
end

% Reconstruct images using iradon
reconstructions_iradon = zeros(size(image_stack));
for i = 1:num_slices
    reconstructions_iradon(:, :, i) = iradon(sinograms(:, :, i), angles, 'linear', 'none', 1);
end

% Reconstruct images using ifanbeam
reconstructions_ifanbeam = zeros(size(image_stack));
for i = 1:num_slices
    reconstructions_ifanbeam(:, :, i) = ifanbeam(sinograms(:, :, i), angles, 'FanSensorGeometry', 'line', 'none', 1);
end

% Display images
slice_indices = [10, 20, 30]; % choose some representative slices
for i = 1:length(slice_indices)
    slice_index = slice_indices(i);
    figure;
    subplot(1, 3, 1);
    imshow(image_stack(:, :, slice_index), []);
    title(sprintf('Original slice %d', slice_index));
    subplot(1, 3, 2);
    imshow(reconstructions_iradon(:, :, slice_index), []);
    title(sprintf('Reconstructed slice %d (iradon)', slice_index));
    subplot(1, 3, 3);
    imshow(reconstructions_ifanbeam(:, :, slice_index), []);
    title(sprintf('Reconstructed slice %d (ifanbeam)', slice_index));
end

% Simulate non-idealised conditions
sinograms_noisy = sinograms .* (0.9 + rand(size(sinograms)) * 0.2);

% Reconstruct noisy images using iradon
reconstructions_iradon_noisy = zeros(size(image_stack));
for i = 1:num_slices
    reconstructions_iradon_noisy(:, :, i) = iradon(sinograms_noisy(:, :, i), angles, 'linear', 'none', 1);
end

% Reconstruct noisy images using ifanbeam
reconstructions_ifanbeam_noisy = zeros(size(image_stack));
for i = 1:num_slices
    reconstructions_ifanbeam_noisy(:, :, i) = ifanbeam(sinograms_noisy(:, :, i), angles, 'FanSensorGeometry', 'line', 'none', 1);
end

% Display noisy images
for i = 1:length(slice_indices)
    slice_index = slice_indices(i);
    figure;
    subplot(1, 3, 1);
    imshow(image_stack(:, :, slice_index), []);
    title(sprintf('Original slice %d', slice_index));
    subplot(1, 3, 2);
    imshow(reconstructions_iradon_noisy(:, :, slice_index), []);
    title(sprintf('Noisy reconstructed slice %d (iradon)', slice_index));
    subplot(1, 3, 3);
    imshow(reconstructions_ifanbeam_noisy(:, :, slice_index), []);
    title(sprintf('Noisy reconstructed slice %d (ifanbeam)', slice_index));
end

% Remove filter
sinograms_unfiltered = fftshift(ifft(ifftshift(sinograms), [], 1), 1); % inverse FFT along rows
reconstructions_unfiltered = zeros(size(image_stack));
for i = 1:num_slices
    reconstruction_unfiltered = fftshift(ifft(ifftshift(sinograms_unfiltered(:,