choose_image_SNR

PURPOSE ^

% CHOOSE_IMAGE_SNR: choose hyperparameter based on image SNR calculation

SYNOPSIS ^

function HP = choose_image_SNR(imdl)

DESCRIPTION ^

% CHOOSE_IMAGE_SNR: choose hyperparameter based on image SNR calculation
 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).

   HP = choose_image_SNR(imdl)

 Output:
   HP   - hyperparameter selected to achieve the desired image SNR

 Input:
   imdl            - inverse model (EIDORS struct)
      imdl.hyperparameter.image_SNR - desired target image SNR

 NOTE
  This script was adapted from the function CHOOSE_NOISE_FIGURE


 See also: CALC_IMAGE_SNR, CHOOSE_NOISE_FIGURE

 Fabian Braun, December 2016

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function HP = choose_image_SNR(imdl)
0002 %% CHOOSE_IMAGE_SNR: choose hyperparameter based on image SNR calculation
0003 % 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 %   HP = choose_image_SNR(imdl)
0008 %
0009 % Output:
0010 %   HP   - hyperparameter selected to achieve the desired image SNR
0011 %
0012 % Input:
0013 %   imdl            - inverse model (EIDORS struct)
0014 %      imdl.hyperparameter.image_SNR - desired target image SNR
0015 %
0016 % NOTE
0017 %  This script was adapted from the function CHOOSE_NOISE_FIGURE
0018 %
0019 %
0020 % See also: CALC_IMAGE_SNR, CHOOSE_NOISE_FIGURE
0021 %
0022 % Fabian Braun, December 2016
0023 %
0024 
0025 % (C) 2016 Fabian Braun. License: GPL version 2 or version 3
0026 % $Id: choose_image_SNR.m 5424 2017-04-25 17:45:19Z aadler $
0027 
0028 
0029 if ~isfield(imdl.hyperparameter, 'SNR_func')
0030     imdl.hyperparameter.SNR_func = @calc_image_SNR;
0031 end
0032 
0033 SnrGoal = imdl.hyperparameter.image_SNR;
0034 assert(SnrGoal > 0, 'SNR is expected to be greater than zero!');
0035 
0036 hpInitValue = 1E-1;
0037 try % remove the value field
0038    hpInitValue = imdl.hyperparameter.value;
0039    imdl.hyperparameter = rmfield(imdl.hyperparameter,'value');
0040 end
0041 
0042 copt.boost_priority = 3;
0043 copt.fstr = 'choose_image_SNR';
0044 HP = eidors_cache(@HP_for_SNR_search,{SnrGoal,imdl,hpInitValue},copt);
0045 
0046 
0047 
0048 function HP = HP_for_SNR_search(SnrGoal,imdl,hpInitValue)
0049    llv = eidors_msg('log_level');
0050    if llv==3; eidors_msg('log_level',2); end % HIDE A LOT OF MESSAGES
0051 
0052    if true
0053        % new version first finding maximal possible SNR
0054        % first of all look for maximal possible SNR
0055        % we can't go higher, for once "the sky is the limit" does not apply :-)
0056        fms_opts.MaxIter = 15;
0057        fms_opts.TolFun = 0.01 * SnrGoal;   % we don't need to be more accurate
0058        [HpMaxSnr, NegSnrMax] = fminsearch(@(x) -calc_image_SNR(imdl, x), hpInitValue, fms_opts);
0059        % TODO: how to ensure we won't fall into a local extremum?
0060        
0061        assert(SnrGoal <= abs(NegSnrMax), ...
0062             sprintf('Desired SNR (%0.4d) exceeds maximal possible SNR of %0.4d!', SnrGoal, abs(NegSnrMax)));
0063    else
0064        % old version with very dummy initialization
0065        HpMaxSnr = hpInitValue;
0066    end
0067    
0068    % now go and start from this to find best one
0069    [hp, x, y] = search1(SnrGoal, imdl, -log10(HpMaxSnr));
0070    eidors_msg('Crossed target value, now refining search in this region...', 2);
0071 
0072    dx= hp-linspace(-0.7,0.7,5);
0073    [hp, Snr] = search2(SnrGoal, imdl, dx);
0074 
0075    dx= hp-linspace(-0.2,0.2,5);
0076    [hp, Snr] = search2(SnrGoal,imdl, dx);
0077 
0078    dx= hp-linspace(-0.1,0.1,5); 
0079    [hp, Snr] = search2(SnrGoal,imdl, dx);
0080 
0081    dx= hp-linspace(-0.05,0.05,21);
0082    [hp, Snr] = search2(SnrGoal,imdl, dx);
0083    
0084    % TODO: verify that SNR = f(lambda) has positive derivative
0085    % else we're in the wrong range (after SNR max)
0086 
0087    HP= 10^-hp;
0088    assert(HP < HpMaxSnr, 'Something went wrong we''re on the wrong side of the curve!');
0089    eidors_msg('log_level',llv); %% Reset log_level
0090    
0091    % TODO: check that we achieved our goal without too much deviation
0092   
0093 function [hp, x, y] = search1(SnrGoal, imdl, hp)
0094   [Snr,SE]=imdl.hyperparameter.SNR_func(imdl, 10^(-hp));
0095   hpInitial = hp;
0096   SnrInitial = Snr;
0097   y = Snr;
0098   x = 10^(-hp);
0099    if     Snr > SnrGoal; 
0100        dir = 1;   % we're too high >> decrease hp
0101    elseif Snr < SnrGoal; 
0102        dir = -1;  % we're too low >> increase hp
0103        warning('this should not happen!!!');
0104    else   dir = 0; end
0105    
0106 %    slope = dir;
0107    while  (dir*Snr+3*SE > dir*SnrGoal)% || (sign(slope) ~= dir)
0108      hp= hp+0.5*dir;
0109      [Snr,SE]=imdl.hyperparameter.SNR_func( imdl, 10^(-hp));
0110      y = [y Snr];
0111      x = [x 10^(-hp)];
0112 %      slope = (Snr - SnrInitial)/(10^(-hp) - 10^(-hpInitial));
0113      
0114      if Snr < 0     % TODO: verify that this approach works in all situations
0115        % NO GOOD: we're getting in the negative t-Value zone
0116        % --> fall back to initial hp and inverse direction
0117        hp = hpInitial;
0118        dir = -dir;
0119        Snr = inf;
0120        eidors_msg('Falling back to initial hp and reversing search direction due to negative SNRs.', 2);
0121      end
0122    end
0123 
0124 function [hp, Snr] =search2(SnrGoal, imdl, dx)
0125    hp = [];
0126    it = 0;
0127    Snr = zeros(1,length(dx));
0128    while isempty(hp) && it < 2  % TODO: one could increase the number of iterations and reduce Nnoise in calc_image_SNR
0129        it = it+1;
0130        %if it > 1 keyboard, end
0131        for k=1:length(dx)
0132            Snr(k)=Snr(k)+imdl.hyperparameter.SNR_func( imdl, 10^-dx(k));
0133        end
0134        if 0 % TODO: validate which behaviour is more accurate!
0135          log_snr = log10(Snr/it);
0136          p= polyfit( dx, log_snr-log10(SnrGoal), 3);
0137        else
0138          p= polyfit( dx, (Snr/it)-SnrGoal, 3);
0139        end
0140        hp = roots(p);
0141        %eliminate complex roots
0142        hp = hp(imag(hp)==0);
0143        %eliminate out of range roots
0144        hp = hp( hp<max(dx) & hp>min(dx) );  %USE if poly>1
0145        % pick the root with the smallest derivative
0146        if numel(hp) >1
0147            p2 = p.*(numel(p)-1:-1:0); p2(end) = [];
0148            [jnk, pos] = min(abs(polyval(p2,hp)));
0149            hp = hp(pos);
0150        end
0151        
0152        if 0  % visualize for debug purposes
0153         figure(); subplot(121);
0154         plot(dx, Snr-SnrGoal, ':.')
0155         hold on;
0156         xs = linspace(min(dx), max(dx), 100);
0157         plot(xs, polyval(p, xs), 'r')
0158         
0159         subplot(122);
0160         plot(dx, log10(Snr/it)-log10(SnrGoal))
0161        end
0162    end
0163    if isempty(hp)
0164        %fallback
0165        [jnk,idx] = min(abs((Snr/it)-SnrGoal));
0166 %        [jnk,idx] = min(abs(log_snr-log10(SnrGoal)));
0167        hp = dx(idx);
0168    end
Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005