GREIT3D_distribution

 GREIT3D_distribution: create target distributions for 3D GREIT

 The basic usage is part of creating a 3D GREIT model, 
  given a forward model fmdl
    vopt.imgsz = [32 32 19]; % and other parameters
    opt.noise_figure = NF;   % and other GREIT options
    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
    imdl = mk_GREIT_model(imdl, 0.20, [], opt);

 INPUTS:
  fmdl => EIDORS forward model structure
  vopt => Volume options (see mk_voxel_volume help)

 Options for vopt
   1: Simple
        vopt.imgsz = [32 32 19];
        vopt.square_pixels = true;
   2: Specify in detail the z planes (e.g.)
        vopt.imgsz = [32 32];
        vopt.zvec = linspace( 0.5,3.5,15);
          (Vertical size is: length(vopt.zvec)-1
   3: Specify x,y and z planes (e.g.)
        vopt.xvec = linspace(-.95,.95,32);
        vopt.yvec = linspace(-.95,.95,32);
        vopt.zvec = linspace( .1,3.9,11);

 Reducing the density of target distributions:
   By default, GREIT3D_distribution places one target in each voxel. The
   density can be decreased by downsampling:
       vopt.downsample = [N, PHASE]; 
   for uniform downsampling, or
       vopt.downsample = [xN, xPHASE; yN, yPHASE; zN, zPHASE];
   for different factors in each dimension. PHASE can be omitted and 
   defaults to 0.
 Downsampling returns every Nth element, starting from 1+PHASE.

 NOTE: mk_voxel_volume is slow and uses lots of memory
   To reduce the memory footprint, try
   vopt.save_memory = 1; %(or try 10)

 See also MK_VOXEL_VOLUME, MK_GREIT_MODEL, DOWNSAMPLE

0001 function [imdl,distr] = GREIT3D_distribution(fmdl, vopt) 
0002 % GREIT3D_distribution: create target distributions for 3D GREIT
0003 %
0004 % The basic usage is part of creating a 3D GREIT model,
0005 %  given a forward model fmdl
0006 %    vopt.imgsz = [32 32 19]; % and other parameters
0007 %    opt.noise_figure = NF;   % and other GREIT options
0008 %    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0009 %    imdl = mk_GREIT_model(imdl, 0.20, [], opt);
0010 %
0011 % INPUTS:
0012 %  fmdl => EIDORS forward model structure
0013 %  vopt => Volume options (see mk_voxel_volume help)
0014 %
0015 % Options for vopt
0016 %   1: Simple
0017 %        vopt.imgsz = [32 32 19];
0018 %        vopt.square_pixels = true;
0019 %   2: Specify in detail the z planes (e.g.)
0020 %        vopt.imgsz = [32 32];
0021 %        vopt.zvec = linspace( 0.5,3.5,15);
0022 %          (Vertical size is: length(vopt.zvec)-1
0023 %   3: Specify x,y and z planes (e.g.)
0024 %        vopt.xvec = linspace(-.95,.95,32);
0025 %        vopt.yvec = linspace(-.95,.95,32);
0026 %        vopt.zvec = linspace( .1,3.9,11);
0027 %
0028 % Reducing the density of target distributions:
0029 %   By default, GREIT3D_distribution places one target in each voxel. The
0030 %   density can be decreased by downsampling:
0031 %       vopt.downsample = [N, PHASE];
0032 %   for uniform downsampling, or
0033 %       vopt.downsample = [xN, xPHASE; yN, yPHASE; zN, zPHASE];
0034 %   for different factors in each dimension. PHASE can be omitted and
0035 %   defaults to 0.
0036 % Downsampling returns every Nth element, starting from 1+PHASE.
0037 %
0038 % NOTE: mk_voxel_volume is slow and uses lots of memory
0039 %   To reduce the memory footprint, try
0040 %   vopt.save_memory = 1; %(or try 10)
0041 %
0042 % See also MK_VOXEL_VOLUME, MK_GREIT_MODEL, DOWNSAMPLE
0043 
0044 % (C) 2015-2018 Andy Adler and Bartlomiej Grychtol
0045 % License: GPL version 2 or version 3
0046 % $Id: GREIT3D_distribution.m 7013 2024-11-26 20:31:32Z bgrychtol $
0047 
0048 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0049 
0050    imdl = select_imdl(fmdl,{'Basic GN dif'});
0051    imdl = mk_voxel_volume(imdl,vopt);
0052    msm = imdl.rec_model.mdl_slice_mapper;
0053    imdl.rec_model.mdl_slice_mapper.y_pts = ...
0054        fliplr(imdl.rec_model.mdl_slice_mapper.y_pts);
0055    
0056    [xpts, ypts, zpts] = downsample_distr(msm.x_pts, msm.y_pts, msm.z_pts, vopt);
0057 
0058    [x, y, z] = ndgrid(xpts, ypts, zpts);
0059 
0060 
0061    IN = point_in_vox(imdl.rec_model,[x(:),y(:),z(:)]);
0062    INdesity = nnz(IN)/numel(IN);
0063    if INdesity==0
0064       error('GREIT3D_distribution: grid doesn''t match model');
0065    elseif INdesity<0.5
0066       eidors_msg('GREIT3D_distribution: grid only matches %3.1f%% of model', INdesity*100, 1);
0067    end
0068    distr = [x(IN),y(IN),z(IN)]';
0069 
0070 function out = point_in_vox(fmdl,points, epsilon)
0071 
0072    if nargin < 3
0073       epsilon = 0;
0074    end
0075 
0076    levels = unique(fmdl.nodes(:,3));
0077    out = false(numel(points)/3,1);
0078 
0079    opt.faces = true;
0080    fmdl = fix_model(fmdl, opt);
0081 
0082    jnk.nodes = fmdl.nodes(:,1:2);
0083    jnk.type = 'fwd_model';
0084 
0085    for i = 1:numel(levels)-1
0086       idx = (points(:,3) - levels(i)) > epsilon;
0087       idx = idx & ((levels(i+1) - points(:,3)) > epsilon);
0088       
0089       nidx = fmdl.nodes(:,3) == levels(i);
0090       F = all(nidx(fmdl.faces),2);
0091       jnk.elems = fmdl.faces(F,:);
0092       bnd = find_boundary(jnk);
0093       bnd = order_loop(jnk.nodes(unique(bnd(:)),:));
0094       out(idx) = out(idx) | inpolygon(points(idx,1),points(idx,2),bnd(:,1),bnd(:,2));
0095    end
0096 
0097 function [xpts, ypts, zpts] = downsample_distr(xpts, ypts, zpts, vopt)
0098    if ~isfield(vopt, 'downsample')
0099       return
0100    end
0101 
0102    P = zeros(3,1); % phase
0103    N = ones(3,1);  % downsample factor
0104    switch numel(vopt.downsample)
0105       case {1, 2} % same for all dimensions
0106          N(:) = vopt.downsample(1);
0107          try P(:) = vopt.downsample(2); end
0108       case 3
0109          N(:) = vopt.downsample(:);
0110       case 6
0111          N(:) = vopt.downsample(:,1);
0112          P(:) = vopt.downsample(:,2);
0113       otherwise
0114          error('EIDORS:WrongInput','vopt.downsample is wrong size.')
0115    end
0116    xpts = downsample(xpts,N(1),P(1));
0117    ypts = downsample(ypts,N(2),P(2));
0118    zpts = downsample(zpts,N(3),P(3));
0119 
0120       
0121 function pts = downsample(pts, N, P)
0122 % Implements downsample, which is part of the signal processing toolbox
0123 % Coded independently, without consulting Mathwork's sourcecode
0124 if nargin < 3, P = 0; end
0125 pts = pts(1+P:N:end);
0126         
0127 
0128 
0129 function do_small_test
0130    fmdl= ng_mk_cyl_models([4,1],[16,1.5,2.5],[0.10]);
0131    fmdl= mystim_square(fmdl);
0132    vopt.imgsz = [16 16];
0133    vopt.zvec = linspace( 0.5,3.5,4);
0134    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0135    unit_test_cmp('Coarse1 size:', size(opt.distr),[3,672]);
0136    unit_test_cmp('Corase1 vals:', opt.distr(:,1),[-0.4375;-0.9375;1],1e-7);
0137 
0138    vopt.zvec = linspace(-1,1,3);
0139    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0140    unit_test_cmp('Coarse2 size:', size(opt.distr),[3,224]);
0141    unit_test_cmp('Corase2 vals:', opt.distr(:,1),[-0.4375;-0.9375;0.5],1e-7);
0142 
0143    vopt.zvec = linspace(-2,0,3);
0144    try
0145        get_error = 0;
0146        [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0147    catch
0148        get_error = 1;
0149    end
0150    unit_test_cmp('Coarse3: got error', get_error,1)
0151 
0152 function do_unit_test
0153    do_small_test
0154    [vh,vi] = test_fwd_solutions;
0155    % inverse geometry
0156    fmdl= ng_mk_cyl_models([4,1,.5],[16,1.5,2.5],[0.05]);
0157    fmdl= mystim_square(fmdl);
0158 
0159    vopt.imgsz = [32 32];
0160    vopt.zvec = linspace( 0.5,3.5,7);
0161    vopt.save_memory = 1;
0162 %  opt.noise_figure = 2;
0163    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0164    n_weight = 18.1608; % approx NF=2
0165    imdl = mk_GREIT_model(imdl, 0.2, n_weight, opt);
0166 
0167    unit_test_cmp('RM size', size(imdl.solve_use_matrix.RM), [5136,928]);
0168    tst = imdl.solve_use_matrix.RM(1,1:2);
0169    % numbers below depend on the version of Netgen
0170    tst_= [8.572908305490031 -18.600780994066596]; % Linux, 2022
0171    tst_= [8.572908283730975 -18.600781037986948]; % ng 5.3, Windows, 2024
0172    unit_test_cmp('RM', tst, tst_, 1e-7);
0173 
0174    img = inv_solve(imdl, vh, vi);
0175    unit_test_cmp('img size', size(img.elem_data), [5136,5]);
0176    [mm,ll] =max(img.elem_data(:,1));
0177    tst_ = [0.582158646727019 1308]; % Linux, 2022
0178    tst_ = [0.582158646971746 1308]; % ng 5.3, Windows, 2024
0179    unit_test_cmp('img', [mm,ll], ...
0180        [0.582158646727019, 1308], 1e-8);
0181 
0182    img.show_slices.img_cols = 1;
0183    subplot(131); show_fem(fmdl); title 'fmdl'
0184    subplot(132); show_slices(img,[inf,inf,2]); title 'centre';
0185    subplot(133); show_slices(img,[inf,inf,1.5]); title 'bottom';
0186    
0187 
0188 function fmdl = mystim_square(fmdl);
0189    [~,fmdl] = elec_rearrange([16,2],'square', fmdl);
0190    [fmdl.stimulation, fmdl.meas_select] =  ...
0191        mk_stim_patterns(32,1,[0,5],[0,5],{},1);
0192 
0193 function [vh,vi] = test_fwd_solutions;
0194    posns= linspace(1.0,3.0,5);
0195    str=''; for i=1:length(posns);
0196       extra{i} = sprintf('ball%03d',round(posns(i)*100));
0197       str = [str,sprintf('solid %s=sphere(0.5,0,%f;0.1); ', extra{i}, posns(i))];
0198    end
0199    extra{i+1} = str;
0200    fmdl= ng_mk_cyl_models([4,1,.2],[16,1.5,2.5],[0.05],extra); 
0201    fmdl = mystim_square(fmdl);
0202    
0203    img= mk_image(fmdl,1);
0204    vh = fwd_solve(img);
0205    %vh = add_noise(2000, vh);
0206    for i=1:length(posns);
0207       img= mk_image(fmdl,1);
0208       img.elem_data(fmdl.mat_idx{i+1}) = 2;
0209       vi{i} = fwd_solve(img);
0210    end;
0211    vi = [vi{:}];