define_ROIs

 DEFINE_ROIS: define ROI regions in image

 [ROIs, imgR, imdl] = define_ROIs(imdl, xcuts, ycuts, ..., opts)
   imdl: inv_ or fwd_model structure
   xcuts, ycuts: x and y cuts of ROI grid
   opts: options structure
     opts.reorder = reorder defined ROIs
     opts.normalize = set ROI of each =1
     opts.usec2f    = true if use fwd_model.coarse2fine
 
   ROIs: matrix [n_ROIs x n_elems]
   imgR: image of ROIs (mostly for debugging)
   imdl: inv_model which reconstructs onto ROIs (if possible)

 Example:
  fmdl = mk_library_model('neonate_16el_lungs');
  fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
  imdl= select_imdl(fmdl,'GREIT:NF=0.500 32x32');
  opts.reorder = reshape(1:16,4,4)';
  [ROIs,imgR] = define_ROIs(imdl,-1:0.5:1,1:-0.5:-1,opts);

0001 function [ROIs, imgR, imdl] = define_ROIs(imdl, varargin)
0002 % DEFINE_ROIS: define ROI regions in image
0003 %
0004 % [ROIs, imgR, imdl] = define_ROIs(imdl, xcuts, ycuts, ..., opts)
0005 %   imdl: inv_ or fwd_model structure
0006 %   xcuts, ycuts: x and y cuts of ROI grid
0007 %   opts: options structure
0008 %     opts.reorder = reorder defined ROIs
0009 %     opts.normalize = set ROI of each =1
0010 %     opts.usec2f    = true if use fwd_model.coarse2fine
0011 %
0012 %   ROIs: matrix [n_ROIs x n_elems]
0013 %   imgR: image of ROIs (mostly for debugging)
0014 %   imdl: inv_model which reconstructs onto ROIs (if possible)
0015 %
0016 % Example:
0017 %  fmdl = mk_library_model('neonate_16el_lungs');
0018 %  fmdl.stimulation = mk_stim_patterns(16,1,[0,1],[0,1],{},1);
0019 %  imdl= select_imdl(fmdl,'GREIT:NF=0.500 32x32');
0020 %  opts.reorder = reshape(1:16,4,4)';
0021 %  [ROIs,imgR] = define_ROIs(imdl,-1:0.5:1,1:-0.5:-1,opts);
0022 
0023 % TODO:
0024 %   - create imdl
0025 %   - add 3D
0026 %   - create rec_model if not there
0027 %   - can it work for non-solve_use_matrix?
0028 %
0029 
0030 if ischar(imdl) && strcmp(imdl,'UNIT_TEST'); do_unit_test; return; end
0031 
0032 if isstruct(varargin{end})
0033    opts = varargin{end};
0034    last = 1;
0035 else
0036    opts = struct([]);
0037    last = 0;
0038 end
0039 % Set defaults
0040 %     opts.reorder = reorder defined ROIs
0041 %     opts.normalize = set ROI of each =1
0042 %     opts.usec2f    = true if use fwd_model.coarse2fine
0043 opts = mergestructs(struct( ...
0044     'reorder',[],'normalize',false,'usec2f',false), opts);
0045 
0046 recc2f = 1;
0047 switch imdl.type
0048    case 'fwd_model';
0049       fmdl = imdl;
0050    case 'inv_model';
0051       if isfield(imdl,'rec_model');
0052          fmdl = imdl.rec_model;
0053          recc2f = fmdl.coarse2fine;
0054       else
0055          fmdl = imdl.fwd_model;
0056       end
0057 
0058    otherwise;
0059       error('define_ROIs: require fwd_ or inv_model');
0060 end
0061 if opts.usec2f
0062   recc2f = fmdl.coarse2fine;
0063 end
0064 
0065 
0066 rmdl = mk_grid_model([],varargin{1:end-last});
0067 
0068 ROIa = mk_approx_c2f(fmdl,rmdl);
0069 % Octave ./ doesn't expand dimensions for sparse for v9.2
0070 % ROIs = ROIa.' * recc2f./sum(recc2f,1);
0071 % Bug report submitted to octave: 2024-12-19
0072 ROIs = ROIa.' * bsxfun(@rdivide, recc2f, sum(recc2f,1) );
0073 if ~isempty(opts.reorder)
0074    ROIs = ROIs(opts.reorder(:),:);
0075 end
0076 if opts.normalize
0077    ROIs = ROIs ./ sum(ROIs,2);
0078 end
0079 imgR= mk_image(fmdl,ROIs.');
0080 imgR.calc_colours.ref_level = 0;
0081 
0082 % Can only replace in this case
0083 if strcmp(imdl.type,'inv_model') && ...
0084    isfield(imdl,'rec_model') && ...
0085    isfield(imdl,'solve_use_matrix')
0086 
0087    imdl.rec_model.coarse2fine = ROIa;
0088    imdl.solve_use_matrix.RM = ROIs*imdl.solve_use_matrix.RM;
0089 end
0090 
0091 function do_unit_test
0092    do_GREIT_test
0093 %  do_simple_test
0094 %  do_3D_test
0095 
0096 function do_GREIT_test
0097    imdl = mk_common_model('n3r2',[16,2]);
0098    [ROIs,imgR] = define_ROIs(imdl,-1:0.5:1,-1:0.5:1);
0099    unit_test_cmp('n3r2: 2D 1',     size(ROIs), [16,828]);
0100    unit_test_cmp('n3r2: 2D 2',     ROIs(:,1:9), ...
0101        [0;0;0;0;0;0;0;0;0;0;0;1;0;0;0;0]*ones(1,9));
0102    unit_test_cmp('n3r2: 2D 2',     ROIs(:,10:11), ...
0103        [0,0;0,0;0,0;0,0;0,0;0,0;0,0;0,0;
0104         0,0;0,0;0,0;95,80;0,0;0,0;0,0;5,20]/100,1e-14);
0105    [ROIs,imgR] = define_ROIs(imdl,-1:1,-1:1,0:2);
0106    unit_test_cmp('n3r2: 3D 1',     size(ROIs), [8,828]);
0107    unit_test_cmp('n3r2: 3D 2',     ROIs(:,1:9), ...
0108        [0;0;0;1;0;0;0;0]*ones(1,9),1e-14);
0109    unit_test_cmp('n3r2: 3D 3',     ROIs(:,61:69), ...
0110        [0;0;1;0;0;0;0;0]*ones(1,9),1e-14);
0111 
0112    fmdl = mk_library_model('neonate_32el_lungs');
0113    [stim,msel] = mk_stim_patterns(32,1,[0,5],[0,5],{'no_meas_current_next2'},1);
0114    fmdl.stimulation = stim;
0115    fmdl.meas_select = msel;
0116    imdl= select_imdl(fmdl,'GREIT:NF=0.500 32x32');
0117    [ROIs,imgR] = define_ROIs(imdl,-1:0.5:1,-1:0.5:1);
0118    subplot(221); show_slices(imgR)
0119    imdl= select_imdl(fmdl,'GREIT:NF=0.500 16x16');
0120    [ROIs,imgR] = define_ROIs(imdl,-1:0.5:1,-1:0.5:1);
0121    subplot(222); show_slices(imgR)
0122    opts.reorder = reshape(1:16,4,4)';
0123    [ROIs,imgR] = define_ROIs(imdl,-1:0.5:1,-1:0.5:1,opts);
0124    subplot(223); show_slices(imgR)
0125    [ROIs,imgR] = define_ROIs(imdl,-1:0.5:1,1:-0.5:-1,opts);
0126    subplot(224); show_slices(imgR)