calc_image_SNR

% CALC_IMAGE_SNR: Calculates the signal-to-noise ratio (SNR) in the image 
 domain as proposed by Braun et al. in:
 F Braun et al., A Versatile Noise Performance Metric for Electrical
 Impedance Tomography Algorithms, IEEE Trans. Biomed. Eng. 2017 (submitted).

   [SNRmean, SE, debug] = calc_image_SNR(imdl, hyperparameter, doPlot)

 Output:
   SNRmean   - mean of all SNR values as in equ. (9) in publication below
   SE        - std of all SNR values as in equ. (9) in publication below
   debug     - structure holding some information used for debug purposes

 Input:
   imdl            - inverse model (EIDORS struct)
      imdl.hyperparameter.roi_scaling_factor - amount of model shrinking
                                   to determine ROI where seed n_T targets
                                   (DEFAULT: 0.5);
      imdl.hyperparameter.n_targets - number of targets to seed in ROI n_T
                                   (DEFAULT: n_T = 200);
      imdl.hyperparameter.target_radius - relative target radius r_T
                                   (DEFAULT: r_T = 0.05);
      imdl.hyperparameter.xyzr_targets - vector 4 x num targets n_T
                                   specify targets manually [x y z r]
   hyperparameter  - desired hyperparameter value, this will overwrite the
                     imdl.hyperparameter.value and is for compatibility
                     purposes with the function calc_noise_figure()
   doPlot    will enable plotting if set to true (default = false)


 See also: CALC_NOISE_FIGURE

 Fabian Braun, December 2016

 CITATION_REQUEST:
 AUTHOR: F Braun et al.
 TITLE: A Versatile Noise Performance Metric for Electrical Impedance Tomography Algorithms
 JOURNAL: IEEE Transactions on Biomedical Engineering
 YEAR: 2017
 VOL: 64
 NUM: 10
 PAGE: 2321-2330
 DOI: 10.1109/TBME.2017.2659540

0001 function [SNRmean, SE, debug] = calc_image_SNR(imdl, hyperparameter, doPlot)
0002 %% CALC_IMAGE_SNR: Calculates the signal-to-noise ratio (SNR) in the image
0003 % domain as proposed by Braun et al. in:
0004 % F Braun et al., A Versatile Noise Performance Metric for Electrical
0005 % Impedance Tomography Algorithms, IEEE Trans. Biomed. Eng. 2017 (submitted).
0006 %
0007 %   [SNRmean, SE, debug] = calc_image_SNR(imdl, hyperparameter, doPlot)
0008 %
0009 % Output:
0010 %   SNRmean   - mean of all SNR values as in equ. (9) in publication below
0011 %   SE        - std of all SNR values as in equ. (9) in publication below
0012 %   debug     - structure holding some information used for debug purposes
0013 %
0014 % Input:
0015 %   imdl            - inverse model (EIDORS struct)
0016 %      imdl.hyperparameter.roi_scaling_factor - amount of model shrinking
0017 %                                   to determine ROI where seed n_T targets
0018 %                                   (DEFAULT: 0.5);
0019 %      imdl.hyperparameter.n_targets - number of targets to seed in ROI n_T
0020 %                                   (DEFAULT: n_T = 200);
0021 %      imdl.hyperparameter.target_radius - relative target radius r_T
0022 %                                   (DEFAULT: r_T = 0.05);
0023 %      imdl.hyperparameter.xyzr_targets - vector 4 x num targets n_T
0024 %                                   specify targets manually [x y z r]
0025 %   hyperparameter  - desired hyperparameter value, this will overwrite the
0026 %                     imdl.hyperparameter.value and is for compatibility
0027 %                     purposes with the function calc_noise_figure()
0028 %   doPlot    will enable plotting if set to true (default = false)
0029 %
0030 %
0031 % See also: CALC_NOISE_FIGURE
0032 %
0033 % Fabian Braun, December 2016
0034 %
0035 % CITATION_REQUEST:
0036 % AUTHOR: F Braun et al.
0037 % TITLE: A Versatile Noise Performance Metric for Electrical Impedance Tomography Algorithms
0038 % JOURNAL: IEEE Transactions on Biomedical Engineering
0039 % YEAR: 2017
0040 % VOL: 64
0041 % NUM: 10
0042 % PAGE: 2321-2330
0043 % DOI: 10.1109/TBME.2017.2659540
0044 %
0045 
0046 % (C) 2016 Fabian Braun. License: GPL version 2 or version 3
0047 % $Id: calc_image_SNR.m 5646 2017-09-22 06:16:33Z fab-b $
0048 
0049 citeme(mfilename);
0050 
0051 %% Configuration
0052 % model shrunken by this factor delimits the ROI
0053 ROI_SCALING_FACTOR = 0.5; 
0054 % approximate number of targets to uniformly distribute in the ROI
0055 N_DESIRED_TARGETS = 200;
0056 % relative radius of the circular(2D)/spherical(3D) target (this is relative to the outer model radius)
0057 NORM_TGT_RADIUS = 0.05; % original
0058 % TEZ region threshold
0059 TEZ_REGION_THRESH = 0.25;  % as ratio of maximum value : QUATER AMPLITUDE THRESHOLD
0060 
0061 
0062 %% Parse and default inputs
0063 if ~exist('doPlot', 'var') || isempty(doPlot)
0064     doPlot = false;
0065 end
0066 if isfield(imdl.hyperparameter, 'roi_scaling_factor') && ~isempty(imdl.hyperparameter.roi_scaling_factor)
0067     ROI_SCALING_FACTOR = imdl.hyperparameter.roi_scaling_factor;
0068 end
0069 assert(ROI_SCALING_FACTOR < 1, 'ROI must be scaled smaller than effective model');
0070 if isfield(imdl.hyperparameter, 'n_targets') && ~isempty(imdl.hyperparameter.n_targets)
0071     N_DESIRED_TARGETS = imdl.hyperparameter.n_targets;
0072 end
0073 if isfield(imdl.hyperparameter, 'target_radius') && ~isempty(imdl.hyperparameter.target_radius)
0074     NORM_TGT_RADIUS = imdl.hyperparameter.target_radius;
0075 end
0076 if isfield(imdl.hyperparameter, 'xyzr_targets') && ~isempty(imdl.hyperparameter.xyzr_targets)
0077     xyzr_targets = imdl.hyperparameter.xyzr_targets;
0078 else
0079     xyzr_targets = [];
0080 end
0081 
0082 
0083 %% input parsing
0084 if nargin>=2 && numel(hyperparameter) == 1 && ~isempty(hyperparameter)
0085     imdl.hyperparameter.value = hyperparameter;
0086     % Remove function parameter because it will recurse
0087     try; imdl.hyperparameter = rmfield(imdl.hyperparameter,'func'); end
0088 end
0089 
0090 
0091 %% generate targets inside of rec_model
0092 if isfield(imdl, 'rec_model')
0093     RecMdl = imdl.rec_model;   
0094 else
0095     RecMdl = imdl.fwd_model;
0096 end
0097 MdlCtr = mean(RecMdl.nodes,1);
0098 
0099 if isempty(xyzr_targets)
0100     % targets have not been defined, we'll generate them automatically
0101     % by shrinking the model by the desired scaling factor and see the
0102     % descired number of targets inside the shrunk area
0103     %
0104 
0105     % first determine electrode level
0106     if isfield(RecMdl, 'mdl_slice_mapper') && isfield(RecMdl.mdl_slice_mapper, 'level')
0107         ElectrodeLevel = RecMdl.mdl_slice_mapper.level;    
0108     elseif isfield(RecMdl, 'mdl_slice_mapper') && isfield(RecMdl.mdl_slice_mapper, 'z_pts')
0109         ElectrodeLevel = RecMdl.mdl_slice_mapper.z_pts;    
0110         error('3D not supported yet, please seed targets manually via xyzr_targets');
0111     else
0112         if isfield(RecMdl.electrode(1), 'pos')
0113             elec_loc = cell2mat(cellfun(@(x) mean(x)', {RecMdl.electrode.pos}, 'uniformoutput', false))';
0114         elseif isfield(RecMdl.electrode(1), 'nodes')
0115             for i=1:length(RecMdl.electrode)
0116                enodesi = RecMdl.electrode(i).nodes;
0117                elec_loc(i,:) = mean( RecMdl.nodes( enodesi,:),1 );
0118             end  
0119         else
0120             error('not supported!');
0121         end
0122         ElectrodeLevel = mean(elec_loc,1);
0123     end
0124     
0125     % uniformly distribute N_DESIRED_TARGETS targets in the shrunken model
0126     % at the level of the electrodes
0127     
0128     % first, shrink original model
0129     RecMdlShrunk = RecMdl;
0130     RecMdlShrunk.nodes = RecMdlShrunk.nodes - repmat(MdlCtr, size(RecMdlShrunk.nodes, 1), 1);
0131     RecMdlShrunk.nodes = RecMdlShrunk.nodes * ROI_SCALING_FACTOR;
0132     RecMdlShrunk.nodes = RecMdlShrunk.nodes + repmat(MdlCtr, size(RecMdlShrunk.nodes, 1), 1);
0133 
0134     % find boundary of shrunken model
0135     BoundaryNodes = find_boundary(RecMdlShrunk);
0136     % TODO: there's still a little problem here that we don't always
0137     % get a nicely oriented boundary, investigate this!
0138     Boundary = order_loop(RecMdlShrunk.nodes(BoundaryNodes, :));
0139     Boundary = [Boundary; Boundary(1,:)];
0140     % pp = fourier_fit(Boundary, 10); Boundary = fourier_fit(pp, linspace(0,1,20));
0141 
0142     % seed uniformly in rectangle and remove outliers
0143     AreaPoly = polyarea(Boundary(:,1), Boundary(:,2));
0144     Bounds = [max(Boundary); min(Boundary)];
0145     AreaRect = prod([Bounds(1,:) - Bounds(2,:)]);
0146     nUniformTgts = AreaRect * (N_DESIRED_TARGETS / AreaPoly);    
0147 
0148     % Size of the ROI in the two dimensions (x/y)...
0149     RoiSize = abs(diff(Bounds,[], 1));
0150     % ensure uniform spacing and scale according to differences in x/y size
0151     ScaleX = RoiSize(1) / RoiSize(2);
0152     ScaleY = 1./ScaleX;                 
0153     % create ROI centers
0154     [xx, yy] = meshgrid(linspace(Bounds(2,1), Bounds(1,1), ceil(ScaleX * sqrt(nUniformTgts))), ...
0155                         linspace(Bounds(2,2), Bounds(1,2), ceil(ScaleY * sqrt(nUniformTgts))));
0156     IsInside = inpolygon(xx(:), yy(:), Boundary(:,1), Boundary(:,2));
0157     nTgts = sum(IsInside);
0158     assert(abs((nTgts - N_DESIRED_TARGETS)/N_DESIRED_TARGETS) < 0.15, 'Cannot make desired number of targets');
0159     xx = xx(IsInside);
0160     yy = yy(IsInside);    
0161     zz = ones(size(xx))*ElectrodeLevel(3);
0162 else
0163     % targets specified from outside the function, assign them properly
0164     xx = xyzr_targets(1,:)';
0165     yy = xyzr_targets(2,:)';
0166     zz = xyzr_targets(3,:)';
0167 end
0168 
0169 % set target size relative to maximal model radius
0170 BoundsFull = [max(RecMdl.nodes); min(RecMdl.nodes)];
0171 Rmodel = (max(BoundsFull(1,:) - BoundsFull(2,:)))/2;    % maximal model radius
0172 Rtarget = Rmodel * NORM_TGT_RADIUS;
0173 rr = ones(size(xx))*Rtarget;
0174 
0175 if ~isempty(xyzr_targets)
0176     try
0177         rr = xyzr_targets(4,:)';    % overwrite target radii if existing
0178     end
0179 end
0180     
0181 %% generate differential voltages for each conductivity target
0182 img = mk_image(imdl.fwd_model, 1);
0183 if elem_dim(imdl.fwd_model) == 3
0184     xyzr = [xx yy zz rr]';  
0185 elseif elem_dim(imdl.fwd_model) == 2
0186     xyzr = [xx yy rr]';    
0187 else
0188     error('unsupported dimensions');
0189 end
0190 [vh, vi, xyzrOut, c2f] = simulate_movement(img, xyzr);
0191 NotAssigned = ~ismember(xyzr', xyzrOut', 'rows');
0192 assert(sum(NotAssigned) == 0, 'Error: target(s) got missing...');
0193 vd = vi - repmat(vh, 1, size(vi,2));
0194 
0195 
0196 %% get reconstruction matrix
0197 RM = get_RM(imdl);
0198 % calculate volume/area of each element in RecMdl
0199 RecMdlVols = get_elem_volume(RecMdl);
0200 
0201 
0202 %% generate individual target evaluation zones (TEZs)
0203 % first evaluate image response and take quater amplitude pixels as TEZ
0204 
0205 % we calculate it the direct way as we can reuse it again further down!
0206 imgrs = RM*vd;
0207 imgrsNorm = imgrs;
0208 imgrsNorm = imgrsNorm ./ repmat(max(imgrsNorm, [], 1), size(imgrsNorm,1), 1);
0209 TEZs = double((imgrsNorm > TEZ_REGION_THRESH));
0210 
0211 clear imgrsNorm;
0212 
0213 ElemVols = spdiags(RecMdlVols, 0, length(RecMdlVols), length(RecMdlVols));
0214   
0215 % TEZ matrices with proper normalization (as defined in publication)
0216 %      z_i  =  a_i     *  c_i    (in publication)
0217 Z = ElemVols * TEZs;
0218 Vtez = sum(Z,1);                % volume/area of each TEZ
0219 Vt = 4/3 * pi * rr.^3;          % volume of each target
0220 % note: even for the 2D case we take the volume as a target double
0221 % as thick than another but with same area has double the influence
0222 K = Vtez ./ Vt';                % TEZ volume / target volume
0223 % \tilde{z}_i = z_i * k          (in publication)
0224 Ztilde = Z .* repmat(K, size(Z,1), 1);
0225 
0226 % TEZ matrices kept for debugging / bw-compatility reasons
0227 TEZsNonNorm = Z; 
0228 % normalize ROIs
0229 TEZs = TEZs ./ repmat(sum(TEZs,1), size(TEZs,1), 1);  
0230 
0231 
0232 %% calculate target-wise distinguishability
0233 % determine noise covariance matrix
0234 if isfield(imdl.hyperparameter, 'SigmaN') && ~isempty(imdl.hyperparameter.SigmaN)
0235     SigmaN = imdl.hyperparameter.SigmaN;    
0236 else
0237     SigmaN = speye(size(RM,2));        % assuming uniform and uncorrelated noise
0238 end
0239 
0240 Signal = diag(Ztilde' * imgrs);    % numerator in equation (9)
0241 VarPerPixel = diag(RM*SigmaN*RM');      
0242 Noise = sqrt(Z' * VarPerPixel);      % denominator in equation (9)
0243 
0244 % get the average of all SNRs
0245 SNRs = Signal ./ Noise;                 
0246 SNRmean = mean(SNRs); 
0247 SE = std(SNRs);
0248 
0249 if isnan(SNRmean)
0250     SNRmean = -inf;
0251 elseif isinf(SNRmean)
0252     keyboard;
0253 end
0254 
0255 try
0256   eidors_msg('SNR = %e (hp=%e)', SNRmean, imdl.hyperparameter.value, 1);
0257 end
0258 
0259 
0260 %% debug info: fraction of amplitude response (AR) inside each TEZ
0261 if nargout >= 3
0262     % target-wise amplitude response inside TEZ: in terms of energy
0263     ArInFrac = diag((TEZsNonNorm'*abs(RM*vd).^2) ./ (repmat(RecMdlVols, 1, size(TEZs,2))'*abs(RM*vd).^2));    
0264     debug.ArInFrac = ArInFrac;
0265     
0266     debug.SNRs = SNRs;
0267     debug.Signal = Signal;
0268     debug.Noise = Noise;
0269     debug.TEZs = TEZs;
0270     debug.TEZsNoNorm = TEZsNonNorm;
0271     
0272     debug.RoiBounds = Bounds;
0273     debug.RoiBoundary = Boundary;
0274     debug.MdlCtr = MdlCtr;
0275     debug.BoundsFull = BoundsFull;
0276     
0277     imgrsNorm = imgrs;
0278     imgrsNorm = imgrsNorm ./ repmat(max(imgrsNorm, [], 1), size(imgrsNorm,1), 1);
0279     debug.meanInNormImg = diag(TEZs'*imgrsNorm);
0280 end
0281 
0282 
0283 %% visualize for debug purposes
0284 if doPlot
0285     fig = figure;
0286     set(fig, 'position', [  182         700        1531         313]);
0287     
0288     if isfield(imdl.fwd_model, 'coarse2fine')
0289         MapTgts2Img = imdl.fwd_model.coarse2fine'*c2f;
0290     else
0291         try
0292             imdl.fwd_model.coarse2fine = mk_coarse_fine_mapping(imdl.fwd_model, imdl.rec_model);
0293             MapTgts2Img = imdl.fwd_model.coarse2fine'*c2f;
0294         catch
0295             warning('Unable to make c2f mapping');
0296             MapTgts2Img = c2f;
0297         end
0298     end
0299     MapTgts2Img = MapTgts2Img ./ repmat(sum(MapTgts2Img,2), 1, size(MapTgts2Img, 2));
0300     MapTgts2Img(isnan(MapTgts2Img(:))) = 0;
0301     
0302     img = mk_image(RecMdl, nan);
0303     img.elem_data = MapTgts2Img * SNRs;
0304     SnrImg = calc_slices(img);
0305     
0306     img.elem_data = nan(size(img.elem_data));
0307     img.elem_data = MapTgts2Img * Signal;
0308     AmpSens = img.elem_data;
0309     AmpSensImg = calc_slices(img);
0310     
0311     img.elem_data = nan(size(img.elem_data));
0312     img.elem_data = MapTgts2Img * Noise;
0313     NoiseSens = img.elem_data;
0314     NoiseSensImg = calc_slices(img);
0315     
0316     sp1 = subplot(131);
0317     imagescnan(SnrImg);
0318     title(['SNR image: ', num2str(SNRmean, '%0.2d')]);
0319     colorbar; colormap jet;
0320 
0321     % plot signal sensitivity image
0322     sp2 = subplot(132);
0323     imagescnan(AmpSensImg);
0324     title(['Amplitude response: ', num2str(nanmean(AmpSens(:)), '%0.2d')]);
0325     colorbar;
0326 
0327     % plot noise sensitivity image
0328     sp3 = subplot(133);
0329     imagescnan(NoiseSensImg);
0330     title(['Noise sensitivity: ', num2str(nanmean(NoiseSens(:)), '%0.2d')]);
0331     colorbar;
0332 
0333     if doPlot < 2
0334         linkaxes([sp1 sp2 sp3], 'xy');
0335     end
0336    
0337     if doPlot == 2
0338        % some extra debugging
0339        sp1 = subplot(131);
0340        imgTmp = img;
0341        iTgt = round(size(vd,2)/2);
0342        imgTmp.elem_data = RM*vd(:,iTgt);
0343        hh = show_fem(imgTmp);
0344        set(hh, 'edgecolor', 'none');
0345        hold on; 
0346        circle(xyzr(1:2, iTgt), xyzr(3, iTgt), 100, 'k');
0347        circle(xyzr(1:2, iTgt), Rtarget, 100, 'k');
0348        axis equal;
0349     end
0350     
0351     
0352     if doPlot == 3
0353        sp3 = subplot(133); 
0354        cla;
0355        hh = show_fem(RecMdl);
0356        set(hh, 'edgecolor', 'none')
0357        hold on;
0358        plot(Boundary(:,1), Boundary(:,2), '.-k');
0359        plot(xyzr(1,:), xyzr(2,:), '.r');
0360        plot(xyzr(1,:), xyzr(2,:), 'or');
0361     end
0362     
0363     if doPlot == 4
0364         % show pixel-wise noise sensitivity
0365        sp2 = subplot(132);
0366        cla;
0367        imgN = mk_image(imdl,0);
0368        imgN.elem_data = VarPerPixel;
0369        imgN.elem_data(sum(MapTgts2Img,2) == 0) = 0;
0370        imagescnan(calc_slices(imgN));
0371        title(['Pixel-wise NoiseSens: ', num2str(nanmean(imgN.elem_data(:)), '%0.2d')]);
0372        colorbar;        
0373     end
0374     
0375     if doPlot == 5
0376         sp1 = subplot(131);
0377         hist(SNRs); 
0378         xlim([0 max(xlim())]);
0379     end
0380     
0381 end
0382 
0383 end
0384 
0385