Skip to content

Commit

Permalink
Merge pull request #174 from VirtualPhotonics/bug/173-monte-carlo-inv…
Browse files Browse the repository at this point in the history 
…erse-code-fixes

Bug/173 monte carlo inverse code fixes
  • Loading branch information
@hayakawa16
hayakawa16 authored Sep 26, 2024
2 parents 66cd90a + dfd209d commit d7e80eb
Show file tree
Hide file tree
Showing 8 changed files with 137 additions and 94 deletions.
Binary file modified matlab/monte_carlo_inverse/general/SpectralDictionary.xlsx
View file
Binary file not shown.
35 changes: 24 additions & 11 deletions matlab/monte_carlo_inverse/general/get_mua.m
View file
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
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

0 comments on commit d7e80eb

Please sign in to comment.