Description of pdipm

0001 function img=pdipm_abs( inv_model, data);
0002 % PDIPM_ABS  inverse solver for absolute data using Primal/Dual interior point method
0003 % img= pdipm_abs( 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.pdipm_abs.norm_data  1 or 2 (DEFAULT 2)
0009 %  inv_model.pdipm_abs.norm_prior 1 or 2 (DEFAULT 2)
0010 %  inv_model.pdipm_abs.beta     (default 1e-6)
0011 %
0012 % Parameters:
0013 %  max_iter =  inv_model.parameters.max_iteration (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: pdipm_abs.html 2819 2011-09-07 16:43:11Z aadler $
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('pdipm_abs-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       r = sqrt(g.^2 + pp.beta);R= spdiags(r,0,D,D); % S in paper
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*inv(E)*F) * J;
0161       As3 = (speye(D,D) - Y*inv(R)*G) * L;
0162       Ax1 = J'*W;
0163       Ax2 = -E;
0164       Ax3 = sparse(D,M);
0165       Ay1 = L';
0166       Ay2 = sparse(M,D);
0167       Ay3 = -R;
0168       B1  = J'*W*x + L'*y;
0169       B2  = f - E*x;
0170       B3  = g - R*y;
0171 
0172       DD = -[As1,Ax1,Ay1; ...
0173              As2,Ax2,Ay2; ...
0174              As3,Ax3,Ay3] \ [B1;B2;B3];
0175 
0176       ds = DD(1:N);
0177       img = line_optimize(img, ds, d);
0178 
0179       dx = x_update(x, DD(N+(1:M)));
0180       x= x + dx;
0181 
0182       dy = x_update(y, DD(N+M+(1:D)));
0183       y= y + dy;
0184 
0185       loop_display(i)
0186       pp = manage_beta(pp);
0187    end
0188 
0189 % abs(x + dx) must be <= 1
0190 function dx = x_update( x, dx)
0191    dx(dx==0) = eps; % can't have zeros
0192    sx = sign(dx);
0193    % space to limits in direction of x
0194    clr = sx - x;
0195    % how much to multiply by to get to limits
0196    fac = clr./dx;
0197    % choose min amount to get to limits
0198    dx = dx*min(fac);
0199 
0200 function pp = manage_beta(pp);
0201    pp.beta = pp.beta * pp.beta_reduce;
0202    if pp.beta < pp.beta_minimum;
0203       pp.beta = pp.beta_minimum;
0204    end
0205 
0206 function pp= process_parameters(imdl);
0207    try    pp.max_iter = imdl.parameters.max_iterations;
0208    catch  pp.max_iter = 10;
0209    end
0210 
0211    try    pp.min_change = imdl.parameters.min_change;
0212    catch  pp.min_change = 0;
0213    end
0214 
0215    try    pp.beta = imdl.pdipm_abs.beta; 
0216    catch  pp.beta = 1e-6;
0217    end
0218 
0219    pp.beta_reduce = 0.2;
0220    pp.beta_minimum= 1e-16;
0221 
0222    try    pp.norm_data = imdl.pdipm_abs.norm_data;
0223    catch  pp.norm_data = 2;
0224    end
0225 
0226    try    pp.norm_image = imdl.pdipm_abs.norm_image;
0227    catch  pp.norm_image = 2;
0228    end
0229 
0230    if     pp.norm_data==2 && pp.norm_image==2;
0231       pp.fn = @pdipm_2_2;
0232    elseif pp.norm_data==2 && pp.norm_image==1;
0233       pp.fn = @pdipm_2_1;
0234    elseif pp.norm_data==1 && pp.norm_image==2;
0235       pp.fn = @pdipm_1_2;
0236    elseif pp.norm_data==1 && pp.norm_image==1;
0237       pp.fn = @pdipm_1_1;
0238    else
0239       error('norm_data and norm_image should be 1 or 2');
0240    end
0241 
0242 
0243 
0244 % Fit a parabola to the linefit and pick the best point
0245 % This is faster than doing an exhaustive search
0246 function  img = line_optimize(imgk, dx, data1);
0247   flist = [ 0.1,  0.5, 1.0];
0248   clim = mean(imgk.elem_data)/10; % prevent zero and negative conductivity
0249   img = imgk;
0250   for i = 1:length(flist);
0251      img.elem_data = imgk.elem_data + flist(i)*dx;
0252      img.elem_data(img.elem_data <= clim ) = clim;
0253      dv = sim_diff( img, data1);
0254      mlist(i) = norm(dv);
0255   end
0256   pf = polyfit(flist, mlist, 2);
0257   fmin = -pf(2)/pf(1)/2; % poly minimum
0258   fmin(fmin>1) = 1; fmin(fmin<0) = 0;
0259 
0260   img.elem_data = imgk.elem_data + flist(i)*dx;
0261   img.elem_data(img.elem_data <= clim ) = clim;
0262 
0263 function img = homogeneous_estimate( imdl, data );
0264    img = calc_jacobian_bkgnd( imdl );
0265    vs = fwd_solve(img);
0266    data = data_vector( data, imdl );
0267 
0268    pf = polyfit(data,vs.meas,1);
0269 
0270    img.elem_data = img.elem_data*pf(1);
0271 
0272 function data = data_vector( data, imdl );
0273    if isstruct(data)
0274       data = data.meas;
0275    else
0276      meas_select = [];
0277      try
0278         meas_select = imdl.fwd_model.meas_select;
0279      end
0280      if length(data) == length(meas_select)
0281         data = data(meas_select);
0282      end
0283    end
0284 
0285 function dv = sim_diff( img, data1);
0286   vsim = fwd_solve( img );
0287   dv = calc_difference_data( vsim , data1, img.fwd_model);
0288 
0289 function loop_display(i)
0290    fprintf('+');
pdipm_abs

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
