line_optimize

LINE_OPTIMIZE Cheap line optimizer
 [img fmin res] = line_optimize(imgk, dx, data0, opt)
 img     : output image
 fmin    : optimal step size
 res     : value of the objective function
 imgk    : starting image
 dx      : step direction
 data0   : data to fit
 opt     : options structure

 Options:
   opt.perturb  : vector of step sizes to try (default: [0.1 0.5 1.0])
   opt.min_value: lower limit of img values (default: -Inf)
   opt.max_value: upper limit of img values (default: Inf)
   opt.objective_function: @my_objective_fun
       handle to an objective funtion with the following signature
           val = my_objective_fun(data0, data, img0, img, opt)
       where:
           data0   : data to fit
           data    : data simulated on the current image
           img0    : starting image
           img     : current image
           opt     : the entire option structure (can be used to pass
           additional parameters)
       The default objective function is:
           val = norm(calc_difference_data(data0,data, img0.fwd_model));

 Note that the value of fmin will be limited by the maximum and minimum
 given in the opt.perturb vector.

 See also: INV_SOLVE_ABS_GN

0001 function [img fmin res] = line_optimize(imgk, dx, data0, opt)
0002 %LINE_OPTIMIZE Cheap line optimizer
0003 % [img fmin res] = line_optimize(imgk, dx, data0, opt)
0004 % img     : output image
0005 % fmin    : optimal step size
0006 % res     : value of the objective function
0007 % imgk    : starting image
0008 % dx      : step direction
0009 % data0   : data to fit
0010 % opt     : options structure
0011 %
0012 % Options:
0013 %   opt.perturb  : vector of step sizes to try (default: [0.1 0.5 1.0])
0014 %   opt.min_value: lower limit of img values (default: -Inf)
0015 %   opt.max_value: upper limit of img values (default: Inf)
0016 %   opt.objective_function: @my_objective_fun
0017 %       handle to an objective funtion with the following signature
0018 %           val = my_objective_fun(data0, data, img0, img, opt)
0019 %       where:
0020 %           data0   : data to fit
0021 %           data    : data simulated on the current image
0022 %           img0    : starting image
0023 %           img     : current image
0024 %           opt     : the entire option structure (can be used to pass
0025 %           additional parameters)
0026 %       The default objective function is:
0027 %           val = norm(calc_difference_data(data0,data, img0.fwd_model));
0028 %
0029 % Note that the value of fmin will be limited by the maximum and minimum
0030 % given in the opt.perturb vector.
0031 %
0032 % See also: INV_SOLVE_ABS_GN
0033 
0034 % (C) 2010-2013 Copyright Bartlomiej Grychtol, Andy Adler & Nolwenn Lesparre.
0035 % License: GPL version 2 or 3.
0036 % $Id: line_optimize.m 4596 2014-05-22 01:26:53Z alistair_boyle $
0037 
0038 
0039 if nargin <4 || isempty(opt)
0040    opt = struct;
0041 end
0042 
0043 opt = parse_options(opt);
0044 
0045 img = imgk;
0046 for i = 1:length(opt.perturb)
0047    img = calc_perturb(imgk,opt.perturb(i),dx, opt);
0048    vsim = fwd_solve(img);
0049    mlist(i) = feval(opt.objective_func,data0,vsim,imgk,img,opt);
0050 end
0051 
0052 % perform quadratic line fit in log space
0053 p10 = log10(opt.perturb);
0054 
0055 pf = polyfit(p10, mlist, 2);
0056 fmin = -pf(2)/pf(1)/2; % poly minimum for a 2nd order poly
0057 val = polyval(pf, fmin);
0058 
0059 if val > min(mlist) % fit didn't find the minimum, correct
0060    [mlist_o ik] = sort(mlist);
0061    % in case the user provided a 0 perturbation, avoid it
0062    if opt.perturb( ik(1)) ~= 0
0063       fmin = opt.perturb(1);
0064    else
0065       fmin = opt.perturb(2);
0066    end
0067 end
0068 
0069 fmin = 10^fmin; % convert back to linear
0070 
0071 % limit to the values given in perturb
0072 if fmin < min( opt.perturb )
0073    fmin = min( opt.perturb );
0074 end
0075 
0076 if fmin > max( opt.perturb )
0077    fmin = max( opt.perturb );
0078 end
0079 
0080 % RETURN VALUES
0081 img  = calc_perturb(imgk,fmin,dx,opt);
0082 vsim = fwd_solve(img);
0083 res  = feval(opt.objective_func,data0,vsim,imgk,img,opt);
0084 
0085 
0086 
0087 function img = calc_perturb(imgk, p, dx, opt)
0088    imgk = data_mapper(imgk);
0089    img = imgk;
0090    img.elem_data = imgk.elem_data + p*dx;
0091    img = apply_limits(img,opt);
0092    img = data_mapper(img, 1);
0093 
0094 
0095 function img = apply_limits(img,opt)
0096    img.elem_data(img.elem_data > opt.max_value) = opt.max_value;
0097    img.elem_data(img.elem_data < opt.min_value) = opt.min_value;
0098 
0099 
0100 function val = default_obj_fun(data0, data, img0, img, opt)
0101 dv = calc_difference_data(data0,data,img0.fwd_model);
0102 val = norm(dv);
0103 
0104 
0105 function opt = parse_options(opt)
0106 if ~isfield(opt,'perturb')
0107    opt.perturb = [0.1 0.5 1.0];
0108 end
0109 
0110 if ~isfield(opt,'min_value')
0111    opt.min_value = -Inf;
0112 end
0113 
0114 if ~isfield(opt,'max_value')
0115    opt.max_value = Inf;
0116 end
0117 
0118 if ~isfield(opt,'objective_func')
0119     opt.objective_func = @default_obj_fun;
0120 end