Monte Carlo Basic produces erroneous reflectance

[ccampb19]

I'm trying to use VTS to generate some reflectance data in Matlab, but forward-modeled reflectance is way off when using the basic MC solver. Have I made an incorrect assumption?

To reproduce: VTS_MATLAB_v4.3.0Beta Matlab R2017b:

clear
startup

% Generate spectra
rho = 22;
wv = 600:1000;
absorbers.Names =           {'HbO2', 'Hb', 'H2O', 'Fat', 'Melanin'};
absorbers.Concentrations = [ 28.4, 22.4, 0.7, 0, 0.0051 ];
scatterer.Type = 'Intralipid';
scatterer.VolumeFraction = 0.01;
op = VtsSpectroscopy.GetOP(absorbers, scatterer, wv);

% Forward model reflectance with three models
VtsSolvers.SetSolverType('PointSourceSDA');
R_sda = VtsSolvers.ROfRho(op, rho);

VtsSolvers.SetSolverType('MonteCarlo');
R_basic = VtsSolvers.ROfRho(op, rho);

VtsSolvers.SetSolverType('Nurbs');
R_nurbs = VtsSolvers.ROfRho(op, rho);

% plot (log scale to keep all spectra visible)
figure
subplot(1,2,1)
plot(wv,op(:,1))
xlabel('Wavelength (nm)')
title('\mu_A (mm^{-1}')

subplot(1,2,2)
plot(wv,(op(:,2)))
xlabel('Wavelength (nm)')
title('\mu_S'' (mm^{-1}')

figure;semilogy(wv,R_sda,wv,R_basic,wv,R_nurbs)
xlabel('Wavelength')
title('Forward Reflectance')
legend('SDA','MC Basic','Nurbs','Location','Best')

[hayakawa16]

This is a very good catch! I am able to replicate the result. The error occurs in the Vts library and is independent of the Matlab interop code. I can reproduce the error using our WPF GUI Forward Solver panel and checking "use spectral panel inputs" to obtain R vs wavelength results using Basic MC and comparing with any other solver. Important difference: R(rho=10) mua=0.0666/mm mus'=1.2/mm produces a different result than R(rho=10,lambda=1000nm) using skin OPs mua=0.666/mm mus'=1.2/mm.

Okay this took me a little bit to debug, but I got it!! I'm providing details primarily for my fellow developers so please bear with me. The problem was that in the Basic MC code for R(rho) there are two loops, the outer loop loops over OPs and the inner loop loops over rhos. The array RatRhoMC was initialized outside of the outer loop, however it is created by a cumulative sum so if there were multiple OPs (or wavelengths), RatRhoMC was just growing and growing producing erroneous results for all OPs greater than the first. To fix this, I added an "Array.Clear" of RatRhoMC inside the outer loop. This problem also occurred in the code for ROfFx.

I have pushed these fixes to Github. It might take a little bit before we have a new Matlab interop download for you. I'll post when we do.

Thanks for finding this bug and helping us improve our code base! Best, Carole

[hayakawa16]

ccampb19 has also pointed out that another forward solver does not work across wavelength.
1) Distributed Gaussian Source with beam diameter=0, check "use spectral panel inputs", select R(rho,ft), check "allow multi-axis selection", check lambda, leave rho and ft unchecked with rho=1mm ft=0GHz crashes the GUI.

I tried other solution domains and found two other setups that crash the GUI.
2) Distributed Gaussian Source with beam diameter=0, check "use spectral panel inputs", select R(fx,t), check "allow multi-axis selection", check lambda, leave fx and t unchecked with fx=0 and t=0.05ns crashes. 3) Distributed Gaussian Source with beam diameter=0, check "use spectral panel inputs", select R(rho,t), check "allow multi-axis selection", check lambda, leave rho and t unchecked with rho=1 and t=0.05ns crashes.

I'd like to fix these problems before we close this issue. Carole

[hayakawa16]

I determined that the reason the GUI crashes for these selections is that time-dependent solutions for the DistributedGaussianSourceSDA solver are not implemented. I modified the GUI code so that when this forward solver is selected all time-dependent solution domain options (e.g. R(rho,time), R(rho,ft), R(fx,time), R(fx,ft)) are greyed out.

[dcuccia]

Thx Carole. We should have a generic way of greying out non-implemented domains.

On Thu, Sep 20, 2018, 6:42 PM Carole Hayakawa notifications@github.com wrote:

I determined that the reason the GUI crashes for these selections is that time-dependent solutions for the DistributedGaussianSourceSDA solver are not implemented. I modified the GUI code so that when this forward solver is selected all time-dependent solution domain options (e.g. R(rho,time), R(rho,ft), R(fx,time), R(fx,ft)) are greyed out.

— You are receiving this because you are subscribed to this thread. Reply to this email directly, view it on GitHub https://github.com/VirtualPhotonics/VTS/issues/12#issuecomment-423386996, or mute the thread https://github.com/notifications/unsubscribe-auth/AAdRgbHRCeISOo_x_R2r6imO-kQCt08oks5udEQPgaJpZM4Wnh-V .

[hayakawa16]

Hi, I just wanted to close the loop on these issues.
1) In the latest Matlab download, version 4.5 Beta, this scaled Monte Carlo now produces improved results. 2) In the latest WPF GUI release, WPF Version 2.3 Beta, on the Forward panel when the user selects distributed Gaussian source, the Solution Domain options are R(rho) and R(fx), all others are greyed out. Carole