## Contents

- Clear the workspace and close all open figures
- define variables
- read in the image
- Obtatain projections of the image using Matlab's radon function
- calculate the amount of Poisson noise for each ray in each projection
- Correct for infs in sinogram
- reconstruct the image using the noisy projections
- Do a noise only recon to fix problems with matlab (method 2)
- How to adjust No to get some level of "dose" reduction
- Picking an No value to simulate any dose lower than what you started with

%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 clc

## 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'); fclose(fid); %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; end % 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'); xlabel('1/No');ylabel('Variance');

## 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