Description of GREIT_desired_img

0001 function PSF = GREIT_desired_img_sigmoid(xyz,radius, opt)
0002 %GREIT_DESIRED_IMG_SIGMOID sigmoid-decay desired image function for GREIT
0003 % PSF= GREIT_desired_img_sigmoid(xyc, radius, opt)
0004 %   xyz     - array of centers of desired target images
0005 %               [2xN] - xy only
0006 %               [3xN] - xyz or, if radius = [], xyr
0007 %               [4xN] - xyzr (radius is ignored)
0008 %   radius  - the radius of the target on the desired image as a fraction
0009 %             of the model radius (half the larger dimension in xy)
0010 %   opt     - a struct with these fields:
0011 %      .rec_model   a 2D or 3D model as generated by MK_GRID_MODEL that may
0012 %                   have had some pixels/voxels removed (a rec_model that
0013 %                   is not a complete rectangle/parallelepiped).
0014 %      .steepness   [optional] a factor controling the amount of blur, see
0015 %                   below for details. Lower steepness gives more blur, but
0016 %                   if the value is too low, image may not reach the value
0017 %                   of 1 at the center of the target.
0018 %                   May be specified as a scalar, a [1xN] vector, or a
0019 %                   function handle with the signature:
0020 %                        func(pts)
0021 %                   where pts is either [2xN] or [3xN], dependig on the xyz
0022 %                   function input, e.g. @(xyz) 50./(xyz(3,:))
0023 %                   Default: 10./radius
0024 %      .desired_img_radius
0025 %                   [optional] Overwrites the radius input. May be
0026 %                   specified as a scalar, a [1xN] vector or a function
0027 %                   handle with the signature:
0028 %                        func(pts)
0029 %                   where pts is either [2xN] or [3xN], dependig on the xyz
0030 %                   function input, e.g. @(xyz) abs(xyz(3,:))/5
0031 %      .threshold   [optional] voxels where the function is smaller than
0032 %                   threshold will be set 0. Values bigger than 1-threshold
0033 %                   will be 1. The smaller the threshold the more
0034 %                   computationally expensive is the evaluation.
0035 %                   Default: 1e-4;
0036 %
0037 % The desired images approximate in each pixel the area integral of:
0038 %       f(r) = 1 / (1 + exp(s*(|r-r0| - radius)))
0039 % where
0040 %       r   - position vector in 2D/3D space
0041 %       s   - opt.steepness
0042 %       r0  - target center
0043 % For |r-r0| = radius, f(r) = 0.5.
0044 %
0045 % As of 2015-03-29, this is the default desired image function used by
0046 % MK_GREIT_MODEL.
0047 %
0048 % See also: CALC_GREIT_RM, MK_GREIT_MODEL, MK_PIXEL_SLICE
0049 
0050 % (C) 2015 Bartlomiej Grychtol - all rights reserved Swisstom AG.
0051 % License: GPL version 2 or 3
0052 % $Id: GREIT_desired_img_sigmoid.m 6044 2019-12-31 20:07:48Z aadler $
0053 
0054 % >> SWISSTOM CONTRIBUTION <<
0055 
0056 if ischar(xyz) && strcmp(xyz,'UNIT_TEST'), do_unit_test; return, end
0057 
0058 [xyzr, radius, opt] = parse_opt(xyz, radius, opt);
0059 
0060 copt.cache_obj = {xyzr, radius, opt.rec_model.nodes, opt.rec_model.elems, opt.steepness};
0061 copt.fstr = 'GREIT_desired_img_sigmoid';
0062 PSF = eidors_cache(@desired_soln,{xyzr, radius, opt},copt);
0063 
0064 
0065 end
0066 
0067 %-------------------------------------------------------------------------%
0068 % The main function
0069 function PSF = desired_soln(xyz, radius, opt)
0070    num_it = size(xyz,2);
0071     progress_msg('Composing desired images:',0,num_it);
0072     mdl = opt.rec_model;
0073     switch opt.n_dim
0074         case 3
0075             mdl.vox = [mdl.elems(1:6:end,:) mdl.elems(6:6:end,:)];
0076 
0077             min_vox_edge = min( [min(diff(unique(mdl.nodes(:,1)))), ...
0078                                  min(diff(unique(mdl.nodes(:,2)))), ...
0079                                  min(diff(unique(mdl.nodes(:,3))))] );
0080         case 2
0081             mdl.vox = [mdl.elems(1:2:end,:) mdl.elems(2:2:end,:)];
0082             min_vox_edge = min( [min(diff(unique(mdl.nodes(:,1)))), ...
0083                                  min(diff(unique(mdl.nodes(:,2))))] );
0084 
0085     end
0086     [Xnodes,Ynodes,Znodes] = voxnodes(mdl);
0087     
0088 
0089     
0090     warned = false;
0091     interp_elem_new('reset',size(mdl.vox,1),[],opt);
0092     
0093     % estimate the amount of memory needed to store PSF
0094     n_el = size(mdl.vox,1);
0095     n_pt = size(xyz,2);
0096     n_ep = n_el*n_pt;
0097     ll = eidors_msg('log_level',0);
0098     vox_box = sum(get_elem_volume(mdl));
0099     eidors_msg('log_level',ll);
0100     tgt_box = (2*mean(radius+opt.threshold./opt.steepness))^3;
0101     n_nz = ceil(n_ep*tgt_box/vox_box);
0102     SPARSE_STEP = ceil(n_nz/5);
0103     PSF = spalloc(n_el,n_pt,n_nz);
0104     
0105     for i=1:size(xyz,2)
0106         if nnz(PSF) > n_nz
0107            n_nz = nnz(PSF) + SPARSE_STEP;
0108            [r,c,v] = find(PSF);
0109            PSF = sparse(r,c,v,n_el,n_pt,n_nz); 
0110         end
0111         th = opt.threshold/opt.steepness(i);
0112         progress_msg(i,num_it);
0113         farel = far_elems(Xnodes,Ynodes,Znodes, xyz(:,i), radius(i), th);
0114         % also check for elements so close to the center they are approx 1
0115         close_el = false(size(farel));
0116         if radius(i) - th > min_vox_edge
0117            close_el = close_elems(Xnodes,Ynodes,Znodes, xyz(:,i), radius(i), th);
0118         end
0119         el_idx = find(~farel & ~close_el);
0120         % start with an initial size
0121         STEP = 1000;
0122         idx = zeros(STEP,1);
0123         factor = zeros(STEP,1);
0124         X = zeros(STEP,1);
0125         Y = zeros(STEP,1);
0126         Z = zeros(STEP,1);
0127         L = STEP;
0128         N = 1;
0129 
0130         for e = el_idx'
0131             [x,y,z] = interp_elem_new(mdl,e,radius(i), opt);
0132             n = numel(x);
0133             N1 = N; 
0134             N = N+n;
0135             if N > L % expand arrays as needed
0136                 fill = zeros(STEP,1);
0137                 idx = [idx; fill];
0138                 factor = [factor; fill];
0139                 X = [X; fill]; Y = [Y; fill]; Z = [Z; fill];
0140                 L = L + STEP;
0141             end
0142             v = N1:N-1;
0143             idx(v) = e;
0144             factor(v) = 1/n;
0145             X(v) = x; Y(v) = y; Z(v) = z;
0146         end
0147         % remove spare length
0148         v = N:L;
0149         idx(v) = [];
0150         factor(v) = [];
0151         X(v) = []; Y(v) = []; Z(v) = [];
0152         
0153         if isempty(X)
0154             if ~warned
0155                warning('EIDORS:OutsidePoint',...
0156                   'Desired image generation failed for point %d (and maybe others)',i);
0157                warned = true;
0158             end
0159             % PSF(:,i) = 0; % not needed
0160             continue;
0161         end
0162         V = [X Y Z];
0163         D = sqrt(sum(bsxfun(@minus,V(:,1:opt.n_dim),xyz(1:opt.n_dim,i)').^2, 2));
0164         x = D - radius(i);
0165         tmp = 1 ./ (1 + exp( opt.steepness(i) * x));
0166         PSF(:,i) = sparse(idx,1,tmp(:) .* factor,size(mdl.vox,1),1);
0167         PSF(close_el,i) = 1;
0168     end
0169     interp_elem_new('clear');
0170     progress_msg(Inf);
0171 end
0172 
0173 %-------------------------------------------------------------------------%
0174 % Generate internal points in elements
0175 function [X, Y, Z] = voxnodes(mdl)
0176     x = mdl.nodes(:,1); X = x(mdl.vox);
0177     y = mdl.nodes(:,2); Y = y(mdl.vox);
0178     try
0179         z = mdl.nodes(:,3); Z = z(mdl.vox);
0180     catch
0181         Z = [];
0182     end  
0183 end
0184 
0185 function [x,y,z] = interp_elem_new(mdl,e,radius,opt)
0186     persistent N_entries X  Y Z MAP N minnode maxnode done_elems sep
0187     if ischar(mdl) && strcmp(mdl,'reset')
0188         N_entries = 0; X = []; Y = []; Z = []; MAP = [];
0189         N       = ones(1,3);
0190         minnode = zeros(e,opt.n_dim);
0191         maxnode = zeros(e,opt.n_dim);
0192         sep     = zeros(e,opt.n_dim);
0193         done_elems = false(e,1);
0194         return;
0195     end
0196     if ischar(mdl) && strcmp(mdl,'clear')
0197        clear N_entries X  Y Z MAP N minnode maxnode done_elems sep
0198        return;
0199     end
0200     maxsep = radius/5;
0201     
0202     if ~done_elems(e)
0203        minnode(e,:) = min(mdl.nodes(mdl.vox(e,:),:));
0204        maxnode(e,:) = max(mdl.nodes(mdl.vox(e,:),:));
0205        sep(e,:) = maxnode(e,:) - minnode(e,:);
0206        done_elems(e) = true;
0207     end
0208     
0209     N(1:opt.n_dim) = max(3, ceil(sep(e,:)/maxsep)+1);
0210     
0211     try
0212         entry = MAP(N(1),N(2),N(3));
0213         x = X{entry};
0214         y = Y{entry};
0215         z = Z{entry};
0216     catch
0217         entry = N_entries+1;
0218         
0219         vx = linvec(0,1,N(1));
0220         vx = vx + .5*(vx(2) -vx(1));
0221         
0222         vy = linvec(0,1,N(2));
0223         vy = vy + .5*(vy(2) -vy(1));
0224         
0225         
0226         switch opt.n_dim
0227             case 2
0228                 [x, y] = grid2d(vx,vy);
0229                 z = zeros(size(x));
0230             case 3
0231                 vz = linvec(0,1,N(3));
0232                 vz = vz + .5*(vz(2) -vz(1));
0233                 [x, y, z] = grid3d(vx,vy,vz);
0234         end
0235         X{entry} = x;
0236         Y{entry} = y;
0237         Z{entry} = z;
0238         N_entries = entry;
0239         MAP(N(1),N(2),N(3)) = entry;
0240     end
0241     x = x*sep(e,1) + minnode(e,1);
0242     y = y*sep(e,2) + minnode(e,2);
0243     if opt.n_dim == 3
0244        z = z*sep(e,3) + minnode(e,3);
0245     end
0246     
0247 end
0248 
0249 %-------------------------------------------------------------------------%
0250 % Generate internal points in elements
0251 function [x,y,z] = interp_elem(mdl,e,radius, opt)
0252     maxsep = radius/5;
0253 
0254     minnode = min(mdl.nodes(mdl.vox(e,:),:));
0255     maxnode = max(mdl.nodes(mdl.vox(e,:),:));
0256     
0257     sep = maxnode - minnode;
0258     N = max(3, ceil(sep/maxsep)+1);
0259 
0260     vx = linvec(minnode(1),maxnode(1),N(1));
0261     vx = vx + .5*(vx(2) -vx(1));
0262         
0263     vy = linvec(minnode(2),maxnode(2),N(2));
0264     vy = vy + .5*(vy(2) -vy(1));
0265     
0266 
0267     switch opt.n_dim
0268         case 2
0269             [x, y] = grid2d(vx,vy);
0270             z = [];
0271         case 3
0272             vz = linvec(minnode(3),maxnode(3),N(3));
0273             vz = vz + .5*(vz(2) -vz(1));
0274             [x, y, z] = grid3d(vx,vy,vz);
0275     end
0276 end
0277 
0278 
0279 function out = linvec(v1, v2, N)
0280 % equivalent to out= linspace(v1,v2,N); out(end) = []; minus error checking
0281 out = v1 + (0:(N-2)).*(v2-v1)/(N-1);
0282 
0283 end
0284 
0285 function [x, y] = grid2d(vx, vy)
0286 % ndgrid minus overhead, asumes both vx and vy are row vectors
0287     vx = vx.';
0288     x = repmat(vx,size(vy));
0289     y = repmat(vy,size(vx));
0290 end
0291 
0292 function [x, y, z] = grid3d(vx,vy,vz)
0293 % ndgrid minus overhead
0294     sz = [numel(vx) numel(vy) numel(vz)];
0295     x = reshape(vx,[sz(1) 1 1]);
0296     x = repmat(x,[1 sz(2) sz(3)]);
0297     
0298     y = reshape(vy,[1 sz(2) 1]);
0299     y = repmat(y,[sz(1) 1 sz(3)]);
0300     
0301     z = reshape(vz,[1 1 sz(3)]);
0302     z = repmat(z,[sz(1) sz(2) 1]);
0303     
0304 
0305 end
0306 
0307 
0308 %-------------------------------------------------------------------------%
0309 % Find elements where the function value is negligable
0310 function farel = far_elems(Xnodes,Ynodes,Znodes,xyz,radius, th)
0311     farel = false(size(Xnodes,1),1);
0312 
0313     nodes_test = Xnodes < xyz(1) - radius - th;
0314     farel = farel | all(nodes_test,2);
0315     if all(farel), return, end;
0316     nodes_test = Xnodes > xyz(1) + radius + th;
0317     farel = farel | all(nodes_test,2);
0318     if all(farel), return, end;
0319     nodes_test = Ynodes < xyz(2) - radius - th;
0320     farel = farel | all(nodes_test,2);
0321     if all(farel), return, end;
0322     nodes_test = Ynodes > xyz(2) + radius + th;
0323     farel = farel | all(nodes_test,2);
0324     if all(farel), return, end;
0325     if ~isempty(Znodes)
0326         nodes_test  = Znodes > xyz(3) + radius + th;
0327         farel = farel | all(nodes_test,2);
0328         if all(farel), return, end;
0329         nodes_test  = Znodes < xyz(3) - radius - th;
0330         farel = farel | all(nodes_test,2);
0331         if all(farel), return, end;
0332     end
0333     idx = find(~farel);
0334 end
0335 
0336 %-------------------------------------------------------------------------%
0337 % Find elements where the function value is essentially 1
0338 function farel = close_elems(Xnodes,Ynodes,Znodes,xyz,radius, th)
0339   
0340    farel = true(size(Xnodes,1),1);
0341 
0342     nodes_test = Xnodes < xyz(1) + radius - th;
0343     farel = farel & all(nodes_test,2);
0344     if ~any(farel), return, end;
0345     nodes_test = Xnodes > xyz(1) - radius + th;
0346     farel = farel & all(nodes_test,2);
0347     if ~any(farel), return, end;
0348     nodes_test = Ynodes < xyz(2) + radius - th;
0349     farel = farel & all(nodes_test,2);
0350     if ~any(farel), return, end;
0351     nodes_test = Ynodes > xyz(2) - radius + th;
0352     farel = farel & all(nodes_test,2);
0353     if ~any(farel), return, end;
0354     if ~isempty(Znodes)
0355         nodes_test  = Znodes > xyz(3) - radius + th;
0356         farel = farel & all(nodes_test,2);
0357         if ~any(farel), return, end;
0358         nodes_test  = Znodes < xyz(3) + radius - th;
0359         farel = farel & all(nodes_test,2);
0360         if ~any(farel), return, end;
0361     end
0362     idx = find(farel);
0363 end
0364 
0365 
0366 %-------------------------------------------------------------------------%
0367 % Parse options
0368 function [xyz, radius, opt] = parse_opt(xyz, radius, opt)
0369 
0370     scale_radius = false;
0371     if isempty(radius)
0372         radius = xyz(end,:);
0373         scale_radius = true;
0374         xyz(end,:) = [];
0375     end
0376     
0377     if isfield(opt,'desired_img_radius')
0378        scale_radius = false;
0379        if isnumeric(opt.desired_img_radius)
0380           radius = opt.desired_img_radius;
0381        end
0382        if isa(opt.desired_img_radius, 'function_handle')
0383           radius = feval(opt.desired_img_radius,xyz);
0384        end
0385     end
0386         
0387  
0388     
0389     mdl = opt.rec_model; % must exist
0390     opt.n_dim =mdl_dim(mdl);
0391     xyz = xyz(1:opt.n_dim, :); % ingore z if model is 2D
0392 
0393     
0394 
0395     opt.meshsz = [];
0396     try
0397         for i = 1:3
0398             opt.meshsz = [opt.meshsz min(mdl.nodes(:,i)) max(mdl.nodes(:,i))];
0399         end
0400     end
0401     
0402     opt.n_dim = length(opt.meshsz)/2;
0403     opt.meshsz = reshape(opt.meshsz,2,[])';
0404 
0405     opt.minpt = opt.meshsz(:,1);
0406     opt.maxpt = opt.meshsz(:,2);
0407     opt.range = opt.maxpt - opt.minpt;
0408     opt.maxrange = max(opt.range(1:2));
0409 
0410     %rescale points to between -1 and 1 in the x-y plane
0411     xyz = 2*bsxfun(@minus, xyz, mean(opt.meshsz,2))/opt.maxrange;
0412     mdl.nodes = 2*bsxfun(@minus, mdl.nodes, mean(opt.meshsz(1:size(mdl.nodes,2),:),2)')/opt.maxrange;
0413     opt.rec_model = mdl;
0414     if scale_radius
0415         radius = 2*radius / opt.maxrange;
0416     end
0417 
0418     if ~isfield(opt, 'steepness')
0419         opt.steepness = 10./radius;
0420     elseif isa(opt.steepness,'function_handle')
0421         opt.steepness = feval(opt.steepness,xyz);
0422     end
0423     
0424     if numel(opt.steepness) == 1
0425        opt.steepness = ones(1,size(xyz,2)) * opt.steepness;
0426     end
0427     
0428     if numel(radius) == 1
0429         radius = ones(1,size(xyz,2)) * radius;
0430     end
0431     
0432     if opt.n_dim == 2
0433        xyz(3,:) = 0;
0434     end
0435    
0436     if ~isfield(opt, 'threshold')
0437        opt.threshold = 1e-4;
0438     end
0439     
0440     opt.threshold = log(1/opt.threshold - 1); 
0441     
0442 end
0443 
0444 %-------------------------------------------------------------------------%
0445 % UNIT_TEST
0446 function do_unit_test
0447     eidors_cache clear_name GREIT_desired_img_sigmoid
0448     v = linspace(-1,1,32);
0449     mdl= mk_grid_model([],v,v,[0 .7 1:.2:2 2.3 3]);
0450     opt.rec_model = mdl;
0451     opt.steepness = @(xyz) 50./xyz(3,:);
0452     opt.desired_img_radius = @(xyz) xyz(3,:)/5;
0453 %     xyzr = zeros(5,4);
0454     xyzr(:,3) = .5:.5:2.5;
0455     xyzr(:,4) = .25;
0456     sol = GREIT_desired_img_sigmoid(xyzr',[],opt);
0457     img = mk_image(mdl,0);
0458     for i = 1:5
0459         subplot(2,3,i)
0460         img.elem_data= sol(:,i);
0461         show_3d_slices(img,xyzr(i,3),xyzr(i,2),xyzr(i,1));
0462     end
0463     eidors_msg('UNIT_TEST: Showed %d images to verify',i,0);
0464     eidors_cache on GREIT_desired_img_sigmoid
0465 end
GREIT_desired_img_sigmoid

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
