diff --git a/src/generateOmega.m b/src/generateOmega.m
index 8e17366c96c6c32601ebb175bf281040655111d2..c0df53cb0863f62d0585ad07f8a4c60f0aa74b39 100644
--- a/src/generateOmega.m
+++ b/src/generateOmega.m
@@ -1,130 +1,133 @@
-function [P_Omega, P_OmegaC] = generateOmega(OmegaTypes, nxhat, nd, params)
-% INPUT: OmegaTypes - Cell-string-Array; possible strings are 'rectangle',
-% 'constant_height', 'constant_width', 'random', 'random position constant size'
-% nxhat - Number of observation directions/ rows
-% nd - Number of illumination directions/ columns
-% params - 5xlength(OmegaTypes) - array that specifies the
-% shapes of the components of Omega
- [~,nr_of_OmegaTypes] = size(OmegaTypes);
- P_OmegaC = ones(nxhat,nd);
- for iterk = 1:nr_of_OmegaTypes
- OmegaPos_xhat = params(iterk,1);
- OmegaPos_d = params(iterk,2);
- OmegaSize_xhat = params(iterk,3);
- OmegaSize_d = params(iterk,4);
- sigma = params(iterk,5);
- switch lower(OmegaTypes{1,iterk}) % Omega modells unknown (!) matrix entries
- case 'rectangle'
- omega = [OmegaPos_xhat, OmegaPos_xhat + OmegaSize_xhat];
- illuminationDirections = [OmegaPos_d, OmegaPos_d + OmegaSize_d];
- ind_omega_l = round(omega(1)*(nxhat-1)/(2*pi))+1;
- ind_omega_r = round(omega(2)*(nxhat-1)/(2*pi))+1;
- if omega(2) == 2*pi
- ind_omega_r = nxhat;
- end
- if omega(2) == pi
- ind_omega_r = nxhat/2;
- end
- height_omega = ind_omega_r-ind_omega_l+1;
- ind_illumiDir_l = round(illuminationDirections(1)*(nd-1)/(2*pi))+1;
- ind_illumiDir_r = round(illuminationDirections(2)*(nd-1)/(2*pi))+1;
- if illuminationDirections(2) == 2*pi
- ind_illumiDir_r = nd;
- end
- if illuminationDirections(2) == pi
- ind_illumiDir_r = nd/2;
- end
- width_illumiDir = ind_illumiDir_r-ind_illumiDir_l+1;
- clear omega illuminationDirections
- Omega = zeros(height_omega*width_illumiDir,2);
- for iterj = ind_illumiDir_l:ind_illumiDir_r
- help_iterj = iterj-ind_illumiDir_l+1;
- Omega(height_omega*(help_iterj-1)+1:height_omega*help_iterj,:) = [(ind_omega_l:ind_omega_r)', iterj*ones(height_omega,1)];
- end
- P_OmegaC_help = ones(nxhat,nd);
- P_OmegaC_help(sub2ind([nd,nxhat],Omega(:,1),Omega(:,2))) = 0;
- P_OmegaC = P_OmegaC .* P_OmegaC_help;
- clear ind_omega_l ind_omega_r heigth_omega ind_illumiDir_l ind_illumiDir_r width_illumiDir Omega help_iterj P_OmegaC_help
-
- case 'constant_height'
- omega = [OmegaPos_xhat, OmegaPos_xhat + OmegaSize_xhat];
- clear bla start
- ind_omega_l = round(omega(:,1)*(nxhat-1)/(2*pi)) + ones(nd,1);
- ind_omega_r = round(omega(:,2)*(nxhat-1)/(2*pi)) + ones(nd,1);
- P_OmegaC_help = ones(nxhat,nd);
- P_OmegaC_help(ind_omega_l:ind_omega_r,:) = 0;
- P_OmegaC = P_OmegaC .* P_OmegaC_help;
-
- clear ind_omega_l ind_omega_r P_OmegaC_help
-
- case 'constant_width'
- omega = [OmegaPos_d, OmegaPos_d + OmegaSize_d];
- clear bla start
- ind_omega_b = round(omega(:,1)*(nxhat-1)/(2*pi)) + ones(nxhat,1);
- ind_omega_t = round(omega(:,2)*(nxhat-1)/(2*pi)) + ones(nxhat,1);
- P_OmegaC_help = ones(nxhat,nd);
- P_OmegaC_help(:,ind_omega_b:ind_omega_t) = 0;
- P_OmegaC = P_OmegaC .* P_OmegaC_help;
-
- clear ind_omega_b ind_omega_t P_OmegaC_help
-
- case 'random'
- positions = 2*pi*rand(nd,1);
- sizes = OmegaSize_xhat*rand(nd,1);
- Omega = [positions, positions + sizes];
- ind_Omega_l = round(Omega(:,1)*(nxhat-1)/(2*pi)) + ones(nd,1);
- ind_Omega_r = round(Omega(:,2)*(nxhat-1)/(2*pi)) + ones(nd,1);
- Omega = zeros(1,2);
- for iterj = 1:nd
- left = ind_Omega_l(iterj,1);
- right = ind_Omega_r(iterj,1);
- length = right-left+1;
- if right > nxhat
- lengthLeft = length-(nxhat-left+1);
- NewOmega = [iterj*ones(length,1), vertcat((left:nxhat)',(1:lengthLeft)')];
- else
- NewOmega = [iterj*ones(length,1), (left:right)'];
- end
- Omega = vertcat(Omega, NewOmega);
- end
- Omega = Omega(2:end,:);
- P_OmegaC_help = ones(nxhat,nd);
- P_OmegaC_help(sub2ind([nd, nxhat],Omega(:,2),Omega(:,1))) = 0;
- P_OmegaC = P_OmegaC .* P_OmegaC_help;
- clear positions Omega ind_Omega_l ind_Omega_r left right NewOmega length lengthLeft P_OmegaC_help
-
- case 'random position constant size'
- positions = 2*pi*rand(nd,1);
- Omega = horzcat(positions, positions + OmegaSize_xhat * ones(nd,1));
- ind_Omega_l = round(Omega(:,1)*(nxhat-1)/(2*pi)) + ones(nd,1);
- ind_Omega_r = round(Omega(:,2)*(nxhat-1)/(2*pi)) + ones(nd,1);
- Omega = zeros(1,2);
- for iterj = 1:nd
- left = ind_Omega_l(iterj,1);
- right = ind_Omega_r(iterj,1);
- length = right-left+1;
- if right > nxhat
- lengthLeft = length-(nxhat-left+1);
- NewOmega = [iterj*ones(length,1), vertcat((left:nxhat)',(1:lengthLeft)')];
- else
- NewOmega = [iterj*ones(length,1), (left:right)'];
- end
- Omega = vertcat(Omega, NewOmega);
- end
- Omega = Omega(2:end,:);
- P_OmegaC_help = ones(nxhat,nd);
- P_OmegaC_help(sub2ind([nd, nxhat],Omega(:,2),Omega(:,1))) = 0;
- P_OmegaC = P_OmegaC .* P_OmegaC_help;
- clear positions Omega ind_Omega_l ind_Omega_r left right NewOmega length lengthLeft P_OmegaC_help
- case 'random entries'
- s = round(sigma*nxhat*nd); % l0 support of the additive sparse component A
- omega = randperm(nxhat*nd,s)';
- Omega = horzcat(omega-(ceil(omega/nxhat)-1)*nxhat, ceil(omega/nd)); % s*2 array
- P_OmegaC_help = ones(nxhat,nd);
- P_OmegaC_help(sub2ind([nd, nxhat],Omega(:,2),Omega(:,1))) = 0;
- P_OmegaC = P_OmegaC .* P_OmegaC_help;
- clear s Omega omega
- end
- end
- P_Omega = ones(nxhat,nd) - P_OmegaC;
-end
\ No newline at end of file
+function [P_Omega, P_OmegaC] = generateOmega(OmegaTypes, nxhat, nd, params)
+
+% Simulates projection operator corresponding to some missing data segment Omega.
+%
+% INPUT: OmegaTypes Cell-string-Array; possible strings are 'rectangle',
+% 'constant_height', 'constant_width', 'random', 'random position constant size'
+% nxhat - Number of observation directions/ rows
+% nd - Number of illumination directions/ columns
+% params - 5xlength(OmegaTypes) - array that specifies the
+% shapes of the components of Omega
+ [~,nr_of_OmegaTypes] = size(OmegaTypes);
+ P_OmegaC = ones(nxhat,nd);
+ for iterk = 1:nr_of_OmegaTypes
+ OmegaPos_xhat = params(iterk,1);
+ OmegaPos_d = params(iterk,2);
+ OmegaSize_xhat = params(iterk,3);
+ OmegaSize_d = params(iterk,4);
+ sigma = params(iterk,5);
+ switch lower(OmegaTypes{1,iterk}) % Omega modells unknown (!) matrix entries
+ case 'rectangle'
+ omega = [OmegaPos_xhat, OmegaPos_xhat + OmegaSize_xhat];
+ illuminationDirections = [OmegaPos_d, OmegaPos_d + OmegaSize_d];
+ ind_omega_l = round(omega(1)*(nxhat-1)/(2*pi))+1;
+ ind_omega_r = round(omega(2)*(nxhat-1)/(2*pi))+1;
+ if omega(2) == 2*pi
+ ind_omega_r = nxhat;
+ end
+ if omega(2) == pi
+ ind_omega_r = nxhat/2;
+ end
+ height_omega = ind_omega_r-ind_omega_l+1;
+ ind_illumiDir_l = round(illuminationDirections(1)*(nd-1)/(2*pi))+1;
+ ind_illumiDir_r = round(illuminationDirections(2)*(nd-1)/(2*pi))+1;
+ if illuminationDirections(2) == 2*pi
+ ind_illumiDir_r = nd;
+ end
+ if illuminationDirections(2) == pi
+ ind_illumiDir_r = nd/2;
+ end
+ width_illumiDir = ind_illumiDir_r-ind_illumiDir_l+1;
+ clear omega illuminationDirections
+ Omega = zeros(height_omega*width_illumiDir,2);
+ for iterj = ind_illumiDir_l:ind_illumiDir_r
+ help_iterj = iterj-ind_illumiDir_l+1;
+ Omega(height_omega*(help_iterj-1)+1:height_omega*help_iterj,:) = [(ind_omega_l:ind_omega_r)', iterj*ones(height_omega,1)];
+ end
+ P_OmegaC_help = ones(nxhat,nd);
+ P_OmegaC_help(sub2ind([nd,nxhat],Omega(:,1),Omega(:,2))) = 0;
+ P_OmegaC = P_OmegaC .* P_OmegaC_help;
+ clear ind_omega_l ind_omega_r heigth_omega ind_illumiDir_l ind_illumiDir_r width_illumiDir Omega help_iterj P_OmegaC_help
+
+ case 'constant_height'
+ omega = [OmegaPos_xhat, OmegaPos_xhat + OmegaSize_xhat];
+ clear bla start
+ ind_omega_l = round(omega(:,1)*(nxhat-1)/(2*pi)) + ones(nd,1);
+ ind_omega_r = round(omega(:,2)*(nxhat-1)/(2*pi)) + ones(nd,1);
+ P_OmegaC_help = ones(nxhat,nd);
+ P_OmegaC_help(ind_omega_l:ind_omega_r,:) = 0;
+ P_OmegaC = P_OmegaC .* P_OmegaC_help;
+
+ clear ind_omega_l ind_omega_r P_OmegaC_help
+
+ case 'constant_width'
+ omega = [OmegaPos_d, OmegaPos_d + OmegaSize_d];
+ clear bla start
+ ind_omega_b = round(omega(:,1)*(nxhat-1)/(2*pi)) + ones(nxhat,1);
+ ind_omega_t = round(omega(:,2)*(nxhat-1)/(2*pi)) + ones(nxhat,1);
+ P_OmegaC_help = ones(nxhat,nd);
+ P_OmegaC_help(:,ind_omega_b:ind_omega_t) = 0;
+ P_OmegaC = P_OmegaC .* P_OmegaC_help;
+
+ clear ind_omega_b ind_omega_t P_OmegaC_help
+
+ case 'random'
+ positions = 2*pi*rand(nd,1);
+ sizes = OmegaSize_xhat*rand(nd,1);
+ Omega = [positions, positions + sizes];
+ ind_Omega_l = round(Omega(:,1)*(nxhat-1)/(2*pi)) + ones(nd,1);
+ ind_Omega_r = round(Omega(:,2)*(nxhat-1)/(2*pi)) + ones(nd,1);
+ Omega = zeros(1,2);
+ for iterj = 1:nd
+ left = ind_Omega_l(iterj,1);
+ right = ind_Omega_r(iterj,1);
+ length = right-left+1;
+ if right > nxhat
+ lengthLeft = length-(nxhat-left+1);
+ NewOmega = [iterj*ones(length,1), vertcat((left:nxhat)',(1:lengthLeft)')];
+ else
+ NewOmega = [iterj*ones(length,1), (left:right)'];
+ end
+ Omega = vertcat(Omega, NewOmega);
+ end
+ Omega = Omega(2:end,:);
+ P_OmegaC_help = ones(nxhat,nd);
+ P_OmegaC_help(sub2ind([nd, nxhat],Omega(:,2),Omega(:,1))) = 0;
+ P_OmegaC = P_OmegaC .* P_OmegaC_help;
+ clear positions Omega ind_Omega_l ind_Omega_r left right NewOmega length lengthLeft P_OmegaC_help
+
+ case 'random position constant size'
+ positions = 2*pi*rand(nd,1);
+ Omega = horzcat(positions, positions + OmegaSize_xhat * ones(nd,1));
+ ind_Omega_l = round(Omega(:,1)*(nxhat-1)/(2*pi)) + ones(nd,1);
+ ind_Omega_r = round(Omega(:,2)*(nxhat-1)/(2*pi)) + ones(nd,1);
+ Omega = zeros(1,2);
+ for iterj = 1:nd
+ left = ind_Omega_l(iterj,1);
+ right = ind_Omega_r(iterj,1);
+ length = right-left+1;
+ if right > nxhat
+ lengthLeft = length-(nxhat-left+1);
+ NewOmega = [iterj*ones(length,1), vertcat((left:nxhat)',(1:lengthLeft)')];
+ else
+ NewOmega = [iterj*ones(length,1), (left:right)'];
+ end
+ Omega = vertcat(Omega, NewOmega);
+ end
+ Omega = Omega(2:end,:);
+ P_OmegaC_help = ones(nxhat,nd);
+ P_OmegaC_help(sub2ind([nd, nxhat],Omega(:,2),Omega(:,1))) = 0;
+ P_OmegaC = P_OmegaC .* P_OmegaC_help;
+ clear positions Omega ind_Omega_l ind_Omega_r left right NewOmega length lengthLeft P_OmegaC_help
+ case 'random entries'
+ s = round(sigma*nxhat*nd); % l0 support of the additive sparse component A
+ omega = randperm(nxhat*nd,s)';
+ Omega = horzcat(omega-(ceil(omega/nxhat)-1)*nxhat, ceil(omega/nd)); % s*2 array
+ P_OmegaC_help = ones(nxhat,nd);
+ P_OmegaC_help(sub2ind([nd, nxhat],Omega(:,2),Omega(:,1))) = 0;
+ P_OmegaC = P_OmegaC .* P_OmegaC_help;
+ clear s Omega omega
+ end
+ end
+ P_Omega = ones(nxhat,nd) - P_OmegaC;
+end