spot focused Gaussian source distribution

430 views
Skip to first unread message

G Buist

unread,
Aug 2, 2021, 3:20:15 AM8/2/21
to mcx-users
Dear Qianqian,

In the MCX code is there an option available for a spot focused Gaussian beam similar to;
or
?

Best regards,
Gijs Buist

Qianqian Fang

unread,
Aug 3, 2021, 4:41:05 PM8/3/21
to mcx-...@googlegroups.com
On 8/2/21 3:20 AM, 'G Buist' via mcx-users wrote:
Dear Qianqian,

In the MCX code is there an option available for a spot focused Gaussian beam similar to;
or
?


hi Gijs,

mcx supports collimated Gaussian beam (srctype='gaussian'), convergent/divergent gaussian beam (srcdir(4)>0 or <0), or angular Gaussian beam (srctype='zgaussian').

For example, collimated Gaussian beam

http://mcx.space/cloud/?tab=preview&data=N4Ig9gDgLglmB2BnEAuUMDGCA2MBGqIAZglAIYDuApomALZUCsIANOHgFZUZQD62ZAJ5gArlELwwAJzplsrEIgAWYCrypSp05ChAx4UDWR5x4CqEqoNCiCNyhSqCybzIATNzFgI5vCFrspWBpUBxEqNg5aMzQQAGUaRFNUUABJABFCOgwADwUABRUoBB0ARgAGSsq2dLAAWRhEWSgMJVCpcJqwAEExMAAVJUd3ds6QWriyADcqADUwbBFrFDCI8bB8siCYOU2LUbXagDlpWVwALycVjsOwSZmAJSoiVCI5RFv7qgBRHK9X96faZUBJUNwA7AfLrpADCcQhUPWcTsGBEAikBxqVDwIgA5gAxAS4wgKL7pMjkOpkRAAaxJbAA8mIIGJ8acKRIYDAFEyoCyoP1BHZCHk2A8jgBxUHglClABsABYABwAZhVjEYpSVAF82GypBQtjLQP1yqhymx+qVUIwqABaACcNXEKFtjt1IAZ0DAbiuoDioikGD9IEFwt0uLIIkQSTIZjY+TAOgA2iqLWmWOUALo1GAYlDJi0W0o5kCbKRkOjWgsVTN17MJraVgBMqEL9czpfxjgAjuF4BhBOaPekqIYeNI21mPXElGQ7CnQBKpDBjaGyMTZWw4jBLm25RaDywD9PtaXarJ9ClPSvcfow1dVrm6PvD2/SwAZKjwXEWACq8D/FuIB1GCOz6s0hB4IIhgKPMiwMPiMDYFcIAKKBnhkG2oB0CIWEoBaIC4TohGbqUbAxKUuo4Xh5oAHRVORRHRuabBkRRqBUae2rakAA=


a point-focusing Gaussian beam

http://mcx.space/cloud/?tab=preview&data=N4Ig9gDgLglmB2BnEAuUMDGCA2MBGqIAZglAIYDuApomALZUCsIANOHgFZUZQD62ZAJ5gArlELwwAJzplsrEIgAWYCrypSp05ChAx4UDWR5x4CqEqoNCiCNyhSqCybzIATNzFgI5vCFrspWBpUBxEqNg5aMzQQAGUaRFNUUABJABFCOgwADwUABRUoBB0ARgAGSsq2dLAAWRhEWSgMJVCpcJqwAEExMAAVJUd3ds6QWriyADcqADUwbBFrFDCI8bB8siCYOU2LUbXagDlpWVwALycVjsOwSZmAJSoiVCI5RFv7qgBRHK9X96faZUBJUNwA7AfLrpADCcQhUPWcTsGBEAikBxqVDwIgA5gAxAS4wgKL7pMjkOpkRAAaxJbAA8mIIGJ8acKRIYDAFEyoCyoP1BHZCHk2A8jgBxUHglClABsABYABwAZhVjEYpSVAF82GypBQtjLQP1yqhymx+qVUIwqABaACcNXEKFtjt1IAZ0DAbiuoDioikGD9IEFwt0uLIIkQSTIZjY+TAOgA2iqLWmWOUALo1GAYlDJi0W0osNM5kCbKRkOjWgsAJiLmcz5cr1brqELTaL5fxjgAjuF4BhBOaPekqIYeNIO1mPXElGQ7CnQBKpDBjaGyMTZWw4jBLh25Rajywj7PteXarJ9ClPWvcfow1dVrm6Ifjx/ywAZKjwXEWABVeB/h3EA6jBHZ9WaQg8EEQwFHmRYGHxGBsCuEAFHAzwyA7UA6BEHCUAtEB8J0YjtxLEAYlKXU8II80ADoqko0jzTYCi2Goi9tR47UgA

the focused gaussian beam output screenshot is attached.


Qianqian


Best regards,
Gijs Buist
--
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/35508d2d-2155-4f17-8f7b-d3fcf7cd55e5n%40googlegroups.com.
point_focusing_gaussian.png

Buist, G.

unread,
Aug 4, 2021, 9:35:45 AM8/4/21
to mcx-...@googlegroups.com

Hi Qianqian,

 

Does the focussed Gaussian of MCX work via method b) in the figure below?

I wish to simulated a setup with geometry b) where the focus inside the medium.

 

 

 


 

Best regards,

Gijs Buist

--
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/wauKd1IbEJE/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/f747312e-eb93-9f2a-fb02-cc27fb738bdf%40neu.edu.

Fang, Qianqian

unread,
Aug 4, 2021, 12:22:15 PM8/4/21
to mcx-...@googlegroups.com
On 8/4/21 9:35 AM, 'Buist, G.' via mcx-users wrote:

Hi Qianqian,

 

Does the focussed Gaussian of MCX work via method b) in the figure below?

I wish to simulated a setup with geometry b) where the focus inside the medium.


hi Gijs, the one that I showed earlier is method (a). true gaussian beam is not yet supported in mcx.


Buist, G.

unread,
Sep 7, 2021, 9:39:25 AM9/7/21
to mcx-...@googlegroups.com

Hi Qianqian,


I have written code to determine the start positions and directions of photons to be launched in order to create a true Gaussian beam.

It works by realizing that the Gaussian beam is a hyperboloid of one sheet and that such shapes can be fully described by straight lines.

And needs  scrparam1=(w,zf,zr) (radius beam waist at focus, z-position focus, rayleigh range) as input.


I now wish to add this code to the MCX code but am not sure on how exactly.

I believe it should look somewhat like this:


case(MCX_SRC_TRUE_GAUSSIAN): { // Gaussian-beam with spot focus
              //hyperboloid of one sheet, doubly ruled surface
              // scrparam1=(w,zf,zr) beam waist, z-position focus, rayleigh range
              float sphi, cphi;
              float phi=TWO_PI*rand_uniform01(t);
              sincosf(phi,&sphi,&cphi);
                 
            float r; // sampled distance from z-axis at focus
            float x0; // sampled x coordinate photon at focus
            float y0; // sampled y coordinate photon at focus
            
            r=sqrtf(-0.5f*logf(rand_uniform01(t)))*gcfg->srcparam1.x;
            x0=r*cphi;
            y0=r*sphi;
            
            float xi; // x position photon at entry
            float yi; // y position photon at entry
            float zi=0; // z position photon at entry (can be any value but then zf difference between focus position and zi)
            float t = -gcfc->srcparam1.y/gcfc->scrparam1.z;
            
            xi= r*(cphi-t*sphi);
            yi= r*(sphi+tcphi);
            
            float l; //use rsqrtf instead to calculate 1/l ? used to normalize directional vector
            float vx; // normalized x part of directional vector photon
            float vy; // normalized y part of directional vector photon
            float vz; // normalized z part of directional vector photon
            
            l=sqrtf(r*r+gcfc->srcparam1.z*gcfc->srcparam1.z); //use rsqrtf instead to calculate 1/l?
            vx=-r*sphi/l;
            vy=r*cphi/l;
            vz=gcfc->srcparam1.z/l
            
            // How to load xi,yi,zi in photon start position ?
            // How to load vx,vy,vz in photon direction ?           
 
          }

But I am not sure how integrate it with the rest of the MCX code.

Could you tell me how I get xi,yi and zi to be the start position of the photon in the MCX code, as well as how to make vx,vy,vz the directional components?

Many thanks in advance.


Gijs Buist


Van: mcx-...@googlegroups.com <mcx-...@googlegroups.com> namens Fang, Qianqian <q.f...@northeastern.edu>
Verzonden: woensdag 4 augustus 2021 18:22:11
Aan: mcx-...@googlegroups.com
Onderwerp: Re: [mcx-users] spot focused Gaussian source distribution
 

Fang, Qianqian

unread,
Sep 7, 2021, 1:00:38 PM9/7/21
to mcx-...@googlegroups.com
On 9/7/21 9:39 AM, 'Buist, G.' via mcx-users wrote:

Hi Qianqian,


I have written code to determine the start positions and directions of photons to be launched in order to create a true Gaussian beam.

It works by realizing that the Gaussian beam is a hyperboloid of one sheet and that such shapes can be fully described by straight lines.

And needs  scrparam1=(w,zf,zr) (radius beam waist at focus, z-position focus, rayleigh range) as input.



hi Gijs, thanks for contributing.


currently, all mcx sources can be pointed to any direction, but reading your below code, it seems it is only oriented along the +z direction.

so, in order for me to integrate this code, I would like to understand this better and generalize it for all directions.


can you send me offline a drawing/slides showing how is this source implemented? I just want to get some idea how is this

different from the existing point-focusing Gaussian and how to add it with minimal lines of codes. please email that drawing/schematic

in a reply or to my email directly.


Qianqian



Buist, G.

unread,
Sep 10, 2021, 8:12:07 AM9/10/21
to mcx-...@googlegroups.com

Hi Qianqian,

 

I indeed started with the simple case of propagating along the z-direction.

Attached are slides explaining the math.

 

Regards,

Gijs Buist

true_gaussian_MCX_distribution_presi.pptx

Fang, Qianqian

unread,
Sep 10, 2021, 10:42:03 AM9/10/21
to mcx-...@googlegroups.com
On 9/10/21 8:12 AM, 'Buist, G.' via mcx-users wrote:

Hi Qianqian,

 

I indeed started with the simple case of propagating along the z-direction.

Attached are slides explaining the math.


hi Gijs, thank you for your detailed slides! well done and thanks for sharing with the list.

I have a busy schedule today but will take a look during the weekend. I may likely follow up and schedule a meeting with you next week and learn more about your implementation.

Qianqian


Buist, G.

unread,
Sep 21, 2021, 6:19:45 AM9/21/21
to mcx-...@googlegroups.com

Hi Qianqian,

 

If you have any questions I would be happy to schedule a meeting.

Additionally I have been thinking on how to create the gaussians for any direction and I believe the best way is by rotation and translation of the ‘standard’ z-propagation centred at the origin result.

 

Regards,

Gijs

Fang, Qianqian

unread,
Sep 21, 2021, 12:03:58 PM9/21/21
to mcx-...@googlegroups.com, Shijie Yan
On 9/21/21 6:19 AM, 'Buist, G.' via mcx-users wrote:

Hi Qianqian,

 

If you have any questions I would be happy to schedule a meeting.

Additionally I have been thinking on how to create the gaussians for any direction and I believe the best way is by rotation and translation of the ‘standard’ z-propagation centred at the origin result.


that's what I thought too. we've already have a 3-D vector rotation function and have been using it for launching photons

https://github.com/fangq/mcx/blob/master/src/mcx_core.cu#L902
https://github.com/fangq/mcx/blob/master/src/mcx_core.cu#L1201-L1203

can you update your code and add rotation to allow it to point to srcdir?

I am CCing my student Shijie Yan who can answer questions about our part of the code.

Qianqian


Buist, G.

unread,
Sep 24, 2021, 3:23:51 AM9/24/21
to mcx-...@googlegroups.com, Shijie Yan

Dear,

 

I wrote and tested a first draft of code in Python that should give a true gaussian for arbitrary source position and direction (see below).

First the position and direction for source at (0,0,0) with dir (0,0,1) is calculated and then based on the desired direction a rotation around the y-axis and z-axis (in that order) is preformed on the direction and position.

Then the position is shifted by the desired source position.

I did not fully understand the rotatevector of MCX so the code still needs to be translated to .cu

 

Best,

Gijs

 

PYTHON CODE

###########################################################

Import numpy as np

 

# INPUT PARAMETERS NECESSARY FOR PHOTON LAUNCH OF TRUE GAUSSIAN FROM ARBITRARY POSITION AND ARBITRARY BEAM DIRECTION

 

#beam parameters

 

w0 = 1 #radius beam waist 1/e reduction E-field (1/e2 for intensity) at focus

df = -1 # distance between launch plane and focus

zr = 2 #rayleigh range

 

# direction and position beam for rotation and shift respectively

beam_dir = np.array([1,0,0]) # (x,y,z) beam direction

source_pos = np.array([2,3,0]) # (x,y,z) position center beam source

 

# GET SOLUTION FOR SOURCE AT ORIGIN PROPAGATING IN +Z DIRECTION

 

rng = np.random.default_rng()

num_phot=1 # In MCX core it should be a single photon launch function, here arrays of num_phot length

ran1 = rng.uniform(low=0,high=1.0, size=num_phot) # for distance from beam center

ran2 = rng.uniform(low=0,high=1.0, size=num_phot) # for azimuthal angle

 

#r_samps = r = w0*np.sqrt(-0.5*np.log(ran1)) # no sampling to test rotations reliably

r_samps = np.zeros(num_phot)+1.0 # all at r=1.0

thet_samps = 2*np.pi*ran2 # azimuthal angle of photon before rotation

 

x0 = r_samps*np.cos(thet_samps) # x-position at focus

y0 = r_samps*np.sin(thet_samps) # y-position at focus

 

xi,yi = entry_cor(r_samps,thet_samps,zr,df) #positions at launch plane

zi = np.zeros(len(xi))                      # launch at plane z=0

vx,vy,vz = dir_zf0(r_samps,thet_samps,zr)   # components directional vector photon before rotation

 

 

# GET (COSINES AND SINES OF) ROTATION ANGLES TO GO FROM (0,0,1) TO BEAM_DIR ([0],[1],[2])

# Rotation of phi around y-axis and rotation of theta around z-axis in that order

# as rotations around z-axis do not do anything for (0,0,1) vector

 

# theta is counter-clockwise angle with x-axis (azimuthal angle) (from +x axis to +y axis)

# phi is angle with z-axis towards positive x-axis (zenith) (+z axis to +x axis)

 

# to prevent division by zero

if beam_dir[2]==1:

    0 #no rotations necessary

 

else:

    tmp0 = 1/np.sqrt(beam_dir[0]**2 + beam_dir[1]**2) #1/length in xy plane

    tmp1 = 1/np.sqrt(beam_dir[0]**2 + beam_dir[1]**2+beam_dir[2]**2) #1/total_length

   

    ctheta = beam_dir[0]*tmp0

    stheta = beam_dir[1]*tmp0

    cphi = beam_dir[2]*tmp1

    sphi = np.sqrt(beam_dir[0]**2+beam_dir[1]**2)*tmp1

 

    # ROTATE LAUNCH POSITION AND DIRECTION

    # first rotation around y-axis then around z-axis.

    # Rotating first around z-axis does not do anything due to symmetry

    # Replace with the rotatevector of mcx if possible

   

    # rotated launch coordinates, zi part left out as zi=0

    xrot = (xi*ctheta*cphi - yi*stheta) #+zi*ctheta*sphi

    yrot = (xi*stheta*cphi + yi*ctheta) #+ zi*stheta*sphi

    zrot = (-xi*sphi) #+ yi*0 + zi*cphi

   

    # rotated directional vector

    vxrot = vx*ctheta*cphi - vy*stheta +vz*ctheta*sphi

    vyrot = vx*stheta*cphi + vy*ctheta + vz*stheta*sphi

    vzrot = -vx*sphi + vy*0 + vz*cphi

   

# SHIFT LAUNCH POSITION OF PHOTON WITH SOURCE_POS ([0],[1],[2])

 

xrot_shift = xrot + source_pos[0]

yrot_shift = yrot + source_pos[1]

zrot_shift = zrot + source_pos[2]

 

# RETURN FINAL LAUNCH POSITION AND DIRECTION OF PHOTON

 

phot_launch_pos = np.array([xrot_shift,yrot_shift,zrot_shift])

phot_launch_dir = np.array([vxrot,vyrot,vzrot])

 

#########################################################################################

Buist, G.

unread,
Sep 30, 2021, 12:01:19 PM9/30/21
to mcx-...@googlegroups.com, Shijie Yan

Hi,

 

I believe that the code for a focused Gaussian with arbitrary direction and position should look something like the code below.

However, I am not sure if the naming is consistent with the rest of the MCX code nor what to do with the resulting coordinates and directions.

 

Regards,

Gijs

 

case(MCX_SRC_TRUE_GAUSSIAN): { // Gaussian-beam with spot focus

              //hyperboloid of one sheet, doubly ruled surface

              // scrparam1=(w,df,zr) beam waist, distance from source to focus, rayleigh range

              // srcdir = (dx,dy,dz) direction of gausssian beam

              // srcpos = (x,y,z) location of the center of the gaussian source

           

            // Sampling for srcdir = (0,0,1) and srcpos = (0,0,0)

            float sphi, cphi;

            float phi=TWO_PI*rand_uniform01(t);

            sincosf(phi,&sphi,&cphi);

                

            float r; // sampled distance from z-axis at focus

            float x0; // sampled x coordinate photon at focus

            float y0; // sampled y coordinate photon at focus

           

            r=sqrtf(-0.5f*logf(rand_uniform01(t)))*gcfg->srcparam1.x;

            x0=r*cphi;

            y0=r*sphi;

           

            float xi; // x position photon at source

            float yi; // y position photon at source

            float zi=0; // z position photon at source

            float t = -gcfc->srcparam1.y/gcfc->scrparam1.z; //parameter to generate photon path from coordinates at focus (depends focal distance and rayleigh range)

           

            xi= r*(cphi-t*sphi);

            yi= r*(sphi+tcphi);

           

            float l; //use rsqrtf instead to calculate 1/l ? used to normalize directional vector

            float vx; // normalized x part of directional vector photon

            float vy; // normalized y part of directional vector photon

            float vz; // normalized z part of directional vector photon

           

            l=sqrtf(r*r+gcfc->srcparam1.z*gcfc->srcparam1.z); //use rsqrtf instead to calculate 1/l?

            vx=-r*sphi/l;

            vy=r*cphi/l;

            vz=gcfc->srcparam1.z/l;

           

            // rotate photon direction and position if srcdir!=(0,0,1)

            // rotatevector function updates photon direction so cannot be used for position

           

            if( v->z>-1.f+EPS && v->z<1.f-EPS ) {

                float tmp0=1.f-srcdir->z*srcdir->z; // if srcdir unit length then this is transverse direction squared

               

                float ctheta = srcdir->x*rsqrt(tmp0); // cos(phi), phi counter clockwise angle with x-axis (rotation of phi around z-axis)

                float stheta = srcdir->y*rsqrt(tmp0);

               

                float cphi = srcdir->z; // only in case of unit vector srcdir, theta angle with z-axis (rotation towards +x)

                float sphi = sqrt(tmp0); // sqrt(dx*dx+dy*dy)/1

               

                // position photon after rotating (xi,yi,zi=0)=>(xrot,yrot,zrot)

                float xrot = (xi*ctheta*cphi - yi*stheta) 

                float yrot = (xi*stheta*cphi + yi*ctheta)

                float zrot = (-xi*sphi)

                

                // add the shifted source position to the photon

                xrot+=srcpos->x

                yrot+=scrpos->y

                zrot+=srcpos->z

               

                // rotate direction of photon

                // could be done by rotatevector function?

   

                float vxrot = vx*ctheta*cphi - vy*stheta +vz*ctheta*sphi

                float vyrot = vx*stheta*cphi + vy*ctheta + vz*stheta*sphi

                float vzrot = -vx*sphi + vy*0 + vz*cphi

               

                // photon position (xrot,yrot,zrot) now need to be changed in the correct object

                // photon direction (vxrot,vyrott,vzrot) now need to be changed in the correct object

               

            }

            

            // no rotation happened as srcdir==(0,0,1)

            // photon position (xi,yi,zi) now need to be changed in the correct object

            // photon direction (vx,vy,vz) now need to be changed in the correct object

Qianqian Fang

unread,
Sep 30, 2021, 12:12:44 PM9/30/21
to mcx-...@googlegroups.com, Shijie Yan

hi Gijs,

thanks a lot for the patch. unfortunately my coding time becomes extremely limited as my dept service load eats most of my time. can only answer short questions, but for testing, I am CCing my student Shijie -

Shijie, can you please test Gijs's patch? I think this will be a great addition to mcx! please read the slides that Gijs sent previously to understand the algorithm. if everything works well, please create a PR on github. use new sourcetype string 'hyperboloid' and MCX_SRC_HYPERBOLOID_GAUSSIAN macro.


Qianqian

Shijie Yan

unread,
Oct 4, 2021, 8:36:30 PM10/4/21
to Fang, Qianqian, mcx-...@googlegroups.com
Hi Gijs,

After reviewing the code, I have a few questions:
  1. What are x0 and y0 respectively? It seems that they are assigned but not used in the following code.
  2. Do you have a working example in Python or MATLAB? Currently I am not able to run the python script because the functions "entry_cor" and "dir_zf0" are missing.
  3. The function "rotatevector" updates the global direction cosine given the sampled local polar (phi) and azimuthal (theta) angles (see the figure attached). The math is detailed in the section 3.4.7 "Scattering of a Photon Packet" of the textbook "Biomedical Optics: Principles and Imaging (2007)".
Thanks,
Shijie

From: Fang, Qianqian <q.f...@northeastern.edu>
Sent: Thursday, September 30, 2021 12:12 PM
To: mcx-...@googlegroups.com <mcx-...@googlegroups.com>

Shijie Yan

unread,
Oct 4, 2021, 8:37:40 PM10/4/21
to Fang, Qianqian, mcx-...@googlegroups.com
Selection_101.png

Buist, G.

unread,
Oct 6, 2021, 4:45:24 AM10/6/21
to mcx-...@googlegroups.com

Hi Shijie,

 

Attached you will find a complete python script for the hyperboloid source distribution.

 

As for your questions;

  1. x0 and y0 are the x and y position at the focus and are not needed and I forgot to delete them
  2. The working example can be found in the attachments

 

Best,

Gijs

hyperboloid_gaussian.py

Shijie Yan

unread,
Oct 6, 2021, 2:19:30 PM10/6/21
to mcx-...@googlegroups.com
Hi Gijs,

I can run the script and the output is a figure showing "unrotated launch positions". I will take a look and get back to you.

Thanks,
Shijie

From: 'Buist, G.' via mcx-users <mcx-...@googlegroups.com>
Sent: Wednesday, October 6, 2021 4:45 AM
To: mcx-...@googlegroups.com <mcx-...@googlegroups.com>

Buist, G.

unread,
Oct 7, 2021, 3:01:43 AM10/7/21
to mcx-...@googlegroups.com

Hi Shijie,

 

There should be multiple figures after the one you mentioned. I usually use Spyder to execute the code which gives all of the figures, if the Python command is used the next figure will appear after the previous figure is closed. Therefore a different figure should appear after closing “unrotated launch positions”.

 

Additionally about the rotatevector in MCX, if I understand correctly it updates the direction vector based the scattering angles (which are local properties).

In that case I do not think it would be useful for rotating the Gaussian beam as you do not know the local rotation (scattering angles) but the global rotation of the whole beam. I therefore determine the (global) sines and cosines of the rotation from (0,0,1) to an arbitrary vector of unit length. After determining those once they can be used in the rotation matrix RzRy (see Wikipedia) to rotate all the direction vectors. Finally as the hyperboloid is an extended object all the starting coordinates have to be rotated too.

 

Would something work in the MCX code?

 

case(MCX_SRC_TRUE_GAUSSIAN): { // Gaussian-beam with spot focus

              //hyperboloid of one sheet, doubly ruled surface

              // scrparam1=(w,df,zr) beam waist, distance from source to focus, rayleigh range

              // srcdir = (dx,dy,dz) direction of gausssian beam

              // srcpos = (x,y,z) location of the center of the gaussian source

           

            // Sampling for srcdir = (0,0,1) and srcpos = (0,0,0)

            float sphi, cphi;

            float phi=TWO_PI*rand_uniform01(t);

            sincosf(phi,&sphi,&cphi);

                

            float r; // sampled distance from z-axis at focus

           

            r=sqrtf(-0.5f*logf(rand_uniform01(t)))*gcfg->srcparam1.x;

           

            float xi; // x position photon at source

            float yi; // y position photon at source

            float zi=0; // z position photon at source

            float t = -gcfc->srcparam1.y/gcfc->scrparam1.z; //parameter to generate photon path from coordinates at focus (depends focal distance and rayleigh range)

           

            xi= r*(cphi-t*sphi);

            yi= r*(sphi+t*cphi);

            

            float l; //use rsqrtf instead to calculate 1/l ? used to normalize directional vector

            float vx; // normalized x part of directional vector photon

            float vy; // normalized y part of directional vector photon

            float vz; // normalized z part of directional vector photon

           

            l=rsqrtf(r*r+gcfc->srcparam1.z*gcfc->srcparam1.z); //use rsqrtf instead to calculate 1/l?

            vx=-r*sphi*l;

            vy=r*cphi*l;

            vz=gcfc->srcparam1.z*l;

           

            // calculate global rotation matrix and rotate photon direction and position if srcdir!=(0,0,1)         

            if( v->z>-1.f+EPS && v->z<1.f-EPS ) {

                // here the sines and cosines for the rotation matrix are determined

                // these are global properties and thus the same for all photons

                // the way it is now it is calculated every time but maybe there is a good way to change that

                float tmp0=1.f-srcdir->z*srcdir->z; // if srcdir unit length then this is transverse direction squared

               

                float ctheta = srcdir->x*rsqrt(tmp0); // cos(phi), phi counter clockwise angle with x-axis (rotation of phi around z-axis)

                float stheta = srcdir->y*rsqrt(tmp0);

               

                float cphi = srcdir->z; // only in case of unit vector srcdir, theta angle with z-axis (rotation towards +x)

                float sphi = sqrt(tmp0); // sqrt(dx*dx+dy*dy)/1

               

                // here the rotations based on the global rotation matrix are done

                // for the photon launch coordinates and directional vector

               

                // position photon after rotating (xi,yi,zi=0)=>(xrot,yrot,zrot)

                float xrot = (xi*ctheta*cphi - yi*stheta) 

                float yrot = (xi*stheta*cphi + yi*ctheta)

                float zrot = (-xi*sphi)

                

                // add the shifted source position to the photon

                xrot+=srcpos->x

                yrot+=scrpos->y

                zrot+=srcpos->z

               

                // rotate direction of photon

   

                float vxrot = vx*ctheta*cphi - vy*stheta +vz*ctheta*sphi

                float vyrot = vx*stheta*cphi + vy*ctheta + vz*stheta*sphi

                float vzrot = -vx*sphi + vy*0 + vz*cphi

               

                // photon position (xrot,yrot,zrot) now need to be changed in the correct object

                // photon direction (vxrot,vyrott,vzrot) now need to be changed in the correct object

               

            }

            

            // no rotation happened as srcdir==(0,0,1)

            // photon position (xi,yi,zi) now need to be changed in the correct object

            // photon direction (vx,vy,vz) now need to be changed in the correct object

 

           

                                }

 

Best,

Shijie Yan

unread,
Oct 10, 2021, 11:32:57 PM10/10/21
to mcx-...@googlegroups.com
Hi Gijs,

Sorry for the delayed response. It took me a while to walk through your example and understand the algorithm. I have finished the initial implementation for MCX. The CUDA code should now match what you did in your python script. However, when running the python script, I noticed that the cross-section of the one-sheet hyperboloid changes from a circle to an ellipse after rotation. I just want to make sure that this is correct (I suppose that you want to preserve the shape after rotation).

Thanks,
Shijie

From: 'Buist, G.' via mcx-users <mcx-...@googlegroups.com>
Sent: Thursday, October 7, 2021 3:01 AM

Buist, G.

unread,
Oct 11, 2021, 3:26:52 AM10/11/21
to mcx-...@googlegroups.com

Hi Shijie,

 

Thanks for your response and if you have any questions about the code I am happy to answer them. As for the cross-section of the hyperboloid, this should indeed remain circular but I have not extensively tested this. Just looking at the 3d plots seems to be of no help as the perspective makes gauging the precise 3d shape difficult. To check this robustly the distance to the optical axis (at the same t-parameter) needs to be calculated for every photon. At the bottom of the attached script this is done for a beam_dir = (0,0,1), (0,1,0) and (1,0,0) (different lines need to be (un)commented for different beams). For arbitrary beam_dir the coordinates of the optical axis for all the t-parameter values need to be determined for which I have currently no function.

I hope this clears things up a bit.

hyperboloid_gaussian.py

Shijie Yan

unread,
Oct 11, 2021, 10:25:45 AM10/11/21
to mcx-...@googlegroups.com
Hi Gijs,

I am not entirely sure how to use the code attached at the bottom. However, as shown in the script you just sent (where the beam direction is [0 1 0]), the photon positions at the launch plane formed a circle (figure #1), which was transformed to an ellipse after rotation. I assume these positions have the same parameter t. Please correct me if I am wrong.

Thanks,
Shijie

From: 'Buist, G.' via mcx-users <mcx-...@googlegroups.com>
Sent: Monday, October 11, 2021 3:26 AM
Selection_103.png
Selection_104.png

Shijie Yan

unread,
Oct 11, 2021, 10:47:08 AM10/11/21
to mcx-...@googlegroups.com
Hi Gijs,

I guess I figure it out that the plot issue is due to unequal axis (not sure about how to set equal axis in matplotlib though). Sorry for the confusion.

Shijie

Sent: Monday, October 11, 2021 10:25 AM
To: mcx-...@googlegroups.com <mcx-...@googlegroups.com>
Subject: Re: [mcx-users] spot focused Gaussian source distribution
 

Shijie Yan

unread,
Oct 12, 2021, 8:00:21 PM10/12/21
to mcx-...@googlegroups.com
Hi Gijs,

I have completed the initial support of the hyperboloid gaussian source in MCXLAB. The code should match the attached python script (adapted from your python script): 
  1. Edge-cases such as [0 0 1] and [0 0 -1] are added.
  2. Definitions of phi and theta are swapped. 
Please feel free to download and test it from: https://github.com/ShijieYan/mcx/tree/hyperboloid. An example MATLAB script is attached.

Shijie

Sent: Monday, October 11, 2021 10:47 AM
hyperboloid_gaussian_shijie.py
test_true_gaussian.m

Buist, G.

unread,
Oct 13, 2021, 5:49:48 AM10/13/21
to mcx-...@googlegroups.com

Hi Shijie,

 

Thank you for your work in implementing the hypberboloid gaussian.  At the moment I cannot compile from source with my Windows pc’s so testing the github code will have to wait.

Shijie Yan

unread,
Oct 14, 2021, 5:26:15 PM10/14/21
to mcx-...@googlegroups.com
Hi Gijs,

We are merging the code with the main branch (https://github.com/fangq/mcx/pull/127/commits). Once it is done, the new feature should be available in the MCX nightly build

Thanks again for your contribution,
Shijie 
Monte Carlo eXtreme (MCX) - GPU-accelerated photon transport simulator - Initial support of the hyperboloid gaussian source in MCXLAB by ShijieYan · Pull Request #127 · fangq/mcx


From: 'Buist, G.' via mcx-users <mcx-...@googlegroups.com>
Sent: Wednesday, October 13, 2021 5:49 AM

Qianqian Fang

unread,
Oct 14, 2021, 11:56:58 PM10/14/21
to mcx-...@googlegroups.com, Shijie Yan
On 10/14/21 5:26 PM, Shijie Yan wrote:
Hi Gijs,

We are merging the code with the main branch (https://github.com/fangq/mcx/pull/127/commits). Once it is done, the new feature should be available in the MCX nightly build.


hi Gijs, I have accepted Shijie's pull request, and manually built the windows/linux/mac binaries. feel free to go ahead and download it from our nightly built site and see if it works as expected.


thanks again for making this useful feature available on mcx.


Qianqian


Buist, G.

unread,
Oct 20, 2021, 3:04:31 AM10/20/21
to mcx-...@googlegroups.com, Shijie Yan

Hi Qianqian, Shijie

 

I have installed the nightlybuild and the (true) hyberboloid gaussian seems to work as expected.

Thanks for implementing it in MCX.

 

Best,

Gijs

 

From: mcx-...@googlegroups.com <mcx-...@googlegroups.com> On Behalf Of Qianqian Fang
Sent: vrijdag 15 oktober 2021 05:57
To: mcx-...@googlegroups.com; Shijie Yan <yan....@northeastern.edu>
Subject: Re: [mcx-users] spot focused Gaussian source distribution

 

On 10/14/21 5:26 PM, Shijie Yan wrote:

Hi Gijs,

 

We are merging the code with the main branch (https://github.com/fangq/mcx/pull/127/commits). Once it is done, the new feature should be available in the MCX nightly build.

 

hi Gijs, I have accepted Shijie's pull request, and manually built the windows/linux/mac binaries. feel free to go ahead and download it from our nightly built site and see if it works as expected.

 

thanks again for making this useful feature available on mcx.

 

Qianqian

 

 

Thanks again for your contribution,

Shijie 

Image removed by sender.

Qianqian Fang

unread,
Oct 20, 2021, 2:30:19 PM10/20/21
to mcx-...@googlegroups.com, Shijie Yan
On 10/20/21 3:04 AM, 'Buist, G.' via mcx-users wrote:

Hi Qianqian, Shijie

 

I have installed the nightlybuild and the (true) hyberboloid gaussian seems to work as expected.

Thanks for implementing it in MCX.


glad to hear! thank you for the initial patch and Shijie for implementing/debugging!


Qianqian Fang

unread,
Mar 10, 2024, 2:19:58 PM3/10/24
to mcx-...@googlegroups.com, Shijie Yan

hi Gijs,

in an attempt of reducing register counts in mcx kernels, I optimized the hyperboloid gaussian source implementation by reusing some existing registers. see

https://github.com/fangq/mcx/commit/7b6e0e065ebda62952c1a979f97080d2ad7a72c6

can you please download the latest nightly build version (https://mcx.space/nightly/, or recompile on your side) and see if the source still functions as expected?

thanks

Qianqian

G Buist

unread,
Apr 9, 2024, 4:56:26 AM4/9/24
to mcx-users
Hi Qianqian,

I apologize for the late response. I have not have had the time to test the new optimized hyperboloid Gaussian source and I do not foresee when I will have time to test this.
When I manage to find the time to test I will let you know of the results.

Best,
Gijs

Op zondag 10 maart 2024 om 19:19:58 UTC+1 schreef q.fang:

G Buist

unread,
May 28, 2024, 8:06:54 AM5/28/24
to mcx-users
Hi Qianqian,

I have downloaded the latest version of pmcx and the hyperboloid source seems to work as expected.

Best,
Gijs

Op dinsdag 9 april 2024 om 10:56:26 UTC+2 schreef G Buist:
Reply all
Reply to author
Forward
0 new messages