Add files via upload

MIPT-Oulu · Apr 3, 2023 · ac800bd · ac800bd
1 parent 5e193b4
commit ac800bd
Show file tree

Hide file tree

Showing 51 changed files with 1,420 additions and 0 deletions.
diff --git a/CTlab v.1.4/CTlab.v1 funktiot/FBP_2D_opTomo.m b/CTlab v.1.4/CTlab.v1 funktiot/FBP_2D_opTomo.m
@@ -0,0 +1,101 @@
+function rec = FBP_2D_opTomo(proj_data,W,FilterType)
+
+if nargin == 2
+    FilterType = 'ram-lak';
+end
+
+d = 1;
+[p,H] = filterProjections(proj_data', FilterType, d);
+p = p';
+rec = reshape(W'*p(:),W.vol_size);
+
+%% Subfunctions and modules
+
+function [p,H] = filterProjections(p_in, filter, d)
+
+p = p_in;
+
+% Design the filter
+len = size(p,1);
+H = designFilter(filter, len, d);
+
+if strcmpi(filter, 'none')
+    return;
+end
+
+p(length(H),1)=0;  % Zero pad projections
+
+% In the code below, I continuously reuse the array p so as to
+% save memory.  This makes it harder to read, but the comments
+% explain what is going on.
+
+p = fft(p);               % p holds fft of projections
+
+p = bsxfun(@times, p, H); % faster than for-loop
+
+p = ifft(p,'symmetric');  % p is the filtered projections
+
+p(len+1:end,:) = [];      % Truncate the filtered projections
+%----------------------------------------------------------------------
+
+%======================================================================
+function filt = designFilter(filter, len, d)
+% Returns the Fourier Transform of the filter which will be
+% used to filter the projections
+%
+% INPUT ARGS:   filter - either the string specifying the filter
+%               len    - the length of the projections
+%               d      - the fraction of frequencies below the nyquist
+%                        which we want to pass
+%
+% OUTPUT ARGS:  filt   - the filter to use on the projections
+
+
+order = max(64,2^nextpow2(2*len));
+
+if strcmpi(filter, 'none')
+    filt = ones(1, order);
+    return;
+end
+
+% First create a bandlimited ramp filter (Eqn. 61 Chapter 3, Kak and
+% Slaney) - go up to the next highest power of 2.
+
+n = 0:(order/2); % 'order' is always even. 
+filtImpResp = zeros(1,(order/2)+1); % 'filtImpResp' is the bandlimited ramp's impulse response (values for even n are 0)
+filtImpResp(1) = 1/4; % Set the DC term 
+filtImpResp(2:2:end) = -1./((pi*n(2:2:end)).^2); % Set the values for odd n
+filtImpResp = [filtImpResp filtImpResp(end-1:-1:2)]; 
+filt = 2*real(fft(filtImpResp)); 
+filt = filt(1:(order/2)+1);
+
+w = 2*pi*(0:size(filt,2)-1)/order;   % frequency axis up to Nyquist
+
+switch filter
+    case 'rect'
+        filt = filt.^0.5;
+    case 'custom'
+        L = 4.2;%4; % Sample averages
+        H = 1/L*(1-exp(-i*L*w))./(1-exp(-i*w)); % Frequency response of a moving average filter
+        H = abs(H);
+        % Multiply with Hann filter
+        H = H.*(1+cos(w./d)) / 2;
+        filt(2:end) = filt(2:end).*H(2:end);
+        %
+    case 'ram-lak'
+        % Do nothing
+    case 'shepp-logan'
+        % be careful not to divide by 0:
+        filt(2:end) = filt(2:end) .* (sin(w(2:end)/(2*d))./(w(2:end)/(2*d)));
+    case 'cosine'
+        filt(2:end) = filt(2:end) .* cos(w(2:end)/(2*d));
+    case 'hamming'
+        filt(2:end) = filt(2:end) .* (.54 + .46 * cos(w(2:end)/d));
+    case 'hann'
+        filt(2:end) = filt(2:end) .*(1+cos(w(2:end)./d)) / 2;
+    otherwise
+        error(message('images:iradon:invalidFilter'))
+end
+
+filt(w>pi*d) = 0;                      % Crop the frequency response
+filt = [filt' ; filt(end-1:-1:2)'];    % Symmetry of the filter
diff --git a/CTlab v.1.4/CTlab.v1 funktiot/calculateMasks.m b/CTlab v.1.4/CTlab.v1 funktiot/calculateMasks.m
@@ -0,0 +1,6 @@
+function [maskIodine,maskWater,maskAir,maskBone]=calculateMasks(p)
+maskIodine = p < 0.31 & p > 0.29;                   %Iodine mask
+maskWater = p < 0.21 & p > 0.19;                  %Water mask
+maskAir = p < 0.01 & p > -0.01;                   %Air mask
+maskBone = p < 1.01 & p > 0.99;
+end
diff --git a/CTlab v.1.4/CTlab.v1 funktiot/calculateProjections.m b/CTlab v.1.4/CTlab.v1 funktiot/calculateProjections.m
@@ -0,0 +1,60 @@
+function [p,W,I_logsum,I,maskIodine, maskWater, maskAir, maskBone] = calculateProjections(fig,minAngle,stepAngle,maxAngle,imageVolume,...
+                                                                     detectorElementsize,detectorWidth,I0,imageNoise,minE0,maxE0,time,...
+                                                                     geom, SOD, ODD,mA)
+
+d = uiprogressdlg(fig,'Title','Please wait',...
+    'Message','Calculation in progress...','Cancelable','on');
+drawnow
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%                      1.Attenuation data for different energies                       %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Example densities (g/cm^3) for different materials in Shepp-Logan phantom
+boneDensity = 0.96;
+waterDensity = 0.5;
+iodineDensity =0.037;
+airDensity = 0.001225;
+
+[attenuationMatrix,energyMatrix]=loadAttenuationData()
+p=loadPhantom(imageVolume);
+
+
+%Mass attenuation coefficients data of iodine, water, air & bone: NIST
+%(https://www.nist.gov/pml/x-ray-mass-attenuation-coefficients)
+%Import data from a spreadsheet
+
+[u_iodine_E, u_water_E, u_air_E, u_bone_E] = splineFit(energyMatrix,attenuationMatrix,iodineDensity,waterDensity,boneDensity,airDensity);
+
+% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%              2.Parameters for volume & projection geometry and also for noise:                 %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+W = systemMatrix(minAngle,stepAngle,maxAngle,geom,imageVolume,detectorElementsize,detectorWidth);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%                           3.Masks for different tissues and materials:                            %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+[maskIodine,maskWater,maskAir,maskBone]=calculateMasks(p);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%                                  4.Forward project basis materials:                               %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+[projIodine,projWater,projAir,projBone]=forwardprojectPhantomMaterials(W,maskIodine,maskWater,maskAir,maskBone);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%                   5.Forward projections (Beer-Lambert law) for different energies:                 %%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+I=forwardprojectBeerLambert(minE0,maxE0,u_air_E,u_water_E,u_iodine_E,u_bone_E,projAir,...
+                            projWater,projIodine,projBone,imageNoise,mA,time,I0,d);
+
+
+I_logsum = -log(sum(I,3)/sum(I0))/0.1;
+
+close(d)
+findall(0,'type','figure','tag','TMWWaitbar')
+
+end
diff --git a/CTlab v.1.4/CTlab.v1 funktiot/determineReconstruction.m b/CTlab v.1.4/CTlab.v1 funktiot/determineReconstruction.m
@@ -0,0 +1,7 @@
+function reco=determineReconstruction(reconstruction,i,ParametersinfoListBox,BeamgeometryButton)
+            reco=reconstruction;
+            t = datetime('Now')
+            t=datestr(t)
+            ParametersinfoListBox.Items(i) = {strcat(t, ' <Reconstruction algorithm: ',reco,'>')}
+            BeamgeometryButton.Enable ='on'
+end
diff --git a/CTlab v.1.4/CTlab.v1 funktiot/determineSpectrum.m b/CTlab v.1.4/CTlab.v1 funktiot/determineSpectrum.m
@@ -0,0 +1,39 @@
+function [kVp,mA,Cu,Al,E0,I0,maxE0,minE0] = determineSpectrum(q,tubeVoltage,tubeCurrent,copperFilter,aluminiumFilter,ax,recButtonState)
+I0 = q;
+
+%Energies in keV
+E0=1:1:length(I0);
+
+%Eliminate zero values from vector
+E0=E0(I0~=0);
+
+%Max/Min value of energy in the spectrum
+maxE0 = max(E0);
+minE0 = min(E0);
+
+%String which is displayed in text box
+stringforBox1='Min energy: %d keV';
+stringforBox2='Max energy: %d keV';
+stringforBox1 = sprintf(stringforBox1,minE0);
+stringforBox2 = sprintf(stringforBox2,maxE0);
+
+kVp=tubeVoltage;
+mA=tubeCurrent;
+Cu=copperFilter;
+Al=aluminiumFilter;
+
+%Plot spectrum in Spectrum figure
+plot(I0,'LineWidth',2, 'Parent',ax);
+
+limY = max(I0);
+ax.YLim = [0 limY];
+
+%Display min and max energies of the spectrum in text box
+text(ax,70,limY-limY/10,stringforBox1,'Color','w','FontSize',11);
+text(ax,70,limY-2*limY/10,stringforBox2,'Color','w','FontSize',11);
+
+ax.XTickMode = 'auto';
+ax.YTickMode = 'auto';
+recButtonState.Enable = 'on';
+
+end
diff --git a/CTlab v.1.4/CTlab.v1 funktiot/energySampling.m b/CTlab v.1.4/CTlab.v1 funktiot/energySampling.m
@@ -0,0 +1,21 @@
+function [energyWindows, numBins, xl] = energySampling(maxEnergy, minEnergy, bin, ax3)
+
+res=mod(maxEnergy-minEnergy, bin);
+
+energyWindows = minEnergy:bin:maxEnergy;
+
+if res ~= 0
+    energyWindows = [energyWindows max(energyWindows)+res'];
+end
+
+for k=1:length(energyWindows)
+    xl=xline(energyWindows(k),':','Parent',ax3);
+    xl.Color = 'w';
+    hold on
+end
+
+numBins = length(energyWindows);
+
+close
+
+end
diff --git a/CTlab v.1.4/CTlab.v1 funktiot/filteredBackprojection2.m b/CTlab v.1.4/CTlab.v1 funktiot/filteredBackprojection2.m
@@ -0,0 +1,58 @@
+function [rec,minI,maxI,h] = filteredBackprojection2(fig,text,I_logsum,mask_I,mask_water,mask_air,mask_luu,W,axis,axis2)
+
+% h = waitbar(0,'Please wait...');
+d = uiprogressdlg(fig,'Title','Please wait',...
+    'Message','Calculation in progress...','Cancelable','on');
+drawnow
+
+rec=FBP_2D_opTomo(I_logsum,W);
+rec1=mask_I.*rec;
+rec2=mask_water.*rec;
+rec3=mask_air.*rec;
+rec4=mask_luu.*rec;
+
+water=rec2(:);
+water=water(water~=0);
+water=abs(water);
+water = mean(water);
+
+
+% HU_iodine = 1000*((rec1-mean(water))/mean(water));
+% HU_iodine = HU_iodine.*mask_I;
+% HU_water = 1000*((rec2-mean(water))/mean(water));
+% HU_water = HU_water.*mask_water;
+% HU_air = 1000*((rec3-mean(water))/mean(water));
+% HU_air = HU_air.*mask_air;
+% HU_luu = 1000*((rec4-mean(water))/mean(water));
+% HU_luu = HU_luu.*mask_luu;
+
+% rec=HU_iodine+HU_water+HU_air+HU_luu;
+% rec=1000*((rec-mean(water))/mean(water));
+rec=1000*((rec-mean(water))/mean(water));
+%             imshow(rec,[], 'Parent', axis);
+drawnow
+pause(1)
+h = histogram(axis2, rec, 'FaceColor','Black','EdgeColor', 'none');
+% h.EdgeColor = 'r'
+h.FaceColor = [0.50,0.62,0.67];
+h.FaceAlpha = 0.4;
+% minI=min(h.BinEdges);
+% maxI=max(h.BinEdges);
+minI=h.BinLimits(1);
+maxI=h.BinLimits(2);
+I = imshow(rec,[minI-minI/2 maxI+maxI/2], 'Parent', axis, ...
+    'XData', [0 axis.Position(3)], ...
+    'YData', [0 axis.Position(4)]);
+% Set limits of axes
+axis.XLim = [0 I.XData(2)];
+axis.YLim = [0 I.YData(2)];
+%imshow(rec,[], 'Parent', axes);
+d.Value = 1;
+d.Message = sprintf(text);
+pause(1);
+close(d);
+display(minI);
+display(maxI);
+% F = findall(0,'type','figure','tag','TMWWaitbar')
+
+end
diff --git a/CTlab v.1.4/CTlab.v1 funktiot/forwardprojectBeerLambert.m b/CTlab v.1.4/CTlab.v1 funktiot/forwardprojectBeerLambert.m
@@ -0,0 +1,38 @@
+function I=forwardprojectBeerLambert(minE0,maxE0,u_air_E,u_water_E,u_iodine_E,u_bone_E,projAir,...
+                                                                   projWater,projIodine,projBone,imageNoise,mA,time,I0,d)
+
+N1=size(projAir,1);
+N2=size(projAir,2);
+I_air = zeros(N1,N2,maxE0);
+I_iodine = zeros(N1,N2,maxE0);
+I_water = zeros(N1,N2,maxE0);
+I_bone = zeros(N1,N2,maxE0);
+I = zeros(N1,N2,maxE0);
+noiseType='gaussian';                       %Noise type
+
+for i=minE0:1:maxE0 %length(I0)
+    %Beer-Lambert law for each energy of forward projection
+
+    %Add noise for each forward projection
+    I_air(:,:,i)=imnoise(exp(-u_air_E(i)*projAir*0.1),noiseType,0,imageNoise);
+    I_water(:,:,i)=imnoise(exp(-u_water_E(i)*projWater*0.1),noiseType,0,imageNoise);
+    I_iodine(:,:,i)=imnoise(exp(-u_iodine_E(i)*projIodine*0.1),noiseType,0,imageNoise);
+    I_bone(:,:,i)=imnoise(exp(-u_bone_E(i)*projBone*0.1),noiseType,0,imageNoise);
+
+    %Calculate forward projection which includes all materials, energies,
+    %and exposure (mAs = mA*time)
+    I(:,:,i) = mA*time*I0(i)*I_air(:,:,i).*I_iodine(:,:,i).*I_water(:,:,i).*I_bone(:,:,i);
+
+    % Check for Cancel button press
+    if d.CancelRequested
+        break
+    end
+
+    %Update progress, report current estimate
+    d.Value = i/length(I0);
+    d.Message = sprintf('Projections for energy %d keV processed',i);
+end
+
+
+
+end
diff --git a/CTlab v.1.4/CTlab.v1 funktiot/forwardprojectBeerLambert2.m b/CTlab v.1.4/CTlab.v1 funktiot/forwardprojectBeerLambert2.m
@@ -0,0 +1,38 @@
+function I=forwardprojectBeerLambert2(minE0,maxE0,u_air_E,u_water_E,u_iodine_E,u_bone_E,projAir,...
+                                                                   projWater,projIodine,projBone,imageNoise,mA,time,I0,d)
+
+N1=size(projAir,1);
+N2=size(projAir,2);
+I_air = zeros(N1,N2,maxE0);
+I_iodine = zeros(N1,N2,maxE0);
+I_water = zeros(N1,N2,maxE0);
+I_bone = zeros(N1,N2,maxE0);
+I = zeros(N1,N2,maxE0);
+
+
+for i=minE0:1:maxE0 %length(I0)
+    %Beer-Lambert law for each energy of forward projection
+
+    %Add noise for each forward projection
+    I_air(:,:,i)=exp(-u_air_E(i)*projAir*0.1);
+    I_water(:,:,i)=exp(-u_water_E(i)*projWater*0.1);
+    I_iodine(:,:,i)=exp(-u_iodine_E(i)*projIodine*0.1);
+    I_bone(:,:,i)=exp(-u_bone_E(i)*projBone*0.1);
+
+    %Calculate forward projection which includes all materials, energies,
+    %and exposure (mAs = mA*time)
+    I(:,:,i) = mA*time*I0(i)*I_air(:,:,i).*I_iodine(:,:,i).*I_water(:,:,i).*I_bone(:,:,i);
+
+    % Check for Cancel button press
+    if d.CancelRequested
+        break
+    end
+
+    %Update progress, report current estimate
+    d.Value = i/length(I0);
+    d.Message = sprintf('Projections for energy %d keV processed',i);
+end
+
+
+
+end
diff --git a/CTlab v.1.4/CTlab.v1 funktiot/forwardprojectPhantomMaterials.m b/CTlab v.1.4/CTlab.v1 funktiot/forwardprojectPhantomMaterials.m
@@ -0,0 +1,8 @@
+function [projIodine,projWater,projAir,projBone]=forwardprojectPhantomMaterials(W,maskIodine,maskWater,maskAir,maskBone)
+
+projWater = reshape(W*double(maskWater(:)),W.proj_size);        %Water forward projection
+projAir = reshape(W*double(maskAir(:)), W.proj_size);           %Air forward projection
+projIodine = reshape(W*double(maskIodine(:)), W.proj_size);               %Iodine forward projection
+projBone = reshape(W*double(maskBone(:)), W.proj_size);               %Iodine forward projection
+
+end