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 5791 2018-05-22 15:16:02Z 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    distr = [x(IN),y(IN),z(IN)]';
0047 
0048 function out = point_in_vox(fmdl,points, epsilon)
0049 
0050    if nargin < 3
0051       epsilon = 0;
0052    end
0053 
0054    levels = unique(fmdl.nodes(:,3));
0055    out = false(numel(points)/3,1);
0056 
0057    opt.faces = true;
0058    fmdl = fix_model(fmdl, opt);
0059 
0060    jnk.nodes = fmdl.nodes(:,1:2);
0061    jnk.type = 'fwd_model';
0062 
0063    for i = 1:numel(levels)-1
0064       idx = (points(:,3) - levels(i)) > epsilon;
0065       idx = idx & ((levels(i+1) - points(:,3)) > epsilon);
0066       
0067       nidx = fmdl.nodes(:,3) == levels(i);
0068       F = all(nidx(fmdl.faces),2);
0069       jnk.elems = fmdl.faces(F,:);
0070       bnd = find_boundary(jnk);
0071       bnd = order_loop(jnk.nodes(unique(bnd(:)),:));
0072       out(idx) = out(idx) | inpolygon(points(idx,1),points(idx,2),bnd(:,1),bnd(:,2));
0073    end
0074 
0075 function do_unit_test
0076    [vh,vi] = test_fwd_solutions;
0077    % inverse geometry
0078    fmdl= ng_mk_cyl_models([4,1,.5],[16,1.5,2.5],[0.05]);
0079    fmdl= mystim_square(fmdl);
0080 
0081    vopt.imgsz = [32 32];
0082    vopt.zvec = linspace( 0.5,3.5,7);
0083    vopt.save_memory = 1;
0084    opt.noise_figure = 2;
0085    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0086    imdl = mk_GREIT_model(imdl, 0.2, [], opt);
0087 
0088    unit_test_cmp('RM size', size(imdl.solve_use_matrix.RM), [5136,928]);
0089    unit_test_cmp('RM', imdl.solve_use_matrix.RM(1,1:2), ...
0090        [8.573008430710381 -18.601036109804095], 1e-10);
0091 
0092    img = inv_solve(imdl, vh, vi);
0093    unit_test_cmp('img size', size(img.elem_data), [5136,5]);
0094    [mm,ll] =max(img.elem_data(:,1));
0095    unit_test_cmp('img', [mm,ll], ...
0096        [0.582161052175761, 1308], 1e-10);
0097 
0098    img.show_slices.img_cols = 1;
0099    subplot(131); show_fem(fmdl); title 'fmdl'
0100    subplot(132); show_slices(img,[inf,inf,2]); title 'centre';
0101    subplot(133); show_slices(img,[inf,inf,1.5]); title 'bottom';
0102    
0103 
0104 function fmdl = mystim_square(fmdl);
0105    [~,fmdl] = elec_rearrange([16,2],'square', fmdl);
0106    [fmdl.stimulation, fmdl.meas_select] =  ...
0107        mk_stim_patterns(32,1,[0,5],[0,5],{},1);
0108 
0109 function [vh,vi] = test_fwd_solutions;
0110    posns= linspace(1.0,3.0,5);
0111    str=''; for i=1:length(posns);
0112       extra{i} = sprintf('ball%03d',round(posns(i)*100));
0113       str = [str,sprintf('solid %s=sphere(0.5,0,%f;0.1); ', extra{i}, posns(i))];
0114    end
0115    extra{i+1} = str;
0116    fmdl= ng_mk_cyl_models([4,1,.2],[16,1.5,2.5],[0.05],extra); 
0117    fmdl = mystim_square(fmdl);
0118    
0119    img= mk_image(fmdl,1);
0120    vh = fwd_solve(img);
0121    %vh = add_noise(2000, vh);
0122    for i=1:length(posns);
0123       img= mk_image(fmdl,1);
0124       img.elem_data(fmdl.mat_idx{i+1}) = 2;
0125       vi{i} = fwd_solve(img);
0126    end;
0127    vi = [vi{:}];