diff options
Diffstat (limited to 'supp')
| -rw-r--r-- | supp/RMSE.m | 7 | ||||
| -rw-r--r-- | supp/add_wedges.m | 30 | ||||
| -rw-r--r-- | supp/filtersinc.m | 28 | ||||
| -rw-r--r-- | supp/my_red_yellowMAP.mat | bin | 0 -> 1761 bytes | |||
| -rw-r--r-- | supp/ssim_index.m | 181 | ||||
| -rw-r--r-- | supp/subplot_tight.m | 1 | ||||
| -rw-r--r-- | supp/zing_rings_add.m | 84 |
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 Binary files differnew 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 |