inv_solve_d_bar

D-Bar reconstruction
inv_model is an inverse model
 data1 is the inhomogene data
 data2 is backgounrd of 1 S/m

0001 function img = inv_solve_d_bar(inv_model, data1, data2)
0002 %D-Bar reconstruction
0003 %inv_model is an inverse model
0004 % data1 is the inhomogene data
0005 % data2 is backgounrd of 1 S/m
0006 
0007 % (C) 2023 Joeran Rixen. Licensed under GPL version 2 or 3
0008 
0009 %
0010 % CITATION_REQUEST:
0011 % AUTHOR: Rixen et al.
0012 % TITLE: The D-Bar Algorithm Fusing Electrical Impedance Tomography
0013 % with A Priori Radar Data: A Hands-On Analysis
0014 % JOURNAL: Algorithms
0015 % VOL: 16
0016 % NUM: 1
0017 % YEAR: 2023
0018 % LINK: https://doi.org/10.3390/a16010043
0019 % DOI: 10.3390/a16010043
0020 
0021 % TODO:
0022 % - remove elec_area hardcode
0023 % - tests for different stimulation patterns
0024 % - see how to remove stim_pattern_mat
0025 
0026 
0027     if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); test_d_bar; return; end
0028 
0029     citeme(mfilename);
0030 
0031     fmdl = inv_model.fwd_model;
0032     n_elec = num_elecs(fmdl);
0033     
0034     try
0035         trunc_rad = inv_model.hyperparameter.truncation_radius;
0036         grid_size = inv_model.hyperparameter.grid_size;
0037     catch
0038         error('This does not look like a valid d-bar model!')
0039     end
0040 
0041     %% Calculations regarding the Dirichlet to Neumann map
0042     stim_pattern_mat = create_stim_pattern_matrix(fmdl);
0043     stim_pattern_mat = stim_pattern_mat./sqrt(sum(stim_pattern_mat.^2,1));
0044 
0045     %% Calculations regarding the scattering transformation
0046     elec_pos_angle = get_elec_pos_angle(fmdl);
0047     
0048     %compute the grid size
0049     M = 2^grid_size;
0050     %step size between each point in K
0051     stepsize_grid = 2*trunc_rad/(M - 1);
0052     grid_base_vec = [-trunc_rad:stepsize_grid:trunc_rad];
0053     [K_x, K_y] = meshgrid(grid_base_vec);
0054     %creation of K
0055     K = K_x + 1i*K_y;
0056     %boolean mask of which points are actually inside the truncation radius
0057     K_inside_R = K( abs(K) < trunc_rad);
0058     
0059     del_L = dirichlet_neumann_difference(fmdl,data1,data2);
0060     
0061     %computation of the scattering transform
0062     t_aporx = zeros(length(K_inside_R),1);
0063     for i = 1:length(K_inside_R)
0064         k = K_inside_R(i);
0065         c_k = (full(stim_pattern_mat))'*exp(1i*k*exp(1i*elec_pos_angle));%%electrode_pos
0066         d_k = (full(stim_pattern_mat))'*exp(1i*conj(k)*conj(exp(1i*elec_pos_angle)));
0067     
0068         t_aporx(i) = conj(d_k)'*(del_L*c_k);
0069     end
0070     
0071     %% Setting up of the calculations for solving the D-Bar equation
0072     %calculation of the slightly larger Grid Q (comapred to K) for the Greens function
0073     begin_Q = trunc_rad*((2*M - 1)/(M - 1));
0074     grid_vec_Q = [-begin_Q:2*trunc_rad/(M-1):begin_Q];
0075     [Q_x, Q_y] = meshgrid(grid_vec_Q);
0076     Q = Q_x + 1i*Q_y;
0077     
0078     % factor that decides how big the smoothing zone outside of R is
0079     smooth_zone_factor = 1/4;
0080     % calcualte G bar
0081     G_bar = 1./(pi*Q);
0082     % set everything to 0 outisde the smoothing zone
0083     G_bar(abs(Q) > trunc_rad*(2 + smooth_zone_factor)) = 0;
0084     % get index of the smoothing zone
0085     smooth_zone = abs(Q) > trunc_rad*(2) & abs(Q) < trunc_rad*(2 + smooth_zone_factor);
0086     % apply the smoothing (in this case linear) to G_bar
0087     G_bar(smooth_zone) = G_bar(smooth_zone).*((trunc_rad*(2 + smooth_zone_factor) - abs(Q(smooth_zone)))/(trunc_rad*smooth_zone_factor));
0088     
0089     % make the scattering transform compatible to the larger grid Q
0090     scat = zeros(size(K));
0091     rad_idx = abs(K) < trunc_rad;
0092     scat(rad_idx) = t_aporx;
0093     scat = padarray_EIDORS(scat, [M/2 M/2], 0, 'both');
0094     
0095     % use the FFT trick by Mueller and Siltanen
0096     %(due to Matlabs indexing fftshift is needed before apllying fft2)
0097     G_bar_fft = fft2(fftshift(G_bar));
0098     
0099     % calcualte the first part of T (this is independnt of the location z)
0100     T_R_first = scat./(4*pi*conj(Q));
0101     
0102     
0103     % Stuff needed for the computation, this is not directly related to the D-Bar algorithm
0104     % index of which vlaues are inside R
0105     rad_idx_dbar = abs(Q) < trunc_rad;
0106     % number of points inside R
0107     num_pts_in_R = sum(rad_idx_dbar(:));
0108     
0109     % Construct right hand side of the Dbar equation as the initial guess
0110     right_side_eq = [ones(num_pts_in_R,1);zeros(num_pts_in_R,1)];
0111     
0112     recon_cords = interp_mesh(fmdl);
0113     % transform the reconstruction points into complex coordiantes
0114     recon_pts = recon_cords(:, 1) + 1i*recon_cords(:, 2);
0115     num_rec_pts  = length(recon_pts);
0116     
0117     mu_final = ones(num_rec_pts,1);
0118     iniguess = [ones(num_pts_in_R,1);zeros(num_pts_in_R,1)];
0119     
0120     
0121     middle_val = size(T_R_first, 1)/2;
0122     % When you have the parralel computation toolkit feel free to use parfor instead of for
0123     for iii = 1:num_rec_pts
0124         % get positions for recosntruction
0125         z = recon_pts(iii);
0126     
0127         % multiplicator for the D-Bar function
0128         T_R = T_R_first.*(exp(-1i*(Q*z+conj(Q*z))));
0129     
0130         % Solve D-Bar equation; Note that instead of a matrix a function is
0131         % given, the input aargument for that function follow after the
0132         % second []
0133         [sol, ~] = gmres(@DB_oper_demyst, right_side_eq, 40, 1e-4, 450, [], [], iniguess, M*2, rad_idx_dbar, num_pts_in_R, 2*trunc_rad/(M - 1), G_bar_fft, T_R);
0134     
0135         %use sol to get mu for the recosntruction
0136         mu = zeros(size(Q));
0137         mu(rad_idx_dbar) = sol(1:num_pts_in_R) + 1i*sol((num_pts_in_R+1):end);
0138     
0139         % get mu from interpolation of the center
0140         mu_final(iii) = mean([mu(middle_val,middle_val), mu(middle_val+1,middle_val), mu(middle_val, middle_val+1), mu(middle_val+1, middle_val+1)]);
0141     
0142     end
0143     img = struct;
0144     img.name= ['solved by D-Bar'];
0145     img.elem_data = get_conductivity(mu_final);
0146     img.fwd_model= fmdl;
0147 end
0148 
0149 
0150 function stim_pattern_mat = create_stim_pattern_matrix(fmdl)
0151     % converts the stimpattern in the model into matrix form for ease of
0152     % computation
0153     n_elec = length(fmdl.electrode);
0154     stim_pattern_mat = zeros(n_elec, n_elec-1);
0155     
0156     for i = 1:n_elec-1
0157         stim_pattern_mat(:, i) = full(fmdl.stimulation(i).stim_pattern);
0158     end
0159 end
0160 
0161 function elec_pos_angle = get_elec_pos_angle(fmdl)
0162     % calculates the position of the electrodes with respect to it's angle
0163     n_elec = length(fmdl.electrode);
0164     elec_pos_cart = zeros(n_elec, 2);
0165     
0166     for i = 1:n_elec
0167         elec_pos_cart(i, :) = [mean(fmdl.nodes(fmdl.electrode(i).nodes, 1)), mean(fmdl.nodes(fmdl.electrode(i).nodes, 2))];
0168     end
0169     
0170     elec_pos_cart = [elec_pos_cart(:, 1) - mean(elec_pos_cart(:, 1)), elec_pos_cart(:, 2) - mean(elec_pos_cart(:, 2))];
0171     
0172     elec_pos_angle = atan2(elec_pos_cart(:, 2), elec_pos_cart(:, 1));
0173     
0174     elec_pos_angle(elec_pos_angle < 0) = elec_pos_angle(elec_pos_angle < 0) + 2*pi;
0175 end
0176 
0177 function B = padarray_EIDORS(A,padsize,padval,direction)
0178 % padarray function is not availabe in standard matlan.
0179     if nargin < 4
0180         direction = 'both';
0181     end
0182     if nargin < 3
0183         padval = 0;
0184     end
0185     if strcmp(direction,'both')
0186         B = padval*ones(size(A)+2*padsize);
0187         B(padsize(1)+1:end-padsize(1),padsize(2)+1:end-padsize(2)) = A;
0188     elseif strcmp(direction,'pre')
0189         B = padval*ones(size(A)+padsize);
0190         B(padsize(1)+1:end,padsize(2)+1:end) = A;
0191     elseif strcmp(direction,'post')
0192         B = padval*ones(size(A)+padsize);
0193         B(1:size(A,1),1:size(A,2)) = A;
0194     end
0195 end
0196 
0197 function [conductivity] = get_conductivity(mu)
0198     %computes the conductivty from the solution of the D-Bar equation:
0199     
0200     %mu <number of points to reconstruct>:
0201     %solution of the D-Bar equation
0202     
0203     conductivity = real(mu.^2);
0204 end
0205 
0206 function [fmdl, inv_model] = mk_dbar_model(num_electrs)
0207     % creates a forward model compatible to the d-bar algorithm
0208     % and creates a backwards model for the d-bar algorithm
0209     fmdl = ng_mk_cyl_models(.5,[num_electrs,0.25],[0.075]);
0210     [stim, ~]= mk_stim_patterns(num_electrs, 1, '{trig}', '{mono}', {'meas_current'}, 0.1);
0211     fmdl.stimulation = stim;
0212     
0213     fmdl_2d = ng_mk_cyl_models([0,1,0.05],[num_electrs,0],[0.075]);
0214     fmdl_2d.stimulation = stim;
0215 
0216     inv_model = struct;
0217     inv_model.name = 'EIDORS D-Bar model';
0218     inv_model.solve = 'inv_solve_d_bar';
0219     inv_model.RtR_prior = 'eidors_default';
0220 
0221     % set multiple hyperparameters
0222     inv_model.hyperparameter.truncation_radius = 3.0;
0223     inv_model.hyperparameter.grid_size = 5;
0224    
0225     inv_model.jacobian_bkgnd = 1;
0226     inv_model.type = 'inv_model';
0227     inv_model.fwd_model = fmdl_2d;
0228 end
0229 
0230 function result = DB_oper_demyst(sol, M, rad_idx_dbar, num_pts_in_R, h, G_bar_fft, TR)
0231     %rearranges the solution of the D-Bar algorithm to fit the square, perform
0232     %the D-Bar operation and transform it back:
0233 
0234     %sol <(total number of gridpoints)*2>: solution of the D-Bar equation so
0235     %far
0236     
0237     %M <(size of Q)*2>: size of Q
0238     
0239     %rad_idx_dbar <size of Q>*<size of Q>: indices whether grid points are
0240     %inside R or not
0241     
0242     %num_pts_in_R <1>: number of points inside R
0243     
0244     %h <1>: spacing of the grid
0245     
0246     %G_bar_fft <size of Q>*<size of Q>:  FFT trick operator
0247     
0248     
0249     sol_square = zeros(M, M);
0250     sol_square(rad_idx_dbar) = sol(1:num_pts_in_R) + 1i*sol((num_pts_in_R+1):end);
0251     
0252     % Apply real-linear operator %%eqn. 15.65 from "Linear- and non_liner
0253     % inverse problems"
0254     result = sol_square - h^2*ifft2(G_bar_fft .* fft2( TR.*conj(sol_square) ));
0255     
0256     % Construct result as a vector with real and imaginary parts separate
0257     % %%matrl
0258     result = [real(result(rad_idx_dbar)); imag(result(rad_idx_dbar))];
0259 end
0260 
0261 function test_d_bar
0262     [fmdl, inv_model] = mk_dbar_model(16);
0263     
0264     img = mk_image(fmdl, 1);
0265     vh = fwd_solve(img);
0266     img_inh = set_test_target(img);
0267     vi = fwd_solve(img_inh);
0268     
0269     img_solved = inv_solve(inv_model, vi, vh);
0270     show_fem(img_solved);
0271 end
0272 
0273 function img_out = set_test_target(img)
0274     cog = interp_mesh(img.fwd_model, 0);
0275     tar_offset = 0.6;
0276     cog(:,1,:) = cog(:,1,:) - tar_offset;
0277     tar_rad = 0.3;
0278     idx_tar = cog(:,1,:).^2 + cog(:,2,:).^2 < tar_rad.^2;
0279     img.elem_data(idx_tar) = img.elem_data(1)*1.2;
0280     img_out = img;
0281 end