inv_solve_abs_pdipm

PURPOSE ^

INV_SOLVE_ABS_PDIPM inverse solver for absolute data using Primal/Dual interior point method

SYNOPSIS ^

function img=inv_solve_abs_pdipm( inv_model, data);

DESCRIPTION ^

 INV_SOLVE_ABS_PDIPM  inverse solver for absolute data using Primal/Dual interior point method
 img= inv_solve_abs_pdipm( inv_model, data);
 img        => output image (or vector of images)
 inv_model  => inverse model struct
 data       => vector of eit data

  inv_model.inv_solve_abs_pdipm.norm_data  1 or 2 (DEFAULT 2)
  inv_model.inv_solve_abs_pdipm.norm_prior 1 or 2 (DEFAULT 2)
  inv_model.inv_solve_abs_pdipm.beta     (default 1e-6)

 Parameters:
  max_iter =  inv_model.parameters.max_iterations (default 10)
      Max number of iterations before stopping
  min change = inv_model.parameters.min_change   (default 0)
      Min Change in objective fcn (norm(y-Jx)^2 + hp*TV(x)) before stopping
 beta is the parameter that smooths the TV functional

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function img=inv_solve_abs_pdipm( inv_model, data);
0002 % INV_SOLVE_ABS_PDIPM  inverse solver for absolute data using Primal/Dual interior point method
0003 % img= inv_solve_abs_pdipm( inv_model, data);
0004 % img        => output image (or vector of images)
0005 % inv_model  => inverse model struct
0006 % data       => vector of eit data
0007 %
0008 %  inv_model.inv_solve_abs_pdipm.norm_data  1 or 2 (DEFAULT 2)
0009 %  inv_model.inv_solve_abs_pdipm.norm_prior 1 or 2 (DEFAULT 2)
0010 %  inv_model.inv_solve_abs_pdipm.beta     (default 1e-6)
0011 %
0012 % Parameters:
0013 %  max_iter =  inv_model.parameters.max_iterations (default 10)
0014 %      Max number of iterations before stopping
0015 %  min change = inv_model.parameters.min_change   (default 0)
0016 %      Min Change in objective fcn (norm(y-Jx)^2 + hp*TV(x)) before stopping
0017 % beta is the parameter that smooths the TV functional
0018 
0019 % (C) 2010 Andrea Borsic + Andy Adler. License: GPL v2 or v3
0020 % $Id: inv_solve_abs_pdipm.m 5081 2015-06-01 21:20:16Z michaelcrabb30 $
0021 
0022 
0023 pp= process_parameters(inv_model);
0024 
0025 img_bkgnd = homogeneous_estimate( inv_model, data );
0026 
0027 alpha=calc_hyperparameter( inv_model );
0028 L= calc_R_prior( inv_model );
0029 W= calc_meas_icov( inv_model );
0030 
0031 img= feval(pp.fn, img_bkgnd, W,alpha*L,data, pp);
0032 
0033 img.name = sprintf('inv_solve_abs_pdipm-nd%d-ni%d',pp.norm_data,pp.norm_image);
0034 
0035 % This is the Gauss-Newton algorithm
0036 %   for the linear case it is: s= (J'*W*J + L'*L)\J'*W*d;
0037 function img= pdipm_2_2(  img,W,L,d, pp);
0038    img0 = img;
0039    hp2RtR = L'*L;
0040    for i = 1:pp.max_iter
0041      dv = sim_diff( img, d);
0042      J = calc_jacobian( img );
0043 
0044      RDs = hp2RtR*(img0.elem_data - img.elem_data);
0045      ds = (J'*W*J + hp2RtR)\(J'*dv + RDs);
0046 
0047      img = line_optimize(img, ds, d);
0048 
0049      pp = manage_beta(pp);
0050      loop_display(i)
0051    end
0052 
0053 function img= pdipm_1_2( img,W,L,d, pp);
0054    [M]   = size(W,1); % M measurements
0055    [jnk,N] = size(L); % E edges, N parameters
0056    x= zeros( M, 1 ); % dual var - start with zeros
0057 
0058    for loop = 1:pp.max_iter
0059       % Jacobian
0060       J = calc_jacobian( img );
0061       dv = -sim_diff(img, d);
0062       s = img.elem_data;
0063 
0064       % Define variables
0065       f = dv;                  F= spdiags(f,0,M,M);
0066                                X= spdiags(x,0,M,M);
0067       e = sqrt(f.^2 + pp.beta);E= spdiags(e,0,M,M);
0068 
0069       % Define derivatives
0070       dFc_ds = (speye(M,M) - X*inv(E)*F)*J;
0071       dFc_dx = -E;
0072       dFf_ds = L'*L;
0073       dFf_dx = J'*W;
0074 
0075       dsdx = -[dFc_ds, dFc_dx; dFf_ds, dFf_dx] \ ...
0076               [ f-E*x; J'*W*x + L'*L*s ];
0077 
0078       ds = dsdx(1:N);
0079       img = line_optimize(img, ds, d);
0080 
0081       dx = x_update(x, dsdx(N+(1:M)));
0082       x= x + dx;
0083  
0084       pp = manage_beta(pp);
0085       loop_display(i)
0086    end
0087 
0088 function img= pdipm_2_1(img,W,L,d, pp);
0089    [M]   = size(W,1); % M measurements
0090    [G,N] = size(L); % E edges, N parameters
0091    x= zeros( G, 1 ); % dual var - start with zeros
0092 
0093    for loop = 1:pp.max_iter
0094       % Jacobian
0095       J = calc_jacobian( img );
0096       dv = -sim_diff(img, d);
0097       
0098       % Define variables
0099       f = L*img.elem_data;     F= spdiags(f,0,G,G);
0100                                X= spdiags(x,0,G,G);
0101       e = sqrt(f.^2 + pp.beta);E= spdiags(e,0,G,G);
0102 
0103       % Define derivatives
0104       dFc_ds = (speye(G,G) - X*inv(E)*F)*L;
0105       dFc_dx = -E;
0106       dFf_ds = J'*J;
0107       dFf_dx = L';
0108 
0109       dsdx = -[dFc_ds, dFc_dx; dFf_ds, dFf_dx] \ ...
0110               [ f-E*x; J'*dv + L'*x ];
0111 
0112       ds = dsdx(1:N);
0113       img = line_optimize(img, ds, d);
0114 
0115       dx = x_update(x, dsdx(N+(1:G)));
0116       x= x + dx;
0117 
0118       pp = manage_beta(pp);
0119       loop_display(i)
0120    end
0121 
0122 %   img0 = img;
0123 %   hp2RtR = L'*L;
0124 %   for i = 1:pp.max_iter
0125 %     vsim = fwd_solve( img );
0126 %     dv = calc_difference_data( vsim , d, img.fwd_model);
0127 %     J = calc_jacobian( img );
0128 %
0129 %     RDx = hp2RtR*(img0.elem_data - img.elem_data);
0130 %     dx = (J'*W*J + hp2RtR)\(J'*dv + RDx);
0131 %
0132 %     img = line_optimize(img, dx, d);
0133 %
0134 %     loop_display(i)
0135 %   end
0136 
0137 function img= pdipm_1_1( img,W,L,d, pp);
0138    [M]   = size(W,1); % M measurements
0139    [D,N] = size(L); % E edges, N parameters
0140    s= img.elem_data;
0141    y= zeros( D, 1 ); % dual var - start with zeros
0142    x= zeros( M, 1 ); % dual var - start with zeros
0143 
0144    for loop = 1:pp.max_iter
0145       % Jacobian
0146       J = calc_jacobian( img );
0147       dv = -sim_diff(img, d);
0148 
0149       % Define variables
0150       g = L*img.elem_data;     G= spdiags(g,0,D,D);
0151       s = sqrt(g.^2 + pp.beta);S= spdiags(s,0,D,D);
0152                                Y= spdiags(y,0,D,D);
0153 
0154       f = dv;                  F= spdiags(f,0,M,M);
0155       e = sqrt(f.^2 + pp.beta);E= spdiags(e,0,M,M);
0156                                X= spdiags(x,0,M,M);
0157 
0158       % Define derivatives
0159       As1 = sparse(N,N);
0160       As2 = (speye(M,M) - X/E*F) * J; % As2 = (speye(M,M) - X*inv(E)*F) * J;
0161       As3 = (speye(D,D) - Y/S*G) * L; % As3 = (speye(D,D) - Y*inv(S)*G) * L;
0162       Ax1 = J'*W;
0163       Ax2 = -E;
0164       Ax3 = sparse(D,M);
0165       Ay1 = L';
0166       Ay2 = sparse(M,D);
0167       Ay3 = -S;
0168       B1  = J'*W*x + L'*y;
0169       B2  = f - E*x;
0170       B3  = g - S*y;
0171 
0172 %     DD = -[As1,Ax1,Ay1; ...
0173 %            As2,Ax2,Ay2; ...
0174 %            As3,Ax3,Ay3] \ [B1;B2;B3];
0175 
0176 %     dm = DD(1:N); dx = DD(N+(1:M)); dy = DD(N+M+(1:D));
0177 
0178       JtWiE = J'*W/E; LtiS = L'/S;
0179       dm= -(JtWiE*As2 + LtiS*As3)\(JtWiE*f + LtiS*g);
0180       dx= E\(B2 + As2*dm);
0181       dy= S\(B3 + As3*dm);
0182 
0183       img = line_optimize(img, dm, d);
0184 
0185       dx = x_update(x, dx);
0186       x= x + dx;
0187 
0188       dy = x_update(y, dy);
0189       y= y + dy;
0190 
0191       loop_display(i)
0192       pp = manage_beta(pp);
0193    end
0194 
0195 % abs(x + dx) must be <= 1
0196 function dx = x_update( x, dx)
0197    dx(dx==0) = eps; % can't have zeros
0198    sx = sign(dx);
0199    % space to limits in direction of x
0200    clr = sx - x;
0201    % how much to multiply by to get to limits
0202    fac = clr./dx;
0203    % choose min amount to get to limits
0204    dx = dx*min(fac);
0205 
0206 function pp = manage_beta(pp);
0207    pp.beta = pp.beta * pp.beta_reduce;
0208    if pp.beta < pp.beta_minimum;
0209       pp.beta = pp.beta_minimum;
0210    end
0211 
0212 function pp= process_parameters(imdl);
0213    try    pp.max_iter = imdl.parameters.max_iterations;
0214    catch  pp.max_iter = 10;
0215    end
0216 
0217    try    pp.min_change = imdl.parameters.min_change;
0218    catch  pp.min_change = 0;
0219    end
0220 
0221    try    pp.beta = imdl.inv_solve_abs_pdipm.beta; 
0222    catch  pp.beta = 1e-6;
0223    end
0224 
0225    pp.beta_reduce = 0.2;
0226    pp.beta_minimum= 1e-16;
0227 
0228    try    pp.norm_data = imdl.inv_solve_abs_pdipm.norm_data;
0229    catch  pp.norm_data = 2;
0230    end
0231 
0232    try    pp.norm_image = imdl.inv_solve_abs_pdipm.norm_image;
0233    catch  pp.norm_image = 2;
0234    end
0235 
0236    if     pp.norm_data==2 && pp.norm_image==2;
0237       pp.fn = @pdipm_2_2;
0238    elseif pp.norm_data==2 && pp.norm_image==1;
0239       pp.fn = @pdipm_2_1;
0240    elseif pp.norm_data==1 && pp.norm_image==2;
0241       pp.fn = @pdipm_1_2;
0242    elseif pp.norm_data==1 && pp.norm_image==1;
0243       pp.fn = @pdipm_1_1;
0244    else
0245       error('norm_data and norm_image should be 1 or 2');
0246    end
0247 
0248 
0249 
0250 % Fit a parabola to the linefit and pick the best point
0251 % This is faster than doing an exhaustive search
0252 function  img = line_optimize(imgk, dx, data1);
0253   flist = [ 0.1,  0.5, 1.0];
0254   clim = mean(imgk.elem_data)/10; % prevent zero and negative conductivity
0255   img = imgk;
0256   for i = 1:length(flist);
0257      img.elem_data = imgk.elem_data + flist(i)*dx;
0258      img.elem_data(img.elem_data <= clim ) = clim;
0259      dv = sim_diff( img, data1);
0260      mlist(i) = norm(dv);
0261   end
0262   pf = polyfit(flist, mlist, 2);
0263   fmin = -pf(2)/pf(1)/2; % poly minimum
0264   fmin(fmin>1) = 1; fmin(fmin<0) = 0;
0265 
0266   img.elem_data = imgk.elem_data + flist(i)*dx;
0267   img.elem_data(img.elem_data <= clim ) = clim;
0268 
0269 function img = homogeneous_estimate( imdl, data );
0270    img = calc_jacobian_bkgnd( imdl );
0271    vs = fwd_solve(img);
0272    data = data_vector( data, imdl );
0273 
0274    pf = polyfit(data,vs.meas,1);
0275 
0276    img.elem_data = img.elem_data*pf(1);
0277 
0278 function data = data_vector( data, imdl );
0279    if isstruct(data)
0280       data = data.meas;
0281    else
0282      meas_select = [];
0283      try
0284         meas_select = imdl.fwd_model.meas_select;
0285      end
0286      if length(data) == length(meas_select)
0287         data = data(meas_select);
0288      end
0289    end
0290 
0291 function dv = sim_diff( img, data1);
0292   vsim = fwd_solve( img );
0293   dv = calc_difference_data( vsim , data1, img.fwd_model);
0294 
0295 function loop_display(i)
0296    fprintf('+');
Generated on Fri 30-Dec-2022 19:44:54 by m2html © 2005