Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Bug/173 monte carlo inverse code fixes #174

Merged
merged 6 commits into from
Sep 26, 2024
Merged

Bug/173 monte carlo inverse code fixes #174

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
51c73be
Found bug using wvs instead of wv, modified all to wvs since multiple…
hayakawa16 Sep 20, 2024
b862f01
@janakarana helped me to create SpectralDictionary.xlsx with the newl…
hayakawa16 Sep 21, 2024
7f2314c
Modified get_mua to use readtable instead of xlsread which is obsolet…
hayakawa16 Sep 23, 2024
a26f77e
Modified readtable to readmatrix because readmatrix takes away headers.
hayakawa16 Sep 24, 2024
9fcbce3
Found an error in the setting of the measParms in Example 1. Did not…
hayakawa16 Sep 24, 2024
dfd209d
Missed updating Linux comments that I corrected in mc_inverse_Win.m
hayakawa16 Sep 24, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Binary file modified matlab/monte_carlo_inverse/general/SpectralDictionary.xlsx
View file
Open in desktop
Binary file not shown.
35 changes: 24 additions & 11 deletions matlab/monte_carlo_inverse/general/get_mua.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,26 +2,35 @@
function [mua,dmua]=get_mua(absorbers,wavelengths)
persistent data;
if isempty(data)
data=xlsread('SpectralDictionary.xlsx'); % file has 1 tab only
data=readmatrix('SpectralDictionary.xlsx','Range','A2:J1202'); % file has 1 tab only
end
% the following assumes all columns in table have same length
% read Fat data
wv_fat=data(:,1);
spec_fat=data(:,2);
% read H2O data
last_line_h2o=209;
wv_h2o=data(2:last_line_h2o,1);
spec_h2o=data(2:last_line_h2o,2);
wv_h2o=data(:,3);
spec_h2o=data(:,4);
% read Hb data
last_line_hb=376;
wv_hb=data(2:last_line_hb,3);
spec_hb=data(2:last_line_hb,4);
wv_hb=data(:,5);
spec_hb=data(:,6);
% read HbO2 data
last_line_hbo2=376;
wv_hbo2=data(2:last_line_hbo2,5);
spec_hbo2=data(2:last_line_hbo2,6);
wv_hbo2=data(:,7);
spec_hbo2=data(:,8);
% read melanin data
wv_melanin=data(:,9);
spec_melanin=data(:,10);

mua=zeros(1,length(wavelengths));
dmua=zeros(length(wavelengths),length(absorbers.Names));
values=zeros(1,length(wavelengths));
for i=1:length(absorbers.Names)
% determine if spectra is in molar absorption coefficients or fractional abs coefficient
if (strcmp(absorbers.Names(i),'H2O')) % fractional
if (strcmp(absorbers.Names(i),'Fat')) % fractional
values=interp1(wv_fat,spec_fat,wavelengths);
mua=mua+absorbers.Concentrations(i)*values;
dmua(:,i)=dmua(:,i)+values';
elseif (strcmp(absorbers.Names(i),'H2O')) % fractional
values=interp1(wv_h2o,spec_h2o,wavelengths);
mua=mua+absorbers.Concentrations(i)*values;
dmua(:,i)=dmua(:,i)+values';
Expand All @@ -33,6 +42,10 @@
values=interp1(wv_hbo2,spec_hbo2,wavelengths);
mua=mua+log(10)*absorbers.Concentrations(i)*values;
dmua(:,i)=dmua(:,i)+log(10)*values';
elseif (strcmp(absorbers.Names(i),'Melanin')) % fractional
values=interp1(wv_melanin,spec_melanin,wavelengths);
mua=mua+absorbers.Concentrations(i)*values;
dmua(:,i)=dmua(:,i)+values';
end
end
end
57 changes: 33 additions & 24 deletions matlab/monte_carlo_inverse/linux_and_mac/mc_inverse.m
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,14 +2,18 @@
% This script includes three examples:
% 1) Example 1: equivalent to Example 7 in vts_mc_demo.m: use N=1000 in
% infile_pMC_db_gen_template.txt to reproduce results in vts_mc_demo.m
% for multiple rho
% 2) Example 2: inverse solution for chromophore concentrations
% for multiple wavelengths, single rho: use
% wv=500:100:1000 and rho=1, N=10000 in infile_pMC_db_gen_template.txt
% for multiple wavelengths, single rho: use N=10000 in
% infile_pMC_db_gen_template.txt. Example specifies
% wvs=500:100:1000nm and rho=1mm,
% these selections different than those used in vts_solver_demo.m
% 3) Example 3: inverse solution for chromophore concentrations
% for multiple wavelengths, two rho: use
% for multiple wavelengths, two rho: use N=10000 in
% infile_pMC_db_gen_template.txt. Example specifies
% wvs=600:100:900nm and rho=0.1429mm and 1mm.
%
% but does not require *MATLAB interop* code to run
% Note! This script does not require *MATLAB interop* code to run
%%
clear all
clc
Expand All @@ -25,7 +29,7 @@
% 4) normalization of chi2
% 5) optimset options selected
%% read in baseline OPs from database gen infile
x0 = [0.01, 5.0]; % baseline values for database and initial guess [mua, mus]
x0 = [0.01, 5.0]; % baseline values for database and initial guess [mua, mus], i.e. musp=1.0
g = 0.8;
% input rho: in inversion don't use last point since includes tallies from
% beyond
Expand All @@ -51,7 +55,7 @@
lb=[]; ub=[];

% input measData
measParms = [ 0.04 0.95 ]; % mua, musp
measParms = [ 0.04 4.75 ]; % mua, mus, i.e. musp=0.95
measData = [0.0331 0.0099 0.0042 0.0020 0.0010];

% option: divide measured data and forward model by measured data
Expand Down Expand Up @@ -86,9 +90,9 @@
legend('Meas','IG','Converged','Location','Best');
legend boxoff;
grid on;
set(gca,'TickDir','out','FontSize',20);
title('Inverse solution using pMC/dMC');
set(f, 'Name', 'Inverse solution using pMC/dMC');
set(gca,'TickDir','out','FontSize',12);
title('Ex 1: Inverse solution using pMC/dMC');
set(f, 'Name', 'Ex 1: Inverse solution using pMC/dMC');
disp(sprintf('Meas = [%f %5.3f]',measParms(1),measParms(2)));
disp(sprintf('IG = [%f %5.3f] Chi2=%5.3e',x0(1),x0(2),...
(measData-R_ig)*(measData-R_ig)'));
Expand All @@ -104,7 +108,9 @@
% Run a Monte Carlo simulation with pMC post-processing enabled
% Use generated database to solve inverse problem with measured data
% generated using Nurbs and selected concentrations
% define 8 rho, however only 4th used in inverse problem

% define 8 rho, however only 4th bin is used in inverse problem
% following rho specifications needed to fill out infile template
rhostart=0;
rhostop=2;
rhocount=8;
Expand All @@ -125,7 +131,8 @@
g=0.8;
n=1.4;

% ops has dimensions [numwv 4] NOTE: if using Octave type "pkg load io" in Command Window
% ops has dimensions [numwvs 4] 4=[mua,musp,g,n]
% NOTE: if using Octave type "pkg load io" in Command Window
% to load IO package and read SpectralDictionary.xlsx in next command
[ops,dmua,dmusp]=get_optical_properties(absorbers,scatterers,wvs);

Expand Down Expand Up @@ -154,10 +161,11 @@
lb=[]; ub=[];
% input measData taken from vts_solver_demo using Nurbs rho=1mm and
% concentrations
%measParms = [ 70, 30, 0.8 ]; % debug with same as measured
%measParms = [ 70, 30, 0.8 ]; % debug with same as initial guess
%measData = [0.0089 0.0221 0.0346 0.0301 0.0251 0.0198];
measParms = [ 72, 35, 0.6 ];
measData = [0.0082 0.0208 0.0342 0.0299 0.0250 0.0205];

% run lsqcurvefit if have Optimization Toolbox because it makes use of
% dMC differential Monte Carlo predictions
% if don't have Optimization Toolbox, run non-gradient, non-constrained
Expand Down Expand Up @@ -185,13 +193,13 @@
wvs,R_ig,'g-',...
wvs,R_conv,'b:','LineWidth',2);
xlabel('\lambda [nm]');
ylabel('log10[R(\lambda)]');
ylabel('R(\lambda)');
legend('Meas','IG','Converged','Location','Best');
legend boxoff;
grid on;
set(gca,'TickDir','out','FontSize',20);
%title('Inverse solution using pMC/dMC');
set(f, 'Name', 'Inverse solution using pMC/dMC');
set(gca,'TickDir','out','FontSize',12);
title('Ex 2: Inverse solution using pMC/dMC');
set(f, 'Name', 'Ex 2: Inverse solution using pMC/dMC');
disp(sprintf('Meas = [%5.3f %5.3f %5.3f]',measParms(1),measParms(2),measParms(3)));
disp(sprintf('IG = [%5.3f %5.3f %5.3f] Chi2=%5.3e',igConc(1),igConc(2),...
igConc(3),(measData-R_ig)*(measData-R_ig)'));
Expand All @@ -201,9 +209,9 @@
abs(measParms(2)-recoveredOPs(2))/measParms(2),abs(measParms(3)-recoveredOPs(3))/measParms(3)));
%% ======================================================================= %
% Example 3: Inverse solution for R(rho,wavelength) chromophore
% concentrations (Hb, HbO2) and scatterer coefficients (a,b) for
% concentrations (Hb, HbO2) and scatterer coefficients (a,b) for
% 4 wavelengths and 2 rho
% define 8 rho, however only 1st and 4th used in inverse problem
% define 8 rho, however only 1st and 4th bins used in inverse problem
rhostart=0;
rhostop=2;
rhocount=8;
Expand All @@ -226,7 +234,8 @@
n=1.4;
igParms = [absorbers.Concentrations(1) absorbers.Concentrations(2) scatterers.Coefficients(1) scatterers.Coefficients(2)];

% ops has dimensions [numwv 4] NOTE: if using Octave type "pkg load io" in Command Window
% ops has dimensions [numwv 4] 4=[mua,musp,g,n]
% NOTE: if using Octave type "pkg load io" in Command Window
% to load IO package and read SpectralDictionary.xlsx in next command
[ops,dmua,dmusp]=get_optical_properties(absorbers,scatterers,wvs);

Expand Down Expand Up @@ -255,7 +264,7 @@

%% use unconstrained optimization lb=[-inf -inf]; ub=[inf inf];
lb=[]; ub=[];
% input measData taken from vts_solver_demo using Nurbs rhoMidpoint=0.1429,1.0mm and
% input measData taken from vts_solver_demo using Nurbs rhoMidpoints=0.1429,1.0mm and
% parameters { 28.4, 22.4, 1.2, 1.42 }
measParms = [28.4, 22.4, 1.2, 1.42 ];
measData = [0.2931 0.2548 0.2068 0.1711 0.0255 0.0355 0.0320 0.0277];
Expand Down Expand Up @@ -287,14 +296,14 @@
wvs,R_ig(1:length(wvs)),'g-',wvs,R_ig(length(wvs)+1:end),'g--',...
wvs,R_conv(1:length(wvs)),'b-',wvs,R_conv(length(wvs)+1:end),'b--','LineWidth',2);
xlabel('\lambda [nm]');
ylabel('log10[R(\lambda)]');
ylabel('R(\lambda)');
legend('Meas rho=0.1429','Meas rho=1.0','IG rho=0.1429','IG rho=1.0',...
'Converged rho=0.1429','Converged rho=1.0','Location','Best');
legend boxoff;
grid on;
set(gca,'TickDir','out','FontSize',20);
title('Inverse solution using pMC/dMC');
set(f, 'Name', 'Inverse solution using pMC/dMC');
set(gca,'TickDir','out','FontSize',12);
title('Ex 3: Inverse solution using pMC/dMC');
set(f, 'Name', 'Ex 3: Inverse solution using pMC/dMC');
disp(sprintf('Meas = [%5.3f %5.3f %5.3f %5.3f]',measParms(1),measParms(2),...
measParms(3),measParms(4)));
disp(sprintf('IG = [%5.3f %5.3f %5.3f %5.3f] Chi2=%5.3e',igParms(1),igParms(2),...
Expand Down
Loading
Loading