Initial revision

7 files changed, 331 insertions, 0 deletions

diff --git a/supp/RMSE.m b/supp/RMSE.m
new file mode 100644
index 0000000..734e4b8
--- /dev/null
+++ b/supp/RMSE.m
@@ -0,0 +1,7 @@
+function err = RMSE(signal1, signal2)
+%RMSE Root Mean Squared Error
+
+err = sum((signal1 - signal2).^2)/length(signal1);  % MSE
+err = sqrt(err);                                    % RMSE
+
+end
\ No newline at end of file
diff --git a/supp/add_wedges.m b/supp/add_wedges.m
new file mode 100644
index 0000000..8b8f2a7
--- /dev/null
+++ b/supp/add_wedges.m
@@ -0,0 +1,30 @@
+% create a wedge mask to simulate the missing wedge
+
+[rows, columns] = size(sino_zing_rings);
+grayImage = ones(rows, columns, 'uint8');
+xCoords = [0 360 0];
+yCoords = [35 7 7];
+mask = poly2mask(xCoords, yCoords, rows, columns);
+grayImage(mask) = 0; 
+
+xCoords = [727 360 727];
+yCoords = [35 7 7];
+mask = poly2mask(xCoords, yCoords, rows, columns);
+grayImage(mask) = 0; 
+
+xCoords = [0 360 0];
+yCoords = [145 173 173];
+mask = poly2mask(xCoords, yCoords, rows, columns);
+grayImage(mask) = 0; 
+
+xCoords = [727 360 727];
+yCoords = [145 173 173];
+mask = poly2mask(xCoords, yCoords, rows, columns);
+grayImage(mask) = 0; 
+
+grayImage(1:7,:) = 0;
+grayImage(173:end,:) = 0;
+
+%figure; imshow(grayImage, [0 1]);
+MW_sino_artifacts = sino_zing_rings.*double(grayImage);
+%Dweights = Dweights.*double(grayImage);
\ No newline at end of file
diff --git a/supp/filtersinc.m b/supp/filtersinc.m
new file mode 100644
index 0000000..6c29c98
--- /dev/null
+++ b/supp/filtersinc.m
@@ -0,0 +1,28 @@
+function g = filtersinc(PR)
+
+
+% filtersinc.m
+%
+% Written by Waqas Akram
+%
+% "a":	This parameter varies the filter magnitude response.
+%	When "a" is very small (a<<1), the response approximates |w|
+%	As "a" is increased, the filter response starts to 
+%	roll off at high frequencies.
+a = 1;
+
+[Length, Count] = size(PR);
+w = [-pi:(2*pi)/Length:pi-(2*pi)/Length];
+
+rn1 = abs(2/a*sin(a.*w./2));
+rn2 = sin(a.*w./2);
+rd = (a*w)./2;
+r = rn1*(rn2/rd)^2;
+
+f = fftshift(r);
+for i = 1:Count
+        IMG = fft(PR(:,i));
+        fimg = IMG.*f';
+        g(:,i) = ifft(fimg);
+end
+g = real(g);
\ No newline at end of file
diff --git a/supp/my_red_yellowMAP.mat b/supp/my_red_yellowMAP.mat
new file mode 100644
index 0000000..c2a5b87
--- /dev/null
+++ b/supp/my_red_yellowMAP.mat
diff --git a/supp/ssim_index.m b/supp/ssim_index.m
new file mode 100644
index 0000000..4fa7a79
--- /dev/null
+++ b/supp/ssim_index.m
@@ -0,0 +1,181 @@
+function [mssim, ssim_map] = ssim_index(img1, img2, K, window, L)
+
+%========================================================================
+%SSIM Index, Version 1.0
+%Copyright(c) 2003 Zhou Wang
+%All Rights Reserved.
+%
+%This is an implementation of the algorithm for calculating the
+%Structural SIMilarity (SSIM) index between two images. Please refer
+%to the following paper:
+%
+%Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, "Image
+%quality assessment: From error visibility to structural similarity"
+%IEEE Transactios on Image Processing, vol. 13, no. 4, pp.600-612,
+%Apr. 2004.
+%
+%Kindly report any suggestions or corrections to zhouwang@ieee.org
+%
+%----------------------------------------------------------------------
+%
+%Input : (1) img1: the first image being compared
+%        (2) img2: the second image being compared
+%        (3) K: constants in the SSIM index formula (see the above
+%            reference). defualt value: K = [0.01 0.03]
+%        (4) window: local window for statistics (see the above
+%            reference). default widnow is Gaussian given by
+%            window = fspecial('gaussian', 11, 1.5);
+%        (5) L: dynamic range of the images. default: L = 255
+%
+%Output: (1) mssim: the mean SSIM index value between 2 images.
+%            If one of the images being compared is regarded as 
+%            perfect quality, then mssim can be considered as the
+%            quality measure of the other image.
+%            If img1 = img2, then mssim = 1.
+%        (2) ssim_map: the SSIM index map of the test image. The map
+%            has a smaller size than the input images. The actual size:
+%            size(img1) - size(window) + 1.
+%
+%Default Usage:
+%   Given 2 test images img1 and img2, whose dynamic range is 0-255
+%
+%   [mssim ssim_map] = ssim_index(img1, img2);
+%
+%Advanced Usage:
+%   User defined parameters. For example
+%
+%   K = [0.05 0.05];
+%   window = ones(8);
+%   L = 100;
+%   [mssim ssim_map] = ssim_index(img1, img2, K, window, L);
+%
+%See the results:
+%
+%   mssim                        %Gives the mssim value
+%   imshow(max(0, ssim_map).^4)  %Shows the SSIM index map
+%
+%========================================================================
+
+
+if (nargin < 2 | nargin > 5)
+   ssim_index = -Inf;
+   ssim_map = -Inf;
+   return;
+end
+
+if (size(img1) ~= size(img2))
+   ssim_index = -Inf;
+   ssim_map = -Inf;
+   return;
+end
+
+[M N] = size(img1);
+
+if (nargin == 2)
+   if ((M < 11) | (N < 11))   % ͼССû塣
+           ssim_index = -Inf;
+           ssim_map = -Inf;
+      return
+   end
+   window = fspecial('gaussian', 11, 1.5);        % һ׼ƫ1.511*11ĸ˹ͨ˲
+   K(1) = 0.01;                                                                      % default settings
+   K(2) = 0.03;                                                                      %
+   L = 255;                                  %
+end
+
+if (nargin == 3)
+   if ((M < 11) | (N < 11))
+           ssim_index = -Inf;
+           ssim_map = -Inf;
+      return
+   end
+   window = fspecial('gaussian', 11, 1.5);
+   L = 255;
+   if (length(K) == 2)
+      if (K(1) < 0 | K(2) < 0)
+                   ssim_index = -Inf;
+                   ssim_map = -Inf;
+                   return;
+      end
+   else
+           ssim_index = -Inf;
+           ssim_map = -Inf;
+           return;
+   end
+end
+
+if (nargin == 4)
+   [H W] = size(window);
+   if ((H*W) < 4 | (H > M) | (W > N))
+           ssim_index = -Inf;
+           ssim_map = -Inf;
+      return
+   end
+   L = 255;
+   if (length(K) == 2)
+      if (K(1) < 0 | K(2) < 0)
+                   ssim_index = -Inf;
+                   ssim_map = -Inf;
+                   return;
+      end
+   else
+           ssim_index = -Inf;
+           ssim_map = -Inf;
+           return;
+   end
+end
+
+if (nargin == 5)
+   [H W] = size(window);
+   if ((H*W) < 4 | (H > M) | (W > N))
+           ssim_index = -Inf;
+           ssim_map = -Inf;
+      return
+   end
+   if (length(K) == 2)
+      if (K(1) < 0 | K(2) < 0)
+                   ssim_index = -Inf;
+                   ssim_map = -Inf;
+                   return;
+      end
+   else
+           ssim_index = -Inf;
+           ssim_map = -Inf;
+           return;
+   end
+end
+%%
+C1 = (K(1)*L)^2;    % C1Lxyá
+C2 = (K(2)*L)^2;    % C2ԱȶCxyá
+window = window/sum(sum(window)); %˲һ
+img1 = double(img1); 
+img2 = double(img2);
+
+mu1   = filter2(window, img1, 'valid');  % ͼ˲ӼȨ
+mu2   = filter2(window, img2, 'valid');  % ͼ˲ӼȨ
+
+mu1_sq = mu1.*mu1;     % Uxƽֵ
+mu2_sq = mu2.*mu2;     % Uyƽֵ
+mu1_mu2 = mu1.*mu2;    % Ux*Uyֵ
+
+sigma1_sq = filter2(window, img1.*img1, 'valid') - mu1_sq;  % sigmax ׼
+sigma2_sq = filter2(window, img2.*img2, 'valid') - mu2_sq;  % sigmay ׼
+sigma12 = filter2(window, img1.*img2, 'valid') - mu1_mu2;   % sigmaxy׼
+
+if (C1 > 0 & C2 > 0)
+   ssim_map = ((2*mu1_mu2 + C1).*(2*sigma12 + C2))./((mu1_sq + mu2_sq + C1).*(sigma1_sq + sigma2_sq + C2));
+else
+   numerator1 = 2*mu1_mu2 + C1;
+   numerator2 = 2*sigma12 + C2;
+   denominator1 = mu1_sq + mu2_sq + C1;
+   denominator2 = sigma1_sq + sigma2_sq + C2;
+   ssim_map = ones(size(mu1));
+   index = (denominator1.*denominator2 > 0);
+   ssim_map(index) = (numerator1(index).*numerator2(index))./(denominator1(index).*denominator2(index));
+   index = (denominator1 ~= 0) & (denominator2 == 0);
+   ssim_map(index) = numerator1(index)./denominator1(index);
+end
+
+mssim = mean2(ssim_map);
+
+return
\ No newline at end of file
diff --git a/supp/subplot_tight.m b/supp/subplot_tight.m
new file mode 100644
index 0000000..0b0cbd5
--- /dev/null
+++ b/supp/subplot_tight.m
@@ -0,0 +1 @@
+function vargout=subplot_tight(m, n, p, margins, varargin)
%% subplot_tight
% A subplot function substitude with margins user tunabble parameter.
%
%% Syntax
%  h=subplot_tight(m, n, p);
%  h=subplot_tight(m, n, p, margins);
%  h=subplot_tight(m, n, p, margins, subplotArgs...);
%
%% Description
% Our goal is to grant the user the ability to define the margins between neighbouring
%  subplots. Unfotrtunately Matlab subplot function lacks this functionality, and the
%  margins between subplots can reach 40% of figure area, which is pretty lavish. While at
%  the begining the function was implememnted as wrapper function for Matlab function
%  subplot, it was modified due to axes del;etion resulting from what Matlab subplot
%  detected as overlapping. Therefore, the current implmenetation makes no use of Matlab
%  subplot function, using axes instead. This can be problematic, as axis and subplot
%  parameters are quie different. Set isWrapper to "True" to return to wrapper mode, which
%  fully supports subplot format.
%
%% Input arguments (defaults exist):
%   margins- two elements vector [vertical,horizontal] defining the margins between
%        neighbouring axes. Default value is 0.04
%
%% Output arguments
%   same as subplot- none, or axes handle according to function call.
%
%% Issues & Comments
%  - Note that if additional elements are used in order to be passed to subplot, margins
%     parameter must be defined. For default margins value use empty element- [].
%  - 
%
%% Example
% close all;
% img=imread('peppers.png');
% figSubplotH=figure('Name', 'subplot');
% figSubplotTightH=figure('Name', 'subplot_tight');
% nElems=17;
% subplotRows=ceil(sqrt(nElems)-1);
% subplotRows=max(1, subplotRows);
% subplotCols=ceil(nElems/subplotRows);
% for iElem=1:nElems
%    figure(figSubplotH);
%    subplot(subplotRows, subplotCols, iElem);
%    imshow(img);
%    figure(figSubplotTightH);
%    subplot_tight(subplotRows, subplotCols, iElem, [0.0001]);
%    imshow(img);
% end
%
%% See also
%  - subplot
%
%% Revision history
% First version: Nikolay S. 2011-03-29.
% Last update:   Nikolay S. 2012-05-24.
%
% *List of Changes:*
% 2012-05-24
%  Non wrapping mode (based on axes command) added, to deal with an issue of disappearing
%     subplots occuring with massive axes.

%% Default params
isWrapper=false;
if (nargin<4) || isempty(margins)
    margins=[0.04,0.04]; % default margins value- 4% of figure
end
if length(margins)==1
    margins(2)=margins;
end

%note n and m are switched as Matlab indexing is column-wise, while subplot indexing is row-wise :(
[subplot_col,subplot_row]=ind2sub([n,m],p);  


height=(1-(m+1)*margins(1))/m; % single subplot height
width=(1-(n+1)*margins(2))/n;  % single subplot width

% note subplot suppors vector p inputs- so a merged subplot of higher dimentions will be created
subplot_cols=1+max(subplot_col)-min(subplot_col); % number of column elements in merged subplot 
subplot_rows=1+max(subplot_row)-min(subplot_row); % number of row elements in merged subplot   

merged_height=subplot_rows*( height+margins(1) )- margins(1);   % merged subplot height
merged_width= subplot_cols*( width +margins(2) )- margins(2);   % merged subplot width

merged_bottom=(m-max(subplot_row))*(height+margins(1)) +margins(1); % merged subplot bottom position
merged_left=min(subplot_col)*(width+margins(2))-width;              % merged subplot left position
pos=[merged_left, merged_bottom, merged_width, merged_height];


if isWrapper
   h=subplot(m, n, p, varargin{:}, 'Units', 'Normalized', 'Position', pos);
else
   h=axes('Position', pos, varargin{:});
end

if nargout==1
   vargout=h;
end
\ No newline at end of file
diff --git a/supp/zing_rings_add.m b/supp/zing_rings_add.m
new file mode 100644
index 0000000..023ac27
--- /dev/null
+++ b/supp/zing_rings_add.m
@@ -0,0 +1,84 @@
+%  uncomment this part of script to generate data with different noise characterisitcs
+
+fprintf('%s\n', 'Generating Projection Data...');
+multfactor = 1000;
+% Creating RHS (b) - the sinogram (using a strip projection model)
+vol_geom = astra_create_vol_geom(N, N);
+proj_geom = astra_create_proj_geom('parallel', 1.0, P, theta_rad);
+proj_id_temp = astra_create_projector('strip', proj_geom, vol_geom);
+[sinogram_id, sinogramIdeal] = astra_create_sino(phantom./multfactor, proj_id_temp);
+astra_mex_data2d('delete',sinogram_id);
+astra_mex_algorithm('delete',proj_id_temp);
+%
+% % adding Gaussian noise
+% eta = 0.04;  % Relative noise level
+% E = randn(size(sinogram));
+% sinogram = sinogram + eta*norm(sinogram,'fro')*E/norm(E,'fro');  % adding noise to the sinogram
+% sinogram(sinogram<0) = 0;
+% clear E;
+
+%%
+% adding zingers
+val_offset = 0;
+sino_zing = sinogramIdeal;
+vec1 = [60, 80, 80, 70, 70, 90, 90, 40, 130, 145, 155, 125];
+vec2 = [350, 450, 190, 500, 250, 530, 330, 230, 550, 250, 450, 195];
+for jj = 1:length(vec1)
+    for i1 = -2:2
+        for j1 = -2:2
+            sino_zing(vec1(jj)+i1, vec2(jj)+j1) = val_offset;
+        end
+    end
+end
+
+% adding stripes into the signogram
+sino_zing_rings = sino_zing;
+coeff = linspace2(0.01,0.15,180);
+vmax = max(sinogramIdeal(:));
+sino_zing_rings(1:180,120) = sino_zing_rings(1:180,120) + vmax*0.13;
+sino_zing_rings(80:180,209) = sino_zing_rings(80:180,209) + vmax*0.14;
+sino_zing_rings(50:110,210) = sino_zing_rings(50:110,210) + vmax*0.12;
+sino_zing_rings(1:180,211) = sino_zing_rings(1:180,211) + vmax*0.14;
+sino_zing_rings(1:180,300) = sino_zing_rings(1:180,300) + vmax*coeff(:);
+sino_zing_rings(1:180,301) = sino_zing_rings(1:180,301) + vmax*0.14;
+sino_zing_rings(10:100,302) = sino_zing_rings(10:100,302) + vmax*0.15;
+sino_zing_rings(90:180,350) = sino_zing_rings(90:180,350) + vmax*0.11;
+sino_zing_rings(60:140,410) = sino_zing_rings(60:140,410) + vmax*0.12;
+sino_zing_rings(1:180,411) = sino_zing_rings(1:180,411) + vmax*0.14;
+sino_zing_rings(1:180,412) = sino_zing_rings(1:180,412) + vmax*coeff(:);
+sino_zing_rings(1:180,413) = sino_zing_rings(1:180,413) + vmax*coeff(:);
+sino_zing_rings(1:180,500) = sino_zing_rings(1:180,500) - vmax*0.12;
+sino_zing_rings(1:180,501) = sino_zing_rings(1:180,501) - vmax*0.12;
+sino_zing_rings(1:180,550) = sino_zing_rings(1:180,550) + vmax*0.11;
+sino_zing_rings(1:180,551) = sino_zing_rings(1:180,551) + vmax*0.11;
+sino_zing_rings(1:180,552) = sino_zing_rings(1:180,552) + vmax*0.11;
+
+sino_zing_rings(sino_zing_rings < 0) = 0;
+%%
+
+% adding Poisson noise
+dose = 50000;
+dataExp = dose.*exp(-sino_zing_rings); % noiseless raw data
+dataPnoise = astra_add_noise_to_sino(dataExp,2*dose); % pre-log noisy raw data (weights)
+Dweights = dataPnoise; 
+sinogram = log(dose./dataPnoise);  %log corrected data -> sinogram
+sinogram = abs(sinogram);
+clear dataPnoise dataExp
+
+% normalizing
+sinogram = sinogram.*multfactor;
+sino_zing_rings = sinogram;
+Dweights = multfactor./Dweights;
+
+%
+% figure(1);
+% set(gcf, 'Position', get(0,'Screensize'));
+% subplot(1,2,1); imshow(phantom,[0 0.6]); title('Ideal Phantom'); colorbar;
+% subplot(1,2,2); imshow(sinogram,[0 180]);  title('Noisy Sinogram');  colorbar;
+% colormap(cmapnew);
+
+% figure;
+% set(gcf, 'Position', get(0,'Screensize'));
+% subplot(1,2,1); imshow(sinogramIdeal,[0 180]);  title('Ideal Sinogram');  colorbar;
+% imshow(sino_zing_rings,[0 180]); title('Noisy Sinogram with zingers and stripes');  colorbar;
+% colormap(cmapnew);
\ No newline at end of file