phase function implementation with lookup table

[Anouk Post]

Dear Dr. Fang,

I am working on simulations where I would like to be able to use any phase function I could come up with. This can be implemented with a look-up table with the approach from this paper by Peter Naglic. et al. I have implemented this myself in MCML. However, since MCXLab is considerably faster I would like to use MCXLab instead.

I could try to implement it myself, but before I do that I would just like to check - this has not yet been implemented, has it? Or is someone currently already working on it?

Thank you in advance,

Anouk Post

[Qianqian Fang]

hi Anouk,

we had this plan (same method, i.e. user-specified inv(CDF) as a lookup table) since 2015 in this TODO, however, due to other development proprieties, we still haven't implemented this feature yet.

https://github.com/fangq/mcx/issues/13

if you've already done this or know how to do this quickly, we are happy to accept/include your patch once you finish the implementation and testing. My recommendation is to place this look-up table in the constant memory or, even faster, shared memory, because one needs to read it every scattering event. I was also planning to accept a spline/polynomial coefficients instead of a look-up table to greatly reduce the look-up time. Either look-up table or the spline coefficients should be read from the JSON input interface under the "Domain" section, as a 1D vector, something like cfg.phasefun.

otherwise, we will start working on this feature in the development cycle after the v2021.2 release is published.

https://github.com/fangq/mcx#whats-new

Qianqian

Thank you in advance,

Anouk Post

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/93c67727-809a-4036-9ac2-58f1daf5923an%40googlegroups.com.

[Anouk Post]

Dear Dr. Fang,

Thank you for your quick reply!

Since I am currently working on a paper where I need quite some MC simulations (tens of thousands), I am currently trying to make the modifications myself. 

Just to give you an idea, I have copied the code below that I use in my current MC software that is written in C++. The results of this code have been validated with regular MCML. I already have Matlab code that creates a phase function look-up table, that could be stored as cfg.phasefun (which I am happy to share). So all I would need to change in the MCX code is to read in cfg.phasefun and to sample the angle from this. cfg.phasefun would then be an inverse look-up table based on the cumulative distribution function, and it would be a vector with a number of cos(theta) values equal to Ntable. This list of cos(theta) values then needs to be read by the MCX code and used to determine the angle. 

Within my C-code I then have the following few lines to determine the angle of a scattering event. Here, mui is the list of cos(theta) values in cfg.phasefun, Ntable is the length of cfg.phasefun. Finally, the cosine of the scattering angle costheta and the sine of the scattering angle sintheta are calculated and later on used to determine the direction of the photon. If your code uses the angle itself that could be easily modified. 

rnd = RandomGen(1, 0, NULL);

double ind = rnd*(Ntable - 1); /* find index in table (decimal number)*/

int low = floor(ind); /*find lower boundary of table index*/

int high = ceil(ind); /* find higher boundary of table index */

/* linear interpolation*/

double diffcos = mui[high] - mui[low];

double mu_low = mui[low];

double frac = ind - low;

costheta = diffcos*frac + mu_low;

sintheta = sqrt(1.0 - costheta*costheta);   /*sqrt faster than sin()*/

So I do know which changes to make in the C++-code, but I am by no means an experienced programmer in C++ (I mainly work in Matlab). Unfortunately, I have been struggling all afternoon to try to compile the source code, so that I can change the source code as a next step to implement the phase function look-up table. I am working on a Windows computer. I have tried the visual studio approach, MSYS, and the command line to compile but I keep getting errors. I have looked online for solutions to those errors and they seem to work, but I also keep getting new ones. 

So I am currently a bit stuck. Could you give me some tips to compile the software with Windows? Which approach is most likely to succeed?

And since there are quite a few source code files, could you point me to where in the source code I should make the modifications I outlined above?

I am very grateful that you developed MCX and I understand that you must be very busy with the new release, so any help is greatly appreciated. 

Yours sincerely,

Anouk

Op maandag 22 februari 2021 om 14:43:39 UTC+1 schreef fan...@gmail.com:

[Fang, Qianqian]

On 2/22/21 3:56 PM, Anouk Post wrote:

Dear Dr. Fang,

Thank you for your quick reply!

Since I am currently working on a paper where I need quite some MC simulations (tens of thousands), I am currently trying to make the modifications myself. 

Just to give you an idea, I have copied the code below that I use in my current MC software that is written in C++. The results of this code have been validated with regular MCML. I already have Matlab code that creates a phase function look-up table, that could be stored as cfg.phasefun (which I am happy to share). So all I would need to change in the MCX code is to read in cfg.phasefun and to sample the angle from this. cfg.phasefun would then be an inverse look-up table based on the cumulative distribution function, and it would be a vector with a number of cos(theta) values equal to Ntable. This list of cos(theta) values then needs to be read by the MCX code and used to determine the angle. 

Within my C-code I then have the following few lines to determine the angle of a scattering event. Here, mui is the list of cos(theta) values in cfg.phasefun, Ntable is the length of cfg.phasefun. Finally, the cosine of the scattering angle costheta and the sine of the scattering angle sintheta are calculated and later on used to determine the direction of the photon. If your code uses the angle itself that could be easily modified. 

rnd = RandomGen(1, 0, NULL);

double ind = rnd*(Ntable - 1); /* find index in table (decimal number)*/

int low = floor(ind); /*find lower boundary of table index*/

int high = ceil(ind); /* find higher boundary of table index */

/* linear interpolation*/

double diffcos = mui[high] - mui[low];

double mu_low = mui[low];

double frac = ind - low;

costheta = diffcos*frac + mu_low;

sintheta = sqrt(1.0 - costheta*costheta);   /*sqrt faster than sin()*/

So I do know which changes to make in the C++-code, but I am by no means an experienced programmer in C++ (I mainly work in Matlab). Unfortunately, I have been struggling all afternoon to try to compile the source code, so that I can change the source code as a next step to implement the phase function look-up table. I am working on a Windows computer. I have tried the visual studio approach, MSYS, and the command line to compile but I keep getting errors. I have looked online for solutions to those errors and they seem to work, but I also keep getting new ones. 

So I am currently a bit stuck. Could you give me some tips to compile the software with Windows? Which approach is most likely to succeed?

building stuff on Windows is a total nightmare for deploying cross-platform software ... if you have access to a Linux box, building mcx/mcxlab on a Linux machine is orders of magnitude easier - just install a few packages via apt-get (sudo apt-get install gcc nvidia-cuda-toolkit cmake), the Makefile in mcx should work out of box. Even building mcx/mcxlab on a lightweight CI linux container is not that difficult

https://github.com/fangq/mcx/blob/master/.travis.yml#L25-L47

however, if you really have to build mcx on Windows, you need to install

- cygwin64 or MSYS2 (select gcc, make during setup)
- cuda toolkit (nvcc)
- visual studio (cl.exe)
- zlib (to provide zlib.lib and zlib.dll), or use zlib-win-build at https://github.com/kiyolee/zlib-win-build

the biggest issue is that nvcc on windows only supports VS's compiler cl (see this forum thread), and does not support gcc. This makes things dramatically more difficult. Even linking with the most commonly used compression library, zlib, is quite time-consuming (I still haven't figured out how to do static linking with zlibstat ...)

anyways, to get it to work on Windows, you need to install cuda and VS and add nvcc/cl's paths to your env PATH variable. You need to then install cygwin64 or MSYS2 with make/gcc packages (we can't use gcc, but the makefile still needs it). Once this is done, open a cygwin/MSYS2 terminal, you should be able to type nvcc and cl to see version info.

then, you need to download zlib source package from https://zlib.net and build zlib.lib/zlib.dll, or use this one to build https://github.com/kiyolee/zlib-win-build, once this is done, you should copy zlib.lib/zlib.dll to mcx/src and include/{zlib.h,zconf.h} to VS's include folder (C:\Program Files (x86)\Microsoft Visual Studio ??.0\VC\include).

Then, you should be able to cd mcx/src, and type

make static

to generate mcx.exe under mcx/bin. Setting up matlab/mex to build mcxlab requires a lot more steps, so I will skip it for now.

And since there are quite a few source code files, could you point me to where in the source code I should make the modifications I outlined above?

you need to at least modify

- mcx_utils.h/mcx_utils.c to add the cfg subfield phasefun as a float pointer and and phasefunlen as an integer for the length of the buffer, then read in both from a JSON file in the mcx_loadjson() function
- mcx_core.cu: this is the GPU host code and kernel, you need to change a few things,

    - in mcx_run_simulation() host code, you can either pass on the host pointer cfg.phasefun to a constant mem buffer, similar to gproperty, or first to a global memory, then copy to the shared member inside the gpu kernel, and
    - in mcx_main_loop() GPU kernel, do all things necessary to use this buffer for scattering angle calculation, especially here, https://github.com/fangq/mcx/blob/master/src/mcx_core.cu#L1468-L1483, if you use a look-up table, you need to implement a binary search

https://github.com/fangq/mcx/blob/master/src/mcx_core.cu#L2038-L2049
https://github.com/fangq/mcx/blob/master/src/mcx_core.cu#L1339-L1366

the source code was documented in details, but feel free to let me know if you have any further question.

Qianqian

I am very grateful that you developed MCX and I understand that you must be very busy with the new release, so any help is greatly appreciated. 

Yours sincerely,

Anouk

Op maandag 22 februari 2021 om 14:43:39 UTC+1 schreef fan...@gmail.com:

On 2/22/21 7:49 AM, Anouk Post wrote:

Dear Dr. Fang,

I am working on simulations where I would like to be able to use any phase function I could come up with. This can be implemented with a look-up table with the approach from this paper by Peter Naglic. et al. I have implemented this myself in MCML. However, since MCXLab is considerably faster I would like to use MCXLab instead.

I could try to implement it myself, but before I do that I would just like to check - this has not yet been implemented, has it? Or is someone currently already working on it?

hi Anouk,

we had this plan (same method, i.e. user-specified inv(CDF) as a lookup table) since 2015 in this TODO, however, due to other development proprieties, we still haven't implemented this feature yet.

https://github.com/fangq/mcx/issues/13

if you've already done this or know how to do this quickly, we are happy to accept/include your patch once you finish the implementation and testing. My recommendation is to place this look-up table in the constant memory or, even faster, shared memory, because one needs to read it every scattering event. I was also planning to accept a spline/polynomial coefficients instead of a look-up table to greatly reduce the look-up time. Either look-up table or the spline coefficients should be read from the JSON input interface under the "Domain" section, as a 1D vector, something like cfg.phasefun.

otherwise, we will start working on this feature in the development cycle after the v2021.2 release is published.

https://github.com/fangq/mcx#whats-new

Qianqian

Thank you in advance,

Anouk Post

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/93c67727-809a-4036-9ac2-58f1daf5923an%40googlegroups.com.

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/5f7b6522-7ab8-49ab-854b-bf4605bbbcb9n%40googlegroups.com.

[Maxime Baillot]

Just to add my contribution to the discussion.

I worked on doing that (and sorry Dr Fang for not giving you any feedback about that)

Having a look up table is really hard because of multiple problem.

1. You need a number of point big enough so the error is not too big
2. Limitation of the GPU memory

The only solution that I found would be to use the texture memory of the GPU because it is not use in MCX and it do an "automatic" interpolation.

I didn't implemented it because it would take a lot of time and a big modification of the code. I don't have the time and enough knowledge about MCX source code.

The solution that I found, and use, is to modify directly the phase function and did other simple modifications. (I didn't submit any modification on github because it's messy and not really flexible. Dr Fang if you are interested I can send you the source code to take a look at)


If it help you, and you are still trying to do something, the main modification is changing that :

//__constant__ float4 gproperty[MAX_PROP_AND_DETECTORS];
__constant__ Medium property[MAX_PROP];

__constant__ float4 detectors[MAX_DETECTORS];

[Anouk Post]

Dear Maxime,

Thank you!

I originally asked this question, but I haven't had time yet to modify the MCX code myself, since I have some problems getting access to a Linux system. Could you explain a bit more what you modified? Because I thought that MCX used a cumulative distribution function of the phase function to sample the angles, which has an analytical solution for the Henyey-Greenstein. But for other phase functions, there is (as far as I know) no formula that describes their cumulative distribution functions. 

Regarding what you said about look-up tables - a look-up table should be possible for most phase function types. For a Mie phase function you might indeed need a very large look-up table for it to be accurate, but for the Modified Henyey Greenstein, the double Henyey Greenstein, and the Reynolds-McCormick phase functions it should work. I currently have implemented that myself in a modification of MCML. However, that code results in slow simulations, and I can only make two layers, nothing else. That is why I am so interested in MCX. 

Best wishes,

Anouk

You received this message because you are subscribed to a topic in the Google Groups "mcx-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mcx-users/pckAVbf8z7U/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mcx-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/3f2f42c9-cd95-485c-8773-f51c3bc360d4n%40googlegroups.com.

[Maxime Baillot]

There is,in general, formula for the CDF  but no formula for the Inverse of CDF (as far as we know or maybe we just need to read russian (every one know that the russians did everything in maths! hahaha))

Just be careful that even with Modified Henyey Greenstein you can still have problem with lookup table especially with g that are close to 1.

 

So at the end what I did is use the Newton Method  with the f = CDF - Rng (Rng is the random number in [0,1] to find the angle for the inverse CDF).

So to do all of that you need to 

1. Do the modification to change __constant__ float4 gproperty[MAX_PROP_AND_DETECTORS] in two different table __constant__ Medium property[MAX_PROP] and __constant__ float4 detectors[MAX_DETECTORS] (in /src/mcx_core.cu) and do all the modifications that is need because of the modification
2.  In the struct  MCXMedium (in /src/mcx_utils.h) you can add other parameter for your phase function like g, alpha or even the lookup table
3. Add in the function mcx_loadjson (in /src/mcx_utils.c) appropriate modification to read the new parameter in the .json file
4. in the mcx_core.cu do the modification where the costheta is found (you should fin it by searching for Henyey Greenstein in the file)

That's all what I did in a nutshell. 

Because there is other people that are looking to do the same thing. I'm gonna try to do something easy for people to modify in the code and share it on github (or something like that)

 

[Fang, Qianqian]

On 6/3/21 10:06 AM, 'Maxime Baillot' via mcx-users wrote:

There is,in general, formula for the CDF  but no formula for the Inverse of CDF (as far as we know or maybe we just need to read russian (every one know that the russians did everything in maths! hahaha))

Just be careful that even with Modified Henyey Greenstein you can still have problem with lookup table especially with g that are close to 1.

under-sampling could be an issue for very complex phase functions, but I assume for most use cases, sampling inv(CDF) at a few hundred samples should be sufficient.

So at the end what I did is use the Newton Method  with the f = CDF - Rng (Rng is the random number in [0,1] to find the angle for the inverse CDF).

why not use a linear index (equal-spaced sample) for invCDF? there is no look-up cost, just storage.

So to do all of that you need to 

1. Do the modification to change __constant__ float4 gproperty[MAX_PROP_AND_DETECTORS] in two different table __constant__ Medium property[MAX_PROP] and __constant__ float4 detectors[MAX_DETECTORS] (in /src/mcx_core.cu) and do all the modifications that is need because of the modification

maybe I overlook something - I don't see splitting the constant memory array gproperty to two (property/detectors) helps for implementing this feature - can you explain why you need to do this?

I would imagine that the shared memory might be the best place for this lookup table - it is the fastest, and currently under-utilized in mcx.

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/22fa0f3c-79cd-40fd-bb70-0abacd63cbc4n%40googlegroups.com.

[Maxime Baillot]

On Thursday, June 3, 2021 at 10:56:01 AM UTC-4 Fang, Qianqian wrote:

On 6/3/21 10:06 AM, 'Maxime Baillot' via mcx-users wrote:

There is,in general, formula for the CDF  but no formula for the Inverse of CDF (as far as we know or maybe we just need to read russian (every one know that the russians did everything in maths! hahaha))

Just be careful that even with Modified Henyey Greenstein you can still have problem with lookup table especially with g that are close to 1.

under-sampling could be an issue for very complex phase functions, but I assume for most use cases, sampling inv(CDF) at a few hundred samples should be sufficient.

In a lot of cases under sampling is not a problem. But even with Henyey Greenstein I saw some problem for g close to 1 like g = 0.98 or .99 (I did some comparative simulation that is why I speak about that). But for most biological media it's fine I think. I was looking at those values for other reason at first.

 

So at the end what I did is use the Newton Method  with the f = CDF - Rng (Rng is the random number in [0,1] to find the angle for the inverse CDF).

why not use a linear index (equal-spaced sample) for invCDF? there is no look-up cost, just storage.

I am not using a look-up table anymore now. That is why I use the Newton Method now. What I do is I find the corresponding costheta for a value of the CDF on the fly (it take a bit of time but it's ok for me)

So to do all of that you need to 

1. Do the modification to change __constant__ float4 gproperty[MAX_PROP_AND_DETECTORS] in two different table __constant__ Medium property[MAX_PROP] and __constant__ float4 detectors[MAX_DETECTORS] (in /src/mcx_core.cu) and do all the modifications that is need because of the modification

maybe I overlook something - I don't see splitting the constant memory array gproperty to two (property/detectors) helps for implementing this feature - can you explain why you need to do this?

I would imagine that the shared memory might be the best place for this lookup table - it is the fastest, and currently under-utilized in mcx.

I am adding other medium property in the struct MCXMedium so floart4 is not suitable anymore (MCXMedium have more than 4 floats in the structure now) 

[Fang, Qianqian]

hi Maxime and Anouk

since implementing user-defined phase function has been one of the outstanding TODOs on github (I added it back in 2015)

https://github.com/fangq/mcx/issues/13

I decided to take a step and implemented something that followed the lines that I suggested earlier. see my commits here

https://github.com/fangq/mcx/commit/ebf1ea1df701806b240c0571fa89b7461b82debc
https://github.com/fangq/mcx/commit/64c5dd061f81af8546658f42ee52eff33c672dcf

an example was added to show how to use this feature, including externally defined Henyey-Greenstein and Rayleigh scattering phase functions.

https://github.com/fangq/mcx/commit/ebf1ea1df701806b240c0571fa89b7461b82debc#diff-35d87707a598643dbe94356044238841f29b9b14fe208eea0b9ca3e7ed16b9cc

the needed shared memory is quite moderate and can accommodate higher resolution, the overall speed penalty is not noticeable.

please give it a try and let me know if you see any issue.

at this point, only a homogeneous phase function is supported, in the future, it can be extended for heterogeneous phase functions - by serializing cfg.invcdf for each scattering type and using g to store indices to find the corresponding invcdf for the corresponding medium label.

Qianqian

--
You received this message because you are subscribed to the Google Groups "mcx-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to mcx-users+...@googlegroups.com.

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/7d9c4372-7a84-4989-932e-a05c4a93d289n%40googlegroups.com.

[Maxime Baillot]

Hi,

That is great! I will take a look deeper in your change and try to do some tester next week (and also try to update pymcx to use it at the same time)

I could not do that by my self because of my lack of knowledge on how to do it and mcx internals working. 

I did the way I did it because I found it simpler to add the parameter in the medium struct and also because I'm using heteorgeneous phase function (if I understand I found a way to add that to your change I will do it and late you know)

I still think that adding the tabulated invCDF in the texture memory would be great because it is not use and it have an auto interpolation (if i'm not mistaking) but it would take a lot of work and I can't do it unfortunately.

Anyway I will come back to you as soon as I did my testing.

PS: again thanks for the modifications and sorry for not sharing more stuff that I did

[Qianqian Fang]

On 6/4/21 10:44 AM, 'Maxime Baillot' via mcx-users wrote:

Hi,

That is great! I will take a look deeper in your change and try to do some tester next week (and also try to update pymcx to use it at the same time)

thanks, can you also check out one of the recent questions regarding pymcx?

https://groups.google.com/g/mcx-users/c/xPj4Q3WVYwA/m/qt6RDf0oBAAJ

I could not do that by my self because of my lack of knowledge on how to do it and mcx internals working. 

I did the way I did it because I found it simpler to add the parameter in the medium struct and also because I'm using heteorgeneous phase function (if I understand I found a way to add that to your change I will do it and late you know)

I am sure it should work too - just my personal preference, I'd like to keep the float4 data structure whenever possible, because it is often the most optimized data unit in modern hardware (CPU/GPU) and can benefit from various caching mechanisms.

I still think that adding the tabulated invCDF in the texture memory would be great because it is not use and it have an auto interpolation (if i'm not mistaking) but it would take a lot of work and I can't do it unfortunately.

yes, texture memory does automatic interpolation. however, it will bring higher memory overhead compared to shared/constant memory, especially such data is needed in every scattering event. for 1d interpolation, doing interpolation manually is still ok

https://github.com/fangq/mcx/blob/master/src/mcx_core.cu#L1487-L1488

Anyway I will come back to you as soon as I did my testing.

PS: again thanks for the modifications and sorry for not sharing more stuff that I did

absolutely no worries!

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/4881862b-bd3e-4a16-b493-e578af6ae744n%40googlegroups.com.

[Anouk Post]

Dear Prof Fang,

Thank you so much for the modification! 

Unfortunately, I cannot test it, because I cannot compile it myself (I tried it before - I only have access to Windows). I did look at the commits and I was wondering if I understand it correctly. I think right now cfg.invcdf has a stored look-up table in it, with only the values of the inverse cdf, but not the corresponding mu values. The mu values are 0.01:0.01:0.99. Are those always the mu values the MC software will use? Or is there flexibility in choosing other mu values? 

For diffuse measurements this is probably enough, but for subdiffuse measurements I generally have look-up tables in the order of 10.000 mu values. Would it be possible to make it user-defined how much values the look-up table has? Or do you expect that such high values would make the simulations run much slower? With my adapted version of MCML it works on my laptop, but I'm not sure what the influence would be on GPU.

If you would like me to compare the results with the new MXC code to the adapted version of MCML that I currently have, I would happily do that. In that case, I would need a compiled version of the code since I can't seem to compile it on my own computer, unfortunately. 

Best wishes,

Anouk

You received this message because you are subscribed to a topic in the Google Groups "mcx-users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/mcx-users/pckAVbf8z7U/unsubscribe.
To unsubscribe from this group and all its topics, send an email to mcx-users+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/0c00205d-6d60-6691-0297-d3e0006f7cb1%40neu.edu.

[Fang, Qianqian]

On 6/8/21 10:03 AM, Anouk Post wrote:

Dear Prof Fang,

Thank you so much for the modification! 

Unfortunately, I cannot test it, because I cannot compile it myself (I tried it before - I only have access to Windows).

you don't have to recompile, our nightly build is already updated. you can download from here

http://mcx.space/nightly/win64/

I did look at the commits and I was wondering if I understand it correctly. I think right now cfg.invcdf has a stored look-up table in it, with only the values of the inverse cdf, but not the corresponding mu values. The mu values are 0.01:0.01:0.99. Are those always the mu values the MC software will use? Or is there flexibility in choosing other mu values?

it must be in the format of du:du:1-du, but you can decide the sample spacing (du) and 1/du must be an integer

For diffuse measurements this is probably enough, but for subdiffuse measurements I generally have look-up tables in the order of 10.000 mu values.

because typical shared memory is about 48kB on modern NVIDIA GPUs, so, you should be ok by setting du=1/10000, but I don't think it is necessary in most cases

Would it be possible to make it user-defined how much values the look-up table has?

yes, see above

Or do you expect that such high values would make the simulations run much slower? With my adapted version of MCML it works on my laptop, but I'm not sure what the influence would be on GPU.

If you would like me to compare the results with the new MXC code to the adapted version of MCML that I currently have, I would happily do that. In that case, I would need a compiled version of the code since I can't seem to compile it on my own computer, unfortunately.

that would be great, happy to see some cross-validation on this.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/CAOmd5%3D3va2Ux7BU_YhmRi2-ApOBWBfNcbiGUTT8fPnzTohvV5g%40mail.gmail.com.

[Maxime Baillot]

With what you where saying to have nonhomogeneous phase function I think I found a way to do it.

invcdf could be a table of table (pointer of pointer)

float **invcdf;

cfg->invcdf=(float**)calloc(cfg->medianum,sizeof(float*))

cfg->invcdf[idmedia] = (float*)calloc(cfg->nphase,sizeof(float));

something like that should work. I'm gonna try to implement this (I will need to add a lot of change in the lines you already change in the commit but it should be just adding the id of the media or medianum). Using something like that and we can still have the g as the g_1 value of the medium. The only problem that I can see is about memory management of the pointer of pointer.... if someone see where it could go wrong please let me know.

Also Anouk. I think you can use linspace(0,1,nbp) it should work with no problem especially if you want a big number of points. (The edge of the table will have two 0 and 1 but I don't think it's a problem with a big number of points)

[Qianqian Fang]

On 6/10/21 10:22 AM, 'Maxime Baillot' via mcx-users wrote:

With what you where saying to have nonhomogeneous phase function I think I found a way to do it.

invcdf could be a table of table (pointer of pointer)

float **invcdf;

cfg->invcdf=(float**)calloc(cfg->medianum,sizeof(float*))

cfg->invcdf[idmedia] = (float*)calloc(cfg->nphase,sizeof(float));

something like that should work. I'm gonna try to implement this (I will need to add a lot of change in the lines you already change in the commit but it should be just adding the id of the media or medianum). Using something like that and we can still have the g as the g_1 value of the medium. The only problem that I can see is about memory management of the pointer of pointer.... if someone see where it could go wrong please let me know.

use 2D arrays as pointer-of-pointers is only beneficial if you are dealing with different sampling for phase functions in different labels. Unless some of the phase functions are extremely spiky (Mie-like) while others are smooth, I don't see much added benefit to accept label-varying lookup tables - but certainly you can do it.

in either case (label-varying sampling vs uniform sampling), a 2D array should to be serialized into a 1D buffer in the GPU like we did for cfg.vol, cfg.srcpattern etc.

Declaring a float* pointer will cost you 8-byte (sizeof(float*)) register space - a pointer-of-pointers will consume (cfg->mediumnum*8) bytes of register space to store just address information, this is a pretty big overhead in the GPU kernel.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/6b521db0-fdc5-440f-927b-22ee46bb9ccdn%40googlegroups.com.

[Maxime Baillot]

Hello Dr Frang and Anouk,

I did the serialization for the phase function : https://github.com/4D42/mcx/tree/dev_NinvCDF

I did some test and it look like it work.

For the moment it can only be use with a .json file. "InverseCDF" is a list of  N list (N is the number of medium) of invCDF list. 

• In the json file :"Domain": { "InverseCDF": [[, , , ], [, , ,] , [, , ,]] }
• For Pymcx : cfg['Domain']['InverseCDF'] = [ [ ], [ ] ,[ ]]

If you find any mistake please let me know (I did some modification that I'm not 100% certain) 

Meanwhile I will continue to do other test and update the branch if I find bugs or problems.

[Anouk Post]

Dear Prof. Fang,

I had to finish some work for a grant proposal, so I didn't have any time earlier to test the phase function implementation. But I started testing it today I want to compare the reflectance versus radial distance from the source for different phase functions. For that I generally use information from detphoton. However, when I run:

[fluence,detphoton,vol,seed,trajectory]=mcxlab(cfg);

in combination with cfg.invcfg I get the following error:

However, if I use:

flux = mcxlab(cfg);

in combination cfg.invcfg it works fine.

Do you have an idea what could cause this error?

Thank you in advance.

Kind regards,

Anouk

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/2b661a13-5ccf-5343-ecda-3b91d21c848a%40neu.edu.

[Qianqian Fang]

On 7/7/21 6:52 AM, Anouk Post wrote:

Dear Prof. Fang,

I had to finish some work for a grant proposal, so I didn't have any time earlier to test the phase function implementation. But I started testing it today I want to compare the reflectance versus radial distance from the source for different phase functions. For that I generally use information from detphoton. However, when I run:

[fluence,detphoton,vol,seed,trajectory]=mcxlab(cfg);

in combination with cfg.invcfg I get the following error:

hi Anouk,

I was able to reproduce the issue you reported. when the seed output is included in mcxlab's output, the cfg.invcdf data structure caused a misalignment error in the shared memory.

I created the below issue on github explaining the cause

https://github.com/fangq/mcx/issues/118

the cfg.invcdf created in the demo_mcxlab_phasefun.m script has a length of 99, adding the 2 ends, making the invcdf buffer a length of 101, which is an odd number. this makes reading/writing the seed data across byte alignment boundaries.

I just committed a fix to this issue and a quick test shows that the error is gone.

https://github.com/fangq/mcx/commit/53d7ac0f2b26ca7cc59b5a50466f228b9fc0dbe1

please redownload the nightly build and let me know if you see any other issues.

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/CAOmd5%3D3Gf80xcNVaxHCyhwj6tyc7eDpYHHpnVyd9nxvN8_j%2BEA%40mail.gmail.com.

[x.at...@gmail.com]

Dear all,

First off thank you very much for the implementation of the lookup table-based phase function. It is much appreciated. A directly related request would be the implementation of a source type which allows the user to define the angular emission profile of a source. This is relevant for the accurate modelling of multimode optical fibers, which can be represented as a disk source with uniform initial intensity, where each point emits photons following a certain angular distribution. The simplest assumption is that this angular distribution is uniform within the acceptance angle of the fiber. This model can be expanded to more complex forms where the probability of a photon having a certain origin and initial angle can be described by the radiance function (see attached link for a more complete example/explanation).

https://www.osapublishing.org/ao/fulltext.cfm?uri=ao-57-5-1110&id=381396

I believe implementing this lookup table should not significantly affect overall runtime as it would only be evaluated once for each photon packet rather than once for each scattering interaction as is the case for the phase function. 

Would something like this be possible ?

Best regards,

Xavier

[Qianqian Fang]

On 10/19/21 9:45 AM, x.at...@gmail.com wrote:

Dear all,

First off thank you very much for the implementation of the lookup table-based phase function. It is much appreciated. A directly related request would be the implementation of a source type which allows the user to define the angular emission profile of a source. This is relevant for the accurate modelling of multimode optical fibers, which can be represented as a disk source with uniform initial intensity, where each point emits photons following a certain angular distribution. The simplest assumption is that this angular distribution is uniform within the acceptance angle of the fiber. This model can be expanded to more complex forms where the probability of a photon having a certain origin and initial angle can be described by the radiance function (see attached link for a more complete example/explanation).

https://www.osapublishing.org/ao/fulltext.cfm?uri=ao-57-5-1110&id=381396

I believe implementing this lookup table should not significantly affect overall runtime as it would only be evaluated once for each photon packet rather than once for each scattering interaction as is the case for the phase function. 

Would something like this be possible ?

hi Xavier, yes, this feature has come across my mind as well. I've added this todo on github to track our progress

https://github.com/fangq/mcx/issues/129

Qianqian

To view this discussion on the web visit https://groups.google.com/d/msgid/mcx-users/ded5b1e8-eec0-488e-bb71-d29d913d6844n%40googlegroups.com.

[x.at...@gmail.com]

Dear Qianqian,

Thank you for adding the controllable angular distribution to the todo. It will definitely be very useful to our research. Is there anything I can do to assist with regards to this?

Best regards,

Xavier