map_RM_slice

map_RM_slice - map a 3D reconstruction matrix to an arbitrary 2D slice

 IMDL2 = map_RM_slice(IMDL3, LEVEL)
 IMDL2 = map_RM_slice(IMDL3, LEVEL, OPT) 
  Maps the reconstruction matrix in IMDL3.solve_use_matrix.RM onto an 
  arbitrary slice defined by LEVEL. IMDL3 must be a dual-model inverse 
  model providing IMDL3.rec_model and IMDL3.fwd_model.coarse2fine.
  The OPT struct defines the pixel grid of the new reconstruction; see 
  MK_PIXEL_SLICE for details.

 Note: 
 - If the provided imdl.rec_model used medical orientation (decreasing
   imdl.rec_model.mdl_slice_mapper.y_pts), it may be necessary to flip the
   orientation of the resultant imdl.rec_model.
 - For selecting a horizontal cut, consider SELECT_RM_SLICE instead. 
 - At the moment, MAP_RM_SLICE only works for voxel-based rec_models as
   as provided by MK_GRID_MODEL or MK_VOXEL_VOLUME.

 See also MK_PIXEL_SLICE, SELECT_RM_SLICE, MK_VOXEL_VOLUME, MK_GRID_MODEL

0001 function imdl = map_RM_slice(imdl,varargin)
0002 %map_RM_slice - map a 3D reconstruction matrix to an arbitrary 2D slice
0003 %
0004 % IMDL2 = map_RM_slice(IMDL3, LEVEL)
0005 % IMDL2 = map_RM_slice(IMDL3, LEVEL, OPT)
0006 %  Maps the reconstruction matrix in IMDL3.solve_use_matrix.RM onto an
0007 %  arbitrary slice defined by LEVEL. IMDL3 must be a dual-model inverse
0008 %  model providing IMDL3.rec_model and IMDL3.fwd_model.coarse2fine.
0009 %  The OPT struct defines the pixel grid of the new reconstruction; see
0010 %  MK_PIXEL_SLICE for details.
0011 %
0012 % Note:
0013 % - If the provided imdl.rec_model used medical orientation (decreasing
0014 %   imdl.rec_model.mdl_slice_mapper.y_pts), it may be necessary to flip the
0015 %   orientation of the resultant imdl.rec_model.
0016 % - For selecting a horizontal cut, consider SELECT_RM_SLICE instead.
0017 % - At the moment, MAP_RM_SLICE only works for voxel-based rec_models as
0018 %   as provided by MK_GRID_MODEL or MK_VOXEL_VOLUME.
0019 %
0020 % See also MK_PIXEL_SLICE, SELECT_RM_SLICE, MK_VOXEL_VOLUME, MK_GRID_MODEL
0021 
0022 % TODO:
0023 % - Expand support for any 3D rec_model (see 3D-3D dual tutorial), test
0024 % with inv_solve_diff_GN_one_step (factor out get_RM).
0025 
0026 
0027 if ischar(imdl) && strcmp(imdl,'UNIT_TEST'); do_unit_test; return; end
0028 [rec, fwd] = mk_pixel_slice(imdl.rec_model,varargin{:});
0029 
0030 % f = fwd_model elems
0031 % r = rec_model elems
0032 % c = rec_model values
0033 % s = slice values
0034 
0035 c2f = imdl.fwd_model.coarse2fine; % fraction of fwd_model elem in coarse value
0036 c2r = imdl.rec_model.coarse2fine; % fraction of rec_model elem in coarse value
0037 s2r = fwd.coarse2fine;            % fraction of rec_model elem in slice value
0038 
0039 rvol = get_elem_volume(imdl.rec_model,'no_c2f');
0040 r2s = s2r' * spdiag(rvol); % volume of each slice elem built from rec elems
0041 
0042 % we need to divide by the volume of each pixel, but we don't know
0043 % the height, so we estimate from the area and the maximum volume
0044 s_area = get_elem_volume(rec); 
0045 s_hght = max(sum(r2s,2)) ./ s_area; 
0046 svol = s_area .* s_hght;
0047 r2s = r2s ./ svol;
0048 
0049 if 0
0050    c = ones(size(c2f,2),1);
0051    r = c2r * c;
0052    s = r2s * r;
0053    unit_test_cmp('still 1', s(30),1,1e-13)
0054 end
0055 
0056 s2c = r2s * c2r;
0057 RM = s2c * imdl.solve_use_matrix.RM;
0058 
0059 % to map from the slice to the original fwd_model, we need r2c.
0060 % Same approach:
0061 r2c = c2r' * spdiag(rvol);
0062 % again, in general, we don't know the volume of each element (because the
0063 % coarse elements on the rec_model could theoretically extend beyond the
0064 % elements of it, though we don't build such rec_models). We take the max
0065 % of two estimates:
0066 vol_fmdl = get_elem_volume(imdl.fwd_model);
0067 vol_rmdl = get_elem_volume(imdl.rec_model);
0068 rvol = max([vol_rmdl, vol_fmdl],[],2);
0069 r2c = r2c ./ rvol;
0070 
0071 c2f = imdl.fwd_model.coarse2fine * r2c * s2r;
0072 
0073 if 0
0074    fmdl = imdl.fwd_model;
0075    fmdl.coarse2fine = c2f;
0076    img = mk_image(fmdl, ones(size(c2f,2),1));
0077    show_fem(img)
0078 end
0079 
0080 imdl.fwd_model.coarse2fine = c2f;
0081 imdl.solve_use_matrix.RM = RM;
0082 imdl.rec_model = rec;
0083 
0084 function do_unit_test
0085    warning('off','EIDORS:FirstImageOnly');
0086    [vh,vi] = test_fwd_solutions;
0087    % inverse geometry
0088    fmdl= ng_mk_cyl_models([4,1,.5],[16,1.5,2.5],[0.05]);
0089    fmdl= mystim_square(fmdl);
0090 
0091    vopt.imgsz = [32 32];
0092    vopt.zvec = linspace( 0.75,3.25,6);
0093    vopt.save_memory = 1;
0094 %  opt.noise_figure = 2;
0095    weight = 8.42412109211;
0096    [imdl,opt.distr] = GREIT3D_distribution(fmdl, vopt);
0097    imdl = mk_GREIT_model(imdl, 0.2, weight, opt);
0098 
0099 
0100    img = inv_solve(imdl, vh, vi);
0101    unit_test_cmp('img size', size(img.elem_data), [4280,5]);
0102 
0103    img.show_slices.img_cols = 1;
0104    slice1 = calc_slices(img,[inf inf 2]);
0105 
0106    subplot(231); show_fem(fmdl); title 'fmdl'
0107    img0 = img;
0108    img0.elem_data = img0.elem_data(:,3);
0109    subplot(234); show_fem(img0);  title '3D'
0110    subplot(232); show_slices(img,[inf,inf,2]); title '3D rec';
0111    eidors_colourbar(img);
0112 
0113    imd2 = map_RM_slice(imdl,[Inf Inf 2.0],struct('z_depth',0.25));
0114    % flip for medical orientation
0115    imd2.rec_model.mdl_slice_mapper.y_pts = ...
0116       fliplr(imd2.rec_model.mdl_slice_mapper.y_pts);
0117 
0118    unit_test_cmp('RM size', size(imd2.solve_use_matrix.RM), [856,928]);
0119    rm1 = imdl.solve_use_matrix.RM(856*2 + (1:856),:);
0120    rm2 = imd2.solve_use_matrix.RM;
0121    unit_test_cmp('RM values', rm2, rm1, 1e-10);
0122 
0123    img = inv_solve(imd2, vh, vi);
0124    unit_test_cmp('img size', size(img.elem_data), [856,5]);
0125    im1 = img0.elem_data(856*2 + (1:856));
0126    unit_test_cmp ('img values', img.elem_data(:,3), im1, 1e-13);
0127 
0128    img.show_slices.img_cols = 1;
0129    slice2 = calc_slices(img);
0130    subplot(235), show_slices(img); title('2D rec')
0131    eidors_colourbar(img);
0132 
0133    subplot(2,3,[3,6])
0134 
0135    level = struct( 'centre', [0 0 2],...
0136                    'normal_angle',[1 -1 1 0]/sqrt(3)');
0137    opt.z_depth = 0.5;
0138    opt.square_pixes = true;
0139    opt.imgsz = [64 64];
0140    imd2 = map_RM_slice(imdl, level, opt);
0141    img = inv_solve(imd2, vh, vi);
0142    img.elem_data = img.elem_data(:,3);
0143    h1 = show_fem(img);
0144 %    h1.EdgeColor = 0.8*ones(3,1);
0145 
0146    
0147 
0148    xlim([-1 1]); ylim([-1 1]); zlim([0 4]);
0149    hold on
0150    h2 = show_fem(imd2.fwd_model);
0151    h2.EdgeColor = 0.8*ones(3,1);
0152    show_fem(img0.fwd_model);
0153    hold off
0154    view(-20,40)
0155 
0156 
0157 
0158    unit_test_cmp('slice Nan', isnan(slice1), isnan(slice2))
0159    slice1(isnan(slice1))= 0;
0160    slice2(isnan(slice2))= 0;
0161    % special case:
0162    unit_test_cmp('slice values', slice1, slice2, 1e-13)
0163 
0164   
0165 
0166    warning('on','EIDORS:FirstImageOnly');
0167 
0168 function fmdl = mystim_square(fmdl);
0169    [~,fmdl] = elec_rearrange([16,2],'square', fmdl);
0170    [fmdl.stimulation, fmdl.meas_select] =  ...
0171        mk_stim_patterns(32,1,[0,5],[0,5],{},1);
0172 
0173 function [vh,vi] = test_fwd_solutions;
0174    posns= linspace(1.0,3.0,5);
0175    str=''; for i=1:length(posns);
0176       extra{i} = sprintf('ball%03d',round(posns(i)*100));
0177       str = [str,sprintf('solid %s=sphere(0.5,0.3,%f;0.1); ', extra{i}, posns(i))];
0178    end
0179    extra{i+1} = str;
0180    fmdl= ng_mk_cyl_models([4,1,.2],[16,1.5,2.5],[0.05],extra); 
0181    fmdl = mystim_square(fmdl);
0182    
0183    img= mk_image(fmdl,1);
0184    vh = fwd_solve(img);
0185    %vh = add_noise(2000, vh);
0186    for i=1:length(posns);
0187       img= mk_image(fmdl,1);
0188       img.elem_data(fmdl.mat_idx{i+1}) = 2;
0189       vi{i} = fwd_solve(img);
0190    end;
0191    vi = [vi{:}];