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