Add files via upload

Functions/ringingRemovalFilt.m

@@ -0,0 +1,36 @@
+function img_out = ringingRemovalFilt(xin, yin, img_in, c0, f, cutoff, ord)
+%RINGINGREMOVALFILT Filter to Remove Ringing from Waveform Inversion Result
+%
+% img_out = ringingRemovalFilt(xin, yin, img_in, c0, f, cutoff)
+% INPUT:
+%   xin = x (N-element array) grid [in m] for the input image
+%   yin = y (M-element array) grid [in m] for the input image
+%   img_in = input image (M x N array)
+%   c0 = reference sound speed [m/s]
+%   f = frequency [Hz]
+%   cutoff = ranging from 0 (only keep DC) to 1 (up to ringing frequency)
+%   ord = order of radial Butterworth filter cutoff (Inf for sharp cutoff)
+% OUTPUT:
+%   img_out = output image (M x N array)
+
+% Number and spacing of input points (assumes a uniform spacing)
+Nxin = numel(xin); dxin = mean(diff(xin));
+Nyin = numel(xin); dyin = mean(diff(yin));
+
+% K-Space
+kxin = fftshift(((0:Nxin-1)/Nxin)/dxin);
+kxin(kxin>=1/(2*dxin)) = kxin(kxin>=1/(2*dxin)) - 1/dxin;
+kyin = fftshift(((0:Nyin-1)/Nyin)/dyin);
+kyin(kyin>=1/(2*dyin)) = kyin(kyin>=1/(2*dyin)) - 1/dyin;
+[Kxin, Kyin] = meshgrid(kxin, kyin);
+
+% Cutoff Wavenumber
+k = 2*f/c0; kcutoff = k*cutoff;
+
+% Filter Based on Cutoff
+radialFilt = 1./(1+((Kxin.^2 + Kyin.^2)/(kcutoff^2)).^ord);
+
+% Apply Filter
+img_out = real(ifft2(ifftshift(radialFilt.*(fftshift(fft2(img_in))))));
+
+end

Functions/stencilOptParams.m

@@ -0,0 +1,42 @@
+function [b,d,e] = stencilOptParams(vmin,vmax,f,h,g)
+%STENCILOPTPARAMS Optimal Params for 9-Point Stencil 
+%   INPUTS:
+%       vmin = minimum wave velocity [L/T]
+%       vmax = maximum wave velocity [L/T]
+%       f = frequency [1/T]
+%       h = grid spacing in X [L]
+%       g = (grid spacing in Y [L])/(grid spacing in X [L])
+%   OUTPUTS:
+%       b, d, e = optimal params according to Chen/Cheng/Feng/Wu 2013 Paper
+
+l = 100; r = 10;
+Gmin = vmin/(f*h); Gmax = vmax/(f*h);
+
+m = 1:l; n = 1:r;
+theta = (m-1)*pi/(4*(l-1));
+G = 1./(1/Gmax + ((n-1)/(r-1))*(1/Gmin-1/Gmax));
+
+[TH, GG] = meshgrid(theta, G);
+
+P = cos(g*2*pi*cos(TH)./GG);
+Q = cos(2*pi*sin(TH)./GG);
+
+S1 = (1+1/(g^2))*(GG.^2).*(1-P-Q+P.*Q);
+S2 = (pi^2)*(2-P-Q);
+S3 = (2*pi^2)*(1-P.*Q);
+S4 = 2*pi^2 + (GG.^2).*((1+1/(g^2))*P.*Q-P-Q/(g^2));
+
+fixB = true;
+if fixB
+    b = 5/6; % Fix the Value to 5/6 based on Laplacian Derived by Robert E. Lynch
+    A = [S2(:), S3(:)]; y = S4(:)-b*S1(:);
+    params = (A'*A)\(A'*y);
+    d = params(1); e = params(2);
+else
+    A = [S1(:), S2(:), S3(:)]; y = S4(:);
+    params = (A'*A)\(A'*y);
+    b = params(1); d = params(2); e = params(3);
+end
+
+end
+

MultiFrequencyWaveformInvKCI.m

@@ -0,0 +1,298 @@
+clear
+clc
+
+% Add Functions to Path
+addpath(genpath('Functions'));
+
+% Load Extracted Dataset
+filename = 'Malignancy'; % 'BenignCyst', 'Malignancy'
+load(['SampleData/', filename, '.mat'], ...
+    'time', 'transducerPositionsXY', 'full_dataset');
+
+% Assemble Full Dataset
+numElements = size(transducerPositionsXY,2);
+fsa_rf_data = single(full_dataset);
+clearvars full_dataset;
+
+% Extract the Desired Frequency for Waveform Inversion
+fDATA_SoS = (0.3:0.05:1.25)*(1e6); % Frequencies for SoS-only Iterations [Hz]
+fDATA_SoSAtten = (0.325:0.05:1.275)*(1e6); % Frequencies for SoS/Attenuation Iterations [Hz]
+fDATA = [fDATA_SoS, fDATA_SoSAtten]; % All Frequencies [Hz]
+
+% Attenuation Iterations Always Happen After SoS Iterations
+niterSoSPerFreq = [3*ones(size(fDATA_SoS)), 3*ones(size(fDATA_SoSAtten))]; % Number of Sound Speed Iterations Per Frequency
+niterAttenPerFreq = [0*ones(size(fDATA_SoS)), 3*ones(size(fDATA_SoSAtten))]; % Number of Attenuation Iterations Per Frequency
+
+% Discrete-Time Fourier Transform (DTFT) - Not an FFT though!
+DTFT = exp(-1i*2*pi*fDATA'*time)*mean(diff(time));
+
+% Geometric TOF Based on Sound Speed
+c_geom = 1540; % Sound Speed [mm/us] 
+transducerPositionsX = transducerPositionsXY(1,:);
+transducerPositionsY = transducerPositionsXY(2,:);
+geomTOFs = zeros(size(transducerPositionsXY,2));
+for col = 1:size(transducerPositionsXY,2)
+    geomTOFs(:,col) = ...
+        sqrt((transducerPositionsX-transducerPositionsX(col)).^2 + ...
+        (transducerPositionsY-transducerPositionsY(col)).^2)/c_geom;
+end
+clearvars col transducerPositionsXY;
+
+% Window Time Traces - Extract the Frequencies for Waveform Inversion
+REC_DATA = zeros(numElements, numElements, numel(fDATA), 'single');
+sos_perc_change_pre = 0.05; sos_perc_change_post = Inf; 
+twinpre = sos_perc_change_pre*max(geomTOFs(:));
+twinpost = sos_perc_change_post*max(geomTOFs(:));
+for tx_element = 1:numElements
+    times_tx = geomTOFs(tx_element,:);
+    [TIMES_TX, TIME] = meshgrid(times_tx, time);
+    window = exp(-(1/2)*(subplus(TIME-TIMES_TX)/twinpost + ...
+        subplus(TIMES_TX-TIME)/twinpre).^2);
+    REC_DATA(tx_element,:,:) = ...
+        permute(DTFT*(window.*fsa_rf_data(:,:,tx_element)),[2,1]);
+end
+clearvars twin window fsa_rf_data time times_tx TIME TIMES_TX DTFT;
+
+%% Waveform Inversion
+
+% Create Sound Speed Map and Transducer Ring
+dxi = 0.3e-3; xmax = 120e-3;
+xi = -xmax:dxi:xmax; yi = xi;
+Nxi = numel(xi); Nyi = numel(yi);
+[Xi, Yi] = meshgrid(xi, yi);
+x_circ = transducerPositionsX;
+y_circ = transducerPositionsY;
+x_idx = dsearchn(xi(:), x_circ(:));
+y_idx = dsearchn(yi(:), y_circ(:));
+ind = sub2ind([Nyi, Nxi], y_idx, x_idx);
+msk = zeros(Nyi, Nxi); msk(ind) = 1;
+
+% Parameters
+h = dxi; % [m]
+g = 1; % (grid spacing in Y)/(grid spacing in X)
+alphaDB = 0.0; % Attenuation [dB/(MHz*mm)]
+alphaNp = (log(10)/20)*alphaDB*((1e3)/(1e6)); % Attenuation [Np/(Hz*m)]
+ATTEN = alphaNp*ones(Nyi,Nxi); % Make Spatially Varying Attenuation [Np/(Hz m)]
+% Solve Options for Helmholtz Equation
+sign_conv = -1; % Sign Convention
+a0 = 10; % PML Constant
+L_PML = 9.0e-3; % Thickness of PML  
+
+% Conversion of Units for Attenuation Map
+Np2dB = 20/log(10);
+slow2atten = (1e6)/(1e2); % Hz to MHz; m to cm
+
+% Compute Phase Screen to Correct for Discretization of Element Positions
+% Times for Discretized Positions Assuming Constant Sound Speed
+x_circ_disc = xi(x_idx); y_circ_disc = yi(y_idx);
+[x_circ_disc_tx, x_circ_disc_rx] = meshgrid(x_circ_disc, x_circ_disc);
+[y_circ_disc_tx, y_circ_disc_rx] = meshgrid(y_circ_disc, y_circ_disc);
+geomTOFs_disc = sqrt((x_circ_disc_tx-x_circ_disc_rx).^2 + ...
+    (y_circ_disc_tx-y_circ_disc_rx).^2)/c_geom;
+% Time-of-Flight Error for Discretized Positions
+geomTOFs_error = geomTOFs_disc-geomTOFs;
+
+% Phase Screen and Gain Correction
+PS = zeros(numElements, numElements, numel(fDATA));
+for f_idx = 1:numel(fDATA)
+    PS(:,:,f_idx) = exp(1i*sign_conv*2*pi*fDATA(f_idx)*geomTOFs_error);
+end
+REC_DATA = REC_DATA .* PS; % Phase Screen and Gain Corrections
+clearvars geom_distance col PS geomTOFs_error geomTOFs_disc geomTOFs;
+
+% Which Subset of Transmits to Use
+dwnsmp = 1; % can be 1, 2, or 4 (faster but less quality with more downsampling)
+tx_include = 1:dwnsmp:numElements;
+REC_DATA = REC_DATA(tx_include,:,:); 
+REC_DATA(isnan(REC_DATA)) = 0; % Eliminate Blank Channel
+
+% Extract Subset of Times within Acceptance Angle
+numElemLeftRightExcl = 63;
+elemLeftRightExcl = -numElemLeftRightExcl:numElemLeftRightExcl;
+elemInclude = true(numElements, numElements);
+for tx_element = 1:numElements 
+    elemLeftRightExclCurrent = elemLeftRightExcl + tx_element;
+    elemLeftRightExclCurrent(elemLeftRightExclCurrent<1) = numElements + ...
+         elemLeftRightExclCurrent(elemLeftRightExclCurrent<1);
+    elemLeftRightExclCurrent(elemLeftRightExclCurrent>numElements) = ...
+        elemLeftRightExclCurrent(elemLeftRightExclCurrent>numElements) - numElements;
+    elemInclude(tx_element,elemLeftRightExclCurrent) = false;
+end
+
+% Remove Outliers from Observed Signals Prior to Waveform Inversion
+perc_outliers = 0.99; % Confidence Interval Cutoff
+for f_idx = 1:numel(fDATA)
+    REC_DATA_SINGLE_FREQ = REC_DATA(:,:,f_idx); 
+    signalMagnitudes = elemInclude(tx_include,:).*abs(REC_DATA_SINGLE_FREQ);
+    num_outliers = ceil((1-perc_outliers)*numel(signalMagnitudes));
+    [~,idx_outliers] = maxk(signalMagnitudes(:),num_outliers);
+    REC_DATA_SINGLE_FREQ(idx_outliers) = 0;
+    REC_DATA(:,:,f_idx) = REC_DATA_SINGLE_FREQ;
+end
+
+% Initial Constant Sound Speed Map [m/s]
+c_init = 1480; % Initial Homogeneous Sound Speed [m/s] Guess
+VEL_INIT = c_init*ones(Nyi,Nxi); 
+
+% Initial Constant Attenuation [Np/(Hz m)]
+ATTEN_INIT = 0*alphaNp*ones(Nyi,Nxi);
+
+% (Nonlinear) Conjugate Gradient
+search_dir = zeros(Nyi,Nxi); % Conjugate Gradient Direction
+gradient_img_prev = zeros(Nyi,Nxi); % Previous Gradient Image
+VEL_ESTIM = VEL_INIT; ATTEN_ESTIM = ATTEN_INIT;
+SLOW_ESTIM = 1./VEL_ESTIM + ...
+    1i*sign(sign_conv)*ATTEN_ESTIM/(2*pi); % Initial Slowness Image [s/m]
+crange = [1350, 1600]; % For reconstruction display [m/s]
+attenrange = 10*[-1,1]; % For reconstruction display [dB/(cm MHz)]
+c0 = mean(VEL_ESTIM(:)); cutoff = 0.75; ord = Inf; % Parameters for Ringing Removal Filter
+% Values to Save at Each Iteration
+Niter = sum(niterSoSPerFreq)+sum(niterAttenPerFreq);
+VEL_ESTIM_ITER = zeros(Nyi,Nxi,Niter);
+ATTEN_ESTIM_ITER = zeros(Nyi, Nxi, Niter);
+GRAD_IMG_ITER = zeros(Nyi,Nxi,Niter);
+SEARCH_DIR_ITER = zeros(Nyi,Nxi,Niter);
+for f_idx = 1:numel(fDATA)
+    % Iterations at Each Frequency
+    for iter_f_idx = 1:(niterSoSPerFreq(f_idx)+niterAttenPerFreq(f_idx))
+        tic;
+        % Step 1: Accumulate Backprojection Over Each Element
+        iter = iter_f_idx + sum(niterSoSPerFreq(1:f_idx-1)) + ...
+            sum(niterAttenPerFreq(1:f_idx-1));
+        % Reset CG at Each Frequency (SoS and Attenuation)
+        if ((iter_f_idx == 1) || (iter_f_idx == 1+niterSoSPerFreq(f_idx)))
+            search_dir = zeros(Nyi, Nxi); % Conjugate Gradient Direction
+            gradient_img_prev = zeros(Nyi, Nxi); % Previous Gradient Image
+        end
+        % Attenuation Iterations Always Happen After SoS Iterations
+        if iter_f_idx > niterSoSPerFreq(f_idx)
+            updateAttenuation = true;
+        else
+            updateAttenuation = false; 
+        end
+        gradient_img = zeros(Nyi,Nxi); 
+        % Generate Sources
+        SRC = zeros(Nyi, Nxi, numel(tx_include));
+        for elmt_idx = 1:numel(tx_include)
+            % Single Element Source
+            x_idx_src = x_idx(tx_include(elmt_idx)); 
+            y_idx_src = y_idx(tx_include(elmt_idx)); 
+            SRC(y_idx_src, x_idx_src, elmt_idx) = 1; 
+        end
+        % Forward Solve Helmholtz Equation
+        HS = HelmholtzSolver(xi, yi, VEL_ESTIM, ATTEN_ESTIM, ...
+            fDATA(f_idx), sign_conv, a0, L_PML);
+        [WVFIELD, VIRT_SRC] = HS.solve(SRC,false);
+        % Virtual Sources and Virtual (Sound Speed-Only) Wavefield for Attenuation
+        if updateAttenuation % Virtual Sources for Attenuation
+            % Modify Virtual Source for Attenuation
+            VIRT_SRC = 1i*sign(sign_conv)*VIRT_SRC; 
+        end
+        % Build Adjoint Sources
+        scaling = zeros(numel(tx_include), 1);
+        ADJ_SRC = zeros(Nyi, Nxi, numel(tx_include));
+        for elmt_idx = 1:numel(tx_include)
+            WVFIELD_elmt = WVFIELD(:,:,elmt_idx);
+            REC_SIM = WVFIELD_elmt(ind(elemInclude(tx_include(elmt_idx),:))); 
+            REC = REC_DATA(elmt_idx, ...
+                elemInclude(tx_include(elmt_idx),:), f_idx); REC = REC(:);
+            scaling(elmt_idx) = (REC_SIM(:)'*REC(:)) / ...
+                (REC_SIM(:)'*REC_SIM(:)); % Source Scaling
+            ADJ_SRC_elmt = zeros(Nyi, Nxi);
+            ADJ_SRC_elmt(ind(elemInclude(tx_include(elmt_idx),:))) = ...
+                scaling(elmt_idx)*REC_SIM - REC;
+            ADJ_SRC(:,:,elmt_idx) = ADJ_SRC_elmt;
+        end
+        % Backproject Error
+        ADJ_WVFIELD = HS.solve(ADJ_SRC,true);
+        SCALING = repmat(reshape(scaling, [1,1,numel(scaling)]), [Nyi, Nxi, 1]);
+        BACKPROJ = -real(conj(SCALING.*VIRT_SRC).*ADJ_WVFIELD);
+        % Accumulate Gradient Over Each Element
+        for elmt_idx = 1:numel(tx_include)
+            % Accumulate Backprojection
+            gradient_img = gradient_img + BACKPROJ(:,:,elmt_idx);
+            showAnimation = false; % Show Backprojection Animation
+            if showAnimation
+                % Visualize Numerical Solution
+                imagesc(xi,yi,gradient_img)
+                xlabel('Lateral [m]'); ylabel('Axial [m]'); axis image;
+                title(['Backprojection up to Transmit ', num2str(elmt_idx)]); 
+                colorbar; colormap gray; drawnow;
+                % Source and Receiver Positions Included
+                hold on; plot(x_circ(tx_include(elmt_idx)), ...
+                    y_circ(tx_include(elmt_idx)), 'r*')
+                plot(x_circ(elemInclude(tx_include(elmt_idx),:)), ...
+                    y_circ(elemInclude(tx_include(elmt_idx),:)), 'y.'); 
+                drawnow; clf;
+            end
+        end
+        % Remove Ringing from Gradient Image
+        gradient_img = ringingRemovalFilt(xi, yi, gradient_img, c0, fDATA(f_idx), cutoff, ord);
+        % Step 2: Compute New Conjugate Gradient Search Direction from Gradient
+        % Conjugate Gradient Direction Scaling Factor for Updates
+        if ((iter_f_idx == 1) || (iter_f_idx == 1+niterSoSPerFreq(f_idx)))
+            beta = 0; 
+        else 
+            betaPR = (gradient_img(:)'*...
+                (gradient_img(:)-gradient_img_prev(:))) / ...
+                (gradient_img_prev(:)'*gradient_img_prev(:));
+            betaFR = (gradient_img(:)'*gradient_img(:)) / ...
+                (gradient_img_prev(:)'*gradient_img_prev(:));
+            beta = min(max(betaPR,0),betaFR);
+        end
+        search_dir = beta*search_dir-gradient_img;
+        gradient_img_prev = gradient_img;
+        % Step 3: Compute Forward Projection of Current Search Direction
+        PERTURBED_WVFIELD = HS.solve(VIRT_SRC.*search_dir, false); 
+        dREC_SIM = zeros(numel(tx_include), numElements);
+        for elmt_idx = 1:numel(tx_include)
+            % Forward Projection of Search Direction Image
+            PERTURBED_WVFIELD_elmt = PERTURBED_WVFIELD(:,:,elmt_idx);
+            dREC_SIM(elmt_idx,elemInclude(tx_include(elmt_idx),:)) = -permute(scaling(elmt_idx) * ...
+                PERTURBED_WVFIELD_elmt(ind(elemInclude(tx_include(elmt_idx),:))),[2,1]);
+        end
+        % Step 4: Perform a Linear Approximation of Exact Line Search
+        perc_step_size = 1; % (<1/2) Introduced to Improve Compliance with Strong Wolfe Conditions 
+        alpha = -(gradient_img(:)'*search_dir(:))/(dREC_SIM(:)'*dREC_SIM(:));
+        if updateAttenuation % Update Complex Slowness
+            SI = sign(sign_conv) * imag(SLOW_ESTIM) + perc_step_size * alpha * search_dir;
+            SLOW_ESTIM = real(SLOW_ESTIM) + 1i * sign(sign_conv) * SI;
+        else
+            SLOW_ESTIM = SLOW_ESTIM + perc_step_size * alpha * search_dir;
+        end 
+        VEL_ESTIM = 1./real(SLOW_ESTIM); % Wave Velocity Estimate [m/s]
+        ATTEN_ESTIM = 2*pi*imag(SLOW_ESTIM)*sign(sign_conv);
+        % Save Intermediate Results
+        VEL_ESTIM_ITER(:,:,iter) = VEL_ESTIM;
+        ATTEN_ESTIM_ITER(:,:,iter) = ATTEN_ESTIM;
+        GRAD_IMG_ITER(:,:,iter) = gradient_img;
+        SEARCH_DIR_ITER(:,:,iter) = search_dir;
+        % Visualize Numerical Solution
+        subplot(2,2,1); imagesc(xi,yi,VEL_ESTIM,crange);
+        title(['Estimated Wave Velocity ', num2str(iter)]); axis image;
+        xlabel('Lateral [m]'); ylabel('Axial [m]'); colorbar; colormap gray;
+        subplot(2,2,2); imagesc(xi,yi,Np2dB*slow2atten*ATTEN_ESTIM,attenrange);
+        title(['Estimated Attenuation ', num2str(iter)]); axis image;
+        xlabel('Lateral [m]'); ylabel('Axial [m]'); colorbar; colormap gray;
+        subplot(2,2,3); imagesc(xi,yi,search_dir)
+        xlabel('Lateral [m]'); ylabel('Axial [m]'); axis image;
+        title(['Search Direction Iteration ', num2str(iter)]); colorbar; colormap gray; 
+        subplot(2,2,4); imagesc(xi,yi,-gradient_img)
+        xlabel('Lateral [m]'); ylabel('Axial [m]'); axis image;
+        title(['Gradient Iteration ', num2str(iter)]); colorbar; colormap gray; 
+        drawnow; disp(['Iteration ', num2str(iter)]); toc;
+    end
+end
+
+% Plot Final Reconstructions
+subplot(1,2,1); imagesc(xi,yi,VEL_ESTIM,crange);
+title(['Estimated Wave Velocity ', num2str(iter)]); axis image;
+xlabel('Lateral [m]'); ylabel('Axial [m]'); colorbar; colormap gray;
+subplot(1,2,2); imagesc(xi,yi,Np2dB*slow2atten*ATTEN_ESTIM,attenrange); 
+title(['Estimated Attenuation ', num2str(iter)]); axis image;
+xlabel('Lateral [m]'); ylabel('Axial [m]'); colorbar; colormap gray;
+
+% Save the Result to File
+filename_results = ['Results/', filename, '_WaveformInversionResults.mat'];
+save(filename_results, '-v7.3', 'xi', 'yi', 'fDATA', 'niterAttenPerFreq', ...
+    'niterSoSPerFreq', 'VEL_ESTIM_ITER', 'ATTEN_ESTIM_ITER', 'GRAD_IMG_ITER', 'SEARCH_DIR_ITER')