Description of mk_thorax_model

0001 function out = mk_thorax_model_bp3d(varargin)
0002 %MK_THORAX_MODEL_BP3D Torso model based on BodyParts3D
0003 %
0004 % FMDL = MK_THORAX_MODEL_BP3D()
0005 % Returns a torso fwd_model with no electrodes or organs
0006 %
0007 % FMDL = MK_THORAX_MODEL_BP3D(elec_pos, elec_shape)
0008 % FMDL = MK_THORAX_MODEL_BP3D(elec_pos, elec_shape, maxh)
0009 % Allows specifying electrode configuration, see PLACE_ELEC_ON_SURF
0010 %
0011 % FMDL = MK_THORAX_MODEL_BP3D(stimpat)
0012 % A predefined model with round electrodes (radius 1 cm). Accepts a string
0013 % parameter to define a stim pattern:
0014 %    '2x16_planar'   - 2 rings of 16 electrodes, planar stimulation
0015 %    '2x16_odd-even' - 2 rings, each stimulation is cross-plane
0016 %    '2x16_square'   - 2 rings, every second stimulation is cross-plane
0017 %    '2x16_adjacent' - like square, but with adjacent stimulation
0018 %    '1x32_ring'     - single ring of 32 electrodes
0019 %
0020 % IMG = MK_THORAX_MODEL_BP3D(... , organspec)
0021 % Returns a torso model containing organ shapes as an image structure. Can
0022 % be used with all the above electrode defintions.
0023 % ORGANSPEC must be one of:
0024 %    'fine'          - a relatively dense model, ~7 million elems
0025 %    'coarse'        - a coarser model, *NOT IMPLEMENTED YET*
0026 %
0027 %
0028 % The 'fine' model contains the following materials:
0029 %    #   mat_names         Description
0030 %   -----------------------------------------------------------------------
0031 %     1    Body                 Volume between the skin and the internal organs
0032 %     2    LungLeft             Left lung
0033 %     3    LungRight            Right lung
0034 %     4    HeartWall            Heart wall
0035 %     5    Aorta                Aorta (wall)
0036 %     6    SVC                  Superior Vena Cava (wall)
0037 %     7    IVC                  Inferior Vena Cava (wall)
0038 %     8    PA                   Pulmonary Artery (wall)
0039 %     9    LSPV                 Left Superior Pulmonary Vein (wall)
0040 %    10    LIPV                 Left Inferior Pulmonary Vein (wall)
0041 %    11    RSPV                 Right Superior Pulmonary Vein (wall)
0042 %    12    RIPV                 Right Inferior Pulmonary Vein (wall)
0043 %    13    BloodHeartLeft       Blood in left atrium and left ventricle
0044 %    14    BloodHeartRight      Blood in right atrium and right ventricle
0045 %    15    BloodAorta           Aorta (lumen)
0046 %    16    BloodSVC             Superior Vena Cava (lumen)
0047 %    17    BloodIVC             Inferior Vena Cava (lumen)
0048 %    18    BloodPA              Pulmonary Artery (lumen)
0049 %    19    BloodLSPV            Left Superior Pulmonary Vein (lumen)
0050 %    20    BloodLIPV            Left Inferior Pulmonary Vein (lumen)
0051 %    21    BloodRSPV            Right Superior Pulmonary Vein (lumen)
0052 %    22    BloodRIPV            Right Inferior Pulmonary Vein (lumen)
0053 %
0054 %
0055 % CITATION_REQUEST:
0056 % AUTHOR: Grychtol B
0057 % TITLE: Torso Organ Meshes
0058 % JOURNAL: Zenodo
0059 % YEAR: 2024
0060 % DOI: 10.5281/zenodo.11047831
0061 %
0062 % LICENSE INFORMATION:
0063 % The model generated by this function is a derivative of the <a href=
0064 % "matlab: web https://https://doi.org/10.5281/zenodo.11047830 -browser"
0065 % >"Torso Organ Meshes"</a>
0066 % by Bartłomiej Grychtol, itself a derivative of the <a href="matlab: web
0067 % https://dbarchive.biosciencedbc.jp/en/bodyparts3d -browser">BodyParts3D</a> database by
0068 % <a href="matlab: web https://biosciencedbc.jp/ -browser"
0069 % >The Database Center for Life Science</a>, used under the Creative Commons
0070 % Attribution-ShareAlike 2.1 Japan (<a href="matlab: web
0071 % http://creativecommons.org/licenses/by-sa/2.1/jp/deed.en_US
0072 % -browser">CC BY-SA 2.1 JP</a>) license.
0073 % The model must be licensed under CC BY-SA 2.1 JP and attributed
0074 % correctly.
0075 %
0076 % See also MK_THORAX_MODEL, PLACE_ELEC_ON_SURF
0077 
0078 % (C) 2024 Bartek Grychtol. License: GPL v2 or v3
0079 % $Id: mk_thorax_model_bp3d.m 7105 2024-12-26 20:49:00Z bgrychtol $
0080 
0081 citeme(mfilename)
0082 
0083 if nargin == 0 % get a torso with no electrodes
0084    out = get_torso; 
0085 end
0086 organspec = '';
0087 copyright = '';
0088 if nargin > 1 && ischar(varargin{end})
0089    organspec = varargin{end};
0090    varargin(end) = [];
0091 end
0092 
0093 switch length(varargin)
0094    case 0
0095       out = mk_organ_image(out, organspec); % for the sake of completeness
0096    case 1
0097       if isempty(varargin{1}) 
0098          out = get_torso;
0099          if ~isempty(organspec)
0100             out = mk_organ_image(out, organspec);
0101          end
0102       elseif ischar(varargin{1}) 
0103          copyright = '(C) EIDORS Project';
0104          out = build_predef_elecs(varargin{1}, organspec);
0105       else
0106          error('EIDORS:WrongInput','Electrode specification not understood');
0107       end
0108    case {2,3}
0109       out = mk_thorax_model('bp3d',varargin{:});
0110       out = mk_organ_image(out, organspec);
0111    otherwise
0112       error('EIDORS:WrongInput', 'Too many inputs');
0113 end
0114 
0115 out = eidors_obj('set', out); % sort fields
0116 
0117 if ~isempty(copyright)
0118    out.copyright = copyright;
0119 end
0120 
0121 out.attribution = ...
0122    ['Modified from "Torso Organ Meshes" ' ...
0123     '(https://doi.org/10.5281/zenodo.11047831) by Bartłomiej Grychtol ' ...
0124     '/ A derivative of "BodyParts3D" '...
0125     '(https://dbarchive.biosciencedbc.jp/en/bodyparts3d) by '...
0126     'The Database Center for Life Science.'];
0127 
0128 out.license   =  ['Creative Commons Attribution-ShareAlike 2.1 Japan '...;
0129        '(http://creativecommons.org/licenses/by-sa/2.1/jp/deed.en_US)'];
0130 
0131 if strcmp(out.type, 'image')
0132    try out.fwd_model.copyright = out.copyright; end
0133    out.fwd_model.attribution = out.attribution;
0134    out.fwd_model.license = out.license;
0135 end
0136 
0137 ws = warning('off', 'backtrace');
0138 warning('EIDORS:License', ...
0139    ['The bp3d torso model is licensed under CC BY-SA 2.1 JP.\n' ...
0140     'Please make sure to attribute correctly. ' ...
0141     'See the attribution field for details.']);
0142 warning(ws.state, 'backtrace');
0143 
0144 
0145 function out = get_torso()
0146 contrib = 'bp3d-torso'; file =  'bp3d_torso.mat';
0147 load(get_contrib_data(contrib,file),'torso');
0148 eidors_msg('Building bp3d body surface', 2);
0149 out = gmsh_stl2tet(torso);
0150 out = fix_boundary(out);
0151 out.name = 'BP3D Torso Model';
0152 
0153 
0154 function out = build_predef_elecs(stimpat, organspec)
0155 cache_path = get_cache_path;
0156 fname = 'bp3d_predef_torso';
0157 suffix = '';
0158 if ~isempty(organspec)
0159    fname = [fname '_' organspec];
0160    suffix = [' - ' organspec];
0161 end
0162 name = ['BP3D Torso Model' suffix];
0163 fpath = [cache_path filesep fname '.mat'];
0164 if exist(fpath,'file')
0165    eidors_msg('@@ Using stored model.', 3);
0166    load(fpath, 'out');
0167 else
0168    eth16 = 360*cumsum([0.2 0.4*ones(1,7) 0.5 0.4*ones(1,7)])/6.5 - 90; eth16 = eth16';
0169    eth32 = 360*cumsum([0.1 0.2*ones(1,15) 0.25 0.2*ones(1,15)])/6.45 - 90; eth32 = eth32';
0170    ep = eth16; ep(:,2) = 1175;
0171    ep(17:48,1) = eth32; ep(17:48,2) = 1212.5;
0172    ep(49:64,1) = eth16; ep(49:64,2) = 1250;
0173    out = mk_thorax_model('bp3d',ep,[10 0 1],20);
0174    out.name = name;
0175    out = mk_organ_image(out, organspec);
0176    save(fpath,'out');
0177 end
0178 out = set_predef_stim(out, stimpat);
0179 out.name = [out.name ' - ' stimpat];
0180 
0181 
0182 function out = mk_organ_image(out, organspec)
0183 if isempty(organspec), return, end
0184 organs = get_organ_models(organspec);
0185 out = eidors_cache(@add_organs,{out, organs});
0186 out = set_organ_values(out);
0187 out.name = 'BP3D Torso Image';
0188 
0189 
0190 function out = add_organs(torso, organs)
0191 tri.vertices = torso.nodes;
0192 tri.faces = torso.boundary;
0193 N = size(tri.vertices, 1);
0194 tri.vertices = [tri.vertices; organs.nodes];
0195 tri.faces = [tri.faces; N + organs.boundary(:,[1 3 2])];
0196 
0197 eidors_msg('Re-meshing body volume', 2);
0198 shell = gmsh_stl2tet(tri,[],[],4); % lots of optimization
0199 print_mdl_info(shell, 'Body');
0200 
0201 eidors_msg('Restoring electrodes', 2);
0202 if isfield(torso, 'electrode') && ~isempty(torso.electrode)
0203    for i = 1:numel(torso.electrode)
0204       electrode(i).nodes = torso.nodes(torso.electrode(i).nodes,:);
0205    end
0206    shell = restore_electrodes(shell, electrode, 0.01);
0207 end
0208 organs.boundary = organs.show_boundary;
0209 organs = rmfield(organs, 'show_boundary');
0210 
0211 eidors_msg('Merging body and organ volumes', 2);
0212 out = merge_meshes(shell, organs);
0213 out.mat_names = ['Body'; organs.mat_names];
0214  
0215 eidors_msg('Removing trachea', 2);
0216 idx = strcmp(out.mat_names,'Trachea');
0217 out = remove_elems(out, out.mat_idx{idx});
0218 
0219 % nicely re-order materials
0220 idx = [1 5 10 2  17:2:21 15 6 8 11 13 3 4 18:2:22 16 7 9 12 14];
0221 out.mat_idx = out.mat_idx(idx);
0222 out.mat_names = out.mat_names(idx);
0223 
0224 out = eidors_obj('set', out); % order fields
0225 
0226 print_mdl_info(out, 'Final model');
0227 
0228 function img = set_organ_values(torso)
0229 img = mk_image(torso, 1, 'bp3d-torso');
0230 for i = 2:numel(torso.mat_idx)
0231     if strcmp(torso.mat_names{i},'HeartWall')
0232         img.elem_data(torso.mat_idx{i}) = 1.5;
0233 %     elseif strcmp(torso.mat_names{i},'Trachea')
0234 %         img.elem_data(torso.mat_idx{i}) = NaN;
0235     elseif ~isempty(strfind(torso.mat_names{i}, 'Blood'))
0236         img.elem_data(torso.mat_idx{i}) = 2;
0237     elseif ~isempty(strfind(torso.mat_names{i}, 'Lung'))
0238         img.elem_data(torso.mat_idx{i}) = .5;
0239     else % vessel
0240         img.elem_data(torso.mat_idx{i}) = .8;
0241     end
0242 end
0243 img.calc_colours.ref_level = 1;
0244 
0245 function mm = restore_electrodes(mm, electrode, z_contact, thresh)
0246 if nargin < 4
0247    thresh = eps;
0248 end
0249 surf_node_idx = unique(find_boundary(mm));
0250 surf_nodes = mm.nodes(surf_node_idx,:);
0251 % ignore nodes close to the bottom surface
0252 rng = max(surf_nodes(:,3)) - min(surf_nodes(:,3));
0253 idx = surf_nodes(:,3) > min(surf_nodes(:,3)) + 0.02 * rng;
0254 surf_node_idx = surf_node_idx(idx);
0255 surf_nodes = surf_nodes(idx,:);
0256 for e = 1:numel(electrode)
0257    % if we got this far, we should have enough memory for vectorized calcs
0258    dist = (surf_nodes(:,1) - electrode(e).nodes(:,1)').^2;
0259    for d = 2:3
0260       dist = dist + (surf_nodes(:,d) - electrode(e).nodes(:,d)').^2;
0261    end
0262 %    dist = sqrt(dist); % takes a long time and changes nothing
0263    [val, pos] = min(dist);
0264    if any(val > thresh)
0265       warning('Some nodes of electrode %d seem too far. Furthest distance: %f',e,sqrt(max(val)));
0266    end
0267    mm.electrode(e).nodes = surf_node_idx(pos);
0268    mm.electrode(e).z_contact = z_contact;
0269 end
0270 
0271 
0272 function organs = get_organ_models(organspec)
0273 switch organspec
0274    case 'fine'
0275       fhandle = @mk_fine_organs;
0276    case 'coarse'
0277       error('EIDORS:NotImplemented', ['Coarse organ model is not '...
0278             'implemented yet.'])
0279    otherwise
0280       error('EIDORS:WrongInput', 'Organ specification not understood.');
0281 end
0282 cache_path = get_cache_path;
0283 fname = sprintf('bp3d_organs_%s.mat', organspec);
0284 fpath = [cache_path filesep fname];
0285 if ~exist(fpath, 'file')
0286    eidors_msg('This takes a REALLY long time');
0287    organs = feval(fhandle);
0288    save(fpath, 'organs');
0289 else
0290    load(fpath, 'organs');
0291 end
0292 
0293 
0294 function mdl = mk_fine_organs
0295 
0296 start_time = now(); 
0297 organs = mk_individual_organs;
0298 mdl = assemble_organs(organs);
0299 mdl.show_boundary = mdl.boundary;
0300 eidors_msg('Recalculating boundary',2);
0301 mdl = fix_boundary(mdl);
0302 
0303 %  cashed results are not useful outside
0304 eidors_cache('clear_old', start_time);
0305 
0306 
0307 function mdl = assemble_organs(organs)
0308 pcs = struct();
0309 
0310 eidors_msg('Merging vessels', 2);
0311 vessels = {'Aorta','IVC','SVC','PA', 'LSPV', 'LIPV','RSPV','RIPV'};
0312 for i = 1:numel(vessels)
0313     v = vessels{i};
0314     disp(v)
0315     pcs.(v) = merge_meshes(organs.(v), organs.(['Blood' v]),1e-2);
0316     pcs.(v).mat_names = {v; ['Blood' v]};
0317 end
0318 
0319 eidors_msg('Merging left side', 2);
0320 left = merge_meshes(organs.LungLeft, ...            
0321                     pcs.LSPV, ...
0322                     pcs.LIPV, ...
0323                     1e-2);
0324 left.mat_names = vertcat({'LungLeft'},pcs.LSPV.mat_names, pcs.LIPV.mat_names);
0325 
0326 eidors_msg('Merging right side', 2);
0327 right = merge_meshes(organs.LungRight, ...            
0328                      pcs.RSPV, ...
0329                      pcs.RIPV, ...
0330                      2e-1); % need a larger threshold due to mesh mismatch
0331 right.mat_names = vertcat({'LungRight'},pcs.RSPV.mat_names, pcs.RIPV.mat_names);
0332 
0333 eidors_msg('Merging heart', 2);
0334 heart = merge_meshes(organs.HeartWall, ...
0335                      organs.BloodLeftHeart, ...
0336                      organs.BloodRightHeart, ...
0337                      1e-2);
0338 heart.mat_names = {'HeartWall';'BloodHeartLeft'; 'BloodHeartRight'};
0339 
0340 eidors_msg('Merging heart, left and right', 2);
0341 mdl = merge_meshes(heart, left, right, 1e-2);
0342 mdl.mat_names = vertcat(heart.mat_names, left.mat_names, right.mat_names);
0343 
0344 eidors_msg('Merging remaining vessels', 2);
0345 mdl = merge_meshes(mdl, ...
0346                   pcs.PA, ...
0347                   pcs.Aorta, ...
0348                   pcs.SVC, ...
0349                   pcs.IVC, ...
0350                   1e-2);
0351 mdl.mat_names = vertcat(mdl.mat_names, pcs.PA.mat_names,pcs.Aorta.mat_names, ...
0352                         pcs.SVC.mat_names, pcs.IVC.mat_names);
0353 
0354 eidors_msg('Merging trachea', 2);
0355 mdl = merge_meshes(mdl, organs.Trachea, 5e-2);
0356 mdl.mat_names{end+1} = 'Trachea';
0357 
0358 
0359 function organs = mk_individual_organs
0360 org_path = get_contrib_data('bp3d-torso','bp3d_organs.mat');
0361 srf = load(org_path);
0362 
0363 eidors_msg('Meshing vessels', 2);
0364 vessels = {'Aorta','IVC','SVC','PA', 'LSPV', 'LIPV','RSPV','RIPV'};
0365 for i = 1:numel(vessels)
0366     v = vessels{i};
0367     disp(v)
0368     [vessel, blood] = eidors_cache(@mk_vessel, srf.(v));
0369     organs.(v) = vessel;
0370     b = sprintf('Blood%s', v);
0371     organs.(b) = blood;
0372     print_mdl_info(vessel, v);
0373     print_mdl_info(blood, b);
0374 end
0375 
0376 eidors_msg('Meshing trachea and lungs', 2);
0377 organs.Trachea = gmsh_stl2tet(srf.Trachea,[],[],1);
0378 print_mdl_info(organs.Trachea, 'Trachea');
0379 organs.LungLeft = gmsh_stl2tet(srf.LungLeft,[],[],2);
0380 print_mdl_info(organs.LungLeft, 'LungLeft');
0381 organs.LungRight = gmsh_stl2tet(srf.LungRight,[],[],2);
0382 print_mdl_info(organs.LungRight, 'LungRight');
0383 
0384 eidors_msg('Meshing left heart', 2);
0385 left = eidors_cache(@mk_blood_left_surf,srf);
0386 organs.BloodLeftHeart = gmsh_stl2tet(left, [3, 20],[],4);
0387 print_mdl_info(organs.BloodLeftHeart, 'BloodLeftHeart');
0388 
0389 eidors_msg('Meshing right heart', 2);
0390 right = eidors_cache(@mk_blood_right_surf,srf);
0391 organs.BloodRightHeart = gmsh_stl2tet(right, [3, 20],[],4);
0392 print_mdl_info(organs.BloodRightHeart, 'BloodRightHeart');
0393 
0394 eidors_msg('Meshing heart wall', 2);
0395 heart = eidors_cache(@mk_heart_surf, srf);
0396 organs.HeartWall = gmsh_stl2tet(heart,[],4);
0397 print_mdl_info(organs.HeartWall, 'HeartWall');
0398 
0399 
0400 function right = mk_blood_right_surf(srf)
0401 fn = @(x) flip_normals(x);
0402 right = merge_meshes(srf.BloodRight, ...
0403                      fn(srf.SVC.cap), ...
0404                      fn(srf.IVC.cap), ...
0405                      fn(srf.PA.cap));
0406 
0407 
0408 function left = mk_blood_left_surf(srf)
0409 fn = @(x) flip_normals(x);
0410 left = merge_meshes(srf.BloodLeft, ...
0411                     fn(srf.Aorta.cap), ...
0412                     fn(srf.LSPV.cap), ...
0413                     fn(srf.LIPV.cap), ...
0414                     fn(srf.RSPV.cap), ...
0415                     fn(srf.RIPV.cap), ...
0416                     1e-2);
0417 
0418 
0419 function [heart, right, left] = mk_heart_surf(srf)
0420 right = merge_meshes(srf.BloodRight, ...
0421                      srf.SVC.ring, ...
0422                      srf.IVC.ring, ...
0423                      srf.PA.ring, ...
0424                      1e-2);
0425 right = flip_normals(right);
0426 
0427 left = merge_meshes(srf.BloodLeft, ...
0428                     srf.Aorta.ring, ...
0429                     srf.LSPV.ring, ...
0430                     srf.LIPV.ring, ...
0431                     srf.RSPV.ring, ...
0432                     srf.RIPV.ring, ...
0433                     1e-2);
0434 left = flip_normals(left);
0435 
0436 heart = merge_meshes(bd(left), bd(right), bd(srf.HeartShell));
0437 
0438 
0439 function [V, B] = mk_vessel(vstruct)
0440 msh = merge_meshes(vstruct.outer, vstruct.ring, ...
0441                     flip_normals(vstruct.inner),1e-2);
0442 V = gmsh_stl2tet(msh, 1,[],3); % max 1 mm edge (vessels are 2 mm thick
0443 msh = merge_meshes(vstruct.cap, vstruct.inner);
0444 B = gmsh_stl2tet(msh,[],1);
0445 
0446 
0447 function mdl = flip_normals(mdl)
0448 mdl.elems = mdl.elems(:,[1 3 2]);
0449 
0450 
0451 function mdl = bd(mdl)
0452 mdl.boundary = mdl.elems;
0453 
0454 
0455 function print_mdl_info(mdl, name)
0456 fprintf('=> %20s: %6d nodes, %7d elems\n',name, num_nodes(mdl), num_elems(mdl));
0457 
0458 
0459 function out = get_cache_path
0460 global eidors_objects
0461 out = [eidors_objects.model_cache filesep 'bp3d-torso'];
0462 if ~exist(out, 'dir')
0463    mkdir(out);
0464 end
mk_thorax_model_bp3d

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
