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 => Vertical 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);

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

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 => Vertical 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 % NOTE: mk_voxel_volume is slow and uses lots of memory
0029 %   To reduce the memory footprint, try
0030 %   vopt.save_memory = 1; %(or try 10)
0031 
0032 % (C) 2015-2018 Andy Adler and Bartlomiej Grychtol
0033 % License: GPL version 2 or version 3
0034 % $Id: GREIT3D_distribution.m 6327 2022-04-21 13:50:55Z aadler $
0035 
0036 if ischar(fmdl) && strcmp(fmdl,'UNIT_TEST'); do_unit_test; return; end
0037 
0038    imdl = select_imdl(fmdl,{'Basic GN dif'});
0039    imdl = mk_voxel_volume(imdl,vopt);
0040    msm = imdl.rec_model.mdl_slice_mapper;
0041    imdl.rec_model.mdl_slice_mapper.y_pts = ...
0042        fliplr(imdl.rec_model.mdl_slice_mapper.y_pts);
0043    [x, y, z] = ndgrid(msm.x_pts, msm.y_pts, msm.z_pts);
0044 
0045    IN = point_in_vox(imdl.rec_model,[x(:),y(:),z(:)]);
0046    INdesity = nnz(IN)/prod(size(IN));
0047    if INdesity==0
0048       error('GREIT3D_distribution: grid doesn''t match model');
0049    elseif INdesity<0.5
0050       eidors_msg('GREIT3D_distribution: grid only matches %3.1f%% of model', INdesity*100, 1);
0051    end
0052    distr = [x(IN),y(IN),z(IN)]';
0053 
0054 function out = point_in_vox(fmdl,points, epsilon)
0055 
0056    if nargin < 3
0057       epsilon = 0;
0058    end
0059 
0060    levels = unique(fmdl.nodes(:,3));
0061    out = false(numel(points)/3,1);
0062 
0063    opt.faces = true;
0064    fmdl = fix_model(fmdl, opt);
0065 
0066    jnk.nodes = fmdl.nodes(:,1:2);
0067    jnk.type = 'fwd_model';
0068 
0069    for i = 1:numel(levels)-1
0070       idx = (points(:,3) - levels(i)) > epsilon;
0071       idx = idx & ((levels(i+1) - points(:,3)) > epsilon);
0072       
0073       nidx = fmdl.nodes(:,3) == levels(i);
0074       F = all(nidx(fmdl.faces),2);
0075       jnk.elems = fmdl.faces(F,:);
0076       bnd = find_boundary(jnk);
0077       bnd = order_loop(jnk.nodes(unique(bnd(:)),:));
0078       out(idx) = out(idx) | inpolygon(points(idx,1),points(idx,2),bnd(:,1),bnd(:,2));
0079    end
0080 
0081 function do_small_test
0082    fmdl= ng_mk_cyl_models([4,1],[16,1.5,2.5],[0.10]);
0083    fmdl= mystim_square(fmdl);
0084    vopt.imgsz = [16 16];
0085    vopt.zvec = linspace( 0.5,3.5,4);
0086    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0087    unit_test_cmp('Coarse1 size:', size(opt.distr),[3,672]);
0088    unit_test_cmp('Corase1 vals:', opt.distr(:,1),[-0.4375;-0.9375;1],1e-7);
0089 
0090    vopt.zvec = linspace(-1,1,3);
0091    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0092    unit_test_cmp('Coarse2 size:', size(opt.distr),[3,224]);
0093    unit_test_cmp('Corase2 vals:', opt.distr(:,1),[-0.4375;-0.9375;0.5],1e-7);
0094 
0095    vopt.zvec = linspace(-2,0,3);
0096    try
0097        get_error = 0;
0098        [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0099    catch
0100        get_error = 1;
0101    end
0102    unit_test_cmp('Coarse3: got error', get_error,1)
0103 
0104 function do_unit_test
0105    do_small_test
0106    [vh,vi] = test_fwd_solutions;
0107    % inverse geometry
0108    fmdl= ng_mk_cyl_models([4,1,.5],[16,1.5,2.5],[0.05]);
0109    fmdl= mystim_square(fmdl);
0110 
0111    vopt.imgsz = [32 32];
0112    vopt.zvec = linspace( 0.5,3.5,7);
0113    vopt.save_memory = 1;
0114 %  opt.noise_figure = 2;
0115    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0116    n_weight = 18.1608; % approx NF=2
0117    imdl = mk_GREIT_model(imdl, 0.2, n_weight, opt);
0118 
0119    unit_test_cmp('RM size', size(imdl.solve_use_matrix.RM), [5136,928]);
0120    tst = imdl.solve_use_matrix.RM(1,1:2);
0121    tst_= [8.572908305490031 -18.600780994066596];
0122    unit_test_cmp('RM', tst, tst_, 1e-10);
0123 
0124    img = inv_solve(imdl, vh, vi);
0125    unit_test_cmp('img size', size(img.elem_data), [5136,5]);
0126    [mm,ll] =max(img.elem_data(:,1));
0127    unit_test_cmp('img', [mm,ll], ...
0128        [0.582158646727019, 1308], 1e-10);
0129 
0130    img.show_slices.img_cols = 1;
0131    subplot(131); show_fem(fmdl); title 'fmdl'
0132    subplot(132); show_slices(img,[inf,inf,2]); title 'centre';
0133    subplot(133); show_slices(img,[inf,inf,1.5]); title 'bottom';
0134    
0135 
0136 function fmdl = mystim_square(fmdl);
0137    [~,fmdl] = elec_rearrange([16,2],'square', fmdl);
0138    [fmdl.stimulation, fmdl.meas_select] =  ...
0139        mk_stim_patterns(32,1,[0,5],[0,5],{},1);
0140 
0141 function [vh,vi] = test_fwd_solutions;
0142    posns= linspace(1.0,3.0,5);
0143    str=''; for i=1:length(posns);
0144       extra{i} = sprintf('ball%03d',round(posns(i)*100));
0145       str = [str,sprintf('solid %s=sphere(0.5,0,%f;0.1); ', extra{i}, posns(i))];
0146    end
0147    extra{i+1} = str;
0148    fmdl= ng_mk_cyl_models([4,1,.2],[16,1.5,2.5],[0.05],extra); 
0149    fmdl = mystim_square(fmdl);
0150    
0151    img= mk_image(fmdl,1);
0152    vh = fwd_solve(img);
0153    %vh = add_noise(2000, vh);
0154    for i=1:length(posns);
0155       img= mk_image(fmdl,1);
0156       img.elem_data(fmdl.mat_idx{i+1}) = 2;
0157       vi{i} = fwd_solve(img);
0158    end;
0159    vi = [vi{:}];