%This script will add noise to CT images
% It does not model polychromatic energy spectrums, bowtie filters, mA
% modulation, detector electronic noise... There are many "CT simulators"
% that can do all of these things, but in reality, poeple in the clinic may
% have settled on a certain protocol and simply want to see what effect a
% little decrease in dose will have. I have not validated this yet with experimental data, too
% busy, I will soon, so it may totally suck.

Clear the workspace and close all open figures

clear all
close all

define variables

imLoc = '/home/szczykutowic/Desktop/testImagesAddNoise/series2reconFDK0984_0000_0983_000_xy.d';
No = 3.2315E6; % controls the level of noise added to the image
epsilon = 5; % corrects for negative vlaues in the pre logged sinogram data that would otherwise results in infinite projection values
voxelSize = 100/978; % in mm

read in the image

%This assumes the image is <> by <> and has data type = float
% if you can't fiugre out how to get your image in this format, change the
% code or change your image or email me
fid = fopen(imLoc,'rb');
imin = fread(fid,[1024 1024],'float');
%the image may need to be transformed back to attenuation coefficient
%first get images to HU unit --> these two transforms will vary depending
%on the format of the images you are starting with, most likely if you have
%a dicom image you will only need the secound of the following 2 transforms
imin = 2.317*imin - 1001;
%now get images to linear attenuation coef (1/mm)
imin = 0.0227*imin/1000 + 0.0227;%here, 0.0227 is the attenuation coef (1/mm) of water at 50 keV
%plot the image
figure;imagesc(imin',[0.022 0.026]);colormap('gray');axis off;axis equal;title('Original Image');

Obtatain projections of the image using Matlab's radon function

theta = linspace(0,360,900);
%bwelow we have to multiply the projections by the voxel size or else each
%voxel is given a "unit" wieght. This means for the same image simply
%changing the matrix size would change the projections values which is
%wrong, so we must account for this.
projs = radon(imin,theta)*voxelSize;
%plot the projections
figure;imagesc(projs,[0 max(max(projs))]);colormap('gray');axis off;axis equal;title('Original Sinogram');

calculate the amount of Poisson noise for each ray in each projection

N = No*exp(-projs);
N_noise = poissrnd(N);
projs_noise = -log(N_noise/No);
%plot the projections with noise added
figure;imagesc(projs_noise,[0 max(max(projs))]);colormap('gray');axis off;axis equal;title('Sinogram w/ noise added');

Correct for infs in sinogram

since sometimes N_noise is <= 0, projs_noise can be inf! We need to correct for this, I will do so here by finding all inf values and replacing them by -log(epsilon/No) where eplison is some small number >0 that reflects the smallest possible detected photon count

idx = isinf(projs_noise);
projs_noise(idx) = -log(epsilon/No);

reconstruct the image using the noisy projections

imin_noise = iradon(projs_noise,theta,'linear','Ram-Lak',1,size(imin,1))/voxelSize;
%show the image with noise added
figure;imagesc(imin_noise',[0.022 0.026]);colormap('gray');axis off;axis equal;title(['Image w/ noise No = ' int2str(No) ]);

Do a noise only recon to fix problems with matlab (method 2)

since matlab is using a pretty simple forward and backprojector, some edges may get messed up, to remedy this problem, we can do a recon of only the noise in the projection and then add that back to the original image

projs_noise_only = projs_noise - projs;
imin_noise2 = iradon(projs_noise_only,theta,'linear','Ram-Lak',1,size(imin,1))/voxelSize;
%plot the image with noise added
figure;imagesc(imin_noise2' + imin',[0.022 0.026]);colormap('gray');axis off;axis equal;title(['Image w/ noise No = ' int2str(No) ', 2nd method']);

%this method will appear to have more noise than the previous becasue in
%the process of forward projecting and reconstructing the image, a blurring
%effect is seen b/c of the non ideal forward and backprojects used in
%matlab (and b/c of the filtrations step in FBP), to see this more clearly, lets recon the original image from it's
%projections with no noise added and compare

imin_org_recon =  iradon(projs,theta,'linear','Ram-Lak',1,size(imin,1))/voxelSize;
%plot the original image reconstructed with the projections from radon
figure;imagesc(imin_org_recon',[0.022 0.026]);colormap('gray');axis off;axis equal;title(['Original image reconsteucted w/o noise added']);
%plot the difference between imin and imin_org_recon
figure;imagesc(imin'-imin_org_recon',[-0.004 0.005]);colormap('gray');axis off;axis equal;title(['Original image - Reconstructed original image']);
% you can see edges in the difference image, especially of the table.
% you can also see pretty easily how the original image got blurry

How to adjust No to get some level of "dose" reduction

this is not straight forward, we know that picking No = inf should add no noise to the projections and No ~ very small a lot of noise, the trouble is how do we connect No with the mAs/dose of the scan? Here we are assuming no mA modulation, which could be added by simply varying No with some sinosodial function, look in the literature

% if we assume var = var_anatomy + var_quantum + var_electronics and then
% assume we neglect var_anatomy b/c we choose a relativly flat ROI and
% var_electornics b/c we are not operating in a photon starved region of the
% detector, then var \approx 1/Dose. Now we know also that with these same
% assumptions var \approx 1/No. So if we want some decrease in dose, we can
% predict how the var should increase. To figure out what No to use to get
% this increase, first let us try a bunch of No and then plot var vs 1/No
% to get the relationship b/w the two. Then we can easily use this fit to
% solve for the No required for any level of dose reduction (so long as we
% dont go too low to where var_electronic comes into play).

% first pick an ROI to use to measure the variance in
%IMPORTANT when you are done choosing your ROI, you must right click inside
%of it and select create mask
figure;imagesc(imin);colormap('gray');axis off; axis equal;title('Choose your ROI in a flat region, then right click and sel. make mask');
[roi roiXIdx roiYIdx ] = roipoly;
BWprops = regionprops(roi, 'BoundingBox');

% now let us vary No and look at the variance inside our ROI!
indexx = 1E4:1E4:1E5;% you may need to change this for different size phantoms to be relevant (especially if your data is in 1/mm or 1/cm or 1/m or ... units)
for i=1:size(indexx,2)
% this code simply computes the image noise using method 2 from above
NoTest = indexx(i);
NTest = NoTest*exp(-projs);
N_noiseTest = poissrnd(NTest);
projs_noiseTest = -log(N_noiseTest/NoTest);
projs_noise_onlyTest = projs_noiseTest - projs;
% we have to correct for inf values again, like we did above
idx = isinf(projs_noise_onlyTest);
projs_noise_onlyTest(idx) = -log(epsilon/NoTest);

imin_noise2Test = imin + iradon(projs_noise_onlyTest,theta,'linear','Ram-Lak',1,size(imin,1))/voxelSize;

% all this code simply extracts the ROI you just chose and computes the
% variance within that ROI
iminCropped =  imin_noise2Test.* cast(roi, class(imin_noise2Test));
tmp = imcrop(iminCropped, BWprops.BoundingBox);
idx = tmp > 0;
tmp2 = tmp(idx);
varr(i) = std(tmp2(:))^2;


% now we need to fit var vs 1/No
p = polyfit(1./indexx,varr,1);
figure;plot(1./indexx,varr,'-',1./indexx,polyval(p,1./indexx),'--');legend('Data','Fit');title('Var vs. 1/No');

Picking an No value to simulate any dose lower than what you started with

% again, as low a dose as the assumption electronic noise is not a problem

% based on our fit from above, we can calculate the No corresponding to any
% dose we want since for some dose decrease d, we know the var will increase
% by 1/(1-d) (for example if we decrease dose by 10%, the var will increase
% by 10/9 or 1/(1-0.1))

% first figure out the starting noise variance
iminCropped =  imin.* cast(roi, class(imin));
tmp = imcrop(iminCropped, BWprops.BoundingBox);
idx = tmp > 0;
tmp2 = tmp(idx);
startingVar = std(tmp2(:))^2;

doseDecreaseDesired = 0.875;
desiredVar = 1/(1-doseDecreaseDesired) * startingVar;
NoToUse =  p(1)/(desiredVar - p(2));

% now you can use this NoToUse value in the code at the beginning of this
% script and it should give you results similar to actually reducing your
% mAs/Dose by doseDecreaseDesired