0001 function [fwd_mdl, mat_idx_reordered]= ...
0002 ng_mk_fwd_model( ng_vol_filename, centres, ...
0003 name, stim_pattern, z_contact, postprocmesh)
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033 if isempty(name);
0034 name = ['MDL from', ng_vol_filename];
0035 end
0036
0037 if nargin<5
0038 z_contact=0.01;
0039 end
0040
0041
0042 [srf,vtx,fc,bc,simp,edg,mat_ind] = ng_read_mesh(ng_vol_filename);
0043 if isempty(vtx);
0044 error('EIDORS: ng_mk_fwd_model: Netgen meshing failed. Stopping');
0045 end
0046 if nargin>=6
0047 N_elec = max(size(centres));
0048 [srf,vtx,fc,bc,simp,edg,mat_ind] = feval(postprocmesh,...
0049 srf,vtx,fc,bc,simp,edg,mat_ind, N_elec);
0050 end
0051 fwd_mdl= construct_fwd_model(srf,vtx,simp,bc, name, ...
0052 stim_pattern, centres, z_contact,fc);
0053
0054 [fwd_mdl.mat_idx, fwd_mdl.mat_idx_reordered] = mk_mat_indices(mat_ind);
0055
0056 mat_idx_reordered = fwd_mdl.mat_idx_reordered;
0057
0058 if ~isfield(fwd_mdl,'normalize_measurements')
0059 fwd_mdl.normalize_measurements = 0;
0060 end
0061
0062
0063 function fwd_mdl= construct_fwd_model(srf,vtx,simp,bc, name, ...
0064 stim_pattern, centres, z_contact,fc)
0065 mdl.nodes = vtx;
0066 mdl.elems = simp;
0067 mdl.boundary = srf;
0068 mdl.boundary_numbers=fc;
0069 mdl.gnd_node= find_centre_node(vtx);
0070 mdl.np_fwd_solve.perm_sym = '{n}';
0071 mdl.name = name;
0072
0073
0074 if ~isempty(stim_pattern)
0075 mdl.stimulation= stim_pattern;
0076 end
0077
0078 nelec= size(centres,1);
0079 if nelec>0
0080
0081 [elec,sels,electrodes] = ng_tank_find_elec(srf,vtx,bc,centres);
0082 if size(elec,1) ~= nelec
0083 error('Failed to find all the electrodes')
0084 end
0085
0086
0087 z_contact= z_contact.*ones(nelec,1);
0088 for i=1:nelec
0089 electrodes(i).z_contact= z_contact(i);
0090 end
0091
0092 mdl.electrode = electrodes;
0093 end
0094
0095 mdl.solve= 'eidors_default';
0096 mdl.jacobian= 'eidors_default';
0097 mdl.system_mat= 'eidors_default';
0098
0099 fwd_mdl= eidors_obj('fwd_model', mdl);
0100
0101
0102
0103
0104
0105
0106 function [mat_idx, mat_idx_reordered] = mk_mat_indices( mat_ind);
0107
0108
0109
0110 if isempty(mat_ind)
0111 mat_idx = [];
0112 mat_idx_reordered = [];
0113 return
0114 end
0115 mat_indices = unique( mat_ind );
0116 for i= 1:length(mat_indices);
0117 mat_idx{i}= find(mat_ind == mat_indices(i));
0118 end
0119
0120
0121
0122 for i= 1:length(mat_indices);
0123 mat_idx_l(i) = length( mat_idx{i} );
0124 end
0125 [jnk, max_l] = max(mat_idx_l);
0126 new_idx = 1:length(mat_indices); new_idx(max_l) = [];
0127 ver= eidors_obj('interpreter_version');
0128 if ~ver.isoctave && ver.ver < 7
0129
0130 mat_idx_reordered = {};
0131 ii=1;
0132 for i=[max_l, new_idx];
0133 mat_idx_reordered{ii} = mat_idx{i};
0134 ii=ii+1;
0135 end
0136 else
0137 mat_idx_reordered = cell(length(mat_idx),1);
0138 [mat_idx_reordered{:}] = mat_idx{[max_l, new_idx]};
0139 end
0140
0141 function gnd_node= find_centre_node(vtx);
0142
0143 d = sum( vtx.^2, 2);
0144 [jnk,gnd_node] = min(d);
0145 gnd_node= gnd_node(1);