fangq / mcx

Monte Carlo eXtreme (MCX) - GPU-accelerated photon transport simulator
http://mcx.space
Other
132 stars 73 forks source link

How does the detector work? #102

Closed CGLemon closed 4 years ago

CGLemon commented 4 years ago

Qianqian Fang: Great for the MCX. Let photon simulation be easier and flexible.

In our collage lab, we want to use MCX to replace the old program. After testing, we do not have any idea how does the detector work. In a case, we changed the parameter "ma" to a greater number. we guess that detector would detected lower photons because the photon had been absorbed by media. Out of the blue, the detected photons is same as before. It confusing to me. Another, we also want to know that why the detector is sphere. Is the square detector bad?

I will appreciate your answers. Thank you very much!

                                                                                                  -- Hung Tse, Lin
fangq commented 4 years ago

@CGLemon, let me first explain how an mcx detector captures photons.

there are currently 3 ways to capture photons in mcx/mcxcl

  1. using a detector (cfg.detpos = [dx, dy, dz, Rd], use cfg.savedetflag to decide what data to store). A photon is captured by a detector only if both of the below conditions are met: 1) they are trying to escape from the domain (any non-zero voxel) to the background (zero-valued voxel or outside of the bounding box), and 2) the exit position is within a sphere (radius Rd) centered at the nearest detector.

you can see that these two conditions restrict that a captured photon must be located on the bounding box/an interface between zero and non-zero voxels, truncated by a sphere centered at a given detector.

if the above conditions are met, a photon will be captured, regardless of their propagation direction. you can, however, set cfg.issaveexit=1 or append letters 'xv' in cfg.savedetflag to store the exit position detpt.p and exit direction vector detpt.v. If your actual physical detector has an NA that limits the incident photon angles, you can then use detp.{p,v} to figure out if it will be detected by your camera/optics in a post processing.

to show what this means, here is the plot of the detected photon positions when a series of detectors is placed on a curved volume, you can see that they are limited to the voxel boundaries.

curved_bc_det

  1. in the latest github version of mcx/mcxcl, I added a new way to detect photons on the bounding box. One can set digits 7-12 in the cfg.bc='______000000' input to ask mcx to capture all photons escaping from the 6 bounding box facets: x=0,y=0,z=0,x=Nx, y=Ny, z=Nz planes, respectively according to the orders of the letters - a letter '1' indicates that you want to capture all photons escaping from that plane, a letter '0' - do not capture. This way, you can define a rectangle-shaped detector but only on the bounding box.the output data are similar to a spherical detector, where you can use cfg.savedetflag to get various info, including exit position and angle. More details, please see

https://github.com/fangq/mcx/blob/master/README.txt#L244-L250

you can also run this sample script to get a demo: https://github.com/fangq/mcx/blob/master/mcxlab/examples/demo_bc_det.m

bc_det

  1. you can use cfg.issaveref to accumulate total diffuse reflectance within zero-valued voxels immediately adjacent to a non-zero valued voxel. This gives you the distribution of light leaking from the entire surface of the object, as long as you pad a layer of zeros around the object. This accumulated surface light distribution does not store individually detected photon information, instead, their final weight are summed together regardless of their escaping directions. This is more memory efficient, but does not provide the granularity if you want to individually process the detected photon.

Lastly, I want to highlight that both mcx and mmc have the capability of modeling wide-field detectors (i.e. non contact), as long as you discretize the space between the detector and the domain and include that as part of your volume (in mcx) or mesh (in mmc). An example of this can be seen for mmc in Fig. 7 of our wide-field mmc paper - a rectangular detector can be placed outside of the object as long as you tessellate the gap between them.

mmc_widefield_det

Please be aware that any empty space between your widefield source/detector in mcx/mcxcl, if they are rasterized, must be labeled as a non-zero label with an optical property of air, instead of setting it to zero-valued voxel. This is because zero-valued voxels are special in mcx/mcxcl as the background, as I mentioned above, a photon moving from a non-zero voxel to a zero-voxel will be terminated, so they won't make it to the wide-field detector.

let me know if these addressed your questions.

CGLemon commented 4 years ago

Thank you for your reply. So, I got these points.

1). Detector can only exist at non-zero voxel and zero-valued voxel boundary. 2). A photon will be captured when it meets detector and doesn't run away on next step. 3). It is more effective that detector is sphere. 4). You developed a effective way to set a rectangle(wild-field) detector and source in MCX/MMC.

Because I do not specialize in this, I do not understand how hard to implement these functions. However, it looks more powerful than I had thought.

But I still have no idea why it is same photons captured by detector in different case. For example, I set a box in the grid. the
source is below grid. The detector is upper grid. It looks like below.

Screenshot from 2020-08-29 12-09-18

I set different parameters "ua" for box. I thought this would change the number of captured photons. However it is the same number of captured photons. Is it just because that the photon still exist when it is no energy?

fangq commented 4 years ago

@CGLemon, the reason you don't see photon number changes when changing mua is because mcx simulates photon packets instead of individual photons. You must call mcxdetweight.m to get the detected photon weights. It requires you to provide cfg.prop and use the updated mua for this calculation.

please checkout this mailing list post and the associated thread for additional details

https://groups.google.com/g/mcx-users/c/nxP4BKkTSY4/m/YslmOj7gCAAJ

in mcx/mmc, mua only determines weight attenuation and mus only determines photon trajectories, which is different from MCML. This method is referred to as the mBLL method in Angelo Sassaroli's JOSA paper. Based on that paper, mBLL provides better simulations than the Albedo-Weight approach used in MCML.

CGLemon commented 4 years ago

Spending some time to understand it. It solves my questions. Thank you very much.