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.

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

SDA forward solvers need check on optical properties to make sure solution is valid #88

Closed
hayakawa16 opened this issue Aug 22, 2023 · 7 comments · Fixed by #160
Closed
Assignees

Comments

@hayakawa16
Copy link
Member

SDA solutions are only valid when mus'>>mua. This check is not currently performed in these solvers and lead to erroneous results.

The forward solvers that require this check are:
DistributedGaussianSourceSDAForwardSolver, DistributedPointSourceSDAForwardSolver, PointSourceSDAForwardSolver, TwoLayerSDAForwardSolver and possibly SFDDiffusionForwardSolver.

To reproduce this error define mua=1, mus'=1 as optical properties in TwoLayerSDAForwardSolver and view FluenceOfRhoAndZ method results.

Code should be added to throw an ArgumentException after checking optical properties for case of mua==musp. This is the extreme lower bound for mus'/mua. It is difficult to know a higher lower bound that would invalidate each SDA solution.

@hayakawa16 hayakawa16 self-assigned this Aug 22, 2023
@dcuccia
Copy link
Contributor

dcuccia commented Aug 22, 2023

Thanks Carole. I think there's an anomaly at exactly (1, 1) for upper layer properties in the TwoLayer solver, not necessarily in the general case.

@hayakawa16
Copy link
Member Author

Sorry! I posted on other issue and should be posting here.
I return to the purple book section 6.5.1: In previous sections, it was stated but not proven that mua << mus(1-g) is the criterion for the validity of the diffusion approximation....[other criteria]....Finally, since the fluence rate is the most important measurable quantity, one could also require that Eq.(6.15) holds: fluence rate phi_s(r)=4pi L0 where L0 is radiance for perfectly isotropic distribution.
Envoking in code all of these criteria could get tricky because far enough from the source, the light field is isotropic for all OPs. When that happens is difficult to know apriori. And even if I check the OPs, the SDA fails close to the source (since field is not isotropic yet).
Part of me says let users use the code erroneously if they choose to, part of me says they should be alerted.

@dcuccia
Copy link
Contributor

dcuccia commented Aug 22, 2023

Perhaps putting out a one-time warning in the validation check. But I'd err strongly on the side of letting the computation proceed. Hard to know how someone wants to use it. For example, they could use it to show how wrong the SDA is, compared to a more accurate model.

@hayakawa16
Copy link
Member Author

Like your feedback. A similar philosophy could be taken with Issue #13 and we can work on that too.

@scottprahl
Copy link

@hayakawa16 Just noticed this. The diffusion approximation is exact when mus=0. The makes the mua << mus' limitation marginally useful. IIRC the diffusion approximation is worst when mua≈mus' (which just happens to be the case you started this issue with :)

In any case, I agree with @dcuccia that there are use cases for allowing any optical properties.

@hayakawa16
Copy link
Member Author

Thanks for your input @scottprahl! I will create a branch and work on adding a warning but letting the code proceed.

@hayakawa16
Copy link
Member Author

The way that I fixed this is that I added a check in DiffusionParameters:
if (op.Mua >= op.Musp) Console.WriteLine("Warning Mua >= Musp");

This warning is shown in the MATLAB interop output window if a user specifies Mua >= Musp. The warning is output and the code proceeds.

This warning is not seen in the text window of the GUI yet. That fix requires a bit more work.

I welcome any comments or concerns.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging a pull request may close this issue.

3 participants