Skip to content
Snippets Groups Projects
Commit 47f8683f authored by Lisa Schätzle's avatar Lisa Schätzle
Browse files

Update generateOmega.m

parent 604ed21f
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 133 additions and 130 deletions
src/generateOmega.m
+ 133
− 130
View file @ 47f8683f
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

Help | Imprint | Privacy policy | Accessibility | Contact