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 6367 2022-05-05 14:03:03Z 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     Npts = size(xyz,2);
0106     Nfailed = 0;
0107     for i=1:Npts
0108         if nnz(PSF) > n_nz
0109            n_nz = nnz(PSF) + SPARSE_STEP;
0110            [r,c,v] = find(PSF);
0111            PSF = sparse(r,c,v,n_el,n_pt,n_nz); 
0112         end
0113         th = opt.threshold/opt.steepness(i);
0114         progress_msg(i,num_it);
0115         farel = far_elems(Xnodes,Ynodes,Znodes, xyz(:,i), radius(i), th);
0116         % also check for elements so close to the center they are approx 1
0117         close_el = false(size(farel));
0118         if radius(i) - th > min_vox_edge
0119            close_el = close_elems(Xnodes,Ynodes,Znodes, xyz(:,i), radius(i), th);
0120         end
0121         el_idx = find(~farel & ~close_el);
0122         % start with an initial size
0123         STEP = 1000;
0124         idx = zeros(STEP,1);
0125         factor = zeros(STEP,1);
0126         X = zeros(STEP,1);
0127         Y = zeros(STEP,1);
0128         Z = zeros(STEP,1);
0129         L = STEP;
0130         N = 1;
0131 
0132         for e = el_idx'
0133             [x,y,z] = interp_elem_new(mdl,e,radius(i), opt);
0134             n = numel(x);
0135             N1 = N; 
0136             N = N+n;
0137             if N > L % expand arrays as needed
0138                 fill = zeros(STEP,1);
0139                 idx = [idx; fill];
0140                 factor = [factor; fill];
0141                 X = [X; fill]; Y = [Y; fill]; Z = [Z; fill];
0142                 L = L + STEP;
0143             end
0144             v = N1:N-1;
0145             idx(v) = e;
0146             factor(v) = 1/n;
0147             X(v) = x; Y(v) = y; Z(v) = z;
0148         end
0149         % remove spare length
0150         v = N:L;
0151         idx(v) = [];
0152         factor(v) = [];
0153         X(v) = []; Y(v) = []; Z(v) = [];
0154         
0155         if isempty(X)
0156 %  REMOVED: Now only one warning at the end
0157 %           if 0 && ~warned
0158 %              warning('EIDORS:OutsidePoint',...
0159 %                 'Desired image generation failed for point %d (and maybe others)',i);
0160 %              warned = true;
0161 %           end
0162             % PSF(:,i) = 0; % not needed
0163             Nfailed = Nfailed + 1;
0164             continue;
0165         end
0166         V = [X Y Z];
0167         D = sqrt(sum(bsxfun(@minus,V(:,1:opt.n_dim),xyz(1:opt.n_dim,i)').^2, 2));
0168         x = D - radius(i);
0169         tmp = 1 ./ (1 + exp( opt.steepness(i) * x));
0170         PSF(:,i) = sparse(idx,1,tmp(:) .* factor,size(mdl.vox,1),1);
0171         PSF(close_el,i) = 1;
0172     end
0173     interp_elem_new('clear');
0174     if Nfailed==0
0175         progress_msg(Inf);
0176     else
0177         progress_msg(sprintf( ...
0178       'Interpolation failed for %d of %d points', ...
0179                 Nfailed,Npts), Inf);
0180     end
0181   
0182 end
0183 
0184 %-------------------------------------------------------------------------%
0185 % Generate internal points in elements
0186 function [X, Y, Z] = voxnodes(mdl)
0187     x = mdl.nodes(:,1); X = x(mdl.vox);
0188     y = mdl.nodes(:,2); Y = y(mdl.vox);
0189     try
0190         z = mdl.nodes(:,3); Z = z(mdl.vox);
0191     catch
0192         Z = [];
0193     end  
0194 end
0195 
0196 function [x,y,z] = interp_elem_new(mdl,e,radius,opt)
0197     persistent N_entries X  Y Z MAP N minnode maxnode done_elems sep
0198     if ischar(mdl) && strcmp(mdl,'reset')
0199         N_entries = 0; X = []; Y = []; Z = []; MAP = [];
0200         N       = ones(1,3);
0201         minnode = zeros(e,opt.n_dim);
0202         maxnode = zeros(e,opt.n_dim);
0203         sep     = zeros(e,opt.n_dim);
0204         done_elems = false(e,1);
0205         return;
0206     end
0207     if ischar(mdl) && strcmp(mdl,'clear')
0208        clear N_entries X  Y Z MAP N minnode maxnode done_elems sep
0209        return;
0210     end
0211     maxsep = radius/5;
0212     
0213     if ~done_elems(e)
0214        minnode(e,:) = min(mdl.nodes(mdl.vox(e,:),:));
0215        maxnode(e,:) = max(mdl.nodes(mdl.vox(e,:),:));
0216        sep(e,:) = maxnode(e,:) - minnode(e,:);
0217        done_elems(e) = true;
0218     end
0219     
0220     N(1:opt.n_dim) = max(3, ceil(sep(e,:)/maxsep)+1);
0221     
0222     try
0223         entry = MAP(N(1),N(2),N(3));
0224         x = X{entry};
0225         y = Y{entry};
0226         z = Z{entry};
0227     catch
0228         entry = N_entries+1;
0229         
0230         vx = linvec(0,1,N(1));
0231         vx = vx + .5*(vx(2) -vx(1));
0232         
0233         vy = linvec(0,1,N(2));
0234         vy = vy + .5*(vy(2) -vy(1));
0235         
0236         
0237         switch opt.n_dim
0238             case 2
0239                 [x, y] = grid2d(vx,vy);
0240                 z = zeros(size(x));
0241             case 3
0242                 vz = linvec(0,1,N(3));
0243                 vz = vz + .5*(vz(2) -vz(1));
0244                 [x, y, z] = grid3d(vx,vy,vz);
0245         end
0246         X{entry} = x;
0247         Y{entry} = y;
0248         Z{entry} = z;
0249         N_entries = entry;
0250         MAP(N(1),N(2),N(3)) = entry;
0251     end
0252     x = x*sep(e,1) + minnode(e,1);
0253     y = y*sep(e,2) + minnode(e,2);
0254     if opt.n_dim == 3
0255        z = z*sep(e,3) + minnode(e,3);
0256     end
0257     
0258 end
0259 
0260 %-------------------------------------------------------------------------%
0261 % Generate internal points in elements
0262 function [x,y,z] = interp_elem(mdl,e,radius, opt)
0263     maxsep = radius/5;
0264 
0265     minnode = min(mdl.nodes(mdl.vox(e,:),:));
0266     maxnode = max(mdl.nodes(mdl.vox(e,:),:));
0267     
0268     sep = maxnode - minnode;
0269     N = max(3, ceil(sep/maxsep)+1);
0270 
0271     vx = linvec(minnode(1),maxnode(1),N(1));
0272     vx = vx + .5*(vx(2) -vx(1));
0273         
0274     vy = linvec(minnode(2),maxnode(2),N(2));
0275     vy = vy + .5*(vy(2) -vy(1));
0276     
0277 
0278     switch opt.n_dim
0279         case 2
0280             [x, y] = grid2d(vx,vy);
0281             z = [];
0282         case 3
0283             vz = linvec(minnode(3),maxnode(3),N(3));
0284             vz = vz + .5*(vz(2) -vz(1));
0285             [x, y, z] = grid3d(vx,vy,vz);
0286     end
0287 end
0288 
0289 
0290 function out = linvec(v1, v2, N)
0291 % equivalent to out= linspace(v1,v2,N); out(end) = []; minus error checking
0292 out = v1 + (0:(N-2)).*(v2-v1)/(N-1);
0293 
0294 end
0295 
0296 function [x, y] = grid2d(vx, vy)
0297 % ndgrid minus overhead, asumes both vx and vy are row vectors
0298     vx = vx.';
0299     x = repmat(vx,size(vy));
0300     y = repmat(vy,size(vx));
0301 end
0302 
0303 function [x, y, z] = grid3d(vx,vy,vz)
0304 % ndgrid minus overhead
0305     sz = [numel(vx) numel(vy) numel(vz)];
0306     x = reshape(vx,[sz(1) 1 1]);
0307     x = repmat(x,[1 sz(2) sz(3)]);
0308     
0309     y = reshape(vy,[1 sz(2) 1]);
0310     y = repmat(y,[sz(1) 1 sz(3)]);
0311     
0312     z = reshape(vz,[1 1 sz(3)]);
0313     z = repmat(z,[sz(1) sz(2) 1]);
0314     
0315 
0316 end
0317 
0318 
0319 %-------------------------------------------------------------------------%
0320 % Find elements where the function value is negligable
0321 function farel = far_elems(Xnodes,Ynodes,Znodes,xyz,radius, th)
0322     farel = false(size(Xnodes,1),1);
0323 
0324     nodes_test = Xnodes < xyz(1) - radius - th;
0325     farel = farel | all(nodes_test,2);
0326     if all(farel), return, end;
0327     nodes_test = Xnodes > xyz(1) + radius + th;
0328     farel = farel | all(nodes_test,2);
0329     if all(farel), return, end;
0330     nodes_test = Ynodes < xyz(2) - radius - th;
0331     farel = farel | all(nodes_test,2);
0332     if all(farel), return, end;
0333     nodes_test = Ynodes > xyz(2) + radius + th;
0334     farel = farel | all(nodes_test,2);
0335     if all(farel), return, end;
0336     if ~isempty(Znodes)
0337         nodes_test  = Znodes > xyz(3) + radius + th;
0338         farel = farel | all(nodes_test,2);
0339         if all(farel), return, end;
0340         nodes_test  = Znodes < xyz(3) - radius - th;
0341         farel = farel | all(nodes_test,2);
0342         if all(farel), return, end;
0343     end
0344     idx = find(~farel);
0345 end
0346 
0347 %-------------------------------------------------------------------------%
0348 % Find elements where the function value is essentially 1
0349 function farel = close_elems(Xnodes,Ynodes,Znodes,xyz,radius, th)
0350   
0351    farel = true(size(Xnodes,1),1);
0352 
0353     nodes_test = Xnodes < xyz(1) + radius - th;
0354     farel = farel & all(nodes_test,2);
0355     if ~any(farel), return, end;
0356     nodes_test = Xnodes > xyz(1) - radius + th;
0357     farel = farel & all(nodes_test,2);
0358     if ~any(farel), return, end;
0359     nodes_test = Ynodes < xyz(2) + radius - th;
0360     farel = farel & all(nodes_test,2);
0361     if ~any(farel), return, end;
0362     nodes_test = Ynodes > xyz(2) - radius + th;
0363     farel = farel & all(nodes_test,2);
0364     if ~any(farel), return, end;
0365     if ~isempty(Znodes)
0366         nodes_test  = Znodes > xyz(3) - radius + th;
0367         farel = farel & all(nodes_test,2);
0368         if ~any(farel), return, end;
0369         nodes_test  = Znodes < xyz(3) + radius - th;
0370         farel = farel & all(nodes_test,2);
0371         if ~any(farel), return, end;
0372     end
0373     idx = find(farel);
0374 end
0375 
0376 
0377 %-------------------------------------------------------------------------%
0378 % Parse options
0379 function [xyz, radius, opt] = parse_opt(xyz, radius, opt)
0380 
0381     scale_radius = false;
0382     if isempty(radius)
0383         radius = xyz(end,:);
0384         scale_radius = true;
0385         xyz(end,:) = [];
0386     end
0387     
0388     if isfield(opt,'desired_img_radius')
0389        scale_radius = false;
0390        if isnumeric(opt.desired_img_radius)
0391           radius = opt.desired_img_radius;
0392        end
0393        if isa(opt.desired_img_radius, 'function_handle')
0394           radius = feval(opt.desired_img_radius,xyz);
0395        end
0396     end
0397         
0398  
0399     
0400     mdl = opt.rec_model; % must exist
0401     opt.n_dim =mdl_dim(mdl);
0402     xyz = xyz(1:opt.n_dim, :); % ingore z if model is 2D
0403 
0404     
0405 
0406     opt.meshsz = [];
0407     try
0408         for i = 1:3
0409             opt.meshsz = [opt.meshsz min(mdl.nodes(:,i)) max(mdl.nodes(:,i))];
0410         end
0411     end
0412     
0413     opt.n_dim = length(opt.meshsz)/2;
0414     opt.meshsz = reshape(opt.meshsz,2,[])';
0415 
0416     opt.minpt = opt.meshsz(:,1);
0417     opt.maxpt = opt.meshsz(:,2);
0418     opt.range = opt.maxpt - opt.minpt;
0419     opt.maxrange = max(opt.range(1:2));
0420 
0421     %rescale points to between -1 and 1 in the x-y plane
0422     xyz = 2*bsxfun(@minus, xyz, mean(opt.meshsz,2))/opt.maxrange;
0423     mdl.nodes = 2*bsxfun(@minus, mdl.nodes, mean(opt.meshsz(1:size(mdl.nodes,2),:),2)')/opt.maxrange;
0424     opt.rec_model = mdl;
0425     if scale_radius
0426         radius = 2*radius / opt.maxrange;
0427     end
0428 
0429     if ~isfield(opt, 'steepness')
0430         opt.steepness = 10./radius;
0431     elseif isa(opt.steepness,'function_handle')
0432         opt.steepness = feval(opt.steepness,xyz);
0433     end
0434     
0435     if numel(opt.steepness) == 1
0436        opt.steepness = ones(1,size(xyz,2)) * opt.steepness;
0437     end
0438     
0439     if numel(radius) == 1
0440         radius = ones(1,size(xyz,2)) * radius;
0441     end
0442     
0443     if opt.n_dim == 2
0444        xyz(3,:) = 0;
0445     end
0446    
0447     if ~isfield(opt, 'threshold')
0448        opt.threshold = 1e-4;
0449     end
0450     
0451     opt.threshold = log(1/opt.threshold - 1); 
0452     
0453 end
0454 
0455 %-------------------------------------------------------------------------%
0456 % UNIT_TEST
0457 function do_unit_test
0458     eidors_cache clear_name GREIT_desired_img_sigmoid
0459     v = linspace(-1,1,32);
0460     mdl= mk_grid_model([],v,v,[0 .7 1:.2:2 2.3 3]);
0461     opt.rec_model = mdl;
0462     opt.steepness = @(xyz) 50./xyz(3,:);
0463     opt.desired_img_radius = @(xyz) xyz(3,:)/5;
0464 %     xyzr = zeros(5,4);
0465     xyzr(:,3) = .5:.5:2.5;
0466     xyzr(:,4) = .25;
0467     sol = GREIT_desired_img_sigmoid(xyzr',[],opt);
0468     img = mk_image(mdl,0);
0469     for i = 1:5
0470         subplot(2,3,i)
0471         img.elem_data= sol(:,i);
0472         show_3d_slices(img,xyzr(i,3),xyzr(i,2),xyzr(i,1));
0473     end
0474     eidors_msg('UNIT_TEST: Showed %d images to verify',i,0);
0475     eidors_cache on GREIT_desired_img_sigmoid
0476 end
GREIT_desired_img_sigmoid

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
